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

#if @HELP_FUNM

    FUNM

    Purpose: Computes the general matrix function of a square matrix.

    Syntax:  FUNM(a, "fun", tol)

                 a - A square matrix.

             "fun" - A string, the matrix function.

               tol - Optional. A real, the imaginary part tolerence.
                     Defaults to 1e-9.

    Returns: A square matrix.
               
    Example:
             W1: ravel(1..9, 3)
             W2: funm(W1, "sinh")
             W3: funm(W2, "asinh")

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

             W2 == {{559452.165, 1266940.566, 1974428.967},
                    {687407.399, 1556707.149, 2426006.900},
                    {815362.634, 1846473.733, 2877584.832}}

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

             W2 computes the matrix sinh function and W3 recovers the
             original matrix by computing the matrix asinh function.

    Remarks:
             The function string can refer to any unary built-in or
             user defined function that operates on scalar values.

             The imaginary component of the result is set to zero if
             all values are less than TOL and the input was purely
             real.

   See Also:
             Expm
             Logm
             Sqrtm
#endif


/* generic matrix form of a given function */
funm(A, fun, tol)
{
        local V, D, F, scalar = 0;

        if (argc < 3)
        {
                if (argc < 2) error("funm - input matrix and function string required");

                tol = 1e-9;
        }

        if (isscalar(A))
        {
                A      = {A};
                scalar = 1;
        }

        if (numcols(A) != numrows(A))
        {
                error("funm - square matrix required");
        }

        if (not(isstring(fun))) error("funm - function string required");

        /* compute matrix function */

        (V, D) = eig(A);

        F = V *^ diag(feval(fun, diag(D))) /^ V;

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

        if (scalar)
        {
                F = F[1];
        }
        else
        {
                setvunits(F, getvunits(A, 1), -1);
        }

        return(F);
}