Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/ALGLIB/hsschur.cs @ 2575

Last change on this file since 2575 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: 48.7 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 hsschur
32    {
33        /*************************************************************************
34        Subroutine performing  the  Schur  decomposition  of  a  matrix  in  upper
35        Hessenberg form using the QR algorithm with multiple shifts.
36
37        The  source matrix  H  is  represented as  S'*H*S = T, where H - matrix in
38        upper Hessenberg form,  S - orthogonal matrix (Schur vectors),   T - upper
39        quasi-triangular matrix (with blocks of sizes  1x1  and  2x2  on  the main
40        diagonal).
41
42        Input parameters:
43            H   -   matrix to be decomposed.
44                    Array whose indexes range within [1..N, 1..N].
45            N   -   size of H, N>=0.
46
47
48        Output parameters:
49            H   –   contains the matrix T.
50                    Array whose indexes range within [1..N, 1..N].
51                    All elements below the blocks on the main diagonal are equal
52                    to 0.
53            S   -   contains Schur vectors.
54                    Array whose indexes range within [1..N, 1..N].
55
56        Note 1:
57            The block structure of matrix T could be easily recognized: since  all
58            the elements  below  the blocks are zeros, the elements a[i+1,i] which
59            are equal to 0 show the block border.
60
61        Note 2:
62            the algorithm  performance  depends  on  the  value  of  the  internal
63            parameter NS of InternalSchurDecomposition  subroutine  which  defines
64            the number of shifts in the QR algorithm (analog of  the  block  width
65            in block matrix algorithms in linear algebra). If you require  maximum
66            performance  on  your  machine,  it  is  recommended  to  adjust  this
67            parameter manually.
68
69        Result:
70            True, if the algorithm has converged and the parameters H and S contain
71                the result.
72            False, if the algorithm has not converged.
73
74        Algorithm implemented on the basis of subroutine DHSEQR (LAPACK 3.0 library).
75        *************************************************************************/
76        public static bool upperhessenbergschurdecomposition(ref double[,] h,
77            int n,
78            ref double[,] s)
79        {
80            bool result = new bool();
81            double[] wi = new double[0];
82            double[] wr = new double[0];
83            int info = 0;
84
85            internalschurdecomposition(ref h, n, 1, 2, ref wr, ref wi, ref s, ref info);
86            result = info==0;
87            return result;
88        }
89
90
91        public static void internalschurdecomposition(ref double[,] h,
92            int n,
93            int tneeded,
94            int zneeded,
95            ref double[] wr,
96            ref double[] wi,
97            ref double[,] z,
98            ref int info)
99        {
100            double[] work = new double[0];
101            int i = 0;
102            int i1 = 0;
103            int i2 = 0;
104            int ierr = 0;
105            int ii = 0;
106            int itemp = 0;
107            int itn = 0;
108            int its = 0;
109            int j = 0;
110            int k = 0;
111            int l = 0;
112            int maxb = 0;
113            int nr = 0;
114            int ns = 0;
115            int nv = 0;
116            double absw = 0;
117            double ovfl = 0;
118            double smlnum = 0;
119            double tau = 0;
120            double temp = 0;
121            double tst1 = 0;
122            double ulp = 0;
123            double unfl = 0;
124            double[,] s = new double[0,0];
125            double[] v = new double[0];
126            double[] vv = new double[0];
127            double[] workc1 = new double[0];
128            double[] works1 = new double[0];
129            double[] workv3 = new double[0];
130            double[] tmpwr = new double[0];
131            double[] tmpwi = new double[0];
132            bool initz = new bool();
133            bool wantt = new bool();
134            bool wantz = new bool();
135            double cnst = 0;
136            bool failflag = new bool();
137            int p1 = 0;
138            int p2 = 0;
139            double vt = 0;
140            int i_ = 0;
141            int i1_ = 0;
142
143           
144            //
145            // Set the order of the multi-shift QR algorithm to be used.
146            // If you want to tune algorithm, change this values
147            //
148            ns = 12;
149            maxb = 50;
150           
151            //
152            // Now 2 < NS <= MAXB < NH.
153            //
154            maxb = Math.Max(3, maxb);
155            ns = Math.Min(maxb, ns);
156           
157            //
158            // Initialize
159            //
160            cnst = 1.5;
161            work = new double[Math.Max(n, 1)+1];
162            s = new double[ns+1, ns+1];
163            v = new double[ns+1+1];
164            vv = new double[ns+1+1];
165            wr = new double[Math.Max(n, 1)+1];
166            wi = new double[Math.Max(n, 1)+1];
167            workc1 = new double[1+1];
168            works1 = new double[1+1];
169            workv3 = new double[3+1];
170            tmpwr = new double[Math.Max(n, 1)+1];
171            tmpwi = new double[Math.Max(n, 1)+1];
172            System.Diagnostics.Debug.Assert(n>=0, "InternalSchurDecomposition: incorrect N!");
173            System.Diagnostics.Debug.Assert(tneeded==0 | tneeded==1, "InternalSchurDecomposition: incorrect TNeeded!");
174            System.Diagnostics.Debug.Assert(zneeded==0 | zneeded==1 | zneeded==2, "InternalSchurDecomposition: incorrect ZNeeded!");
175            wantt = tneeded==1;
176            initz = zneeded==2;
177            wantz = zneeded!=0;
178            info = 0;
179           
180            //
181            // Initialize Z, if necessary
182            //
183            if( initz )
184            {
185                z = new double[n+1, n+1];
186                for(i=1; i<=n; i++)
187                {
188                    for(j=1; j<=n; j++)
189                    {
190                        if( i==j )
191                        {
192                            z[i,j] = 1;
193                        }
194                        else
195                        {
196                            z[i,j] = 0;
197                        }
198                    }
199                }
200            }
201           
202            //
203            // Quick return if possible
204            //
205            if( n==0 )
206            {
207                return;
208            }
209            if( n==1 )
210            {
211                wr[1] = h[1,1];
212                wi[1] = 0;
213                return;
214            }
215           
216            //
217            // Set rows and columns 1 to N to zero below the first
218            // subdiagonal.
219            //
220            for(j=1; j<=n-2; j++)
221            {
222                for(i=j+2; i<=n; i++)
223                {
224                    h[i,j] = 0;
225                }
226            }
227           
228            //
229            // Test if N is sufficiently small
230            //
231            if( ns<=2 | ns>n | maxb>=n )
232            {
233               
234                //
235                // Use the standard double-shift algorithm
236                //
237                internalauxschur(wantt, wantz, n, 1, n, ref h, ref wr, ref wi, 1, n, ref z, ref work, ref workv3, ref workc1, ref works1, ref info);
238               
239                //
240                // fill entries under diagonal blocks of T with zeros
241                //
242                if( wantt )
243                {
244                    j = 1;
245                    while( j<=n )
246                    {
247                        if( (double)(wi[j])==(double)(0) )
248                        {
249                            for(i=j+1; i<=n; i++)
250                            {
251                                h[i,j] = 0;
252                            }
253                            j = j+1;
254                        }
255                        else
256                        {
257                            for(i=j+2; i<=n; i++)
258                            {
259                                h[i,j] = 0;
260                                h[i,j+1] = 0;
261                            }
262                            j = j+2;
263                        }
264                    }
265                }
266                return;
267            }
268            unfl = AP.Math.MinRealNumber;
269            ovfl = 1/unfl;
270            ulp = 2*AP.Math.MachineEpsilon;
271            smlnum = unfl*(n/ulp);
272           
273            //
274            // I1 and I2 are the indices of the first row and last column of H
275            // to which transformations must be applied. If eigenvalues only are
276            // being computed, I1 and I2 are set inside the main loop.
277            //
278            i1 = 1;
279            i2 = n;
280           
281            //
282            // ITN is the total number of multiple-shift QR iterations allowed.
283            //
284            itn = 30*n;
285           
286            //
287            // The main loop begins here. I is the loop index and decreases from
288            // IHI to ILO in steps of at most MAXB. Each iteration of the loop
289            // works with the active submatrix in rows and columns L to I.
290            // Eigenvalues I+1 to IHI have already converged. Either L = ILO or
291            // H(L,L-1) is negligible so that the matrix splits.
292            //
293            i = n;
294            while( true )
295            {
296                l = 1;
297                if( i<1 )
298                {
299                   
300                    //
301                    // fill entries under diagonal blocks of T with zeros
302                    //
303                    if( wantt )
304                    {
305                        j = 1;
306                        while( j<=n )
307                        {
308                            if( (double)(wi[j])==(double)(0) )
309                            {
310                                for(i=j+1; i<=n; i++)
311                                {
312                                    h[i,j] = 0;
313                                }
314                                j = j+1;
315                            }
316                            else
317                            {
318                                for(i=j+2; i<=n; i++)
319                                {
320                                    h[i,j] = 0;
321                                    h[i,j+1] = 0;
322                                }
323                                j = j+2;
324                            }
325                        }
326                    }
327                   
328                    //
329                    // Exit
330                    //
331                    return;
332                }
333               
334                //
335                // Perform multiple-shift QR iterations on rows and columns ILO to I
336                // until a submatrix of order at most MAXB splits off at the bottom
337                // because a subdiagonal element has become negligible.
338                //
339                failflag = true;
340                for(its=0; its<=itn; its++)
341                {
342                   
343                    //
344                    // Look for a single small subdiagonal element.
345                    //
346                    for(k=i; k>=l+1; k--)
347                    {
348                        tst1 = Math.Abs(h[k-1,k-1])+Math.Abs(h[k,k]);
349                        if( (double)(tst1)==(double)(0) )
350                        {
351                            tst1 = blas.upperhessenberg1norm(ref h, l, i, l, i, ref work);
352                        }
353                        if( (double)(Math.Abs(h[k,k-1]))<=(double)(Math.Max(ulp*tst1, smlnum)) )
354                        {
355                            break;
356                        }
357                    }
358                    l = k;
359                    if( l>1 )
360                    {
361                       
362                        //
363                        // H(L,L-1) is negligible.
364                        //
365                        h[l,l-1] = 0;
366                    }
367                   
368                    //
369                    // Exit from loop if a submatrix of order <= MAXB has split off.
370                    //
371                    if( l>=i-maxb+1 )
372                    {
373                        failflag = false;
374                        break;
375                    }
376                   
377                    //
378                    // Now the active submatrix is in rows and columns L to I. If
379                    // eigenvalues only are being computed, only the active submatrix
380                    // need be transformed.
381                    //
382                    if( its==20 | its==30 )
383                    {
384                       
385                        //
386                        // Exceptional shifts.
387                        //
388                        for(ii=i-ns+1; ii<=i; ii++)
389                        {
390                            wr[ii] = cnst*(Math.Abs(h[ii,ii-1])+Math.Abs(h[ii,ii]));
391                            wi[ii] = 0;
392                        }
393                    }
394                    else
395                    {
396                       
397                        //
398                        // Use eigenvalues of trailing submatrix of order NS as shifts.
399                        //
400                        blas.copymatrix(ref h, i-ns+1, i, i-ns+1, i, ref s, 1, ns, 1, ns);
401                        internalauxschur(false, false, ns, 1, ns, ref s, ref tmpwr, ref tmpwi, 1, ns, ref z, ref work, ref workv3, ref workc1, ref works1, ref ierr);
402                        for(p1=1; p1<=ns; p1++)
403                        {
404                            wr[i-ns+p1] = tmpwr[p1];
405                            wi[i-ns+p1] = tmpwi[p1];
406                        }
407                        if( ierr>0 )
408                        {
409                           
410                            //
411                            // If DLAHQR failed to compute all NS eigenvalues, use the
412                            // unconverged diagonal elements as the remaining shifts.
413                            //
414                            for(ii=1; ii<=ierr; ii++)
415                            {
416                                wr[i-ns+ii] = s[ii,ii];
417                                wi[i-ns+ii] = 0;
418                            }
419                        }
420                    }
421                   
422                    //
423                    // Form the first column of (G-w(1)) (G-w(2)) . . . (G-w(ns))
424                    // where G is the Hessenberg submatrix H(L:I,L:I) and w is
425                    // the vector of shifts (stored in WR and WI). The result is
426                    // stored in the local array V.
427                    //
428                    v[1] = 1;
429                    for(ii=2; ii<=ns+1; ii++)
430                    {
431                        v[ii] = 0;
432                    }
433                    nv = 1;
434                    for(j=i-ns+1; j<=i; j++)
435                    {
436                        if( (double)(wi[j])>=(double)(0) )
437                        {
438                            if( (double)(wi[j])==(double)(0) )
439                            {
440                               
441                                //
442                                // real shift
443                                //
444                                p1 = nv+1;
445                                for(i_=1; i_<=p1;i_++)
446                                {
447                                    vv[i_] = v[i_];
448                                }
449                                blas.matrixvectormultiply(ref h, l, l+nv, l, l+nv-1, false, ref vv, 1, nv, 1.0, ref v, 1, nv+1, -wr[j]);
450                                nv = nv+1;
451                            }
452                            else
453                            {
454                                if( (double)(wi[j])>(double)(0) )
455                                {
456                                   
457                                    //
458                                    // complex conjugate pair of shifts
459                                    //
460                                    p1 = nv+1;
461                                    for(i_=1; i_<=p1;i_++)
462                                    {
463                                        vv[i_] = v[i_];
464                                    }
465                                    blas.matrixvectormultiply(ref h, l, l+nv, l, l+nv-1, false, ref v, 1, nv, 1.0, ref vv, 1, nv+1, -(2*wr[j]));
466                                    itemp = blas.vectoridxabsmax(ref vv, 1, nv+1);
467                                    temp = 1/Math.Max(Math.Abs(vv[itemp]), smlnum);
468                                    p1 = nv+1;
469                                    for(i_=1; i_<=p1;i_++)
470                                    {
471                                        vv[i_] = temp*vv[i_];
472                                    }
473                                    absw = blas.pythag2(wr[j], wi[j]);
474                                    temp = temp*absw*absw;
475                                    blas.matrixvectormultiply(ref h, l, l+nv+1, l, l+nv, false, ref vv, 1, nv+1, 1.0, ref v, 1, nv+2, temp);
476                                    nv = nv+2;
477                                }
478                            }
479                           
480                            //
481                            // Scale V(1:NV) so that max(abs(V(i))) = 1. If V is zero,
482                            // reset it to the unit vector.
483                            //
484                            itemp = blas.vectoridxabsmax(ref v, 1, nv);
485                            temp = Math.Abs(v[itemp]);
486                            if( (double)(temp)==(double)(0) )
487                            {
488                                v[1] = 1;
489                                for(ii=2; ii<=nv; ii++)
490                                {
491                                    v[ii] = 0;
492                                }
493                            }
494                            else
495                            {
496                                temp = Math.Max(temp, smlnum);
497                                vt = 1/temp;
498                                for(i_=1; i_<=nv;i_++)
499                                {
500                                    v[i_] = vt*v[i_];
501                                }
502                            }
503                        }
504                    }
505                   
506                    //
507                    // Multiple-shift QR step
508                    //
509                    for(k=l; k<=i-1; k++)
510                    {
511                       
512                        //
513                        // The first iteration of this loop determines a reflection G
514                        // from the vector V and applies it from left and right to H,
515                        // thus creating a nonzero bulge below the subdiagonal.
516                        //
517                        // Each subsequent iteration determines a reflection G to
518                        // restore the Hessenberg form in the (K-1)th column, and thus
519                        // chases the bulge one step toward the bottom of the active
520                        // submatrix. NR is the order of G.
521                        //
522                        nr = Math.Min(ns+1, i-k+1);
523                        if( k>l )
524                        {
525                            p1 = k-1;
526                            p2 = k+nr-1;
527                            i1_ = (k) - (1);
528                            for(i_=1; i_<=nr;i_++)
529                            {
530                                v[i_] = h[i_+i1_,p1];
531                            }
532                        }
533                        reflections.generatereflection(ref v, nr, ref tau);
534                        if( k>l )
535                        {
536                            h[k,k-1] = v[1];
537                            for(ii=k+1; ii<=i; ii++)
538                            {
539                                h[ii,k-1] = 0;
540                            }
541                        }
542                        v[1] = 1;
543                       
544                        //
545                        // Apply G from the left to transform the rows of the matrix in
546                        // columns K to I2.
547                        //
548                        reflections.applyreflectionfromtheleft(ref h, tau, ref v, k, k+nr-1, k, i2, ref work);
549                       
550                        //
551                        // Apply G from the right to transform the columns of the
552                        // matrix in rows I1 to min(K+NR,I).
553                        //
554                        reflections.applyreflectionfromtheright(ref h, tau, ref v, i1, Math.Min(k+nr, i), k, k+nr-1, ref work);
555                        if( wantz )
556                        {
557                           
558                            //
559                            // Accumulate transformations in the matrix Z
560                            //
561                            reflections.applyreflectionfromtheright(ref z, tau, ref v, 1, n, k, k+nr-1, ref work);
562                        }
563                    }
564                }
565               
566                //
567                // Failure to converge in remaining number of iterations
568                //
569                if( failflag )
570                {
571                    info = i;
572                    return;
573                }
574               
575                //
576                // A submatrix of order <= MAXB in rows and columns L to I has split
577                // off. Use the double-shift QR algorithm to handle it.
578                //
579                internalauxschur(wantt, wantz, n, l, i, ref h, ref wr, ref wi, 1, n, ref z, ref work, ref workv3, ref workc1, ref works1, ref info);
580                if( info>0 )
581                {
582                    return;
583                }
584               
585                //
586                // Decrement number of remaining iterations, and return to start of
587                // the main loop with a new value of I.
588                //
589                itn = itn-its;
590                i = l-1;
591            }
592        }
593
594
595        private static void internalauxschur(bool wantt,
596            bool wantz,
597            int n,
598            int ilo,
599            int ihi,
600            ref double[,] h,
601            ref double[] wr,
602            ref double[] wi,
603            int iloz,
604            int ihiz,
605            ref double[,] z,
606            ref double[] work,
607            ref double[] workv3,
608            ref double[] workc1,
609            ref double[] works1,
610            ref int info)
611        {
612            int i = 0;
613            int i1 = 0;
614            int i2 = 0;
615            int itn = 0;
616            int its = 0;
617            int j = 0;
618            int k = 0;
619            int l = 0;
620            int m = 0;
621            int nh = 0;
622            int nr = 0;
623            int nz = 0;
624            double ave = 0;
625            double cs = 0;
626            double disc = 0;
627            double h00 = 0;
628            double h10 = 0;
629            double h11 = 0;
630            double h12 = 0;
631            double h21 = 0;
632            double h22 = 0;
633            double h33 = 0;
634            double h33s = 0;
635            double h43h34 = 0;
636            double h44 = 0;
637            double h44s = 0;
638            double ovfl = 0;
639            double s = 0;
640            double smlnum = 0;
641            double sn = 0;
642            double sum = 0;
643            double t1 = 0;
644            double t2 = 0;
645            double t3 = 0;
646            double tst1 = 0;
647            double unfl = 0;
648            double v1 = 0;
649            double v2 = 0;
650            double v3 = 0;
651            bool failflag = new bool();
652            double dat1 = 0;
653            double dat2 = 0;
654            int p1 = 0;
655            double him1im1 = 0;
656            double him1i = 0;
657            double hiim1 = 0;
658            double hii = 0;
659            double wrim1 = 0;
660            double wri = 0;
661            double wiim1 = 0;
662            double wii = 0;
663            double ulp = 0;
664
665            info = 0;
666            dat1 = 0.75;
667            dat2 = -0.4375;
668            ulp = AP.Math.MachineEpsilon;
669           
670            //
671            // Quick return if possible
672            //
673            if( n==0 )
674            {
675                return;
676            }
677            if( ilo==ihi )
678            {
679                wr[ilo] = h[ilo,ilo];
680                wi[ilo] = 0;
681                return;
682            }
683            nh = ihi-ilo+1;
684            nz = ihiz-iloz+1;
685           
686            //
687            // Set machine-dependent constants for the stopping criterion.
688            // If norm(H) <= sqrt(OVFL), overflow should not occur.
689            //
690            unfl = AP.Math.MinRealNumber;
691            ovfl = 1/unfl;
692            smlnum = unfl*(nh/ulp);
693           
694            //
695            // I1 and I2 are the indices of the first row and last column of H
696            // to which transformations must be applied. If eigenvalues only are
697            // being computed, I1 and I2 are set inside the main loop.
698            //
699            i1 = 1;
700            i2 = n;
701           
702            //
703            // ITN is the total number of QR iterations allowed.
704            //
705            itn = 30*nh;
706           
707            //
708            // The main loop begins here. I is the loop index and decreases from
709            // IHI to ILO in steps of 1 or 2. Each iteration of the loop works
710            // with the active submatrix in rows and columns L to I.
711            // Eigenvalues I+1 to IHI have already converged. Either L = ILO or
712            // H(L,L-1) is negligible so that the matrix splits.
713            //
714            i = ihi;
715            while( true )
716            {
717                l = ilo;
718                if( i<ilo )
719                {
720                    return;
721                }
722               
723                //
724                // Perform QR iterations on rows and columns ILO to I until a
725                // submatrix of order 1 or 2 splits off at the bottom because a
726                // subdiagonal element has become negligible.
727                //
728                failflag = true;
729                for(its=0; its<=itn; its++)
730                {
731                   
732                    //
733                    // Look for a single small subdiagonal element.
734                    //
735                    for(k=i; k>=l+1; k--)
736                    {
737                        tst1 = Math.Abs(h[k-1,k-1])+Math.Abs(h[k,k]);
738                        if( (double)(tst1)==(double)(0) )
739                        {
740                            tst1 = blas.upperhessenberg1norm(ref h, l, i, l, i, ref work);
741                        }
742                        if( (double)(Math.Abs(h[k,k-1]))<=(double)(Math.Max(ulp*tst1, smlnum)) )
743                        {
744                            break;
745                        }
746                    }
747                    l = k;
748                    if( l>ilo )
749                    {
750                       
751                        //
752                        // H(L,L-1) is negligible
753                        //
754                        h[l,l-1] = 0;
755                    }
756                   
757                    //
758                    // Exit from loop if a submatrix of order 1 or 2 has split off.
759                    //
760                    if( l>=i-1 )
761                    {
762                        failflag = false;
763                        break;
764                    }
765                   
766                    //
767                    // Now the active submatrix is in rows and columns L to I. If
768                    // eigenvalues only are being computed, only the active submatrix
769                    // need be transformed.
770                    //
771                    if( its==10 | its==20 )
772                    {
773                       
774                        //
775                        // Exceptional shift.
776                        //
777                        s = Math.Abs(h[i,i-1])+Math.Abs(h[i-1,i-2]);
778                        h44 = dat1*s+h[i,i];
779                        h33 = h44;
780                        h43h34 = dat2*s*s;
781                    }
782                    else
783                    {
784                       
785                        //
786                        // Prepare to use Francis' double shift
787                        // (i.e. 2nd degree generalized Rayleigh quotient)
788                        //
789                        h44 = h[i,i];
790                        h33 = h[i-1,i-1];
791                        h43h34 = h[i,i-1]*h[i-1,i];
792                        s = h[i-1,i-2]*h[i-1,i-2];
793                        disc = (h33-h44)*0.5;
794                        disc = disc*disc+h43h34;
795                        if( (double)(disc)>(double)(0) )
796                        {
797                           
798                            //
799                            // Real roots: use Wilkinson's shift twice
800                            //
801                            disc = Math.Sqrt(disc);
802                            ave = 0.5*(h33+h44);
803                            if( (double)(Math.Abs(h33)-Math.Abs(h44))>(double)(0) )
804                            {
805                                h33 = h33*h44-h43h34;
806                                h44 = h33/(extschursign(disc, ave)+ave);
807                            }
808                            else
809                            {
810                                h44 = extschursign(disc, ave)+ave;
811                            }
812                            h33 = h44;
813                            h43h34 = 0;
814                        }
815                    }
816                   
817                    //
818                    // Look for two consecutive small subdiagonal elements.
819                    //
820                    for(m=i-2; m>=l; m--)
821                    {
822                       
823                        //
824                        // Determine the effect of starting the double-shift QR
825                        // iteration at row M, and see if this would make H(M,M-1)
826                        // negligible.
827                        //
828                        h11 = h[m,m];
829                        h22 = h[m+1,m+1];
830                        h21 = h[m+1,m];
831                        h12 = h[m,m+1];
832                        h44s = h44-h11;
833                        h33s = h33-h11;
834                        v1 = (h33s*h44s-h43h34)/h21+h12;
835                        v2 = h22-h11-h33s-h44s;
836                        v3 = h[m+2,m+1];
837                        s = Math.Abs(v1)+Math.Abs(v2)+Math.Abs(v3);
838                        v1 = v1/s;
839                        v2 = v2/s;
840                        v3 = v3/s;
841                        workv3[1] = v1;
842                        workv3[2] = v2;
843                        workv3[3] = v3;
844                        if( m==l )
845                        {
846                            break;
847                        }
848                        h00 = h[m-1,m-1];
849                        h10 = h[m,m-1];
850                        tst1 = Math.Abs(v1)*(Math.Abs(h00)+Math.Abs(h11)+Math.Abs(h22));
851                        if( (double)(Math.Abs(h10)*(Math.Abs(v2)+Math.Abs(v3)))<=(double)(ulp*tst1) )
852                        {
853                            break;
854                        }
855                    }
856                   
857                    //
858                    // Double-shift QR step
859                    //
860                    for(k=m; k<=i-1; k++)
861                    {
862                       
863                        //
864                        // The first iteration of this loop determines a reflection G
865                        // from the vector V and applies it from left and right to H,
866                        // thus creating a nonzero bulge below the subdiagonal.
867                        //
868                        // Each subsequent iteration determines a reflection G to
869                        // restore the Hessenberg form in the (K-1)th column, and thus
870                        // chases the bulge one step toward the bottom of the active
871                        // submatrix. NR is the order of G.
872                        //
873                        nr = Math.Min(3, i-k+1);
874                        if( k>m )
875                        {
876                            for(p1=1; p1<=nr; p1++)
877                            {
878                                workv3[p1] = h[k+p1-1,k-1];
879                            }
880                        }
881                        reflections.generatereflection(ref workv3, nr, ref t1);
882                        if( k>m )
883                        {
884                            h[k,k-1] = workv3[1];
885                            h[k+1,k-1] = 0;
886                            if( k<i-1 )
887                            {
888                                h[k+2,k-1] = 0;
889                            }
890                        }
891                        else
892                        {
893                            if( m>l )
894                            {
895                                h[k,k-1] = -h[k,k-1];
896                            }
897                        }
898                        v2 = workv3[2];
899                        t2 = t1*v2;
900                        if( nr==3 )
901                        {
902                            v3 = workv3[3];
903                            t3 = t1*v3;
904                           
905                            //
906                            // Apply G from the left to transform the rows of the matrix
907                            // in columns K to I2.
908                            //
909                            for(j=k; j<=i2; j++)
910                            {
911                                sum = h[k,j]+v2*h[k+1,j]+v3*h[k+2,j];
912                                h[k,j] = h[k,j]-sum*t1;
913                                h[k+1,j] = h[k+1,j]-sum*t2;
914                                h[k+2,j] = h[k+2,j]-sum*t3;
915                            }
916                           
917                            //
918                            // Apply G from the right to transform the columns of the
919                            // matrix in rows I1 to min(K+3,I).
920                            //
921                            for(j=i1; j<=Math.Min(k+3, i); j++)
922                            {
923                                sum = h[j,k]+v2*h[j,k+1]+v3*h[j,k+2];
924                                h[j,k] = h[j,k]-sum*t1;
925                                h[j,k+1] = h[j,k+1]-sum*t2;
926                                h[j,k+2] = h[j,k+2]-sum*t3;
927                            }
928                            if( wantz )
929                            {
930                               
931                                //
932                                // Accumulate transformations in the matrix Z
933                                //
934                                for(j=iloz; j<=ihiz; j++)
935                                {
936                                    sum = z[j,k]+v2*z[j,k+1]+v3*z[j,k+2];
937                                    z[j,k] = z[j,k]-sum*t1;
938                                    z[j,k+1] = z[j,k+1]-sum*t2;
939                                    z[j,k+2] = z[j,k+2]-sum*t3;
940                                }
941                            }
942                        }
943                        else
944                        {
945                            if( nr==2 )
946                            {
947                               
948                                //
949                                // Apply G from the left to transform the rows of the matrix
950                                // in columns K to I2.
951                                //
952                                for(j=k; j<=i2; j++)
953                                {
954                                    sum = h[k,j]+v2*h[k+1,j];
955                                    h[k,j] = h[k,j]-sum*t1;
956                                    h[k+1,j] = h[k+1,j]-sum*t2;
957                                }
958                               
959                                //
960                                // Apply G from the right to transform the columns of the
961                                // matrix in rows I1 to min(K+3,I).
962                                //
963                                for(j=i1; j<=i; j++)
964                                {
965                                    sum = h[j,k]+v2*h[j,k+1];
966                                    h[j,k] = h[j,k]-sum*t1;
967                                    h[j,k+1] = h[j,k+1]-sum*t2;
968                                }
969                                if( wantz )
970                                {
971                                   
972                                    //
973                                    // Accumulate transformations in the matrix Z
974                                    //
975                                    for(j=iloz; j<=ihiz; j++)
976                                    {
977                                        sum = z[j,k]+v2*z[j,k+1];
978                                        z[j,k] = z[j,k]-sum*t1;
979                                        z[j,k+1] = z[j,k+1]-sum*t2;
980                                    }
981                                }
982                            }
983                        }
984                    }
985                }
986                if( failflag )
987                {
988                   
989                    //
990                    // Failure to converge in remaining number of iterations
991                    //
992                    info = i;
993                    return;
994                }
995                if( l==i )
996                {
997                   
998                    //
999                    // H(I,I-1) is negligible: one eigenvalue has converged.
1000                    //
1001                    wr[i] = h[i,i];
1002                    wi[i] = 0;
1003                }
1004                else
1005                {
1006                    if( l==i-1 )
1007                    {
1008                       
1009                        //
1010                        // H(I-1,I-2) is negligible: a pair of eigenvalues have converged.
1011                        //
1012                        //        Transform the 2-by-2 submatrix to standard Schur form,
1013                        //        and compute and store the eigenvalues.
1014                        //
1015                        him1im1 = h[i-1,i-1];
1016                        him1i = h[i-1,i];
1017                        hiim1 = h[i,i-1];
1018                        hii = h[i,i];
1019                        aux2x2schur(ref him1im1, ref him1i, ref hiim1, ref hii, ref wrim1, ref wiim1, ref wri, ref wii, ref cs, ref sn);
1020                        wr[i-1] = wrim1;
1021                        wi[i-1] = wiim1;
1022                        wr[i] = wri;
1023                        wi[i] = wii;
1024                        h[i-1,i-1] = him1im1;
1025                        h[i-1,i] = him1i;
1026                        h[i,i-1] = hiim1;
1027                        h[i,i] = hii;
1028                        if( wantt )
1029                        {
1030                           
1031                            //
1032                            // Apply the transformation to the rest of H.
1033                            //
1034                            if( i2>i )
1035                            {
1036                                workc1[1] = cs;
1037                                works1[1] = sn;
1038                                rotations.applyrotationsfromtheleft(true, i-1, i, i+1, i2, ref workc1, ref works1, ref h, ref work);
1039                            }
1040                            workc1[1] = cs;
1041                            works1[1] = sn;
1042                            rotations.applyrotationsfromtheright(true, i1, i-2, i-1, i, ref workc1, ref works1, ref h, ref work);
1043                        }
1044                        if( wantz )
1045                        {
1046                           
1047                            //
1048                            // Apply the transformation to Z.
1049                            //
1050                            workc1[1] = cs;
1051                            works1[1] = sn;
1052                            rotations.applyrotationsfromtheright(true, iloz, iloz+nz-1, i-1, i, ref workc1, ref works1, ref z, ref work);
1053                        }
1054                    }
1055                }
1056               
1057                //
1058                // Decrement number of remaining iterations, and return to start of
1059                // the main loop with new value of I.
1060                //
1061                itn = itn-its;
1062                i = l-1;
1063            }
1064        }
1065
1066
1067        private static void aux2x2schur(ref double a,
1068            ref double b,
1069            ref double c,
1070            ref double d,
1071            ref double rt1r,
1072            ref double rt1i,
1073            ref double rt2r,
1074            ref double rt2i,
1075            ref double cs,
1076            ref double sn)
1077        {
1078            double multpl = 0;
1079            double aa = 0;
1080            double bb = 0;
1081            double bcmax = 0;
1082            double bcmis = 0;
1083            double cc = 0;
1084            double cs1 = 0;
1085            double dd = 0;
1086            double eps = 0;
1087            double p = 0;
1088            double sab = 0;
1089            double sac = 0;
1090            double scl = 0;
1091            double sigma = 0;
1092            double sn1 = 0;
1093            double tau = 0;
1094            double temp = 0;
1095            double z = 0;
1096
1097            multpl = 4.0;
1098            eps = AP.Math.MachineEpsilon;
1099            if( (double)(c)==(double)(0) )
1100            {
1101                cs = 1;
1102                sn = 0;
1103            }
1104            else
1105            {
1106                if( (double)(b)==(double)(0) )
1107                {
1108                   
1109                    //
1110                    // Swap rows and columns
1111                    //
1112                    cs = 0;
1113                    sn = 1;
1114                    temp = d;
1115                    d = a;
1116                    a = temp;
1117                    b = -c;
1118                    c = 0;
1119                }
1120                else
1121                {
1122                    if( (double)(a-d)==(double)(0) & extschursigntoone(b)!=extschursigntoone(c) )
1123                    {
1124                        cs = 1;
1125                        sn = 0;
1126                    }
1127                    else
1128                    {
1129                        temp = a-d;
1130                        p = 0.5*temp;
1131                        bcmax = Math.Max(Math.Abs(b), Math.Abs(c));
1132                        bcmis = Math.Min(Math.Abs(b), Math.Abs(c))*extschursigntoone(b)*extschursigntoone(c);
1133                        scl = Math.Max(Math.Abs(p), bcmax);
1134                        z = p/scl*p+bcmax/scl*bcmis;
1135                       
1136                        //
1137                        // If Z is of the order of the machine accuracy, postpone the
1138                        // decision on the nature of eigenvalues
1139                        //
1140                        if( (double)(z)>=(double)(multpl*eps) )
1141                        {
1142                           
1143                            //
1144                            // Real eigenvalues. Compute A and D.
1145                            //
1146                            z = p+extschursign(Math.Sqrt(scl)*Math.Sqrt(z), p);
1147                            a = d+z;
1148                            d = d-bcmax/z*bcmis;
1149                           
1150                            //
1151                            // Compute B and the rotation matrix
1152                            //
1153                            tau = blas.pythag2(c, z);
1154                            cs = z/tau;
1155                            sn = c/tau;
1156                            b = b-c;
1157                            c = 0;
1158                        }
1159                        else
1160                        {
1161                           
1162                            //
1163                            // Complex eigenvalues, or real (almost) equal eigenvalues.
1164                            // Make diagonal elements equal.
1165                            //
1166                            sigma = b+c;
1167                            tau = blas.pythag2(sigma, temp);
1168                            cs = Math.Sqrt(0.5*(1+Math.Abs(sigma)/tau));
1169                            sn = -(p/(tau*cs)*extschursign(1, sigma));
1170                           
1171                            //
1172                            // Compute [ AA  BB ] = [ A  B ] [ CS -SN ]
1173                            //         [ CC  DD ]   [ C  D ] [ SN  CS ]
1174                            //
1175                            aa = a*cs+b*sn;
1176                            bb = -(a*sn)+b*cs;
1177                            cc = c*cs+d*sn;
1178                            dd = -(c*sn)+d*cs;
1179                           
1180                            //
1181                            // Compute [ A  B ] = [ CS  SN ] [ AA  BB ]
1182                            //         [ C  D ]   [-SN  CS ] [ CC  DD ]
1183                            //
1184                            a = aa*cs+cc*sn;
1185                            b = bb*cs+dd*sn;
1186                            c = -(aa*sn)+cc*cs;
1187                            d = -(bb*sn)+dd*cs;
1188                            temp = 0.5*(a+d);
1189                            a = temp;
1190                            d = temp;
1191                            if( (double)(c)!=(double)(0) )
1192                            {
1193                                if( (double)(b)!=(double)(0) )
1194                                {
1195                                    if( extschursigntoone(b)==extschursigntoone(c) )
1196                                    {
1197                                       
1198                                        //
1199                                        // Real eigenvalues: reduce to upper triangular form
1200                                        //
1201                                        sab = Math.Sqrt(Math.Abs(b));
1202                                        sac = Math.Sqrt(Math.Abs(c));
1203                                        p = extschursign(sab*sac, c);
1204                                        tau = 1/Math.Sqrt(Math.Abs(b+c));
1205                                        a = temp+p;
1206                                        d = temp-p;
1207                                        b = b-c;
1208                                        c = 0;
1209                                        cs1 = sab*tau;
1210                                        sn1 = sac*tau;
1211                                        temp = cs*cs1-sn*sn1;
1212                                        sn = cs*sn1+sn*cs1;
1213                                        cs = temp;
1214                                    }
1215                                }
1216                                else
1217                                {
1218                                    b = -c;
1219                                    c = 0;
1220                                    temp = cs;
1221                                    cs = -sn;
1222                                    sn = temp;
1223                                }
1224                            }
1225                        }
1226                    }
1227                }
1228            }
1229           
1230            //
1231            // Store eigenvalues in (RT1R,RT1I) and (RT2R,RT2I).
1232            //
1233            rt1r = a;
1234            rt2r = d;
1235            if( (double)(c)==(double)(0) )
1236            {
1237                rt1i = 0;
1238                rt2i = 0;
1239            }
1240            else
1241            {
1242                rt1i = Math.Sqrt(Math.Abs(b))*Math.Sqrt(Math.Abs(c));
1243                rt2i = -rt1i;
1244            }
1245        }
1246
1247
1248        private static double extschursign(double a,
1249            double b)
1250        {
1251            double result = 0;
1252
1253            if( (double)(b)>=(double)(0) )
1254            {
1255                result = Math.Abs(a);
1256            }
1257            else
1258            {
1259                result = -Math.Abs(a);
1260            }
1261            return result;
1262        }
1263
1264
1265        private static int extschursigntoone(double b)
1266        {
1267            int result = 0;
1268
1269            if( (double)(b)>=(double)(0) )
1270            {
1271                result = 1;
1272            }
1273            else
1274            {
1275                result = -1;
1276            }
1277            return result;
1278        }
1279    }
1280}
Note: See TracBrowser for help on using the repository browser.