/*****************************************************************************
*                                                                            *
*   COND.SPL     Copyright (C) 2001-2011 DSP Development Corporation         *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Estimates the condition number of a matrix                  *
*                                                                            *
*   Revisions:   24 Jul 2001  RRR  Creation                                  *
*                 6 May 2005  RRR  override default_math_value to return inf *
*                26 Aug 2011  RRR  p norm                                    *
*                                                                            *
*****************************************************************************/


#if @HELP_COND

    COND

    Purpose: Estimates the condition number of a matrix

    Syntax:  COND(a, p)

                  a - input matrix, defaults to the current series

                  p - optional integer, the p-norm. Defaults to 2,
                      singular value method.

    Returns: A real, the estimated condition number.

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

             W2: {{1, 2, 3},
                  {4, 5, 6},
                  {7, 8, 0}}


             cond(w1) == 8.579581E+016

             cond(w2) == 35.105870

    Remarks:
             COND uses SVD to compute the singular values of a matrix. The
             result is the ratio of the maximum and minimum singular values.

             The condition number is infinite for a singular matrix,
             thus a large condition number indicates an ill conditioned
             matrix.

             If p is specified, the condition number is computed as:

                 c = norm(a, p) * norm(inv(a), p)

             See RCOND to compute the reciprocal condition number.

    See Also:
             Norm
             Rcond
             Svd
#endif


/* condition number */
cond(a, p)
{
        local s, c, defmv;

        if (argc < 2)
        {
                if (argc < 1) a = refseries(1);

                p = 2;
        }

        /* if not an array, cast */
        if (not(isarray(a))) a = {a};

        if (isscalar(p) && castint(p) == 2)
        {
                /* singular values */
                s = svd(a);

                /* current default math value */
                defmv = getconf("default_math_value");

                /* override to return inf for singular matrices */
                setconf("default_math_value", "inf");

                /* condition value */
                c = max(s) / min(s);

                /* restore */
                setconf("default_math_value", defmv);
        }
        else
        {
                c = norm(a, p) * norm(inv(a), p);
        }

        return(c);

}


/* if no matrix inverse, return inf */
cond_error()
{
        return(inf);
}