/*****************************************************************************
*                                                                            *
*   SQRTM.SPL    Copyright (C) 2011 DSP Development Corporation              *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Computes the matrix square root                             *
*                                                                            *
*   Revisions:   20 Apr 2011  RRR  Creation                                  *
*                                                                            *
*****************************************************************************/


#if @HELP_SQRTM

    SQRTM

    Purpose: Computes the principal square root of a square matrix.

    Syntax:  SQRTM(a)

             (x, res) = SQRTM(a)

              a - input array, must be square


    Returns: An N x N array.
               
             (x, res) = SQRTM(a) returns the square root matrix and residual
             where res = norm(a - x *^ x, "fro") / norm(a, "fro").

    Example:
             W1: ravel(1..9, 3)
             W2: sqrtm(W1);
             W3: real(W2 *^ W2)

             W1 contains the array: 

              {{1, 4, 7}.
               {2, 5, 8},
               {3, 6, 9}}

             W2 == 

              {{0.4498 + 0.7623i, 1.0185 + 0.0842i, 1.5873 - 0.5940i},
               {0.5526 + 0.2068i, 1.2515 + 0.0229i, 1.9503 - 0.1611i},
               {0.6555 - 0.3487i, 1.4844 - 0.0385i, 2.3134 + 0.2717i}}

             W3 ==

              {{1, 4, 7},
               {2, 5, 8},
               {3, 6, 9}}

    Remarks:
             SQRTM(A) returns an array X such that X *^ X == A.

             Not all matrices have a square root and some matrices have
             multiple square roots.

    See Also:
             Norm
             Eig
             Expm
             Funm
             Logm
#endif


/* matrix square root */
sqrtm(a, tol)
{
        local nr, nc, x, res;

        if (argc < 2)
        {
                if (argc < 1) error("sqrtm - matrix required");

                tol = 10 * eps;
        }

        (nr, nc) = size(a);

        if (nr != nc) error("sqrtm - matrix must be square");

        /* sqrt */
        x = a ^^ 0.5;

        /* check if purely real */
        if (isreal(a) && all(all(abs(imag(x)) < tol)))
        {
                x = real(x);
        }

        if (outargc > 1)
        {
                /* return residual */
                res = norm(x *^ x - a, "fro") / norm(a, "fro");
                return(x, res);
        }
        else
        {
                return(x);
        }
}