Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/ALGLIB/bdsvd.cs @ 2636

Last change on this file since 2636 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: 50.6 KB
Line 
1/*************************************************************************
2Copyright (c) 1992-2007 The University of Tennessee.  All rights reserved.
3
4Contributors:
5    * Sergey Bochkanov (ALGLIB project). Translation from FORTRAN to
6      pseudocode.
7
8See subroutines comments for additional copyrights.
9
10>>> SOURCE LICENSE >>>
11This program is free software; you can redistribute it and/or modify
12it under the terms of the GNU General Public License as published by
13the Free Software Foundation (www.fsf.org); either version 2 of the
14License, or (at your option) any later version.
15
16This program is distributed in the hope that it will be useful,
17but WITHOUT ANY WARRANTY; without even the implied warranty of
18MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
19GNU General Public License for more details.
20
21A copy of the GNU General Public License is available at
22http://www.fsf.org/licensing/licenses
23
24>>> END OF LICENSE >>>
25*************************************************************************/
26
27using System;
28
29namespace alglib
30{
31    public class bdsvd
32    {
33        /*************************************************************************
34        Singular value decomposition of a bidiagonal matrix (extended algorithm)
35
36        The algorithm performs the singular value decomposition  of  a  bidiagonal
37        matrix B (upper or lower) representing it as B = Q*S*P^T, where Q and  P -
38        orthogonal matrices, S - diagonal matrix with non-negative elements on the
39        main diagonal, in descending order.
40
41        The  algorithm  finds  singular  values.  In  addition,  the algorithm can
42        calculate  matrices  Q  and P (more precisely, not the matrices, but their
43        product  with  given  matrices U and VT - U*Q and (P^T)*VT)).  Of  course,
44        matrices U and VT can be of any type, including identity. Furthermore, the
45        algorithm can calculate Q'*C (this product is calculated more  effectively
46        than U*Q,  because  this calculation operates with rows instead  of matrix
47        columns).
48
49        The feature of the algorithm is its ability to find  all  singular  values
50        including those which are arbitrarily close to 0  with  relative  accuracy
51        close to  machine precision. If the parameter IsFractionalAccuracyRequired
52        is set to True, all singular values will have high relative accuracy close
53        to machine precision. If the parameter is set to False, only  the  biggest
54        singular value will have relative accuracy  close  to  machine  precision.
55        The absolute error of other singular values is equal to the absolute error
56        of the biggest singular value.
57
58        Input parameters:
59            D       -   main diagonal of matrix B.
60                        Array whose index ranges within [0..N-1].
61            E       -   superdiagonal (or subdiagonal) of matrix B.
62                        Array whose index ranges within [0..N-2].
63            N       -   size of matrix B.
64            IsUpper -   True, if the matrix is upper bidiagonal.
65            IsFractionalAccuracyRequired -
66                        accuracy to search singular values with.
67            U       -   matrix to be multiplied by Q.
68                        Array whose indexes range within [0..NRU-1, 0..N-1].
69                        The matrix can be bigger, in that case only the  submatrix
70                        [0..NRU-1, 0..N-1] will be multiplied by Q.
71            NRU     -   number of rows in matrix U.
72            C       -   matrix to be multiplied by Q'.
73                        Array whose indexes range within [0..N-1, 0..NCC-1].
74                        The matrix can be bigger, in that case only the  submatrix
75                        [0..N-1, 0..NCC-1] will be multiplied by Q'.
76            NCC     -   number of columns in matrix C.
77            VT      -   matrix to be multiplied by P^T.
78                        Array whose indexes range within [0..N-1, 0..NCVT-1].
79                        The matrix can be bigger, in that case only the  submatrix
80                        [0..N-1, 0..NCVT-1] will be multiplied by P^T.
81            NCVT    -   number of columns in matrix VT.
82
83        Output parameters:
84            D       -   singular values of matrix B in descending order.
85            U       -   if NRU>0, contains matrix U*Q.
86            VT      -   if NCVT>0, contains matrix (P^T)*VT.
87            C       -   if NCC>0, contains matrix Q'*C.
88
89        Result:
90            True, if the algorithm has converged.
91            False, if the algorithm hasn't converged (rare case).
92
93        Additional information:
94            The type of convergence is controlled by the internal  parameter  TOL.
95            If the parameter is greater than 0, the singular values will have
96            relative accuracy TOL. If TOL<0, the singular values will have
97            absolute accuracy ABS(TOL)*norm(B).
98            By default, |TOL| falls within the range of 10*Epsilon and 100*Epsilon,
99            where Epsilon is the machine precision. It is not  recommended  to  use
100            TOL less than 10*Epsilon since this will  considerably  slow  down  the
101            algorithm and may not lead to error decreasing.
102        History:
103            * 31 March, 2007.
104                changed MAXITR from 6 to 12.
105
106          -- LAPACK routine (version 3.0) --
107             Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
108             Courant Institute, Argonne National Lab, and Rice University
109             October 31, 1999.
110        *************************************************************************/
111        public static bool rmatrixbdsvd(ref double[] d,
112            double[] e,
113            int n,
114            bool isupper,
115            bool isfractionalaccuracyrequired,
116            ref double[,] u,
117            int nru,
118            ref double[,] c,
119            int ncc,
120            ref double[,] vt,
121            int ncvt)
122        {
123            bool result = new bool();
124            double[] d1 = new double[0];
125            double[] e1 = new double[0];
126            int i_ = 0;
127            int i1_ = 0;
128
129            e = (double[])e.Clone();
130
131            d1 = new double[n+1];
132            i1_ = (0) - (1);
133            for(i_=1; i_<=n;i_++)
134            {
135                d1[i_] = d[i_+i1_];
136            }
137            if( n>1 )
138            {
139                e1 = new double[n-1+1];
140                i1_ = (0) - (1);
141                for(i_=1; i_<=n-1;i_++)
142                {
143                    e1[i_] = e[i_+i1_];
144                }
145            }
146            result = bidiagonalsvddecompositioninternal(ref d1, e1, n, isupper, isfractionalaccuracyrequired, ref u, 0, nru, ref c, 0, ncc, ref vt, 0, ncvt);
147            i1_ = (1) - (0);
148            for(i_=0; i_<=n-1;i_++)
149            {
150                d[i_] = d1[i_+i1_];
151            }
152            return result;
153        }
154
155
156        public static bool bidiagonalsvddecomposition(ref double[] d,
157            double[] e,
158            int n,
159            bool isupper,
160            bool isfractionalaccuracyrequired,
161            ref double[,] u,
162            int nru,
163            ref double[,] c,
164            int ncc,
165            ref double[,] vt,
166            int ncvt)
167        {
168            bool result = new bool();
169
170            e = (double[])e.Clone();
171
172            result = bidiagonalsvddecompositioninternal(ref d, e, n, isupper, isfractionalaccuracyrequired, ref u, 1, nru, ref c, 1, ncc, ref vt, 1, ncvt);
173            return result;
174        }
175
176
177        /*************************************************************************
178        Internal working subroutine for bidiagonal decomposition
179        *************************************************************************/
180        private static bool bidiagonalsvddecompositioninternal(ref double[] d,
181            double[] e,
182            int n,
183            bool isupper,
184            bool isfractionalaccuracyrequired,
185            ref double[,] u,
186            int ustart,
187            int nru,
188            ref double[,] c,
189            int cstart,
190            int ncc,
191            ref double[,] vt,
192            int vstart,
193            int ncvt)
194        {
195            bool result = new bool();
196            int i = 0;
197            int idir = 0;
198            int isub = 0;
199            int iter = 0;
200            int j = 0;
201            int ll = 0;
202            int lll = 0;
203            int m = 0;
204            int maxit = 0;
205            int oldll = 0;
206            int oldm = 0;
207            double abse = 0;
208            double abss = 0;
209            double cosl = 0;
210            double cosr = 0;
211            double cs = 0;
212            double eps = 0;
213            double f = 0;
214            double g = 0;
215            double h = 0;
216            double mu = 0;
217            double oldcs = 0;
218            double oldsn = 0;
219            double r = 0;
220            double shift = 0;
221            double sigmn = 0;
222            double sigmx = 0;
223            double sinl = 0;
224            double sinr = 0;
225            double sll = 0;
226            double smax = 0;
227            double smin = 0;
228            double sminl = 0;
229            double sminlo = 0;
230            double sminoa = 0;
231            double sn = 0;
232            double thresh = 0;
233            double tol = 0;
234            double tolmul = 0;
235            double unfl = 0;
236            double[] work0 = new double[0];
237            double[] work1 = new double[0];
238            double[] work2 = new double[0];
239            double[] work3 = new double[0];
240            int maxitr = 0;
241            bool matrixsplitflag = new bool();
242            bool iterflag = new bool();
243            double[] utemp = new double[0];
244            double[] vttemp = new double[0];
245            double[] ctemp = new double[0];
246            double[] etemp = new double[0];
247            bool rightside = new bool();
248            bool fwddir = new bool();
249            double tmp = 0;
250            int mm1 = 0;
251            int mm0 = 0;
252            bool bchangedir = new bool();
253            int uend = 0;
254            int cend = 0;
255            int vend = 0;
256            int i_ = 0;
257
258            e = (double[])e.Clone();
259
260            result = true;
261            if( n==0 )
262            {
263                return result;
264            }
265            if( n==1 )
266            {
267                if( (double)(d[1])<(double)(0) )
268                {
269                    d[1] = -d[1];
270                    if( ncvt>0 )
271                    {
272                        for(i_=vstart; i_<=vstart+ncvt-1;i_++)
273                        {
274                            vt[vstart,i_] = -1*vt[vstart,i_];
275                        }
276                    }
277                }
278                return result;
279            }
280           
281            //
282            // init
283            //
284            work0 = new double[n-1+1];
285            work1 = new double[n-1+1];
286            work2 = new double[n-1+1];
287            work3 = new double[n-1+1];
288            uend = ustart+Math.Max(nru-1, 0);
289            vend = vstart+Math.Max(ncvt-1, 0);
290            cend = cstart+Math.Max(ncc-1, 0);
291            utemp = new double[uend+1];
292            vttemp = new double[vend+1];
293            ctemp = new double[cend+1];
294            maxitr = 12;
295            rightside = true;
296            fwddir = true;
297           
298            //
299            // resize E from N-1 to N
300            //
301            etemp = new double[n+1];
302            for(i=1; i<=n-1; i++)
303            {
304                etemp[i] = e[i];
305            }
306            e = new double[n+1];
307            for(i=1; i<=n-1; i++)
308            {
309                e[i] = etemp[i];
310            }
311            e[n] = 0;
312            idir = 0;
313           
314            //
315            // Get machine constants
316            //
317            eps = AP.Math.MachineEpsilon;
318            unfl = AP.Math.MinRealNumber;
319           
320            //
321            // If matrix lower bidiagonal, rotate to be upper bidiagonal
322            // by applying Givens rotations on the left
323            //
324            if( !isupper )
325            {
326                for(i=1; i<=n-1; i++)
327                {
328                    rotations.generaterotation(d[i], e[i], ref cs, ref sn, ref r);
329                    d[i] = r;
330                    e[i] = sn*d[i+1];
331                    d[i+1] = cs*d[i+1];
332                    work0[i] = cs;
333                    work1[i] = sn;
334                }
335               
336                //
337                // Update singular vectors if desired
338                //
339                if( nru>0 )
340                {
341                    rotations.applyrotationsfromtheright(fwddir, ustart, uend, 1+ustart-1, n+ustart-1, ref work0, ref work1, ref u, ref utemp);
342                }
343                if( ncc>0 )
344                {
345                    rotations.applyrotationsfromtheleft(fwddir, 1+cstart-1, n+cstart-1, cstart, cend, ref work0, ref work1, ref c, ref ctemp);
346                }
347            }
348           
349            //
350            // Compute singular values to relative accuracy TOL
351            // (By setting TOL to be negative, algorithm will compute
352            // singular values to absolute accuracy ABS(TOL)*norm(input matrix))
353            //
354            tolmul = Math.Max(10, Math.Min(100, Math.Pow(eps, -0.125)));
355            tol = tolmul*eps;
356            if( !isfractionalaccuracyrequired )
357            {
358                tol = -tol;
359            }
360           
361            //
362            // Compute approximate maximum, minimum singular values
363            //
364            smax = 0;
365            for(i=1; i<=n; i++)
366            {
367                smax = Math.Max(smax, Math.Abs(d[i]));
368            }
369            for(i=1; i<=n-1; i++)
370            {
371                smax = Math.Max(smax, Math.Abs(e[i]));
372            }
373            sminl = 0;
374            if( (double)(tol)>=(double)(0) )
375            {
376               
377                //
378                // Relative accuracy desired
379                //
380                sminoa = Math.Abs(d[1]);
381                if( (double)(sminoa)!=(double)(0) )
382                {
383                    mu = sminoa;
384                    for(i=2; i<=n; i++)
385                    {
386                        mu = Math.Abs(d[i])*(mu/(mu+Math.Abs(e[i-1])));
387                        sminoa = Math.Min(sminoa, mu);
388                        if( (double)(sminoa)==(double)(0) )
389                        {
390                            break;
391                        }
392                    }
393                }
394                sminoa = sminoa/Math.Sqrt(n);
395                thresh = Math.Max(tol*sminoa, maxitr*n*n*unfl);
396            }
397            else
398            {
399               
400                //
401                // Absolute accuracy desired
402                //
403                thresh = Math.Max(Math.Abs(tol)*smax, maxitr*n*n*unfl);
404            }
405           
406            //
407            // Prepare for main iteration loop for the singular values
408            // (MAXIT is the maximum number of passes through the inner
409            // loop permitted before nonconvergence signalled.)
410            //
411            maxit = maxitr*n*n;
412            iter = 0;
413            oldll = -1;
414            oldm = -1;
415           
416            //
417            // M points to last element of unconverged part of matrix
418            //
419            m = n;
420           
421            //
422            // Begin main iteration loop
423            //
424            while( true )
425            {
426               
427                //
428                // Check for convergence or exceeding iteration count
429                //
430                if( m<=1 )
431                {
432                    break;
433                }
434                if( iter>maxit )
435                {
436                    result = false;
437                    return result;
438                }
439               
440                //
441                // Find diagonal block of matrix to work on
442                //
443                if( (double)(tol)<(double)(0) & (double)(Math.Abs(d[m]))<=(double)(thresh) )
444                {
445                    d[m] = 0;
446                }
447                smax = Math.Abs(d[m]);
448                smin = smax;
449                matrixsplitflag = false;
450                for(lll=1; lll<=m-1; lll++)
451                {
452                    ll = m-lll;
453                    abss = Math.Abs(d[ll]);
454                    abse = Math.Abs(e[ll]);
455                    if( (double)(tol)<(double)(0) & (double)(abss)<=(double)(thresh) )
456                    {
457                        d[ll] = 0;
458                    }
459                    if( (double)(abse)<=(double)(thresh) )
460                    {
461                        matrixsplitflag = true;
462                        break;
463                    }
464                    smin = Math.Min(smin, abss);
465                    smax = Math.Max(smax, Math.Max(abss, abse));
466                }
467                if( !matrixsplitflag )
468                {
469                    ll = 0;
470                }
471                else
472                {
473                   
474                    //
475                    // Matrix splits since E(LL) = 0
476                    //
477                    e[ll] = 0;
478                    if( ll==m-1 )
479                    {
480                       
481                        //
482                        // Convergence of bottom singular value, return to top of loop
483                        //
484                        m = m-1;
485                        continue;
486                    }
487                }
488                ll = ll+1;
489               
490                //
491                // E(LL) through E(M-1) are nonzero, E(LL-1) is zero
492                //
493                if( ll==m-1 )
494                {
495                   
496                    //
497                    // 2 by 2 block, handle separately
498                    //
499                    svdv2x2(d[m-1], e[m-1], d[m], ref sigmn, ref sigmx, ref sinr, ref cosr, ref sinl, ref cosl);
500                    d[m-1] = sigmx;
501                    e[m-1] = 0;
502                    d[m] = sigmn;
503                   
504                    //
505                    // Compute singular vectors, if desired
506                    //
507                    if( ncvt>0 )
508                    {
509                        mm0 = m+(vstart-1);
510                        mm1 = m-1+(vstart-1);
511                        for(i_=vstart; i_<=vend;i_++)
512                        {
513                            vttemp[i_] = cosr*vt[mm1,i_];
514                        }
515                        for(i_=vstart; i_<=vend;i_++)
516                        {
517                            vttemp[i_] = vttemp[i_] + sinr*vt[mm0,i_];
518                        }
519                        for(i_=vstart; i_<=vend;i_++)
520                        {
521                            vt[mm0,i_] = cosr*vt[mm0,i_];
522                        }
523                        for(i_=vstart; i_<=vend;i_++)
524                        {
525                            vt[mm0,i_] = vt[mm0,i_] - sinr*vt[mm1,i_];
526                        }
527                        for(i_=vstart; i_<=vend;i_++)
528                        {
529                            vt[mm1,i_] = vttemp[i_];
530                        }
531                    }
532                    if( nru>0 )
533                    {
534                        mm0 = m+ustart-1;
535                        mm1 = m-1+ustart-1;
536                        for(i_=ustart; i_<=uend;i_++)
537                        {
538                            utemp[i_] = cosl*u[i_,mm1];
539                        }
540                        for(i_=ustart; i_<=uend;i_++)
541                        {
542                            utemp[i_] = utemp[i_] + sinl*u[i_,mm0];
543                        }
544                        for(i_=ustart; i_<=uend;i_++)
545                        {
546                            u[i_,mm0] = cosl*u[i_,mm0];
547                        }
548                        for(i_=ustart; i_<=uend;i_++)
549                        {
550                            u[i_,mm0] = u[i_,mm0] - sinl*u[i_,mm1];
551                        }
552                        for(i_=ustart; i_<=uend;i_++)
553                        {
554                            u[i_,mm1] = utemp[i_];
555                        }
556                    }
557                    if( ncc>0 )
558                    {
559                        mm0 = m+cstart-1;
560                        mm1 = m-1+cstart-1;
561                        for(i_=cstart; i_<=cend;i_++)
562                        {
563                            ctemp[i_] = cosl*c[mm1,i_];
564                        }
565                        for(i_=cstart; i_<=cend;i_++)
566                        {
567                            ctemp[i_] = ctemp[i_] + sinl*c[mm0,i_];
568                        }
569                        for(i_=cstart; i_<=cend;i_++)
570                        {
571                            c[mm0,i_] = cosl*c[mm0,i_];
572                        }
573                        for(i_=cstart; i_<=cend;i_++)
574                        {
575                            c[mm0,i_] = c[mm0,i_] - sinl*c[mm1,i_];
576                        }
577                        for(i_=cstart; i_<=cend;i_++)
578                        {
579                            c[mm1,i_] = ctemp[i_];
580                        }
581                    }
582                    m = m-2;
583                    continue;
584                }
585               
586                //
587                // If working on new submatrix, choose shift direction
588                // (from larger end diagonal element towards smaller)
589                //
590                // Previously was
591                //     "if (LL>OLDM) or (M<OLDLL) then"
592                // fixed thanks to Michael Rolle < m@rolle.name >
593                // Very strange that LAPACK still contains it.
594                //
595                bchangedir = false;
596                if( idir==1 & (double)(Math.Abs(d[ll]))<(double)(1.0E-3*Math.Abs(d[m])) )
597                {
598                    bchangedir = true;
599                }
600                if( idir==2 & (double)(Math.Abs(d[m]))<(double)(1.0E-3*Math.Abs(d[ll])) )
601                {
602                    bchangedir = true;
603                }
604                if( ll!=oldll | m!=oldm | bchangedir )
605                {
606                    if( (double)(Math.Abs(d[ll]))>=(double)(Math.Abs(d[m])) )
607                    {
608                       
609                        //
610                        // Chase bulge from top (big end) to bottom (small end)
611                        //
612                        idir = 1;
613                    }
614                    else
615                    {
616                       
617                        //
618                        // Chase bulge from bottom (big end) to top (small end)
619                        //
620                        idir = 2;
621                    }
622                }
623               
624                //
625                // Apply convergence tests
626                //
627                if( idir==1 )
628                {
629                   
630                    //
631                    // Run convergence test in forward direction
632                    // First apply standard test to bottom of matrix
633                    //
634                    if( (double)(Math.Abs(e[m-1]))<=(double)(Math.Abs(tol)*Math.Abs(d[m])) | (double)(tol)<(double)(0) & (double)(Math.Abs(e[m-1]))<=(double)(thresh) )
635                    {
636                        e[m-1] = 0;
637                        continue;
638                    }
639                    if( (double)(tol)>=(double)(0) )
640                    {
641                       
642                        //
643                        // If relative accuracy desired,
644                        // apply convergence criterion forward
645                        //
646                        mu = Math.Abs(d[ll]);
647                        sminl = mu;
648                        iterflag = false;
649                        for(lll=ll; lll<=m-1; lll++)
650                        {
651                            if( (double)(Math.Abs(e[lll]))<=(double)(tol*mu) )
652                            {
653                                e[lll] = 0;
654                                iterflag = true;
655                                break;
656                            }
657                            sminlo = sminl;
658                            mu = Math.Abs(d[lll+1])*(mu/(mu+Math.Abs(e[lll])));
659                            sminl = Math.Min(sminl, mu);
660                        }
661                        if( iterflag )
662                        {
663                            continue;
664                        }
665                    }
666                }
667                else
668                {
669                   
670                    //
671                    // Run convergence test in backward direction
672                    // First apply standard test to top of matrix
673                    //
674                    if( (double)(Math.Abs(e[ll]))<=(double)(Math.Abs(tol)*Math.Abs(d[ll])) | (double)(tol)<(double)(0) & (double)(Math.Abs(e[ll]))<=(double)(thresh) )
675                    {
676                        e[ll] = 0;
677                        continue;
678                    }
679                    if( (double)(tol)>=(double)(0) )
680                    {
681                       
682                        //
683                        // If relative accuracy desired,
684                        // apply convergence criterion backward
685                        //
686                        mu = Math.Abs(d[m]);
687                        sminl = mu;
688                        iterflag = false;
689                        for(lll=m-1; lll>=ll; lll--)
690                        {
691                            if( (double)(Math.Abs(e[lll]))<=(double)(tol*mu) )
692                            {
693                                e[lll] = 0;
694                                iterflag = true;
695                                break;
696                            }
697                            sminlo = sminl;
698                            mu = Math.Abs(d[lll])*(mu/(mu+Math.Abs(e[lll])));
699                            sminl = Math.Min(sminl, mu);
700                        }
701                        if( iterflag )
702                        {
703                            continue;
704                        }
705                    }
706                }
707                oldll = ll;
708                oldm = m;
709               
710                //
711                // Compute shift.  First, test if shifting would ruin relative
712                // accuracy, and if so set the shift to zero.
713                //
714                if( (double)(tol)>=(double)(0) & (double)(n*tol*(sminl/smax))<=(double)(Math.Max(eps, 0.01*tol)) )
715                {
716                   
717                    //
718                    // Use a zero shift to avoid loss of relative accuracy
719                    //
720                    shift = 0;
721                }
722                else
723                {
724                   
725                    //
726                    // Compute the shift from 2-by-2 block at end of matrix
727                    //
728                    if( idir==1 )
729                    {
730                        sll = Math.Abs(d[ll]);
731                        svd2x2(d[m-1], e[m-1], d[m], ref shift, ref r);
732                    }
733                    else
734                    {
735                        sll = Math.Abs(d[m]);
736                        svd2x2(d[ll], e[ll], d[ll+1], ref shift, ref r);
737                    }
738                   
739                    //
740                    // Test if shift negligible, and if so set to zero
741                    //
742                    if( (double)(sll)>(double)(0) )
743                    {
744                        if( (double)(AP.Math.Sqr(shift/sll))<(double)(eps) )
745                        {
746                            shift = 0;
747                        }
748                    }
749                }
750               
751                //
752                // Increment iteration count
753                //
754                iter = iter+m-ll;
755               
756                //
757                // If SHIFT = 0, do simplified QR iteration
758                //
759                if( (double)(shift)==(double)(0) )
760                {
761                    if( idir==1 )
762                    {
763                       
764                        //
765                        // Chase bulge from top to bottom
766                        // Save cosines and sines for later singular vector updates
767                        //
768                        cs = 1;
769                        oldcs = 1;
770                        for(i=ll; i<=m-1; i++)
771                        {
772                            rotations.generaterotation(d[i]*cs, e[i], ref cs, ref sn, ref r);
773                            if( i>ll )
774                            {
775                                e[i-1] = oldsn*r;
776                            }
777                            rotations.generaterotation(oldcs*r, d[i+1]*sn, ref oldcs, ref oldsn, ref tmp);
778                            d[i] = tmp;
779                            work0[i-ll+1] = cs;
780                            work1[i-ll+1] = sn;
781                            work2[i-ll+1] = oldcs;
782                            work3[i-ll+1] = oldsn;
783                        }
784                        h = d[m]*cs;
785                        d[m] = h*oldcs;
786                        e[m-1] = h*oldsn;
787                       
788                        //
789                        // Update singular vectors
790                        //
791                        if( ncvt>0 )
792                        {
793                            rotations.applyrotationsfromtheleft(fwddir, ll+vstart-1, m+vstart-1, vstart, vend, ref work0, ref work1, ref vt, ref vttemp);
794                        }
795                        if( nru>0 )
796                        {
797                            rotations.applyrotationsfromtheright(fwddir, ustart, uend, ll+ustart-1, m+ustart-1, ref work2, ref work3, ref u, ref utemp);
798                        }
799                        if( ncc>0 )
800                        {
801                            rotations.applyrotationsfromtheleft(fwddir, ll+cstart-1, m+cstart-1, cstart, cend, ref work2, ref work3, ref c, ref ctemp);
802                        }
803                       
804                        //
805                        // Test convergence
806                        //
807                        if( (double)(Math.Abs(e[m-1]))<=(double)(thresh) )
808                        {
809                            e[m-1] = 0;
810                        }
811                    }
812                    else
813                    {
814                       
815                        //
816                        // Chase bulge from bottom to top
817                        // Save cosines and sines for later singular vector updates
818                        //
819                        cs = 1;
820                        oldcs = 1;
821                        for(i=m; i>=ll+1; i--)
822                        {
823                            rotations.generaterotation(d[i]*cs, e[i-1], ref cs, ref sn, ref r);
824                            if( i<m )
825                            {
826                                e[i] = oldsn*r;
827                            }
828                            rotations.generaterotation(oldcs*r, d[i-1]*sn, ref oldcs, ref oldsn, ref tmp);
829                            d[i] = tmp;
830                            work0[i-ll] = cs;
831                            work1[i-ll] = -sn;
832                            work2[i-ll] = oldcs;
833                            work3[i-ll] = -oldsn;
834                        }
835                        h = d[ll]*cs;
836                        d[ll] = h*oldcs;
837                        e[ll] = h*oldsn;
838                       
839                        //
840                        // Update singular vectors
841                        //
842                        if( ncvt>0 )
843                        {
844                            rotations.applyrotationsfromtheleft(!fwddir, ll+vstart-1, m+vstart-1, vstart, vend, ref work2, ref work3, ref vt, ref vttemp);
845                        }
846                        if( nru>0 )
847                        {
848                            rotations.applyrotationsfromtheright(!fwddir, ustart, uend, ll+ustart-1, m+ustart-1, ref work0, ref work1, ref u, ref utemp);
849                        }
850                        if( ncc>0 )
851                        {
852                            rotations.applyrotationsfromtheleft(!fwddir, ll+cstart-1, m+cstart-1, cstart, cend, ref work0, ref work1, ref c, ref ctemp);
853                        }
854                       
855                        //
856                        // Test convergence
857                        //
858                        if( (double)(Math.Abs(e[ll]))<=(double)(thresh) )
859                        {
860                            e[ll] = 0;
861                        }
862                    }
863                }
864                else
865                {
866                   
867                    //
868                    // Use nonzero shift
869                    //
870                    if( idir==1 )
871                    {
872                       
873                        //
874                        // Chase bulge from top to bottom
875                        // Save cosines and sines for later singular vector updates
876                        //
877                        f = (Math.Abs(d[ll])-shift)*(extsignbdsqr(1, d[ll])+shift/d[ll]);
878                        g = e[ll];
879                        for(i=ll; i<=m-1; i++)
880                        {
881                            rotations.generaterotation(f, g, ref cosr, ref sinr, ref r);
882                            if( i>ll )
883                            {
884                                e[i-1] = r;
885                            }
886                            f = cosr*d[i]+sinr*e[i];
887                            e[i] = cosr*e[i]-sinr*d[i];
888                            g = sinr*d[i+1];
889                            d[i+1] = cosr*d[i+1];
890                            rotations.generaterotation(f, g, ref cosl, ref sinl, ref r);
891                            d[i] = r;
892                            f = cosl*e[i]+sinl*d[i+1];
893                            d[i+1] = cosl*d[i+1]-sinl*e[i];
894                            if( i<m-1 )
895                            {
896                                g = sinl*e[i+1];
897                                e[i+1] = cosl*e[i+1];
898                            }
899                            work0[i-ll+1] = cosr;
900                            work1[i-ll+1] = sinr;
901                            work2[i-ll+1] = cosl;
902                            work3[i-ll+1] = sinl;
903                        }
904                        e[m-1] = f;
905                       
906                        //
907                        // Update singular vectors
908                        //
909                        if( ncvt>0 )
910                        {
911                            rotations.applyrotationsfromtheleft(fwddir, ll+vstart-1, m+vstart-1, vstart, vend, ref work0, ref work1, ref vt, ref vttemp);
912                        }
913                        if( nru>0 )
914                        {
915                            rotations.applyrotationsfromtheright(fwddir, ustart, uend, ll+ustart-1, m+ustart-1, ref work2, ref work3, ref u, ref utemp);
916                        }
917                        if( ncc>0 )
918                        {
919                            rotations.applyrotationsfromtheleft(fwddir, ll+cstart-1, m+cstart-1, cstart, cend, ref work2, ref work3, ref c, ref ctemp);
920                        }
921                       
922                        //
923                        // Test convergence
924                        //
925                        if( (double)(Math.Abs(e[m-1]))<=(double)(thresh) )
926                        {
927                            e[m-1] = 0;
928                        }
929                    }
930                    else
931                    {
932                       
933                        //
934                        // Chase bulge from bottom to top
935                        // Save cosines and sines for later singular vector updates
936                        //
937                        f = (Math.Abs(d[m])-shift)*(extsignbdsqr(1, d[m])+shift/d[m]);
938                        g = e[m-1];
939                        for(i=m; i>=ll+1; i--)
940                        {
941                            rotations.generaterotation(f, g, ref cosr, ref sinr, ref r);
942                            if( i<m )
943                            {
944                                e[i] = r;
945                            }
946                            f = cosr*d[i]+sinr*e[i-1];
947                            e[i-1] = cosr*e[i-1]-sinr*d[i];
948                            g = sinr*d[i-1];
949                            d[i-1] = cosr*d[i-1];
950                            rotations.generaterotation(f, g, ref cosl, ref sinl, ref r);
951                            d[i] = r;
952                            f = cosl*e[i-1]+sinl*d[i-1];
953                            d[i-1] = cosl*d[i-1]-sinl*e[i-1];
954                            if( i>ll+1 )
955                            {
956                                g = sinl*e[i-2];
957                                e[i-2] = cosl*e[i-2];
958                            }
959                            work0[i-ll] = cosr;
960                            work1[i-ll] = -sinr;
961                            work2[i-ll] = cosl;
962                            work3[i-ll] = -sinl;
963                        }
964                        e[ll] = f;
965                       
966                        //
967                        // Test convergence
968                        //
969                        if( (double)(Math.Abs(e[ll]))<=(double)(thresh) )
970                        {
971                            e[ll] = 0;
972                        }
973                       
974                        //
975                        // Update singular vectors if desired
976                        //
977                        if( ncvt>0 )
978                        {
979                            rotations.applyrotationsfromtheleft(!fwddir, ll+vstart-1, m+vstart-1, vstart, vend, ref work2, ref work3, ref vt, ref vttemp);
980                        }
981                        if( nru>0 )
982                        {
983                            rotations.applyrotationsfromtheright(!fwddir, ustart, uend, ll+ustart-1, m+ustart-1, ref work0, ref work1, ref u, ref utemp);
984                        }
985                        if( ncc>0 )
986                        {
987                            rotations.applyrotationsfromtheleft(!fwddir, ll+cstart-1, m+cstart-1, cstart, cend, ref work0, ref work1, ref c, ref ctemp);
988                        }
989                    }
990                }
991               
992                //
993                // QR iteration finished, go back and check convergence
994                //
995                continue;
996            }
997           
998            //
999            // All singular values converged, so make them positive
1000            //
1001            for(i=1; i<=n; i++)
1002            {
1003                if( (double)(d[i])<(double)(0) )
1004                {
1005                    d[i] = -d[i];
1006                   
1007                    //
1008                    // Change sign of singular vectors, if desired
1009                    //
1010                    if( ncvt>0 )
1011                    {
1012                        for(i_=vstart; i_<=vend;i_++)
1013                        {
1014                            vt[i+vstart-1,i_] = -1*vt[i+vstart-1,i_];
1015                        }
1016                    }
1017                }
1018            }
1019           
1020            //
1021            // Sort the singular values into decreasing order (insertion sort on
1022            // singular values, but only one transposition per singular vector)
1023            //
1024            for(i=1; i<=n-1; i++)
1025            {
1026               
1027                //
1028                // Scan for smallest D(I)
1029                //
1030                isub = 1;
1031                smin = d[1];
1032                for(j=2; j<=n+1-i; j++)
1033                {
1034                    if( (double)(d[j])<=(double)(smin) )
1035                    {
1036                        isub = j;
1037                        smin = d[j];
1038                    }
1039                }
1040                if( isub!=n+1-i )
1041                {
1042                   
1043                    //
1044                    // Swap singular values and vectors
1045                    //
1046                    d[isub] = d[n+1-i];
1047                    d[n+1-i] = smin;
1048                    if( ncvt>0 )
1049                    {
1050                        j = n+1-i;
1051                        for(i_=vstart; i_<=vend;i_++)
1052                        {
1053                            vttemp[i_] = vt[isub+vstart-1,i_];
1054                        }
1055                        for(i_=vstart; i_<=vend;i_++)
1056                        {
1057                            vt[isub+vstart-1,i_] = vt[j+vstart-1,i_];
1058                        }
1059                        for(i_=vstart; i_<=vend;i_++)
1060                        {
1061                            vt[j+vstart-1,i_] = vttemp[i_];
1062                        }
1063                    }
1064                    if( nru>0 )
1065                    {
1066                        j = n+1-i;
1067                        for(i_=ustart; i_<=uend;i_++)
1068                        {
1069                            utemp[i_] = u[i_,isub+ustart-1];
1070                        }
1071                        for(i_=ustart; i_<=uend;i_++)
1072                        {
1073                            u[i_,isub+ustart-1] = u[i_,j+ustart-1];
1074                        }
1075                        for(i_=ustart; i_<=uend;i_++)
1076                        {
1077                            u[i_,j+ustart-1] = utemp[i_];
1078                        }
1079                    }
1080                    if( ncc>0 )
1081                    {
1082                        j = n+1-i;
1083                        for(i_=cstart; i_<=cend;i_++)
1084                        {
1085                            ctemp[i_] = c[isub+cstart-1,i_];
1086                        }
1087                        for(i_=cstart; i_<=cend;i_++)
1088                        {
1089                            c[isub+cstart-1,i_] = c[j+cstart-1,i_];
1090                        }
1091                        for(i_=cstart; i_<=cend;i_++)
1092                        {
1093                            c[j+cstart-1,i_] = ctemp[i_];
1094                        }
1095                    }
1096                }
1097            }
1098            return result;
1099        }
1100
1101
1102        private static double extsignbdsqr(double a,
1103            double b)
1104        {
1105            double result = 0;
1106
1107            if( (double)(b)>=(double)(0) )
1108            {
1109                result = Math.Abs(a);
1110            }
1111            else
1112            {
1113                result = -Math.Abs(a);
1114            }
1115            return result;
1116        }
1117
1118
1119        private static void svd2x2(double f,
1120            double g,
1121            double h,
1122            ref double ssmin,
1123            ref double ssmax)
1124        {
1125            double aas = 0;
1126            double at = 0;
1127            double au = 0;
1128            double c = 0;
1129            double fa = 0;
1130            double fhmn = 0;
1131            double fhmx = 0;
1132            double ga = 0;
1133            double ha = 0;
1134
1135            fa = Math.Abs(f);
1136            ga = Math.Abs(g);
1137            ha = Math.Abs(h);
1138            fhmn = Math.Min(fa, ha);
1139            fhmx = Math.Max(fa, ha);
1140            if( (double)(fhmn)==(double)(0) )
1141            {
1142                ssmin = 0;
1143                if( (double)(fhmx)==(double)(0) )
1144                {
1145                    ssmax = ga;
1146                }
1147                else
1148                {
1149                    ssmax = Math.Max(fhmx, ga)*Math.Sqrt(1+AP.Math.Sqr(Math.Min(fhmx, ga)/Math.Max(fhmx, ga)));
1150                }
1151            }
1152            else
1153            {
1154                if( (double)(ga)<(double)(fhmx) )
1155                {
1156                    aas = 1+fhmn/fhmx;
1157                    at = (fhmx-fhmn)/fhmx;
1158                    au = AP.Math.Sqr(ga/fhmx);
1159                    c = 2/(Math.Sqrt(aas*aas+au)+Math.Sqrt(at*at+au));
1160                    ssmin = fhmn*c;
1161                    ssmax = fhmx/c;
1162                }
1163                else
1164                {
1165                    au = fhmx/ga;
1166                    if( (double)(au)==(double)(0) )
1167                    {
1168                       
1169                        //
1170                        // Avoid possible harmful underflow if exponent range
1171                        // asymmetric (true SSMIN may not underflow even if
1172                        // AU underflows)
1173                        //
1174                        ssmin = fhmn*fhmx/ga;
1175                        ssmax = ga;
1176                    }
1177                    else
1178                    {
1179                        aas = 1+fhmn/fhmx;
1180                        at = (fhmx-fhmn)/fhmx;
1181                        c = 1/(Math.Sqrt(1+AP.Math.Sqr(aas*au))+Math.Sqrt(1+AP.Math.Sqr(at*au)));
1182                        ssmin = fhmn*c*au;
1183                        ssmin = ssmin+ssmin;
1184                        ssmax = ga/(c+c);
1185                    }
1186                }
1187            }
1188        }
1189
1190
1191        private static void svdv2x2(double f,
1192            double g,
1193            double h,
1194            ref double ssmin,
1195            ref double ssmax,
1196            ref double snr,
1197            ref double csr,
1198            ref double snl,
1199            ref double csl)
1200        {
1201            bool gasmal = new bool();
1202            bool swp = new bool();
1203            int pmax = 0;
1204            double a = 0;
1205            double clt = 0;
1206            double crt = 0;
1207            double d = 0;
1208            double fa = 0;
1209            double ft = 0;
1210            double ga = 0;
1211            double gt = 0;
1212            double ha = 0;
1213            double ht = 0;
1214            double l = 0;
1215            double m = 0;
1216            double mm = 0;
1217            double r = 0;
1218            double s = 0;
1219            double slt = 0;
1220            double srt = 0;
1221            double t = 0;
1222            double temp = 0;
1223            double tsign = 0;
1224            double tt = 0;
1225            double v = 0;
1226
1227            ft = f;
1228            fa = Math.Abs(ft);
1229            ht = h;
1230            ha = Math.Abs(h);
1231           
1232            //
1233            // PMAX points to the maximum absolute element of matrix
1234            //  PMAX = 1 if F largest in absolute values
1235            //  PMAX = 2 if G largest in absolute values
1236            //  PMAX = 3 if H largest in absolute values
1237            //
1238            pmax = 1;
1239            swp = (double)(ha)>(double)(fa);
1240            if( swp )
1241            {
1242               
1243                //
1244                // Now FA .ge. HA
1245                //
1246                pmax = 3;
1247                temp = ft;
1248                ft = ht;
1249                ht = temp;
1250                temp = fa;
1251                fa = ha;
1252                ha = temp;
1253            }
1254            gt = g;
1255            ga = Math.Abs(gt);
1256            if( (double)(ga)==(double)(0) )
1257            {
1258               
1259                //
1260                // Diagonal matrix
1261                //
1262                ssmin = ha;
1263                ssmax = fa;
1264                clt = 1;
1265                crt = 1;
1266                slt = 0;
1267                srt = 0;
1268            }
1269            else
1270            {
1271                gasmal = true;
1272                if( (double)(ga)>(double)(fa) )
1273                {
1274                    pmax = 2;
1275                    if( (double)(fa/ga)<(double)(AP.Math.MachineEpsilon) )
1276                    {
1277                       
1278                        //
1279                        // Case of very large GA
1280                        //
1281                        gasmal = false;
1282                        ssmax = ga;
1283                        if( (double)(ha)>(double)(1) )
1284                        {
1285                            v = ga/ha;
1286                            ssmin = fa/v;
1287                        }
1288                        else
1289                        {
1290                            v = fa/ga;
1291                            ssmin = v*ha;
1292                        }
1293                        clt = 1;
1294                        slt = ht/gt;
1295                        srt = 1;
1296                        crt = ft/gt;
1297                    }
1298                }
1299                if( gasmal )
1300                {
1301                   
1302                    //
1303                    // Normal case
1304                    //
1305                    d = fa-ha;
1306                    if( (double)(d)==(double)(fa) )
1307                    {
1308                        l = 1;
1309                    }
1310                    else
1311                    {
1312                        l = d/fa;
1313                    }
1314                    m = gt/ft;
1315                    t = 2-l;
1316                    mm = m*m;
1317                    tt = t*t;
1318                    s = Math.Sqrt(tt+mm);
1319                    if( (double)(l)==(double)(0) )
1320                    {
1321                        r = Math.Abs(m);
1322                    }
1323                    else
1324                    {
1325                        r = Math.Sqrt(l*l+mm);
1326                    }
1327                    a = 0.5*(s+r);
1328                    ssmin = ha/a;
1329                    ssmax = fa*a;
1330                    if( (double)(mm)==(double)(0) )
1331                    {
1332                       
1333                        //
1334                        // Note that M is very tiny
1335                        //
1336                        if( (double)(l)==(double)(0) )
1337                        {
1338                            t = extsignbdsqr(2, ft)*extsignbdsqr(1, gt);
1339                        }
1340                        else
1341                        {
1342                            t = gt/extsignbdsqr(d, ft)+m/t;
1343                        }
1344                    }
1345                    else
1346                    {
1347                        t = (m/(s+t)+m/(r+l))*(1+a);
1348                    }
1349                    l = Math.Sqrt(t*t+4);
1350                    crt = 2/l;
1351                    srt = t/l;
1352                    clt = (crt+srt*m)/a;
1353                    v = ht/ft;
1354                    slt = v*srt/a;
1355                }
1356            }
1357            if( swp )
1358            {
1359                csl = srt;
1360                snl = crt;
1361                csr = slt;
1362                snr = clt;
1363            }
1364            else
1365            {
1366                csl = clt;
1367                snl = slt;
1368                csr = crt;
1369                snr = srt;
1370            }
1371           
1372            //
1373            // Correct signs of SSMAX and SSMIN
1374            //
1375            if( pmax==1 )
1376            {
1377                tsign = extsignbdsqr(1, csr)*extsignbdsqr(1, csl)*extsignbdsqr(1, f);
1378            }
1379            if( pmax==2 )
1380            {
1381                tsign = extsignbdsqr(1, snr)*extsignbdsqr(1, csl)*extsignbdsqr(1, g);
1382            }
1383            if( pmax==3 )
1384            {
1385                tsign = extsignbdsqr(1, snr)*extsignbdsqr(1, snl)*extsignbdsqr(1, h);
1386            }
1387            ssmax = extsignbdsqr(ssmax, tsign);
1388            ssmin = extsignbdsqr(ssmin, tsign*extsignbdsqr(1, f)*extsignbdsqr(1, h));
1389        }
1390    }
1391}
Note: See TracBrowser for help on using the repository browser.