/*****************************************************************************
*                                                                            *
*   SVD.SPL      Copyright (C) 1999-2012 DSP Development Corporation         *
*                                All Rights Reserved                         *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Calculates the singular value decomposition of an array     *
*                                                                            *
*   Revisions:   13 Jun 1999  RRR  Creation                                  *
*                18 Jan 2010  RRR  compact form                              *
*                13 Dec 2012  RRR  economy form                              *
*                                                                            *
*****************************************************************************/

#include 

#if @HELP_SVD

    SVD

    Purpose: Calculates the singular value decomposition of an array

    Syntax:  (u, w, v) = SVD(a, mode)

               a    - the input array

               mode - an optional integer, type of output matrices:
                  
                       0 - Compact form.
                       1 - Standard form.

    Alternative Syntax:  

             SVD(a, mode)

               a    - the input array

               mode - an optional integer, type of output matrix:
                 
                       00 - W, Singular values as a column (default).
                       01 - V, Right singular value matrix.
                       10 - U, Left singular value matrix.
                       11 - Combined U W V matrix.

    Returns:

             (U, W, V) = SVD(A)

             returns U, W and V as three separate matrices. For input
             matrix A of size NxM, matrix U is size NxN, W is size

             SVD(A) returns a series containing the singular values of A.

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

             (u, w, v) = svd(w1)

             u == {{-0.1525, -0.8226, -0.3945, -0.3800},
                   {-0.3499, -0.4214,  0.2428,  0.8007},
                   {-0.5474, -0.0201,  0.6979, -0.4614},
                   {-0.7448,  0.3812, -0.5462,  0.0407}}

             w == {{14.2691,  0.0000},
                   { 0.0000,  0.6268},
                   { 0.0000,  0.0000},
                   { 0.0000,  0.0000}}

             v == {{-0.6414, -0.7672},
                   {-0.7672, -0.6414}}

             u *^ w *^ v' == {{1, 2},
                              {3, 4},
                              {5, 6},
                              {7, 8}}


             SVD computes the decomposition in standard form.

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

             (u, w, v) = svd(w1, 0)

             u == {{-0.1525, -0.8226},
                   {-0.3499, -0.4214},
                   {-0.5474, -0.0201},
                   {-0.7448,  0.3812}}

             w == {{14.2691,  0.0000},
                   { 0.0000,  0.6268}}

             v == {{-0.6414, -0.7672},
                   {-0.7672, -0.6414}}

             u *^ w *^ v' == {{1, 2},
                              {3, 4},
                              {5, 6},
                              {7, 8}}


             SVD computes the decomposition in compact form.

    Example:
             W1: {{1, 2},
                  {3, 4},
                  {5, 6},
                  {7, 8}}
     
             W2: svd(w1)

             W2 == {14.2691, 
                     0.6268}


             W2 contains the singular values of W1 as a single column
             series.  The singular values are ranked from largest to
             smallest.

    Remarks:
             The input matrix A is decomposed into a left singular value
             matrix U, a diagonal matrix W, and a right singular matrix
             V such that:

             A = U *^ W *^ V'

             The inverse of the matrix can be calculated as:

             V *^ (diag(1/diag(W)) *^ U'


             For input matrix A of size NxM, matrix U is size NxN, W is
             size NxM and V is size MxM.

             (u, w, v) = svd(a, 0)

             returns the components in compact form where U is size NxM, 
             W is size MxM and V is size MxM.

    See Also:
             Linfit
             QR
             Svddiv

    References:

             [1] Press, Flannery, Teukolsky, Vetterling
                 Numerical Recipes in C
                 Cambridge Press, 1988
                 pp. 407-552
#endif


/* singular value decomposition - wrapper */
svd(a, mode)
{
        local u, s, v;

        if (outargc > 1)
        {
                /* generate full components */
                if (argc < 2) mode = 11;

                if (isstring(mode))
                {
                        /* economy form */
                        mode = -1;
                }
                else if (mode > 0)
                {
                        /* standard form */
                        mode = 11;
                }
                else if (mode < 0)
                {
                        /* economy form, note mode == 0 implies compact form */
                        mode = -1;
                }

                /* svd calculation */
                (u, s, v) = _svd(a, mode);

                /* return all components */
                return(u, s, v);
        }
        else if (argc > 1)
        {
                /* return a single component */
                return(_svd(a, mode));
        }
        else
        {
                /* return singular values as a column */
                return(_svd(a));
        }
}