source: trunk/sources/ALGLIB/svd.cs @ 2430

/*************************************************************************
Copyright (c) 2005-2007, Sergey Bochkanov (ALGLIB project).

>>> SOURCE LICENSE >>>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation (www.fsf.org); either version 2 of the 
License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

A copy of the GNU General Public License is available at
http://www.fsf.org/licensing/licenses

>>> END OF LICENSE >>>
*************************************************************************/

using System;

namespace alglib
{
    public class svd
    {
        /*************************************************************************
        Singular value decomposition of a rectangular matrix.

        The algorithm calculates the singular value decomposition of a matrix of
        size MxN: A = U * S * V^T

        The algorithm finds the singular values and, optionally, matrices U and V^T.
        The algorithm can find both first min(M,N) columns of matrix U and rows of
        matrix V^T (singular vectors), and matrices U and V^T wholly (of sizes MxM
        and NxN respectively).

        Take into account that the subroutine does not return matrix V but V^T.

        Input parameters:
            A           -   matrix to be decomposed.
                            Array whose indexes range within [0..M-1, 0..N-1].
            M           -   number of rows in matrix A.
            N           -   number of columns in matrix A.
            UNeeded     -   0, 1 or 2. See the description of the parameter U.
            VTNeeded    -   0, 1 or 2. See the description of the parameter VT.
            AdditionalMemory -
                            If the parameter:
                             * equals 0, the algorithm doesn’t use additional
                               memory (lower requirements, lower performance).
                             * equals 1, the algorithm uses additional
                               memory of size min(M,N)*min(M,N) of real numbers.
                               It often speeds up the algorithm.
                             * equals 2, the algorithm uses additional
                               memory of size M*min(M,N) of real numbers.
                               It allows to get a maximum performance.
                            The recommended value of the parameter is 2.

        Output parameters:
            W           -   contains singular values in descending order.
            U           -   if UNeeded=0, U isn't changed, the left singular vectors
                            are not calculated.
                            if Uneeded=1, U contains left singular vectors (first
                            min(M,N) columns of matrix U). Array whose indexes range
                            within [0..M-1, 0..Min(M,N)-1].
                            if UNeeded=2, U contains matrix U wholly. Array whose
                            indexes range within [0..M-1, 0..M-1].
            VT          -   if VTNeeded=0, VT isn’t changed, the right singular vectors
                            are not calculated.
                            if VTNeeded=1, VT contains right singular vectors (first
                            min(M,N) rows of matrix V^T). Array whose indexes range
                            within [0..min(M,N)-1, 0..N-1].
                            if VTNeeded=2, VT contains matrix V^T wholly. Array whose
                            indexes range within [0..N-1, 0..N-1].

          -- ALGLIB --
             Copyright 2005 by Bochkanov Sergey
        *************************************************************************/
        public static bool rmatrixsvd(double[,] a,
            int m,
            int n,
            int uneeded,
            int vtneeded,
            int additionalmemory,
            ref double[] w,
            ref double[,] u,
            ref double[,] vt)
        {
            bool result = new bool();
            double[] tauq = new double[0];
            double[] taup = new double[0];
            double[] tau = new double[0];
            double[] e = new double[0];
            double[] work = new double[0];
            double[,] t2 = new double[0,0];
            bool isupper = new bool();
            int minmn = 0;
            int ncu = 0;
            int nrvt = 0;
            int nru = 0;
            int ncvt = 0;
            int i = 0;
            int j = 0;
            int im1 = 0;
            double sm = 0;

            a = (double[,])a.Clone();

            result = true;
            if( m==0 | n==0 )
            {
                return result;
            }
            System.Diagnostics.Debug.Assert(uneeded>=0 & uneeded<=2, "SVDDecomposition: wrong parameters!");
            System.Diagnostics.Debug.Assert(vtneeded>=0 & vtneeded<=2, "SVDDecomposition: wrong parameters!");
            System.Diagnostics.Debug.Assert(additionalmemory>=0 & additionalmemory<=2, "SVDDecomposition: wrong parameters!");
            
            //
            // initialize
            //
            minmn = Math.Min(m, n);
            w = new double[minmn+1];
            ncu = 0;
            nru = 0;
            if( uneeded==1 )
            {
                nru = m;
                ncu = minmn;
                u = new double[nru-1+1, ncu-1+1];
            }
            if( uneeded==2 )
            {
                nru = m;
                ncu = m;
                u = new double[nru-1+1, ncu-1+1];
            }
            nrvt = 0;
            ncvt = 0;
            if( vtneeded==1 )
            {
                nrvt = minmn;
                ncvt = n;
                vt = new double[nrvt-1+1, ncvt-1+1];
            }
            if( vtneeded==2 )
            {
                nrvt = n;
                ncvt = n;
                vt = new double[nrvt-1+1, ncvt-1+1];
            }
            
            //
            // M much larger than N
            // Use bidiagonal reduction with QR-decomposition
            //
            if( m>1.6*n )
            {
                if( uneeded==0 )
                {
                    
                    //
                    // No left singular vectors to be computed
                    //
                    qr.rmatrixqr(ref a, m, n, ref tau);
                    for(i=0; i<=n-1; i++)
                    {
                        for(j=0; j<=i-1; j++)
                        {
                            a[i,j] = 0;
                        }
                    }
                    bidiagonal.rmatrixbd(ref a, n, n, ref tauq, ref taup);
                    bidiagonal.rmatrixbdunpackpt(ref a, n, n, ref taup, nrvt, ref vt);
                    bidiagonal.rmatrixbdunpackdiagonals(ref a, n, n, ref isupper, ref w, ref e);
                    result = bdsvd.rmatrixbdsvd(ref w, e, n, isupper, false, ref u, 0, ref a, 0, ref vt, ncvt);
                    return result;
                }
                else
                {
                    
                    //
                    // Left singular vectors (may be full matrix U) to be computed
                    //
                    qr.rmatrixqr(ref a, m, n, ref tau);
                    qr.rmatrixqrunpackq(ref a, m, n, ref tau, ncu, ref u);
                    for(i=0; i<=n-1; i++)
                    {
                        for(j=0; j<=i-1; j++)
                        {
                            a[i,j] = 0;
                        }
                    }
                    bidiagonal.rmatrixbd(ref a, n, n, ref tauq, ref taup);
                    bidiagonal.rmatrixbdunpackpt(ref a, n, n, ref taup, nrvt, ref vt);
                    bidiagonal.rmatrixbdunpackdiagonals(ref a, n, n, ref isupper, ref w, ref e);
                    if( additionalmemory<1 )
                    {
                        
                        //
                        // No additional memory can be used
                        //
                        bidiagonal.rmatrixbdmultiplybyq(ref a, n, n, ref tauq, ref u, m, n, true, false);
                        result = bdsvd.rmatrixbdsvd(ref w, e, n, isupper, false, ref u, m, ref a, 0, ref vt, ncvt);
                    }
                    else
                    {
                        
                        //
                        // Large U. Transforming intermediate matrix T2
                        //
                        work = new double[Math.Max(m, n)+1];
                        bidiagonal.rmatrixbdunpackq(ref a, n, n, ref tauq, n, ref t2);
                        blas.copymatrix(ref u, 0, m-1, 0, n-1, ref a, 0, m-1, 0, n-1);
                        blas.inplacetranspose(ref t2, 0, n-1, 0, n-1, ref work);
                        result = bdsvd.rmatrixbdsvd(ref w, e, n, isupper, false, ref u, 0, ref t2, n, ref vt, ncvt);
                        blas.matrixmatrixmultiply(ref a, 0, m-1, 0, n-1, false, ref t2, 0, n-1, 0, n-1, true, 1.0, ref u, 0, m-1, 0, n-1, 0.0, ref work);
                    }
                    return result;
                }
            }
            
            //
            // N much larger than M
            // Use bidiagonal reduction with LQ-decomposition
            //
            if( n>1.6*m )
            {
                if( vtneeded==0 )
                {
                    
                    //
                    // No right singular vectors to be computed
                    //
                    lq.rmatrixlq(ref a, m, n, ref tau);
                    for(i=0; i<=m-1; i++)
                    {
                        for(j=i+1; j<=m-1; j++)
                        {
                            a[i,j] = 0;
                        }
                    }
                    bidiagonal.rmatrixbd(ref a, m, m, ref tauq, ref taup);
                    bidiagonal.rmatrixbdunpackq(ref a, m, m, ref tauq, ncu, ref u);
                    bidiagonal.rmatrixbdunpackdiagonals(ref a, m, m, ref isupper, ref w, ref e);
                    work = new double[m+1];
                    blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
                    result = bdsvd.rmatrixbdsvd(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref vt, 0);
                    blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
                    return result;
                }
                else
                {
                    
                    //
                    // Right singular vectors (may be full matrix VT) to be computed
                    //
                    lq.rmatrixlq(ref a, m, n, ref tau);
                    lq.rmatrixlqunpackq(ref a, m, n, ref tau, nrvt, ref vt);
                    for(i=0; i<=m-1; i++)
                    {
                        for(j=i+1; j<=m-1; j++)
                        {
                            a[i,j] = 0;
                        }
                    }
                    bidiagonal.rmatrixbd(ref a, m, m, ref tauq, ref taup);
                    bidiagonal.rmatrixbdunpackq(ref a, m, m, ref tauq, ncu, ref u);
                    bidiagonal.rmatrixbdunpackdiagonals(ref a, m, m, ref isupper, ref w, ref e);
                    work = new double[Math.Max(m, n)+1];
                    blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
                    if( additionalmemory<1 )
                    {
                        
                        //
                        // No additional memory available
                        //
                        bidiagonal.rmatrixbdmultiplybyp(ref a, m, m, ref taup, ref vt, m, n, false, true);
                        result = bdsvd.rmatrixbdsvd(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref vt, n);
                    }
                    else
                    {
                        
                        //
                        // Large VT. Transforming intermediate matrix T2
                        //
                        bidiagonal.rmatrixbdunpackpt(ref a, m, m, ref taup, m, ref t2);
                        result = bdsvd.rmatrixbdsvd(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref t2, m);
                        blas.copymatrix(ref vt, 0, m-1, 0, n-1, ref a, 0, m-1, 0, n-1);
                        blas.matrixmatrixmultiply(ref t2, 0, m-1, 0, m-1, false, ref a, 0, m-1, 0, n-1, false, 1.0, ref vt, 0, m-1, 0, n-1, 0.0, ref work);
                    }
                    blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
                    return result;
                }
            }
            
            //
            // M<=N
            // We can use inplace transposition of U to get rid of columnwise operations
            //
            if( m<=n )
            {
                bidiagonal.rmatrixbd(ref a, m, n, ref tauq, ref taup);
                bidiagonal.rmatrixbdunpackq(ref a, m, n, ref tauq, ncu, ref u);
                bidiagonal.rmatrixbdunpackpt(ref a, m, n, ref taup, nrvt, ref vt);
                bidiagonal.rmatrixbdunpackdiagonals(ref a, m, n, ref isupper, ref w, ref e);
                work = new double[m+1];
                blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
                result = bdsvd.rmatrixbdsvd(ref w, e, minmn, isupper, false, ref a, 0, ref u, nru, ref vt, ncvt);
                blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
                return result;
            }
            
            //
            // Simple bidiagonal reduction
            //
            bidiagonal.rmatrixbd(ref a, m, n, ref tauq, ref taup);
            bidiagonal.rmatrixbdunpackq(ref a, m, n, ref tauq, ncu, ref u);
            bidiagonal.rmatrixbdunpackpt(ref a, m, n, ref taup, nrvt, ref vt);
            bidiagonal.rmatrixbdunpackdiagonals(ref a, m, n, ref isupper, ref w, ref e);
            if( additionalmemory<2 | uneeded==0 )
            {
                
                //
                // We cant use additional memory or there is no need in such operations
                //
                result = bdsvd.rmatrixbdsvd(ref w, e, minmn, isupper, false, ref u, nru, ref a, 0, ref vt, ncvt);
            }
            else
            {
                
                //
                // We can use additional memory
                //
                t2 = new double[minmn-1+1, m-1+1];
                blas.copyandtranspose(ref u, 0, m-1, 0, minmn-1, ref t2, 0, minmn-1, 0, m-1);
                result = bdsvd.rmatrixbdsvd(ref w, e, minmn, isupper, false, ref u, 0, ref t2, m, ref vt, ncvt);
                blas.copyandtranspose(ref t2, 0, minmn-1, 0, m-1, ref u, 0, m-1, 0, minmn-1);
            }
            return result;
        }


        /*************************************************************************
        Obsolete 1-based subroutine.
        See RMatrixSVD for 0-based replacement.
        *************************************************************************/
        public static bool svddecomposition(double[,] a,
            int m,
            int n,
            int uneeded,
            int vtneeded,
            int additionalmemory,
            ref double[] w,
            ref double[,] u,
            ref double[,] vt)
        {
            bool result = new bool();
            double[] tauq = new double[0];
            double[] taup = new double[0];
            double[] tau = new double[0];
            double[] e = new double[0];
            double[] work = new double[0];
            double[,] t2 = new double[0,0];
            bool isupper = new bool();
            int minmn = 0;
            int ncu = 0;
            int nrvt = 0;
            int nru = 0;
            int ncvt = 0;
            int i = 0;
            int j = 0;
            int im1 = 0;
            double sm = 0;

            a = (double[,])a.Clone();

            result = true;
            if( m==0 | n==0 )
            {
                return result;
            }
            System.Diagnostics.Debug.Assert(uneeded>=0 & uneeded<=2, "SVDDecomposition: wrong parameters!");
            System.Diagnostics.Debug.Assert(vtneeded>=0 & vtneeded<=2, "SVDDecomposition: wrong parameters!");
            System.Diagnostics.Debug.Assert(additionalmemory>=0 & additionalmemory<=2, "SVDDecomposition: wrong parameters!");
            
            //
            // initialize
            //
            minmn = Math.Min(m, n);
            w = new double[minmn+1];
            ncu = 0;
            nru = 0;
            if( uneeded==1 )
            {
                nru = m;
                ncu = minmn;
                u = new double[nru+1, ncu+1];
            }
            if( uneeded==2 )
            {
                nru = m;
                ncu = m;
                u = new double[nru+1, ncu+1];
            }
            nrvt = 0;
            ncvt = 0;
            if( vtneeded==1 )
            {
                nrvt = minmn;
                ncvt = n;
                vt = new double[nrvt+1, ncvt+1];
            }
            if( vtneeded==2 )
            {
                nrvt = n;
                ncvt = n;
                vt = new double[nrvt+1, ncvt+1];
            }
            
            //
            // M much larger than N
            // Use bidiagonal reduction with QR-decomposition
            //
            if( m>1.6*n )
            {
                if( uneeded==0 )
                {
                    
                    //
                    // No left singular vectors to be computed
                    //
                    qr.qrdecomposition(ref a, m, n, ref tau);
                    for(i=2; i<=n; i++)
                    {
                        for(j=1; j<=i-1; j++)
                        {
                            a[i,j] = 0;
                        }
                    }
                    bidiagonal.tobidiagonal(ref a, n, n, ref tauq, ref taup);
                    bidiagonal.unpackptfrombidiagonal(ref a, n, n, ref taup, nrvt, ref vt);
                    bidiagonal.unpackdiagonalsfrombidiagonal(ref a, n, n, ref isupper, ref w, ref e);
                    result = bdsvd.bidiagonalsvddecomposition(ref w, e, n, isupper, false, ref u, 0, ref a, 0, ref vt, ncvt);
                    return result;
                }
                else
                {
                    
                    //
                    // Left singular vectors (may be full matrix U) to be computed
                    //
                    qr.qrdecomposition(ref a, m, n, ref tau);
                    qr.unpackqfromqr(ref a, m, n, ref tau, ncu, ref u);
                    for(i=2; i<=n; i++)
                    {
                        for(j=1; j<=i-1; j++)
                        {
                            a[i,j] = 0;
                        }
                    }
                    bidiagonal.tobidiagonal(ref a, n, n, ref tauq, ref taup);
                    bidiagonal.unpackptfrombidiagonal(ref a, n, n, ref taup, nrvt, ref vt);
                    bidiagonal.unpackdiagonalsfrombidiagonal(ref a, n, n, ref isupper, ref w, ref e);
                    if( additionalmemory<1 )
                    {
                        
                        //
                        // No additional memory can be used
                        //
                        bidiagonal.multiplybyqfrombidiagonal(ref a, n, n, ref tauq, ref u, m, n, true, false);
                        result = bdsvd.bidiagonalsvddecomposition(ref w, e, n, isupper, false, ref u, m, ref a, 0, ref vt, ncvt);
                    }
                    else
                    {
                        
                        //
                        // Large U. Transforming intermediate matrix T2
                        //
                        work = new double[Math.Max(m, n)+1];
                        bidiagonal.unpackqfrombidiagonal(ref a, n, n, ref tauq, n, ref t2);
                        blas.copymatrix(ref u, 1, m, 1, n, ref a, 1, m, 1, n);
                        blas.inplacetranspose(ref t2, 1, n, 1, n, ref work);
                        result = bdsvd.bidiagonalsvddecomposition(ref w, e, n, isupper, false, ref u, 0, ref t2, n, ref vt, ncvt);
                        blas.matrixmatrixmultiply(ref a, 1, m, 1, n, false, ref t2, 1, n, 1, n, true, 1.0, ref u, 1, m, 1, n, 0.0, ref work);
                    }
                    return result;
                }
            }
            
            //
            // N much larger than M
            // Use bidiagonal reduction with LQ-decomposition
            //
            if( n>1.6*m )
            {
                if( vtneeded==0 )
                {
                    
                    //
                    // No right singular vectors to be computed
                    //
                    lq.lqdecomposition(ref a, m, n, ref tau);
                    for(i=1; i<=m-1; i++)
                    {
                        for(j=i+1; j<=m; j++)
                        {
                            a[i,j] = 0;
                        }
                    }
                    bidiagonal.tobidiagonal(ref a, m, m, ref tauq, ref taup);
                    bidiagonal.unpackqfrombidiagonal(ref a, m, m, ref tauq, ncu, ref u);
                    bidiagonal.unpackdiagonalsfrombidiagonal(ref a, m, m, ref isupper, ref w, ref e);
                    work = new double[m+1];
                    blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
                    result = bdsvd.bidiagonalsvddecomposition(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref vt, 0);
                    blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
                    return result;
                }
                else
                {
                    
                    //
                    // Right singular vectors (may be full matrix VT) to be computed
                    //
                    lq.lqdecomposition(ref a, m, n, ref tau);
                    lq.unpackqfromlq(ref a, m, n, ref tau, nrvt, ref vt);
                    for(i=1; i<=m-1; i++)
                    {
                        for(j=i+1; j<=m; j++)
                        {
                            a[i,j] = 0;
                        }
                    }
                    bidiagonal.tobidiagonal(ref a, m, m, ref tauq, ref taup);
                    bidiagonal.unpackqfrombidiagonal(ref a, m, m, ref tauq, ncu, ref u);
                    bidiagonal.unpackdiagonalsfrombidiagonal(ref a, m, m, ref isupper, ref w, ref e);
                    work = new double[Math.Max(m, n)+1];
                    blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
                    if( additionalmemory<1 )
                    {
                        
                        //
                        // No additional memory available
                        //
                        bidiagonal.multiplybypfrombidiagonal(ref a, m, m, ref taup, ref vt, m, n, false, true);
                        result = bdsvd.bidiagonalsvddecomposition(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref vt, n);
                    }
                    else
                    {
                        
                        //
                        // Large VT. Transforming intermediate matrix T2
                        //
                        bidiagonal.unpackptfrombidiagonal(ref a, m, m, ref taup, m, ref t2);
                        result = bdsvd.bidiagonalsvddecomposition(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref t2, m);
                        blas.copymatrix(ref vt, 1, m, 1, n, ref a, 1, m, 1, n);
                        blas.matrixmatrixmultiply(ref t2, 1, m, 1, m, false, ref a, 1, m, 1, n, false, 1.0, ref vt, 1, m, 1, n, 0.0, ref work);
                    }
                    blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
                    return result;
                }
            }
            
            //
            // M<=N
            // We can use inplace transposition of U to get rid of columnwise operations
            //
            if( m<=n )
            {
                bidiagonal.tobidiagonal(ref a, m, n, ref tauq, ref taup);
                bidiagonal.unpackqfrombidiagonal(ref a, m, n, ref tauq, ncu, ref u);
                bidiagonal.unpackptfrombidiagonal(ref a, m, n, ref taup, nrvt, ref vt);
                bidiagonal.unpackdiagonalsfrombidiagonal(ref a, m, n, ref isupper, ref w, ref e);
                work = new double[m+1];
                blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
                result = bdsvd.bidiagonalsvddecomposition(ref w, e, minmn, isupper, false, ref a, 0, ref u, nru, ref vt, ncvt);
                blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
                return result;
            }
            
            //
            // Simple bidiagonal reduction
            //
            bidiagonal.tobidiagonal(ref a, m, n, ref tauq, ref taup);
            bidiagonal.unpackqfrombidiagonal(ref a, m, n, ref tauq, ncu, ref u);
            bidiagonal.unpackptfrombidiagonal(ref a, m, n, ref taup, nrvt, ref vt);
            bidiagonal.unpackdiagonalsfrombidiagonal(ref a, m, n, ref isupper, ref w, ref e);
            if( additionalmemory<2 | uneeded==0 )
            {
                
                //
                // We cant use additional memory or there is no need in such operations
                //
                result = bdsvd.bidiagonalsvddecomposition(ref w, e, minmn, isupper, false, ref u, nru, ref a, 0, ref vt, ncvt);
            }
            else
            {
                
                //
                // We can use additional memory
                //
                t2 = new double[minmn+1, m+1];
                blas.copyandtranspose(ref u, 1, m, 1, minmn, ref t2, 1, minmn, 1, m);
                result = bdsvd.bidiagonalsvddecomposition(ref w, e, minmn, isupper, false, ref u, 0, ref t2, m, ref vt, ncvt);
                blas.copyandtranspose(ref t2, 1, minmn, 1, m, ref u, 1, m, 1, minmn);
            }
            return result;
        }
    }
}