Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/ALGLIB/nsevd.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: 88.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 nsevd
32    {
33        /*************************************************************************
34        Finding eigenvalues and eigenvectors of a general matrix
35
36        The algorithm finds eigenvalues and eigenvectors of a general matrix by
37        using the QR algorithm with multiple shifts. The algorithm can find
38        eigenvalues and both left and right eigenvectors.
39
40        The right eigenvector is a vector x such that A*x = w*x, and the left
41        eigenvector is a vector y such that y'*A = w*y' (here y' implies a complex
42        conjugate transposition of vector y).
43
44        Input parameters:
45            A       -   matrix. Array whose indexes range within [0..N-1, 0..N-1].
46            N       -   size of matrix A.
47            VNeeded -   flag controlling whether eigenvectors are needed or not.
48                        If VNeeded is equal to:
49                         * 0, eigenvectors are not returned;
50                         * 1, right eigenvectors are returned;
51                         * 2, left eigenvectors are returned;
52                         * 3, both left and right eigenvectors are returned.
53
54        Output parameters:
55            WR      -   real parts of eigenvalues.
56                        Array whose index ranges within [0..N-1].
57            WR      -   imaginary parts of eigenvalues.
58                        Array whose index ranges within [0..N-1].
59            VL, VR  -   arrays of left and right eigenvectors (if they are needed).
60                        If WI[i]=0, the respective eigenvalue is a real number,
61                        and it corresponds to the column number I of matrices VL/VR.
62                        If WI[i]>0, we have a pair of complex conjugate numbers with
63                        positive and negative imaginary parts:
64                            the first eigenvalue WR[i] + sqrt(-1)*WI[i];
65                            the second eigenvalue WR[i+1] + sqrt(-1)*WI[i+1];
66                            WI[i]>0
67                            WI[i+1] = -WI[i] < 0
68                        In that case, the eigenvector  corresponding to the first
69                        eigenvalue is located in i and i+1 columns of matrices
70                        VL/VR (the column number i contains the real part, and the
71                        column number i+1 contains the imaginary part), and the vector
72                        corresponding to the second eigenvalue is a complex conjugate to
73                        the first vector.
74                        Arrays whose indexes range within [0..N-1, 0..N-1].
75
76        Result:
77            True, if the algorithm has converged.
78            False, if the algorithm has not converged.
79
80        Note 1:
81            Some users may ask the following question: what if WI[N-1]>0?
82            WI[N] must contain an eigenvalue which is complex conjugate to the
83            N-th eigenvalue, but the array has only size N?
84            The answer is as follows: such a situation cannot occur because the
85            algorithm finds a pairs of eigenvalues, therefore, if WI[i]>0, I is
86            strictly less than N-1.
87
88        Note 2:
89            The algorithm performance depends on the value of the internal parameter
90            NS of the InternalSchurDecomposition subroutine which defines the number
91            of shifts in the QR algorithm (similarly to the block width in block-matrix
92            algorithms of linear algebra). If you require maximum performance
93            on your machine, it is recommended to adjust this parameter manually.
94
95
96        See also the InternalTREVC subroutine.
97
98        The algorithm is based on the LAPACK 3.0 library.
99        *************************************************************************/
100        public static bool rmatrixevd(double[,] a,
101            int n,
102            int vneeded,
103            ref double[] wr,
104            ref double[] wi,
105            ref double[,] vl,
106            ref double[,] vr)
107        {
108            bool result = new bool();
109            double[,] a1 = new double[0,0];
110            double[,] vl1 = new double[0,0];
111            double[,] vr1 = new double[0,0];
112            double[] wr1 = new double[0];
113            double[] wi1 = new double[0];
114            int i = 0;
115            double mx = 0;
116            int i_ = 0;
117            int i1_ = 0;
118
119            a = (double[,])a.Clone();
120
121            System.Diagnostics.Debug.Assert(vneeded>=0 & vneeded<=3, "RMatrixEVD: incorrect VNeeded!");
122            a1 = new double[n+1, n+1];
123            for(i=1; i<=n; i++)
124            {
125                i1_ = (0) - (1);
126                for(i_=1; i_<=n;i_++)
127                {
128                    a1[i,i_] = a[i-1,i_+i1_];
129                }
130            }
131            result = nonsymmetricevd(a1, n, vneeded, ref wr1, ref wi1, ref vl1, ref vr1);
132            if( result )
133            {
134                wr = new double[n-1+1];
135                wi = new double[n-1+1];
136                i1_ = (1) - (0);
137                for(i_=0; i_<=n-1;i_++)
138                {
139                    wr[i_] = wr1[i_+i1_];
140                }
141                i1_ = (1) - (0);
142                for(i_=0; i_<=n-1;i_++)
143                {
144                    wi[i_] = wi1[i_+i1_];
145                }
146                if( vneeded==2 | vneeded==3 )
147                {
148                    vl = new double[n-1+1, n-1+1];
149                    for(i=0; i<=n-1; i++)
150                    {
151                        i1_ = (1) - (0);
152                        for(i_=0; i_<=n-1;i_++)
153                        {
154                            vl[i,i_] = vl1[i+1,i_+i1_];
155                        }
156                    }
157                }
158                if( vneeded==1 | vneeded==3 )
159                {
160                    vr = new double[n-1+1, n-1+1];
161                    for(i=0; i<=n-1; i++)
162                    {
163                        i1_ = (1) - (0);
164                        for(i_=0; i_<=n-1;i_++)
165                        {
166                            vr[i,i_] = vr1[i+1,i_+i1_];
167                        }
168                    }
169                }
170            }
171            return result;
172        }
173
174
175        public static bool nonsymmetricevd(double[,] a,
176            int n,
177            int vneeded,
178            ref double[] wr,
179            ref double[] wi,
180            ref double[,] vl,
181            ref double[,] vr)
182        {
183            bool result = new bool();
184            double[,] s = new double[0,0];
185            double[] tau = new double[0];
186            bool[] sel = new bool[0];
187            int i = 0;
188            int info = 0;
189            int m = 0;
190            int i_ = 0;
191
192            a = (double[,])a.Clone();
193
194            System.Diagnostics.Debug.Assert(vneeded>=0 & vneeded<=3, "NonSymmetricEVD: incorrect VNeeded!");
195            if( vneeded==0 )
196            {
197               
198                //
199                // Eigen values only
200                //
201                hessenberg.toupperhessenberg(ref a, n, ref tau);
202                hsschur.internalschurdecomposition(ref a, n, 0, 0, ref wr, ref wi, ref s, ref info);
203                result = info==0;
204                return result;
205            }
206           
207            //
208            // Eigen values and vectors
209            //
210            hessenberg.toupperhessenberg(ref a, n, ref tau);
211            hessenberg.unpackqfromupperhessenberg(ref a, n, ref tau, ref s);
212            hsschur.internalschurdecomposition(ref a, n, 1, 1, ref wr, ref wi, ref s, ref info);
213            result = info==0;
214            if( !result )
215            {
216                return result;
217            }
218            if( vneeded==1 | vneeded==3 )
219            {
220                vr = new double[n+1, n+1];
221                for(i=1; i<=n; i++)
222                {
223                    for(i_=1; i_<=n;i_++)
224                    {
225                        vr[i,i_] = s[i,i_];
226                    }
227                }
228            }
229            if( vneeded==2 | vneeded==3 )
230            {
231                vl = new double[n+1, n+1];
232                for(i=1; i<=n; i++)
233                {
234                    for(i_=1; i_<=n;i_++)
235                    {
236                        vl[i,i_] = s[i,i_];
237                    }
238                }
239            }
240            internaltrevc(ref a, n, vneeded, 1, sel, ref vl, ref vr, ref m, ref info);
241            result = info==0;
242            return result;
243        }
244
245
246        private static void internaltrevc(ref double[,] t,
247            int n,
248            int side,
249            int howmny,
250            bool[] vselect,
251            ref double[,] vl,
252            ref double[,] vr,
253            ref int m,
254            ref int info)
255        {
256            bool allv = new bool();
257            bool bothv = new bool();
258            bool leftv = new bool();
259            bool over = new bool();
260            bool pair = new bool();
261            bool rightv = new bool();
262            bool somev = new bool();
263            int i = 0;
264            int ierr = 0;
265            int ii = 0;
266            int ip = 0;
267            int iis = 0;
268            int j = 0;
269            int j1 = 0;
270            int j2 = 0;
271            int jnxt = 0;
272            int k = 0;
273            int ki = 0;
274            int n2 = 0;
275            double beta = 0;
276            double bignum = 0;
277            double emax = 0;
278            double ovfl = 0;
279            double rec = 0;
280            double remax = 0;
281            double scl = 0;
282            double smin = 0;
283            double smlnum = 0;
284            double ulp = 0;
285            double unfl = 0;
286            double vcrit = 0;
287            double vmax = 0;
288            double wi = 0;
289            double wr = 0;
290            double xnorm = 0;
291            double[,] x = new double[0,0];
292            double[] work = new double[0];
293            double[] temp = new double[0];
294            double[,] temp11 = new double[0,0];
295            double[,] temp22 = new double[0,0];
296            double[,] temp11b = new double[0,0];
297            double[,] temp21b = new double[0,0];
298            double[,] temp12b = new double[0,0];
299            double[,] temp22b = new double[0,0];
300            bool skipflag = new bool();
301            int k1 = 0;
302            int k2 = 0;
303            int k3 = 0;
304            int k4 = 0;
305            double vt = 0;
306            bool[] rswap4 = new bool[0];
307            bool[] zswap4 = new bool[0];
308            int[,] ipivot44 = new int[0,0];
309            double[] civ4 = new double[0];
310            double[] crv4 = new double[0];
311            int i_ = 0;
312            int i1_ = 0;
313
314            vselect = (bool[])vselect.Clone();
315
316            x = new double[2+1, 2+1];
317            temp11 = new double[1+1, 1+1];
318            temp11b = new double[1+1, 1+1];
319            temp21b = new double[2+1, 1+1];
320            temp12b = new double[1+1, 2+1];
321            temp22b = new double[2+1, 2+1];
322            temp22 = new double[2+1, 2+1];
323            work = new double[3*n+1];
324            temp = new double[n+1];
325            rswap4 = new bool[4+1];
326            zswap4 = new bool[4+1];
327            ipivot44 = new int[4+1, 4+1];
328            civ4 = new double[4+1];
329            crv4 = new double[4+1];
330            if( howmny!=1 )
331            {
332                if( side==1 | side==3 )
333                {
334                    vr = new double[n+1, n+1];
335                }
336                if( side==2 | side==3 )
337                {
338                    vl = new double[n+1, n+1];
339                }
340            }
341           
342            //
343            // Decode and test the input parameters
344            //
345            bothv = side==3;
346            rightv = side==1 | bothv;
347            leftv = side==2 | bothv;
348            allv = howmny==2;
349            over = howmny==1;
350            somev = howmny==3;
351            info = 0;
352            if( n<0 )
353            {
354                info = -2;
355                return;
356            }
357            if( !rightv & !leftv )
358            {
359                info = -3;
360                return;
361            }
362            if( !allv & !over & !somev )
363            {
364                info = -4;
365                return;
366            }
367           
368            //
369            // Set M to the number of columns required to store the selected
370            // eigenvectors, standardize the array SELECT if necessary, and
371            // test MM.
372            //
373            if( somev )
374            {
375                m = 0;
376                pair = false;
377                for(j=1; j<=n; j++)
378                {
379                    if( pair )
380                    {
381                        pair = false;
382                        vselect[j] = false;
383                    }
384                    else
385                    {
386                        if( j<n )
387                        {
388                            if( (double)(t[j+1,j])==(double)(0) )
389                            {
390                                if( vselect[j] )
391                                {
392                                    m = m+1;
393                                }
394                            }
395                            else
396                            {
397                                pair = true;
398                                if( vselect[j] | vselect[j+1] )
399                                {
400                                    vselect[j] = true;
401                                    m = m+2;
402                                }
403                            }
404                        }
405                        else
406                        {
407                            if( vselect[n] )
408                            {
409                                m = m+1;
410                            }
411                        }
412                    }
413                }
414            }
415            else
416            {
417                m = n;
418            }
419           
420            //
421            // Quick return if possible.
422            //
423            if( n==0 )
424            {
425                return;
426            }
427           
428            //
429            // Set the constants to control overflow.
430            //
431            unfl = AP.Math.MinRealNumber;
432            ovfl = 1/unfl;
433            ulp = AP.Math.MachineEpsilon;
434            smlnum = unfl*(n/ulp);
435            bignum = (1-ulp)/smlnum;
436           
437            //
438            // Compute 1-norm of each column of strictly upper triangular
439            // part of T to control overflow in triangular solver.
440            //
441            work[1] = 0;
442            for(j=2; j<=n; j++)
443            {
444                work[j] = 0;
445                for(i=1; i<=j-1; i++)
446                {
447                    work[j] = work[j]+Math.Abs(t[i,j]);
448                }
449            }
450           
451            //
452            // Index IP is used to specify the real or complex eigenvalue:
453            // IP = 0, real eigenvalue,
454            //      1, first of conjugate complex pair: (wr,wi)
455            //     -1, second of conjugate complex pair: (wr,wi)
456            //
457            n2 = 2*n;
458            if( rightv )
459            {
460               
461                //
462                // Compute right eigenvectors.
463                //
464                ip = 0;
465                iis = m;
466                for(ki=n; ki>=1; ki--)
467                {
468                    skipflag = false;
469                    if( ip==1 )
470                    {
471                        skipflag = true;
472                    }
473                    else
474                    {
475                        if( ki!=1 )
476                        {
477                            if( (double)(t[ki,ki-1])!=(double)(0) )
478                            {
479                                ip = -1;
480                            }
481                        }
482                        if( somev )
483                        {
484                            if( ip==0 )
485                            {
486                                if( !vselect[ki] )
487                                {
488                                    skipflag = true;
489                                }
490                            }
491                            else
492                            {
493                                if( !vselect[ki-1] )
494                                {
495                                    skipflag = true;
496                                }
497                            }
498                        }
499                    }
500                    if( !skipflag )
501                    {
502                       
503                        //
504                        // Compute the KI-th eigenvalue (WR,WI).
505                        //
506                        wr = t[ki,ki];
507                        wi = 0;
508                        if( ip!=0 )
509                        {
510                            wi = Math.Sqrt(Math.Abs(t[ki,ki-1]))*Math.Sqrt(Math.Abs(t[ki-1,ki]));
511                        }
512                        smin = Math.Max(ulp*(Math.Abs(wr)+Math.Abs(wi)), smlnum);
513                        if( ip==0 )
514                        {
515                           
516                            //
517                            // Real right eigenvector
518                            //
519                            work[ki+n] = 1;
520                           
521                            //
522                            // Form right-hand side
523                            //
524                            for(k=1; k<=ki-1; k++)
525                            {
526                                work[k+n] = -t[k,ki];
527                            }
528                           
529                            //
530                            // Solve the upper quasi-triangular system:
531                            //   (T(1:KI-1,1:KI-1) - WR)*X = SCALE*WORK.
532                            //
533                            jnxt = ki-1;
534                            for(j=ki-1; j>=1; j--)
535                            {
536                                if( j>jnxt )
537                                {
538                                    continue;
539                                }
540                                j1 = j;
541                                j2 = j;
542                                jnxt = j-1;
543                                if( j>1 )
544                                {
545                                    if( (double)(t[j,j-1])!=(double)(0) )
546                                    {
547                                        j1 = j-1;
548                                        jnxt = j-2;
549                                    }
550                                }
551                                if( j1==j2 )
552                                {
553                                   
554                                    //
555                                    // 1-by-1 diagonal block
556                                    //
557                                    temp11[1,1] = t[j,j];
558                                    temp11b[1,1] = work[j+n];
559                                    internalhsevdlaln2(false, 1, 1, smin, 1, ref temp11, 1.0, 1.0, ref temp11b, wr, 0.0, ref rswap4, ref zswap4, ref ipivot44, ref civ4, ref crv4, ref x, ref scl, ref xnorm, ref ierr);
560                                   
561                                    //
562                                    // Scale X(1,1) to avoid overflow when updating
563                                    // the right-hand side.
564                                    //
565                                    if( (double)(xnorm)>(double)(1) )
566                                    {
567                                        if( (double)(work[j])>(double)(bignum/xnorm) )
568                                        {
569                                            x[1,1] = x[1,1]/xnorm;
570                                            scl = scl/xnorm;
571                                        }
572                                    }
573                                   
574                                    //
575                                    // Scale if necessary
576                                    //
577                                    if( (double)(scl)!=(double)(1) )
578                                    {
579                                        k1 = n+1;
580                                        k2 = n+ki;
581                                        for(i_=k1; i_<=k2;i_++)
582                                        {
583                                            work[i_] = scl*work[i_];
584                                        }
585                                    }
586                                    work[j+n] = x[1,1];
587                                   
588                                    //
589                                    // Update right-hand side
590                                    //
591                                    k1 = 1+n;
592                                    k2 = j-1+n;
593                                    k3 = j-1;
594                                    vt = -x[1,1];
595                                    i1_ = (1) - (k1);
596                                    for(i_=k1; i_<=k2;i_++)
597                                    {
598                                        work[i_] = work[i_] + vt*t[i_+i1_,j];
599                                    }
600                                }
601                                else
602                                {
603                                   
604                                    //
605                                    // 2-by-2 diagonal block
606                                    //
607                                    temp22[1,1] = t[j-1,j-1];
608                                    temp22[1,2] = t[j-1,j];
609                                    temp22[2,1] = t[j,j-1];
610                                    temp22[2,2] = t[j,j];
611                                    temp21b[1,1] = work[j-1+n];
612                                    temp21b[2,1] = work[j+n];
613                                    internalhsevdlaln2(false, 2, 1, smin, 1.0, ref temp22, 1.0, 1.0, ref temp21b, wr, 0, ref rswap4, ref zswap4, ref ipivot44, ref civ4, ref crv4, ref x, ref scl, ref xnorm, ref ierr);
614                                   
615                                    //
616                                    // Scale X(1,1) and X(2,1) to avoid overflow when
617                                    // updating the right-hand side.
618                                    //
619                                    if( (double)(xnorm)>(double)(1) )
620                                    {
621                                        beta = Math.Max(work[j-1], work[j]);
622                                        if( (double)(beta)>(double)(bignum/xnorm) )
623                                        {
624                                            x[1,1] = x[1,1]/xnorm;
625                                            x[2,1] = x[2,1]/xnorm;
626                                            scl = scl/xnorm;
627                                        }
628                                    }
629                                   
630                                    //
631                                    // Scale if necessary
632                                    //
633                                    if( (double)(scl)!=(double)(1) )
634                                    {
635                                        k1 = 1+n;
636                                        k2 = ki+n;
637                                        for(i_=k1; i_<=k2;i_++)
638                                        {
639                                            work[i_] = scl*work[i_];
640                                        }
641                                    }
642                                    work[j-1+n] = x[1,1];
643                                    work[j+n] = x[2,1];
644                                   
645                                    //
646                                    // Update right-hand side
647                                    //
648                                    k1 = 1+n;
649                                    k2 = j-2+n;
650                                    k3 = j-2;
651                                    k4 = j-1;
652                                    vt = -x[1,1];
653                                    i1_ = (1) - (k1);
654                                    for(i_=k1; i_<=k2;i_++)
655                                    {
656                                        work[i_] = work[i_] + vt*t[i_+i1_,k4];
657                                    }
658                                    vt = -x[2,1];
659                                    i1_ = (1) - (k1);
660                                    for(i_=k1; i_<=k2;i_++)
661                                    {
662                                        work[i_] = work[i_] + vt*t[i_+i1_,j];
663                                    }
664                                }
665                            }
666                           
667                            //
668                            // Copy the vector x or Q*x to VR and normalize.
669                            //
670                            if( !over )
671                            {
672                                k1 = 1+n;
673                                k2 = ki+n;
674                                i1_ = (k1) - (1);
675                                for(i_=1; i_<=ki;i_++)
676                                {
677                                    vr[i_,iis] = work[i_+i1_];
678                                }
679                                ii = blas.columnidxabsmax(ref vr, 1, ki, iis);
680                                remax = 1/Math.Abs(vr[ii,iis]);
681                                for(i_=1; i_<=ki;i_++)
682                                {
683                                    vr[i_,iis] = remax*vr[i_,iis];
684                                }
685                                for(k=ki+1; k<=n; k++)
686                                {
687                                    vr[k,iis] = 0;
688                                }
689                            }
690                            else
691                            {
692                                if( ki>1 )
693                                {
694                                    for(i_=1; i_<=n;i_++)
695                                    {
696                                        temp[i_] = vr[i_,ki];
697                                    }
698                                    blas.matrixvectormultiply(ref vr, 1, n, 1, ki-1, false, ref work, 1+n, ki-1+n, 1.0, ref temp, 1, n, work[ki+n]);
699                                    for(i_=1; i_<=n;i_++)
700                                    {
701                                        vr[i_,ki] = temp[i_];
702                                    }
703                                }
704                                ii = blas.columnidxabsmax(ref vr, 1, n, ki);
705                                remax = 1/Math.Abs(vr[ii,ki]);
706                                for(i_=1; i_<=n;i_++)
707                                {
708                                    vr[i_,ki] = remax*vr[i_,ki];
709                                }
710                            }
711                        }
712                        else
713                        {
714                           
715                            //
716                            // Complex right eigenvector.
717                            //
718                            // Initial solve
719                            //     [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0.
720                            //     [ (T(KI,KI-1)   T(KI,KI)   )               ]
721                            //
722                            if( (double)(Math.Abs(t[ki-1,ki]))>=(double)(Math.Abs(t[ki,ki-1])) )
723                            {
724                                work[ki-1+n] = 1;
725                                work[ki+n2] = wi/t[ki-1,ki];
726                            }
727                            else
728                            {
729                                work[ki-1+n] = -(wi/t[ki,ki-1]);
730                                work[ki+n2] = 1;
731                            }
732                            work[ki+n] = 0;
733                            work[ki-1+n2] = 0;
734                           
735                            //
736                            // Form right-hand side
737                            //
738                            for(k=1; k<=ki-2; k++)
739                            {
740                                work[k+n] = -(work[ki-1+n]*t[k,ki-1]);
741                                work[k+n2] = -(work[ki+n2]*t[k,ki]);
742                            }
743                           
744                            //
745                            // Solve upper quasi-triangular system:
746                            // (T(1:KI-2,1:KI-2) - (WR+i*WI))*X = SCALE*(WORK+i*WORK2)
747                            //
748                            jnxt = ki-2;
749                            for(j=ki-2; j>=1; j--)
750                            {
751                                if( j>jnxt )
752                                {
753                                    continue;
754                                }
755                                j1 = j;
756                                j2 = j;
757                                jnxt = j-1;
758                                if( j>1 )
759                                {
760                                    if( (double)(t[j,j-1])!=(double)(0) )
761                                    {
762                                        j1 = j-1;
763                                        jnxt = j-2;
764                                    }
765                                }
766                                if( j1==j2 )
767                                {
768                                   
769                                    //
770                                    // 1-by-1 diagonal block
771                                    //
772                                    temp11[1,1] = t[j,j];
773                                    temp12b[1,1] = work[j+n];
774                                    temp12b[1,2] = work[j+n+n];
775                                    internalhsevdlaln2(false, 1, 2, smin, 1.0, ref temp11, 1.0, 1.0, ref temp12b, wr, wi, ref rswap4, ref zswap4, ref ipivot44, ref civ4, ref crv4, ref x, ref scl, ref xnorm, ref ierr);
776                                   
777                                    //
778                                    // Scale X(1,1) and X(1,2) to avoid overflow when
779                                    // updating the right-hand side.
780                                    //
781                                    if( (double)(xnorm)>(double)(1) )
782                                    {
783                                        if( (double)(work[j])>(double)(bignum/xnorm) )
784                                        {
785                                            x[1,1] = x[1,1]/xnorm;
786                                            x[1,2] = x[1,2]/xnorm;
787                                            scl = scl/xnorm;
788                                        }
789                                    }
790                                   
791                                    //
792                                    // Scale if necessary
793                                    //
794                                    if( (double)(scl)!=(double)(1) )
795                                    {
796                                        k1 = 1+n;
797                                        k2 = ki+n;
798                                        for(i_=k1; i_<=k2;i_++)
799                                        {
800                                            work[i_] = scl*work[i_];
801                                        }
802                                        k1 = 1+n2;
803                                        k2 = ki+n2;
804                                        for(i_=k1; i_<=k2;i_++)
805                                        {
806                                            work[i_] = scl*work[i_];
807                                        }
808                                    }
809                                    work[j+n] = x[1,1];
810                                    work[j+n2] = x[1,2];
811                                   
812                                    //
813                                    // Update the right-hand side
814                                    //
815                                    k1 = 1+n;
816                                    k2 = j-1+n;
817                                    k3 = 1;
818                                    k4 = j-1;
819                                    vt = -x[1,1];
820                                    i1_ = (k3) - (k1);
821                                    for(i_=k1; i_<=k2;i_++)
822                                    {
823                                        work[i_] = work[i_] + vt*t[i_+i1_,j];
824                                    }
825                                    k1 = 1+n2;
826                                    k2 = j-1+n2;
827                                    k3 = 1;
828                                    k4 = j-1;
829                                    vt = -x[1,2];
830                                    i1_ = (k3) - (k1);
831                                    for(i_=k1; i_<=k2;i_++)
832                                    {
833                                        work[i_] = work[i_] + vt*t[i_+i1_,j];
834                                    }
835                                }
836                                else
837                                {
838                                   
839                                    //
840                                    // 2-by-2 diagonal block
841                                    //
842                                    temp22[1,1] = t[j-1,j-1];
843                                    temp22[1,2] = t[j-1,j];
844                                    temp22[2,1] = t[j,j-1];
845                                    temp22[2,2] = t[j,j];
846                                    temp22b[1,1] = work[j-1+n];
847                                    temp22b[1,2] = work[j-1+n+n];
848                                    temp22b[2,1] = work[j+n];
849                                    temp22b[2,2] = work[j+n+n];
850                                    internalhsevdlaln2(false, 2, 2, smin, 1.0, ref temp22, 1.0, 1.0, ref temp22b, wr, wi, ref rswap4, ref zswap4, ref ipivot44, ref civ4, ref crv4, ref x, ref scl, ref xnorm, ref ierr);
851                                   
852                                    //
853                                    // Scale X to avoid overflow when updating
854                                    // the right-hand side.
855                                    //
856                                    if( (double)(xnorm)>(double)(1) )
857                                    {
858                                        beta = Math.Max(work[j-1], work[j]);
859                                        if( (double)(beta)>(double)(bignum/xnorm) )
860                                        {
861                                            rec = 1/xnorm;
862                                            x[1,1] = x[1,1]*rec;
863                                            x[1,2] = x[1,2]*rec;
864                                            x[2,1] = x[2,1]*rec;
865                                            x[2,2] = x[2,2]*rec;
866                                            scl = scl*rec;
867                                        }
868                                    }
869                                   
870                                    //
871                                    // Scale if necessary
872                                    //
873                                    if( (double)(scl)!=(double)(1) )
874                                    {
875                                        for(i_=1+n; i_<=ki+n;i_++)
876                                        {
877                                            work[i_] = scl*work[i_];
878                                        }
879                                        for(i_=1+n2; i_<=ki+n2;i_++)
880                                        {
881                                            work[i_] = scl*work[i_];
882                                        }
883                                    }
884                                    work[j-1+n] = x[1,1];
885                                    work[j+n] = x[2,1];
886                                    work[j-1+n2] = x[1,2];
887                                    work[j+n2] = x[2,2];
888                                   
889                                    //
890                                    // Update the right-hand side
891                                    //
892                                    vt = -x[1,1];
893                                    i1_ = (1) - (n+1);
894                                    for(i_=n+1; i_<=n+j-2;i_++)
895                                    {
896                                        work[i_] = work[i_] + vt*t[i_+i1_,j-1];
897                                    }
898                                    vt = -x[2,1];
899                                    i1_ = (1) - (n+1);
900                                    for(i_=n+1; i_<=n+j-2;i_++)
901                                    {
902                                        work[i_] = work[i_] + vt*t[i_+i1_,j];
903                                    }
904                                    vt = -x[1,2];
905                                    i1_ = (1) - (n2+1);
906                                    for(i_=n2+1; i_<=n2+j-2;i_++)
907                                    {
908                                        work[i_] = work[i_] + vt*t[i_+i1_,j-1];
909                                    }
910                                    vt = -x[2,2];
911                                    i1_ = (1) - (n2+1);
912                                    for(i_=n2+1; i_<=n2+j-2;i_++)
913                                    {
914                                        work[i_] = work[i_] + vt*t[i_+i1_,j];
915                                    }
916                                }
917                            }
918                           
919                            //
920                            // Copy the vector x or Q*x to VR and normalize.
921                            //
922                            if( !over )
923                            {
924                                i1_ = (n+1) - (1);
925                                for(i_=1; i_<=ki;i_++)
926                                {
927                                    vr[i_,iis-1] = work[i_+i1_];
928                                }
929                                i1_ = (n2+1) - (1);
930                                for(i_=1; i_<=ki;i_++)
931                                {
932                                    vr[i_,iis] = work[i_+i1_];
933                                }
934                                emax = 0;
935                                for(k=1; k<=ki; k++)
936                                {
937                                    emax = Math.Max(emax, Math.Abs(vr[k,iis-1])+Math.Abs(vr[k,iis]));
938                                }
939                                remax = 1/emax;
940                                for(i_=1; i_<=ki;i_++)
941                                {
942                                    vr[i_,iis-1] = remax*vr[i_,iis-1];
943                                }
944                                for(i_=1; i_<=ki;i_++)
945                                {
946                                    vr[i_,iis] = remax*vr[i_,iis];
947                                }
948                                for(k=ki+1; k<=n; k++)
949                                {
950                                    vr[k,iis-1] = 0;
951                                    vr[k,iis] = 0;
952                                }
953                            }
954                            else
955                            {
956                                if( ki>2 )
957                                {
958                                    for(i_=1; i_<=n;i_++)
959                                    {
960                                        temp[i_] = vr[i_,ki-1];
961                                    }
962                                    blas.matrixvectormultiply(ref vr, 1, n, 1, ki-2, false, ref work, 1+n, ki-2+n, 1.0, ref temp, 1, n, work[ki-1+n]);
963                                    for(i_=1; i_<=n;i_++)
964                                    {
965                                        vr[i_,ki-1] = temp[i_];
966                                    }
967                                    for(i_=1; i_<=n;i_++)
968                                    {
969                                        temp[i_] = vr[i_,ki];
970                                    }
971                                    blas.matrixvectormultiply(ref vr, 1, n, 1, ki-2, false, ref work, 1+n2, ki-2+n2, 1.0, ref temp, 1, n, work[ki+n2]);
972                                    for(i_=1; i_<=n;i_++)
973                                    {
974                                        vr[i_,ki] = temp[i_];
975                                    }
976                                }
977                                else
978                                {
979                                    vt = work[ki-1+n];
980                                    for(i_=1; i_<=n;i_++)
981                                    {
982                                        vr[i_,ki-1] = vt*vr[i_,ki-1];
983                                    }
984                                    vt = work[ki+n2];
985                                    for(i_=1; i_<=n;i_++)
986                                    {
987                                        vr[i_,ki] = vt*vr[i_,ki];
988                                    }
989                                }
990                                emax = 0;
991                                for(k=1; k<=n; k++)
992                                {
993                                    emax = Math.Max(emax, Math.Abs(vr[k,ki-1])+Math.Abs(vr[k,ki]));
994                                }
995                                remax = 1/emax;
996                                for(i_=1; i_<=n;i_++)
997                                {
998                                    vr[i_,ki-1] = remax*vr[i_,ki-1];
999                                }
1000                                for(i_=1; i_<=n;i_++)
1001                                {
1002                                    vr[i_,ki] = remax*vr[i_,ki];
1003                                }
1004                            }
1005                        }
1006                        iis = iis-1;
1007                        if( ip!=0 )
1008                        {
1009                            iis = iis-1;
1010                        }
1011                    }
1012                    if( ip==1 )
1013                    {
1014                        ip = 0;
1015                    }
1016                    if( ip==-1 )
1017                    {
1018                        ip = 1;
1019                    }
1020                }
1021            }
1022            if( leftv )
1023            {
1024               
1025                //
1026                // Compute left eigenvectors.
1027                //
1028                ip = 0;
1029                iis = 1;
1030                for(ki=1; ki<=n; ki++)
1031                {
1032                    skipflag = false;
1033                    if( ip==-1 )
1034                    {
1035                        skipflag = true;
1036                    }
1037                    else
1038                    {
1039                        if( ki!=n )
1040                        {
1041                            if( (double)(t[ki+1,ki])!=(double)(0) )
1042                            {
1043                                ip = 1;
1044                            }
1045                        }
1046                        if( somev )
1047                        {
1048                            if( !vselect[ki] )
1049                            {
1050                                skipflag = true;
1051                            }
1052                        }
1053                    }
1054                    if( !skipflag )
1055                    {
1056                       
1057                        //
1058                        // Compute the KI-th eigenvalue (WR,WI).
1059                        //
1060                        wr = t[ki,ki];
1061                        wi = 0;
1062                        if( ip!=0 )
1063                        {
1064                            wi = Math.Sqrt(Math.Abs(t[ki,ki+1]))*Math.Sqrt(Math.Abs(t[ki+1,ki]));
1065                        }
1066                        smin = Math.Max(ulp*(Math.Abs(wr)+Math.Abs(wi)), smlnum);
1067                        if( ip==0 )
1068                        {
1069                           
1070                            //
1071                            // Real left eigenvector.
1072                            //
1073                            work[ki+n] = 1;
1074                           
1075                            //
1076                            // Form right-hand side
1077                            //
1078                            for(k=ki+1; k<=n; k++)
1079                            {
1080                                work[k+n] = -t[ki,k];
1081                            }
1082                           
1083                            //
1084                            // Solve the quasi-triangular system:
1085                            // (T(KI+1:N,KI+1:N) - WR)'*X = SCALE*WORK
1086                            //
1087                            vmax = 1;
1088                            vcrit = bignum;
1089                            jnxt = ki+1;
1090                            for(j=ki+1; j<=n; j++)
1091                            {
1092                                if( j<jnxt )
1093                                {
1094                                    continue;
1095                                }
1096                                j1 = j;
1097                                j2 = j;
1098                                jnxt = j+1;
1099                                if( j<n )
1100                                {
1101                                    if( (double)(t[j+1,j])!=(double)(0) )
1102                                    {
1103                                        j2 = j+1;
1104                                        jnxt = j+2;
1105                                    }
1106                                }
1107                                if( j1==j2 )
1108                                {
1109                                   
1110                                    //
1111                                    // 1-by-1 diagonal block
1112                                    //
1113                                    // Scale if necessary to avoid overflow when forming
1114                                    // the right-hand side.
1115                                    //
1116                                    if( (double)(work[j])>(double)(vcrit) )
1117                                    {
1118                                        rec = 1/vmax;
1119                                        for(i_=ki+n; i_<=n+n;i_++)
1120                                        {
1121                                            work[i_] = rec*work[i_];
1122                                        }
1123                                        vmax = 1;
1124                                        vcrit = bignum;
1125                                    }
1126                                    i1_ = (ki+1+n)-(ki+1);
1127                                    vt = 0.0;
1128                                    for(i_=ki+1; i_<=j-1;i_++)
1129                                    {
1130                                        vt += t[i_,j]*work[i_+i1_];
1131                                    }
1132                                    work[j+n] = work[j+n]-vt;
1133                                   
1134                                    //
1135                                    // Solve (T(J,J)-WR)'*X = WORK
1136                                    //
1137                                    temp11[1,1] = t[j,j];
1138                                    temp11b[1,1] = work[j+n];
1139                                    internalhsevdlaln2(false, 1, 1, smin, 1.0, ref temp11, 1.0, 1.0, ref temp11b, wr, 0, ref rswap4, ref zswap4, ref ipivot44, ref civ4, ref crv4, ref x, ref scl, ref xnorm, ref ierr);
1140                                   
1141                                    //
1142                                    // Scale if necessary
1143                                    //
1144                                    if( (double)(scl)!=(double)(1) )
1145                                    {
1146                                        for(i_=ki+n; i_<=n+n;i_++)
1147                                        {
1148                                            work[i_] = scl*work[i_];
1149                                        }
1150                                    }
1151                                    work[j+n] = x[1,1];
1152                                    vmax = Math.Max(Math.Abs(work[j+n]), vmax);
1153                                    vcrit = bignum/vmax;
1154                                }
1155                                else
1156                                {
1157                                   
1158                                    //
1159                                    // 2-by-2 diagonal block
1160                                    //
1161                                    // Scale if necessary to avoid overflow when forming
1162                                    // the right-hand side.
1163                                    //
1164                                    beta = Math.Max(work[j], work[j+1]);
1165                                    if( (double)(beta)>(double)(vcrit) )
1166                                    {
1167                                        rec = 1/vmax;
1168                                        for(i_=ki+n; i_<=n+n;i_++)
1169                                        {
1170                                            work[i_] = rec*work[i_];
1171                                        }
1172                                        vmax = 1;
1173                                        vcrit = bignum;
1174                                    }
1175                                    i1_ = (ki+1+n)-(ki+1);
1176                                    vt = 0.0;
1177                                    for(i_=ki+1; i_<=j-1;i_++)
1178                                    {
1179                                        vt += t[i_,j]*work[i_+i1_];
1180                                    }
1181                                    work[j+n] = work[j+n]-vt;
1182                                    i1_ = (ki+1+n)-(ki+1);
1183                                    vt = 0.0;
1184                                    for(i_=ki+1; i_<=j-1;i_++)
1185                                    {
1186                                        vt += t[i_,j+1]*work[i_+i1_];
1187                                    }
1188                                    work[j+1+n] = work[j+1+n]-vt;
1189                                   
1190                                    //
1191                                    // Solve
1192                                    //    [T(J,J)-WR   T(J,J+1)     ]'* X = SCALE*( WORK1 )
1193                                    //    [T(J+1,J)    T(J+1,J+1)-WR]             ( WORK2 )
1194                                    //
1195                                    temp22[1,1] = t[j,j];
1196                                    temp22[1,2] = t[j,j+1];
1197                                    temp22[2,1] = t[j+1,j];
1198                                    temp22[2,2] = t[j+1,j+1];
1199                                    temp21b[1,1] = work[j+n];
1200                                    temp21b[2,1] = work[j+1+n];
1201                                    internalhsevdlaln2(true, 2, 1, smin, 1.0, ref temp22, 1.0, 1.0, ref temp21b, wr, 0, ref rswap4, ref zswap4, ref ipivot44, ref civ4, ref crv4, ref x, ref scl, ref xnorm, ref ierr);
1202                                   
1203                                    //
1204                                    // Scale if necessary
1205                                    //
1206                                    if( (double)(scl)!=(double)(1) )
1207                                    {
1208                                        for(i_=ki+n; i_<=n+n;i_++)
1209                                        {
1210                                            work[i_] = scl*work[i_];
1211                                        }
1212                                    }
1213                                    work[j+n] = x[1,1];
1214                                    work[j+1+n] = x[2,1];
1215                                    vmax = Math.Max(Math.Abs(work[j+n]), Math.Max(Math.Abs(work[j+1+n]), vmax));
1216                                    vcrit = bignum/vmax;
1217                                }
1218                            }
1219                           
1220                            //
1221                            // Copy the vector x or Q*x to VL and normalize.
1222                            //
1223                            if( !over )
1224                            {
1225                                i1_ = (ki+n) - (ki);
1226                                for(i_=ki; i_<=n;i_++)
1227                                {
1228                                    vl[i_,iis] = work[i_+i1_];
1229                                }
1230                                ii = blas.columnidxabsmax(ref vl, ki, n, iis);
1231                                remax = 1/Math.Abs(vl[ii,iis]);
1232                                for(i_=ki; i_<=n;i_++)
1233                                {
1234                                    vl[i_,iis] = remax*vl[i_,iis];
1235                                }
1236                                for(k=1; k<=ki-1; k++)
1237                                {
1238                                    vl[k,iis] = 0;
1239                                }
1240                            }
1241                            else
1242                            {
1243                                if( ki<n )
1244                                {
1245                                    for(i_=1; i_<=n;i_++)
1246                                    {
1247                                        temp[i_] = vl[i_,ki];
1248                                    }
1249                                    blas.matrixvectormultiply(ref vl, 1, n, ki+1, n, false, ref work, ki+1+n, n+n, 1.0, ref temp, 1, n, work[ki+n]);
1250                                    for(i_=1; i_<=n;i_++)
1251                                    {
1252                                        vl[i_,ki] = temp[i_];
1253                                    }
1254                                }
1255                                ii = blas.columnidxabsmax(ref vl, 1, n, ki);
1256                                remax = 1/Math.Abs(vl[ii,ki]);
1257                                for(i_=1; i_<=n;i_++)
1258                                {
1259                                    vl[i_,ki] = remax*vl[i_,ki];
1260                                }
1261                            }
1262                        }
1263                        else
1264                        {
1265                           
1266                            //
1267                            // Complex left eigenvector.
1268                            //
1269                            // Initial solve:
1270                            //   ((T(KI,KI)    T(KI,KI+1) )' - (WR - I* WI))*X = 0.
1271                            //   ((T(KI+1,KI) T(KI+1,KI+1))                )
1272                            //
1273                            if( (double)(Math.Abs(t[ki,ki+1]))>=(double)(Math.Abs(t[ki+1,ki])) )
1274                            {
1275                                work[ki+n] = wi/t[ki,ki+1];
1276                                work[ki+1+n2] = 1;
1277                            }
1278                            else
1279                            {
1280                                work[ki+n] = 1;
1281                                work[ki+1+n2] = -(wi/t[ki+1,ki]);
1282                            }
1283                            work[ki+1+n] = 0;
1284                            work[ki+n2] = 0;
1285                           
1286                            //
1287                            // Form right-hand side
1288                            //
1289                            for(k=ki+2; k<=n; k++)
1290                            {
1291                                work[k+n] = -(work[ki+n]*t[ki,k]);
1292                                work[k+n2] = -(work[ki+1+n2]*t[ki+1,k]);
1293                            }
1294                           
1295                            //
1296                            // Solve complex quasi-triangular system:
1297                            // ( T(KI+2,N:KI+2,N) - (WR-i*WI) )*X = WORK1+i*WORK2
1298                            //
1299                            vmax = 1;
1300                            vcrit = bignum;
1301                            jnxt = ki+2;
1302                            for(j=ki+2; j<=n; j++)
1303                            {
1304                                if( j<jnxt )
1305                                {
1306                                    continue;
1307                                }
1308                                j1 = j;
1309                                j2 = j;
1310                                jnxt = j+1;
1311                                if( j<n )
1312                                {
1313                                    if( (double)(t[j+1,j])!=(double)(0) )
1314                                    {
1315                                        j2 = j+1;
1316                                        jnxt = j+2;
1317                                    }
1318                                }
1319                                if( j1==j2 )
1320                                {
1321                                   
1322                                    //
1323                                    // 1-by-1 diagonal block
1324                                    //
1325                                    // Scale if necessary to avoid overflow when
1326                                    // forming the right-hand side elements.
1327                                    //
1328                                    if( (double)(work[j])>(double)(vcrit) )
1329                                    {
1330                                        rec = 1/vmax;
1331                                        for(i_=ki+n; i_<=n+n;i_++)
1332                                        {
1333                                            work[i_] = rec*work[i_];
1334                                        }
1335                                        for(i_=ki+n2; i_<=n+n2;i_++)
1336                                        {
1337                                            work[i_] = rec*work[i_];
1338                                        }
1339                                        vmax = 1;
1340                                        vcrit = bignum;
1341                                    }
1342                                    i1_ = (ki+2+n)-(ki+2);
1343                                    vt = 0.0;
1344                                    for(i_=ki+2; i_<=j-1;i_++)
1345                                    {
1346                                        vt += t[i_,j]*work[i_+i1_];
1347                                    }
1348                                    work[j+n] = work[j+n]-vt;
1349                                    i1_ = (ki+2+n2)-(ki+2);
1350                                    vt = 0.0;
1351                                    for(i_=ki+2; i_<=j-1;i_++)
1352                                    {
1353                                        vt += t[i_,j]*work[i_+i1_];
1354                                    }
1355                                    work[j+n2] = work[j+n2]-vt;
1356                                   
1357                                    //
1358                                    // Solve (T(J,J)-(WR-i*WI))*(X11+i*X12)= WK+I*WK2
1359                                    //
1360                                    temp11[1,1] = t[j,j];
1361                                    temp12b[1,1] = work[j+n];
1362                                    temp12b[1,2] = work[j+n+n];
1363                                    internalhsevdlaln2(false, 1, 2, smin, 1.0, ref temp11, 1.0, 1.0, ref temp12b, wr, -wi, ref rswap4, ref zswap4, ref ipivot44, ref civ4, ref crv4, ref x, ref scl, ref xnorm, ref ierr);
1364                                   
1365                                    //
1366                                    // Scale if necessary
1367                                    //
1368                                    if( (double)(scl)!=(double)(1) )
1369                                    {
1370                                        for(i_=ki+n; i_<=n+n;i_++)
1371                                        {
1372                                            work[i_] = scl*work[i_];
1373                                        }
1374                                        for(i_=ki+n2; i_<=n+n2;i_++)
1375                                        {
1376                                            work[i_] = scl*work[i_];
1377                                        }
1378                                    }
1379                                    work[j+n] = x[1,1];
1380                                    work[j+n2] = x[1,2];
1381                                    vmax = Math.Max(Math.Abs(work[j+n]), Math.Max(Math.Abs(work[j+n2]), vmax));
1382                                    vcrit = bignum/vmax;
1383                                }
1384                                else
1385                                {
1386                                   
1387                                    //
1388                                    // 2-by-2 diagonal block
1389                                    //
1390                                    // Scale if necessary to avoid overflow when forming
1391                                    // the right-hand side elements.
1392                                    //
1393                                    beta = Math.Max(work[j], work[j+1]);
1394                                    if( (double)(beta)>(double)(vcrit) )
1395                                    {
1396                                        rec = 1/vmax;
1397                                        for(i_=ki+n; i_<=n+n;i_++)
1398                                        {
1399                                            work[i_] = rec*work[i_];
1400                                        }
1401                                        for(i_=ki+n2; i_<=n+n2;i_++)
1402                                        {
1403                                            work[i_] = rec*work[i_];
1404                                        }
1405                                        vmax = 1;
1406                                        vcrit = bignum;
1407                                    }
1408                                    i1_ = (ki+2+n)-(ki+2);
1409                                    vt = 0.0;
1410                                    for(i_=ki+2; i_<=j-1;i_++)
1411                                    {
1412                                        vt += t[i_,j]*work[i_+i1_];
1413                                    }
1414                                    work[j+n] = work[j+n]-vt;
1415                                    i1_ = (ki+2+n2)-(ki+2);
1416                                    vt = 0.0;
1417                                    for(i_=ki+2; i_<=j-1;i_++)
1418                                    {
1419                                        vt += t[i_,j]*work[i_+i1_];
1420                                    }
1421                                    work[j+n2] = work[j+n2]-vt;
1422                                    i1_ = (ki+2+n)-(ki+2);
1423                                    vt = 0.0;
1424                                    for(i_=ki+2; i_<=j-1;i_++)
1425                                    {
1426                                        vt += t[i_,j+1]*work[i_+i1_];
1427                                    }
1428                                    work[j+1+n] = work[j+1+n]-vt;
1429                                    i1_ = (ki+2+n2)-(ki+2);
1430                                    vt = 0.0;
1431                                    for(i_=ki+2; i_<=j-1;i_++)
1432                                    {
1433                                        vt += t[i_,j+1]*work[i_+i1_];
1434                                    }
1435                                    work[j+1+n2] = work[j+1+n2]-vt;
1436                                   
1437                                    //
1438                                    // Solve 2-by-2 complex linear equation
1439                                    //   ([T(j,j)   T(j,j+1)  ]'-(wr-i*wi)*I)*X = SCALE*B
1440                                    //   ([T(j+1,j) T(j+1,j+1)]             )
1441                                    //
1442                                    temp22[1,1] = t[j,j];
1443                                    temp22[1,2] = t[j,j+1];
1444                                    temp22[2,1] = t[j+1,j];
1445                                    temp22[2,2] = t[j+1,j+1];
1446                                    temp22b[1,1] = work[j+n];
1447                                    temp22b[1,2] = work[j+n+n];
1448                                    temp22b[2,1] = work[j+1+n];
1449                                    temp22b[2,2] = work[j+1+n+n];
1450                                    internalhsevdlaln2(true, 2, 2, smin, 1.0, ref temp22, 1.0, 1.0, ref temp22b, wr, -wi, ref rswap4, ref zswap4, ref ipivot44, ref civ4, ref crv4, ref x, ref scl, ref xnorm, ref ierr);
1451                                   
1452                                    //
1453                                    // Scale if necessary
1454                                    //
1455                                    if( (double)(scl)!=(double)(1) )
1456                                    {
1457                                        for(i_=ki+n; i_<=n+n;i_++)
1458                                        {
1459                                            work[i_] = scl*work[i_];
1460                                        }
1461                                        for(i_=ki+n2; i_<=n+n2;i_++)
1462                                        {
1463                                            work[i_] = scl*work[i_];
1464                                        }
1465                                    }
1466                                    work[j+n] = x[1,1];
1467                                    work[j+n2] = x[1,2];
1468                                    work[j+1+n] = x[2,1];
1469                                    work[j+1+n2] = x[2,2];
1470                                    vmax = Math.Max(Math.Abs(x[1,1]), vmax);
1471                                    vmax = Math.Max(Math.Abs(x[1,2]), vmax);
1472                                    vmax = Math.Max(Math.Abs(x[2,1]), vmax);
1473                                    vmax = Math.Max(Math.Abs(x[2,2]), vmax);
1474                                    vcrit = bignum/vmax;
1475                                }
1476                            }
1477                           
1478                            //
1479                            // Copy the vector x or Q*x to VL and normalize.
1480                            //
1481                            if( !over )
1482                            {
1483                                i1_ = (ki+n) - (ki);
1484                                for(i_=ki; i_<=n;i_++)
1485                                {
1486                                    vl[i_,iis] = work[i_+i1_];
1487                                }
1488                                i1_ = (ki+n2) - (ki);
1489                                for(i_=ki; i_<=n;i_++)
1490                                {
1491                                    vl[i_,iis+1] = work[i_+i1_];
1492                                }
1493                                emax = 0;
1494                                for(k=ki; k<=n; k++)
1495                                {
1496                                    emax = Math.Max(emax, Math.Abs(vl[k,iis])+Math.Abs(vl[k,iis+1]));
1497                                }
1498                                remax = 1/emax;
1499                                for(i_=ki; i_<=n;i_++)
1500                                {
1501                                    vl[i_,iis] = remax*vl[i_,iis];
1502                                }
1503                                for(i_=ki; i_<=n;i_++)
1504                                {
1505                                    vl[i_,iis+1] = remax*vl[i_,iis+1];
1506                                }
1507                                for(k=1; k<=ki-1; k++)
1508                                {
1509                                    vl[k,iis] = 0;
1510                                    vl[k,iis+1] = 0;
1511                                }
1512                            }
1513                            else
1514                            {
1515                                if( ki<n-1 )
1516                                {
1517                                    for(i_=1; i_<=n;i_++)
1518                                    {
1519                                        temp[i_] = vl[i_,ki];
1520                                    }
1521                                    blas.matrixvectormultiply(ref vl, 1, n, ki+2, n, false, ref work, ki+2+n, n+n, 1.0, ref temp, 1, n, work[ki+n]);
1522                                    for(i_=1; i_<=n;i_++)
1523                                    {
1524                                        vl[i_,ki] = temp[i_];
1525                                    }
1526                                    for(i_=1; i_<=n;i_++)
1527                                    {
1528                                        temp[i_] = vl[i_,ki+1];
1529                                    }
1530                                    blas.matrixvectormultiply(ref vl, 1, n, ki+2, n, false, ref work, ki+2+n2, n+n2, 1.0, ref temp, 1, n, work[ki+1+n2]);
1531                                    for(i_=1; i_<=n;i_++)
1532                                    {
1533                                        vl[i_,ki+1] = temp[i_];
1534                                    }
1535                                }
1536                                else
1537                                {
1538                                    vt = work[ki+n];
1539                                    for(i_=1; i_<=n;i_++)
1540                                    {
1541                                        vl[i_,ki] = vt*vl[i_,ki];
1542                                    }
1543                                    vt = work[ki+1+n2];
1544                                    for(i_=1; i_<=n;i_++)
1545                                    {
1546                                        vl[i_,ki+1] = vt*vl[i_,ki+1];
1547                                    }
1548                                }
1549                                emax = 0;
1550                                for(k=1; k<=n; k++)
1551                                {
1552                                    emax = Math.Max(emax, Math.Abs(vl[k,ki])+Math.Abs(vl[k,ki+1]));
1553                                }
1554                                remax = 1/emax;
1555                                for(i_=1; i_<=n;i_++)
1556                                {
1557                                    vl[i_,ki] = remax*vl[i_,ki];
1558                                }
1559                                for(i_=1; i_<=n;i_++)
1560                                {
1561                                    vl[i_,ki+1] = remax*vl[i_,ki+1];
1562                                }
1563                            }
1564                        }
1565                        iis = iis+1;
1566                        if( ip!=0 )
1567                        {
1568                            iis = iis+1;
1569                        }
1570                    }
1571                    if( ip==-1 )
1572                    {
1573                        ip = 0;
1574                    }
1575                    if( ip==1 )
1576                    {
1577                        ip = -1;
1578                    }
1579                }
1580            }
1581        }
1582
1583
1584        private static void internalhsevdlaln2(bool ltrans,
1585            int na,
1586            int nw,
1587            double smin,
1588            double ca,
1589            ref double[,] a,
1590            double d1,
1591            double d2,
1592            ref double[,] b,
1593            double wr,
1594            double wi,
1595            ref bool[] rswap4,
1596            ref bool[] zswap4,
1597            ref int[,] ipivot44,
1598            ref double[] civ4,
1599            ref double[] crv4,
1600            ref double[,] x,
1601            ref double scl,
1602            ref double xnorm,
1603            ref int info)
1604        {
1605            int icmax = 0;
1606            int j = 0;
1607            double bbnd = 0;
1608            double bi1 = 0;
1609            double bi2 = 0;
1610            double bignum = 0;
1611            double bnorm = 0;
1612            double br1 = 0;
1613            double br2 = 0;
1614            double ci21 = 0;
1615            double ci22 = 0;
1616            double cmax = 0;
1617            double cnorm = 0;
1618            double cr21 = 0;
1619            double cr22 = 0;
1620            double csi = 0;
1621            double csr = 0;
1622            double li21 = 0;
1623            double lr21 = 0;
1624            double smini = 0;
1625            double smlnum = 0;
1626            double temp = 0;
1627            double u22abs = 0;
1628            double ui11 = 0;
1629            double ui11r = 0;
1630            double ui12 = 0;
1631            double ui12s = 0;
1632            double ui22 = 0;
1633            double ur11 = 0;
1634            double ur11r = 0;
1635            double ur12 = 0;
1636            double ur12s = 0;
1637            double ur22 = 0;
1638            double xi1 = 0;
1639            double xi2 = 0;
1640            double xr1 = 0;
1641            double xr2 = 0;
1642            double tmp1 = 0;
1643            double tmp2 = 0;
1644
1645            zswap4[1] = false;
1646            zswap4[2] = false;
1647            zswap4[3] = true;
1648            zswap4[4] = true;
1649            rswap4[1] = false;
1650            rswap4[2] = true;
1651            rswap4[3] = false;
1652            rswap4[4] = true;
1653            ipivot44[1,1] = 1;
1654            ipivot44[2,1] = 2;
1655            ipivot44[3,1] = 3;
1656            ipivot44[4,1] = 4;
1657            ipivot44[1,2] = 2;
1658            ipivot44[2,2] = 1;
1659            ipivot44[3,2] = 4;
1660            ipivot44[4,2] = 3;
1661            ipivot44[1,3] = 3;
1662            ipivot44[2,3] = 4;
1663            ipivot44[3,3] = 1;
1664            ipivot44[4,3] = 2;
1665            ipivot44[1,4] = 4;
1666            ipivot44[2,4] = 3;
1667            ipivot44[3,4] = 2;
1668            ipivot44[4,4] = 1;
1669            smlnum = 2*AP.Math.MinRealNumber;
1670            bignum = 1/smlnum;
1671            smini = Math.Max(smin, smlnum);
1672           
1673            //
1674            // Don't check for input errors
1675            //
1676            info = 0;
1677           
1678            //
1679            // Standard Initializations
1680            //
1681            scl = 1;
1682            if( na==1 )
1683            {
1684               
1685                //
1686                // 1 x 1  (i.e., scalar) system   C X = B
1687                //
1688                if( nw==1 )
1689                {
1690                   
1691                    //
1692                    // Real 1x1 system.
1693                    //
1694                    // C = ca A - w D
1695                    //
1696                    csr = ca*a[1,1]-wr*d1;
1697                    cnorm = Math.Abs(csr);
1698                   
1699                    //
1700                    // If | C | < SMINI, use C = SMINI
1701                    //
1702                    if( (double)(cnorm)<(double)(smini) )
1703                    {
1704                        csr = smini;
1705                        cnorm = smini;
1706                        info = 1;
1707                    }
1708                   
1709                    //
1710                    // Check scaling for  X = B / C
1711                    //
1712                    bnorm = Math.Abs(b[1,1]);
1713                    if( (double)(cnorm)<(double)(1) & (double)(bnorm)>(double)(1) )
1714                    {
1715                        if( (double)(bnorm)>(double)(bignum*cnorm) )
1716                        {
1717                            scl = 1/bnorm;
1718                        }
1719                    }
1720                   
1721                    //
1722                    // Compute X
1723                    //
1724                    x[1,1] = b[1,1]*scl/csr;
1725                    xnorm = Math.Abs(x[1,1]);
1726                }
1727                else
1728                {
1729                   
1730                    //
1731                    // Complex 1x1 system (w is complex)
1732                    //
1733                    // C = ca A - w D
1734                    //
1735                    csr = ca*a[1,1]-wr*d1;
1736                    csi = -(wi*d1);
1737                    cnorm = Math.Abs(csr)+Math.Abs(csi);
1738                   
1739                    //
1740                    // If | C | < SMINI, use C = SMINI
1741                    //
1742                    if( (double)(cnorm)<(double)(smini) )
1743                    {
1744                        csr = smini;
1745                        csi = 0;
1746                        cnorm = smini;
1747                        info = 1;
1748                    }
1749                   
1750                    //
1751                    // Check scaling for  X = B / C
1752                    //
1753                    bnorm = Math.Abs(b[1,1])+Math.Abs(b[1,2]);
1754                    if( (double)(cnorm)<(double)(1) & (double)(bnorm)>(double)(1) )
1755                    {
1756                        if( (double)(bnorm)>(double)(bignum*cnorm) )
1757                        {
1758                            scl = 1/bnorm;
1759                        }
1760                    }
1761                   
1762                    //
1763                    // Compute X
1764                    //
1765                    internalhsevdladiv(scl*b[1,1], scl*b[1,2], csr, csi, ref tmp1, ref tmp2);
1766                    x[1,1] = tmp1;
1767                    x[1,2] = tmp2;
1768                    xnorm = Math.Abs(x[1,1])+Math.Abs(x[1,2]);
1769                }
1770            }
1771            else
1772            {
1773               
1774                //
1775                // 2x2 System
1776                //
1777                // Compute the real part of  C = ca A - w D  (or  ca A' - w D )
1778                //
1779                crv4[1+0] = ca*a[1,1]-wr*d1;
1780                crv4[2+2] = ca*a[2,2]-wr*d2;
1781                if( ltrans )
1782                {
1783                    crv4[1+2] = ca*a[2,1];
1784                    crv4[2+0] = ca*a[1,2];
1785                }
1786                else
1787                {
1788                    crv4[2+0] = ca*a[2,1];
1789                    crv4[1+2] = ca*a[1,2];
1790                }
1791                if( nw==1 )
1792                {
1793                   
1794                    //
1795                    // Real 2x2 system  (w is real)
1796                    //
1797                    // Find the largest element in C
1798                    //
1799                    cmax = 0;
1800                    icmax = 0;
1801                    for(j=1; j<=4; j++)
1802                    {
1803                        if( (double)(Math.Abs(crv4[j]))>(double)(cmax) )
1804                        {
1805                            cmax = Math.Abs(crv4[j]);
1806                            icmax = j;
1807                        }
1808                    }
1809                   
1810                    //
1811                    // If norm(C) < SMINI, use SMINI*identity.
1812                    //
1813                    if( (double)(cmax)<(double)(smini) )
1814                    {
1815                        bnorm = Math.Max(Math.Abs(b[1,1]), Math.Abs(b[2,1]));
1816                        if( (double)(smini)<(double)(1) & (double)(bnorm)>(double)(1) )
1817                        {
1818                            if( (double)(bnorm)>(double)(bignum*smini) )
1819                            {
1820                                scl = 1/bnorm;
1821                            }
1822                        }
1823                        temp = scl/smini;
1824                        x[1,1] = temp*b[1,1];
1825                        x[2,1] = temp*b[2,1];
1826                        xnorm = temp*bnorm;
1827                        info = 1;
1828                        return;
1829                    }
1830                   
1831                    //
1832                    // Gaussian elimination with complete pivoting.
1833                    //
1834                    ur11 = crv4[icmax];
1835                    cr21 = crv4[ipivot44[2,icmax]];
1836                    ur12 = crv4[ipivot44[3,icmax]];
1837                    cr22 = crv4[ipivot44[4,icmax]];
1838                    ur11r = 1/ur11;
1839                    lr21 = ur11r*cr21;
1840                    ur22 = cr22-ur12*lr21;
1841                   
1842                    //
1843                    // If smaller pivot < SMINI, use SMINI
1844                    //
1845                    if( (double)(Math.Abs(ur22))<(double)(smini) )
1846                    {
1847                        ur22 = smini;
1848                        info = 1;
1849                    }
1850                    if( rswap4[icmax] )
1851                    {
1852                        br1 = b[2,1];
1853                        br2 = b[1,1];
1854                    }
1855                    else
1856                    {
1857                        br1 = b[1,1];
1858                        br2 = b[2,1];
1859                    }
1860                    br2 = br2-lr21*br1;
1861                    bbnd = Math.Max(Math.Abs(br1*(ur22*ur11r)), Math.Abs(br2));
1862                    if( (double)(bbnd)>(double)(1) & (double)(Math.Abs(ur22))<(double)(1) )
1863                    {
1864                        if( (double)(bbnd)>=(double)(bignum*Math.Abs(ur22)) )
1865                        {
1866                            scl = 1/bbnd;
1867                        }
1868                    }
1869                    xr2 = br2*scl/ur22;
1870                    xr1 = scl*br1*ur11r-xr2*(ur11r*ur12);
1871                    if( zswap4[icmax] )
1872                    {
1873                        x[1,1] = xr2;
1874                        x[2,1] = xr1;
1875                    }
1876                    else
1877                    {
1878                        x[1,1] = xr1;
1879                        x[2,1] = xr2;
1880                    }
1881                    xnorm = Math.Max(Math.Abs(xr1), Math.Abs(xr2));
1882                   
1883                    //
1884                    // Further scaling if  norm(A) norm(X) > overflow
1885                    //
1886                    if( (double)(xnorm)>(double)(1) & (double)(cmax)>(double)(1) )
1887                    {
1888                        if( (double)(xnorm)>(double)(bignum/cmax) )
1889                        {
1890                            temp = cmax/bignum;
1891                            x[1,1] = temp*x[1,1];
1892                            x[2,1] = temp*x[2,1];
1893                            xnorm = temp*xnorm;
1894                            scl = temp*scl;
1895                        }
1896                    }
1897                }
1898                else
1899                {
1900                   
1901                    //
1902                    // Complex 2x2 system  (w is complex)
1903                    //
1904                    // Find the largest element in C
1905                    //
1906                    civ4[1+0] = -(wi*d1);
1907                    civ4[2+0] = 0;
1908                    civ4[1+2] = 0;
1909                    civ4[2+2] = -(wi*d2);
1910                    cmax = 0;
1911                    icmax = 0;
1912                    for(j=1; j<=4; j++)
1913                    {
1914                        if( (double)(Math.Abs(crv4[j])+Math.Abs(civ4[j]))>(double)(cmax) )
1915                        {
1916                            cmax = Math.Abs(crv4[j])+Math.Abs(civ4[j]);
1917                            icmax = j;
1918                        }
1919                    }
1920                   
1921                    //
1922                    // If norm(C) < SMINI, use SMINI*identity.
1923                    //
1924                    if( (double)(cmax)<(double)(smini) )
1925                    {
1926                        bnorm = Math.Max(Math.Abs(b[1,1])+Math.Abs(b[1,2]), Math.Abs(b[2,1])+Math.Abs(b[2,2]));
1927                        if( (double)(smini)<(double)(1) & (double)(bnorm)>(double)(1) )
1928                        {
1929                            if( (double)(bnorm)>(double)(bignum*smini) )
1930                            {
1931                                scl = 1/bnorm;
1932                            }
1933                        }
1934                        temp = scl/smini;
1935                        x[1,1] = temp*b[1,1];
1936                        x[2,1] = temp*b[2,1];
1937                        x[1,2] = temp*b[1,2];
1938                        x[2,2] = temp*b[2,2];
1939                        xnorm = temp*bnorm;
1940                        info = 1;
1941                        return;
1942                    }
1943                   
1944                    //
1945                    // Gaussian elimination with complete pivoting.
1946                    //
1947                    ur11 = crv4[icmax];
1948                    ui11 = civ4[icmax];
1949                    cr21 = crv4[ipivot44[2,icmax]];
1950                    ci21 = civ4[ipivot44[2,icmax]];
1951                    ur12 = crv4[ipivot44[3,icmax]];
1952                    ui12 = civ4[ipivot44[3,icmax]];
1953                    cr22 = crv4[ipivot44[4,icmax]];
1954                    ci22 = civ4[ipivot44[4,icmax]];
1955                    if( icmax==1 | icmax==4 )
1956                    {
1957                       
1958                        //
1959                        // Code when off-diagonals of pivoted C are real
1960                        //
1961                        if( (double)(Math.Abs(ur11))>(double)(Math.Abs(ui11)) )
1962                        {
1963                            temp = ui11/ur11;
1964                            ur11r = 1/(ur11*(1+AP.Math.Sqr(temp)));
1965                            ui11r = -(temp*ur11r);
1966                        }
1967                        else
1968                        {
1969                            temp = ur11/ui11;
1970                            ui11r = -(1/(ui11*(1+AP.Math.Sqr(temp))));
1971                            ur11r = -(temp*ui11r);
1972                        }
1973                        lr21 = cr21*ur11r;
1974                        li21 = cr21*ui11r;
1975                        ur12s = ur12*ur11r;
1976                        ui12s = ur12*ui11r;
1977                        ur22 = cr22-ur12*lr21;
1978                        ui22 = ci22-ur12*li21;
1979                    }
1980                    else
1981                    {
1982                       
1983                        //
1984                        // Code when diagonals of pivoted C are real
1985                        //
1986                        ur11r = 1/ur11;
1987                        ui11r = 0;
1988                        lr21 = cr21*ur11r;
1989                        li21 = ci21*ur11r;
1990                        ur12s = ur12*ur11r;
1991                        ui12s = ui12*ur11r;
1992                        ur22 = cr22-ur12*lr21+ui12*li21;
1993                        ui22 = -(ur12*li21)-ui12*lr21;
1994                    }
1995                    u22abs = Math.Abs(ur22)+Math.Abs(ui22);
1996                   
1997                    //
1998                    // If smaller pivot < SMINI, use SMINI
1999                    //
2000                    if( (double)(u22abs)<(double)(smini) )
2001                    {
2002                        ur22 = smini;
2003                        ui22 = 0;
2004                        info = 1;
2005                    }
2006                    if( rswap4[icmax] )
2007                    {
2008                        br2 = b[1,1];
2009                        br1 = b[2,1];
2010                        bi2 = b[1,2];
2011                        bi1 = b[2,2];
2012                    }
2013                    else
2014                    {
2015                        br1 = b[1,1];
2016                        br2 = b[2,1];
2017                        bi1 = b[1,2];
2018                        bi2 = b[2,2];
2019                    }
2020                    br2 = br2-lr21*br1+li21*bi1;
2021                    bi2 = bi2-li21*br1-lr21*bi1;
2022                    bbnd = Math.Max((Math.Abs(br1)+Math.Abs(bi1))*(u22abs*(Math.Abs(ur11r)+Math.Abs(ui11r))), Math.Abs(br2)+Math.Abs(bi2));
2023                    if( (double)(bbnd)>(double)(1) & (double)(u22abs)<(double)(1) )
2024                    {
2025                        if( (double)(bbnd)>=(double)(bignum*u22abs) )
2026                        {
2027                            scl = 1/bbnd;
2028                            br1 = scl*br1;
2029                            bi1 = scl*bi1;
2030                            br2 = scl*br2;
2031                            bi2 = scl*bi2;
2032                        }
2033                    }
2034                    internalhsevdladiv(br2, bi2, ur22, ui22, ref xr2, ref xi2);
2035                    xr1 = ur11r*br1-ui11r*bi1-ur12s*xr2+ui12s*xi2;
2036                    xi1 = ui11r*br1+ur11r*bi1-ui12s*xr2-ur12s*xi2;
2037                    if( zswap4[icmax] )
2038                    {
2039                        x[1,1] = xr2;
2040                        x[2,1] = xr1;
2041                        x[1,2] = xi2;
2042                        x[2,2] = xi1;
2043                    }
2044                    else
2045                    {
2046                        x[1,1] = xr1;
2047                        x[2,1] = xr2;
2048                        x[1,2] = xi1;
2049                        x[2,2] = xi2;
2050                    }
2051                    xnorm = Math.Max(Math.Abs(xr1)+Math.Abs(xi1), Math.Abs(xr2)+Math.Abs(xi2));
2052                   
2053                    //
2054                    // Further scaling if  norm(A) norm(X) > overflow
2055                    //
2056                    if( (double)(xnorm)>(double)(1) & (double)(cmax)>(double)(1) )
2057                    {
2058                        if( (double)(xnorm)>(double)(bignum/cmax) )
2059                        {
2060                            temp = cmax/bignum;
2061                            x[1,1] = temp*x[1,1];
2062                            x[2,1] = temp*x[2,1];
2063                            x[1,2] = temp*x[1,2];
2064                            x[2,2] = temp*x[2,2];
2065                            xnorm = temp*xnorm;
2066                            scl = temp*scl;
2067                        }
2068                    }
2069                }
2070            }
2071        }
2072
2073
2074        private static void internalhsevdladiv(double a,
2075            double b,
2076            double c,
2077            double d,
2078            ref double p,
2079            ref double q)
2080        {
2081            double e = 0;
2082            double f = 0;
2083
2084            if( (double)(Math.Abs(d))<(double)(Math.Abs(c)) )
2085            {
2086                e = d/c;
2087                f = c+d*e;
2088                p = (a+b*e)/f;
2089                q = (b-a*e)/f;
2090            }
2091            else
2092            {
2093                e = c/d;
2094                f = d+c*e;
2095                p = (b+a*e)/f;
2096                q = (-a+b*e)/f;
2097            }
2098        }
2099    }
2100}
Note: See TracBrowser for help on using the repository browser.