/*****************************************************************************
*                                                                            *
*   LU.SPL       Copyright (C) 2011 DSP Development Corporation              *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    LU decomposition of a matrix, wrapper to built-in _lu       *
*                                                                            *
*   Revisions:   24 Aug 2011  RRR  Creation                                  *
*                                                                            *
*****************************************************************************/

#include 

#if @HELP_LU

    Purpose: Computes the LU decomposition of a square matrix.

    Syntax:  LU(a, type, permute)

             (l, u) = LU(a)

             (l, u, p) = LU(a, shape)

                    a -  Real or complex matrix.
 
                 type -  Optional. An integer specifying the type of matrix.
                         Defaults to 2. Valid arguments are:

                           0: Lower LU decomposition matrix.

                           1: Upper LU decomposition matrix (default).

                           2: Complete LU decomposition matrix.

              permute -  Optional. An integer specifying whether or not to
                         permute. 

                           0: No permutation

                           1: permute (default).

                shape -  Optional string if returning the permutation value.
                         "vector" returns the permutation values as a single
                         column series. Defaults to a matrix of 1 and 0.


    Returns: A matrix, the full LU decomposition or an individual LU
             component.

             (l, u) = LU(a) returns the lower and upper components such
             that l *^ u == a.

             (l, u, p) = LU(a) returns the lower, upper and permutation
             matrices such that l *^ u == p *^ a.

    Example:
              a = {{1, 2, 3},
                   {4, 5, 6},
                   {7, 8, 10}};

              m = lu(a);

              (l, u) = lu(a);


              m == {{7,        8,       10},
                    {0.142857, 0.857143, 1.57143},
                    {0.571429, 0.5,     -0.5}}

              l == {{0.142857, 1,        0},
                    {0.571429, 0.5,      1},
                    {1,        0,        0}}

              u == {{7,        8,       10},
                    {0,        0.857143, 1.57143},
                    {0,        0,       -0.5}}

              
              l *^ u == {{1, 2, 3},
                         {4, 5, 6},
                         {7, 8, 10}}

             l is a permuted lower triangle matrix and u is the corresponding
             upper triangle matrix such that l *^ y == a.

    Example:
             (l1, u1, p) = lu(a);

             b = l1 *^ u1;
             c = p  *^ a; 

             l == {{1,         0,        0},
                   {0.142857,  1,        0}
                   {0.571429,  0.5,      1}

             u == {{7,         8,       10},
                   {0,         0.857143, 1.57143},
                   {0,         0,       -0.5}}

             p == {{0,         0,        1},
                   {1,         0,        0},
                   {0,         1,        0}}

              
             b == c == {{7, 8, 10},
                        {1, 2,  3},
                        {4, 5,  6}}


             l is now a true lower triangular matrix and u is the same as
             the previous example. The permutation matrix, p, is formed such
             that l *^ u == p *^ a, a permutation of the original matrix.

    Remarks:
             The built in \^ operator uses LU decomposition for system of
             equations with a square matrix. In this case, the equation
             A *^ x = b can be solved with:

               y = l \^ b
               x = u \^ y

             where the solutions to y and x are computed by simple back
             substitution.

             INV invokes LU decompostion to compute the inverse of a square
             matrix with:

              inv(u) *^ inv(l) == inv(a)

    See Also:
             \^
             *^
             det
             inv
             qr
#endif


/* LU decomposition */
lu(a, type, permute)
{
        if (argc < 3)
        {
                if (argc < 2)
                {
                        if (argc < 1) error("lu - input matrix required");

                        /* full matrix if returing a single matrix */
                        type = 2;
                }

                permute = 1;
        }

        if (outargc > 1)
        {
                if (outargc > 2)
                {
                        /* return permutation matrix or vector */
                        (l, u, p) = _lu(a, type);

                        return(l, u, p);
                }
                else
                {
                        /* return each component */
                        (l, u) = _lu(a, type);

                        return(l, u);
                }
        }
        else
        {
                /* return full matrix or a single component */
                return(_lu(a, type, permute));
        }
}