Free cookie consent management tool by TermsFeed Policy Generator

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

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

Updated ALGLIB to latest version. #751 (Plugin for for data-modeling with ANN (integrated into CEDMA))

File size: 26.2 KB
RevLine 
[2154]1/*************************************************************************
2Copyright (c) 2005-2007, Sergey Bochkanov (ALGLIB project).
3
[2430]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.
[2154]9
[2430]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.
[2154]14
[2430]15A copy of the GNU General Public License is available at
16http://www.fsf.org/licensing/licenses
[2154]17
[2430]18>>> END OF LICENSE >>>
[2154]19*************************************************************************/
20
21using System;
22
[2430]23namespace alglib
[2154]24{
[2430]25    public class svd
26    {
27        /*************************************************************************
28        Singular value decomposition of a rectangular matrix.
[2154]29
[2430]30        The algorithm calculates the singular value decomposition of a matrix of
31        size MxN: A = U * S * V^T
[2154]32
[2430]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).
[2154]37
[2430]38        Take into account that the subroutine does not return matrix V but V^T.
[2154]39
[2430]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.
[2154]58
[2430]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].
[2154]75
[2430]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;
[2154]106
[2430]107            a = (double[,])a.Clone();
[2154]108
[2430]109            result = true;
110            if( m==0 | n==0 )
[2154]111            {
112                return result;
113            }
[2430]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 )
[2154]126            {
[2430]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            //
[2563]156            if( (double)(m)>(double)(1.6*n) )
[2430]157            {
158                if( uneeded==0 )
[2154]159                {
160                   
161                    //
[2430]162                    // No left singular vectors to be computed
[2154]163                    //
[2430]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;
[2154]177                }
178                else
179                {
180                   
181                    //
[2430]182                    // Left singular vectors (may be full matrix U) to be computed
[2154]183                    //
[2430]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++)
[2154]187                    {
[2430]188                        for(j=0; j<=i-1; j++)
189                        {
190                            a[i,j] = 0;
191                        }
[2154]192                    }
[2430]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;
[2154]219                }
220            }
[2430]221           
222            //
223            // N much larger than M
224            // Use bidiagonal reduction with LQ-decomposition
225            //
[2563]226            if( (double)(n)>(double)(1.6*m) )
[2154]227            {
[2430]228                if( vtneeded==0 )
[2154]229                {
230                   
231                    //
[2430]232                    // No right singular vectors to be computed
[2154]233                    //
[2430]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;
[2154]250                }
251                else
252                {
253                   
254                    //
[2430]255                    // Right singular vectors (may be full matrix VT) to be computed
[2154]256                    //
[2430]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;
[2154]293                }
[2430]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];
[2154]307                blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
[2430]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);
[2154]310                return result;
311            }
[2430]312           
313            //
314            // Simple bidiagonal reduction
315            //
[2154]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);
[2430]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            }
[2154]339            return result;
340        }
[2430]341
342
343        public static bool svddecomposition(double[,] a,
344            int m,
345            int n,
346            int uneeded,
347            int vtneeded,
348            int additionalmemory,
349            ref double[] w,
350            ref double[,] u,
351            ref double[,] vt)
[2154]352        {
[2430]353            bool result = new bool();
354            double[] tauq = new double[0];
355            double[] taup = new double[0];
356            double[] tau = new double[0];
357            double[] e = new double[0];
358            double[] work = new double[0];
359            double[,] t2 = new double[0,0];
360            bool isupper = new bool();
361            int minmn = 0;
362            int ncu = 0;
363            int nrvt = 0;
364            int nru = 0;
365            int ncvt = 0;
366            int i = 0;
367            int j = 0;
368            int im1 = 0;
369            double sm = 0;
370
371            a = (double[,])a.Clone();
372
373            result = true;
374            if( m==0 | n==0 )
375            {
376                return result;
377            }
378            System.Diagnostics.Debug.Assert(uneeded>=0 & uneeded<=2, "SVDDecomposition: wrong parameters!");
379            System.Diagnostics.Debug.Assert(vtneeded>=0 & vtneeded<=2, "SVDDecomposition: wrong parameters!");
380            System.Diagnostics.Debug.Assert(additionalmemory>=0 & additionalmemory<=2, "SVDDecomposition: wrong parameters!");
[2154]381           
382            //
[2430]383            // initialize
[2154]384            //
[2430]385            minmn = Math.Min(m, n);
386            w = new double[minmn+1];
387            ncu = 0;
388            nru = 0;
389            if( uneeded==1 )
390            {
391                nru = m;
392                ncu = minmn;
393                u = new double[nru+1, ncu+1];
394            }
395            if( uneeded==2 )
396            {
397                nru = m;
398                ncu = m;
399                u = new double[nru+1, ncu+1];
400            }
401            nrvt = 0;
402            ncvt = 0;
403            if( vtneeded==1 )
404            {
405                nrvt = minmn;
406                ncvt = n;
407                vt = new double[nrvt+1, ncvt+1];
408            }
409            if( vtneeded==2 )
410            {
411                nrvt = n;
412                ncvt = n;
413                vt = new double[nrvt+1, ncvt+1];
414            }
[2154]415           
416            //
[2430]417            // M much larger than N
418            // Use bidiagonal reduction with QR-decomposition
[2154]419            //
[2563]420            if( (double)(m)>(double)(1.6*n) )
[2154]421            {
[2430]422                if( uneeded==0 )
[2154]423                {
[2430]424                   
425                    //
426                    // No left singular vectors to be computed
427                    //
428                    qr.qrdecomposition(ref a, m, n, ref tau);
429                    for(i=2; i<=n; i++)
[2154]430                    {
[2430]431                        for(j=1; j<=i-1; j++)
432                        {
433                            a[i,j] = 0;
434                        }
[2154]435                    }
[2430]436                    bidiagonal.tobidiagonal(ref a, n, n, ref tauq, ref taup);
437                    bidiagonal.unpackptfrombidiagonal(ref a, n, n, ref taup, nrvt, ref vt);
438                    bidiagonal.unpackdiagonalsfrombidiagonal(ref a, n, n, ref isupper, ref w, ref e);
439                    result = bdsvd.bidiagonalsvddecomposition(ref w, e, n, isupper, false, ref u, 0, ref a, 0, ref vt, ncvt);
440                    return result;
[2154]441                }
[2430]442                else
[2154]443                {
[2430]444                   
445                    //
446                    // Left singular vectors (may be full matrix U) to be computed
447                    //
448                    qr.qrdecomposition(ref a, m, n, ref tau);
449                    qr.unpackqfromqr(ref a, m, n, ref tau, ncu, ref u);
450                    for(i=2; i<=n; i++)
[2154]451                    {
[2430]452                        for(j=1; j<=i-1; j++)
453                        {
454                            a[i,j] = 0;
455                        }
[2154]456                    }
[2430]457                    bidiagonal.tobidiagonal(ref a, n, n, ref tauq, ref taup);
458                    bidiagonal.unpackptfrombidiagonal(ref a, n, n, ref taup, nrvt, ref vt);
459                    bidiagonal.unpackdiagonalsfrombidiagonal(ref a, n, n, ref isupper, ref w, ref e);
460                    if( additionalmemory<1 )
461                    {
462                       
463                        //
464                        // No additional memory can be used
465                        //
466                        bidiagonal.multiplybyqfrombidiagonal(ref a, n, n, ref tauq, ref u, m, n, true, false);
467                        result = bdsvd.bidiagonalsvddecomposition(ref w, e, n, isupper, false, ref u, m, ref a, 0, ref vt, ncvt);
468                    }
469                    else
470                    {
471                       
472                        //
473                        // Large U. Transforming intermediate matrix T2
474                        //
475                        work = new double[Math.Max(m, n)+1];
476                        bidiagonal.unpackqfrombidiagonal(ref a, n, n, ref tauq, n, ref t2);
477                        blas.copymatrix(ref u, 1, m, 1, n, ref a, 1, m, 1, n);
478                        blas.inplacetranspose(ref t2, 1, n, 1, n, ref work);
479                        result = bdsvd.bidiagonalsvddecomposition(ref w, e, n, isupper, false, ref u, 0, ref t2, n, ref vt, ncvt);
480                        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);
481                    }
482                    return result;
[2154]483                }
[2430]484            }
485           
486            //
487            // N much larger than M
488            // Use bidiagonal reduction with LQ-decomposition
489            //
[2563]490            if( (double)(n)>(double)(1.6*m) )
[2430]491            {
492                if( vtneeded==0 )
[2154]493                {
494                   
495                    //
[2430]496                    // No right singular vectors to be computed
[2154]497                    //
[2430]498                    lq.lqdecomposition(ref a, m, n, ref tau);
499                    for(i=1; i<=m-1; i++)
500                    {
501                        for(j=i+1; j<=m; j++)
502                        {
503                            a[i,j] = 0;
504                        }
505                    }
506                    bidiagonal.tobidiagonal(ref a, m, m, ref tauq, ref taup);
507                    bidiagonal.unpackqfrombidiagonal(ref a, m, m, ref tauq, ncu, ref u);
508                    bidiagonal.unpackdiagonalsfrombidiagonal(ref a, m, m, ref isupper, ref w, ref e);
509                    work = new double[m+1];
510                    blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
511                    result = bdsvd.bidiagonalsvddecomposition(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref vt, 0);
512                    blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
513                    return result;
[2154]514                }
515                else
516                {
517                   
518                    //
[2430]519                    // Right singular vectors (may be full matrix VT) to be computed
[2154]520                    //
[2430]521                    lq.lqdecomposition(ref a, m, n, ref tau);
522                    lq.unpackqfromlq(ref a, m, n, ref tau, nrvt, ref vt);
523                    for(i=1; i<=m-1; i++)
524                    {
525                        for(j=i+1; j<=m; j++)
526                        {
527                            a[i,j] = 0;
528                        }
529                    }
530                    bidiagonal.tobidiagonal(ref a, m, m, ref tauq, ref taup);
531                    bidiagonal.unpackqfrombidiagonal(ref a, m, m, ref tauq, ncu, ref u);
532                    bidiagonal.unpackdiagonalsfrombidiagonal(ref a, m, m, ref isupper, ref w, ref e);
[2154]533                    work = new double[Math.Max(m, n)+1];
[2430]534                    blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
535                    if( additionalmemory<1 )
536                    {
537                       
538                        //
539                        // No additional memory available
540                        //
541                        bidiagonal.multiplybypfrombidiagonal(ref a, m, m, ref taup, ref vt, m, n, false, true);
542                        result = bdsvd.bidiagonalsvddecomposition(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref vt, n);
543                    }
544                    else
545                    {
546                       
547                        //
548                        // Large VT. Transforming intermediate matrix T2
549                        //
550                        bidiagonal.unpackptfrombidiagonal(ref a, m, m, ref taup, m, ref t2);
551                        result = bdsvd.bidiagonalsvddecomposition(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref t2, m);
552                        blas.copymatrix(ref vt, 1, m, 1, n, ref a, 1, m, 1, n);
553                        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);
554                    }
555                    blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
556                    return result;
[2154]557                }
558            }
[2430]559           
560            //
561            // M<=N
562            // We can use inplace transposition of U to get rid of columnwise operations
563            //
564            if( m<=n )
[2154]565            {
[2430]566                bidiagonal.tobidiagonal(ref a, m, n, ref tauq, ref taup);
567                bidiagonal.unpackqfrombidiagonal(ref a, m, n, ref tauq, ncu, ref u);
568                bidiagonal.unpackptfrombidiagonal(ref a, m, n, ref taup, nrvt, ref vt);
569                bidiagonal.unpackdiagonalsfrombidiagonal(ref a, m, n, ref isupper, ref w, ref e);
[2154]570                work = new double[m+1];
571                blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
[2430]572                result = bdsvd.bidiagonalsvddecomposition(ref w, e, minmn, isupper, false, ref a, 0, ref u, nru, ref vt, ncvt);
[2154]573                blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
574                return result;
575            }
[2430]576           
577            //
578            // Simple bidiagonal reduction
579            //
580            bidiagonal.tobidiagonal(ref a, m, n, ref tauq, ref taup);
581            bidiagonal.unpackqfrombidiagonal(ref a, m, n, ref tauq, ncu, ref u);
582            bidiagonal.unpackptfrombidiagonal(ref a, m, n, ref taup, nrvt, ref vt);
583            bidiagonal.unpackdiagonalsfrombidiagonal(ref a, m, n, ref isupper, ref w, ref e);
584            if( additionalmemory<2 | uneeded==0 )
585            {
586               
587                //
588                // We cant use additional memory or there is no need in such operations
589                //
590                result = bdsvd.bidiagonalsvddecomposition(ref w, e, minmn, isupper, false, ref u, nru, ref a, 0, ref vt, ncvt);
591            }
[2154]592            else
593            {
594               
595                //
[2430]596                // We can use additional memory
[2154]597                //
[2430]598                t2 = new double[minmn+1, m+1];
599                blas.copyandtranspose(ref u, 1, m, 1, minmn, ref t2, 1, minmn, 1, m);
600                result = bdsvd.bidiagonalsvddecomposition(ref w, e, minmn, isupper, false, ref u, 0, ref t2, m, ref vt, ncvt);
601                blas.copyandtranspose(ref t2, 1, minmn, 1, m, ref u, 1, m, 1, minmn);
[2154]602            }
603            return result;
604        }
605    }
606}
Note: See TracBrowser for help on using the repository browser.