Free cookie consent management tool by TermsFeed Policy Generator

Ignore:
Timestamp:
10/15/09 13:22:07 (15 years ago)
Author:
gkronber
Message:

Moved ALGLIB code into a separate plugin. #783

Location:
trunk/sources/ALGLIB
Files:
1 added
1 copied

Legend:

Unmodified
Added
Removed

  • trunk/sources/ALGLIB/svd.cs

    r2426 r2430  
    22Copyright (c) 2005-2007, Sergey Bochkanov (ALGLIB project).
    33
    4 Redistribution and use in source and binary forms, with or without
    5 modification, are permitted provided that the following conditions are
    6 met:
    7 
    8 - Redistributions of source code must retain the above copyright
    9   notice, this list of conditions and the following disclaimer.
    10 
    11 - Redistributions in binary form must reproduce the above copyright
    12   notice, this list of conditions and the following disclaimer listed
    13   in this license in the documentation and/or other materials
    14   provided with the distribution.
    15 
    16 - Neither the name of the copyright holders nor the names of its
    17   contributors may be used to endorse or promote products derived from
    18   this software without specific prior written permission.
    19 
    20 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
    21 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
    22 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
    23 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
    24 OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
    25 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
    26 LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
    27 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
    28 THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
    29 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
    30 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
     4>>> SOURCE LICENSE >>>
     5This program is free software; you can redistribute it and/or modify
     6it under the terms of the GNU General Public License as published by
     7the Free Software Foundation (www.fsf.org); either version 2 of the
     8License, or (at your option) any later version.
     9
     10This program is distributed in the hope that it will be useful,
     11but WITHOUT ANY WARRANTY; without even the implied warranty of
     12MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
     13GNU General Public License for more details.
     14
     15A copy of the GNU General Public License is available at
     16http://www.fsf.org/licensing/licenses
     17
     18>>> END OF LICENSE >>>
    3119*************************************************************************/
    3220
    3321using System;
    3422
    35 class svd
     23namespace alglib
    3624{
    37     /*************************************************************************
    38     Singular value decomposition of a rectangular matrix.
    39 
    40     The algorithm calculates the singular value decomposition of a matrix of
    41     size MxN: A = U * S * V^T
    42 
    43     The algorithm finds the singular values and, optionally, matrices U and V^T.
    44     The algorithm can find both first min(M,N) columns of matrix U and rows of
    45     matrix V^T (singular vectors), and matrices U and V^T wholly (of sizes MxM
    46     and NxN respectively).
    47 
    48     Take into account that the subroutine does not return matrix V but V^T.
    49 
    50     Input parameters:
    51         A           -   matrix to be decomposed.
    52                         Array whose indexes range within [0..M-1, 0..N-1].
    53         M           -   number of rows in matrix A.
    54         N           -   number of columns in matrix A.
    55         UNeeded     -   0, 1 or 2. See the description of the parameter U.
    56         VTNeeded    -   0, 1 or 2. See the description of the parameter VT.
    57         AdditionalMemory -
    58                         If the parameter:
    59                          * equals 0, the algorithm doesn’t use additional
    60                            memory (lower requirements, lower performance).
    61                          * equals 1, the algorithm uses additional
    62                            memory of size min(M,N)*min(M,N) of real numbers.
    63                            It often speeds up the algorithm.
    64                          * equals 2, the algorithm uses additional
    65                            memory of size M*min(M,N) of real numbers.
    66                            It allows to get a maximum performance.
    67                         The recommended value of the parameter is 2.
    68 
    69     Output parameters:
    70         W           -   contains singular values in descending order.
    71         U           -   if UNeeded=0, U isn't changed, the left singular vectors
    72                         are not calculated.
    73                         if Uneeded=1, U contains left singular vectors (first
    74                         min(M,N) columns of matrix U). Array whose indexes range
    75                         within [0..M-1, 0..Min(M,N)-1].
    76                         if UNeeded=2, U contains matrix U wholly. Array whose
    77                         indexes range within [0..M-1, 0..M-1].
    78         VT          -   if VTNeeded=0, VT isn’t changed, the right singular vectors
    79                         are not calculated.
    80                         if VTNeeded=1, VT contains right singular vectors (first
    81                         min(M,N) rows of matrix V^T). Array whose indexes range
    82                         within [0..min(M,N)-1, 0..N-1].
    83                         if VTNeeded=2, VT contains matrix V^T wholly. Array whose
    84                         indexes range within [0..N-1, 0..N-1].
    85 
    86       -- ALGLIB --
    87          Copyright 2005 by Bochkanov Sergey
    88     *************************************************************************/
    89     public static bool rmatrixsvd(double[,] a,
    90         int m,
    91         int n,
    92         int uneeded,
    93         int vtneeded,
    94         int additionalmemory,
    95         ref double[] w,
    96         ref double[,] u,
    97         ref double[,] vt)
     25    public class svd
    9826    {
    99         bool result = new bool();
    100         double[] tauq = new double[0];
    101         double[] taup = new double[0];
    102         double[] tau = new double[0];
    103         double[] e = new double[0];
    104         double[] work = new double[0];
    105         double[,] t2 = new double[0,0];
    106         bool isupper = new bool();
    107         int minmn = 0;
    108         int ncu = 0;
    109         int nrvt = 0;
    110         int nru = 0;
    111         int ncvt = 0;
    112         int i = 0;
    113         int j = 0;
    114         int im1 = 0;
    115         double sm = 0;
    116 
    117         a = (double[,])a.Clone();
    118 
    119         result = true;
    120         if( m==0 | n==0 )
     27        /*************************************************************************
     28        Singular value decomposition of a rectangular matrix.
     29
     30        The algorithm calculates the singular value decomposition of a matrix of
     31        size MxN: A = U * S * V^T
     32
     33        The algorithm finds the singular values and, optionally, matrices U and V^T.
     34        The algorithm can find both first min(M,N) columns of matrix U and rows of
     35        matrix V^T (singular vectors), and matrices U and V^T wholly (of sizes MxM
     36        and NxN respectively).
     37
     38        Take into account that the subroutine does not return matrix V but V^T.
     39
     40        Input parameters:
     41            A           -   matrix to be decomposed.
     42                            Array whose indexes range within [0..M-1, 0..N-1].
     43            M           -   number of rows in matrix A.
     44            N           -   number of columns in matrix A.
     45            UNeeded     -   0, 1 or 2. See the description of the parameter U.
     46            VTNeeded    -   0, 1 or 2. See the description of the parameter VT.
     47            AdditionalMemory -
     48                            If the parameter:
     49                             * equals 0, the algorithm doesn’t use additional
     50                               memory (lower requirements, lower performance).
     51                             * equals 1, the algorithm uses additional
     52                               memory of size min(M,N)*min(M,N) of real numbers.
     53                               It often speeds up the algorithm.
     54                             * equals 2, the algorithm uses additional
     55                               memory of size M*min(M,N) of real numbers.
     56                               It allows to get a maximum performance.
     57                            The recommended value of the parameter is 2.
     58
     59        Output parameters:
     60            W           -   contains singular values in descending order.
     61            U           -   if UNeeded=0, U isn't changed, the left singular vectors
     62                            are not calculated.
     63                            if Uneeded=1, U contains left singular vectors (first
     64                            min(M,N) columns of matrix U). Array whose indexes range
     65                            within [0..M-1, 0..Min(M,N)-1].
     66                            if UNeeded=2, U contains matrix U wholly. Array whose
     67                            indexes range within [0..M-1, 0..M-1].
     68            VT          -   if VTNeeded=0, VT isn’t changed, the right singular vectors
     69                            are not calculated.
     70                            if VTNeeded=1, VT contains right singular vectors (first
     71                            min(M,N) rows of matrix V^T). Array whose indexes range
     72                            within [0..min(M,N)-1, 0..N-1].
     73                            if VTNeeded=2, VT contains matrix V^T wholly. Array whose
     74                            indexes range within [0..N-1, 0..N-1].
     75
     76          -- ALGLIB --
     77             Copyright 2005 by Bochkanov Sergey
     78        *************************************************************************/
     79        public static bool rmatrixsvd(double[,] a,
     80            int m,
     81            int n,
     82            int uneeded,
     83            int vtneeded,
     84            int additionalmemory,
     85            ref double[] w,
     86            ref double[,] u,
     87            ref double[,] vt)
    12188        {
    122             return result;
    123         }
    124         System.Diagnostics.Debug.Assert(uneeded>=0 & uneeded<=2, "SVDDecomposition: wrong parameters!");
    125         System.Diagnostics.Debug.Assert(vtneeded>=0 & vtneeded<=2, "SVDDecomposition: wrong parameters!");
    126         System.Diagnostics.Debug.Assert(additionalmemory>=0 & additionalmemory<=2, "SVDDecomposition: wrong parameters!");
    127        
    128         //
    129         // initialize
    130         //
    131         minmn = Math.Min(m, n);
    132         w = new double[minmn+1];
    133         ncu = 0;
    134         nru = 0;
    135         if( uneeded==1 )
    136         {
    137             nru = m;
    138             ncu = minmn;
    139             u = new double[nru-1+1, ncu-1+1];
    140         }
    141         if( uneeded==2 )
    142         {
    143             nru = m;
    144             ncu = m;
    145             u = new double[nru-1+1, ncu-1+1];
    146         }
    147         nrvt = 0;
    148         ncvt = 0;
    149         if( vtneeded==1 )
    150         {
    151             nrvt = minmn;
    152             ncvt = n;
    153             vt = new double[nrvt-1+1, ncvt-1+1];
    154         }
    155         if( vtneeded==2 )
    156         {
    157             nrvt = n;
    158             ncvt = n;
    159             vt = new double[nrvt-1+1, ncvt-1+1];
    160         }
    161        
    162         //
    163         // M much larger than N
    164         // Use bidiagonal reduction with QR-decomposition
    165         //
    166         if( m>1.6*n )
    167         {
    168             if( uneeded==0 )
    169             {
    170                
    171                 //
    172                 // No left singular vectors to be computed
    173                 //
    174                 qr.rmatrixqr(ref a, m, n, ref tau);
    175                 for(i=0; i<=n-1; i++)
    176                 {
    177                     for(j=0; j<=i-1; j++)
    178                     {
    179                         a[i,j] = 0;
    180                     }
    181                 }
    182                 bidiagonal.rmatrixbd(ref a, n, n, ref tauq, ref taup);
    183                 bidiagonal.rmatrixbdunpackpt(ref a, n, n, ref taup, nrvt, ref vt);
    184                 bidiagonal.rmatrixbdunpackdiagonals(ref a, n, n, ref isupper, ref w, ref e);
    185                 result = bdsvd.rmatrixbdsvd(ref w, e, n, isupper, false, ref u, 0, ref a, 0, ref vt, ncvt);
     89            bool result = new bool();
     90            double[] tauq = new double[0];
     91            double[] taup = new double[0];
     92            double[] tau = new double[0];
     93            double[] e = new double[0];
     94            double[] work = new double[0];
     95            double[,] t2 = new double[0,0];
     96            bool isupper = new bool();
     97            int minmn = 0;
     98            int ncu = 0;
     99            int nrvt = 0;
     100            int nru = 0;
     101            int ncvt = 0;
     102            int i = 0;
     103            int j = 0;
     104            int im1 = 0;
     105            double sm = 0;
     106
     107            a = (double[,])a.Clone();
     108
     109            result = true;
     110            if( m==0 | n==0 )
     111            {
    186112                return result;
    187113            }
    188             else
    189             {
    190                
    191                 //
    192                 // Left singular vectors (may be full matrix U) to be computed
    193                 //
    194                 qr.rmatrixqr(ref a, m, n, ref tau);
    195                 qr.rmatrixqrunpackq(ref a, m, n, ref tau, ncu, ref u);
    196                 for(i=0; i<=n-1; i++)
    197                 {
    198                     for(j=0; j<=i-1; j++)
    199                     {
    200                         a[i,j] = 0;
    201                     }
    202                 }
    203                 bidiagonal.rmatrixbd(ref a, n, n, ref tauq, ref taup);
    204                 bidiagonal.rmatrixbdunpackpt(ref a, n, n, ref taup, nrvt, ref vt);
    205                 bidiagonal.rmatrixbdunpackdiagonals(ref a, n, n, ref isupper, ref w, ref e);
    206                 if( additionalmemory<1 )
    207                 {
    208                    
    209                     //
    210                     // No additional memory can be used
    211                     //
    212                     bidiagonal.rmatrixbdmultiplybyq(ref a, n, n, ref tauq, ref u, m, n, true, false);
    213                     result = bdsvd.rmatrixbdsvd(ref w, e, n, isupper, false, ref u, m, ref a, 0, ref vt, ncvt);
     114            System.Diagnostics.Debug.Assert(uneeded>=0 & uneeded<=2, "SVDDecomposition: wrong parameters!");
     115            System.Diagnostics.Debug.Assert(vtneeded>=0 & vtneeded<=2, "SVDDecomposition: wrong parameters!");
     116            System.Diagnostics.Debug.Assert(additionalmemory>=0 & additionalmemory<=2, "SVDDecomposition: wrong parameters!");
     117           
     118            //
     119            // initialize
     120            //
     121            minmn = Math.Min(m, n);
     122            w = new double[minmn+1];
     123            ncu = 0;
     124            nru = 0;
     125            if( uneeded==1 )
     126            {
     127                nru = m;
     128                ncu = minmn;
     129                u = new double[nru-1+1, ncu-1+1];
     130            }
     131            if( uneeded==2 )
     132            {
     133                nru = m;
     134                ncu = m;
     135                u = new double[nru-1+1, ncu-1+1];
     136            }
     137            nrvt = 0;
     138            ncvt = 0;
     139            if( vtneeded==1 )
     140            {
     141                nrvt = minmn;
     142                ncvt = n;
     143                vt = new double[nrvt-1+1, ncvt-1+1];
     144            }
     145            if( vtneeded==2 )
     146            {
     147                nrvt = n;
     148                ncvt = n;
     149                vt = new double[nrvt-1+1, ncvt-1+1];
     150            }
     151           
     152            //
     153            // M much larger than N
     154            // Use bidiagonal reduction with QR-decomposition
     155            //
     156            if( m>1.6*n )
     157            {
     158                if( uneeded==0 )
     159                {
     160                   
     161                    //
     162                    // No left singular vectors to be computed
     163                    //
     164                    qr.rmatrixqr(ref a, m, n, ref tau);
     165                    for(i=0; i<=n-1; i++)
     166                    {
     167                        for(j=0; j<=i-1; j++)
     168                        {
     169                            a[i,j] = 0;
     170                        }
     171                    }
     172                    bidiagonal.rmatrixbd(ref a, n, n, ref tauq, ref taup);
     173                    bidiagonal.rmatrixbdunpackpt(ref a, n, n, ref taup, nrvt, ref vt);
     174                    bidiagonal.rmatrixbdunpackdiagonals(ref a, n, n, ref isupper, ref w, ref e);
     175                    result = bdsvd.rmatrixbdsvd(ref w, e, n, isupper, false, ref u, 0, ref a, 0, ref vt, ncvt);
     176                    return result;
    214177                }
    215178                else
    … …  
    217180                   
    218181                    //
    219                     // Large U. Transforming intermediate matrix T2
    220                     //
     182                    // Left singular vectors (may be full matrix U) to be computed
     183                    //
     184                    qr.rmatrixqr(ref a, m, n, ref tau);
     185                    qr.rmatrixqrunpackq(ref a, m, n, ref tau, ncu, ref u);
     186                    for(i=0; i<=n-1; i++)
     187                    {
     188                        for(j=0; j<=i-1; j++)
     189                        {
     190                            a[i,j] = 0;
     191                        }
     192                    }
     193                    bidiagonal.rmatrixbd(ref a, n, n, ref tauq, ref taup);
     194                    bidiagonal.rmatrixbdunpackpt(ref a, n, n, ref taup, nrvt, ref vt);
     195                    bidiagonal.rmatrixbdunpackdiagonals(ref a, n, n, ref isupper, ref w, ref e);
     196                    if( additionalmemory<1 )
     197                    {
     198                       
     199                        //
     200                        // No additional memory can be used
     201                        //
     202                        bidiagonal.rmatrixbdmultiplybyq(ref a, n, n, ref tauq, ref u, m, n, true, false);
     203                        result = bdsvd.rmatrixbdsvd(ref w, e, n, isupper, false, ref u, m, ref a, 0, ref vt, ncvt);
     204                    }
     205                    else
     206                    {
     207                       
     208                        //
     209                        // Large U. Transforming intermediate matrix T2
     210                        //
     211                        work = new double[Math.Max(m, n)+1];
     212                        bidiagonal.rmatrixbdunpackq(ref a, n, n, ref tauq, n, ref t2);
     213                        blas.copymatrix(ref u, 0, m-1, 0, n-1, ref a, 0, m-1, 0, n-1);
     214                        blas.inplacetranspose(ref t2, 0, n-1, 0, n-1, ref work);
     215                        result = bdsvd.rmatrixbdsvd(ref w, e, n, isupper, false, ref u, 0, ref t2, n, ref vt, ncvt);
     216                        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);
     217                    }
     218                    return result;
     219                }
     220            }
     221           
     222            //
     223            // N much larger than M
     224            // Use bidiagonal reduction with LQ-decomposition
     225            //
     226            if( n>1.6*m )
     227            {
     228                if( vtneeded==0 )
     229                {
     230                   
     231                    //
     232                    // No right singular vectors to be computed
     233                    //
     234                    lq.rmatrixlq(ref a, m, n, ref tau);
     235                    for(i=0; i<=m-1; i++)
     236                    {
     237                        for(j=i+1; j<=m-1; j++)
     238                        {
     239                            a[i,j] = 0;
     240                        }
     241                    }
     242                    bidiagonal.rmatrixbd(ref a, m, m, ref tauq, ref taup);
     243                    bidiagonal.rmatrixbdunpackq(ref a, m, m, ref tauq, ncu, ref u);
     244                    bidiagonal.rmatrixbdunpackdiagonals(ref a, m, m, ref isupper, ref w, ref e);
     245                    work = new double[m+1];
     246                    blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
     247                    result = bdsvd.rmatrixbdsvd(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref vt, 0);
     248                    blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
     249                    return result;
     250                }
     251                else
     252                {
     253                   
     254                    //
     255                    // Right singular vectors (may be full matrix VT) to be computed
     256                    //
     257                    lq.rmatrixlq(ref a, m, n, ref tau);
     258                    lq.rmatrixlqunpackq(ref a, m, n, ref tau, nrvt, ref vt);
     259                    for(i=0; i<=m-1; i++)
     260                    {
     261                        for(j=i+1; j<=m-1; j++)
     262                        {
     263                            a[i,j] = 0;
     264                        }
     265                    }
     266                    bidiagonal.rmatrixbd(ref a, m, m, ref tauq, ref taup);
     267                    bidiagonal.rmatrixbdunpackq(ref a, m, m, ref tauq, ncu, ref u);
     268                    bidiagonal.rmatrixbdunpackdiagonals(ref a, m, m, ref isupper, ref w, ref e);
    221269                    work = new double[Math.Max(m, n)+1];
    222                     bidiagonal.rmatrixbdunpackq(ref a, n, n, ref tauq, n, ref t2);
    223                     blas.copymatrix(ref u, 0, m-1, 0, n-1, ref a, 0, m-1, 0, n-1);
    224                     blas.inplacetranspose(ref t2, 0, n-1, 0, n-1, ref work);
    225                     result = bdsvd.rmatrixbdsvd(ref w, e, n, isupper, false, ref u, 0, ref t2, n, ref vt, ncvt);
    226                     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);
    227                 }
    228                 return result;
    229             }
    230         }
    231        
    232         //
    233         // N much larger than M
    234         // Use bidiagonal reduction with LQ-decomposition
    235         //
    236         if( n>1.6*m )
    237         {
    238             if( vtneeded==0 )
    239             {
    240                
    241                 //
    242                 // No right singular vectors to be computed
    243                 //
    244                 lq.rmatrixlq(ref a, m, n, ref tau);
    245                 for(i=0; i<=m-1; i++)
    246                 {
    247                     for(j=i+1; j<=m-1; j++)
    248                     {
    249                         a[i,j] = 0;
    250                     }
    251                 }
    252                 bidiagonal.rmatrixbd(ref a, m, m, ref tauq, ref taup);
    253                 bidiagonal.rmatrixbdunpackq(ref a, m, m, ref tauq, ncu, ref u);
    254                 bidiagonal.rmatrixbdunpackdiagonals(ref a, m, m, ref isupper, ref w, ref e);
     270                    blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
     271                    if( additionalmemory<1 )
     272                    {
     273                       
     274                        //
     275                        // No additional memory available
     276                        //
     277                        bidiagonal.rmatrixbdmultiplybyp(ref a, m, m, ref taup, ref vt, m, n, false, true);
     278                        result = bdsvd.rmatrixbdsvd(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref vt, n);
     279                    }
     280                    else
     281                    {
     282                       
     283                        //
     284                        // Large VT. Transforming intermediate matrix T2
     285                        //
     286                        bidiagonal.rmatrixbdunpackpt(ref a, m, m, ref taup, m, ref t2);
     287                        result = bdsvd.rmatrixbdsvd(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref t2, m);
     288                        blas.copymatrix(ref vt, 0, m-1, 0, n-1, ref a, 0, m-1, 0, n-1);
     289                        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);
     290                    }
     291                    blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
     292                    return result;
     293                }
     294            }
     295           
     296            //
     297            // M<=N
     298            // We can use inplace transposition of U to get rid of columnwise operations
     299            //
     300            if( m<=n )
     301            {
     302                bidiagonal.rmatrixbd(ref a, m, n, ref tauq, ref taup);
     303                bidiagonal.rmatrixbdunpackq(ref a, m, n, ref tauq, ncu, ref u);
     304                bidiagonal.rmatrixbdunpackpt(ref a, m, n, ref taup, nrvt, ref vt);
     305                bidiagonal.rmatrixbdunpackdiagonals(ref a, m, n, ref isupper, ref w, ref e);
    255306                work = new double[m+1];
    256307                blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
    257                 result = bdsvd.rmatrixbdsvd(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref vt, 0);
     308                result = bdsvd.rmatrixbdsvd(ref w, e, minmn, isupper, false, ref a, 0, ref u, nru, ref vt, ncvt);
    258309                blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
    259310                return result;
    260311            }
    261             else
    262             {
    263                
    264                 //
    265                 // Right singular vectors (may be full matrix VT) to be computed
    266                 //
    267                 lq.rmatrixlq(ref a, m, n, ref tau);
    268                 lq.rmatrixlqunpackq(ref a, m, n, ref tau, nrvt, ref vt);
    269                 for(i=0; i<=m-1; i++)
    270                 {
    271                     for(j=i+1; j<=m-1; j++)
    272                     {
    273                         a[i,j] = 0;
    274                     }
    275                 }
    276                 bidiagonal.rmatrixbd(ref a, m, m, ref tauq, ref taup);
    277                 bidiagonal.rmatrixbdunpackq(ref a, m, m, ref tauq, ncu, ref u);
    278                 bidiagonal.rmatrixbdunpackdiagonals(ref a, m, m, ref isupper, ref w, ref e);
    279                 work = new double[Math.Max(m, n)+1];
    280                 blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
    281                 if( additionalmemory<1 )
    282                 {
    283                    
    284                     //
    285                     // No additional memory available
    286                     //
    287                     bidiagonal.rmatrixbdmultiplybyp(ref a, m, m, ref taup, ref vt, m, n, false, true);
    288                     result = bdsvd.rmatrixbdsvd(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref vt, n);
    289                 }
    290                 else
    291                 {
    292                    
    293                     //
    294                     // Large VT. Transforming intermediate matrix T2
    295                     //
    296                     bidiagonal.rmatrixbdunpackpt(ref a, m, m, ref taup, m, ref t2);
    297                     result = bdsvd.rmatrixbdsvd(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref t2, m);
    298                     blas.copymatrix(ref vt, 0, m-1, 0, n-1, ref a, 0, m-1, 0, n-1);
    299                     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);
    300                 }
    301                 blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
    302                 return result;
    303             }
    304         }
    305        
    306         //
    307         // M<=N
    308         // We can use inplace transposition of U to get rid of columnwise operations
    309         //
    310         if( m<=n )
    311         {
     312           
     313            //
     314            // Simple bidiagonal reduction
     315            //
    312316            bidiagonal.rmatrixbd(ref a, m, n, ref tauq, ref taup);
    313317            bidiagonal.rmatrixbdunpackq(ref a, m, n, ref tauq, ncu, ref u);
    314318            bidiagonal.rmatrixbdunpackpt(ref a, m, n, ref taup, nrvt, ref vt);
    315319            bidiagonal.rmatrixbdunpackdiagonals(ref a, m, n, ref isupper, ref w, ref e);
    316             work = new double[m+1];
    317             blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
    318             result = bdsvd.rmatrixbdsvd(ref w, e, minmn, isupper, false, ref a, 0, ref u, nru, ref vt, ncvt);
    319             blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
     320            if( additionalmemory<2 | uneeded==0 )
     321            {
     322               
     323                //
     324                // We cant use additional memory or there is no need in such operations
     325                //
     326                result = bdsvd.rmatrixbdsvd(ref w, e, minmn, isupper, false, ref u, nru, ref a, 0, ref vt, ncvt);
     327            }
     328            else
     329            {
     330               
     331                //
     332                // We can use additional memory
     333                //
     334                t2 = new double[minmn-1+1, m-1+1];
     335                blas.copyandtranspose(ref u, 0, m-1, 0, minmn-1, ref t2, 0, minmn-1, 0, m-1);
     336                result = bdsvd.rmatrixbdsvd(ref w, e, minmn, isupper, false, ref u, 0, ref t2, m, ref vt, ncvt);
     337                blas.copyandtranspose(ref t2, 0, minmn-1, 0, m-1, ref u, 0, m-1, 0, minmn-1);
     338            }
    320339            return result;
    321340        }
    322        
    323         //
    324         // Simple bidiagonal reduction
    325         //
    326         bidiagonal.rmatrixbd(ref a, m, n, ref tauq, ref taup);
    327         bidiagonal.rmatrixbdunpackq(ref a, m, n, ref tauq, ncu, ref u);
    328         bidiagonal.rmatrixbdunpackpt(ref a, m, n, ref taup, nrvt, ref vt);
    329         bidiagonal.rmatrixbdunpackdiagonals(ref a, m, n, ref isupper, ref w, ref e);
    330         if( additionalmemory<2 | uneeded==0 )
     341
     342
     343        /*************************************************************************
     344        Obsolete 1-based subroutine.
     345        See RMatrixSVD for 0-based replacement.
     346        *************************************************************************/
     347        public static bool svddecomposition(double[,] a,
     348            int m,
     349            int n,
     350            int uneeded,
     351            int vtneeded,
     352            int additionalmemory,
     353            ref double[] w,
     354            ref double[,] u,
     355            ref double[,] vt)
    331356        {
    332            
    333             //
    334             // We cant use additional memory or there is no need in such operations
    335             //
    336             result = bdsvd.rmatrixbdsvd(ref w, e, minmn, isupper, false, ref u, nru, ref a, 0, ref vt, ncvt);
    337         }
    338         else
    339         {
    340            
    341             //
    342             // We can use additional memory
    343             //
    344             t2 = new double[minmn-1+1, m-1+1];
    345             blas.copyandtranspose(ref u, 0, m-1, 0, minmn-1, ref t2, 0, minmn-1, 0, m-1);
    346             result = bdsvd.rmatrixbdsvd(ref w, e, minmn, isupper, false, ref u, 0, ref t2, m, ref vt, ncvt);
    347             blas.copyandtranspose(ref t2, 0, minmn-1, 0, m-1, ref u, 0, m-1, 0, minmn-1);
    348         }
    349         return result;
    350     }
    351 
    352 
    353     /*************************************************************************
    354     Obsolete 1-based subroutine.
    355     See RMatrixSVD for 0-based replacement.
    356     *************************************************************************/
    357     public static bool svddecomposition(double[,] a,
    358         int m,
    359         int n,
    360         int uneeded,
    361         int vtneeded,
    362         int additionalmemory,
    363         ref double[] w,
    364         ref double[,] u,
    365         ref double[,] vt)
    366     {
    367         bool result = new bool();
    368         double[] tauq = new double[0];
    369         double[] taup = new double[0];
    370         double[] tau = new double[0];
    371         double[] e = new double[0];
    372         double[] work = new double[0];
    373         double[,] t2 = new double[0,0];
    374         bool isupper = new bool();
    375         int minmn = 0;
    376         int ncu = 0;
    377         int nrvt = 0;
    378         int nru = 0;
    379         int ncvt = 0;
    380         int i = 0;
    381         int j = 0;
    382         int im1 = 0;
    383         double sm = 0;
    384 
    385         a = (double[,])a.Clone();
    386 
    387         result = true;
    388         if( m==0 | n==0 )
    389         {
    390             return result;
    391         }
    392         System.Diagnostics.Debug.Assert(uneeded>=0 & uneeded<=2, "SVDDecomposition: wrong parameters!");
    393         System.Diagnostics.Debug.Assert(vtneeded>=0 & vtneeded<=2, "SVDDecomposition: wrong parameters!");
    394         System.Diagnostics.Debug.Assert(additionalmemory>=0 & additionalmemory<=2, "SVDDecomposition: wrong parameters!");
    395        
    396         //
    397         // initialize
    398         //
    399         minmn = Math.Min(m, n);
    400         w = new double[minmn+1];
    401         ncu = 0;
    402         nru = 0;
    403         if( uneeded==1 )
    404         {
    405             nru = m;
    406             ncu = minmn;
    407             u = new double[nru+1, ncu+1];
    408         }
    409         if( uneeded==2 )
    410         {
    411             nru = m;
    412             ncu = m;
    413             u = new double[nru+1, ncu+1];
    414         }
    415         nrvt = 0;
    416         ncvt = 0;
    417         if( vtneeded==1 )
    418         {
    419             nrvt = minmn;
    420             ncvt = n;
    421             vt = new double[nrvt+1, ncvt+1];
    422         }
    423         if( vtneeded==2 )
    424         {
    425             nrvt = n;
    426             ncvt = n;
    427             vt = new double[nrvt+1, ncvt+1];
    428         }
    429        
    430         //
    431         // M much larger than N
    432         // Use bidiagonal reduction with QR-decomposition
    433         //
    434         if( m>1.6*n )
    435         {
    436             if( uneeded==0 )
    437             {
    438                
    439                 //
    440                 // No left singular vectors to be computed
    441                 //
    442                 qr.qrdecomposition(ref a, m, n, ref tau);
    443                 for(i=2; i<=n; i++)
    444                 {
    445                     for(j=1; j<=i-1; j++)
    446                     {
    447                         a[i,j] = 0;
    448                     }
    449                 }
    450                 bidiagonal.tobidiagonal(ref a, n, n, ref tauq, ref taup);
    451                 bidiagonal.unpackptfrombidiagonal(ref a, n, n, ref taup, nrvt, ref vt);
    452                 bidiagonal.unpackdiagonalsfrombidiagonal(ref a, n, n, ref isupper, ref w, ref e);
    453                 result = bdsvd.bidiagonalsvddecomposition(ref w, e, n, isupper, false, ref u, 0, ref a, 0, ref vt, ncvt);
     357            bool result = new bool();
     358            double[] tauq = new double[0];
     359            double[] taup = new double[0];
     360            double[] tau = new double[0];
     361            double[] e = new double[0];
     362            double[] work = new double[0];
     363            double[,] t2 = new double[0,0];
     364            bool isupper = new bool();
     365            int minmn = 0;
     366            int ncu = 0;
     367            int nrvt = 0;
     368            int nru = 0;
     369            int ncvt = 0;
     370            int i = 0;
     371            int j = 0;
     372            int im1 = 0;
     373            double sm = 0;
     374
     375            a = (double[,])a.Clone();
     376
     377            result = true;
     378            if( m==0 | n==0 )
     379            {
    454380                return result;
    455381            }
    456             else
    457             {
    458                
    459                 //
    460                 // Left singular vectors (may be full matrix U) to be computed
    461                 //
    462                 qr.qrdecomposition(ref a, m, n, ref tau);
    463                 qr.unpackqfromqr(ref a, m, n, ref tau, ncu, ref u);
    464                 for(i=2; i<=n; i++)
    465                 {
    466                     for(j=1; j<=i-1; j++)
    467                     {
    468                         a[i,j] = 0;
    469                     }
    470                 }
    471                 bidiagonal.tobidiagonal(ref a, n, n, ref tauq, ref taup);
    472                 bidiagonal.unpackptfrombidiagonal(ref a, n, n, ref taup, nrvt, ref vt);
    473                 bidiagonal.unpackdiagonalsfrombidiagonal(ref a, n, n, ref isupper, ref w, ref e);
    474                 if( additionalmemory<1 )
    475                 {
    476                    
    477                     //
    478                     // No additional memory can be used
    479                     //
    480                     bidiagonal.multiplybyqfrombidiagonal(ref a, n, n, ref tauq, ref u, m, n, true, false);
    481                     result = bdsvd.bidiagonalsvddecomposition(ref w, e, n, isupper, false, ref u, m, ref a, 0, ref vt, ncvt);
     382            System.Diagnostics.Debug.Assert(uneeded>=0 & uneeded<=2, "SVDDecomposition: wrong parameters!");
     383            System.Diagnostics.Debug.Assert(vtneeded>=0 & vtneeded<=2, "SVDDecomposition: wrong parameters!");
     384            System.Diagnostics.Debug.Assert(additionalmemory>=0 & additionalmemory<=2, "SVDDecomposition: wrong parameters!");
     385           
     386            //
     387            // initialize
     388            //
     389            minmn = Math.Min(m, n);
     390            w = new double[minmn+1];
     391            ncu = 0;
     392            nru = 0;
     393            if( uneeded==1 )
     394            {
     395                nru = m;
     396                ncu = minmn;
     397                u = new double[nru+1, ncu+1];
     398            }
     399            if( uneeded==2 )
     400            {
     401                nru = m;
     402                ncu = m;
     403                u = new double[nru+1, ncu+1];
     404            }
     405            nrvt = 0;
     406            ncvt = 0;
     407            if( vtneeded==1 )
     408            {
     409                nrvt = minmn;
     410                ncvt = n;
     411                vt = new double[nrvt+1, ncvt+1];
     412            }
     413            if( vtneeded==2 )
     414            {
     415                nrvt = n;
     416                ncvt = n;
     417                vt = new double[nrvt+1, ncvt+1];
     418            }
     419           
     420            //
     421            // M much larger than N
     422            // Use bidiagonal reduction with QR-decomposition
     423            //
     424            if( m>1.6*n )
     425            {
     426                if( uneeded==0 )
     427                {
     428                   
     429                    //
     430                    // No left singular vectors to be computed
     431                    //
     432                    qr.qrdecomposition(ref a, m, n, ref tau);
     433                    for(i=2; i<=n; i++)
     434                    {
     435                        for(j=1; j<=i-1; j++)
     436                        {
     437                            a[i,j] = 0;
     438                        }
     439                    }
     440                    bidiagonal.tobidiagonal(ref a, n, n, ref tauq, ref taup);
     441                    bidiagonal.unpackptfrombidiagonal(ref a, n, n, ref taup, nrvt, ref vt);
     442                    bidiagonal.unpackdiagonalsfrombidiagonal(ref a, n, n, ref isupper, ref w, ref e);
     443                    result = bdsvd.bidiagonalsvddecomposition(ref w, e, n, isupper, false, ref u, 0, ref a, 0, ref vt, ncvt);
     444                    return result;
    482445                }
    483446                else
    … …  
    485448                   
    486449                    //
    487                     // Large U. Transforming intermediate matrix T2
    488                     //
     450                    // Left singular vectors (may be full matrix U) to be computed
     451                    //
     452                    qr.qrdecomposition(ref a, m, n, ref tau);
     453                    qr.unpackqfromqr(ref a, m, n, ref tau, ncu, ref u);
     454                    for(i=2; i<=n; i++)
     455                    {
     456                        for(j=1; j<=i-1; j++)
     457                        {
     458                            a[i,j] = 0;
     459                        }
     460                    }
     461                    bidiagonal.tobidiagonal(ref a, n, n, ref tauq, ref taup);
     462                    bidiagonal.unpackptfrombidiagonal(ref a, n, n, ref taup, nrvt, ref vt);
     463                    bidiagonal.unpackdiagonalsfrombidiagonal(ref a, n, n, ref isupper, ref w, ref e);
     464                    if( additionalmemory<1 )
     465                    {
     466                       
     467                        //
     468                        // No additional memory can be used
     469                        //
     470                        bidiagonal.multiplybyqfrombidiagonal(ref a, n, n, ref tauq, ref u, m, n, true, false);
     471                        result = bdsvd.bidiagonalsvddecomposition(ref w, e, n, isupper, false, ref u, m, ref a, 0, ref vt, ncvt);
     472                    }
     473                    else
     474                    {
     475                       
     476                        //
     477                        // Large U. Transforming intermediate matrix T2
     478                        //
     479                        work = new double[Math.Max(m, n)+1];
     480                        bidiagonal.unpackqfrombidiagonal(ref a, n, n, ref tauq, n, ref t2);
     481                        blas.copymatrix(ref u, 1, m, 1, n, ref a, 1, m, 1, n);
     482                        blas.inplacetranspose(ref t2, 1, n, 1, n, ref work);
     483                        result = bdsvd.bidiagonalsvddecomposition(ref w, e, n, isupper, false, ref u, 0, ref t2, n, ref vt, ncvt);
     484                        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);
     485                    }
     486                    return result;
     487                }
     488            }
     489           
     490            //
     491            // N much larger than M
     492            // Use bidiagonal reduction with LQ-decomposition
     493            //
     494            if( n>1.6*m )
     495            {
     496                if( vtneeded==0 )
     497                {
     498                   
     499                    //
     500                    // No right singular vectors to be computed
     501                    //
     502                    lq.lqdecomposition(ref a, m, n, ref tau);
     503                    for(i=1; i<=m-1; i++)
     504                    {
     505                        for(j=i+1; j<=m; j++)
     506                        {
     507                            a[i,j] = 0;
     508                        }
     509                    }
     510                    bidiagonal.tobidiagonal(ref a, m, m, ref tauq, ref taup);
     511                    bidiagonal.unpackqfrombidiagonal(ref a, m, m, ref tauq, ncu, ref u);
     512                    bidiagonal.unpackdiagonalsfrombidiagonal(ref a, m, m, ref isupper, ref w, ref e);
     513                    work = new double[m+1];
     514                    blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
     515                    result = bdsvd.bidiagonalsvddecomposition(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref vt, 0);
     516                    blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
     517                    return result;
     518                }
     519                else
     520                {
     521                   
     522                    //
     523                    // Right singular vectors (may be full matrix VT) to be computed
     524                    //
     525                    lq.lqdecomposition(ref a, m, n, ref tau);
     526                    lq.unpackqfromlq(ref a, m, n, ref tau, nrvt, ref vt);
     527                    for(i=1; i<=m-1; i++)
     528                    {
     529                        for(j=i+1; j<=m; j++)
     530                        {
     531                            a[i,j] = 0;
     532                        }
     533                    }
     534                    bidiagonal.tobidiagonal(ref a, m, m, ref tauq, ref taup);
     535                    bidiagonal.unpackqfrombidiagonal(ref a, m, m, ref tauq, ncu, ref u);
     536                    bidiagonal.unpackdiagonalsfrombidiagonal(ref a, m, m, ref isupper, ref w, ref e);
    489537                    work = new double[Math.Max(m, n)+1];
    490                     bidiagonal.unpackqfrombidiagonal(ref a, n, n, ref tauq, n, ref t2);
    491                     blas.copymatrix(ref u, 1, m, 1, n, ref a, 1, m, 1, n);
    492                     blas.inplacetranspose(ref t2, 1, n, 1, n, ref work);
    493                     result = bdsvd.bidiagonalsvddecomposition(ref w, e, n, isupper, false, ref u, 0, ref t2, n, ref vt, ncvt);
    494                     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);
    495                 }
    496                 return result;
    497             }
    498         }
    499        
    500         //
    501         // N much larger than M
    502         // Use bidiagonal reduction with LQ-decomposition
    503         //
    504         if( n>1.6*m )
    505         {
    506             if( vtneeded==0 )
    507             {
    508                
    509                 //
    510                 // No right singular vectors to be computed
    511                 //
    512                 lq.lqdecomposition(ref a, m, n, ref tau);
    513                 for(i=1; i<=m-1; i++)
    514                 {
    515                     for(j=i+1; j<=m; j++)
    516                     {
    517                         a[i,j] = 0;
    518                     }
    519                 }
    520                 bidiagonal.tobidiagonal(ref a, m, m, ref tauq, ref taup);
    521                 bidiagonal.unpackqfrombidiagonal(ref a, m, m, ref tauq, ncu, ref u);
    522                 bidiagonal.unpackdiagonalsfrombidiagonal(ref a, m, m, ref isupper, ref w, ref e);
     538                    blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
     539                    if( additionalmemory<1 )
     540                    {
     541                       
     542                        //
     543                        // No additional memory available
     544                        //
     545                        bidiagonal.multiplybypfrombidiagonal(ref a, m, m, ref taup, ref vt, m, n, false, true);
     546                        result = bdsvd.bidiagonalsvddecomposition(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref vt, n);
     547                    }
     548                    else
     549                    {
     550                       
     551                        //
     552                        // Large VT. Transforming intermediate matrix T2
     553                        //
     554                        bidiagonal.unpackptfrombidiagonal(ref a, m, m, ref taup, m, ref t2);
     555                        result = bdsvd.bidiagonalsvddecomposition(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref t2, m);
     556                        blas.copymatrix(ref vt, 1, m, 1, n, ref a, 1, m, 1, n);
     557                        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);
     558                    }
     559                    blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
     560                    return result;
     561                }
     562            }
     563           
     564            //
     565            // M<=N
     566            // We can use inplace transposition of U to get rid of columnwise operations
     567            //
     568            if( m<=n )
     569            {
     570                bidiagonal.tobidiagonal(ref a, m, n, ref tauq, ref taup);
     571                bidiagonal.unpackqfrombidiagonal(ref a, m, n, ref tauq, ncu, ref u);
     572                bidiagonal.unpackptfrombidiagonal(ref a, m, n, ref taup, nrvt, ref vt);
     573                bidiagonal.unpackdiagonalsfrombidiagonal(ref a, m, n, ref isupper, ref w, ref e);
    523574                work = new double[m+1];
    524575                blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
    525                 result = bdsvd.bidiagonalsvddecomposition(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref vt, 0);
     576                result = bdsvd.bidiagonalsvddecomposition(ref w, e, minmn, isupper, false, ref a, 0, ref u, nru, ref vt, ncvt);
    526577                blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
    527578                return result;
    528579            }
    529             else
    530             {
    531                
    532                 //
    533                 // Right singular vectors (may be full matrix VT) to be computed
    534                 //
    535                 lq.lqdecomposition(ref a, m, n, ref tau);
    536                 lq.unpackqfromlq(ref a, m, n, ref tau, nrvt, ref vt);
    537                 for(i=1; i<=m-1; i++)
    538                 {
    539                     for(j=i+1; j<=m; j++)
    540                     {
    541                         a[i,j] = 0;
    542                     }
    543                 }
    544                 bidiagonal.tobidiagonal(ref a, m, m, ref tauq, ref taup);
    545                 bidiagonal.unpackqfrombidiagonal(ref a, m, m, ref tauq, ncu, ref u);
    546                 bidiagonal.unpackdiagonalsfrombidiagonal(ref a, m, m, ref isupper, ref w, ref e);
    547                 work = new double[Math.Max(m, n)+1];
    548                 blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
    549                 if( additionalmemory<1 )
    550                 {
    551                    
    552                     //
    553                     // No additional memory available
    554                     //
    555                     bidiagonal.multiplybypfrombidiagonal(ref a, m, m, ref taup, ref vt, m, n, false, true);
    556                     result = bdsvd.bidiagonalsvddecomposition(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref vt, n);
    557                 }
    558                 else
    559                 {
    560                    
    561                     //
    562                     // Large VT. Transforming intermediate matrix T2
    563                     //
    564                     bidiagonal.unpackptfrombidiagonal(ref a, m, m, ref taup, m, ref t2);
    565                     result = bdsvd.bidiagonalsvddecomposition(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref t2, m);
    566                     blas.copymatrix(ref vt, 1, m, 1, n, ref a, 1, m, 1, n);
    567                     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);
    568                 }
    569                 blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
    570                 return result;
    571             }
    572         }
    573        
    574         //
    575         // M<=N
    576         // We can use inplace transposition of U to get rid of columnwise operations
    577         //
    578         if( m<=n )
    579         {
     580           
     581            //
     582            // Simple bidiagonal reduction
     583            //
    580584            bidiagonal.tobidiagonal(ref a, m, n, ref tauq, ref taup);
    581585            bidiagonal.unpackqfrombidiagonal(ref a, m, n, ref tauq, ncu, ref u);
    582586            bidiagonal.unpackptfrombidiagonal(ref a, m, n, ref taup, nrvt, ref vt);
    583587            bidiagonal.unpackdiagonalsfrombidiagonal(ref a, m, n, ref isupper, ref w, ref e);
    584             work = new double[m+1];
    585             blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
    586             result = bdsvd.bidiagonalsvddecomposition(ref w, e, minmn, isupper, false, ref a, 0, ref u, nru, ref vt, ncvt);
    587             blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
     588            if( additionalmemory<2 | uneeded==0 )
     589            {
     590               
     591                //
     592                // We cant use additional memory or there is no need in such operations
     593                //
     594                result = bdsvd.bidiagonalsvddecomposition(ref w, e, minmn, isupper, false, ref u, nru, ref a, 0, ref vt, ncvt);
     595            }
     596            else
     597            {
     598               
     599                //
     600                // We can use additional memory
     601                //
     602                t2 = new double[minmn+1, m+1];
     603                blas.copyandtranspose(ref u, 1, m, 1, minmn, ref t2, 1, minmn, 1, m);
     604                result = bdsvd.bidiagonalsvddecomposition(ref w, e, minmn, isupper, false, ref u, 0, ref t2, m, ref vt, ncvt);
     605                blas.copyandtranspose(ref t2, 1, minmn, 1, m, ref u, 1, m, 1, minmn);
     606            }
    588607            return result;
    589608        }
    590        
    591         //
    592         // Simple bidiagonal reduction
    593         //
    594         bidiagonal.tobidiagonal(ref a, m, n, ref tauq, ref taup);
    595         bidiagonal.unpackqfrombidiagonal(ref a, m, n, ref tauq, ncu, ref u);
    596         bidiagonal.unpackptfrombidiagonal(ref a, m, n, ref taup, nrvt, ref vt);
    597         bidiagonal.unpackdiagonalsfrombidiagonal(ref a, m, n, ref isupper, ref w, ref e);
    598         if( additionalmemory<2 | uneeded==0 )
    599         {
    600            
    601             //
    602             // We cant use additional memory or there is no need in such operations
    603             //
    604             result = bdsvd.bidiagonalsvddecomposition(ref w, e, minmn, isupper, false, ref u, nru, ref a, 0, ref vt, ncvt);
    605         }
    606         else
    607         {
    608            
    609             //
    610             // We can use additional memory
    611             //
    612             t2 = new double[minmn+1, m+1];
    613             blas.copyandtranspose(ref u, 1, m, 1, minmn, ref t2, 1, minmn, 1, m);
    614             result = bdsvd.bidiagonalsvddecomposition(ref w, e, minmn, isupper, false, ref u, 0, ref t2, m, ref vt, ncvt);
    615             blas.copyandtranspose(ref t2, 1, minmn, 1, m, ref u, 1, m, 1, minmn);
    616         }
    617         return result;
    618609    }
    619610}
Note: See TracChangeset for help on using the changeset viewer.