Free cookie consent management tool by TermsFeed Policy Generator

source: branches/3.2/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.1.2.2591/ALGLIB-2.1.2.2591/tdevd.cs @ 3580

Last change on this file since 3580 was 2645, checked in by mkommend, 15 years ago

extracted external libraries and adapted dependent plugins (ticket #837)
File size: 43.9 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 tdevd
32    {
33        /*************************************************************************
34        Finding the eigenvalues and eigenvectors of a tridiagonal symmetric matrix
35
36        The algorithm finds the eigen pairs of a tridiagonal symmetric matrix by
37        using an QL/QR algorithm with implicit shifts.
38
39        Input parameters:
40            D       -   the main diagonal of a tridiagonal matrix.
41                        Array whose index ranges within [0..N-1].
42            E       -   the secondary diagonal of a tridiagonal matrix.
43                        Array whose index ranges within [0..N-2].
44            N       -   size of matrix A.
45            ZNeeded -   flag controlling whether the eigenvectors are needed or not.
46                        If ZNeeded is equal to:
47                         * 0, the eigenvectors are not needed;
48                         * 1, the eigenvectors of a tridiagonal matrix
49                           are multiplied by the square matrix Z. It is used if the
50                           tridiagonal matrix is obtained by the similarity
51                           transformation of a symmetric matrix;
52                         * 2, the eigenvectors of a tridiagonal matrix replace the
53                           square matrix Z;
54                         * 3, matrix Z contains the first row of the eigenvectors
55                           matrix.
56            Z       -   if ZNeeded=1, Z contains the square matrix by which the
57                        eigenvectors are multiplied.
58                        Array whose indexes range within [0..N-1, 0..N-1].
59
60        Output parameters:
61            D       -   eigenvalues in ascending order.
62                        Array whose index ranges within [0..N-1].
63            Z       -   if ZNeeded is equal to:
64                         * 0, Z hasn’t changed;
65                         * 1, Z contains the product of a given matrix (from the left)
66                           and the eigenvectors matrix (from the right);
67                         * 2, Z contains the eigenvectors.
68                         * 3, Z contains the first row of the eigenvectors matrix.
69                        If ZNeeded<3, Z is the array whose indexes range within [0..N-1, 0..N-1].
70                        In that case, the eigenvectors are stored in the matrix columns.
71                        If ZNeeded=3, Z is the array whose indexes range within [0..0, 0..N-1].
72
73        Result:
74            True, if the algorithm has converged.
75            False, if the algorithm hasn't converged.
76
77          -- LAPACK routine (version 3.0) --
78             Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
79             Courant Institute, Argonne National Lab, and Rice University
80             September 30, 1994
81        *************************************************************************/
82        public static bool smatrixtdevd(ref double[] d,
83            double[] e,
84            int n,
85            int zneeded,
86            ref double[,] z)
87        {
88            bool result = new bool();
89            double[] d1 = new double[0];
90            double[] e1 = new double[0];
91            double[,] z1 = new double[0,0];
92            int i = 0;
93            int i_ = 0;
94            int i1_ = 0;
95
96            e = (double[])e.Clone();
97
98           
99            //
100            // Prepare 1-based task
101            //
102            d1 = new double[n+1];
103            e1 = new double[n+1];
104            i1_ = (0) - (1);
105            for(i_=1; i_<=n;i_++)
106            {
107                d1[i_] = d[i_+i1_];
108            }
109            if( n>1 )
110            {
111                i1_ = (0) - (1);
112                for(i_=1; i_<=n-1;i_++)
113                {
114                    e1[i_] = e[i_+i1_];
115                }
116            }
117            if( zneeded==1 )
118            {
119                z1 = new double[n+1, n+1];
120                for(i=1; i<=n; i++)
121                {
122                    i1_ = (0) - (1);
123                    for(i_=1; i_<=n;i_++)
124                    {
125                        z1[i,i_] = z[i-1,i_+i1_];
126                    }
127                }
128            }
129           
130            //
131            // Solve 1-based task
132            //
133            result = tridiagonalevd(ref d1, e1, n, zneeded, ref z1);
134            if( !result )
135            {
136                return result;
137            }
138           
139            //
140            // Convert back to 0-based result
141            //
142            i1_ = (1) - (0);
143            for(i_=0; i_<=n-1;i_++)
144            {
145                d[i_] = d1[i_+i1_];
146            }
147            if( zneeded!=0 )
148            {
149                if( zneeded==1 )
150                {
151                    for(i=1; i<=n; i++)
152                    {
153                        i1_ = (1) - (0);
154                        for(i_=0; i_<=n-1;i_++)
155                        {
156                            z[i-1,i_] = z1[i,i_+i1_];
157                        }
158                    }
159                    return result;
160                }
161                if( zneeded==2 )
162                {
163                    z = new double[n-1+1, n-1+1];
164                    for(i=1; i<=n; i++)
165                    {
166                        i1_ = (1) - (0);
167                        for(i_=0; i_<=n-1;i_++)
168                        {
169                            z[i-1,i_] = z1[i,i_+i1_];
170                        }
171                    }
172                    return result;
173                }
174                if( zneeded==3 )
175                {
176                    z = new double[0+1, n-1+1];
177                    i1_ = (1) - (0);
178                    for(i_=0; i_<=n-1;i_++)
179                    {
180                        z[0,i_] = z1[1,i_+i1_];
181                    }
182                    return result;
183                }
184                System.Diagnostics.Debug.Assert(false, "SMatrixTDEVD: Incorrect ZNeeded!");
185            }
186            return result;
187        }
188
189
190        public static bool tridiagonalevd(ref double[] d,
191            double[] e,
192            int n,
193            int zneeded,
194            ref double[,] z)
195        {
196            bool result = new bool();
197            int maxit = 0;
198            int i = 0;
199            int icompz = 0;
200            int ii = 0;
201            int iscale = 0;
202            int j = 0;
203            int jtot = 0;
204            int k = 0;
205            int t = 0;
206            int l = 0;
207            int l1 = 0;
208            int lend = 0;
209            int lendm1 = 0;
210            int lendp1 = 0;
211            int lendsv = 0;
212            int lm1 = 0;
213            int lsv = 0;
214            int m = 0;
215            int mm = 0;
216            int mm1 = 0;
217            int nm1 = 0;
218            int nmaxit = 0;
219            int tmpint = 0;
220            double anorm = 0;
221            double b = 0;
222            double c = 0;
223            double eps = 0;
224            double eps2 = 0;
225            double f = 0;
226            double g = 0;
227            double p = 0;
228            double r = 0;
229            double rt1 = 0;
230            double rt2 = 0;
231            double s = 0;
232            double safmax = 0;
233            double safmin = 0;
234            double ssfmax = 0;
235            double ssfmin = 0;
236            double tst = 0;
237            double tmp = 0;
238            double[] work1 = new double[0];
239            double[] work2 = new double[0];
240            double[] workc = new double[0];
241            double[] works = new double[0];
242            double[] wtemp = new double[0];
243            bool gotoflag = new bool();
244            int zrows = 0;
245            bool wastranspose = new bool();
246            int i_ = 0;
247
248            e = (double[])e.Clone();
249
250            System.Diagnostics.Debug.Assert(zneeded>=0 & zneeded<=3, "TridiagonalEVD: Incorrent ZNeeded");
251           
252            //
253            // Quick return if possible
254            //
255            if( zneeded<0 | zneeded>3 )
256            {
257                result = false;
258                return result;
259            }
260            result = true;
261            if( n==0 )
262            {
263                return result;
264            }
265            if( n==1 )
266            {
267                if( zneeded==2 | zneeded==3 )
268                {
269                    z = new double[1+1, 1+1];
270                    z[1,1] = 1;
271                }
272                return result;
273            }
274            maxit = 30;
275           
276            //
277            // Initialize arrays
278            //
279            wtemp = new double[n+1];
280            work1 = new double[n-1+1];
281            work2 = new double[n-1+1];
282            workc = new double[n+1];
283            works = new double[n+1];
284           
285            //
286            // Determine the unit roundoff and over/underflow thresholds.
287            //
288            eps = AP.Math.MachineEpsilon;
289            eps2 = AP.Math.Sqr(eps);
290            safmin = AP.Math.MinRealNumber;
291            safmax = AP.Math.MaxRealNumber;
292            ssfmax = Math.Sqrt(safmax)/3;
293            ssfmin = Math.Sqrt(safmin)/eps2;
294           
295            //
296            // Prepare Z
297            //
298            // Here we are using transposition to get rid of column operations
299            //
300            //
301            wastranspose = false;
302            if( zneeded==0 )
303            {
304                zrows = 0;
305            }
306            if( zneeded==1 )
307            {
308                zrows = n;
309            }
310            if( zneeded==2 )
311            {
312                zrows = n;
313            }
314            if( zneeded==3 )
315            {
316                zrows = 1;
317            }
318            if( zneeded==1 )
319            {
320                wastranspose = true;
321                blas.inplacetranspose(ref z, 1, n, 1, n, ref wtemp);
322            }
323            if( zneeded==2 )
324            {
325                wastranspose = true;
326                z = new double[n+1, n+1];
327                for(i=1; i<=n; i++)
328                {
329                    for(j=1; j<=n; j++)
330                    {
331                        if( i==j )
332                        {
333                            z[i,j] = 1;
334                        }
335                        else
336                        {
337                            z[i,j] = 0;
338                        }
339                    }
340                }
341            }
342            if( zneeded==3 )
343            {
344                wastranspose = false;
345                z = new double[1+1, n+1];
346                for(j=1; j<=n; j++)
347                {
348                    if( j==1 )
349                    {
350                        z[1,j] = 1;
351                    }
352                    else
353                    {
354                        z[1,j] = 0;
355                    }
356                }
357            }
358            nmaxit = n*maxit;
359            jtot = 0;
360           
361            //
362            // Determine where the matrix splits and choose QL or QR iteration
363            // for each block, according to whether top or bottom diagonal
364            // element is smaller.
365            //
366            l1 = 1;
367            nm1 = n-1;
368            while( true )
369            {
370                if( l1>n )
371                {
372                    break;
373                }
374                if( l1>1 )
375                {
376                    e[l1-1] = 0;
377                }
378                gotoflag = false;
379                if( l1<=nm1 )
380                {
381                    for(m=l1; m<=nm1; m++)
382                    {
383                        tst = Math.Abs(e[m]);
384                        if( (double)(tst)==(double)(0) )
385                        {
386                            gotoflag = true;
387                            break;
388                        }
389                        if( (double)(tst)<=(double)(Math.Sqrt(Math.Abs(d[m]))*Math.Sqrt(Math.Abs(d[m+1]))*eps) )
390                        {
391                            e[m] = 0;
392                            gotoflag = true;
393                            break;
394                        }
395                    }
396                }
397                if( !gotoflag )
398                {
399                    m = n;
400                }
401               
402                //
403                // label 30:
404                //
405                l = l1;
406                lsv = l;
407                lend = m;
408                lendsv = lend;
409                l1 = m+1;
410                if( lend==l )
411                {
412                    continue;
413                }
414               
415                //
416                // Scale submatrix in rows and columns L to LEND
417                //
418                if( l==lend )
419                {
420                    anorm = Math.Abs(d[l]);
421                }
422                else
423                {
424                    anorm = Math.Max(Math.Abs(d[l])+Math.Abs(e[l]), Math.Abs(e[lend-1])+Math.Abs(d[lend]));
425                    for(i=l+1; i<=lend-1; i++)
426                    {
427                        anorm = Math.Max(anorm, Math.Abs(d[i])+Math.Abs(e[i])+Math.Abs(e[i-1]));
428                    }
429                }
430                iscale = 0;
431                if( (double)(anorm)==(double)(0) )
432                {
433                    continue;
434                }
435                if( (double)(anorm)>(double)(ssfmax) )
436                {
437                    iscale = 1;
438                    tmp = ssfmax/anorm;
439                    tmpint = lend-1;
440                    for(i_=l; i_<=lend;i_++)
441                    {
442                        d[i_] = tmp*d[i_];
443                    }
444                    for(i_=l; i_<=tmpint;i_++)
445                    {
446                        e[i_] = tmp*e[i_];
447                    }
448                }
449                if( (double)(anorm)<(double)(ssfmin) )
450                {
451                    iscale = 2;
452                    tmp = ssfmin/anorm;
453                    tmpint = lend-1;
454                    for(i_=l; i_<=lend;i_++)
455                    {
456                        d[i_] = tmp*d[i_];
457                    }
458                    for(i_=l; i_<=tmpint;i_++)
459                    {
460                        e[i_] = tmp*e[i_];
461                    }
462                }
463               
464                //
465                // Choose between QL and QR iteration
466                //
467                if( (double)(Math.Abs(d[lend]))<(double)(Math.Abs(d[l])) )
468                {
469                    lend = lsv;
470                    l = lendsv;
471                }
472                if( lend>l )
473                {
474                   
475                    //
476                    // QL Iteration
477                    //
478                    // Look for small subdiagonal element.
479                    //
480                    while( true )
481                    {
482                        gotoflag = false;
483                        if( l!=lend )
484                        {
485                            lendm1 = lend-1;
486                            for(m=l; m<=lendm1; m++)
487                            {
488                                tst = AP.Math.Sqr(Math.Abs(e[m]));
489                                if( (double)(tst)<=(double)(eps2*Math.Abs(d[m])*Math.Abs(d[m+1])+safmin) )
490                                {
491                                    gotoflag = true;
492                                    break;
493                                }
494                            }
495                        }
496                        if( !gotoflag )
497                        {
498                            m = lend;
499                        }
500                        if( m<lend )
501                        {
502                            e[m] = 0;
503                        }
504                        p = d[l];
505                        if( m!=l )
506                        {
507                           
508                            //
509                            // If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
510                            // to compute its eigensystem.
511                            //
512                            if( m==l+1 )
513                            {
514                                if( zneeded>0 )
515                                {
516                                    tdevdev2(d[l], e[l], d[l+1], ref rt1, ref rt2, ref c, ref s);
517                                    work1[l] = c;
518                                    work2[l] = s;
519                                    workc[1] = work1[l];
520                                    works[1] = work2[l];
521                                    if( !wastranspose )
522                                    {
523                                        rotations.applyrotationsfromtheright(false, 1, zrows, l, l+1, ref workc, ref works, ref z, ref wtemp);
524                                    }
525                                    else
526                                    {
527                                        rotations.applyrotationsfromtheleft(false, l, l+1, 1, zrows, ref workc, ref works, ref z, ref wtemp);
528                                    }
529                                }
530                                else
531                                {
532                                    tdevde2(d[l], e[l], d[l+1], ref rt1, ref rt2);
533                                }
534                                d[l] = rt1;
535                                d[l+1] = rt2;
536                                e[l] = 0;
537                                l = l+2;
538                                if( l<=lend )
539                                {
540                                    continue;
541                                }
542                               
543                                //
544                                // GOTO 140
545                                //
546                                break;
547                            }
548                            if( jtot==nmaxit )
549                            {
550                               
551                                //
552                                // GOTO 140
553                                //
554                                break;
555                            }
556                            jtot = jtot+1;
557                           
558                            //
559                            // Form shift.
560                            //
561                            g = (d[l+1]-p)/(2*e[l]);
562                            r = tdevdpythag(g, 1);
563                            g = d[m]-p+e[l]/(g+tdevdextsign(r, g));
564                            s = 1;
565                            c = 1;
566                            p = 0;
567                           
568                            //
569                            // Inner loop
570                            //
571                            mm1 = m-1;
572                            for(i=mm1; i>=l; i--)
573                            {
574                                f = s*e[i];
575                                b = c*e[i];
576                                rotations.generaterotation(g, f, ref c, ref s, ref r);
577                                if( i!=m-1 )
578                                {
579                                    e[i+1] = r;
580                                }
581                                g = d[i+1]-p;
582                                r = (d[i]-g)*s+2*c*b;
583                                p = s*r;
584                                d[i+1] = g+p;
585                                g = c*r-b;
586                               
587                                //
588                                // If eigenvectors are desired, then save rotations.
589                                //
590                                if( zneeded>0 )
591                                {
592                                    work1[i] = c;
593                                    work2[i] = -s;
594                                }
595                            }
596                           
597                            //
598                            // If eigenvectors are desired, then apply saved rotations.
599                            //
600                            if( zneeded>0 )
601                            {
602                                for(i=l; i<=m-1; i++)
603                                {
604                                    workc[i-l+1] = work1[i];
605                                    works[i-l+1] = work2[i];
606                                }
607                                if( !wastranspose )
608                                {
609                                    rotations.applyrotationsfromtheright(false, 1, zrows, l, m, ref workc, ref works, ref z, ref wtemp);
610                                }
611                                else
612                                {
613                                    rotations.applyrotationsfromtheleft(false, l, m, 1, zrows, ref workc, ref works, ref z, ref wtemp);
614                                }
615                            }
616                            d[l] = d[l]-p;
617                            e[l] = g;
618                            continue;
619                        }
620                       
621                        //
622                        // Eigenvalue found.
623                        //
624                        d[l] = p;
625                        l = l+1;
626                        if( l<=lend )
627                        {
628                            continue;
629                        }
630                        break;
631                    }
632                }
633                else
634                {
635                   
636                    //
637                    // QR Iteration
638                    //
639                    // Look for small superdiagonal element.
640                    //
641                    while( true )
642                    {
643                        gotoflag = false;
644                        if( l!=lend )
645                        {
646                            lendp1 = lend+1;
647                            for(m=l; m>=lendp1; m--)
648                            {
649                                tst = AP.Math.Sqr(Math.Abs(e[m-1]));
650                                if( (double)(tst)<=(double)(eps2*Math.Abs(d[m])*Math.Abs(d[m-1])+safmin) )
651                                {
652                                    gotoflag = true;
653                                    break;
654                                }
655                            }
656                        }
657                        if( !gotoflag )
658                        {
659                            m = lend;
660                        }
661                        if( m>lend )
662                        {
663                            e[m-1] = 0;
664                        }
665                        p = d[l];
666                        if( m!=l )
667                        {
668                           
669                            //
670                            // If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
671                            // to compute its eigensystem.
672                            //
673                            if( m==l-1 )
674                            {
675                                if( zneeded>0 )
676                                {
677                                    tdevdev2(d[l-1], e[l-1], d[l], ref rt1, ref rt2, ref c, ref s);
678                                    work1[m] = c;
679                                    work2[m] = s;
680                                    workc[1] = c;
681                                    works[1] = s;
682                                    if( !wastranspose )
683                                    {
684                                        rotations.applyrotationsfromtheright(true, 1, zrows, l-1, l, ref workc, ref works, ref z, ref wtemp);
685                                    }
686                                    else
687                                    {
688                                        rotations.applyrotationsfromtheleft(true, l-1, l, 1, zrows, ref workc, ref works, ref z, ref wtemp);
689                                    }
690                                }
691                                else
692                                {
693                                    tdevde2(d[l-1], e[l-1], d[l], ref rt1, ref rt2);
694                                }
695                                d[l-1] = rt1;
696                                d[l] = rt2;
697                                e[l-1] = 0;
698                                l = l-2;
699                                if( l>=lend )
700                                {
701                                    continue;
702                                }
703                                break;
704                            }
705                            if( jtot==nmaxit )
706                            {
707                                break;
708                            }
709                            jtot = jtot+1;
710                           
711                            //
712                            // Form shift.
713                            //
714                            g = (d[l-1]-p)/(2*e[l-1]);
715                            r = tdevdpythag(g, 1);
716                            g = d[m]-p+e[l-1]/(g+tdevdextsign(r, g));
717                            s = 1;
718                            c = 1;
719                            p = 0;
720                           
721                            //
722                            // Inner loop
723                            //
724                            lm1 = l-1;
725                            for(i=m; i<=lm1; i++)
726                            {
727                                f = s*e[i];
728                                b = c*e[i];
729                                rotations.generaterotation(g, f, ref c, ref s, ref r);
730                                if( i!=m )
731                                {
732                                    e[i-1] = r;
733                                }
734                                g = d[i]-p;
735                                r = (d[i+1]-g)*s+2*c*b;
736                                p = s*r;
737                                d[i] = g+p;
738                                g = c*r-b;
739                               
740                                //
741                                // If eigenvectors are desired, then save rotations.
742                                //
743                                if( zneeded>0 )
744                                {
745                                    work1[i] = c;
746                                    work2[i] = s;
747                                }
748                            }
749                           
750                            //
751                            // If eigenvectors are desired, then apply saved rotations.
752                            //
753                            if( zneeded>0 )
754                            {
755                                mm = l-m+1;
756                                for(i=m; i<=l-1; i++)
757                                {
758                                    workc[i-m+1] = work1[i];
759                                    works[i-m+1] = work2[i];
760                                }
761                                if( !wastranspose )
762                                {
763                                    rotations.applyrotationsfromtheright(true, 1, zrows, m, l, ref workc, ref works, ref z, ref wtemp);
764                                }
765                                else
766                                {
767                                    rotations.applyrotationsfromtheleft(true, m, l, 1, zrows, ref workc, ref works, ref z, ref wtemp);
768                                }
769                            }
770                            d[l] = d[l]-p;
771                            e[lm1] = g;
772                            continue;
773                        }
774                       
775                        //
776                        // Eigenvalue found.
777                        //
778                        d[l] = p;
779                        l = l-1;
780                        if( l>=lend )
781                        {
782                            continue;
783                        }
784                        break;
785                    }
786                }
787               
788                //
789                // Undo scaling if necessary
790                //
791                if( iscale==1 )
792                {
793                    tmp = anorm/ssfmax;
794                    tmpint = lendsv-1;
795                    for(i_=lsv; i_<=lendsv;i_++)
796                    {
797                        d[i_] = tmp*d[i_];
798                    }
799                    for(i_=lsv; i_<=tmpint;i_++)
800                    {
801                        e[i_] = tmp*e[i_];
802                    }
803                }
804                if( iscale==2 )
805                {
806                    tmp = anorm/ssfmin;
807                    tmpint = lendsv-1;
808                    for(i_=lsv; i_<=lendsv;i_++)
809                    {
810                        d[i_] = tmp*d[i_];
811                    }
812                    for(i_=lsv; i_<=tmpint;i_++)
813                    {
814                        e[i_] = tmp*e[i_];
815                    }
816                }
817               
818                //
819                // Check for no convergence to an eigenvalue after a total
820                // of N*MAXIT iterations.
821                //
822                if( jtot>=nmaxit )
823                {
824                    result = false;
825                    if( wastranspose )
826                    {
827                        blas.inplacetranspose(ref z, 1, n, 1, n, ref wtemp);
828                    }
829                    return result;
830                }
831            }
832           
833            //
834            // Order eigenvalues and eigenvectors.
835            //
836            if( zneeded==0 )
837            {
838               
839                //
840                // Sort
841                //
842                if( n==1 )
843                {
844                    return result;
845                }
846                if( n==2 )
847                {
848                    if( (double)(d[1])>(double)(d[2]) )
849                    {
850                        tmp = d[1];
851                        d[1] = d[2];
852                        d[2] = tmp;
853                    }
854                    return result;
855                }
856                i = 2;
857                do
858                {
859                    t = i;
860                    while( t!=1 )
861                    {
862                        k = t/2;
863                        if( (double)(d[k])>=(double)(d[t]) )
864                        {
865                            t = 1;
866                        }
867                        else
868                        {
869                            tmp = d[k];
870                            d[k] = d[t];
871                            d[t] = tmp;
872                            t = k;
873                        }
874                    }
875                    i = i+1;
876                }
877                while( i<=n );
878                i = n-1;
879                do
880                {
881                    tmp = d[i+1];
882                    d[i+1] = d[1];
883                    d[+1] = tmp;
884                    t = 1;
885                    while( t!=0 )
886                    {
887                        k = 2*t;
888                        if( k>i )
889                        {
890                            t = 0;
891                        }
892                        else
893                        {
894                            if( k<i )
895                            {
896                                if( (double)(d[k+1])>(double)(d[k]) )
897                                {
898                                    k = k+1;
899                                }
900                            }
901                            if( (double)(d[t])>=(double)(d[k]) )
902                            {
903                                t = 0;
904                            }
905                            else
906                            {
907                                tmp = d[k];
908                                d[k] = d[t];
909                                d[t] = tmp;
910                                t = k;
911                            }
912                        }
913                    }
914                    i = i-1;
915                }
916                while( i>=1 );
917            }
918            else
919            {
920               
921                //
922                // Use Selection Sort to minimize swaps of eigenvectors
923                //
924                for(ii=2; ii<=n; ii++)
925                {
926                    i = ii-1;
927                    k = i;
928                    p = d[i];
929                    for(j=ii; j<=n; j++)
930                    {
931                        if( (double)(d[j])<(double)(p) )
932                        {
933                            k = j;
934                            p = d[j];
935                        }
936                    }
937                    if( k!=i )
938                    {
939                        d[k] = d[i];
940                        d[i] = p;
941                        if( wastranspose )
942                        {
943                            for(i_=1; i_<=n;i_++)
944                            {
945                                wtemp[i_] = z[i,i_];
946                            }
947                            for(i_=1; i_<=n;i_++)
948                            {
949                                z[i,i_] = z[k,i_];
950                            }
951                            for(i_=1; i_<=n;i_++)
952                            {
953                                z[k,i_] = wtemp[i_];
954                            }
955                        }
956                        else
957                        {
958                            for(i_=1; i_<=zrows;i_++)
959                            {
960                                wtemp[i_] = z[i_,i];
961                            }
962                            for(i_=1; i_<=zrows;i_++)
963                            {
964                                z[i_,i] = z[i_,k];
965                            }
966                            for(i_=1; i_<=zrows;i_++)
967                            {
968                                z[i_,k] = wtemp[i_];
969                            }
970                        }
971                    }
972                }
973                if( wastranspose )
974                {
975                    blas.inplacetranspose(ref z, 1, n, 1, n, ref wtemp);
976                }
977            }
978            return result;
979        }
980
981
982        /*************************************************************************
983        DLAE2  computes the eigenvalues of a 2-by-2 symmetric matrix
984           [  A   B  ]
985           [  B   C  ].
986        On return, RT1 is the eigenvalue of larger absolute value, and RT2
987        is the eigenvalue of smaller absolute value.
988
989          -- LAPACK auxiliary routine (version 3.0) --
990             Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
991             Courant Institute, Argonne National Lab, and Rice University
992             October 31, 1992
993        *************************************************************************/
994        private static void tdevde2(double a,
995            double b,
996            double c,
997            ref double rt1,
998            ref double rt2)
999        {
1000            double ab = 0;
1001            double acmn = 0;
1002            double acmx = 0;
1003            double adf = 0;
1004            double df = 0;
1005            double rt = 0;
1006            double sm = 0;
1007            double tb = 0;
1008
1009            sm = a+c;
1010            df = a-c;
1011            adf = Math.Abs(df);
1012            tb = b+b;
1013            ab = Math.Abs(tb);
1014            if( (double)(Math.Abs(a))>(double)(Math.Abs(c)) )
1015            {
1016                acmx = a;
1017                acmn = c;
1018            }
1019            else
1020            {
1021                acmx = c;
1022                acmn = a;
1023            }
1024            if( (double)(adf)>(double)(ab) )
1025            {
1026                rt = adf*Math.Sqrt(1+AP.Math.Sqr(ab/adf));
1027            }
1028            else
1029            {
1030                if( (double)(adf)<(double)(ab) )
1031                {
1032                    rt = ab*Math.Sqrt(1+AP.Math.Sqr(adf/ab));
1033                }
1034                else
1035                {
1036                   
1037                    //
1038                    // Includes case AB=ADF=0
1039                    //
1040                    rt = ab*Math.Sqrt(2);
1041                }
1042            }
1043            if( (double)(sm)<(double)(0) )
1044            {
1045                rt1 = 0.5*(sm-rt);
1046               
1047                //
1048                // Order of execution important.
1049                // To get fully accurate smaller eigenvalue,
1050                // next line needs to be executed in higher precision.
1051                //
1052                rt2 = acmx/rt1*acmn-b/rt1*b;
1053            }
1054            else
1055            {
1056                if( (double)(sm)>(double)(0) )
1057                {
1058                    rt1 = 0.5*(sm+rt);
1059                   
1060                    //
1061                    // Order of execution important.
1062                    // To get fully accurate smaller eigenvalue,
1063                    // next line needs to be executed in higher precision.
1064                    //
1065                    rt2 = acmx/rt1*acmn-b/rt1*b;
1066                }
1067                else
1068                {
1069                   
1070                    //
1071                    // Includes case RT1 = RT2 = 0
1072                    //
1073                    rt1 = 0.5*rt;
1074                    rt2 = -(0.5*rt);
1075                }
1076            }
1077        }
1078
1079
1080        /*************************************************************************
1081        DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix
1082
1083           [  A   B  ]
1084           [  B   C  ].
1085
1086        On return, RT1 is the eigenvalue of larger absolute value, RT2 is the
1087        eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right
1088        eigenvector for RT1, giving the decomposition
1089
1090           [ CS1  SN1 ] [  A   B  ] [ CS1 -SN1 ]  =  [ RT1  0  ]
1091           [-SN1  CS1 ] [  B   C  ] [ SN1  CS1 ]     [  0  RT2 ].
1092
1093
1094          -- LAPACK auxiliary routine (version 3.0) --
1095             Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1096             Courant Institute, Argonne National Lab, and Rice University
1097             October 31, 1992
1098        *************************************************************************/
1099        private static void tdevdev2(double a,
1100            double b,
1101            double c,
1102            ref double rt1,
1103            ref double rt2,
1104            ref double cs1,
1105            ref double sn1)
1106        {
1107            int sgn1 = 0;
1108            int sgn2 = 0;
1109            double ab = 0;
1110            double acmn = 0;
1111            double acmx = 0;
1112            double acs = 0;
1113            double adf = 0;
1114            double cs = 0;
1115            double ct = 0;
1116            double df = 0;
1117            double rt = 0;
1118            double sm = 0;
1119            double tb = 0;
1120            double tn = 0;
1121
1122           
1123            //
1124            // Compute the eigenvalues
1125            //
1126            sm = a+c;
1127            df = a-c;
1128            adf = Math.Abs(df);
1129            tb = b+b;
1130            ab = Math.Abs(tb);
1131            if( (double)(Math.Abs(a))>(double)(Math.Abs(c)) )
1132            {
1133                acmx = a;
1134                acmn = c;
1135            }
1136            else
1137            {
1138                acmx = c;
1139                acmn = a;
1140            }
1141            if( (double)(adf)>(double)(ab) )
1142            {
1143                rt = adf*Math.Sqrt(1+AP.Math.Sqr(ab/adf));
1144            }
1145            else
1146            {
1147                if( (double)(adf)<(double)(ab) )
1148                {
1149                    rt = ab*Math.Sqrt(1+AP.Math.Sqr(adf/ab));
1150                }
1151                else
1152                {
1153                   
1154                    //
1155                    // Includes case AB=ADF=0
1156                    //
1157                    rt = ab*Math.Sqrt(2);
1158                }
1159            }
1160            if( (double)(sm)<(double)(0) )
1161            {
1162                rt1 = 0.5*(sm-rt);
1163                sgn1 = -1;
1164               
1165                //
1166                // Order of execution important.
1167                // To get fully accurate smaller eigenvalue,
1168                // next line needs to be executed in higher precision.
1169                //
1170                rt2 = acmx/rt1*acmn-b/rt1*b;
1171            }
1172            else
1173            {
1174                if( (double)(sm)>(double)(0) )
1175                {
1176                    rt1 = 0.5*(sm+rt);
1177                    sgn1 = 1;
1178                   
1179                    //
1180                    // Order of execution important.
1181                    // To get fully accurate smaller eigenvalue,
1182                    // next line needs to be executed in higher precision.
1183                    //
1184                    rt2 = acmx/rt1*acmn-b/rt1*b;
1185                }
1186                else
1187                {
1188                   
1189                    //
1190                    // Includes case RT1 = RT2 = 0
1191                    //
1192                    rt1 = 0.5*rt;
1193                    rt2 = -(0.5*rt);
1194                    sgn1 = 1;
1195                }
1196            }
1197           
1198            //
1199            // Compute the eigenvector
1200            //
1201            if( (double)(df)>=(double)(0) )
1202            {
1203                cs = df+rt;
1204                sgn2 = 1;
1205            }
1206            else
1207            {
1208                cs = df-rt;
1209                sgn2 = -1;
1210            }
1211            acs = Math.Abs(cs);
1212            if( (double)(acs)>(double)(ab) )
1213            {
1214                ct = -(tb/cs);
1215                sn1 = 1/Math.Sqrt(1+ct*ct);
1216                cs1 = ct*sn1;
1217            }
1218            else
1219            {
1220                if( (double)(ab)==(double)(0) )
1221                {
1222                    cs1 = 1;
1223                    sn1 = 0;
1224                }
1225                else
1226                {
1227                    tn = -(cs/tb);
1228                    cs1 = 1/Math.Sqrt(1+tn*tn);
1229                    sn1 = tn*cs1;
1230                }
1231            }
1232            if( sgn1==sgn2 )
1233            {
1234                tn = cs1;
1235                cs1 = -sn1;
1236                sn1 = tn;
1237            }
1238        }
1239
1240
1241        /*************************************************************************
1242        Internal routine
1243        *************************************************************************/
1244        private static double tdevdpythag(double a,
1245            double b)
1246        {
1247            double result = 0;
1248
1249            if( (double)(Math.Abs(a))<(double)(Math.Abs(b)) )
1250            {
1251                result = Math.Abs(b)*Math.Sqrt(1+AP.Math.Sqr(a/b));
1252            }
1253            else
1254            {
1255                result = Math.Abs(a)*Math.Sqrt(1+AP.Math.Sqr(b/a));
1256            }
1257            return result;
1258        }
1259
1260
1261        /*************************************************************************
1262        Internal routine
1263        *************************************************************************/
1264        private static double tdevdextsign(double a,
1265            double b)
1266        {
1267            double result = 0;
1268
1269            if( (double)(b)>=(double)(0) )
1270            {
1271                result = Math.Abs(a);
1272            }
1273            else
1274            {
1275                result = -Math.Abs(a);
1276            }
1277            return result;
1278        }
1279    }
1280}
Note: See TracBrowser for help on using the repository browser.