Free cookie consent management tool by TermsFeed Policy Generator

source: branches/Persistence Test/ALGLIB/svd.cs @ 2527

Last change on this file since 2527 was 2430, checked in by gkronber, 15 years ago

Moved ALGLIB code into a separate plugin. #783
File size: 26.4 KB
Line 
1/*************************************************************************
2Copyright (c) 2005-2007, Sergey Bochkanov (ALGLIB project).
3
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 >>>
19*************************************************************************/
20
21using System;
22
23namespace alglib
24{
25    public class svd
26    {
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)
88        {
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            {
112                return result;
113            }
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;
177                }
178                else
179                {
180                   
181                    //
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);
269                    work = new double[Math.Max(m, n)+1];
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);
306                work = new double[m+1];
307                blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
308                result = bdsvd.rmatrixbdsvd(ref w, e, minmn, isupper, false, ref a, 0, ref u, nru, ref vt, ncvt);
309                blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
310                return result;
311            }
312           
313            //
314            // Simple bidiagonal reduction
315            //
316            bidiagonal.rmatrixbd(ref a, m, n, ref tauq, ref taup);
317            bidiagonal.rmatrixbdunpackq(ref a, m, n, ref tauq, ncu, ref u);
318            bidiagonal.rmatrixbdunpackpt(ref a, m, n, ref taup, nrvt, ref vt);
319            bidiagonal.rmatrixbdunpackdiagonals(ref a, m, n, ref isupper, ref w, ref e);
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            }
339            return result;
340        }
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)
356        {
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            {
380                return result;
381            }
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;
445                }
446                else
447                {
448                   
449                    //
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);
537                    work = new double[Math.Max(m, n)+1];
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);
574                work = new double[m+1];
575                blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
576                result = bdsvd.bidiagonalsvddecomposition(ref w, e, minmn, isupper, false, ref a, 0, ref u, nru, ref vt, ncvt);
577                blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
578                return result;
579            }
580           
581            //
582            // Simple bidiagonal reduction
583            //
584            bidiagonal.tobidiagonal(ref a, m, n, ref tauq, ref taup);
585            bidiagonal.unpackqfrombidiagonal(ref a, m, n, ref tauq, ncu, ref u);
586            bidiagonal.unpackptfrombidiagonal(ref a, m, n, ref taup, nrvt, ref vt);
587            bidiagonal.unpackdiagonalsfrombidiagonal(ref a, m, n, ref isupper, ref w, ref e);
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            }
607            return result;
608        }
609    }
610}
Note: See TracBrowser for help on using the repository browser.