source: branches/profile/EISPACK/eispack.cxx

Last change on this file was 5390, checked in by westram, 12 years ago
  • TAB-Ex
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 481.8 KB
Line 
1/* eispack.f -- translated by f2c (version 19950110).
2   You must link the resulting object file with the libraries:
3        -lf2c -lm   (in that order)
4*/
5
6#ifdef __cplusplus
7extern "C" {
8#endif
9#include "f2c.h"
10
11/* Table of constant values */
12
13static doublereal c_b141 = 1.;
14static doublereal c_b550 = 0.;
15
16/* Subroutine */ int cdiv_(doublereal *ar, doublereal *ai, doublereal *br, 
17        doublereal *bi, doublereal *cr, doublereal *ci)
18{
19    /* System generated locals */
20    doublereal d_1, d_2;
21
22    /* Local variables */
23    static doublereal s, ais, bis, ars, brs;
24
25
26/*     COMPLEX DIVISION, (CR,CI) = (AR,AI)/(BR,BI) */
27
28    s = abs(*br) + abs(*bi);
29    ars = *ar / s;
30    ais = *ai / s;
31    brs = *br / s;
32    bis = *bi / s;
33/* Computing 2nd power */
34    d_1 = brs;
35/* Computing 2nd power */
36    d_2 = bis;
37    s = d_1 * d_1 + d_2 * d_2;
38    *cr = (ars * brs + ais * bis) / s;
39    *ci = (ais * brs - ars * bis) / s;
40    return 0;
41} /* cdiv_ */
42
43/* Subroutine */ int csroot_(doublereal *xr, doublereal *xi, doublereal *yr, 
44        doublereal *yi)
45{
46    /* Builtin functions */
47    double sqrt(doublereal);
48
49    /* Local variables */
50    static doublereal s, ti, tr;
51    extern doublereal pythag_(doublereal *, doublereal *);
52
53
54/*     (YR,YI) = COMPLEX DSQRT(XR,XI) */
55/*     BRANCH CHOSEN SO THAT YR .GE. 0.0 AND SIGN(YI) .EQ. SIGN(XI) */
56
57    tr = *xr;
58    ti = *xi;
59    s = sqrt((pythag_(&tr, &ti) + abs(tr)) * .5);
60    if (tr >= 0.) {
61        *yr = s;
62    }
63    if (ti < 0.) {
64        s = -s;
65    }
66    if (tr <= 0.) {
67        *yi = s;
68    }
69    if (tr < 0.) {
70        *yr = ti / *yi * .5;
71    }
72    if (tr > 0.) {
73        *yi = ti / *yr * .5;
74    }
75    return 0;
76} /* csroot_ */
77
78doublereal epslon_(doublereal *x)
79{
80    /* System generated locals */
81    doublereal ret_val, d_1;
82
83    /* Local variables */
84    static doublereal a, b, c, eps;
85
86
87/*     ESTIMATE UNIT ROUNDOFF IN QUANTITIES OF SIZE X. */
88
89
90/*     THIS PROGRAM SHOULD FUNCTION PROPERLY ON ALL SYSTEMS */
91/*     SATISFYING THE FOLLOWING TWO ASSUMPTIONS, */
92/*        1.  THE BASE USED IN REPRESENTING FLOATING POINT */
93/*            NUMBERS IS NOT A POWER OF THREE. */
94/*        2.  THE QUANTITY  A  IN STATEMENT 10 IS REPRESENTED TO */
95/*            THE ACCURACY USED IN FLOATING POINT VARIABLES */
96/*            THAT ARE STORED IN MEMORY. */
97/*     THE STATEMENT NUMBER 10 AND THE GO TO 10 ARE INTENDED TO */
98/*     FORCE OPTIMIZING COMPILERS TO GENERATE CODE SATISFYING */
99/*     ASSUMPTION 2. */
100/*     UNDER THESE ASSUMPTIONS, IT SHOULD BE TRUE THAT, */
101/*            A  IS NOT EXACTLY EQUAL TO FOUR-THIRDS, */
102/*            B  HAS A ZERO FOR ITS LAST BIT OR DIGIT, */
103/*            C  IS NOT EXACTLY EQUAL TO ONE, */
104/*            EPS  MEASURES THE SEPARATION OF 1.0 FROM */
105/*                 THE NEXT LARGER FLOATING POINT NUMBER. */
106/*     THE DEVELOPERS OF EISPACK WOULD APPRECIATE BEING INFORMED */
107/*     ABOUT ANY SYSTEMS WHERE THESE ASSUMPTIONS DO NOT HOLD. */
108
109/*     THIS VERSION DATED 4/6/83. */
110
111    a = 1.3333333333333333;
112L10:
113    b = a - 1.;
114    c = b + b + b;
115    eps = (d_1 = c - 1., abs(d_1));
116    if (eps == 0.) {
117        goto L10;
118    }
119    ret_val = eps * abs(*x);
120    return ret_val;
121} /* epslon_ */
122
123doublereal pythag_(doublereal *a, doublereal *b)
124{
125    /* System generated locals */
126    doublereal ret_val, d_1, d_2, d_3;
127
128    /* Local variables */
129    static doublereal p, r, s, t, u;
130
131
132/*     FINDS DSQRT(A**2+B**2) WITHOUT OVERFLOW OR DESTRUCTIVE UNDERFLOW */
133
134/* Computing MAX */
135    d_1 = abs(*a), d_2 = abs(*b);
136    p = max(d_1,d_2);
137    if (p == 0.) {
138        goto L20;
139    }
140/* Computing MIN */
141    d_2 = abs(*a), d_3 = abs(*b);
142/* Computing 2nd power */
143    d_1 = min(d_2,d_3) / p;
144    r = d_1 * d_1;
145L10:
146    t = r + 4.;
147    if (t == 4.) {
148        goto L20;
149    }
150    s = r / t;
151    u = s * 2. + 1.;
152    p = u * p;
153/* Computing 2nd power */
154    d_1 = s / u;
155    r = d_1 * d_1 * r;
156    goto L10;
157L20:
158    ret_val = p;
159    return ret_val;
160} /* pythag_ */
161
162/* Subroutine */ int bakvec_(integer *nm, integer *n, doublereal *t, 
163        doublereal *e, integer *m, doublereal *z, integer *ierr)
164{
165    /* System generated locals */
166    integer t_dim1, t_offset, z_dim1, z_offset, i_1, i_2;
167
168    /* Local variables */
169    static integer i, j;
170
171
172
173/*     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A NONSYMMETRIC */
174/*     TRIDIAGONAL MATRIX BY BACK TRANSFORMING THOSE OF THE */
175/*     CORRESPONDING SYMMETRIC MATRIX DETERMINED BY  FIGI. */
176
177/*     ON INPUT */
178
179/*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
180/*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
181/*          DIMENSION STATEMENT. */
182
183/*        N IS THE ORDER OF THE MATRIX. */
184
185/*        T CONTAINS THE NONSYMMETRIC MATRIX.  ITS SUBDIAGONAL IS */
186/*          STORED IN THE LAST N-1 POSITIONS OF THE FIRST COLUMN, */
187/*          ITS DIAGONAL IN THE N POSITIONS OF THE SECOND COLUMN, */
188/*          AND ITS SUPERDIAGONAL IN THE FIRST N-1 POSITIONS OF */
189/*          THE THIRD COLUMN.  T(1,1) AND T(N,3) ARE ARBITRARY. */
190
191/*        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE SYMMETRIC */
192/*          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY. */
193
194/*        M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED. */
195
196/*        Z CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED */
197/*          IN ITS FIRST M COLUMNS. */
198
199/*     ON OUTPUT */
200
201/*        T IS UNALTERED. */
202
203/*        E IS DESTROYED. */
204
205/*        Z CONTAINS THE TRANSFORMED EIGENVECTORS */
206/*          IN ITS FIRST M COLUMNS. */
207
208/*        IERR IS SET TO */
209/*          ZERO       FOR NORMAL RETURN, */
210/*          2*N+I      IF E(I) IS ZERO WITH T(I,1) OR T(I-1,3) NON-ZERO.
211*/
212/*                     IN THIS CASE, THE SYMMETRIC MATRIX IS NOT SIMILAR
213*/
214/*                     TO THE ORIGINAL MATRIX, AND THE EIGENVECTORS */
215/*                     CANNOT BE FOUND BY THIS PROGRAM. */
216
217/*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
218/*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
219*/
220
221/*     THIS VERSION DATED AUGUST 1983. */
222
223/*     ------------------------------------------------------------------
224*/
225
226    /* Parameter adjustments */
227    t_dim1 = *nm;
228    t_offset = t_dim1 + 1;
229    t -= t_offset;
230    --e;
231    z_dim1 = *nm;
232    z_offset = z_dim1 + 1;
233    z -= z_offset;
234
235    /* Function Body */
236    *ierr = 0;
237    if (*m == 0) {
238        goto L1001;
239    }
240    e[1] = 1.;
241    if (*n == 1) {
242        goto L1001;
243    }
244
245    i_1 = *n;
246    for (i = 2; i <= i_1; ++i) {
247        if (e[i] != 0.) {
248            goto L80;
249        }
250        if (t[i + t_dim1] != 0. || t[i - 1 + t_dim1 * 3] != 0.) {
251            goto L1000;
252        }
253        e[i] = 1.;
254        goto L100;
255L80:
256        e[i] = e[i - 1] * e[i] / t[i - 1 + t_dim1 * 3];
257L100:
258        ;
259    }
260
261    i_1 = *m;
262    for (j = 1; j <= i_1; ++j) {
263
264        i_2 = *n;
265        for (i = 2; i <= i_2; ++i) {
266            z[i + j * z_dim1] *= e[i];
267/* L120: */
268        }
269    }
270
271    goto L1001;
272/*     .......... SET ERROR -- EIGENVECTORS CANNOT BE */
273/*                FOUND BY THIS PROGRAM .......... */
274L1000:
275    *ierr = (*n << 1) + i;
276L1001:
277    return 0;
278} /* bakvec_ */
279
280/* Subroutine */ int balanc_(integer *nm, integer *n, doublereal *a, integer *
281        low, integer *igh, doublereal *scale)
282{
283    /* System generated locals */
284    integer a_dim1, a_offset, i_1, i_2;
285    doublereal d_1;
286
287    /* Local variables */
288    static integer iexc;
289    static doublereal c, f, g;
290    static integer i, j, k, l, m;
291    static doublereal r, s, radix, b2;
292    static integer jj;
293    static logical noconv;
294
295
296
297/*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BALANCE, */
298/*     NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH. */
299/*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971). */
300
301/*     THIS SUBROUTINE BALANCES A REAL MATRIX AND ISOLATES */
302/*     EIGENVALUES WHENEVER POSSIBLE. */
303
304/*     ON INPUT */
305
306/*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
307/*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
308/*          DIMENSION STATEMENT. */
309
310/*        N IS THE ORDER OF THE MATRIX. */
311
312/*        A CONTAINS THE INPUT MATRIX TO BE BALANCED. */
313
314/*     ON OUTPUT */
315
316/*        A CONTAINS THE BALANCED MATRIX. */
317
318/*        LOW AND IGH ARE TWO INTEGERS SUCH THAT A(I,J) */
319/*          IS EQUAL TO ZERO IF */
320/*           (1) I IS GREATER THAN J AND */
321/*           (2) J=1,...,LOW-1 OR I=IGH+1,...,N. */
322
323/*        SCALE CONTAINS INFORMATION DETERMINING THE */
324/*           PERMUTATIONS AND SCALING FACTORS USED. */
325
326/*     SUPPOSE THAT THE PRINCIPAL SUBMATRIX IN ROWS LOW THROUGH IGH */
327/*     HAS BEEN BALANCED, THAT P(J) DENOTES THE INDEX INTERCHANGED */
328/*     WITH J DURING THE PERMUTATION STEP, AND THAT THE ELEMENTS */
329/*     OF THE DIAGONAL MATRIX USED ARE DENOTED BY D(I,J).  THEN */
330/*        SCALE(J) = P(J),    FOR J = 1,...,LOW-1 */
331/*                 = D(J,J),      J = LOW,...,IGH */
332/*                 = P(J)         J = IGH+1,...,N. */
333/*     THE ORDER IN WHICH THE INTERCHANGES ARE MADE IS N TO IGH+1, */
334/*     THEN 1 TO LOW-1. */
335
336/*     NOTE THAT 1 IS RETURNED FOR IGH IF IGH IS ZERO FORMALLY. */
337
338/*     THE ALGOL PROCEDURE EXC CONTAINED IN BALANCE APPEARS IN */
339/*     BALANC  IN LINE.  (NOTE THAT THE ALGOL ROLES OF IDENTIFIERS */
340/*     K,L HAVE BEEN REVERSED.) */
341
342/*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
343/*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
344*/
345
346/*     THIS VERSION DATED AUGUST 1983. */
347
348/*     ------------------------------------------------------------------
349*/
350
351    /* Parameter adjustments */
352    --scale;
353    a_dim1 = *nm;
354    a_offset = a_dim1 + 1;
355    a -= a_offset;
356
357    /* Function Body */
358    radix = 16.;
359
360    b2 = radix * radix;
361    k = 1;
362    l = *n;
363    goto L100;
364/*     .......... IN-LINE PROCEDURE FOR ROW AND */
365/*                COLUMN EXCHANGE .......... */
366L20:
367    scale[m] = (doublereal) j;
368    if (j == m) {
369        goto L50;
370    }
371
372    i_1 = l;
373    for (i = 1; i <= i_1; ++i) {
374        f = a[i + j * a_dim1];
375        a[i + j * a_dim1] = a[i + m * a_dim1];
376        a[i + m * a_dim1] = f;
377/* L30: */
378    }
379
380    i_1 = *n;
381    for (i = k; i <= i_1; ++i) {
382        f = a[j + i * a_dim1];
383        a[j + i * a_dim1] = a[m + i * a_dim1];
384        a[m + i * a_dim1] = f;
385/* L40: */
386    }
387
388L50:
389    switch (iexc) {
390        case 1:  goto L80;
391        case 2:  goto L130;
392    }
393/*     .......... SEARCH FOR ROWS ISOLATING AN EIGENVALUE */
394/*                AND PUSH THEM DOWN .......... */
395L80:
396    if (l == 1) {
397        goto L280;
398    }
399    --l;
400/*     .......... FOR J=L STEP -1 UNTIL 1 DO -- .......... */
401L100:
402    i_1 = l;
403    for (jj = 1; jj <= i_1; ++jj) {
404        j = l + 1 - jj;
405
406        i_2 = l;
407        for (i = 1; i <= i_2; ++i) {
408            if (i == j) {
409                goto L110;
410            }
411            if (a[j + i * a_dim1] != 0.) {
412                goto L120;
413            }
414L110:
415            ;
416        }
417
418        m = l;
419        iexc = 1;
420        goto L20;
421L120:
422        ;
423    }
424
425    goto L140;
426/*     .......... SEARCH FOR COLUMNS ISOLATING AN EIGENVALUE */
427/*                AND PUSH THEM LEFT .......... */
428L130:
429    ++k;
430
431L140:
432    i_1 = l;
433    for (j = k; j <= i_1; ++j) {
434
435        i_2 = l;
436        for (i = k; i <= i_2; ++i) {
437            if (i == j) {
438                goto L150;
439            }
440            if (a[i + j * a_dim1] != 0.) {
441                goto L170;
442            }
443L150:
444            ;
445        }
446
447        m = k;
448        iexc = 2;
449        goto L20;
450L170:
451        ;
452    }
453/*     .......... NOW BALANCE THE SUBMATRIX IN ROWS K TO L .......... */
454    i_1 = l;
455    for (i = k; i <= i_1; ++i) {
456/* L180: */
457        scale[i] = 1.;
458    }
459/*     .......... ITERATIVE LOOP FOR NORM REDUCTION .......... */
460L190:
461    noconv = FALSE_;
462
463    i_1 = l;
464    for (i = k; i <= i_1; ++i) {
465        c = 0.;
466        r = 0.;
467
468        i_2 = l;
469        for (j = k; j <= i_2; ++j) {
470            if (j == i) {
471                goto L200;
472            }
473            c += (d_1 = a[j + i * a_dim1], abs(d_1));
474            r += (d_1 = a[i + j * a_dim1], abs(d_1));
475L200:
476            ;
477        }
478/*     .......... GUARD AGAINST ZERO C OR R DUE TO UNDERFLOW .........
479. */
480        if (c == 0. || r == 0.) {
481            goto L270;
482        }
483        g = r / radix;
484        f = 1.;
485        s = c + r;
486L210:
487        if (c >= g) {
488            goto L220;
489        }
490        f *= radix;
491        c *= b2;
492        goto L210;
493L220:
494        g = r * radix;
495L230:
496        if (c < g) {
497            goto L240;
498        }
499        f /= radix;
500        c /= b2;
501        goto L230;
502/*     .......... NOW BALANCE .......... */
503L240:
504        if ((c + r) / f >= s * .95) {
505            goto L270;
506        }
507        g = 1. / f;
508        scale[i] *= f;
509        noconv = TRUE_;
510
511        i_2 = *n;
512        for (j = k; j <= i_2; ++j) {
513/* L250: */
514            a[i + j * a_dim1] *= g;
515        }
516
517        i_2 = l;
518        for (j = 1; j <= i_2; ++j) {
519/* L260: */
520            a[j + i * a_dim1] *= f;
521        }
522
523L270:
524        ;
525    }
526
527    if (noconv) {
528        goto L190;
529    }
530
531L280:
532    *low = k;
533    *igh = l;
534    return 0;
535} /* balanc_ */
536
537/* Subroutine */ int balbak_(integer *nm, integer *n, integer *low, integer *
538        igh, doublereal *scale, integer *m, doublereal *z)
539{
540    /* System generated locals */
541    integer z_dim1, z_offset, i_1, i_2;
542
543    /* Local variables */
544    static integer i, j, k;
545    static doublereal s;
546    static integer ii;
547
548
549
550/*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BALBAK, */
551/*     NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH. */
552/*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971). */
553
554/*     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A REAL GENERAL */
555/*     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING */
556/*     BALANCED MATRIX DETERMINED BY  BALANC. */
557
558/*     ON INPUT */
559
560/*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
561/*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
562/*          DIMENSION STATEMENT. */
563
564/*        N IS THE ORDER OF THE MATRIX. */
565
566/*        LOW AND IGH ARE INTEGERS DETERMINED BY  BALANC. */
567
568/*        SCALE CONTAINS INFORMATION DETERMINING THE PERMUTATIONS */
569/*          AND SCALING FACTORS USED BY  BALANC. */
570
571/*        M IS THE NUMBER OF COLUMNS OF Z TO BE BACK TRANSFORMED. */
572
573/*        Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGEN- */
574/*          VECTORS TO BE BACK TRANSFORMED IN ITS FIRST M COLUMNS. */
575
576/*     ON OUTPUT */
577
578/*        Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE */
579/*          TRANSFORMED EIGENVECTORS IN ITS FIRST M COLUMNS. */
580
581/*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
582/*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
583*/
584
585/*     THIS VERSION DATED AUGUST 1983. */
586
587/*     ------------------------------------------------------------------
588*/
589
590    /* Parameter adjustments */
591    --scale;
592    z_dim1 = *nm;
593    z_offset = z_dim1 + 1;
594    z -= z_offset;
595
596    /* Function Body */
597    if (*m == 0) {
598        goto L200;
599    }
600    if (*igh == *low) {
601        goto L120;
602    }
603
604    i_1 = *igh;
605    for (i = *low; i <= i_1; ++i) {
606        s = scale[i];
607/*     .......... LEFT HAND EIGENVECTORS ARE BACK TRANSFORMED */
608/*                IF THE FOREGOING STATEMENT IS REPLACED BY */
609/*                S=1.0D0/SCALE(I). .......... */
610        i_2 = *m;
611        for (j = 1; j <= i_2; ++j) {
612/* L100: */
613            z[i + j * z_dim1] *= s;
614        }
615
616/* L110: */
617    }
618/*     ......... FOR I=LOW-1 STEP -1 UNTIL 1, */
619/*               IGH+1 STEP 1 UNTIL N DO -- .......... */
620L120:
621    i_1 = *n;
622    for (ii = 1; ii <= i_1; ++ii) {
623        i = ii;
624        if (i >= *low && i <= *igh) {
625            goto L140;
626        }
627        if (i < *low) {
628            i = *low - ii;
629        }
630        k = (integer) scale[i];
631        if (k == i) {
632            goto L140;
633        }
634
635        i_2 = *m;
636        for (j = 1; j <= i_2; ++j) {
637            s = z[i + j * z_dim1];
638            z[i + j * z_dim1] = z[k + j * z_dim1];
639            z[k + j * z_dim1] = s;
640/* L130: */
641        }
642
643L140:
644        ;
645    }
646
647L200:
648    return 0;
649} /* balbak_ */
650
651/* Subroutine */ int bandr_(integer *nm, integer *n, integer *mb, doublereal *
652        a, doublereal *d, doublereal *e, doublereal *e2, logical *matz, 
653        doublereal *z)
654{
655    /* System generated locals */
656    integer a_dim1, a_offset, z_dim1, z_offset, i_1, i_2, i_3, i_4, i_5, 
657            i_6;
658    doublereal d_1;
659
660    /* Builtin functions */
661    double sqrt(doublereal);
662
663    /* Local variables */
664    static doublereal dmin_;
665    static integer maxl, maxr;
666    static doublereal g;
667    static integer j, k, l, r;
668    static doublereal u, b1, b2, c2, f1, f2;
669    static integer i1, i2, j1, j2, m1, n2, r1;
670    static doublereal s2;
671    static integer kr, mr;
672    static doublereal dminrt;
673    static integer ugl;
674
675
676
677/*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BANDRD, */
678/*     NUM. MATH. 12, 231-241(1968) BY SCHWARZ. */
679/*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 273-283(1971). */
680
681/*     THIS SUBROUTINE REDUCES A REAL SYMMETRIC BAND MATRIX */
682/*     TO A SYMMETRIC TRIDIAGONAL MATRIX USING AND OPTIONALLY */
683/*     ACCUMULATING ORTHOGONAL SIMILARITY TRANSFORMATIONS. */
684
685/*     ON INPUT */
686
687/*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
688/*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
689/*          DIMENSION STATEMENT. */
690
691/*        N IS THE ORDER OF THE MATRIX. */
692
693/*        MB IS THE (HALF) BAND WIDTH OF THE MATRIX, DEFINED AS THE */
694/*          NUMBER OF ADJACENT DIAGONALS, INCLUDING THE PRINCIPAL */
695/*          DIAGONAL, REQUIRED TO SPECIFY THE NON-ZERO PORTION OF THE */
696/*          LOWER TRIANGLE OF THE MATRIX. */
697
698/*        A CONTAINS THE LOWER TRIANGLE OF THE SYMMETRIC BAND INPUT */
699/*          MATRIX STORED AS AN N BY MB ARRAY.  ITS LOWEST SUBDIAGONAL */
700/*          IS STORED IN THE LAST N+1-MB POSITIONS OF THE FIRST COLUMN, */
701/*          ITS NEXT SUBDIAGONAL IN THE LAST N+2-MB POSITIONS OF THE */
702/*          SECOND COLUMN, FURTHER SUBDIAGONALS SIMILARLY, AND FINALLY */
703/*          ITS PRINCIPAL DIAGONAL IN THE N POSITIONS OF THE LAST COLUMN.
704*/
705/*          CONTENTS OF STORAGES NOT PART OF THE MATRIX ARE ARBITRARY. */
706
707/*        MATZ SHOULD BE SET TO .TRUE. IF THE TRANSFORMATION MATRIX IS */
708/*          TO BE ACCUMULATED, AND TO .FALSE. OTHERWISE. */
709
710/*     ON OUTPUT */
711
712/*        A HAS BEEN DESTROYED, EXCEPT FOR ITS LAST TWO COLUMNS WHICH */
713/*          CONTAIN A COPY OF THE TRIDIAGONAL MATRIX. */
714
715/*        D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX. */
716
717/*        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL */
718/*          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO. */
719
720/*        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. */
721/*          E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED. */
722
723/*        Z CONTAINS THE ORTHOGONAL TRANSFORMATION MATRIX PRODUCED IN */
724/*          THE REDUCTION IF MATZ HAS BEEN SET TO .TRUE.  OTHERWISE, Z */
725/*          IS NOT REFERENCED. */
726
727/*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
728/*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
729*/
730
731/*     THIS VERSION DATED AUGUST 1983. */
732
733/*     ------------------------------------------------------------------
734*/
735
736    /* Parameter adjustments */
737    z_dim1 = *nm;
738    z_offset = z_dim1 + 1;
739    z -= z_offset;
740    --e2;
741    --e;
742    --d;
743    a_dim1 = *nm;
744    a_offset = a_dim1 + 1;
745    a -= a_offset;
746
747    /* Function Body */
748    dmin_ = 5.4210108624275222e-20;
749    dminrt = 2.3283064365386963e-10;
750/*     .......... INITIALIZE DIAGONAL SCALING MATRIX .......... */
751    i_1 = *n;
752    for (j = 1; j <= i_1; ++j) {
753/* L30: */
754        d[j] = 1.;
755    }
756
757    if (! (*matz)) {
758        goto L60;
759    }
760
761    i_1 = *n;
762    for (j = 1; j <= i_1; ++j) {
763
764        i_2 = *n;
765        for (k = 1; k <= i_2; ++k) {
766/* L40: */
767            z[j + k * z_dim1] = 0.;
768        }
769
770        z[j + j * z_dim1] = 1.;
771/* L50: */
772    }
773
774L60:
775    m1 = *mb - 1;
776    if ((i_1 = m1 - 1) < 0) {
777        goto L900;
778    } else if (i_1 == 0) {
779        goto L800;
780    } else {
781        goto L70;
782    }
783L70:
784    n2 = *n - 2;
785
786    i_1 = n2;
787    for (k = 1; k <= i_1; ++k) {
788/* Computing MIN */
789        i_2 = m1, i_3 = *n - k;
790        maxr = min(i_2,i_3);
791/*     .......... FOR R=MAXR STEP -1 UNTIL 2 DO -- .......... */
792        i_2 = maxr;
793        for (r1 = 2; r1 <= i_2; ++r1) {
794            r = maxr + 2 - r1;
795            kr = k + r;
796            mr = *mb - r;
797            g = a[kr + mr * a_dim1];
798            a[kr - 1 + a_dim1] = a[kr - 1 + (mr + 1) * a_dim1];
799            ugl = k;
800
801            i_3 = *n;
802            i_4 = m1;
803            for (j = kr; i_4 < 0 ? j >= i_3 : j <= i_3; j += i_4) {
804                j1 = j - 1;
805                j2 = j1 - 1;
806                if (g == 0.) {
807                    goto L600;
808                }
809                b1 = a[j1 + a_dim1] / g;
810                b2 = b1 * d[j1] / d[j];
811                s2 = 1. / (b1 * b2 + 1.);
812                if (s2 >= .5) {
813                    goto L450;
814                }
815                b1 = g / a[j1 + a_dim1];
816                b2 = b1 * d[j] / d[j1];
817                c2 = 1. - s2;
818                d[j1] = c2 * d[j1];
819                d[j] = c2 * d[j];
820                f1 = a[j + m1 * a_dim1] * 2.;
821                f2 = b1 * a[j1 + *mb * a_dim1];
822                a[j + m1 * a_dim1] = -b2 * (b1 * a[j + m1 * a_dim1] - a[j + *
823                        mb * a_dim1]) - f2 + a[j + m1 * a_dim1];
824                a[j1 + *mb * a_dim1] = b2 * (b2 * a[j + *mb * a_dim1] + f1) + 
825                        a[j1 + *mb * a_dim1];
826                a[j + *mb * a_dim1] = b1 * (f2 - f1) + a[j + *mb * a_dim1];
827
828                i_5 = j2;
829                for (l = ugl; l <= i_5; ++l) {
830                    i2 = *mb - j + l;
831                    u = a[j1 + (i2 + 1) * a_dim1] + b2 * a[j + i2 * a_dim1];
832                    a[j + i2 * a_dim1] = -b1 * a[j1 + (i2 + 1) * a_dim1] + a[
833                            j + i2 * a_dim1];
834                    a[j1 + (i2 + 1) * a_dim1] = u;
835/* L200: */
836                }
837
838                ugl = j;
839                a[j1 + a_dim1] += b2 * g;
840                if (j == *n) {
841                    goto L350;
842                }
843/* Computing MIN */
844                i_5 = m1, i_6 = *n - j1;
845                maxl = min(i_5,i_6);
846
847                i_5 = maxl;
848                for (l = 2; l <= i_5; ++l) {
849                    i1 = j1 + l;
850                    i2 = *mb - l;
851                    u = a[i1 + i2 * a_dim1] + b2 * a[i1 + (i2 + 1) * a_dim1];
852                    a[i1 + (i2 + 1) * a_dim1] = -b1 * a[i1 + i2 * a_dim1] + a[
853                            i1 + (i2 + 1) * a_dim1];
854                    a[i1 + i2 * a_dim1] = u;
855/* L300: */
856                }
857
858                i1 = j + m1;
859                if (i1 > *n) {
860                    goto L350;
861                }
862                g = b2 * a[i1 + a_dim1];
863L350:
864                if (! (*matz)) {
865                    goto L500;
866                }
867
868                i_5 = *n;
869                for (l = 1; l <= i_5; ++l) {
870                    u = z[l + j1 * z_dim1] + b2 * z[l + j * z_dim1];
871                    z[l + j * z_dim1] = -b1 * z[l + j1 * z_dim1] + z[l + j * 
872                            z_dim1];
873                    z[l + j1 * z_dim1] = u;
874/* L400: */
875                }
876
877                goto L500;
878
879L450:
880                u = d[j1];
881                d[j1] = s2 * d[j];
882                d[j] = s2 * u;
883                f1 = a[j + m1 * a_dim1] * 2.;
884                f2 = b1 * a[j + *mb * a_dim1];
885                u = b1 * (f2 - f1) + a[j1 + *mb * a_dim1];
886                a[j + m1 * a_dim1] = b2 * (b1 * a[j + m1 * a_dim1] - a[j1 + *
887                        mb * a_dim1]) + f2 - a[j + m1 * a_dim1];
888                a[j1 + *mb * a_dim1] = b2 * (b2 * a[j1 + *mb * a_dim1] + f1) 
889                        + a[j + *mb * a_dim1];
890                a[j + *mb * a_dim1] = u;
891
892                i_5 = j2;
893                for (l = ugl; l <= i_5; ++l) {
894                    i2 = *mb - j + l;
895                    u = b2 * a[j1 + (i2 + 1) * a_dim1] + a[j + i2 * a_dim1];
896                    a[j + i2 * a_dim1] = -a[j1 + (i2 + 1) * a_dim1] + b1 * a[
897                            j + i2 * a_dim1];
898                    a[j1 + (i2 + 1) * a_dim1] = u;
899/* L460: */
900                }
901
902                ugl = j;
903                a[j1 + a_dim1] = b2 * a[j1 + a_dim1] + g;
904                if (j == *n) {
905                    goto L480;
906                }
907/* Computing MIN */
908                i_5 = m1, i_6 = *n - j1;
909                maxl = min(i_5,i_6);
910
911                i_5 = maxl;
912                for (l = 2; l <= i_5; ++l) {
913                    i1 = j1 + l;
914                    i2 = *mb - l;
915                    u = b2 * a[i1 + i2 * a_dim1] + a[i1 + (i2 + 1) * a_dim1];
916                    a[i1 + (i2 + 1) * a_dim1] = -a[i1 + i2 * a_dim1] + b1 * a[
917                            i1 + (i2 + 1) * a_dim1];
918                    a[i1 + i2 * a_dim1] = u;
919/* L470: */
920                }
921
922                i1 = j + m1;
923                if (i1 > *n) {
924                    goto L480;
925                }
926                g = a[i1 + a_dim1];
927                a[i1 + a_dim1] = b1 * a[i1 + a_dim1];
928L480:
929                if (! (*matz)) {
930                    goto L500;
931                }
932
933                i_5 = *n;
934                for (l = 1; l <= i_5; ++l) {
935                    u = b2 * z[l + j1 * z_dim1] + z[l + j * z_dim1];
936                    z[l + j * z_dim1] = -z[l + j1 * z_dim1] + b1 * z[l + j * 
937                            z_dim1];
938                    z[l + j1 * z_dim1] = u;
939/* L490: */
940                }
941
942L500:
943                ;
944            }
945
946L600:
947            ;
948        }
949
950        if (k % 64 != 0) {
951            goto L700;
952        }
953/*     .......... RESCALE TO AVOID UNDERFLOW OR OVERFLOW .......... */
954        i_2 = *n;
955        for (j = k; j <= i_2; ++j) {
956            if (d[j] >= dmin_) {
957                goto L650;
958            }
959/* Computing MAX */
960            i_4 = 1, i_3 = *mb + 1 - j;
961            maxl = max(i_4,i_3);
962
963            i_4 = m1;
964            for (l = maxl; l <= i_4; ++l) {
965/* L610: */
966                a[j + l * a_dim1] = dminrt * a[j + l * a_dim1];
967            }
968
969            if (j == *n) {
970                goto L630;
971            }
972/* Computing MIN */
973            i_4 = m1, i_3 = *n - j;
974            maxl = min(i_4,i_3);
975
976            i_4 = maxl;
977            for (l = 1; l <= i_4; ++l) {
978                i1 = j + l;
979                i2 = *mb - l;
980                a[i1 + i2 * a_dim1] = dminrt * a[i1 + i2 * a_dim1];
981/* L620: */
982            }
983
984L630:
985            if (! (*matz)) {
986                goto L645;
987            }
988
989            i_4 = *n;
990            for (l = 1; l <= i_4; ++l) {
991/* L640: */
992                z[l + j * z_dim1] = dminrt * z[l + j * z_dim1];
993            }
994
995L645:
996            a[j + *mb * a_dim1] = dmin_ * a[j + *mb * a_dim1];
997            d[j] /= dmin_;
998L650:
999            ;
1000        }
1001
1002L700:
1003        ;
1004    }
1005/*     .......... FORM SQUARE ROOT OF SCALING MATRIX .......... */
1006L800:
1007    i_1 = *n;
1008    for (j = 2; j <= i_1; ++j) {
1009/* L810: */
1010        e[j] = sqrt(d[j]);
1011    }
1012
1013    if (! (*matz)) {
1014        goto L840;
1015    }
1016
1017    i_1 = *n;
1018    for (j = 1; j <= i_1; ++j) {
1019
1020        i_2 = *n;
1021        for (k = 2; k <= i_2; ++k) {
1022/* L820: */
1023            z[j + k * z_dim1] = e[k] * z[j + k * z_dim1];
1024        }
1025
1026/* L830: */
1027    }
1028
1029L840:
1030    u = 1.;
1031
1032    i_1 = *n;
1033    for (j = 2; j <= i_1; ++j) {
1034        a[j + m1 * a_dim1] = u * e[j] * a[j + m1 * a_dim1];
1035        u = e[j];
1036/* Computing 2nd power */
1037        d_1 = a[j + m1 * a_dim1];
1038        e2[j] = d_1 * d_1;
1039        a[j + *mb * a_dim1] = d[j] * a[j + *mb * a_dim1];
1040        d[j] = a[j + *mb * a_dim1];
1041        e[j] = a[j + m1 * a_dim1];
1042/* L850: */
1043    }
1044
1045    d[1] = a[*mb * a_dim1 + 1];
1046    e[1] = 0.;
1047    e2[1] = 0.;
1048    goto L1001;
1049
1050L900:
1051    i_1 = *n;
1052    for (j = 1; j <= i_1; ++j) {
1053        d[j] = a[j + *mb * a_dim1];
1054        e[j] = 0.;
1055        e2[j] = 0.;
1056/* L950: */
1057    }
1058
1059L1001:
1060    return 0;
1061} /* bandr_ */
1062
1063/* Subroutine */ int bandv_(integer *nm, integer *n, integer *mbw, doublereal
1064        *a, doublereal *e21, integer *m, doublereal *w, doublereal *z, 
1065        integer *ierr, integer */*nv*/, doublereal *rv, doublereal *rv6)
1066{
1067    /* System generated locals */
1068    integer a_dim1, a_offset, z_dim1, z_offset, i_1, i_2, i_3, i_4, i_5;
1069    doublereal d_1;
1070
1071    /* Builtin functions */
1072    double sqrt(doublereal), d_sign(doublereal *, doublereal *);
1073
1074    /* Local variables */
1075    static integer maxj, maxk;
1076    static doublereal norm;
1077    static integer i, j, k, r;
1078    static doublereal u, v, order;
1079    static integer group, m1;
1080    static doublereal x0, x1;
1081    static integer mb, m21, ii, ij, jj, kj;
1082    static doublereal uk, xu;
1083    extern doublereal pythag_(doublereal *, doublereal *), epslon_(doublereal
1084            *);
1085    static integer ij1, kj1, its;
1086    static doublereal eps2, eps3, eps4;
1087
1088
1089
1090/*     THIS SUBROUTINE FINDS THOSE EIGENVECTORS OF A REAL SYMMETRIC */
1091/*     BAND MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES, USING INVERSE
1092*/
1093/*     ITERATION.  THE SUBROUTINE MAY ALSO BE USED TO SOLVE SYSTEMS */
1094/*     OF LINEAR EQUATIONS WITH A SYMMETRIC OR NON-SYMMETRIC BAND */
1095/*     COEFFICIENT MATRIX. */
1096
1097/*     ON INPUT */
1098
1099/*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
1100/*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
1101/*          DIMENSION STATEMENT. */
1102
1103/*        N IS THE ORDER OF THE MATRIX. */
1104
1105/*        MBW IS THE NUMBER OF COLUMNS OF THE ARRAY A USED TO STORE THE */
1106/*          BAND MATRIX.  IF THE MATRIX IS SYMMETRIC, MBW IS ITS (HALF) */
1107/*          BAND WIDTH, DENOTED MB AND DEFINED AS THE NUMBER OF ADJACENT
1108*/
1109/*          DIAGONALS, INCLUDING THE PRINCIPAL DIAGONAL, REQUIRED TO */
1110/*          SPECIFY THE NON-ZERO PORTION OF THE LOWER TRIANGLE OF THE */
1111/*          MATRIX.  IF THE SUBROUTINE IS BEING USED TO SOLVE SYSTEMS */
1112/*          OF LINEAR EQUATIONS AND THE COEFFICIENT MATRIX IS NOT */
1113/*          SYMMETRIC, IT MUST HOWEVER HAVE THE SAME NUMBER OF ADJACENT */
1114/*          DIAGONALS ABOVE THE MAIN DIAGONAL AS BELOW, AND IN THIS */
1115/*          CASE, MBW=2*MB-1. */
1116
1117/*        A CONTAINS THE LOWER TRIANGLE OF THE SYMMETRIC BAND INPUT */
1118/*          MATRIX STORED AS AN N BY MB ARRAY.  ITS LOWEST SUBDIAGONAL */
1119/*          IS STORED IN THE LAST N+1-MB POSITIONS OF THE FIRST COLUMN, */
1120/*          ITS NEXT SUBDIAGONAL IN THE LAST N+2-MB POSITIONS OF THE */
1121/*          SECOND COLUMN, FURTHER SUBDIAGONALS SIMILARLY, AND FINALLY */
1122/*          ITS PRINCIPAL DIAGONAL IN THE N POSITIONS OF COLUMN MB. */
1123/*          IF THE SUBROUTINE IS BEING USED TO SOLVE SYSTEMS OF LINEAR */
1124/*          EQUATIONS AND THE COEFFICIENT MATRIX IS NOT SYMMETRIC, A IS */
1125/*          N BY 2*MB-1 INSTEAD WITH LOWER TRIANGLE AS ABOVE AND WITH */
1126/*          ITS FIRST SUPERDIAGONAL STORED IN THE FIRST N-1 POSITIONS OF
1127*/
1128/*          COLUMN MB+1, ITS SECOND SUPERDIAGONAL IN THE FIRST N-2 */
1129/*          POSITIONS OF COLUMN MB+2, FURTHER SUPERDIAGONALS SIMILARLY, */
1130/*          AND FINALLY ITS HIGHEST SUPERDIAGONAL IN THE FIRST N+1-MB */
1131/*          POSITIONS OF THE LAST COLUMN. */
1132/*          CONTENTS OF STORAGES NOT PART OF THE MATRIX ARE ARBITRARY. */
1133
1134/*        E21 SPECIFIES THE ORDERING OF THE EIGENVALUES AND CONTAINS */
1135/*            0.0D0 IF THE EIGENVALUES ARE IN ASCENDING ORDER, OR */
1136/*            2.0D0 IF THE EIGENVALUES ARE IN DESCENDING ORDER. */
1137/*          IF THE SUBROUTINE IS BEING USED TO SOLVE SYSTEMS OF LINEAR */
1138/*          EQUATIONS, E21 SHOULD BE SET TO 1.0D0 IF THE COEFFICIENT */
1139/*          MATRIX IS SYMMETRIC AND TO -1.0D0 IF NOT. */
1140
1141/*        M IS THE NUMBER OF SPECIFIED EIGENVALUES OR THE NUMBER OF */
1142/*          SYSTEMS OF LINEAR EQUATIONS. */
1143
1144/*        W CONTAINS THE M EIGENVALUES IN ASCENDING OR DESCENDING ORDER.
1145*/
1146/*          IF THE SUBROUTINE IS BEING USED TO SOLVE SYSTEMS OF LINEAR */
1147/*          EQUATIONS (A-W(R)*I)*X(R)=B(R), WHERE I IS THE IDENTITY */
1148/*          MATRIX, W(R) SHOULD BE SET ACCORDINGLY, FOR R=1,2,...,M. */
1149
1150/*        Z CONTAINS THE CONSTANT MATRIX COLUMNS (B(R),R=1,2,...,M), IF */
1151/*          THE SUBROUTINE IS USED TO SOLVE SYSTEMS OF LINEAR EQUATIONS.
1152*/
1153
1154/*        NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER RV */
1155/*          AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT. */
1156
1157/*     ON OUTPUT */
1158
1159/*        A AND W ARE UNALTERED. */
1160
1161/*        Z CONTAINS THE ASSOCIATED SET OF ORTHOGONAL EIGENVECTORS. */
1162/*          ANY VECTOR WHICH FAILS TO CONVERGE IS SET TO ZERO.  IF THE */
1163/*          SUBROUTINE IS USED TO SOLVE SYSTEMS OF LINEAR EQUATIONS, */
1164/*          Z CONTAINS THE SOLUTION MATRIX COLUMNS (X(R),R=1,2,...,M). */
1165
1166/*        IERR IS SET TO */
1167/*          ZERO       FOR NORMAL RETURN, */
1168/*          -R         IF THE EIGENVECTOR CORRESPONDING TO THE R-TH */
1169/*                     EIGENVALUE FAILS TO CONVERGE, OR IF THE R-TH */
1170/*                     SYSTEM OF LINEAR EQUATIONS IS NEARLY SINGULAR. */
1171
1172/*        RV AND RV6 ARE TEMPORARY STORAGE ARRAYS.  NOTE THAT RV IS */
1173/*          OF DIMENSION AT LEAST N*(2*MB-1).  IF THE SUBROUTINE */
1174/*          IS BEING USED TO SOLVE SYSTEMS OF LINEAR EQUATIONS, THE */
1175/*          DETERMINANT (UP TO SIGN) OF A-W(M)*I IS AVAILABLE, UPON */
1176/*          RETURN, AS THE PRODUCT OF THE FIRST N ELEMENTS OF RV. */
1177
1178/*     CALLS PYTHAG FOR  DSQRT(A*A + B*B) . */
1179
1180/*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
1181/*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
1182*/
1183
1184/*     THIS VERSION DATED AUGUST 1983. */
1185
1186/*     ------------------------------------------------------------------
1187*/
1188
1189    /* Parameter adjustments */
1190    --rv6;
1191    a_dim1 = *nm;
1192    a_offset = a_dim1 + 1;
1193    a -= a_offset;
1194    z_dim1 = *nm;
1195    z_offset = z_dim1 + 1;
1196    z -= z_offset;
1197    --w;
1198    --rv;
1199
1200    /* Function Body */
1201    *ierr = 0;
1202    if (*m == 0) {
1203        goto L1001;
1204    }
1205    mb = *mbw;
1206    if (*e21 < 0.) {
1207        mb = (*mbw + 1) / 2;
1208    }
1209    m1 = mb - 1;
1210    m21 = m1 + mb;
1211    order = 1. - abs(*e21);
1212/*     .......... FIND VECTORS BY INVERSE ITERATION .......... */
1213    i_1 = *m;
1214    for (r = 1; r <= i_1; ++r) {
1215        its = 1;
1216        x1 = w[r];
1217        if (r != 1) {
1218            goto L100;
1219        }
1220/*     .......... COMPUTE NORM OF MATRIX .......... */
1221        norm = 0.;
1222
1223        i_2 = mb;
1224        for (j = 1; j <= i_2; ++j) {
1225            jj = mb + 1 - j;
1226            kj = jj + m1;
1227            ij = 1;
1228            v = 0.;
1229
1230            i_3 = *n;
1231            for (i = jj; i <= i_3; ++i) {
1232                v += (d_1 = a[i + j * a_dim1], abs(d_1));
1233                if (*e21 >= 0.) {
1234                    goto L40;
1235                }
1236                v += (d_1 = a[ij + kj * a_dim1], abs(d_1));
1237                ++ij;
1238L40:
1239                ;
1240            }
1241
1242            norm = max(norm,v);
1243/* L60: */
1244        }
1245
1246        if (*e21 < 0.) {
1247            norm *= .5;
1248        }
1249/*     .......... EPS2 IS THE CRITERION FOR GROUPING, */
1250/*                EPS3 REPLACES ZERO PIVOTS AND EQUAL */
1251/*                ROOTS ARE MODIFIED BY EPS3, */
1252/*                EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW .........
1253. */
1254        if (norm == 0.) {
1255            norm = 1.;
1256        }
1257        eps2 = norm * .001 * abs(order);
1258        eps3 = epslon_(&norm);
1259        uk = (doublereal) (*n);
1260        uk = sqrt(uk);
1261        eps4 = uk * eps3;
1262L80:
1263        group = 0;
1264        goto L120;
1265/*     .......... LOOK FOR CLOSE OR COINCIDENT ROOTS .......... */
1266L100:
1267        if ((d_1 = x1 - x0, abs(d_1)) >= eps2) {
1268            goto L80;
1269        }
1270        ++group;
1271        if (order * (x1 - x0) <= 0.) {
1272            x1 = x0 + order * eps3;
1273        }
1274/*     .......... EXPAND MATRIX, SUBTRACT EIGENVALUE, */
1275/*                AND INITIALIZE VECTOR .......... */
1276L120:
1277        i_2 = *n;
1278        for (i = 1; i <= i_2; ++i) {
1279/* Computing MIN */
1280            i_3 = 0, i_4 = i - m1;
1281            ij = i + min(i_3,i_4) * *n;
1282            kj = ij + mb * *n;
1283            ij1 = kj + m1 * *n;
1284            if (m1 == 0) {
1285                goto L180;
1286            }
1287
1288            i_3 = m1;
1289            for (j = 1; j <= i_3; ++j) {
1290                if (ij > m1) {
1291                    goto L125;
1292                }
1293                if (ij > 0) {
1294                    goto L130;
1295                }
1296                rv[ij1] = 0.;
1297                ij1 += *n;
1298                goto L130;
1299L125:
1300                rv[ij] = a[i + j * a_dim1];
1301L130:
1302                ij += *n;
1303                ii = i + j;
1304                if (ii > *n) {
1305                    goto L150;
1306                }
1307                jj = mb - j;
1308                if (*e21 >= 0.) {
1309                    goto L140;
1310                }
1311                ii = i;
1312                jj = mb + j;
1313L140:
1314                rv[kj] = a[ii + jj * a_dim1];
1315                kj += *n;
1316L150:
1317                ;
1318            }
1319
1320L180:
1321            rv[ij] = a[i + mb * a_dim1] - x1;
1322            rv6[i] = eps4;
1323            if (order == 0.) {
1324                rv6[i] = z[i + r * z_dim1];
1325            }
1326/* L200: */
1327        }
1328
1329        if (m1 == 0) {
1330            goto L600;
1331        }
1332/*     .......... ELIMINATION WITH INTERCHANGES .......... */
1333        i_2 = *n;
1334        for (i = 1; i <= i_2; ++i) {
1335            ii = i + 1;
1336/* Computing MIN */
1337            i_3 = i + m1 - 1;
1338            maxk = min(i_3,*n);
1339/* Computing MIN */
1340            i_3 = *n - i, i_4 = m21 - 2;
1341            maxj = min(i_3,i_4) * *n;
1342
1343            i_3 = maxk;
1344            for (k = i; k <= i_3; ++k) {
1345                kj1 = k;
1346                j = kj1 + *n;
1347                jj = j + maxj;
1348
1349                i_4 = jj;
1350                i_5 = *n;
1351                for (kj = j; i_5 < 0 ? kj >= i_4 : kj <= i_4; kj += i_5) {
1352                    rv[kj1] = rv[kj];
1353                    kj1 = kj;
1354/* L340: */
1355                }
1356
1357                rv[kj1] = 0.;
1358/* L360: */
1359            }
1360
1361            if (i == *n) {
1362                goto L580;
1363            }
1364            u = 0.;
1365/* Computing MIN */
1366            i_3 = i + m1;
1367            maxk = min(i_3,*n);
1368/* Computing MIN */
1369            i_3 = *n - ii, i_5 = m21 - 2;
1370            maxj = min(i_3,i_5) * *n;
1371
1372            i_3 = maxk;
1373            for (j = i; j <= i_3; ++j) {
1374                if ((d_1 = rv[j], abs(d_1)) < abs(u)) {
1375                    goto L450;
1376                }
1377                u = rv[j];
1378                k = j;
1379L450:
1380                ;
1381            }
1382
1383            j = i + *n;
1384            jj = j + maxj;
1385            if (k == i) {
1386                goto L520;
1387            }
1388            kj = k;
1389
1390            i_3 = jj;
1391            i_5 = *n;
1392            for (ij = i; i_5 < 0 ? ij >= i_3 : ij <= i_3; ij += i_5) {
1393                v = rv[ij];
1394                rv[ij] = rv[kj];
1395                rv[kj] = v;
1396                kj += *n;
1397/* L500: */
1398            }
1399
1400            if (order != 0.) {
1401                goto L520;
1402            }
1403            v = rv6[i];
1404            rv6[i] = rv6[k];
1405            rv6[k] = v;
1406L520:
1407            if (u == 0.) {
1408                goto L580;
1409            }
1410
1411            i_5 = maxk;
1412            for (k = ii; k <= i_5; ++k) {
1413                v = rv[k] / u;
1414                kj = k;
1415
1416                i_3 = jj;
1417                i_4 = *n;
1418                for (ij = j; i_4 < 0 ? ij >= i_3 : ij <= i_3; ij += i_4) {
1419                    kj += *n;
1420                    rv[kj] -= v * rv[ij];
1421/* L540: */
1422                }
1423
1424                if (order == 0.) {
1425                    rv6[k] -= v * rv6[i];
1426                }
1427/* L560: */
1428            }
1429
1430L580:
1431            ;
1432        }
1433/*     .......... BACK SUBSTITUTION */
1434/*                FOR I=N STEP -1 UNTIL 1 DO -- .......... */
1435L600:
1436        i_2 = *n;
1437        for (ii = 1; ii <= i_2; ++ii) {
1438            i = *n + 1 - ii;
1439            maxj = min(ii,m21);
1440            if (maxj == 1) {
1441                goto L620;
1442            }
1443            ij1 = i;
1444            j = ij1 + *n;
1445            jj = j + (maxj - 2) * *n;
1446
1447            i_5 = jj;
1448            i_4 = *n;
1449            for (ij = j; i_4 < 0 ? ij >= i_5 : ij <= i_5; ij += i_4) {
1450                ++ij1;
1451                rv6[i] -= rv[ij] * rv6[ij1];
1452/* L610: */
1453            }
1454
1455L620:
1456            v = rv[i];
1457            if (abs(v) >= eps3) {
1458                goto L625;
1459            }
1460/*     .......... SET ERROR -- NEARLY SINGULAR LINEAR SYSTEM .....
1461..... */
1462            if (order == 0.) {
1463                *ierr = -r;
1464            }
1465            v = d_sign(&eps3, &v);
1466L625:
1467            rv6[i] /= v;
1468/* L630: */
1469        }
1470
1471        xu = 1.;
1472        if (order == 0.) {
1473            goto L870;
1474        }
1475/*     .......... ORTHOGONALIZE WITH RESPECT TO PREVIOUS */
1476/*                MEMBERS OF GROUP .......... */
1477        if (group == 0) {
1478            goto L700;
1479        }
1480
1481        i_2 = group;
1482        for (jj = 1; jj <= i_2; ++jj) {
1483            j = r - group - 1 + jj;
1484            xu = 0.;
1485
1486            i_4 = *n;
1487            for (i = 1; i <= i_4; ++i) {
1488/* L640: */
1489                xu += rv6[i] * z[i + j * z_dim1];
1490            }
1491
1492            i_4 = *n;
1493            for (i = 1; i <= i_4; ++i) {
1494/* L660: */
1495                rv6[i] -= xu * z[i + j * z_dim1];
1496            }
1497
1498/* L680: */
1499        }
1500
1501L700:
1502        norm = 0.;
1503
1504        i_2 = *n;
1505        for (i = 1; i <= i_2; ++i) {
1506/* L720: */
1507            norm += (d_1 = rv6[i], abs(d_1));
1508        }
1509
1510        if (norm >= .1) {
1511            goto L840;
1512        }
1513/*     .......... IN-LINE PROCEDURE FOR CHOOSING */
1514/*                A NEW STARTING VECTOR .......... */
1515        if (its >= *n) {
1516            goto L830;
1517        }
1518        ++its;
1519        xu = eps4 / (uk + 1.);
1520        rv6[1] = eps4;
1521
1522        i_2 = *n;
1523        for (i = 2; i <= i_2; ++i) {
1524/* L760: */
1525            rv6[i] = xu;
1526        }
1527
1528        rv6[its] -= eps4 * uk;
1529        goto L600;
1530/*     .......... SET ERROR -- NON-CONVERGED EIGENVECTOR .......... */
1531L830:
1532        *ierr = -r;
1533        xu = 0.;
1534        goto L870;
1535/*     .......... NORMALIZE SO THAT SUM OF SQUARES IS */
1536/*                1 AND EXPAND TO FULL ORDER .......... */
1537L840:
1538        u = 0.;
1539
1540        i_2 = *n;
1541        for (i = 1; i <= i_2; ++i) {
1542/* L860: */
1543            u = pythag_(&u, &rv6[i]);
1544        }
1545
1546        xu = 1. / u;
1547
1548L870:
1549        i_2 = *n;
1550        for (i = 1; i <= i_2; ++i) {
1551/* L900: */
1552            z[i + r * z_dim1] = rv6[i] * xu;
1553        }
1554
1555        x0 = x1;
1556/* L920: */
1557    }
1558
1559L1001:
1560    return 0;
1561} /* bandv_ */
1562
1563/* Subroutine */ int bisect_(integer *n, doublereal *eps1, doublereal *d, 
1564        doublereal *e, doublereal *e2, doublereal *lb, doublereal *ub, 
1565        integer *mm, integer *m, doublereal *w, integer *ind, integer *ierr, 
1566        doublereal *rv4, doublereal *rv5)
1567{
1568    /* System generated locals */
1569    integer i_1, i_2;
1570    doublereal d_1, d_2, d_3;
1571
1572    /* Local variables */
1573    static integer i, j, k, l, p, q, r, s;
1574    static doublereal u, v;
1575    static integer m1, m2;
1576    static doublereal t1, t2, x0, x1;
1577    static integer ii;
1578    static doublereal xu;
1579    extern doublereal epslon_(doublereal *);
1580    static integer isturm, tag;
1581    static doublereal tst1, tst2;
1582
1583
1584
1585/*     THIS SUBROUTINE IS A TRANSLATION OF THE BISECTION TECHNIQUE */
1586/*     IN THE ALGOL PROCEDURE TRISTURM BY PETERS AND WILKINSON. */
1587/*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971). */
1588
1589/*     THIS SUBROUTINE FINDS THOSE EIGENVALUES OF A TRIDIAGONAL */
1590/*     SYMMETRIC MATRIX WHICH LIE IN A SPECIFIED INTERVAL, */
1591/*     USING BISECTION. */
1592
1593/*     ON INPUT */
1594
1595/*        N IS THE ORDER OF THE MATRIX. */
1596
1597/*        EPS1 IS AN ABSOLUTE ERROR TOLERANCE FOR THE COMPUTED */
1598/*          EIGENVALUES.  IF THE INPUT EPS1 IS NON-POSITIVE, */
1599/*          IT IS RESET FOR EACH SUBMATRIX TO A DEFAULT VALUE, */
1600/*          NAMELY, MINUS THE PRODUCT OF THE RELATIVE MACHINE */
1601/*          PRECISION AND THE 1-NORM OF THE SUBMATRIX. */
1602
1603/*        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. */
1604
1605/*        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX */
1606/*          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY. */
1607
1608/*        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. */
1609/*          E2(1) IS ARBITRARY. */
1610
1611/*        LB AND UB DEFINE THE INTERVAL TO BE SEARCHED FOR EIGENVALUES. */
1612/*          IF LB IS NOT LESS THAN UB, NO EIGENVALUES WILL BE FOUND. */
1613
1614/*        MM SHOULD BE SET TO AN UPPER BOUND FOR THE NUMBER OF */
1615/*          EIGENVALUES IN THE INTERVAL.  WARNING. IF MORE THAN */
1616/*          MM EIGENVALUES ARE DETERMINED TO LIE IN THE INTERVAL, */
1617/*          AN ERROR RETURN IS MADE WITH NO EIGENVALUES FOUND. */
1618
1619/*     ON OUTPUT */
1620
1621/*        EPS1 IS UNALTERED UNLESS IT HAS BEEN RESET TO ITS */
1622/*          (LAST) DEFAULT VALUE. */
1623
1624/*        D AND E ARE UNALTERED. */
1625
1626/*        ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED */
1627/*          AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE */
1628/*          MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES. */
1629/*          E2(1) IS ALSO SET TO ZERO. */
1630
1631/*        M IS THE NUMBER OF EIGENVALUES DETERMINED TO LIE IN (LB,UB). */
1632
1633/*        W CONTAINS THE M EIGENVALUES IN ASCENDING ORDER. */
1634
1635/*        IND CONTAINS IN ITS FIRST M POSITIONS THE SUBMATRIX INDICES */
1636/*          ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W -- */
1637/*          1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX FROM */
1638/*          THE TOP, 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC..
1639*/
1640
1641/*        IERR IS SET TO */
1642/*          ZERO       FOR NORMAL RETURN, */
1643/*          3*N+1      IF M EXCEEDS MM. */
1644
1645/*        RV4 AND RV5 ARE TEMPORARY STORAGE ARRAYS. */
1646
1647/*     THE ALGOL PROCEDURE STURMCNT CONTAINED IN TRISTURM */
1648/*     APPEARS IN BISECT IN-LINE. */
1649
1650/*     NOTE THAT SUBROUTINE TQL1 OR IMTQL1 IS GENERALLY FASTER THAN */
1651/*     BISECT, IF MORE THAN N/4 EIGENVALUES ARE TO BE FOUND. */
1652
1653/*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
1654/*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
1655*/
1656
1657/*     THIS VERSION DATED AUGUST 1983. */
1658
1659/*     ------------------------------------------------------------------
1660*/
1661
1662    /* Parameter adjustments */
1663    --rv5;
1664    --rv4;
1665    --e2;
1666    --e;
1667    --d;
1668    --ind;
1669    --w;
1670
1671    /* Function Body */
1672    *ierr = 0;
1673    tag = 0;
1674    t1 = *lb;
1675    t2 = *ub;
1676/*     .......... LOOK FOR SMALL SUB-DIAGONAL ENTRIES .......... */
1677    i_1 = *n;
1678    for (i = 1; i <= i_1; ++i) {
1679        if (i == 1) {
1680            goto L20;
1681        }
1682        tst1 = (d_1 = d[i], abs(d_1)) + (d_2 = d[i - 1], abs(d_2));
1683        tst2 = tst1 + (d_1 = e[i], abs(d_1));
1684        if (tst2 > tst1) {
1685            goto L40;
1686        }
1687L20:
1688        e2[i] = 0.;
1689L40:
1690        ;
1691    }
1692/*     .......... DETERMINE THE NUMBER OF EIGENVALUES */
1693/*                IN THE INTERVAL .......... */
1694    p = 1;
1695    q = *n;
1696    x1 = *ub;
1697    isturm = 1;
1698    goto L320;
1699L60:
1700    *m = s;
1701    x1 = *lb;
1702    isturm = 2;
1703    goto L320;
1704L80:
1705    *m -= s;
1706    if (*m > *mm) {
1707        goto L980;
1708    }
1709    q = 0;
1710    r = 0;
1711/*     .......... ESTABLISH AND PROCESS NEXT SUBMATRIX, REFINING */
1712/*                INTERVAL BY THE GERSCHGORIN BOUNDS .......... */
1713L100:
1714    if (r == *m) {
1715        goto L1001;
1716    }
1717    ++tag;
1718    p = q + 1;
1719    xu = d[p];
1720    x0 = d[p];
1721    u = 0.;
1722
1723    i_1 = *n;
1724    for (q = p; q <= i_1; ++q) {
1725        x1 = u;
1726        u = 0.;
1727        v = 0.;
1728        if (q == *n) {
1729            goto L110;
1730        }
1731        u = (d_1 = e[q + 1], abs(d_1));
1732        v = e2[q + 1];
1733L110:
1734/* Computing MIN */
1735        d_1 = d[q] - (x1 + u);
1736        xu = min(d_1,xu);
1737/* Computing MAX */
1738        d_1 = d[q] + (x1 + u);
1739        x0 = max(d_1,x0);
1740        if (v == 0.) {
1741            goto L140;
1742        }
1743/* L120: */
1744    }
1745
1746L140:
1747/* Computing MAX */
1748    d_2 = abs(xu), d_3 = abs(x0);
1749    d_1 = max(d_2,d_3);
1750    x1 = epslon_(&d_1);
1751    if (*eps1 <= 0.) {
1752        *eps1 = -x1;
1753    }
1754    if (p != q) {
1755        goto L180;
1756    }
1757/*     .......... CHECK FOR ISOLATED ROOT WITHIN INTERVAL .......... */
1758    if (t1 > d[p] || d[p] >= t2) {
1759        goto L940;
1760    }
1761    m1 = p;
1762    m2 = p;
1763    rv5[p] = d[p];
1764    goto L900;
1765L180:
1766    x1 *= q - p + 1;
1767/* Computing MAX */
1768    d_1 = t1, d_2 = xu - x1;
1769    *lb = max(d_1,d_2);
1770/* Computing MIN */
1771    d_1 = t2, d_2 = x0 + x1;
1772    *ub = min(d_1,d_2);
1773    x1 = *lb;
1774    isturm = 3;
1775    goto L320;
1776L200:
1777    m1 = s + 1;
1778    x1 = *ub;
1779    isturm = 4;
1780    goto L320;
1781L220:
1782    m2 = s;
1783    if (m1 > m2) {
1784        goto L940;
1785    }
1786/*     .......... FIND ROOTS BY BISECTION .......... */
1787    x0 = *ub;
1788    isturm = 5;
1789
1790    i_1 = m2;
1791    for (i = m1; i <= i_1; ++i) {
1792        rv5[i] = *ub;
1793        rv4[i] = *lb;
1794/* L240: */
1795    }
1796/*     .......... LOOP FOR K-TH EIGENVALUE */
1797/*                FOR K=M2 STEP -1 UNTIL M1 DO -- */
1798/*                (-DO- NOT USED TO LEGALIZE -COMPUTED GO TO-) ..........
1799*/
1800    k = m2;
1801L250:
1802    xu = *lb;
1803/*     .......... FOR I=K STEP -1 UNTIL M1 DO -- .......... */
1804    i_1 = k;
1805    for (ii = m1; ii <= i_1; ++ii) {
1806        i = m1 + k - ii;
1807        if (xu >= rv4[i]) {
1808            goto L260;
1809        }
1810        xu = rv4[i];
1811        goto L280;
1812L260:
1813        ;
1814    }
1815
1816L280:
1817    if (x0 > rv5[k]) {
1818        x0 = rv5[k];
1819    }
1820/*     .......... NEXT BISECTION STEP .......... */
1821L300:
1822    x1 = (xu + x0) * .5;
1823    if (x0 - xu <= abs(*eps1)) {
1824        goto L420;
1825    }
1826    tst1 = (abs(xu) + abs(x0)) * 2.;
1827    tst2 = tst1 + (x0 - xu);
1828    if (tst2 == tst1) {
1829        goto L420;
1830    }
1831/*     .......... IN-LINE PROCEDURE FOR STURM SEQUENCE .......... */
1832L320:
1833    s = p - 1;
1834    u = 1.;
1835
1836    i_1 = q;
1837    for (i = p; i <= i_1; ++i) {
1838        if (u != 0.) {
1839            goto L325;
1840        }
1841        v = (d_1 = e[i], abs(d_1)) / epslon_(&c_b141);
1842        if (e2[i] == 0.) {
1843            v = 0.;
1844        }
1845        goto L330;
1846L325:
1847        v = e2[i] / u;
1848L330:
1849        u = d[i] - x1 - v;
1850        if (u < 0.) {
1851            ++s;
1852        }
1853/* L340: */
1854    }
1855
1856    switch (isturm) {
1857        case 1:  goto L60;
1858        case 2:  goto L80;
1859        case 3:  goto L200;
1860        case 4:  goto L220;
1861        case 5:  goto L360;
1862    }
1863/*     .......... REFINE INTERVALS .......... */
1864L360:
1865    if (s >= k) {
1866        goto L400;
1867    }
1868    xu = x1;
1869    if (s >= m1) {
1870        goto L380;
1871    }
1872    rv4[m1] = x1;
1873    goto L300;
1874L380:
1875    rv4[s + 1] = x1;
1876    if (rv5[s] > x1) {
1877        rv5[s] = x1;
1878    }
1879    goto L300;
1880L400:
1881    x0 = x1;
1882    goto L300;
1883/*     .......... K-TH EIGENVALUE FOUND .......... */
1884L420:
1885    rv5[k] = x1;
1886    --k;
1887    if (k >= m1) {
1888        goto L250;
1889    }
1890/*     .......... ORDER EIGENVALUES TAGGED WITH THEIR */
1891/*                SUBMATRIX ASSOCIATIONS .......... */
1892L900:
1893    s = r;
1894    r = r + m2 - m1 + 1;
1895    j = 1;
1896    k = m1;
1897
1898    i_1 = r;
1899    for (l = 1; l <= i_1; ++l) {
1900        if (j > s) {
1901            goto L910;
1902        }
1903        if (k > m2) {
1904            goto L940;
1905        }
1906        if (rv5[k] >= w[l]) {
1907            goto L915;
1908        }
1909
1910        i_2 = s;
1911        for (ii = j; ii <= i_2; ++ii) {
1912            i = l + s - ii;
1913            w[i + 1] = w[i];
1914            ind[i + 1] = ind[i];
1915/* L905: */
1916        }
1917
1918L910:
1919        w[l] = rv5[k];
1920        ind[l] = tag;
1921        ++k;
1922        goto L920;
1923L915:
1924        ++j;
1925L920:
1926        ;
1927    }
1928
1929L940:
1930    if (q < *n) {
1931        goto L100;
1932    }
1933    goto L1001;
1934/*     .......... SET ERROR -- UNDERESTIMATE OF NUMBER OF */
1935/*                EIGENVALUES IN INTERVAL .......... */
1936L980:
1937    *ierr = *n * 3 + 1;
1938L1001:
1939    *lb = t1;
1940    *ub = t2;
1941    return 0;
1942} /* bisect_ */
1943
1944/* Subroutine */ int bqr_(integer *nm, integer *n, integer *mb, doublereal *a,
1945         doublereal *t, doublereal *r, integer *ierr, integer */*nv*/, doublereal
1946        *rv)
1947{
1948    /* System generated locals */
1949    integer a_dim1, a_offset, i_1, i_2, i_3;
1950    doublereal d_1;
1951
1952    /* Builtin functions */
1953    double d_sign(doublereal *, doublereal *), sqrt(doublereal);
1954
1955    /* Local variables */
1956    static doublereal f, g;
1957    static integer i, j, k, l, m;
1958    static doublereal q, s, scale;
1959    static integer imult, m1, m2, m3, m4, m21, m31, ii, ik, jk, kj, jm, kk, 
1960            km, ll, mk, mn, ni, mz;
1961    extern doublereal pythag_(doublereal *, doublereal *);
1962    static integer kj1, its;
1963    static doublereal tst1, tst2;
1964
1965
1966
1967/*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BQR, */
1968/*     NUM. MATH. 16, 85-92(1970) BY MARTIN, REINSCH, AND WILKINSON. */
1969/*     HANDBOOK FOR AUTO. COMP., VOL II-LINEAR ALGEBRA, 266-272(1971). */
1970
1971/*     THIS SUBROUTINE FINDS THE EIGENVALUE OF SMALLEST (USUALLY) */
1972/*     MAGNITUDE OF A REAL SYMMETRIC BAND MATRIX USING THE */
1973/*     QR ALGORITHM WITH SHIFTS OF ORIGIN.  CONSECUTIVE CALLS */
1974/*     CAN BE MADE TO FIND FURTHER EIGENVALUES. */
1975
1976/*     ON INPUT */
1977
1978/*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
1979/*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
1980/*          DIMENSION STATEMENT. */
1981
1982/*        N IS THE ORDER OF THE MATRIX. */
1983
1984/*        MB IS THE (HALF) BAND WIDTH OF THE MATRIX, DEFINED AS THE */
1985/*          NUMBER OF ADJACENT DIAGONALS, INCLUDING THE PRINCIPAL */
1986/*          DIAGONAL, REQUIRED TO SPECIFY THE NON-ZERO PORTION OF THE */
1987/*          LOWER TRIANGLE OF THE MATRIX. */
1988
1989/*        A CONTAINS THE LOWER TRIANGLE OF THE SYMMETRIC BAND INPUT */
1990/*          MATRIX STORED AS AN N BY MB ARRAY.  ITS LOWEST SUBDIAGONAL */
1991/*          IS STORED IN THE LAST N+1-MB POSITIONS OF THE FIRST COLUMN, */
1992/*          ITS NEXT SUBDIAGONAL IN THE LAST N+2-MB POSITIONS OF THE */
1993/*          SECOND COLUMN, FURTHER SUBDIAGONALS SIMILARLY, AND FINALLY */
1994/*          ITS PRINCIPAL DIAGONAL IN THE N POSITIONS OF THE LAST COLUMN.
1995*/
1996/*          CONTENTS OF STORAGES NOT PART OF THE MATRIX ARE ARBITRARY. */
1997/*          ON A SUBSEQUENT CALL, ITS OUTPUT CONTENTS FROM THE PREVIOUS */
1998/*          CALL SHOULD BE PASSED. */
1999
2000/*        T SPECIFIES THE SHIFT (OF EIGENVALUES) APPLIED TO THE DIAGONAL
2001*/
2002/*          OF A IN FORMING THE INPUT MATRIX. WHAT IS ACTUALLY DETERMINED
2003*/
2004/*          IS THE EIGENVALUE OF A+TI (I IS THE IDENTITY MATRIX) NEAREST
2005*/
2006/*          TO T.  ON A SUBSEQUENT CALL, THE OUTPUT VALUE OF T FROM THE */
2007/*          PREVIOUS CALL SHOULD BE PASSED IF THE NEXT NEAREST EIGENVALUE
2008*/
2009/*          IS SOUGHT. */
2010
2011/*        R SHOULD BE SPECIFIED AS ZERO ON THE FIRST CALL, AND AS ITS */
2012/*          OUTPUT VALUE FROM THE PREVIOUS CALL ON A SUBSEQUENT CALL. */
2013/*          IT IS USED TO DETERMINE WHEN THE LAST ROW AND COLUMN OF */
2014/*          THE TRANSFORMED BAND MATRIX CAN BE REGARDED AS NEGLIGIBLE. */
2015
2016/*        NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER RV */
2017/*          AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT. */
2018
2019/*     ON OUTPUT */
2020
2021/*        A CONTAINS THE TRANSFORMED BAND MATRIX.  THE MATRIX A+TI */
2022/*          DERIVED FROM THE OUTPUT PARAMETERS IS SIMILAR TO THE */
2023/*          INPUT A+TI TO WITHIN ROUNDING ERRORS.  ITS LAST ROW AND */
2024/*          COLUMN ARE NULL (IF IERR IS ZERO). */
2025
2026/*        T CONTAINS THE COMPUTED EIGENVALUE OF A+TI (IF IERR IS ZERO). */
2027
2028/*        R CONTAINS THE MAXIMUM OF ITS INPUT VALUE AND THE NORM OF THE */
2029/*          LAST COLUMN OF THE INPUT MATRIX A. */
2030
2031/*        IERR IS SET TO */
2032/*          ZERO       FOR NORMAL RETURN, */
2033/*          N          IF THE EIGENVALUE HAS NOT BEEN */
2034/*                     DETERMINED AFTER 30 ITERATIONS. */
2035
2036/*        RV IS A TEMPORARY STORAGE ARRAY OF DIMENSION AT LEAST */
2037/*          (2*MB**2+4*MB-3).  THE FIRST (3*MB-2) LOCATIONS CORRESPOND */
2038/*          TO THE ALGOL ARRAY B, THE NEXT (2*MB-1) LOCATIONS CORRESPOND
2039*/
2040/*          TO THE ALGOL ARRAY H, AND THE FINAL (2*MB**2-MB) LOCATIONS */
2041/*          CORRESPOND TO THE MB BY (2*MB-1) ALGOL ARRAY U. */
2042
2043/*     NOTE. FOR A SUBSEQUENT CALL, N SHOULD BE REPLACED BY N-1, BUT */
2044/*     MB SHOULD NOT BE ALTERED EVEN WHEN IT EXCEEDS THE CURRENT N. */
2045
2046/*     CALLS PYTHAG FOR  DSQRT(A*A + B*B) . */
2047
2048/*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
2049/*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
2050*/
2051
2052/*     THIS VERSION DATED AUGUST 1983. */
2053
2054/*     ------------------------------------------------------------------
2055*/
2056
2057    /* Parameter adjustments */
2058    a_dim1 = *nm;
2059    a_offset = a_dim1 + 1;
2060    a -= a_offset;
2061    --rv;
2062
2063    /* Function Body */
2064    *ierr = 0;
2065    m1 = min(*mb,*n);
2066    m = m1 - 1;
2067    m2 = m + m;
2068    m21 = m2 + 1;
2069    m3 = m21 + m;
2070    m31 = m3 + 1;
2071    m4 = m31 + m2;
2072    mn = m + *n;
2073    mz = *mb - m1;
2074    its = 0;
2075/*     .......... TEST FOR CONVERGENCE .......... */
2076L40:
2077    g = a[*n + *mb * a_dim1];
2078    if (m == 0) {
2079        goto L360;
2080    }
2081    f = 0.;
2082
2083    i_1 = m;
2084    for (k = 1; k <= i_1; ++k) {
2085        mk = k + mz;
2086        f += (d_1 = a[*n + mk * a_dim1], abs(d_1));
2087/* L50: */
2088    }
2089
2090    if (its == 0 && f > *r) {
2091        *r = f;
2092    }
2093    tst1 = *r;
2094    tst2 = tst1 + f;
2095    if (tst2 <= tst1) {
2096        goto L360;
2097    }
2098    if (its == 30) {
2099        goto L1000;
2100    }
2101    ++its;
2102/*     .......... FORM SHIFT FROM BOTTOM 2 BY 2 MINOR .......... */
2103    if (f > *r * .25 && its < 5) {
2104        goto L90;
2105    }
2106    f = a[*n + (*mb - 1) * a_dim1];
2107    if (f == 0.) {
2108        goto L70;
2109    }
2110    q = (a[*n - 1 + *mb * a_dim1] - g) / (f * 2.);
2111    s = pythag_(&q, &c_b141);
2112    g -= f / (q + d_sign(&s, &q));
2113L70:
2114    *t += g;
2115
2116    i_1 = *n;
2117    for (i = 1; i <= i_1; ++i) {
2118/* L80: */
2119        a[i + *mb * a_dim1] -= g;
2120    }
2121
2122L90:
2123    i_1 = m4;
2124    for (k = m31; k <= i_1; ++k) {
2125/* L100: */
2126        rv[k] = 0.;
2127    }
2128
2129    i_1 = mn;
2130    for (ii = 1; ii <= i_1; ++ii) {
2131        i = ii - m;
2132        ni = *n - ii;
2133        if (ni < 0) {
2134            goto L230;
2135        }
2136/*     .......... FORM COLUMN OF SHIFTED MATRIX A-G*I .......... */
2137/* Computing MAX */
2138        i_2 = 1, i_3 = 2 - i;
2139        l = max(i_2,i_3);
2140
2141        i_2 = m3;
2142        for (k = 1; k <= i_2; ++k) {
2143/* L110: */
2144            rv[k] = 0.;
2145        }
2146
2147        i_2 = m1;
2148        for (k = l; k <= i_2; ++k) {
2149            km = k + m;
2150            mk = k + mz;
2151            rv[km] = a[ii + mk * a_dim1];
2152/* L120: */
2153        }
2154
2155        ll = min(m,ni);
2156        if (ll == 0) {
2157            goto L135;
2158        }
2159
2160        i_2 = ll;
2161        for (k = 1; k <= i_2; ++k) {
2162            km = k + m21;
2163            ik = ii + k;
2164            mk = *mb - k;
2165            rv[km] = a[ik + mk * a_dim1];
2166/* L130: */
2167        }
2168/*     .......... PRE-MULTIPLY WITH HOUSEHOLDER REFLECTIONS ..........
2169 */
2170L135:
2171        ll = m2;
2172        imult = 0;
2173/*     .......... MULTIPLICATION PROCEDURE .......... */
2174L140:
2175        kj = m4 - m1;
2176
2177        i_2 = ll;
2178        for (j = 1; j <= i_2; ++j) {
2179            kj += m1;
2180            jm = j + m3;
2181            if (rv[jm] == 0.) {
2182                goto L170;
2183            }
2184            f = 0.;
2185
2186            i_3 = m1;
2187            for (k = 1; k <= i_3; ++k) {
2188                ++kj;
2189                jk = j + k - 1;
2190                f += rv[kj] * rv[jk];
2191/* L150: */
2192            }
2193
2194            f /= rv[jm];
2195            kj -= m1;
2196
2197            i_3 = m1;
2198            for (k = 1; k <= i_3; ++k) {
2199                ++kj;
2200                jk = j + k - 1;
2201                rv[jk] -= rv[kj] * f;
2202/* L160: */
2203            }
2204
2205            kj -= m1;
2206L170:
2207            ;
2208        }
2209
2210        if (imult != 0) {
2211            goto L280;
2212        }
2213/*     .......... HOUSEHOLDER REFLECTION .......... */
2214        f = rv[m21];
2215        s = 0.;
2216        rv[m4] = 0.;
2217        scale = 0.;
2218
2219        i_2 = m3;
2220        for (k = m21; k <= i_2; ++k) {
2221/* L180: */
2222            scale += (d_1 = rv[k], abs(d_1));
2223        }
2224
2225        if (scale == 0.) {
2226            goto L210;
2227        }
2228
2229        i_2 = m3;
2230        for (k = m21; k <= i_2; ++k) {
2231/* L190: */
2232/* Computing 2nd power */
2233            d_1 = rv[k] / scale;
2234            s += d_1 * d_1;
2235        }
2236
2237        s = scale * scale * s;
2238        d_1 = sqrt(s);
2239        g = -d_sign(&d_1, &f);
2240        rv[m21] = g;
2241        rv[m4] = s - f * g;
2242        kj = m4 + m2 * m1 + 1;
2243        rv[kj] = f - g;
2244
2245        i_2 = m1;
2246        for (k = 2; k <= i_2; ++k) {
2247            ++kj;
2248            km = k + m2;
2249            rv[kj] = rv[km];
2250/* L200: */
2251        }
2252/*     .......... SAVE COLUMN OF TRIANGULAR FACTOR R .......... */
2253L210:
2254        i_2 = m1;
2255        for (k = l; k <= i_2; ++k) {
2256            km = k + m;
2257            mk = k + mz;
2258            a[ii + mk * a_dim1] = rv[km];
2259/* L220: */
2260        }
2261
2262L230:
2263/* Computing MAX */
2264        i_2 = 1, i_3 = m1 + 1 - i;
2265        l = max(i_2,i_3);
2266        if (i <= 0) {
2267            goto L300;
2268        }
2269/*     .......... PERFORM ADDITIONAL STEPS .......... */
2270        i_2 = m21;
2271        for (k = 1; k <= i_2; ++k) {
2272/* L240: */
2273            rv[k] = 0.;
2274        }
2275
2276/* Computing MIN */
2277        i_2 = m1, i_3 = ni + m1;
2278        ll = min(i_2,i_3);
2279/*     .......... GET ROW OF TRIANGULAR FACTOR R .......... */
2280        i_2 = ll;
2281        for (kk = 1; kk <= i_2; ++kk) {
2282            k = kk - 1;
2283            km = k + m1;
2284            ik = i + k;
2285            mk = *mb - k;
2286            rv[km] = a[ik + mk * a_dim1];
2287/* L250: */
2288        }
2289/*     .......... POST-MULTIPLY WITH HOUSEHOLDER REFLECTIONS .........
2290. */
2291        ll = m1;
2292        imult = 1;
2293        goto L140;
2294/*     .......... STORE COLUMN OF NEW A MATRIX .......... */
2295L280:
2296        i_2 = m1;
2297        for (k = l; k <= i_2; ++k) {
2298            mk = k + mz;
2299            a[i + mk * a_dim1] = rv[k];
2300/* L290: */
2301        }
2302/*     .......... UPDATE HOUSEHOLDER REFLECTIONS .......... */
2303L300:
2304        if (l > 1) {
2305            --l;
2306        }
2307        kj1 = m4 + l * m1;
2308
2309        i_2 = m2;
2310        for (j = l; j <= i_2; ++j) {
2311            jm = j + m3;
2312            rv[jm] = rv[jm + 1];
2313
2314            i_3 = m1;
2315            for (k = 1; k <= i_3; ++k) {
2316                ++kj1;
2317                kj = kj1 - m1;
2318                rv[kj] = rv[kj1];
2319/* L320: */
2320            }
2321        }
2322
2323/* L350: */
2324    }
2325
2326    goto L40;
2327/*     .......... CONVERGENCE .......... */
2328L360:
2329    *t += g;
2330
2331    i_1 = *n;
2332    for (i = 1; i <= i_1; ++i) {
2333/* L380: */
2334        a[i + *mb * a_dim1] -= g;
2335    }
2336
2337    i_1 = m1;
2338    for (k = 1; k <= i_1; ++k) {
2339        mk = k + mz;
2340        a[*n + mk * a_dim1] = 0.;
2341/* L400: */
2342    }
2343
2344    goto L1001;
2345/*     .......... SET ERROR -- NO CONVERGENCE TO */
2346/*                EIGENVALUE AFTER 30 ITERATIONS .......... */
2347L1000:
2348    *ierr = *n;
2349L1001:
2350    return 0;
2351} /* bqr_ */
2352
2353/* Subroutine */ int cbabk2_(integer *nm, integer *n, integer *low, integer *
2354        igh, doublereal *scale, integer *m, doublereal *zr, doublereal *zi)
2355{
2356    /* System generated locals */
2357    integer zr_dim1, zr_offset, zi_dim1, zi_offset, i_1, i_2;
2358
2359    /* Local variables */
2360    static integer i, j, k;
2361    static doublereal s;
2362    static integer ii;
2363
2364
2365
2366/*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE */
2367/*     CBABK2, WHICH IS A COMPLEX VERSION OF BALBAK, */
2368/*     NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH. */
2369/*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971). */
2370
2371/*     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX GENERAL */
2372/*     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING */
2373/*     BALANCED MATRIX DETERMINED BY  CBAL. */
2374
2375/*     ON INPUT */
2376
2377/*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
2378/*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
2379/*          DIMENSION STATEMENT. */
2380
2381/*        N IS THE ORDER OF THE MATRIX. */
2382
2383/*        LOW AND IGH ARE INTEGERS DETERMINED BY  CBAL. */
2384
2385/*        SCALE CONTAINS INFORMATION DETERMINING THE PERMUTATIONS */
2386/*          AND SCALING FACTORS USED BY  CBAL. */
2387
2388/*        M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED. */
2389
2390/*        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, */
2391/*          RESPECTIVELY, OF THE EIGENVECTORS TO BE */
2392/*          BACK TRANSFORMED IN THEIR FIRST M COLUMNS. */
2393
2394/*     ON OUTPUT */
2395
2396/*        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, */
2397/*          RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS */
2398/*          IN THEIR FIRST M COLUMNS. */
2399
2400/*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
2401/*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
2402*/
2403
2404/*     THIS VERSION DATED AUGUST 1983. */
2405
2406/*     ------------------------------------------------------------------
2407*/
2408
2409    /* Parameter adjustments */
2410    --scale;
2411    zi_dim1 = *nm;
2412    zi_offset = zi_dim1 + 1;
2413    zi -= zi_offset;
2414    zr_dim1 = *nm;
2415    zr_offset = zr_dim1 + 1;
2416    zr -= zr_offset;
2417
2418    /* Function Body */
2419    if (*m == 0) {
2420        goto L200;
2421    }
2422    if (*igh == *low) {
2423        goto L120;
2424    }
2425
2426    i_1 = *igh;
2427    for (i = *low; i <= i_1; ++i) {
2428        s = scale[i];
2429/*     .......... LEFT HAND EIGENVECTORS ARE BACK TRANSFORMED */
2430/*                IF THE FOREGOING STATEMENT IS REPLACED BY */
2431/*                S=1.0D0/SCALE(I). .......... */
2432        i_2 = *m;
2433        for (j = 1; j <= i_2; ++j) {
2434            zr[i + j * zr_dim1] *= s;
2435            zi[i + j * zi_dim1] *= s;
2436/* L100: */
2437        }
2438
2439/* L110: */
2440    }
2441/*     .......... FOR I=LOW-1 STEP -1 UNTIL 1, */
2442/*                IGH+1 STEP 1 UNTIL N DO -- .......... */
2443L120:
2444    i_1 = *n;
2445    for (ii = 1; ii <= i_1; ++ii) {
2446        i = ii;
2447        if (i >= *low && i <= *igh) {
2448            goto L140;
2449        }
2450        if (i < *low) {
2451            i = *low - ii;
2452        }
2453        k = (integer) scale[i];
2454        if (k == i) {
2455            goto L140;
2456        }
2457
2458        i_2 = *m;
2459        for (j = 1; j <= i_2; ++j) {
2460            s = zr[i + j * zr_dim1];
2461            zr[i + j * zr_dim1] = zr[k + j * zr_dim1];
2462            zr[k + j * zr_dim1] = s;
2463            s = zi[i + j * zi_dim1];
2464            zi[i + j * zi_dim1] = zi[k + j * zi_dim1];
2465            zi[k + j * zi_dim1] = s;
2466/* L130: */
2467        }
2468
2469L140:
2470        ;
2471    }
2472
2473L200:
2474    return 0;
2475} /* cbabk2_ */
2476
2477/* Subroutine */ int cbal_(integer *nm, integer *n, doublereal *ar, 
2478        doublereal *ai, integer *low, integer *igh, doublereal *scale)
2479{
2480    /* System generated locals */
2481    integer ar_dim1, ar_offset, ai_dim1, ai_offset, i_1, i_2;
2482    doublereal d_1, d_2;
2483
2484    /* Local variables */
2485    static integer iexc;
2486    static doublereal c, f, g;
2487    static integer i, j, k, l, m;
2488    static doublereal r, s, radix, b2;
2489    static integer jj;
2490    static logical noconv;
2491
2492
2493
2494/*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE */
2495/*     CBALANCE, WHICH IS A COMPLEX VERSION OF BALANCE, */
2496/*     NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH. */
2497/*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971). */
2498
2499/*     THIS SUBROUTINE BALANCES A COMPLEX MATRIX AND ISOLATES */
2500/*     EIGENVALUES WHENEVER POSSIBLE. */
2501
2502/*     ON INPUT */
2503
2504/*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
2505/*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
2506/*          DIMENSION STATEMENT. */
2507
2508/*        N IS THE ORDER OF THE MATRIX. */
2509
2510/*        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */
2511/*          RESPECTIVELY, OF THE COMPLEX MATRIX TO BE BALANCED. */
2512
2513/*     ON OUTPUT */
2514
2515/*        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */
2516/*          RESPECTIVELY, OF THE BALANCED MATRIX. */
2517
2518/*        LOW AND IGH ARE TWO INTEGERS SUCH THAT AR(I,J) AND AI(I,J) */
2519/*          ARE EQUAL TO ZERO IF */
2520/*           (1) I IS GREATER THAN J AND */
2521/*           (2) J=1,...,LOW-1 OR I=IGH+1,...,N. */
2522
2523/*        SCALE CONTAINS INFORMATION DETERMINING THE */
2524/*           PERMUTATIONS AND SCALING FACTORS USED. */
2525
2526/*     SUPPOSE THAT THE PRINCIPAL SUBMATRIX IN ROWS LOW THROUGH IGH */
2527/*     HAS BEEN BALANCED, THAT P(J) DENOTES THE INDEX INTERCHANGED */
2528/*     WITH J DURING THE PERMUTATION STEP, AND THAT THE ELEMENTS */
2529/*     OF THE DIAGONAL MATRIX USED ARE DENOTED BY D(I,J).  THEN */
2530/*        SCALE(J) = P(J),    FOR J = 1,...,LOW-1 */
2531/*                 = D(J,J)       J = LOW,...,IGH */
2532/*                 = P(J)         J = IGH+1,...,N. */
2533/*     THE ORDER IN WHICH THE INTERCHANGES ARE MADE IS N TO IGH+1, */
2534/*     THEN 1 TO LOW-1. */
2535
2536/*     NOTE THAT 1 IS RETURNED FOR IGH IF IGH IS ZERO FORMALLY. */
2537
2538/*     THE ALGOL PROCEDURE EXC CONTAINED IN CBALANCE APPEARS IN */
2539/*     CBAL  IN LINE.  (NOTE THAT THE ALGOL ROLES OF IDENTIFIERS */
2540/*     K,L HAVE BEEN REVERSED.) */
2541
2542/*     ARITHMETIC IS REAL THROUGHOUT. */
2543
2544/*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
2545/*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
2546*/
2547
2548/*     THIS VERSION DATED AUGUST 1983. */
2549
2550/*     ------------------------------------------------------------------
2551*/
2552
2553    /* Parameter adjustments */
2554    --scale;
2555    ai_dim1 = *nm;
2556    ai_offset = ai_dim1 + 1;
2557    ai -= ai_offset;
2558    ar_dim1 = *nm;
2559    ar_offset = ar_dim1 + 1;
2560    ar -= ar_offset;
2561
2562    /* Function Body */
2563    radix = 16.;
2564
2565    b2 = radix * radix;
2566    k = 1;
2567    l = *n;
2568    goto L100;
2569/*     .......... IN-LINE PROCEDURE FOR ROW AND */
2570/*                COLUMN EXCHANGE .......... */
2571L20:
2572    scale[m] = (doublereal) j;
2573    if (j == m) {
2574        goto L50;
2575    }
2576
2577    i_1 = l;
2578    for (i = 1; i <= i_1; ++i) {
2579        f = ar[i + j * ar_dim1];
2580        ar[i + j * ar_dim1] = ar[i + m * ar_dim1];
2581        ar[i + m * ar_dim1] = f;
2582        f = ai[i + j * ai_dim1];
2583        ai[i + j * ai_dim1] = ai[i + m * ai_dim1];
2584        ai[i + m * ai_dim1] = f;
2585/* L30: */
2586    }
2587
2588    i_1 = *n;
2589    for (i = k; i <= i_1; ++i) {
2590        f = ar[j + i * ar_dim1];
2591        ar[j + i * ar_dim1] = ar[m + i * ar_dim1];
2592        ar[m + i * ar_dim1] = f;
2593        f = ai[j + i * ai_dim1];
2594        ai[j + i * ai_dim1] = ai[m + i * ai_dim1];
2595        ai[m + i * ai_dim1] = f;
2596/* L40: */
2597    }
2598
2599L50:
2600    switch (iexc) {
2601        case 1:  goto L80;
2602        case 2:  goto L130;
2603    }
2604/*     .......... SEARCH FOR ROWS ISOLATING AN EIGENVALUE */
2605/*                AND PUSH THEM DOWN .......... */
2606L80:
2607    if (l == 1) {
2608        goto L280;
2609    }
2610    --l;
2611/*     .......... FOR J=L STEP -1 UNTIL 1 DO -- .......... */
2612L100:
2613    i_1 = l;
2614    for (jj = 1; jj <= i_1; ++jj) {
2615        j = l + 1 - jj;
2616
2617        i_2 = l;
2618        for (i = 1; i <= i_2; ++i) {
2619            if (i == j) {
2620                goto L110;
2621            }
2622            if (ar[j + i * ar_dim1] != 0. || ai[j + i * ai_dim1] != 0.) {
2623                goto L120;
2624            }
2625L110:
2626            ;
2627        }
2628
2629        m = l;
2630        iexc = 1;
2631        goto L20;
2632L120:
2633        ;
2634    }
2635
2636    goto L140;
2637/*     .......... SEARCH FOR COLUMNS ISOLATING AN EIGENVALUE */
2638/*                AND PUSH THEM LEFT .......... */
2639L130:
2640    ++k;
2641
2642L140:
2643    i_1 = l;
2644    for (j = k; j <= i_1; ++j) {
2645
2646        i_2 = l;
2647        for (i = k; i <= i_2; ++i) {
2648            if (i == j) {
2649                goto L150;
2650            }
2651            if (ar[i + j * ar_dim1] != 0. || ai[i + j * ai_dim1] != 0.) {
2652                goto L170;
2653            }
2654L150:
2655            ;
2656        }
2657
2658        m = k;
2659        iexc = 2;
2660        goto L20;
2661L170:
2662        ;
2663    }
2664/*     .......... NOW BALANCE THE SUBMATRIX IN ROWS K TO L .......... */
2665    i_1 = l;
2666    for (i = k; i <= i_1; ++i) {
2667/* L180: */
2668        scale[i] = 1.;
2669    }
2670/*     .......... ITERATIVE LOOP FOR NORM REDUCTION .......... */
2671L190:
2672    noconv = FALSE_;
2673
2674    i_1 = l;
2675    for (i = k; i <= i_1; ++i) {
2676        c = 0.;
2677        r = 0.;
2678
2679        i_2 = l;
2680        for (j = k; j <= i_2; ++j) {
2681            if (j == i) {
2682                goto L200;
2683            }
2684            c = c + (d_1 = ar[j + i * ar_dim1], abs(d_1)) + (d_2 = ai[j + 
2685                    i * ai_dim1], abs(d_2));
2686            r = r + (d_1 = ar[i + j * ar_dim1], abs(d_1)) + (d_2 = ai[i + 
2687                    j * ai_dim1], abs(d_2));
2688L200:
2689            ;
2690        }
2691/*     .......... GUARD AGAINST ZERO C OR R DUE TO UNDERFLOW .........
2692. */
2693        if (c == 0. || r == 0.) {
2694            goto L270;
2695        }
2696        g = r / radix;
2697        f = 1.;
2698        s = c + r;
2699L210:
2700        if (c >= g) {
2701            goto L220;
2702        }
2703        f *= radix;
2704        c *= b2;
2705        goto L210;
2706L220:
2707        g = r * radix;
2708L230:
2709        if (c < g) {
2710            goto L240;
2711        }
2712        f /= radix;
2713        c /= b2;
2714        goto L230;
2715/*     .......... NOW BALANCE .......... */
2716L240:
2717        if ((c + r) / f >= s * .95) {
2718            goto L270;
2719        }
2720        g = 1. / f;
2721        scale[i] *= f;
2722        noconv = TRUE_;
2723
2724        i_2 = *n;
2725        for (j = k; j <= i_2; ++j) {
2726            ar[i + j * ar_dim1] *= g;
2727            ai[i + j * ai_dim1] *= g;
2728/* L250: */
2729        }
2730
2731        i_2 = l;
2732        for (j = 1; j <= i_2; ++j) {
2733            ar[j + i * ar_dim1] *= f;
2734            ai[j + i * ai_dim1] *= f;
2735/* L260: */
2736        }
2737
2738L270:
2739        ;
2740    }
2741
2742    if (noconv) {
2743        goto L190;
2744    }
2745
2746L280:
2747    *low = k;
2748    *igh = l;
2749    return 0;
2750} /* cbal_ */
2751
2752/* Subroutine */ int cg_(integer *nm, integer *n, doublereal *ar, doublereal *
2753        ai, doublereal *wr, doublereal *wi, integer *matz, doublereal *zr, 
2754        doublereal *zi, doublereal *fv1, doublereal *fv2, doublereal *fv3, 
2755        integer *ierr)
2756{
2757    /* System generated locals */
2758    integer ar_dim1, ar_offset, ai_dim1, ai_offset, zr_dim1, zr_offset, 
2759            zi_dim1, zi_offset;
2760
2761    /* Local variables */
2762    extern /* Subroutine */ int cbal_(integer *, integer *, doublereal *, 
2763            doublereal *, integer *, integer *, doublereal *), corth_(integer
2764            *, integer *, integer *, integer *, doublereal *, doublereal *, 
2765            doublereal *, doublereal *), comqr_(integer *, integer *, integer
2766            *, integer *, doublereal *, doublereal *, doublereal *, 
2767            doublereal *, integer *), cbabk2_(integer *, integer *, integer *,
2768             integer *, doublereal *, integer *, doublereal *, doublereal *), 
2769            comqr2_(integer *, integer *, integer *, integer *, doublereal *, 
2770            doublereal *, doublereal *, doublereal *, doublereal *, 
2771            doublereal *, doublereal *, doublereal *, integer *);
2772    static integer is1, is2;
2773
2774
2775
2776/*     THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF */
2777/*     SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK) */
2778/*     TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED) */
2779/*     OF A COMPLEX GENERAL MATRIX. */
2780
2781/*     ON INPUT */
2782
2783/*        NM  MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL */
2784/*        ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
2785/*        DIMENSION STATEMENT. */
2786
2787/*        N  IS THE ORDER OF THE MATRIX  A=(AR,AI). */
2788
2789/*        AR  AND  AI  CONTAIN THE REAL AND IMAGINARY PARTS, */
2790/*        RESPECTIVELY, OF THE COMPLEX GENERAL MATRIX. */
2791
2792/*        MATZ  IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF */
2793/*        ONLY EIGENVALUES ARE DESIRED.  OTHERWISE IT IS SET TO */
2794/*        ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS. */
2795
2796/*     ON OUTPUT */
2797
2798/*        WR  AND  WI  CONTAIN THE REAL AND IMAGINARY PARTS, */
2799/*        RESPECTIVELY, OF THE EIGENVALUES. */
2800
2801/*        ZR  AND  ZI  CONTAIN THE REAL AND IMAGINARY PARTS, */
2802/*        RESPECTIVELY, OF THE EIGENVECTORS IF MATZ IS NOT ZERO. */
2803
2804/*        IERR  IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR */
2805/*           COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR COMQR */
2806/*           AND COMQR2.  THE NORMAL COMPLETION CODE IS ZERO. */
2807
2808/*        FV1, FV2, AND  FV3  ARE TEMPORARY STORAGE ARRAYS. */
2809
2810/*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
2811/*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
2812*/
2813
2814/*     THIS VERSION DATED AUGUST 1983. */
2815
2816/*     ------------------------------------------------------------------
2817*/
2818
2819    /* Parameter adjustments */
2820    --fv3;
2821    --fv2;
2822    --fv1;
2823    zi_dim1 = *nm;
2824    zi_offset = zi_dim1 + 1;
2825    zi -= zi_offset;
2826    zr_dim1 = *nm;
2827    zr_offset = zr_dim1 + 1;
2828    zr -= zr_offset;
2829    --wi;
2830    --wr;
2831    ai_dim1 = *nm;
2832    ai_offset = ai_dim1 + 1;
2833    ai -= ai_offset;
2834    ar_dim1 = *nm;
2835    ar_offset = ar_dim1 + 1;
2836    ar -= ar_offset;
2837
2838    /* Function Body */
2839    if (*n <= *nm) {
2840        goto L10;
2841    }
2842    *ierr = *n * 10;
2843    goto L50;
2844
2845L10:
2846    cbal_(nm, n, &ar[ar_offset], &ai[ai_offset], &is1, &is2, &fv1[1]);
2847    corth_(nm, n, &is1, &is2, &ar[ar_offset], &ai[ai_offset], &fv2[1], &fv3[1]
2848            );
2849    if (*matz != 0) {
2850        goto L20;
2851    }
2852/*     .......... FIND EIGENVALUES ONLY .......... */
2853    comqr_(nm, n, &is1, &is2, &ar[ar_offset], &ai[ai_offset], &wr[1], &wi[1], 
2854            ierr);
2855    goto L50;
2856/*     .......... FIND BOTH EIGENVALUES AND EIGENVECTORS .......... */
2857L20:
2858    comqr2_(nm, n, &is1, &is2, &fv2[1], &fv3[1], &ar[ar_offset], &ai[
2859            ai_offset], &wr[1], &wi[1], &zr[zr_offset], &zi[zi_offset], ierr);
2860    if (*ierr != 0) {
2861        goto L50;
2862    }
2863    cbabk2_(nm, n, &is1, &is2, &fv1[1], n, &zr[zr_offset], &zi[zi_offset]);
2864L50:
2865    return 0;
2866} /* cg_ */
2867
2868/* Subroutine */ int ch_(integer *nm, integer *n, doublereal *ar, doublereal *
2869        ai, doublereal *w, integer *matz, doublereal *zr, doublereal *zi, 
2870        doublereal *fv1, doublereal *fv2, doublereal *fm1, integer *ierr)
2871{
2872    /* System generated locals */
2873    integer ar_dim1, ar_offset, ai_dim1, ai_offset, zr_dim1, zr_offset, 
2874            zi_dim1, zi_offset, i_1, i_2;
2875
2876    /* Local variables */
2877    static integer i, j;
2878    extern /* Subroutine */ int htridi_(integer *, integer *, doublereal *, 
2879            doublereal *, doublereal *, doublereal *, doublereal *, 
2880            doublereal *), htribk_(integer *, integer *, doublereal *, 
2881            doublereal *, doublereal *, integer *, doublereal *, doublereal *)
2882            , tqlrat_(integer *, doublereal *, doublereal *, integer *), 
2883            tql2_(integer *, integer *, doublereal *, doublereal *, 
2884            doublereal *, integer *);
2885
2886
2887
2888/*     THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF */
2889/*     SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK) */
2890/*     TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED) */
2891/*     OF A COMPLEX HERMITIAN MATRIX. */
2892
2893/*     ON INPUT */
2894
2895/*        NM  MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL */
2896/*        ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
2897/*        DIMENSION STATEMENT. */
2898
2899/*        N  IS THE ORDER OF THE MATRIX  A=(AR,AI). */
2900
2901/*        AR  AND  AI  CONTAIN THE REAL AND IMAGINARY PARTS, */
2902/*        RESPECTIVELY, OF THE COMPLEX HERMITIAN MATRIX. */
2903
2904/*        MATZ  IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF */
2905/*        ONLY EIGENVALUES ARE DESIRED.  OTHERWISE IT IS SET TO */
2906/*        ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS. */
2907
2908/*     ON OUTPUT */
2909
2910/*        W  CONTAINS THE EIGENVALUES IN ASCENDING ORDER. */
2911
2912/*        ZR  AND  ZI  CONTAIN THE REAL AND IMAGINARY PARTS, */
2913/*        RESPECTIVELY, OF THE EIGENVECTORS IF MATZ IS NOT ZERO. */
2914
2915/*        IERR  IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR */
2916/*           COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR TQLRAT */
2917/*           AND TQL2.  THE NORMAL COMPLETION CODE IS ZERO. */
2918
2919/*        FV1, FV2, AND  FM1  ARE TEMPORARY STORAGE ARRAYS. */
2920
2921/*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
2922/*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
2923*/
2924
2925/*     THIS VERSION DATED AUGUST 1983. */
2926
2927/*     ------------------------------------------------------------------
2928*/
2929
2930    /* Parameter adjustments */
2931    fm1 -= 3;
2932    --fv2;
2933    --fv1;
2934    zi_dim1 = *nm;
2935    zi_offset = zi_dim1 + 1;
2936    zi -= zi_offset;
2937    zr_dim1 = *nm;
2938    zr_offset = zr_dim1 + 1;
2939    zr -= zr_offset;
2940    --w;
2941    ai_dim1 = *nm;
2942    ai_offset = ai_dim1 + 1;
2943    ai -= ai_offset;
2944    ar_dim1 = *nm;
2945    ar_offset = ar_dim1 + 1;
2946    ar -= ar_offset;
2947
2948    /* Function Body */
2949    if (*n <= *nm) {
2950        goto L10;
2951    }
2952    *ierr = *n * 10;
2953    goto L50;
2954
2955L10:
2956    htridi_(nm, n, &ar[ar_offset], &ai[ai_offset], &w[1], &fv1[1], &fv2[1], &
2957            fm1[3]);
2958    if (*matz != 0) {
2959        goto L20;
2960    }
2961/*     .......... FIND EIGENVALUES ONLY .......... */
2962    tqlrat_(n, &w[1], &fv2[1], ierr);
2963    goto L50;
2964/*     .......... FIND BOTH EIGENVALUES AND EIGENVECTORS .......... */
2965L20:
2966    i_1 = *n;
2967    for (i = 1; i <= i_1; ++i) {
2968
2969        i_2 = *n;
2970        for (j = 1; j <= i_2; ++j) {
2971            zr[j + i * zr_dim1] = 0.;
2972/* L30: */
2973        }
2974
2975        zr[i + i * zr_dim1] = 1.;
2976/* L40: */
2977    }
2978
2979    tql2_(nm, n, &w[1], &fv1[1], &zr[zr_offset], ierr);
2980    if (*ierr != 0) {
2981        goto L50;
2982    }
2983    htribk_(nm, n, &ar[ar_offset], &ai[ai_offset], &fm1[3], n, &zr[zr_offset],
2984             &zi[zi_offset]);
2985L50:
2986    return 0;
2987} /* ch_ */
2988
2989/* Subroutine */ int cinvit_(integer *nm, integer *n, doublereal *ar, 
2990        doublereal *ai, doublereal *wr, doublereal *wi, logical *select, 
2991        integer *mm, integer *m, doublereal *zr, doublereal *zi, integer *
2992        ierr, doublereal *rm1, doublereal *rm2, doublereal *rv1, doublereal *
2993        rv2)
2994{
2995    /* System generated locals */
2996    integer ar_dim1, ar_offset, ai_dim1, ai_offset, zr_dim1, zr_offset, 
2997            zi_dim1, zi_offset, rm1_dim1, rm1_offset, rm2_dim1, rm2_offset, 
2998            i_1, i_2, i_3;
2999    doublereal d_1, d_2;
3000
3001    /* Builtin functions */
3002    double sqrt(doublereal);
3003
3004    /* Local variables */
3005    extern /* Subroutine */ int cdiv_(doublereal *, doublereal *, doublereal *
3006            , doublereal *, doublereal *, doublereal *);
3007    static doublereal norm;
3008    static integer i, j, k, s;
3009    static doublereal x, y, normv;
3010    static integer ii;
3011    static doublereal ilambd;
3012    static integer mp, uk;
3013    static doublereal rlambd;
3014    extern doublereal pythag_(doublereal *, doublereal *), epslon_(doublereal
3015            *);
3016    static integer km1, ip1;
3017    static doublereal growto, ukroot;
3018    static integer its;
3019    static doublereal eps3;
3020
3021
3022
3023/*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE CX INVIT */
3024/*     BY PETERS AND WILKINSON. */
3025/*     HANDBOOK FOR AUTO. COMP. VOL.II-LINEAR ALGEBRA, 418-439(1971). */
3026
3027/*     THIS SUBROUTINE FINDS THOSE EIGENVECTORS OF A COMPLEX UPPER */
3028/*     HESSENBERG MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES, */
3029/*     USING INVERSE ITERATION. */
3030
3031/*     ON INPUT */
3032
3033/*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
3034/*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
3035/*          DIMENSION STATEMENT. */
3036
3037/*        N IS THE ORDER OF THE MATRIX. */
3038
3039/*        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */
3040/*          RESPECTIVELY, OF THE HESSENBERG MATRIX. */
3041
3042/*        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, RESPECTIVELY, */
3043/*          OF THE EIGENVALUES OF THE MATRIX.  THE EIGENVALUES MUST BE */
3044/*          STORED IN A MANNER IDENTICAL TO THAT OF SUBROUTINE  COMLR, */
3045/*          WHICH RECOGNIZES POSSIBLE SPLITTING OF THE MATRIX. */
3046
3047/*        SELECT SPECIFIES THE EIGENVECTORS TO BE FOUND.  THE */
3048/*          EIGENVECTOR CORRESPONDING TO THE J-TH EIGENVALUE IS */
3049/*          SPECIFIED BY SETTING SELECT(J) TO .TRUE.. */
3050
3051/*        MM SHOULD BE SET TO AN UPPER BOUND FOR THE NUMBER OF */
3052/*          EIGENVECTORS TO BE FOUND. */
3053
3054/*     ON OUTPUT */
3055
3056/*        AR, AI, WI, AND SELECT ARE UNALTERED. */
3057
3058/*        WR MAY HAVE BEEN ALTERED SINCE CLOSE EIGENVALUES ARE PERTURBED
3059*/
3060/*          SLIGHTLY IN SEARCHING FOR INDEPENDENT EIGENVECTORS. */
3061
3062/*        M IS THE NUMBER OF EIGENVECTORS ACTUALLY FOUND. */
3063
3064/*        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, RESPECTIVELY, */
3065/*          OF THE EIGENVECTORS.  THE EIGENVECTORS ARE NORMALIZED */
3066/*          SO THAT THE COMPONENT OF LARGEST MAGNITUDE IS 1. */
3067/*          ANY VECTOR WHICH FAILS THE ACCEPTANCE TEST IS SET TO ZERO. */
3068
3069/*        IERR IS SET TO */
3070/*          ZERO       FOR NORMAL RETURN, */
3071/*          -(2*N+1)   IF MORE THAN MM EIGENVECTORS HAVE BEEN SPECIFIED,
3072*/
3073/*          -K         IF THE ITERATION CORRESPONDING TO THE K-TH */
3074/*                     VALUE FAILS, */
3075/*          -(N+K)     IF BOTH ERROR SITUATIONS OCCUR. */
3076
3077/*        RM1, RM2, RV1, AND RV2 ARE TEMPORARY STORAGE ARRAYS. */
3078
3079/*     THE ALGOL PROCEDURE GUESSVEC APPEARS IN CINVIT IN LINE. */
3080
3081/*     CALLS CDIV FOR COMPLEX DIVISION. */
3082/*     CALLS PYTHAG FOR  DSQRT(A*A + B*B) . */
3083
3084/*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
3085/*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
3086*/
3087
3088/*     THIS VERSION DATED AUGUST 1983. */
3089
3090/*     ------------------------------------------------------------------
3091*/
3092
3093    /* Parameter adjustments */
3094    --rv2;
3095    --rv1;
3096    rm2_dim1 = *n;
3097    rm2_offset = rm2_dim1 + 1;
3098    rm2 -= rm2_offset;
3099    rm1_dim1 = *n;
3100    rm1_offset = rm1_dim1 + 1;
3101    rm1 -= rm1_offset;
3102    --select;
3103    --wi;
3104    --wr;
3105    ai_dim1 = *nm;
3106    ai_offset = ai_dim1 + 1;
3107    ai -= ai_offset;
3108    ar_dim1 = *nm;
3109    ar_offset = ar_dim1 + 1;
3110    ar -= ar_offset;
3111    zi_dim1 = *nm;
3112    zi_offset = zi_dim1 + 1;
3113    zi -= zi_offset;
3114    zr_dim1 = *nm;
3115    zr_offset = zr_dim1 + 1;
3116    zr -= zr_offset;
3117
3118    /* Function Body */
3119    *ierr = 0;
3120    uk = 0;
3121    s = 1;
3122
3123    i_1 = *n;
3124    for (k = 1; k <= i_1; ++k) {
3125        if (! select[k]) {
3126            goto L980;
3127        }
3128        if (s > *mm) {
3129            goto L1000;
3130        }
3131        if (uk >= k) {
3132            goto L200;
3133        }
3134/*     .......... CHECK FOR POSSIBLE SPLITTING .......... */
3135        i_2 = *n;
3136        for (uk = k; uk <= i_2; ++uk) {
3137            if (uk == *n) {
3138                goto L140;
3139            }
3140            if (ar[uk + 1 + uk * ar_dim1] == 0. && ai[uk + 1 + uk * ai_dim1] 
3141                    == 0.) {
3142                goto L140;
3143            }
3144/* L120: */
3145        }
3146/*     .......... COMPUTE INFINITY NORM OF LEADING UK BY UK */
3147/*                (HESSENBERG) MATRIX .......... */
3148L140:
3149        norm = 0.;
3150        mp = 1;
3151
3152        i_2 = uk;
3153        for (i = 1; i <= i_2; ++i) {
3154            x = 0.;
3155
3156            i_3 = uk;
3157            for (j = mp; j <= i_3; ++j) {
3158/* L160: */
3159                x += pythag_(&ar[i + j * ar_dim1], &ai[i + j * ai_dim1]);
3160            }
3161
3162            if (x > norm) {
3163                norm = x;
3164            }
3165            mp = i;
3166/* L180: */
3167        }
3168/*     .......... EPS3 REPLACES ZERO PIVOT IN DECOMPOSITION */
3169/*                AND CLOSE ROOTS ARE MODIFIED BY EPS3 .......... */
3170        if (norm == 0.) {
3171            norm = 1.;
3172        }
3173        eps3 = epslon_(&norm);
3174/*     .......... GROWTO IS THE CRITERION FOR GROWTH .......... */
3175        ukroot = (doublereal) uk;
3176        ukroot = sqrt(ukroot);
3177        growto = .1 / ukroot;
3178L200:
3179        rlambd = wr[k];
3180        ilambd = wi[k];
3181        if (k == 1) {
3182            goto L280;
3183        }
3184        km1 = k - 1;
3185        goto L240;
3186/*     .......... PERTURB EIGENVALUE IF IT IS CLOSE */
3187/*                TO ANY PREVIOUS EIGENVALUE .......... */
3188L220:
3189        rlambd += eps3;
3190/*     .......... FOR I=K-1 STEP -1 UNTIL 1 DO -- .......... */
3191L240:
3192        i_2 = km1;
3193        for (ii = 1; ii <= i_2; ++ii) {
3194            i = k - ii;
3195            if (select[i] && (d_1 = wr[i] - rlambd, abs(d_1)) < eps3 && (
3196                    d_2 = wi[i] - ilambd, abs(d_2)) < eps3) {
3197                goto L220;
3198            }
3199/* L260: */
3200        }
3201
3202        wr[k] = rlambd;
3203/*     .......... FORM UPPER HESSENBERG (AR,AI)-(RLAMBD,ILAMBD)*I */
3204/*                AND INITIAL COMPLEX VECTOR .......... */
3205L280:
3206        mp = 1;
3207
3208        i_2 = uk;
3209        for (i = 1; i <= i_2; ++i) {
3210
3211            i_3 = uk;
3212            for (j = mp; j <= i_3; ++j) {
3213                rm1[i + j * rm1_dim1] = ar[i + j * ar_dim1];
3214                rm2[i + j * rm2_dim1] = ai[i + j * ai_dim1];
3215/* L300: */
3216            }
3217
3218            rm1[i + i * rm1_dim1] -= rlambd;
3219            rm2[i + i * rm2_dim1] -= ilambd;
3220            mp = i;
3221            rv1[i] = eps3;
3222/* L320: */
3223        }
3224/*     .......... TRIANGULAR DECOMPOSITION WITH INTERCHANGES, */
3225/*                REPLACING ZERO PIVOTS BY EPS3 .......... */
3226        if (uk == 1) {
3227            goto L420;
3228        }
3229
3230        i_2 = uk;
3231        for (i = 2; i <= i_2; ++i) {
3232            mp = i - 1;
3233            if (pythag_(&rm1[i + mp * rm1_dim1], &rm2[i + mp * rm2_dim1]) <= 
3234                    pythag_(&rm1[mp + mp * rm1_dim1], &rm2[mp + mp * rm2_dim1]
3235                    )) {
3236                goto L360;
3237            }
3238
3239            i_3 = uk;
3240            for (j = mp; j <= i_3; ++j) {
3241                y = rm1[i + j * rm1_dim1];
3242                rm1[i + j * rm1_dim1] = rm1[mp + j * rm1_dim1];
3243                rm1[mp + j * rm1_dim1] = y;
3244                y = rm2[i + j * rm2_dim1];
3245                rm2[i + j * rm2_dim1] = rm2[mp + j * rm2_dim1];
3246                rm2[mp + j * rm2_dim1] = y;
3247/* L340: */
3248            }
3249
3250L360:
3251            if (rm1[mp + mp * rm1_dim1] == 0. && rm2[mp + mp * rm2_dim1] == 
3252                    0.) {
3253                rm1[mp + mp * rm1_dim1] = eps3;
3254            }
3255            cdiv_(&rm1[i + mp * rm1_dim1], &rm2[i + mp * rm2_dim1], &rm1[mp + 
3256                    mp * rm1_dim1], &rm2[mp + mp * rm2_dim1], &x, &y);
3257            if (x == 0. && y == 0.) {
3258                goto L400;
3259            }
3260
3261            i_3 = uk;
3262            for (j = i; j <= i_3; ++j) {
3263                rm1[i + j * rm1_dim1] = rm1[i + j * rm1_dim1] - x * rm1[mp + 
3264                        j * rm1_dim1] + y * rm2[mp + j * rm2_dim1];
3265                rm2[i + j * rm2_dim1] = rm2[i + j * rm2_dim1] - x * rm2[mp + 
3266                        j * rm2_dim1] - y * rm1[mp + j * rm1_dim1];
3267/* L380: */
3268            }
3269
3270L400:
3271            ;
3272        }
3273
3274L420:
3275        if (rm1[uk + uk * rm1_dim1] == 0. && rm2[uk + uk * rm2_dim1] == 0.) {
3276            rm1[uk + uk * rm1_dim1] = eps3;
3277        }
3278        its = 0;
3279/*     .......... BACK SUBSTITUTION */
3280/*                FOR I=UK STEP -1 UNTIL 1 DO -- .......... */
3281L660:
3282        i_2 = uk;
3283        for (ii = 1; ii <= i_2; ++ii) {
3284            i = uk + 1 - ii;
3285            x = rv1[i];
3286            y = 0.;
3287            if (i == uk) {
3288                goto L700;
3289            }
3290            ip1 = i + 1;
3291
3292            i_3 = uk;
3293            for (j = ip1; j <= i_3; ++j) {
3294                x = x - rm1[i + j * rm1_dim1] * rv1[j] + rm2[i + j * rm2_dim1]
3295                         * rv2[j];
3296                y = y - rm1[i + j * rm1_dim1] * rv2[j] - rm2[i + j * rm2_dim1]
3297                         * rv1[j];
3298/* L680: */
3299            }
3300
3301L700:
3302            cdiv_(&x, &y, &rm1[i + i * rm1_dim1], &rm2[i + i * rm2_dim1], &
3303                    rv1[i], &rv2[i]);
3304/* L720: */
3305        }
3306/*     .......... ACCEPTANCE TEST FOR EIGENVECTOR */
3307/*                AND NORMALIZATION .......... */
3308        ++its;
3309        norm = 0.;
3310        normv = 0.;
3311
3312        i_2 = uk;
3313        for (i = 1; i <= i_2; ++i) {
3314            x = pythag_(&rv1[i], &rv2[i]);
3315            if (normv >= x) {
3316                goto L760;
3317            }
3318            normv = x;
3319            j = i;
3320L760:
3321            norm += x;
3322/* L780: */
3323        }
3324
3325        if (norm < growto) {
3326            goto L840;
3327        }
3328/*     .......... ACCEPT VECTOR .......... */
3329        x = rv1[j];
3330        y = rv2[j];
3331
3332        i_2 = uk;
3333        for (i = 1; i <= i_2; ++i) {
3334            cdiv_(&rv1[i], &rv2[i], &x, &y, &zr[i + s * zr_dim1], &zi[i + s * 
3335                    zi_dim1]);
3336/* L820: */
3337        }
3338
3339        if (uk == *n) {
3340            goto L940;
3341        }
3342        j = uk + 1;
3343        goto L900;
3344/*     .......... IN-LINE PROCEDURE FOR CHOOSING */
3345/*                A NEW STARTING VECTOR .......... */
3346L840:
3347        if (its >= uk) {
3348            goto L880;
3349        }
3350        x = ukroot;
3351        y = eps3 / (x + 1.);
3352        rv1[1] = eps3;
3353
3354        i_2 = uk;
3355        for (i = 2; i <= i_2; ++i) {
3356/* L860: */
3357            rv1[i] = y;
3358        }
3359
3360        j = uk - its + 1;
3361        rv1[j] -= eps3 * x;
3362        goto L660;
3363/*     .......... SET ERROR -- UNACCEPTED EIGENVECTOR .......... */
3364L880:
3365        j = 1;
3366        *ierr = -k;
3367/*     .......... SET REMAINING VECTOR COMPONENTS TO ZERO ..........
3368*/
3369L900:
3370        i_2 = *n;
3371        for (i = j; i <= i_2; ++i) {
3372            zr[i + s * zr_dim1] = 0.;
3373            zi[i + s * zi_dim1] = 0.;
3374/* L920: */
3375        }
3376
3377L940:
3378        ++s;
3379L980:
3380        ;
3381    }
3382
3383    goto L1001;
3384/*     .......... SET ERROR -- UNDERESTIMATE OF EIGENVECTOR */
3385/*                SPACE REQUIRED .......... */
3386L1000:
3387    if (*ierr != 0) {
3388        *ierr -= *n;
3389    }
3390    if (*ierr == 0) {
3391        *ierr = -((*n << 1) + 1);
3392    }
3393L1001:
3394    *m = s - 1;
3395    return 0;
3396} /* cinvit_ */
3397
3398/* Subroutine */ int combak_(integer *nm, integer *low, integer *igh, 
3399        doublereal *ar, doublereal *ai, integer *int_, integer *m, 
3400        doublereal *zr, doublereal *zi)
3401{
3402    /* System generated locals */
3403    integer ar_dim1, ar_offset, ai_dim1, ai_offset, zr_dim1, zr_offset, 
3404            zi_dim1, zi_offset, i_1, i_2, i_3;
3405
3406    /* Local variables */
3407    static integer i, j, la, mm, mp;
3408    static doublereal xi, xr;
3409    static integer kp1, mp1;
3410
3411
3412
3413/*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE COMBAK, */
3414/*     NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON. */
3415/*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). */
3416
3417/*     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX GENERAL */
3418/*     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING */
3419/*     UPPER HESSENBERG MATRIX DETERMINED BY  COMHES. */
3420
3421/*     ON INPUT */
3422
3423/*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
3424/*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
3425/*          DIMENSION STATEMENT. */
3426
3427/*        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
3428/*          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED, */
3429/*          SET LOW=1 AND IGH EQUAL TO THE ORDER OF THE MATRIX. */
3430
3431/*        AR AND AI CONTAIN THE MULTIPLIERS WHICH WERE USED IN THE */
3432/*          REDUCTION BY  COMHES  IN THEIR LOWER TRIANGLES */
3433/*          BELOW THE SUBDIAGONAL. */
3434
3435/*        INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS */
3436/*          INTERCHANGED IN THE REDUCTION BY  COMHES. */
3437/*          ONLY ELEMENTS LOW THROUGH IGH ARE USED. */
3438
3439/*        M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED. */
3440
3441/*        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, */
3442/*          RESPECTIVELY, OF THE EIGENVECTORS TO BE */
3443/*          BACK TRANSFORMED IN THEIR FIRST M COLUMNS. */
3444
3445/*     ON OUTPUT */
3446
3447/*        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, */
3448/*          RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS */
3449/*          IN THEIR FIRST M COLUMNS. */
3450
3451/*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
3452/*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
3453*/
3454
3455/*     THIS VERSION DATED AUGUST 1983. */
3456
3457/*     ------------------------------------------------------------------
3458*/
3459
3460    /* Parameter adjustments */
3461    --int_;
3462    ai_dim1 = *nm;
3463    ai_offset = ai_dim1 + 1;
3464    ai -= ai_offset;
3465    ar_dim1 = *nm;
3466    ar_offset = ar_dim1 + 1;
3467    ar -= ar_offset;
3468    zi_dim1 = *nm;
3469    zi_offset = zi_dim1 + 1;
3470    zi -= zi_offset;
3471    zr_dim1 = *nm;
3472    zr_offset = zr_dim1 + 1;
3473    zr -= zr_offset;
3474
3475    /* Function Body */
3476    if (*m == 0) {
3477        goto L200;
3478    }
3479    la = *igh - 1;
3480    kp1 = *low + 1;
3481    if (la < kp1) {
3482        goto L200;
3483    }
3484/*     .......... FOR MP=IGH-1 STEP -1 UNTIL LOW+1 DO -- .......... */
3485    i_1 = la;
3486    for (mm = kp1; mm <= i_1; ++mm) {
3487        mp = *low + *igh - mm;
3488        mp1 = mp + 1;
3489
3490        i_2 = *igh;
3491        for (i = mp1; i <= i_2; ++i) {
3492            xr = ar[i + (mp - 1) * ar_dim1];
3493            xi = ai[i + (mp - 1) * ai_dim1];
3494            if (xr == 0. && xi == 0.) {
3495                goto L110;
3496            }
3497
3498            i_3 = *m;
3499            for (j = 1; j <= i_3; ++j) {
3500                zr[i + j * zr_dim1] = zr[i + j * zr_dim1] + xr * zr[mp + j * 
3501                        zr_dim1] - xi * zi[mp + j * zi_dim1];
3502                zi[i + j * zi_dim1] = zi[i + j * zi_dim1] + xr * zi[mp + j * 
3503                        zi_dim1] + xi * zr[mp + j * zr_dim1];
3504/* L100: */
3505            }
3506
3507L110:
3508            ;
3509        }
3510
3511        i = int_[mp];
3512        if (i == mp) {
3513            goto L140;
3514        }
3515
3516        i_2 = *m;
3517        for (j = 1; j <= i_2; ++j) {
3518            xr = zr[i + j * zr_dim1];
3519            zr[i + j * zr_dim1] = zr[mp + j * zr_dim1];
3520            zr[mp + j * zr_dim1] = xr;
3521            xi = zi[i + j * zi_dim1];
3522            zi[i + j * zi_dim1] = zi[mp + j * zi_dim1];
3523            zi[mp + j * zi_dim1] = xi;
3524/* L130: */
3525        }
3526
3527L140:
3528        ;
3529    }
3530
3531L200:
3532    return 0;
3533} /* combak_ */
3534
3535/* Subroutine */ int comhes_(integer *nm, integer *n, integer *low, integer *
3536        igh, doublereal *ar, doublereal *ai, integer *int_)
3537{
3538    /* System generated locals */
3539    integer ar_dim1, ar_offset, ai_dim1, ai_offset, i_1, i_2, i_3;
3540    doublereal d_1, d_2;
3541
3542    /* Local variables */
3543    extern /* Subroutine */ int cdiv_(doublereal *, doublereal *, doublereal *
3544            , doublereal *, doublereal *, doublereal *);
3545    static integer i, j, m, la;
3546    static doublereal xi, yi, xr, yr;
3547    static integer mm1, kp1, mp1;
3548
3549
3550
3551/*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE COMHES, */
3552/*     NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON. */
3553/*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). */
3554
3555/*     GIVEN A COMPLEX GENERAL MATRIX, THIS SUBROUTINE */
3556/*     REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS */
3557/*     LOW THROUGH IGH TO UPPER HESSENBERG FORM BY */
3558/*     STABILIZED ELEMENTARY SIMILARITY TRANSFORMATIONS. */
3559
3560/*     ON INPUT */
3561
3562/*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
3563/*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
3564/*          DIMENSION STATEMENT. */
3565
3566/*        N IS THE ORDER OF THE MATRIX. */
3567
3568/*        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
3569/*          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED, */
3570/*          SET LOW=1, IGH=N. */
3571
3572/*        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */
3573/*          RESPECTIVELY, OF THE COMPLEX INPUT MATRIX. */
3574
3575/*     ON OUTPUT */
3576
3577/*        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */
3578/*          RESPECTIVELY, OF THE HESSENBERG MATRIX.  THE */
3579/*          MULTIPLIERS WHICH WERE USED IN THE REDUCTION */
3580/*          ARE STORED IN THE REMAINING TRIANGLES UNDER THE */
3581/*          HESSENBERG MATRIX. */
3582
3583/*        INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS */
3584/*          INTERCHANGED IN THE REDUCTION. */
3585/*          ONLY ELEMENTS LOW THROUGH IGH ARE USED. */
3586
3587/*     CALLS CDIV FOR COMPLEX DIVISION. */
3588
3589/*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
3590/*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
3591*/
3592
3593/*     THIS VERSION DATED AUGUST 1983. */
3594
3595/*     ------------------------------------------------------------------
3596*/
3597
3598    /* Parameter adjustments */
3599    ai_dim1 = *nm;
3600    ai_offset = ai_dim1 + 1;
3601    ai -= ai_offset;
3602    ar_dim1 = *nm;
3603    ar_offset = ar_dim1 + 1;
3604    ar -= ar_offset;
3605    --int_;
3606
3607    /* Function Body */
3608    la = *igh - 1;
3609    kp1 = *low + 1;
3610    if (la < kp1) {
3611        goto L200;
3612    }
3613
3614    i_1 = la;
3615    for (m = kp1; m <= i_1; ++m) {
3616        mm1 = m - 1;
3617        xr = 0.;
3618        xi = 0.;
3619        i = m;
3620
3621        i_2 = *igh;
3622        for (j = m; j <= i_2; ++j) {
3623            if ((d_1 = ar[j + mm1 * ar_dim1], abs(d_1)) + (d_2 = ai[j + 
3624                    mm1 * ai_dim1], abs(d_2)) <= abs(xr) + abs(xi)) {
3625                goto L100;
3626            }
3627            xr = ar[j + mm1 * ar_dim1];
3628            xi = ai[j + mm1 * ai_dim1];
3629            i = j;
3630L100:
3631            ;
3632        }
3633
3634        int_[m] = i;
3635        if (i == m) {
3636            goto L130;
3637        }
3638/*     .......... INTERCHANGE ROWS AND COLUMNS OF AR AND AI ..........
3639 */
3640        i_2 = *n;
3641        for (j = mm1; j <= i_2; ++j) {
3642            yr = ar[i + j * ar_dim1];
3643            ar[i + j * ar_dim1] = ar[m + j * ar_dim1];
3644            ar[m + j * ar_dim1] = yr;
3645            yi = ai[i + j * ai_dim1];
3646            ai[i + j * ai_dim1] = ai[m + j * ai_dim1];
3647            ai[m + j * ai_dim1] = yi;
3648/* L110: */
3649        }
3650
3651        i_2 = *igh;
3652        for (j = 1; j <= i_2; ++j) {
3653            yr = ar[j + i * ar_dim1];
3654            ar[j + i * ar_dim1] = ar[j + m * ar_dim1];
3655            ar[j + m * ar_dim1] = yr;
3656            yi = ai[j + i * ai_dim1];
3657            ai[j + i * ai_dim1] = ai[j + m * ai_dim1];
3658            ai[j + m * ai_dim1] = yi;
3659/* L120: */
3660        }
3661/*     .......... END INTERCHANGE .......... */
3662L130:
3663        if (xr == 0. && xi == 0.) {
3664            goto L180;
3665        }
3666        mp1 = m + 1;
3667
3668        i_2 = *igh;
3669        for (i = mp1; i <= i_2; ++i) {
3670            yr = ar[i + mm1 * ar_dim1];
3671            yi = ai[i + mm1 * ai_dim1];
3672            if (yr == 0. && yi == 0.) {
3673                goto L160;
3674            }
3675            cdiv_(&yr, &yi, &xr, &xi, &yr, &yi);
3676            ar[i + mm1 * ar_dim1] = yr;
3677            ai[i + mm1 * ai_dim1] = yi;
3678
3679            i_3 = *n;
3680            for (j = m; j <= i_3; ++j) {
3681                ar[i + j * ar_dim1] = ar[i + j * ar_dim1] - yr * ar[m + j * 
3682                        ar_dim1] + yi * ai[m + j * ai_dim1];
3683                ai[i + j * ai_dim1] = ai[i + j * ai_dim1] - yr * ai[m + j * 
3684                        ai_dim1] - yi * ar[m + j * ar_dim1];
3685/* L140: */
3686            }
3687
3688            i_3 = *igh;
3689            for (j = 1; j <= i_3; ++j) {
3690                ar[j + m * ar_dim1] = ar[j + m * ar_dim1] + yr * ar[j + i * 
3691                        ar_dim1] - yi * ai[j + i * ai_dim1];
3692                ai[j + m * ai_dim1] = ai[j + m * ai_dim1] + yr * ai[j + i * 
3693                        ai_dim1] + yi * ar[j + i * ar_dim1];
3694/* L150: */
3695            }
3696
3697L160:
3698            ;
3699        }
3700
3701L180:
3702        ;
3703    }
3704
3705L200:
3706    return 0;
3707} /* comhes_ */
3708
3709/* Subroutine */ int comlr_(integer *nm, integer *n, integer *low, integer *
3710        igh, doublereal *hr, doublereal *hi, doublereal *wr, doublereal *wi, 
3711        integer *ierr)
3712{
3713    /* System generated locals */
3714    integer hr_dim1, hr_offset, hi_dim1, hi_offset, i_1, i_2;
3715    doublereal d_1, d_2, d_3, d_4;
3716
3717    /* Local variables */
3718    extern /* Subroutine */ int cdiv_(doublereal *, doublereal *, doublereal *
3719            , doublereal *, doublereal *, doublereal *);
3720    static integer i, j, l, m, en, ll, mm;
3721    static doublereal si, ti, xi, yi, sr, tr, xr, yr;
3722    static integer im1;
3723    extern /* Subroutine */ int csroot_(doublereal *, doublereal *, 
3724            doublereal *, doublereal *);
3725    static integer mp1, itn, its;
3726    static doublereal zzi, zzr;
3727    static integer enm1;
3728    static doublereal tst1, tst2;
3729
3730
3731
3732/*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE COMLR, */
3733/*     NUM. MATH. 12, 369-376(1968) BY MARTIN AND WILKINSON. */
3734/*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 396-403(1971). */
3735
3736/*     THIS SUBROUTINE FINDS THE EIGENVALUES OF A COMPLEX */
3737/*     UPPER HESSENBERG MATRIX BY THE MODIFIED LR METHOD. */
3738
3739/*     ON INPUT */
3740
3741/*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
3742/*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
3743/*          DIMENSION STATEMENT. */
3744
3745/*        N IS THE ORDER OF THE MATRIX. */
3746
3747/*        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
3748/*          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED, */
3749/*          SET LOW=1, IGH=N. */
3750
3751/*        HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS, */
3752/*          RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX. */
3753/*          THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN THE */
3754/*          MULTIPLIERS WHICH WERE USED IN THE REDUCTION BY  COMHES, */
3755/*          IF PERFORMED. */
3756
3757/*     ON OUTPUT */
3758
3759/*        THE UPPER HESSENBERG PORTIONS OF HR AND HI HAVE BEEN */
3760/*          DESTROYED.  THEREFORE, THEY MUST BE SAVED BEFORE */
3761/*          CALLING  COMLR  IF SUBSEQUENT CALCULATION OF */
3762/*          EIGENVECTORS IS TO BE PERFORMED. */
3763
3764/*        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, */
3765/*          RESPECTIVELY, OF THE EIGENVALUES.  IF AN ERROR */
3766/*          EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT */
3767/*          FOR INDICES IERR+1,...,N. */
3768
3769/*        IERR IS SET TO */
3770/*          ZERO       FOR NORMAL RETURN, */
3771/*          J          IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED */
3772/*                     WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. */
3773
3774/*     CALLS CDIV FOR COMPLEX DIVISION. */
3775/*     CALLS CSROOT FOR COMPLEX SQUARE ROOT. */
3776
3777/*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
3778/*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
3779*/
3780
3781/*     THIS VERSION DATED AUGUST 1983. */
3782
3783/*     ------------------------------------------------------------------
3784*/
3785
3786    /* Parameter adjustments */
3787    --wi;
3788    --wr;
3789    hi_dim1 = *nm;
3790    hi_offset = hi_dim1 + 1;
3791    hi -= hi_offset;
3792    hr_dim1 = *nm;
3793    hr_offset = hr_dim1 + 1;
3794    hr -= hr_offset;
3795
3796    /* Function Body */
3797    *ierr = 0;
3798/*     .......... STORE ROOTS ISOLATED BY CBAL .......... */
3799    i_1 = *n;
3800    for (i = 1; i <= i_1; ++i) {
3801        if (i >= *low && i <= *igh) {
3802            goto L200;
3803        }
3804        wr[i] = hr[i + i * hr_dim1];
3805        wi[i] = hi[i + i * hi_dim1];
3806L200:
3807        ;
3808    }
3809
3810    en = *igh;
3811    tr = 0.;
3812    ti = 0.;
3813    itn = *n * 30;
3814/*     .......... SEARCH FOR NEXT EIGENVALUE .......... */
3815L220:
3816    if (en < *low) {
3817        goto L1001;
3818    }
3819    its = 0;
3820    enm1 = en - 1;
3821/*     .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT */
3822/*                FOR L=EN STEP -1 UNTIL LOW D0 -- .......... */
3823L240:
3824    i_1 = en;
3825    for (ll = *low; ll <= i_1; ++ll) {
3826        l = en + *low - ll;
3827        if (l == *low) {
3828            goto L300;
3829        }
3830        tst1 = (d_1 = hr[l - 1 + (l - 1) * hr_dim1], abs(d_1)) + (d_2 = hi[
3831                l - 1 + (l - 1) * hi_dim1], abs(d_2)) + (d_3 = hr[l + l * 
3832                hr_dim1], abs(d_3)) + (d_4 = hi[l + l * hi_dim1], abs(d_4))
3833                ;
3834        tst2 = tst1 + (d_1 = hr[l + (l - 1) * hr_dim1], abs(d_1)) + (d_2 = 
3835                hi[l + (l - 1) * hi_dim1], abs(d_2));
3836        if (tst2 == tst1) {
3837            goto L300;
3838        }
3839/* L260: */
3840    }
3841/*     .......... FORM SHIFT .......... */
3842L300:
3843    if (l == en) {
3844        goto L660;
3845    }
3846    if (itn == 0) {
3847        goto L1000;
3848    }
3849    if (its == 10 || its == 20) {
3850        goto L320;
3851    }
3852    sr = hr[en + en * hr_dim1];
3853    si = hi[en + en * hi_dim1];
3854    xr = hr[enm1 + en * hr_dim1] * hr[en + enm1 * hr_dim1] - hi[enm1 + en * 
3855            hi_dim1] * hi[en + enm1 * hi_dim1];
3856    xi = hr[enm1 + en * hr_dim1] * hi[en + enm1 * hi_dim1] + hi[enm1 + en * 
3857            hi_dim1] * hr[en + enm1 * hr_dim1];
3858    if (xr == 0. && xi == 0.) {
3859        goto L340;
3860    }
3861    yr = (hr[enm1 + enm1 * hr_dim1] - sr) / 2.;
3862    yi = (hi[enm1 + enm1 * hi_dim1] - si) / 2.;
3863/* Computing 2nd power */
3864    d_2 = yr;
3865/* Computing 2nd power */
3866    d_3 = yi;
3867    d_1 = d_2 * d_2 - d_3 * d_3 + xr;
3868    d_4 = yr * 2. * yi + xi;
3869    csroot_(&d_1, &d_4, &zzr, &zzi);
3870    if (yr * zzr + yi * zzi >= 0.) {
3871        goto L310;
3872    }
3873    zzr = -zzr;
3874    zzi = -zzi;
3875L310:
3876    d_1 = yr + zzr;
3877    d_2 = yi + zzi;
3878    cdiv_(&xr, &xi, &d_1, &d_2, &xr, &xi);
3879    sr -= xr;
3880    si -= xi;
3881    goto L340;
3882/*     .......... FORM EXCEPTIONAL SHIFT .......... */
3883L320:
3884    sr = (d_1 = hr[en + enm1 * hr_dim1], abs(d_1)) + (d_2 = hr[enm1 + (en
3885            - 2) * hr_dim1], abs(d_2));
3886    si = (d_1 = hi[en + enm1 * hi_dim1], abs(d_1)) + (d_2 = hi[enm1 + (en
3887            - 2) * hi_dim1], abs(d_2));
3888
3889L340:
3890    i_1 = en;
3891    for (i = *low; i <= i_1; ++i) {
3892        hr[i + i * hr_dim1] -= sr;
3893        hi[i + i * hi_dim1] -= si;
3894/* L360: */
3895    }
3896
3897    tr += sr;
3898    ti += si;
3899    ++its;
3900    --itn;
3901/*     .......... LOOK FOR TWO CONSECUTIVE SMALL */
3902/*                SUB-DIAGONAL ELEMENTS .......... */
3903    xr = (d_1 = hr[enm1 + enm1 * hr_dim1], abs(d_1)) + (d_2 = hi[enm1 + 
3904            enm1 * hi_dim1], abs(d_2));
3905    yr = (d_1 = hr[en + enm1 * hr_dim1], abs(d_1)) + (d_2 = hi[en + enm1 * 
3906            hi_dim1], abs(d_2));
3907    zzr = (d_1 = hr[en + en * hr_dim1], abs(d_1)) + (d_2 = hi[en + en * 
3908            hi_dim1], abs(d_2));
3909/*     .......... FOR M=EN-1 STEP -1 UNTIL L DO -- .......... */
3910    i_1 = enm1;
3911    for (mm = l; mm <= i_1; ++mm) {
3912        m = enm1 + l - mm;
3913        if (m == l) {
3914            goto L420;
3915        }
3916        yi = yr;
3917        yr = (d_1 = hr[m + (m - 1) * hr_dim1], abs(d_1)) + (d_2 = hi[m + (
3918                m - 1) * hi_dim1], abs(d_2));
3919        xi = zzr;
3920        zzr = xr;
3921        xr = (d_1 = hr[m - 1 + (m - 1) * hr_dim1], abs(d_1)) + (d_2 = hi[m
3922                - 1 + (m - 1) * hi_dim1], abs(d_2));
3923        tst1 = zzr / yi * (zzr + xr + xi);
3924        tst2 = tst1 + yr;
3925        if (tst2 == tst1) {
3926            goto L420;
3927        }
3928/* L380: */
3929    }
3930/*     .......... TRIANGULAR DECOMPOSITION H=L*R .......... */
3931L420:
3932    mp1 = m + 1;
3933
3934    i_1 = en;
3935    for (i = mp1; i <= i_1; ++i) {
3936        im1 = i - 1;
3937        xr = hr[im1 + im1 * hr_dim1];
3938        xi = hi[im1 + im1 * hi_dim1];
3939        yr = hr[i + im1 * hr_dim1];
3940        yi = hi[i + im1 * hi_dim1];
3941        if (abs(xr) + abs(xi) >= abs(yr) + abs(yi)) {
3942            goto L460;
3943        }
3944/*     .......... INTERCHANGE ROWS OF HR AND HI .......... */
3945        i_2 = en;
3946        for (j = im1; j <= i_2; ++j) {
3947            zzr = hr[im1 + j * hr_dim1];
3948            hr[im1 + j * hr_dim1] = hr[i + j * hr_dim1];
3949            hr[i + j * hr_dim1] = zzr;
3950            zzi = hi[im1 + j * hi_dim1];
3951            hi[im1 + j * hi_dim1] = hi[i + j * hi_dim1];
3952            hi[i + j * hi_dim1] = zzi;
3953/* L440: */
3954        }
3955
3956        cdiv_(&xr, &xi, &yr, &yi, &zzr, &zzi);
3957        wr[i] = 1.;
3958        goto L480;
3959L460:
3960        cdiv_(&yr, &yi, &xr, &xi, &zzr, &zzi);
3961        wr[i] = -1.;
3962L480:
3963        hr[i + im1 * hr_dim1] = zzr;
3964        hi[i + im1 * hi_dim1] = zzi;
3965
3966        i_2 = en;
3967        for (j = i; j <= i_2; ++j) {
3968            hr[i + j * hr_dim1] = hr[i + j * hr_dim1] - zzr * hr[im1 + j * 
3969                    hr_dim1] + zzi * hi[im1 + j * hi_dim1];
3970            hi[i + j * hi_dim1] = hi[i + j * hi_dim1] - zzr * hi[im1 + j * 
3971                    hi_dim1] - zzi * hr[im1 + j * hr_dim1];
3972/* L500: */
3973        }
3974
3975/* L520: */
3976    }
3977/*     .......... COMPOSITION R*L=H .......... */
3978    i_1 = en;
3979    for (j = mp1; j <= i_1; ++j) {
3980        xr = hr[j + (j - 1) * hr_dim1];
3981        xi = hi[j + (j - 1) * hi_dim1];
3982        hr[j + (j - 1) * hr_dim1] = 0.;
3983        hi[j + (j - 1) * hi_dim1] = 0.;
3984/*     .......... INTERCHANGE COLUMNS OF HR AND HI, */
3985/*                IF NECESSARY .......... */
3986        if (wr[j] <= 0.) {
3987            goto L580;
3988        }
3989
3990        i_2 = j;
3991        for (i = l; i <= i_2; ++i) {
3992            zzr = hr[i + (j - 1) * hr_dim1];
3993            hr[i + (j - 1) * hr_dim1] = hr[i + j * hr_dim1];
3994            hr[i + j * hr_dim1] = zzr;
3995            zzi = hi[i + (j - 1) * hi_dim1];
3996            hi[i + (j - 1) * hi_dim1] = hi[i + j * hi_dim1];
3997            hi[i + j * hi_dim1] = zzi;
3998/* L540: */
3999        }
4000
4001L580:
4002        i_2 = j;
4003        for (i = l; i <= i_2; ++i) {
4004            hr[i + (j - 1) * hr_dim1] = hr[i + (j - 1) * hr_dim1] + xr * hr[i
4005                    + j * hr_dim1] - xi * hi[i + j * hi_dim1];
4006            hi[i + (j - 1) * hi_dim1] = hi[i + (j - 1) * hi_dim1] + xr * hi[i
4007                    + j * hi_dim1] + xi * hr[i + j * hr_dim1];
4008/* L600: */
4009        }
4010
4011/* L640: */
4012    }
4013
4014    goto L240;
4015/*     .......... A ROOT FOUND .......... */
4016L660:
4017    wr[en] = hr[en + en * hr_dim1] + tr;
4018    wi[en] = hi[en + en * hi_dim1] + ti;
4019    en = enm1;
4020    goto L220;
4021/*     .......... SET ERROR -- ALL EIGENVALUES HAVE NOT */
4022/*                CONVERGED AFTER 30*N ITERATIONS .......... */
4023L1000:
4024    *ierr = en;
4025L1001:
4026    return 0;
4027} /* comlr_ */
4028
4029/* Subroutine */ int comlr2_(integer *nm, integer *n, integer *low, integer *
4030        igh, integer *int_, doublereal *hr, doublereal *hi, doublereal *wr, 
4031        doublereal *wi, doublereal *zr, doublereal *zi, integer *ierr)
4032{
4033    /* System generated locals */
4034    integer hr_dim1, hr_offset, hi_dim1, hi_offset, zr_dim1, zr_offset, 
4035            zi_dim1, zi_offset, i_1, i_2, i_3;
4036    doublereal d_1, d_2, d_3, d_4;
4037
4038    /* Local variables */
4039    static integer iend;
4040    extern /* Subroutine */ int cdiv_(doublereal *, doublereal *, doublereal *
4041            , doublereal *, doublereal *, doublereal *);
4042    static doublereal norm;
4043    static integer i, j, k, l, m, ii, en, jj, ll, mm, nn;
4044    static doublereal si, ti, xi, yi, sr, tr, xr, yr;
4045    static integer im1;
4046    extern /* Subroutine */ int csroot_(doublereal *, doublereal *, 
4047            doublereal *, doublereal *);
4048    static integer ip1, mp1, itn, its;
4049    static doublereal zzi, zzr;
4050    static integer enm1;
4051    static doublereal tst1, tst2;
4052
4053
4054
4055/*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE COMLR2, */
4056/*     NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON. */
4057/*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971). */
4058
4059/*     THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS */
4060/*     OF A COMPLEX UPPER HESSENBERG MATRIX BY THE MODIFIED LR */
4061/*     METHOD.  THE EIGENVECTORS OF A COMPLEX GENERAL MATRIX */
4062/*     CAN ALSO BE FOUND IF  COMHES  HAS BEEN USED TO REDUCE */
4063/*     THIS GENERAL MATRIX TO HESSENBERG FORM. */
4064
4065/*     ON INPUT */
4066
4067/*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
4068/*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
4069/*          DIMENSION STATEMENT. */
4070
4071/*        N IS THE ORDER OF THE MATRIX. */
4072
4073/*        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
4074/*          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED, */
4075/*          SET LOW=1, IGH=N. */
4076
4077/*        INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS INTERCHANGED */
4078/*          IN THE REDUCTION BY  COMHES, IF PERFORMED.  ONLY ELEMENTS */
4079/*          LOW THROUGH IGH ARE USED.  IF THE EIGENVECTORS OF THE HESSEN-
4080*/
4081/*          BERG MATRIX ARE DESIRED, SET INT(J)=J FOR THESE ELEMENTS. */
4082
4083/*        HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS, */
4084/*          RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX. */
4085/*          THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN THE */
4086/*          MULTIPLIERS WHICH WERE USED IN THE REDUCTION BY  COMHES, */
4087/*          IF PERFORMED.  IF THE EIGENVECTORS OF THE HESSENBERG */
4088/*          MATRIX ARE DESIRED, THESE ELEMENTS MUST BE SET TO ZERO. */
4089
4090/*     ON OUTPUT */
4091
4092/*        THE UPPER HESSENBERG PORTIONS OF HR AND HI HAVE BEEN */
4093/*          DESTROYED, BUT THE LOCATION HR(1,1) CONTAINS THE NORM */
4094/*          OF THE TRIANGULARIZED MATRIX. */
4095
4096/*        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, */
4097/*          RESPECTIVELY, OF THE EIGENVALUES.  IF AN ERROR */
4098/*          EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT */
4099/*          FOR INDICES IERR+1,...,N. */
4100
4101/*        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, */
4102/*          RESPECTIVELY, OF THE EIGENVECTORS.  THE EIGENVECTORS */
4103/*          ARE UNNORMALIZED.  IF AN ERROR EXIT IS MADE, NONE OF */
4104/*          THE EIGENVECTORS HAS BEEN FOUND. */
4105
4106/*        IERR IS SET TO */
4107/*          ZERO       FOR NORMAL RETURN, */
4108/*          J          IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED */
4109/*                     WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. */
4110
4111
4112/*     CALLS CDIV FOR COMPLEX DIVISION. */
4113/*     CALLS CSROOT FOR COMPLEX SQUARE ROOT. */
4114
4115/*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
4116/*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
4117*/
4118
4119/*     THIS VERSION DATED AUGUST 1983. */
4120
4121/*     ------------------------------------------------------------------
4122*/
4123
4124    /* Parameter adjustments */
4125    zi_dim1 = *nm;
4126    zi_offset = zi_dim1 + 1;
4127    zi -= zi_offset;
4128    zr_dim1 = *nm;
4129    zr_offset = zr_dim1 + 1;
4130    zr -= zr_offset;
4131    --wi;
4132    --wr;
4133    hi_dim1 = *nm;
4134    hi_offset = hi_dim1 + 1;
4135    hi -= hi_offset;
4136    hr_dim1 = *nm;
4137    hr_offset = hr_dim1 + 1;
4138    hr -= hr_offset;
4139    --int_;
4140
4141    /* Function Body */
4142    *ierr = 0;
4143/*     .......... INITIALIZE EIGENVECTOR MATRIX .......... */
4144    i_1 = *n;
4145    for (i = 1; i <= i_1; ++i) {
4146
4147        i_2 = *n;
4148        for (j = 1; j <= i_2; ++j) {
4149            zr[i + j * zr_dim1] = 0.;
4150            zi[i + j * zi_dim1] = 0.;
4151            if (i == j) {
4152                zr[i + j * zr_dim1] = 1.;
4153            }
4154/* L100: */
4155        }
4156    }
4157/*     .......... FORM THE MATRIX OF ACCUMULATED TRANSFORMATIONS */
4158/*                FROM THE INFORMATION LEFT BY COMHES .......... */
4159    iend = *igh - *low - 1;
4160    if (iend <= 0) {
4161        goto L180;
4162    }
4163/*     .......... FOR I=IGH-1 STEP -1 UNTIL LOW+1 DO -- .......... */
4164    i_2 = iend;
4165    for (ii = 1; ii <= i_2; ++ii) {
4166        i = *igh - ii;
4167        ip1 = i + 1;
4168
4169        i_1 = *igh;
4170        for (k = ip1; k <= i_1; ++k) {
4171            zr[k + i * zr_dim1] = hr[k + (i - 1) * hr_dim1];
4172            zi[k + i * zi_dim1] = hi[k + (i - 1) * hi_dim1];
4173/* L120: */
4174        }
4175
4176        j = int_[i];
4177        if (i == j) {
4178            goto L160;
4179        }
4180
4181        i_1 = *igh;
4182        for (k = i; k <= i_1; ++k) {
4183            zr[i + k * zr_dim1] = zr[j + k * zr_dim1];
4184            zi[i + k * zi_dim1] = zi[j + k * zi_dim1];
4185            zr[j + k * zr_dim1] = 0.;
4186            zi[j + k * zi_dim1] = 0.;
4187/* L140: */
4188        }
4189
4190        zr[j + i * zr_dim1] = 1.;
4191L160:
4192        ;
4193    }
4194/*     .......... STORE ROOTS ISOLATED BY CBAL .......... */
4195L180:
4196    i_2 = *n;
4197    for (i = 1; i <= i_2; ++i) {
4198        if (i >= *low && i <= *igh) {
4199            goto L200;
4200        }
4201        wr[i] = hr[i + i * hr_dim1];
4202        wi[i] = hi[i + i * hi_dim1];
4203L200:
4204        ;
4205    }
4206
4207    en = *igh;
4208    tr = 0.;
4209    ti = 0.;
4210    itn = *n * 30;
4211/*     .......... SEARCH FOR NEXT EIGENVALUE .......... */
4212L220:
4213    if (en < *low) {
4214        goto L680;
4215    }
4216    its = 0;
4217    enm1 = en - 1;
4218/*     .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT */
4219/*                FOR L=EN STEP -1 UNTIL LOW DO -- .......... */
4220L240:
4221    i_2 = en;
4222    for (ll = *low; ll <= i_2; ++ll) {
4223        l = en + *low - ll;
4224        if (l == *low) {
4225            goto L300;
4226        }
4227        tst1 = (d_1 = hr[l - 1 + (l - 1) * hr_dim1], abs(d_1)) + (d_2 = hi[
4228                l - 1 + (l - 1) * hi_dim1], abs(d_2)) + (d_3 = hr[l + l * 
4229                hr_dim1], abs(d_3)) + (d_4 = hi[l + l * hi_dim1], abs(d_4))
4230                ;
4231        tst2 = tst1 + (d_1 = hr[l + (l - 1) * hr_dim1], abs(d_1)) + (d_2 = 
4232                hi[l + (l - 1) * hi_dim1], abs(d_2));
4233        if (tst2 == tst1) {
4234            goto L300;
4235        }
4236/* L260: */
4237    }
4238/*     .......... FORM SHIFT .......... */
4239L300:
4240    if (l == en) {
4241        goto L660;
4242    }
4243    if (itn == 0) {
4244        goto L1000;
4245    }
4246    if (its == 10 || its == 20) {
4247        goto L320;
4248    }
4249    sr = hr[en + en * hr_dim1];
4250    si = hi[en + en * hi_dim1];
4251    xr = hr[enm1 + en * hr_dim1] * hr[en + enm1 * hr_dim1] - hi[enm1 + en * 
4252            hi_dim1] * hi[en + enm1 * hi_dim1];
4253    xi = hr[enm1 + en * hr_dim1] * hi[en + enm1 * hi_dim1] + hi[enm1 + en * 
4254            hi_dim1] * hr[en + enm1 * hr_dim1];
4255    if (xr == 0. && xi == 0.) {
4256        goto L340;
4257    }
4258    yr = (hr[enm1 + enm1 * hr_dim1] - sr) / 2.;
4259    yi = (hi[enm1 + enm1 * hi_dim1] - si) / 2.;
4260/* Computing 2nd power */
4261    d_2 = yr;
4262/* Computing 2nd power */
4263    d_3 = yi;
4264    d_1 = d_2 * d_2 - d_3 * d_3 + xr;
4265    d_4 = yr * 2. * yi + xi;
4266    csroot_(&d_1, &d_4, &zzr, &zzi);
4267    if (yr * zzr + yi * zzi >= 0.) {
4268        goto L310;
4269    }
4270    zzr = -zzr;
4271    zzi = -zzi;
4272L310:
4273    d_1 = yr + zzr;
4274    d_2 = yi + zzi;
4275    cdiv_(&xr, &xi, &d_1, &d_2, &xr, &xi);
4276    sr -= xr;
4277    si -= xi;
4278    goto L340;
4279/*     .......... FORM EXCEPTIONAL SHIFT .......... */
4280L320:
4281    sr = (d_1 = hr[en + enm1 * hr_dim1], abs(d_1)) + (d_2 = hr[enm1 + (en
4282            - 2) * hr_dim1], abs(d_2));
4283    si = (d_1 = hi[en + enm1 * hi_dim1], abs(d_1)) + (d_2 = hi[enm1 + (en
4284            - 2) * hi_dim1], abs(d_2));
4285
4286L340:
4287    i_2 = en;
4288    for (i = *low; i <= i_2; ++i) {
4289        hr[i + i * hr_dim1] -= sr;
4290        hi[i + i * hi_dim1] -= si;
4291/* L360: */
4292    }
4293
4294    tr += sr;
4295    ti += si;
4296    ++its;
4297    --itn;
4298/*     .......... LOOK FOR TWO CONSECUTIVE SMALL */
4299/*                SUB-DIAGONAL ELEMENTS .......... */
4300    xr = (d_1 = hr[enm1 + enm1 * hr_dim1], abs(d_1)) + (d_2 = hi[enm1 + 
4301            enm1 * hi_dim1], abs(d_2));
4302    yr = (d_1 = hr[en + enm1 * hr_dim1], abs(d_1)) + (d_2 = hi[en + enm1 * 
4303            hi_dim1], abs(d_2));
4304    zzr = (d_1 = hr[en + en * hr_dim1], abs(d_1)) + (d_2 = hi[en + en * 
4305            hi_dim1], abs(d_2));
4306/*     .......... FOR M=EN-1 STEP -1 UNTIL L DO -- .......... */
4307    i_2 = enm1;
4308    for (mm = l; mm <= i_2; ++mm) {
4309        m = enm1 + l - mm;
4310        if (m == l) {
4311            goto L420;
4312        }
4313        yi = yr;
4314        yr = (d_1 = hr[m + (m - 1) * hr_dim1], abs(d_1)) + (d_2 = hi[m + (
4315                m - 1) * hi_dim1], abs(d_2));
4316        xi = zzr;
4317        zzr = xr;
4318        xr = (d_1 = hr[m - 1 + (m - 1) * hr_dim1], abs(d_1)) + (d_2 = hi[m
4319                - 1 + (m - 1) * hi_dim1], abs(d_2));
4320        tst1 = zzr / yi * (zzr + xr + xi);
4321        tst2 = tst1 + yr;
4322        if (tst2 == tst1) {
4323            goto L420;
4324        }
4325/* L380: */
4326    }
4327/*     .......... TRIANGULAR DECOMPOSITION H=L*R .......... */
4328L420:
4329    mp1 = m + 1;
4330
4331    i_2 = en;
4332    for (i = mp1; i <= i_2; ++i) {
4333        im1 = i - 1;
4334        xr = hr[im1 + im1 * hr_dim1];
4335        xi = hi[im1 + im1 * hi_dim1];
4336        yr = hr[i + im1 * hr_dim1];
4337        yi = hi[i + im1 * hi_dim1];
4338        if (abs(xr) + abs(xi) >= abs(yr) + abs(yi)) {
4339            goto L460;
4340        }
4341/*     .......... INTERCHANGE ROWS OF HR AND HI .......... */
4342        i_1 = *n;
4343        for (j = im1; j <= i_1; ++j) {
4344            zzr = hr[im1 + j * hr_dim1];
4345            hr[im1 + j * hr_dim1] = hr[i + j * hr_dim1];
4346            hr[i + j * hr_dim1] = zzr;
4347            zzi = hi[im1 + j * hi_dim1];
4348            hi[im1 + j * hi_dim1] = hi[i + j * hi_dim1];
4349            hi[i + j * hi_dim1] = zzi;
4350/* L440: */
4351        }
4352
4353        cdiv_(&xr, &xi, &yr, &yi, &zzr, &zzi);
4354        wr[i] = 1.;
4355        goto L480;
4356L460:
4357        cdiv_(&yr, &yi, &xr, &xi, &zzr, &zzi);
4358        wr[i] = -1.;
4359L480:
4360        hr[i + im1 * hr_dim1] = zzr;
4361        hi[i + im1 * hi_dim1] = zzi;
4362
4363        i_1 = *n;
4364        for (j = i; j <= i_1; ++j) {
4365            hr[i + j * hr_dim1] = hr[i + j * hr_dim1] - zzr * hr[im1 + j * 
4366                    hr_dim1] + zzi * hi[im1 + j * hi_dim1];
4367            hi[i + j * hi_dim1] = hi[i + j * hi_dim1] - zzr * hi[im1 + j * 
4368                    hi_dim1] - zzi * hr[im1 + j * hr_dim1];
4369/* L500: */
4370        }
4371
4372/* L520: */
4373    }
4374/*     .......... COMPOSITION R*L=H .......... */
4375    i_2 = en;
4376    for (j = mp1; j <= i_2; ++j) {
4377        xr = hr[j + (j - 1) * hr_dim1];
4378        xi = hi[j + (j - 1) * hi_dim1];
4379        hr[j + (j - 1) * hr_dim1] = 0.;
4380        hi[j + (j - 1) * hi_dim1] = 0.;
4381/*     .......... INTERCHANGE COLUMNS OF HR, HI, ZR, AND ZI, */
4382/*                IF NECESSARY .......... */
4383        if (wr[j] <= 0.) {
4384            goto L580;
4385        }
4386
4387        i_1 = j;
4388        for (i = 1; i <= i_1; ++i) {
4389            zzr = hr[i + (j - 1) * hr_dim1];
4390            hr[i + (j - 1) * hr_dim1] = hr[i + j * hr_dim1];
4391            hr[i + j * hr_dim1] = zzr;
4392            zzi = hi[i + (j - 1) * hi_dim1];
4393            hi[i + (j - 1) * hi_dim1] = hi[i + j * hi_dim1];
4394            hi[i + j * hi_dim1] = zzi;
4395/* L540: */
4396        }
4397
4398        i_1 = *igh;
4399        for (i = *low; i <= i_1; ++i) {
4400            zzr = zr[i + (j - 1) * zr_dim1];
4401            zr[i + (j - 1) * zr_dim1] = zr[i + j * zr_dim1];
4402            zr[i + j * zr_dim1] = zzr;
4403            zzi = zi[i + (j - 1) * zi_dim1];
4404            zi[i + (j - 1) * zi_dim1] = zi[i + j * zi_dim1];
4405            zi[i + j * zi_dim1] = zzi;
4406/* L560: */
4407        }
4408
4409L580:
4410        i_1 = j;
4411        for (i = 1; i <= i_1; ++i) {
4412            hr[i + (j - 1) * hr_dim1] = hr[i + (j - 1) * hr_dim1] + xr * hr[i
4413                    + j * hr_dim1] - xi * hi[i + j * hi_dim1];
4414            hi[i + (j - 1) * hi_dim1] = hi[i + (j - 1) * hi_dim1] + xr * hi[i
4415                    + j * hi_dim1] + xi * hr[i + j * hr_dim1];
4416/* L600: */
4417        }
4418/*     .......... ACCUMULATE TRANSFORMATIONS .......... */
4419        i_1 = *igh;
4420        for (i = *low; i <= i_1; ++i) {
4421            zr[i + (j - 1) * zr_dim1] = zr[i + (j - 1) * zr_dim1] + xr * zr[i
4422                    + j * zr_dim1] - xi * zi[i + j * zi_dim1];
4423            zi[i + (j - 1) * zi_dim1] = zi[i + (j - 1) * zi_dim1] + xr * zi[i
4424                    + j * zi_dim1] + xi * zr[i + j * zr_dim1];
4425/* L620: */
4426        }
4427
4428/* L640: */
4429    }
4430
4431    goto L240;
4432/*     .......... A ROOT FOUND .......... */
4433L660:
4434    hr[en + en * hr_dim1] += tr;
4435    wr[en] = hr[en + en * hr_dim1];
4436    hi[en + en * hi_dim1] += ti;
4437    wi[en] = hi[en + en * hi_dim1];
4438    en = enm1;
4439    goto L220;
4440/*     .......... ALL ROOTS FOUND.  BACKSUBSTITUTE TO FIND */
4441/*                VECTORS OF UPPER TRIANGULAR FORM .......... */
4442L680:
4443    norm = 0.;
4444
4445    i_2 = *n;
4446    for (i = 1; i <= i_2; ++i) {
4447
4448        i_1 = *n;
4449        for (j = i; j <= i_1; ++j) {
4450            tr = (d_1 = hr[i + j * hr_dim1], abs(d_1)) + (d_2 = hi[i + j * 
4451                    hi_dim1], abs(d_2));
4452            if (tr > norm) {
4453                norm = tr;
4454            }
4455/* L720: */
4456        }
4457    }
4458
4459    hr[hr_dim1 + 1] = norm;
4460    if (*n == 1 || norm == 0.) {
4461        goto L1001;
4462    }
4463/*     .......... FOR EN=N STEP -1 UNTIL 2 DO -- .......... */
4464    i_1 = *n;
4465    for (nn = 2; nn <= i_1; ++nn) {
4466        en = *n + 2 - nn;
4467        xr = wr[en];
4468        xi = wi[en];
4469        hr[en + en * hr_dim1] = 1.;
4470        hi[en + en * hi_dim1] = 0.;
4471        enm1 = en - 1;
4472/*     .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- .......... */
4473        i_2 = enm1;
4474        for (ii = 1; ii <= i_2; ++ii) {
4475            i = en - ii;
4476            zzr = 0.;
4477            zzi = 0.;
4478            ip1 = i + 1;
4479
4480            i_3 = en;
4481            for (j = ip1; j <= i_3; ++j) {
4482                zzr = zzr + hr[i + j * hr_dim1] * hr[j + en * hr_dim1] - hi[i
4483                        + j * hi_dim1] * hi[j + en * hi_dim1];
4484                zzi = zzi + hr[i + j * hr_dim1] * hi[j + en * hi_dim1] + hi[i
4485                        + j * hi_dim1] * hr[j + en * hr_dim1];
4486/* L740: */
4487            }
4488
4489            yr = xr - wr[i];
4490            yi = xi - wi[i];
4491            if (yr != 0. || yi != 0.) {
4492                goto L765;
4493            }
4494            tst1 = norm;
4495            yr = tst1;
4496L760:
4497            yr *= .01;
4498            tst2 = norm + yr;
4499            if (tst2 > tst1) {
4500                goto L760;
4501            }
4502L765:
4503            cdiv_(&zzr, &zzi, &yr, &yi, &hr[i + en * hr_dim1], &hi[i + en * 
4504                    hi_dim1]);
4505/*     .......... OVERFLOW CONTROL .......... */
4506            tr = (d_1 = hr[i + en * hr_dim1], abs(d_1)) + (d_2 = hi[i + en
4507                    * hi_dim1], abs(d_2));
4508            if (tr == 0.) {
4509                goto L780;
4510            }
4511            tst1 = tr;
4512            tst2 = tst1 + 1. / tst1;
4513            if (tst2 > tst1) {
4514                goto L780;
4515            }
4516            i_3 = en;
4517            for (j = i; j <= i_3; ++j) {
4518                hr[j + en * hr_dim1] /= tr;
4519                hi[j + en * hi_dim1] /= tr;
4520/* L770: */
4521            }
4522
4523L780:
4524            ;
4525        }
4526
4527/* L800: */
4528    }
4529/*     .......... END BACKSUBSTITUTION .......... */
4530    enm1 = *n - 1;
4531/*     .......... VECTORS OF ISOLATED ROOTS .......... */
4532    i_1 = enm1;
4533    for (i = 1; i <= i_1; ++i) {
4534        if (i >= *low && i <= *igh) {
4535            goto L840;
4536        }
4537        ip1 = i + 1;
4538
4539        i_2 = *n;
4540        for (j = ip1; j <= i_2; ++j) {
4541            zr[i + j * zr_dim1] = hr[i + j * hr_dim1];
4542            zi[i + j * zi_dim1] = hi[i + j * hi_dim1];
4543/* L820: */
4544        }
4545
4546L840:
4547        ;
4548    }
4549/*     .......... MULTIPLY BY TRANSFORMATION MATRIX TO GIVE */
4550/*                VECTORS OF ORIGINAL FULL MATRIX. */
4551/*                FOR J=N STEP -1 UNTIL LOW+1 DO -- .......... */
4552    i_1 = enm1;
4553    for (jj = *low; jj <= i_1; ++jj) {
4554        j = *n + *low - jj;
4555        m = min(j,*igh);
4556
4557        i_2 = *igh;
4558        for (i = *low; i <= i_2; ++i) {
4559            zzr = 0.;
4560            zzi = 0.;
4561
4562            i_3 = m;
4563            for (k = *low; k <= i_3; ++k) {
4564                zzr = zzr + zr[i + k * zr_dim1] * hr[k + j * hr_dim1] - zi[i
4565                        + k * zi_dim1] * hi[k + j * hi_dim1];
4566                zzi = zzi + zr[i + k * zr_dim1] * hi[k + j * hi_dim1] + zi[i
4567                        + k * zi_dim1] * hr[k + j * hr_dim1];
4568/* L860: */
4569            }
4570
4571            zr[i + j * zr_dim1] = zzr;
4572            zi[i + j * zi_dim1] = zzi;
4573/* L880: */
4574        }
4575    }
4576
4577    goto L1001;
4578/*     .......... SET ERROR -- ALL EIGENVALUES HAVE NOT */
4579/*                CONVERGED AFTER 30*N ITERATIONS .......... */
4580L1000:
4581    *ierr = en;
4582L1001:
4583    return 0;
4584} /* comlr2_ */
4585
4586/* Subroutine */ int comqr_(integer *nm, integer *n, integer *low, integer *
4587        igh, doublereal *hr, doublereal *hi, doublereal *wr, doublereal *wi, 
4588        integer *ierr)
4589{
4590    /* System generated locals */
4591    integer hr_dim1, hr_offset, hi_dim1, hi_offset, i_1, i_2;
4592    doublereal d_1, d_2, d_3, d_4;
4593
4594    /* Local variables */
4595    extern /* Subroutine */ int cdiv_(doublereal *, doublereal *, doublereal *
4596            , doublereal *, doublereal *, doublereal *);
4597    static doublereal norm;
4598    static integer i, j, l, en, ll;
4599    static doublereal si, ti, xi, yi, sr, tr, xr, yr;
4600    extern doublereal pythag_(doublereal *, doublereal *);
4601    extern /* Subroutine */ int csroot_(doublereal *, doublereal *, 
4602            doublereal *, doublereal *);
4603    static integer lp1, itn, its;
4604    static doublereal zzi, zzr;
4605    static integer enm1;
4606    static doublereal tst1, tst2;
4607
4608
4609
4610/*     THIS SUBROUTINE IS A TRANSLATION OF A UNITARY ANALOGUE OF THE */
4611/*     ALGOL PROCEDURE  COMLR, NUM. MATH. 12, 369-376(1968) BY MARTIN */
4612/*     AND WILKINSON. */
4613/*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 396-403(1971). */
4614/*     THE UNITARY ANALOGUE SUBSTITUTES THE QR ALGORITHM OF FRANCIS */
4615/*     (COMP. JOUR. 4, 332-345(1962)) FOR THE LR ALGORITHM. */
4616
4617/*     THIS SUBROUTINE FINDS THE EIGENVALUES OF A COMPLEX */
4618/*     UPPER HESSENBERG MATRIX BY THE QR METHOD. */
4619
4620/*     ON INPUT */
4621
4622/*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
4623/*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
4624/*          DIMENSION STATEMENT. */
4625
4626/*        N IS THE ORDER OF THE MATRIX. */
4627
4628/*        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
4629/*          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED, */
4630/*          SET LOW=1, IGH=N. */
4631
4632/*        HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS, */
4633/*          RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX. */
4634/*          THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN */
4635/*          INFORMATION ABOUT THE UNITARY TRANSFORMATIONS USED IN */
4636/*          THE REDUCTION BY  CORTH, IF PERFORMED. */
4637
4638/*     ON OUTPUT */
4639
4640/*        THE UPPER HESSENBERG PORTIONS OF HR AND HI HAVE BEEN */
4641/*          DESTROYED.  THEREFORE, THEY MUST BE SAVED BEFORE */
4642/*          CALLING  COMQR  IF SUBSEQUENT CALCULATION OF */
4643/*          EIGENVECTORS IS TO BE PERFORMED. */
4644
4645/*        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, */
4646/*          RESPECTIVELY, OF THE EIGENVALUES.  IF AN ERROR */
4647/*          EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT */
4648/*          FOR INDICES IERR+1,...,N. */
4649
4650/*        IERR IS SET TO */
4651/*          ZERO       FOR NORMAL RETURN, */
4652/*          J          IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED */
4653/*                     WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. */
4654
4655/*     CALLS CDIV FOR COMPLEX DIVISION. */
4656/*     CALLS CSROOT FOR COMPLEX SQUARE ROOT. */
4657/*     CALLS PYTHAG FOR  DSQRT(A*A + B*B) . */
4658
4659/*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
4660/*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
4661*/
4662
4663/*     THIS VERSION DATED AUGUST 1983. */
4664
4665/*     ------------------------------------------------------------------
4666*/
4667
4668    /* Parameter adjustments */
4669    --wi;
4670    --wr;
4671    hi_dim1 = *nm;
4672    hi_offset = hi_dim1 + 1;
4673    hi -= hi_offset;
4674    hr_dim1 = *nm;
4675    hr_offset = hr_dim1 + 1;
4676    hr -= hr_offset;
4677
4678    /* Function Body */
4679    *ierr = 0;
4680    if (*low == *igh) {
4681        goto L180;
4682    }
4683/*     .......... CREATE REAL SUBDIAGONAL ELEMENTS .......... */
4684    l = *low + 1;
4685
4686    i_1 = *igh;
4687    for (i = l; i <= i_1; ++i) {
4688/* Computing MIN */
4689        i_2 = i + 1;
4690        ll = min(i_2,*igh);
4691        if (hi[i + (i - 1) * hi_dim1] == 0.) {
4692            goto L170;
4693        }
4694        norm = pythag_(&hr[i + (i - 1) * hr_dim1], &hi[i + (i - 1) * hi_dim1])
4695                ;
4696        yr = hr[i + (i - 1) * hr_dim1] / norm;
4697        yi = hi[i + (i - 1) * hi_dim1] / norm;
4698        hr[i + (i - 1) * hr_dim1] = norm;
4699        hi[i + (i - 1) * hi_dim1] = 0.;
4700
4701        i_2 = *igh;
4702        for (j = i; j <= i_2; ++j) {
4703            si = yr * hi[i + j * hi_dim1] - yi * hr[i + j * hr_dim1];
4704            hr[i + j * hr_dim1] = yr * hr[i + j * hr_dim1] + yi * hi[i + j * 
4705                    hi_dim1];
4706            hi[i + j * hi_dim1] = si;
4707/* L155: */
4708        }
4709
4710        i_2 = ll;
4711        for (j = *low; j <= i_2; ++j) {
4712            si = yr * hi[j + i * hi_dim1] + yi * hr[j + i * hr_dim1];
4713            hr[j + i * hr_dim1] = yr * hr[j + i * hr_dim1] - yi * hi[j + i * 
4714                    hi_dim1];
4715            hi[j + i * hi_dim1] = si;
4716/* L160: */
4717        }
4718
4719L170:
4720        ;
4721    }
4722/*     .......... STORE ROOTS ISOLATED BY CBAL .......... */
4723L180:
4724    i_1 = *n;
4725    for (i = 1; i <= i_1; ++i) {
4726        if (i >= *low && i <= *igh) {
4727            goto L200;
4728        }
4729        wr[i] = hr[i + i * hr_dim1];
4730        wi[i] = hi[i + i * hi_dim1];
4731L200:
4732        ;
4733    }
4734
4735    en = *igh;
4736    tr = 0.;
4737    ti = 0.;
4738    itn = *n * 30;
4739/*     .......... SEARCH FOR NEXT EIGENVALUE .......... */
4740L220:
4741    if (en < *low) {
4742        goto L1001;
4743    }
4744    its = 0;
4745    enm1 = en - 1;
4746/*     .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT */
4747/*                FOR L=EN STEP -1 UNTIL LOW D0 -- .......... */
4748L240:
4749    i_1 = en;
4750    for (ll = *low; ll <= i_1; ++ll) {
4751        l = en + *low - ll;
4752        if (l == *low) {
4753            goto L300;
4754        }
4755        tst1 = (d_1 = hr[l - 1 + (l - 1) * hr_dim1], abs(d_1)) + (d_2 = hi[
4756                l - 1 + (l - 1) * hi_dim1], abs(d_2)) + (d_3 = hr[l + l * 
4757                hr_dim1], abs(d_3)) + (d_4 = hi[l + l * hi_dim1], abs(d_4))
4758                ;
4759        tst2 = tst1 + (d_1 = hr[l + (l - 1) * hr_dim1], abs(d_1));
4760        if (tst2 == tst1) {
4761            goto L300;
4762        }
4763/* L260: */
4764    }
4765/*     .......... FORM SHIFT .......... */
4766L300:
4767    if (l == en) {
4768        goto L660;
4769    }
4770    if (itn == 0) {
4771        goto L1000;
4772    }
4773    if (its == 10 || its == 20) {
4774        goto L320;
4775    }
4776    sr = hr[en + en * hr_dim1];
4777    si = hi[en + en * hi_dim1];
4778    xr = hr[enm1 + en * hr_dim1] * hr[en + enm1 * hr_dim1];
4779    xi = hi[enm1 + en * hi_dim1] * hr[en + enm1 * hr_dim1];
4780    if (xr == 0. && xi == 0.) {
4781        goto L340;
4782    }
4783    yr = (hr[enm1 + enm1 * hr_dim1] - sr) / 2.;
4784    yi = (hi[enm1 + enm1 * hi_dim1] - si) / 2.;
4785/* Computing 2nd power */
4786    d_2 = yr;
4787/* Computing 2nd power */
4788    d_3 = yi;
4789    d_1 = d_2 * d_2 - d_3 * d_3 + xr;
4790    d_4 = yr * 2. * yi + xi;
4791    csroot_(&d_1, &d_4, &zzr, &zzi);
4792    if (yr * zzr + yi * zzi >= 0.) {
4793        goto L310;
4794    }
4795    zzr = -zzr;
4796    zzi = -zzi;
4797L310:
4798    d_1 = yr + zzr;
4799    d_2 = yi + zzi;
4800    cdiv_(&xr, &xi, &d_1, &d_2, &xr, &xi);
4801    sr -= xr;
4802    si -= xi;
4803    goto L340;
4804/*     .......... FORM EXCEPTIONAL SHIFT .......... */
4805L320:
4806    sr = (d_1 = hr[en + enm1 * hr_dim1], abs(d_1)) + (d_2 = hr[enm1 + (en
4807            - 2) * hr_dim1], abs(d_2));
4808    si = 0.;
4809
4810L340:
4811    i_1 = en;
4812    for (i = *low; i <= i_1; ++i) {
4813        hr[i + i * hr_dim1] -= sr;
4814        hi[i + i * hi_dim1] -= si;
4815/* L360: */
4816    }
4817
4818    tr += sr;
4819    ti += si;
4820    ++its;
4821    --itn;
4822/*     .......... REDUCE TO TRIANGLE (ROWS) .......... */
4823    lp1 = l + 1;
4824
4825    i_1 = en;
4826    for (i = lp1; i <= i_1; ++i) {
4827        sr = hr[i + (i - 1) * hr_dim1];
4828        hr[i + (i - 1) * hr_dim1] = 0.;
4829        d_1 = pythag_(&hr[i - 1 + (i - 1) * hr_dim1], &hi[i - 1 + (i - 1) * 
4830                hi_dim1]);
4831        norm = pythag_(&d_1, &sr);
4832        xr = hr[i - 1 + (i - 1) * hr_dim1] / norm;
4833        wr[i - 1] = xr;
4834        xi = hi[i - 1 + (i - 1) * hi_dim1] / norm;
4835        wi[i - 1] = xi;
4836        hr[i - 1 + (i - 1) * hr_dim1] = norm;
4837        hi[i - 1 + (i - 1) * hi_dim1] = 0.;
4838        hi[i + (i - 1) * hi_dim1] = sr / norm;
4839
4840        i_2 = en;
4841        for (j = i; j <= i_2; ++j) {
4842            yr = hr[i - 1 + j * hr_dim1];
4843            yi = hi[i - 1 + j * hi_dim1];
4844            zzr = hr[i + j * hr_dim1];
4845            zzi = hi[i + j * hi_dim1];
4846            hr[i - 1 + j * hr_dim1] = xr * yr + xi * yi + hi[i + (i - 1) * 
4847                    hi_dim1] * zzr;
4848            hi[i - 1 + j * hi_dim1] = xr * yi - xi * yr + hi[i + (i - 1) * 
4849                    hi_dim1] * zzi;
4850            hr[i + j * hr_dim1] = xr * zzr - xi * zzi - hi[i + (i - 1) * 
4851                    hi_dim1] * yr;
4852            hi[i + j * hi_dim1] = xr * zzi + xi * zzr - hi[i + (i - 1) * 
4853                    hi_dim1] * yi;
4854/* L490: */
4855        }
4856
4857/* L500: */
4858    }
4859
4860    si = hi[en + en * hi_dim1];
4861    if (si == 0.) {
4862        goto L540;
4863    }
4864    norm = pythag_(&hr[en + en * hr_dim1], &si);
4865    sr = hr[en + en * hr_dim1] / norm;
4866    si /= norm;
4867    hr[en + en * hr_dim1] = norm;
4868    hi[en + en * hi_dim1] = 0.;
4869/*     .......... INVERSE OPERATION (COLUMNS) .......... */
4870L540:
4871    i_1 = en;
4872    for (j = lp1; j <= i_1; ++j) {
4873        xr = wr[j - 1];
4874        xi = wi[j - 1];
4875
4876        i_2 = j;
4877        for (i = l; i <= i_2; ++i) {
4878            yr = hr[i + (j - 1) * hr_dim1];
4879            yi = 0.;
4880            zzr = hr[i + j * hr_dim1];
4881            zzi = hi[i + j * hi_dim1];
4882            if (i == j) {
4883                goto L560;
4884            }
4885            yi = hi[i + (j - 1) * hi_dim1];
4886            hi[i + (j - 1) * hi_dim1] = xr * yi + xi * yr + hi[j + (j - 1) * 
4887                    hi_dim1] * zzi;
4888L560:
4889            hr[i + (j - 1) * hr_dim1] = xr * yr - xi * yi + hi[j + (j - 1) * 
4890                    hi_dim1] * zzr;
4891            hr[i + j * hr_dim1] = xr * zzr + xi * zzi - hi[j + (j - 1) * 
4892                    hi_dim1] * yr;
4893            hi[i + j * hi_dim1] = xr * zzi - xi * zzr - hi[j + (j - 1) * 
4894                    hi_dim1] * yi;
4895/* L580: */
4896        }
4897
4898/* L600: */
4899    }
4900
4901    if (si == 0.) {
4902        goto L240;
4903    }
4904
4905    i_1 = en;
4906    for (i = l; i <= i_1; ++i) {
4907        yr = hr[i + en * hr_dim1];
4908        yi = hi[i + en * hi_dim1];
4909        hr[i + en * hr_dim1] = sr * yr - si * yi;
4910        hi[i + en * hi_dim1] = sr * yi + si * yr;
4911/* L630: */
4912    }
4913
4914    goto L240;
4915/*     .......... A ROOT FOUND .......... */
4916L660:
4917    wr[en] = hr[en + en * hr_dim1] + tr;
4918    wi[en] = hi[en + en * hi_dim1] + ti;
4919    en = enm1;
4920    goto L220;
4921/*     .......... SET ERROR -- ALL EIGENVALUES HAVE NOT */
4922/*                CONVERGED AFTER 30*N ITERATIONS .......... */
4923L1000:
4924    *ierr = en;
4925L1001:
4926    return 0;
4927} /* comqr_ */
4928
4929/* Subroutine */ int comqr2_(integer *nm, integer *n, integer *low, integer *
4930        igh, doublereal *ortr, doublereal *orti, doublereal *hr, doublereal *
4931        hi, doublereal *wr, doublereal *wi, doublereal *zr, doublereal *zi, 
4932        integer *ierr)
4933{
4934    /* System generated locals */
4935    integer hr_dim1, hr_offset, hi_dim1, hi_offset, zr_dim1, zr_offset, 
4936            zi_dim1, zi_offset, i_1, i_2, i_3;
4937    doublereal d_1, d_2, d_3, d_4;
4938
4939    /* Local variables */
4940    static integer iend;
4941    extern /* Subroutine */ int cdiv_(doublereal *, doublereal *, doublereal *
4942            , doublereal *, doublereal *, doublereal *);
4943    static doublereal norm;
4944    static integer i, j, k, l, m, ii, en, jj, ll, nn;
4945    static doublereal si, ti, xi, yi, sr, tr, xr, yr;
4946    extern doublereal pythag_(doublereal *, doublereal *);
4947    extern /* Subroutine */ int csroot_(doublereal *, doublereal *, 
4948            doublereal *, doublereal *);
4949    static integer ip1, lp1, itn, its;
4950    static doublereal zzi, zzr;
4951    static integer enm1;
4952    static doublereal tst1, tst2;
4953
4954
4955
4956/*     THIS SUBROUTINE IS A TRANSLATION OF A UNITARY ANALOGUE OF THE */
4957/*     ALGOL PROCEDURE  COMLR2, NUM. MATH. 16, 181-204(1970) BY PETERS */
4958/*     AND WILKINSON. */
4959/*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971). */
4960/*     THE UNITARY ANALOGUE SUBSTITUTES THE QR ALGORITHM OF FRANCIS */
4961/*     (COMP. JOUR. 4, 332-345(1962)) FOR THE LR ALGORITHM. */
4962
4963/*     THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS */
4964/*     OF A COMPLEX UPPER HESSENBERG MATRIX BY THE QR */
4965/*     METHOD.  THE EIGENVECTORS OF A COMPLEX GENERAL MATRIX */
4966/*     CAN ALSO BE FOUND IF  CORTH  HAS BEEN USED TO REDUCE */
4967/*     THIS GENERAL MATRIX TO HESSENBERG FORM. */
4968
4969/*     ON INPUT */
4970
4971/*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
4972/*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
4973/*          DIMENSION STATEMENT. */
4974
4975/*        N IS THE ORDER OF THE MATRIX. */
4976
4977/*        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
4978/*          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED, */
4979/*          SET LOW=1, IGH=N. */
4980
4981/*        ORTR AND ORTI CONTAIN INFORMATION ABOUT THE UNITARY TRANS- */
4982/*          FORMATIONS USED IN THE REDUCTION BY  CORTH, IF PERFORMED. */
4983/*          ONLY ELEMENTS LOW THROUGH IGH ARE USED.  IF THE EIGENVECTORS
4984*/
4985/*          OF THE HESSENBERG MATRIX ARE DESIRED, SET ORTR(J) AND */
4986/*          ORTI(J) TO 0.0D0 FOR THESE ELEMENTS. */
4987
4988/*        HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS, */
4989/*          RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX. */
4990/*          THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN FURTHER */
4991/*          INFORMATION ABOUT THE TRANSFORMATIONS WHICH WERE USED IN THE
4992*/
4993/*          REDUCTION BY  CORTH, IF PERFORMED.  IF THE EIGENVECTORS OF */
4994/*          THE HESSENBERG MATRIX ARE DESIRED, THESE ELEMENTS MAY BE */
4995/*          ARBITRARY. */
4996
4997/*     ON OUTPUT */
4998
4999/*        ORTR, ORTI, AND THE UPPER HESSENBERG PORTIONS OF HR AND HI */
5000/*          HAVE BEEN DESTROYED. */
5001
5002/*        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, */
5003/*          RESPECTIVELY, OF THE EIGENVALUES.  IF AN ERROR */
5004/*          EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT */
5005/*          FOR INDICES IERR+1,...,N. */
5006
5007/*        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, */
5008/*          RESPECTIVELY, OF THE EIGENVECTORS.  THE EIGENVECTORS */
5009/*          ARE UNNORMALIZED.  IF AN ERROR EXIT IS MADE, NONE OF */
5010/*          THE EIGENVECTORS HAS BEEN FOUND. */
5011
5012/*        IERR IS SET TO */
5013/*          ZERO       FOR NORMAL RETURN, */
5014/*          J          IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED */
5015/*                     WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. */
5016
5017/*     CALLS CDIV FOR COMPLEX DIVISION. */
5018/*     CALLS CSROOT FOR COMPLEX SQUARE ROOT. */
5019/*     CALLS PYTHAG FOR  DSQRT(A*A + B*B) . */
5020
5021/*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
5022/*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
5023*/
5024
5025/*     THIS VERSION DATED AUGUST 1983. */
5026
5027/*     ------------------------------------------------------------------
5028*/
5029
5030    /* Parameter adjustments */
5031    zi_dim1 = *nm;
5032    zi_offset = zi_dim1 + 1;
5033    zi -= zi_offset;
5034    zr_dim1 = *nm;
5035    zr_offset = zr_dim1 + 1;
5036    zr -= zr_offset;
5037    --wi;
5038    --wr;
5039    hi_dim1 = *nm;
5040    hi_offset = hi_dim1 + 1;
5041    hi -= hi_offset;
5042    hr_dim1 = *nm;
5043    hr_offset = hr_dim1 + 1;
5044    hr -= hr_offset;
5045    --orti;
5046    --ortr;
5047
5048    /* Function Body */
5049    *ierr = 0;
5050/*     .......... INITIALIZE EIGENVECTOR MATRIX .......... */
5051    i_1 = *n;
5052    for (j = 1; j <= i_1; ++j) {
5053
5054        i_2 = *n;
5055        for (i = 1; i <= i_2; ++i) {
5056            zr[i + j * zr_dim1] = 0.;
5057            zi[i + j * zi_dim1] = 0.;
5058/* L100: */
5059        }
5060        zr[j + j * zr_dim1] = 1.;
5061/* L101: */
5062    }
5063/*     .......... FORM THE MATRIX OF ACCUMULATED TRANSFORMATIONS */
5064/*                FROM THE INFORMATION LEFT BY CORTH .......... */
5065    iend = *igh - *low - 1;
5066    if (iend < 0) {
5067        goto L180;
5068    } else if (iend == 0) {
5069        goto L150;
5070    } else {
5071        goto L105;
5072    }
5073/*     .......... FOR I=IGH-1 STEP -1 UNTIL LOW+1 DO -- .......... */
5074L105:
5075    i_1 = iend;
5076    for (ii = 1; ii <= i_1; ++ii) {
5077        i = *igh - ii;
5078        if (ortr[i] == 0. && orti[i] == 0.) {
5079            goto L140;
5080        }
5081        if (hr[i + (i - 1) * hr_dim1] == 0. && hi[i + (i - 1) * hi_dim1] == 
5082                0.) {
5083            goto L140;
5084        }
5085/*     .......... NORM BELOW IS NEGATIVE OF H FORMED IN CORTH ........
5086.. */
5087        norm = hr[i + (i - 1) * hr_dim1] * ortr[i] + hi[i + (i - 1) * hi_dim1]
5088                 * orti[i];
5089        ip1 = i + 1;
5090
5091        i_2 = *igh;
5092        for (k = ip1; k <= i_2; ++k) {
5093            ortr[k] = hr[k + (i - 1) * hr_dim1];
5094            orti[k] = hi[k + (i - 1) * hi_dim1];
5095/* L110: */
5096        }
5097
5098        i_2 = *igh;
5099        for (j = i; j <= i_2; ++j) {
5100            sr = 0.;
5101            si = 0.;
5102
5103            i_3 = *igh;
5104            for (k = i; k <= i_3; ++k) {
5105                sr = sr + ortr[k] * zr[k + j * zr_dim1] + orti[k] * zi[k + j *
5106                         zi_dim1];
5107                si = si + ortr[k] * zi[k + j * zi_dim1] - orti[k] * zr[k + j *
5108                         zr_dim1];
5109/* L115: */
5110            }
5111
5112            sr /= norm;
5113            si /= norm;
5114
5115            i_3 = *igh;
5116            for (k = i; k <= i_3; ++k) {
5117                zr[k + j * zr_dim1] = zr[k + j * zr_dim1] + sr * ortr[k] - si
5118                        * orti[k];
5119                zi[k + j * zi_dim1] = zi[k + j * zi_dim1] + sr * orti[k] + si
5120                        * ortr[k];
5121/* L120: */
5122            }
5123
5124/* L130: */
5125        }
5126
5127L140:
5128        ;
5129    }
5130/*     .......... CREATE REAL SUBDIAGONAL ELEMENTS .......... */
5131L150:
5132    l = *low + 1;
5133
5134    i_1 = *igh;
5135    for (i = l; i <= i_1; ++i) {
5136/* Computing MIN */
5137        i_2 = i + 1;
5138        ll = min(i_2,*igh);
5139        if (hi[i + (i - 1) * hi_dim1] == 0.) {
5140            goto L170;
5141        }
5142        norm = pythag_(&hr[i + (i - 1) * hr_dim1], &hi[i + (i - 1) * hi_dim1])
5143                ;
5144        yr = hr[i + (i - 1) * hr_dim1] / norm;
5145        yi = hi[i + (i - 1) * hi_dim1] / norm;
5146        hr[i + (i - 1) * hr_dim1] = norm;
5147        hi[i + (i - 1) * hi_dim1] = 0.;
5148
5149        i_2 = *n;
5150        for (j = i; j <= i_2; ++j) {
5151            si = yr * hi[i + j * hi_dim1] - yi * hr[i + j * hr_dim1];
5152            hr[i + j * hr_dim1] = yr * hr[i + j * hr_dim1] + yi * hi[i + j * 
5153                    hi_dim1];
5154            hi[i + j * hi_dim1] = si;
5155/* L155: */
5156        }
5157
5158        i_2 = ll;
5159        for (j = 1; j <= i_2; ++j) {
5160            si = yr * hi[j + i * hi_dim1] + yi * hr[j + i * hr_dim1];
5161            hr[j + i * hr_dim1] = yr * hr[j + i * hr_dim1] - yi * hi[j + i * 
5162                    hi_dim1];
5163            hi[j + i * hi_dim1] = si;
5164/* L160: */
5165        }
5166
5167        i_2 = *igh;
5168        for (j = *low; j <= i_2; ++j) {
5169            si = yr * zi[j + i * zi_dim1] + yi * zr[j + i * zr_dim1];
5170            zr[j + i * zr_dim1] = yr * zr[j + i * zr_dim1] - yi * zi[j + i * 
5171                    zi_dim1];
5172            zi[j + i * zi_dim1] = si;
5173/* L165: */
5174        }
5175
5176L170:
5177        ;
5178    }
5179/*     .......... STORE ROOTS ISOLATED BY CBAL .......... */
5180L180:
5181    i_1 = *n;
5182    for (i = 1; i <= i_1; ++i) {
5183        if (i >= *low && i <= *igh) {
5184            goto L200;
5185        }
5186        wr[i] = hr[i + i * hr_dim1];
5187        wi[i] = hi[i + i * hi_dim1];
5188L200:
5189        ;
5190    }
5191
5192    en = *igh;
5193    tr = 0.;
5194    ti = 0.;
5195    itn = *n * 30;
5196/*     .......... SEARCH FOR NEXT EIGENVALUE .......... */
5197L220:
5198    if (en < *low) {
5199        goto L680;
5200    }
5201    its = 0;
5202    enm1 = en - 1;
5203/*     .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT */
5204/*                FOR L=EN STEP -1 UNTIL LOW DO -- .......... */
5205L240:
5206    i_1 = en;
5207    for (ll = *low; ll <= i_1; ++ll) {
5208        l = en + *low - ll;
5209        if (l == *low) {
5210            goto L300;
5211        }
5212        tst1 = (d_1 = hr[l - 1 + (l - 1) * hr_dim1], abs(d_1)) + (d_2 = hi[
5213                l - 1 + (l - 1) * hi_dim1], abs(d_2)) + (d_3 = hr[l + l * 
5214                hr_dim1], abs(d_3)) + (d_4 = hi[l + l * hi_dim1], abs(d_4))
5215                ;
5216        tst2 = tst1 + (d_1 = hr[l + (l - 1) * hr_dim1], abs(d_1));
5217        if (tst2 == tst1) {
5218            goto L300;
5219        }
5220/* L260: */
5221    }
5222/*     .......... FORM SHIFT .......... */
5223L300:
5224    if (l == en) {
5225        goto L660;
5226    }
5227    if (itn == 0) {
5228        goto L1000;
5229    }
5230    if (its == 10 || its == 20) {
5231        goto L320;
5232    }
5233    sr = hr[en + en * hr_dim1];
5234    si = hi[en + en * hi_dim1];
5235    xr = hr[enm1 + en * hr_dim1] * hr[en + enm1 * hr_dim1];
5236    xi = hi[enm1 + en * hi_dim1] * hr[en + enm1 * hr_dim1];
5237    if (xr == 0. && xi == 0.) {
5238        goto L340;
5239    }
5240    yr = (hr[enm1 + enm1 * hr_dim1] - sr) / 2.;
5241    yi = (hi[enm1 + enm1 * hi_dim1] - si) / 2.;
5242/* Computing 2nd power */
5243    d_2 = yr;
5244/* Computing 2nd power */
5245    d_3 = yi;
5246    d_1 = d_2 * d_2 - d_3 * d_3 + xr;
5247    d_4 = yr * 2. * yi + xi;
5248    csroot_(&d_1, &d_4, &zzr, &zzi);
5249    if (yr * zzr + yi * zzi >= 0.) {
5250        goto L310;
5251    }
5252    zzr = -zzr;
5253    zzi = -zzi;
5254L310:
5255    d_1 = yr + zzr;
5256    d_2 = yi + zzi;
5257    cdiv_(&xr, &xi, &d_1, &d_2, &xr, &xi);
5258    sr -= xr;
5259    si -= xi;
5260    goto L340;
5261/*     .......... FORM EXCEPTIONAL SHIFT .......... */
5262L320:
5263    sr = (d_1 = hr[en + enm1 * hr_dim1], abs(d_1)) + (d_2 = hr[enm1 + (en
5264            - 2) * hr_dim1], abs(d_2));
5265    si = 0.;
5266
5267L340:
5268    i_1 = en;
5269    for (i = *low; i <= i_1; ++i) {
5270        hr[i + i * hr_dim1] -= sr;
5271        hi[i + i * hi_dim1] -= si;
5272/* L360: */
5273    }
5274
5275    tr += sr;
5276    ti += si;
5277    ++its;
5278    --itn;
5279/*     .......... REDUCE TO TRIANGLE (ROWS) .......... */
5280    lp1 = l + 1;
5281
5282    i_1 = en;
5283    for (i = lp1; i <= i_1; ++i) {
5284        sr = hr[i + (i - 1) * hr_dim1];
5285        hr[i + (i - 1) * hr_dim1] = 0.;
5286        d_1 = pythag_(&hr[i - 1 + (i - 1) * hr_dim1], &hi[i - 1 + (i - 1) * 
5287                hi_dim1]);
5288        norm = pythag_(&d_1, &sr);
5289        xr = hr[i - 1 + (i - 1) * hr_dim1] / norm;
5290        wr[i - 1] = xr;
5291        xi = hi[i - 1 + (i - 1) * hi_dim1] / norm;
5292        wi[i - 1] = xi;
5293        hr[i - 1 + (i - 1) * hr_dim1] = norm;
5294        hi[i - 1 + (i - 1) * hi_dim1] = 0.;
5295        hi[i + (i - 1) * hi_dim1] = sr / norm;
5296
5297        i_2 = *n;
5298        for (j = i; j <= i_2; ++j) {
5299            yr = hr[i - 1 + j * hr_dim1];
5300            yi = hi[i - 1 + j * hi_dim1];
5301            zzr = hr[i + j * hr_dim1];
5302            zzi = hi[i + j * hi_dim1];
5303            hr[i - 1 + j * hr_dim1] = xr * yr + xi * yi + hi[i + (i - 1) * 
5304                    hi_dim1] * zzr;
5305            hi[i - 1 + j * hi_dim1] = xr * yi - xi * yr + hi[i + (i - 1) * 
5306                    hi_dim1] * zzi;
5307            hr[i + j * hr_dim1] = xr * zzr - xi * zzi - hi[i + (i - 1) * 
5308                    hi_dim1] * yr;
5309            hi[i + j * hi_dim1] = xr * zzi + xi * zzr - hi[i + (i - 1) * 
5310                    hi_dim1] * yi;
5311/* L490: */
5312        }
5313
5314/* L500: */
5315    }
5316
5317    si = hi[en + en * hi_dim1];
5318    if (si == 0.) {
5319        goto L540;
5320    }
5321    norm = pythag_(&hr[en + en * hr_dim1], &si);
5322    sr = hr[en + en * hr_dim1] / norm;
5323    si /= norm;
5324    hr[en + en * hr_dim1] = norm;
5325    hi[en + en * hi_dim1] = 0.;
5326    if (en == *n) {
5327        goto L540;
5328    }
5329    ip1 = en + 1;
5330
5331    i_1 = *n;
5332    for (j = ip1; j <= i_1; ++j) {
5333        yr = hr[en + j * hr_dim1];
5334        yi = hi[en + j * hi_dim1];
5335        hr[en + j * hr_dim1] = sr * yr + si * yi;
5336        hi[en + j * hi_dim1] = sr * yi - si * yr;
5337/* L520: */
5338    }
5339/*     .......... INVERSE OPERATION (COLUMNS) .......... */
5340L540:
5341    i_1 = en;
5342    for (j = lp1; j <= i_1; ++j) {
5343        xr = wr[j - 1];
5344        xi = wi[j - 1];
5345
5346        i_2 = j;
5347        for (i = 1; i <= i_2; ++i) {
5348            yr = hr[i + (j - 1) * hr_dim1];
5349            yi = 0.;
5350            zzr = hr[i + j * hr_dim1];
5351            zzi = hi[i + j * hi_dim1];
5352            if (i == j) {
5353                goto L560;
5354            }
5355            yi = hi[i + (j - 1) * hi_dim1];
5356            hi[i + (j - 1) * hi_dim1] = xr * yi + xi * yr + hi[j + (j - 1) * 
5357                    hi_dim1] * zzi;
5358L560:
5359            hr[i + (j - 1) * hr_dim1] = xr * yr - xi * yi + hi[j + (j - 1) * 
5360                    hi_dim1] * zzr;
5361            hr[i + j * hr_dim1] = xr * zzr + xi * zzi - hi[j + (j - 1) * 
5362                    hi_dim1] * yr;
5363            hi[i + j * hi_dim1] = xr * zzi - xi * zzr - hi[j + (j - 1) * 
5364                    hi_dim1] * yi;
5365/* L580: */
5366        }
5367
5368        i_2 = *igh;
5369        for (i = *low; i <= i_2; ++i) {
5370            yr = zr[i + (j - 1) * zr_dim1];
5371            yi = zi[i + (j - 1) * zi_dim1];
5372            zzr = zr[i + j * zr_dim1];
5373            zzi = zi[i + j * zi_dim1];
5374            zr[i + (j - 1) * zr_dim1] = xr * yr - xi * yi + hi[j + (j - 1) * 
5375                    hi_dim1] * zzr;
5376            zi[i + (j - 1) * zi_dim1] = xr * yi + xi * yr + hi[j + (j - 1) * 
5377                    hi_dim1] * zzi;
5378            zr[i + j * zr_dim1] = xr * zzr + xi * zzi - hi[j + (j - 1) * 
5379                    hi_dim1] * yr;
5380            zi[i + j * zi_dim1] = xr * zzi - xi * zzr - hi[j + (j - 1) * 
5381                    hi_dim1] * yi;
5382/* L590: */
5383        }
5384
5385/* L600: */
5386    }
5387
5388    if (si == 0.) {
5389        goto L240;
5390    }
5391
5392    i_1 = en;
5393    for (i = 1; i <= i_1; ++i) {
5394        yr = hr[i + en * hr_dim1];
5395        yi = hi[i + en * hi_dim1];
5396        hr[i + en * hr_dim1] = sr * yr - si * yi;
5397        hi[i + en * hi_dim1] = sr * yi + si * yr;
5398/* L630: */
5399    }
5400
5401    i_1 = *igh;
5402    for (i = *low; i <= i_1; ++i) {
5403        yr = zr[i + en * zr_dim1];
5404        yi = zi[i + en * zi_dim1];
5405        zr[i + en * zr_dim1] = sr * yr - si * yi;
5406        zi[i + en * zi_dim1] = sr * yi + si * yr;
5407/* L640: */
5408    }
5409
5410    goto L240;
5411/*     .......... A ROOT FOUND .......... */
5412L660:
5413    hr[en + en * hr_dim1] += tr;
5414    wr[en] = hr[en + en * hr_dim1];
5415    hi[en + en * hi_dim1] += ti;
5416    wi[en] = hi[en + en * hi_dim1];
5417    en = enm1;
5418    goto L220;
5419/*     .......... ALL ROOTS FOUND.  BACKSUBSTITUTE TO FIND */
5420/*                VECTORS OF UPPER TRIANGULAR FORM .......... */
5421L680:
5422    norm = 0.;
5423
5424    i_1 = *n;
5425    for (i = 1; i <= i_1; ++i) {
5426
5427        i_2 = *n;
5428        for (j = i; j <= i_2; ++j) {
5429            tr = (d_1 = hr[i + j * hr_dim1], abs(d_1)) + (d_2 = hi[i + j * 
5430                    hi_dim1], abs(d_2));
5431            if (tr > norm) {
5432                norm = tr;
5433            }
5434/* L720: */
5435        }
5436    }
5437
5438    if (*n == 1 || norm == 0.) {
5439        goto L1001;
5440    }
5441/*     .......... FOR EN=N STEP -1 UNTIL 2 DO -- .......... */
5442    i_2 = *n;
5443    for (nn = 2; nn <= i_2; ++nn) {
5444        en = *n + 2 - nn;
5445        xr = wr[en];
5446        xi = wi[en];
5447        hr[en + en * hr_dim1] = 1.;
5448        hi[en + en * hi_dim1] = 0.;
5449        enm1 = en - 1;
5450/*     .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- .......... */
5451        i_1 = enm1;
5452        for (ii = 1; ii <= i_1; ++ii) {
5453            i = en - ii;
5454            zzr = 0.;
5455            zzi = 0.;
5456            ip1 = i + 1;
5457
5458            i_3 = en;
5459            for (j = ip1; j <= i_3; ++j) {
5460                zzr = zzr + hr[i + j * hr_dim1] * hr[j + en * hr_dim1] - hi[i
5461                        + j * hi_dim1] * hi[j + en * hi_dim1];
5462                zzi = zzi + hr[i + j * hr_dim1] * hi[j + en * hi_dim1] + hi[i
5463                        + j * hi_dim1] * hr[j + en * hr_dim1];
5464/* L740: */
5465            }
5466
5467            yr = xr - wr[i];
5468            yi = xi - wi[i];
5469            if (yr != 0. || yi != 0.) {
5470                goto L765;
5471            }
5472            tst1 = norm;
5473            yr = tst1;
5474L760:
5475            yr *= .01;
5476            tst2 = norm + yr;
5477            if (tst2 > tst1) {
5478                goto L760;
5479            }
5480L765:
5481            cdiv_(&zzr, &zzi, &yr, &yi, &hr[i + en * hr_dim1], &hi[i + en * 
5482                    hi_dim1]);
5483/*     .......... OVERFLOW CONTROL .......... */
5484            tr = (d_1 = hr[i + en * hr_dim1], abs(d_1)) + (d_2 = hi[i + en
5485                    * hi_dim1], abs(d_2));
5486            if (tr == 0.) {
5487                goto L780;
5488            }
5489            tst1 = tr;
5490            tst2 = tst1 + 1. / tst1;
5491            if (tst2 > tst1) {
5492                goto L780;
5493            }
5494            i_3 = en;
5495            for (j = i; j <= i_3; ++j) {
5496                hr[j + en * hr_dim1] /= tr;
5497                hi[j + en * hi_dim1] /= tr;
5498/* L770: */
5499            }
5500
5501L780:
5502            ;
5503        }
5504
5505/* L800: */
5506    }
5507/*     .......... END BACKSUBSTITUTION .......... */
5508    enm1 = *n - 1;
5509/*     .......... VECTORS OF ISOLATED ROOTS .......... */
5510    i_2 = enm1;
5511    for (i = 1; i <= i_2; ++i) {
5512        if (i >= *low && i <= *igh) {
5513            goto L840;
5514        }
5515        ip1 = i + 1;
5516
5517        i_1 = *n;
5518        for (j = ip1; j <= i_1; ++j) {
5519            zr[i + j * zr_dim1] = hr[i + j * hr_dim1];
5520            zi[i + j * zi_dim1] = hi[i + j * hi_dim1];
5521/* L820: */
5522        }
5523
5524L840:
5525        ;
5526    }
5527/*     .......... MULTIPLY BY TRANSFORMATION MATRIX TO GIVE */
5528/*                VECTORS OF ORIGINAL FULL MATRIX. */
5529/*                FOR J=N STEP -1 UNTIL LOW+1 DO -- .......... */
5530    i_2 = enm1;
5531    for (jj = *low; jj <= i_2; ++jj) {
5532        j = *n + *low - jj;
5533        m = min(j,*igh);
5534
5535        i_1 = *igh;
5536        for (i = *low; i <= i_1; ++i) {
5537            zzr = 0.;
5538            zzi = 0.;
5539
5540            i_3 = m;
5541            for (k = *low; k <= i_3; ++k) {
5542                zzr = zzr + zr[i + k * zr_dim1] * hr[k + j * hr_dim1] - zi[i
5543                        + k * zi_dim1] * hi[k + j * hi_dim1];
5544                zzi = zzi + zr[i + k * zr_dim1] * hi[k + j * hi_dim1] + zi[i
5545                        + k * zi_dim1] * hr[k + j * hr_dim1];
5546/* L860: */
5547            }
5548
5549            zr[i + j * zr_dim1] = zzr;
5550            zi[i + j * zi_dim1] = zzi;
5551/* L880: */
5552        }
5553    }
5554
5555    goto L1001;
5556/*     .......... SET ERROR -- ALL EIGENVALUES HAVE NOT */
5557/*                CONVERGED AFTER 30*N ITERATIONS .......... */
5558L1000:
5559    *ierr = en;
5560L1001:
5561    return 0;
5562} /* comqr2_ */
5563
5564/* Subroutine */ int cortb_(integer *nm, integer *low, integer *igh, 
5565        doublereal *ar, doublereal *ai, doublereal *ortr, doublereal *orti, 
5566        integer *m, doublereal *zr, doublereal *zi)
5567{
5568    /* System generated locals */
5569    integer ar_dim1, ar_offset, ai_dim1, ai_offset, zr_dim1, zr_offset, 
5570            zi_dim1, zi_offset, i_1, i_2, i_3;
5571
5572    /* Local variables */
5573    static doublereal h;
5574    static integer i, j, la;
5575    static doublereal gi, gr;
5576    static integer mm, mp, kp1, mp1;
5577
5578
5579
5580/*     THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF */
5581/*     THE ALGOL PROCEDURE ORTBAK, NUM. MATH. 12, 349-368(1968) */
5582/*     BY MARTIN AND WILKINSON. */
5583/*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). */
5584
5585/*     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX GENERAL */
5586/*     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING */
5587/*     UPPER HESSENBERG MATRIX DETERMINED BY  CORTH. */
5588
5589/*     ON INPUT */
5590
5591/*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
5592/*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
5593/*          DIMENSION STATEMENT. */
5594
5595/*        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
5596/*          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED, */
5597/*          SET LOW=1 AND IGH EQUAL TO THE ORDER OF THE MATRIX. */
5598
5599/*        AR AND AI CONTAIN INFORMATION ABOUT THE UNITARY */
5600/*          TRANSFORMATIONS USED IN THE REDUCTION BY  CORTH */
5601/*          IN THEIR STRICT LOWER TRIANGLES. */
5602
5603/*        ORTR AND ORTI CONTAIN FURTHER INFORMATION ABOUT THE */
5604/*          TRANSFORMATIONS USED IN THE REDUCTION BY  CORTH. */
5605/*          ONLY ELEMENTS LOW THROUGH IGH ARE USED. */
5606
5607/*        M IS THE NUMBER OF COLUMNS OF ZR AND ZI TO BE BACK TRANSFORMED.
5608*/
5609
5610/*        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, */
5611/*          RESPECTIVELY, OF THE EIGENVECTORS TO BE */
5612/*          BACK TRANSFORMED IN THEIR FIRST M COLUMNS. */
5613
5614/*     ON OUTPUT */
5615
5616/*        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, */
5617/*          RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS */
5618/*          IN THEIR FIRST M COLUMNS. */
5619
5620/*        ORTR AND ORTI HAVE BEEN ALTERED. */
5621
5622/*     NOTE THAT CORTB PRESERVES VECTOR EUCLIDEAN NORMS. */
5623
5624/*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
5625/*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
5626*/
5627
5628/*     THIS VERSION DATED AUGUST 1983. */
5629
5630/*     ------------------------------------------------------------------
5631*/
5632
5633    /* Parameter adjustments */
5634    --orti;
5635    --ortr;
5636    ai_dim1 = *nm;
5637    ai_offset = ai_dim1 + 1;
5638    ai -= ai_offset;
5639    ar_dim1 = *nm;
5640    ar_offset = ar_dim1 + 1;
5641    ar -= ar_offset;
5642    zi_dim1 = *nm;
5643    zi_offset = zi_dim1 + 1;
5644    zi -= zi_offset;
5645    zr_dim1 = *nm;
5646    zr_offset = zr_dim1 + 1;
5647    zr -= zr_offset;
5648
5649    /* Function Body */
5650    if (*m == 0) {
5651        goto L200;
5652    }
5653    la = *igh - 1;
5654    kp1 = *low + 1;
5655    if (la < kp1) {
5656        goto L200;
5657    }
5658/*     .......... FOR MP=IGH-1 STEP -1 UNTIL LOW+1 DO -- .......... */
5659    i_1 = la;
5660    for (mm = kp1; mm <= i_1; ++mm) {
5661        mp = *low + *igh - mm;
5662        if (ar[mp + (mp - 1) * ar_dim1] == 0. && ai[mp + (mp - 1) * ai_dim1] 
5663                == 0.) {
5664            goto L140;
5665        }
5666/*     .......... H BELOW IS NEGATIVE OF H FORMED IN CORTH ..........
5667*/
5668        h = ar[mp + (mp - 1) * ar_dim1] * ortr[mp] + ai[mp + (mp - 1) * 
5669                ai_dim1] * orti[mp];
5670        mp1 = mp + 1;
5671
5672        i_2 = *igh;
5673        for (i = mp1; i <= i_2; ++i) {
5674            ortr[i] = ar[i + (mp - 1) * ar_dim1];
5675            orti[i] = ai[i + (mp - 1) * ai_dim1];
5676/* L100: */
5677        }
5678
5679        i_2 = *m;
5680        for (j = 1; j <= i_2; ++j) {
5681            gr = 0.;
5682            gi = 0.;
5683
5684            i_3 = *igh;
5685            for (i = mp; i <= i_3; ++i) {
5686                gr = gr + ortr[i] * zr[i + j * zr_dim1] + orti[i] * zi[i + j *
5687                         zi_dim1];
5688                gi = gi + ortr[i] * zi[i + j * zi_dim1] - orti[i] * zr[i + j *
5689                         zr_dim1];
5690/* L110: */
5691            }
5692
5693            gr /= h;
5694            gi /= h;
5695
5696            i_3 = *igh;
5697            for (i = mp; i <= i_3; ++i) {
5698                zr[i + j * zr_dim1] = zr[i + j * zr_dim1] + gr * ortr[i] - gi
5699                        * orti[i];
5700                zi[i + j * zi_dim1] = zi[i + j * zi_dim1] + gr * orti[i] + gi
5701                        * ortr[i];
5702/* L120: */
5703            }
5704
5705/* L130: */
5706        }
5707
5708L140:
5709        ;
5710    }
5711
5712L200:
5713    return 0;
5714} /* cortb_ */
5715
5716/* Subroutine */ int corth_(integer *nm, integer *n, integer *low, integer *
5717        igh, doublereal *ar, doublereal *ai, doublereal *ortr, doublereal *
5718        orti)
5719{
5720    /* System generated locals */
5721    integer ar_dim1, ar_offset, ai_dim1, ai_offset, i_1, i_2, i_3;
5722    doublereal d_1, d_2;
5723
5724    /* Builtin functions */
5725    double sqrt(doublereal);
5726
5727    /* Local variables */
5728    static doublereal f, g, h;
5729    static integer i, j, m;
5730    static doublereal scale;
5731    static integer la;
5732    static doublereal fi;
5733    static integer ii, jj;
5734    static doublereal fr;
5735    static integer mp;
5736    extern doublereal pythag_(doublereal *, doublereal *);
5737    static integer kp1;
5738
5739
5740
5741/*     THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF */
5742/*     THE ALGOL PROCEDURE ORTHES, NUM. MATH. 12, 349-368(1968) */
5743/*     BY MARTIN AND WILKINSON. */
5744/*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). */
5745
5746/*     GIVEN A COMPLEX GENERAL MATRIX, THIS SUBROUTINE */
5747/*     REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS */
5748/*     LOW THROUGH IGH TO UPPER HESSENBERG FORM BY */
5749/*     UNITARY SIMILARITY TRANSFORMATIONS. */
5750
5751/*     ON INPUT */
5752
5753/*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
5754/*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
5755/*          DIMENSION STATEMENT. */
5756
5757/*        N IS THE ORDER OF THE MATRIX. */
5758
5759/*        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
5760/*          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED, */
5761/*          SET LOW=1, IGH=N. */
5762
5763/*        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */
5764/*          RESPECTIVELY, OF THE COMPLEX INPUT MATRIX. */
5765
5766/*     ON OUTPUT */
5767
5768/*        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */
5769/*          RESPECTIVELY, OF THE HESSENBERG MATRIX.  INFORMATION */
5770/*          ABOUT THE UNITARY TRANSFORMATIONS USED IN THE REDUCTION */
5771/*          IS STORED IN THE REMAINING TRIANGLES UNDER THE */
5772/*          HESSENBERG MATRIX. */
5773
5774/*        ORTR AND ORTI CONTAIN FURTHER INFORMATION ABOUT THE */
5775/*          TRANSFORMATIONS.  ONLY ELEMENTS LOW THROUGH IGH ARE USED. */
5776
5777/*     CALLS PYTHAG FOR  DSQRT(A*A + B*B) . */
5778
5779/*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
5780/*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
5781*/
5782
5783/*     THIS VERSION DATED AUGUST 1983. */
5784
5785/*     ------------------------------------------------------------------
5786*/
5787
5788    /* Parameter adjustments */
5789    ai_dim1 = *nm;
5790    ai_offset = ai_dim1 + 1;
5791    ai -= ai_offset;
5792    ar_dim1 = *nm;
5793    ar_offset = ar_dim1 + 1;
5794    ar -= ar_offset;
5795    --orti;
5796    --ortr;
5797
5798    /* Function Body */
5799    la = *igh - 1;
5800    kp1 = *low + 1;
5801    if (la < kp1) {
5802        goto L200;
5803    }
5804
5805    i_1 = la;
5806    for (m = kp1; m <= i_1; ++m) {
5807        h = 0.;
5808        ortr[m] = 0.;
5809        orti[m] = 0.;
5810        scale = 0.;
5811/*     .......... SCALE COLUMN (ALGOL TOL THEN NOT NEEDED) ..........
5812*/
5813        i_2 = *igh;
5814        for (i = m; i <= i_2; ++i) {
5815/* L90: */
5816            scale = scale + (d_1 = ar[i + (m - 1) * ar_dim1], abs(d_1)) + (
5817                    d_2 = ai[i + (m - 1) * ai_dim1], abs(d_2));
5818        }
5819
5820        if (scale == 0.) {
5821            goto L180;
5822        }
5823        mp = m + *igh;
5824/*     .......... FOR I=IGH STEP -1 UNTIL M DO -- .......... */
5825        i_2 = *igh;
5826        for (ii = m; ii <= i_2; ++ii) {
5827            i = mp - ii;
5828            ortr[i] = ar[i + (m - 1) * ar_dim1] / scale;
5829            orti[i] = ai[i + (m - 1) * ai_dim1] / scale;
5830            h = h + ortr[i] * ortr[i] + orti[i] * orti[i];
5831/* L100: */
5832        }
5833
5834        g = sqrt(h);
5835        f = pythag_(&ortr[m], &orti[m]);
5836        if (f == 0.) {
5837            goto L103;
5838        }
5839        h += f * g;
5840        g /= f;
5841        ortr[m] = (g + 1.) * ortr[m];
5842        orti[m] = (g + 1.) * orti[m];
5843        goto L105;
5844
5845L103:
5846        ortr[m] = g;
5847        ar[m + (m - 1) * ar_dim1] = scale;
5848/*     .......... FORM (I-(U*UT)/H) * A .......... */
5849L105:
5850        i_2 = *n;
5851        for (j = m; j <= i_2; ++j) {
5852            fr = 0.;
5853            fi = 0.;
5854/*     .......... FOR I=IGH STEP -1 UNTIL M DO -- .......... */
5855            i_3 = *igh;
5856            for (ii = m; ii <= i_3; ++ii) {
5857                i = mp - ii;
5858                fr = fr + ortr[i] * ar[i + j * ar_dim1] + orti[i] * ai[i + j *
5859                         ai_dim1];
5860                fi = fi + ortr[i] * ai[i + j * ai_dim1] - orti[i] * ar[i + j *
5861                         ar_dim1];
5862/* L110: */
5863            }
5864
5865            fr /= h;
5866            fi /= h;
5867
5868            i_3 = *igh;
5869            for (i = m; i <= i_3; ++i) {
5870                ar[i + j * ar_dim1] = ar[i + j * ar_dim1] - fr * ortr[i] + fi
5871                        * orti[i];
5872                ai[i + j * ai_dim1] = ai[i + j * ai_dim1] - fr * orti[i] - fi
5873                        * ortr[i];
5874/* L120: */
5875            }
5876
5877/* L130: */
5878        }
5879/*     .......... FORM (I-(U*UT)/H)*A*(I-(U*UT)/H) .......... */
5880        i_2 = *igh;
5881        for (i = 1; i <= i_2; ++i) {
5882            fr = 0.;
5883            fi = 0.;
5884/*     .......... FOR J=IGH STEP -1 UNTIL M DO -- .......... */
5885            i_3 = *igh;
5886            for (jj = m; jj <= i_3; ++jj) {
5887                j = mp - jj;
5888                fr = fr + ortr[j] * ar[i + j * ar_dim1] - orti[j] * ai[i + j *
5889                         ai_dim1];
5890                fi = fi + ortr[j] * ai[i + j * ai_dim1] + orti[j] * ar[i + j *
5891                         ar_dim1];
5892/* L140: */
5893            }
5894
5895            fr /= h;
5896            fi /= h;
5897
5898            i_3 = *igh;
5899            for (j = m; j <= i_3; ++j) {
5900                ar[i + j * ar_dim1] = ar[i + j * ar_dim1] - fr * ortr[j] - fi
5901                        * orti[j];
5902                ai[i + j * ai_dim1] = ai[i + j * ai_dim1] + fr * orti[j] - fi
5903                        * ortr[j];
5904/* L150: */
5905            }
5906
5907/* L160: */
5908        }
5909
5910        ortr[m] = scale * ortr[m];
5911        orti[m] = scale * orti[m];
5912        ar[m + (m - 1) * ar_dim1] = -g * ar[m + (m - 1) * ar_dim1];
5913        ai[m + (m - 1) * ai_dim1] = -g * ai[m + (m - 1) * ai_dim1];
5914L180:
5915        ;
5916    }
5917
5918L200:
5919    return 0;
5920} /* corth_ */
5921
5922/* Subroutine */ int elmbak_(integer *nm, integer *low, integer *igh, 
5923        doublereal *a, integer *int_, integer *m, doublereal *z)
5924{
5925    /* System generated locals */
5926    integer a_dim1, a_offset, z_dim1, z_offset, i_1, i_2, i_3;
5927
5928    /* Local variables */
5929    static integer i, j;
5930    static doublereal x;
5931    static integer la, mm, mp, kp1, mp1;
5932
5933
5934
5935/*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ELMBAK, */
5936/*     NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON. */
5937/*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). */
5938
5939/*     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A REAL GENERAL */
5940/*     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING */
5941/*     UPPER HESSENBERG MATRIX DETERMINED BY  ELMHES. */
5942
5943/*     ON INPUT */
5944
5945/*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
5946/*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
5947/*          DIMENSION STATEMENT. */
5948
5949/*        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
5950/*          SUBROUTINE  BALANC.  IF  BALANC  HAS NOT BEEN USED, */
5951/*          SET LOW=1 AND IGH EQUAL TO THE ORDER OF THE MATRIX. */
5952
5953/*        A CONTAINS THE MULTIPLIERS WHICH WERE USED IN THE */
5954/*          REDUCTION BY  ELMHES  IN ITS LOWER TRIANGLE */
5955/*          BELOW THE SUBDIAGONAL. */
5956
5957/*        INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS */
5958/*          INTERCHANGED IN THE REDUCTION BY  ELMHES. */
5959/*          ONLY ELEMENTS LOW THROUGH IGH ARE USED. */
5960
5961/*        M IS THE NUMBER OF COLUMNS OF Z TO BE BACK TRANSFORMED. */
5962
5963/*        Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGEN- */
5964/*          VECTORS TO BE BACK TRANSFORMED IN ITS FIRST M COLUMNS. */
5965
5966/*     ON OUTPUT */
5967
5968/*        Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE */
5969/*          TRANSFORMED EIGENVECTORS IN ITS FIRST M COLUMNS. */
5970
5971/*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
5972/*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
5973*/
5974
5975/*     THIS VERSION DATED AUGUST 1983. */
5976
5977/*     ------------------------------------------------------------------
5978*/
5979
5980    /* Parameter adjustments */
5981    --int_;
5982    a_dim1 = *nm;
5983    a_offset = a_dim1 + 1;
5984    a -= a_offset;
5985    z_dim1 = *nm;
5986    z_offset = z_dim1 + 1;
5987    z -= z_offset;
5988
5989    /* Function Body */
5990    if (*m == 0) {
5991        goto L200;
5992    }
5993    la = *igh - 1;
5994    kp1 = *low + 1;
5995    if (la < kp1) {
5996        goto L200;
5997    }
5998/*     .......... FOR MP=IGH-1 STEP -1 UNTIL LOW+1 DO -- .......... */
5999    i_1 = la;
6000    for (mm = kp1; mm <= i_1; ++mm) {
6001        mp = *low + *igh - mm;
6002        mp1 = mp + 1;
6003
6004        i_2 = *igh;
6005        for (i = mp1; i <= i_2; ++i) {
6006            x = a[i + (mp - 1) * a_dim1];
6007            if (x == 0.) {
6008                goto L110;
6009            }
6010
6011            i_3 = *m;
6012            for (j = 1; j <= i_3; ++j) {
6013/* L100: */
6014                z[i + j * z_dim1] += x * z[mp + j * z_dim1];
6015            }
6016
6017L110:
6018            ;
6019        }
6020
6021        i = int_[mp];
6022        if (i == mp) {
6023            goto L140;
6024        }
6025
6026        i_2 = *m;
6027        for (j = 1; j <= i_2; ++j) {
6028            x = z[i + j * z_dim1];
6029            z[i + j * z_dim1] = z[mp + j * z_dim1];
6030            z[mp + j * z_dim1] = x;
6031/* L130: */
6032        }
6033
6034L140:
6035        ;
6036    }
6037
6038L200:
6039    return 0;
6040} /* elmbak_ */
6041
6042/* Subroutine */ int elmhes_(integer *nm, integer *n, integer *low, integer *
6043        igh, doublereal *a, integer *int_)
6044{
6045    /* System generated locals */
6046    integer a_dim1, a_offset, i_1, i_2, i_3;
6047    doublereal d_1;
6048
6049    /* Local variables */
6050    static integer i, j, m;
6051    static doublereal x, y;
6052    static integer la, mm1, kp1, mp1;
6053
6054
6055
6056/*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ELMHES, */
6057/*     NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON. */
6058/*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). */
6059
6060/*     GIVEN A REAL GENERAL MATRIX, THIS SUBROUTINE */
6061/*     REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS */
6062/*     LOW THROUGH IGH TO UPPER HESSENBERG FORM BY */
6063/*     STABILIZED ELEMENTARY SIMILARITY TRANSFORMATIONS. */
6064
6065/*     ON INPUT */
6066
6067/*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
6068/*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
6069/*          DIMENSION STATEMENT. */
6070
6071/*        N IS THE ORDER OF THE MATRIX. */
6072
6073/*        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
6074/*          SUBROUTINE  BALANC.  IF  BALANC  HAS NOT BEEN USED, */
6075/*          SET LOW=1, IGH=N. */
6076
6077/*        A CONTAINS THE INPUT MATRIX. */
6078
6079/*     ON OUTPUT */
6080
6081/*        A CONTAINS THE HESSENBERG MATRIX.  THE MULTIPLIERS */
6082/*          WHICH WERE USED IN THE REDUCTION ARE STORED IN THE */
6083/*          REMAINING TRIANGLE UNDER THE HESSENBERG MATRIX. */
6084
6085/*        INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS */
6086/*          INTERCHANGED IN THE REDUCTION. */
6087/*          ONLY ELEMENTS LOW THROUGH IGH ARE USED. */
6088
6089/*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
6090/*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
6091*/
6092
6093/*     THIS VERSION DATED AUGUST 1983. */
6094
6095/*     ------------------------------------------------------------------
6096*/
6097
6098    /* Parameter adjustments */
6099    a_dim1 = *nm;
6100    a_offset = a_dim1 + 1;
6101    a -= a_offset;
6102    --int_;
6103
6104    /* Function Body */
6105    la = *igh - 1;
6106    kp1 = *low + 1;
6107    if (la < kp1) {
6108        goto L200;
6109    }
6110
6111    i_1 = la;
6112    for (m = kp1; m <= i_1; ++m) {
6113        mm1 = m - 1;
6114        x = 0.;
6115        i = m;
6116
6117        i_2 = *igh;
6118        for (j = m; j <= i_2; ++j) {
6119            if ((d_1 = a[j + mm1 * a_dim1], abs(d_1)) <= abs(x)) {
6120                goto L100;
6121            }
6122            x = a[j + mm1 * a_dim1];
6123            i = j;
6124L100:
6125            ;
6126        }
6127
6128        int_[m] = i;
6129        if (i == m) {
6130            goto L130;
6131        }
6132/*     .......... INTERCHANGE ROWS AND COLUMNS OF A .......... */
6133        i_2 = *n;
6134        for (j = mm1; j <= i_2; ++j) {
6135            y = a[i + j * a_dim1];
6136            a[i + j * a_dim1] = a[m + j * a_dim1];
6137            a[m + j * a_dim1] = y;
6138/* L110: */
6139        }
6140
6141        i_2 = *igh;
6142        for (j = 1; j <= i_2; ++j) {
6143            y = a[j + i * a_dim1];
6144            a[j + i * a_dim1] = a[j + m * a_dim1];
6145            a[j + m * a_dim1] = y;
6146/* L120: */
6147        }
6148/*     .......... END INTERCHANGE .......... */
6149L130:
6150        if (x == 0.) {
6151            goto L180;
6152        }
6153        mp1 = m + 1;
6154
6155        i_2 = *igh;
6156        for (i = mp1; i <= i_2; ++i) {
6157            y = a[i + mm1 * a_dim1];
6158            if (y == 0.) {
6159                goto L160;
6160            }
6161            y /= x;
6162            a[i + mm1 * a_dim1] = y;
6163
6164            i_3 = *n;
6165            for (j = m; j <= i_3; ++j) {
6166/* L140: */
6167                a[i + j * a_dim1] -= y * a[m + j * a_dim1];
6168            }
6169
6170            i_3 = *igh;
6171            for (j = 1; j <= i_3; ++j) {
6172/* L150: */
6173                a[j + m * a_dim1] += y * a[j + i * a_dim1];
6174            }
6175
6176L160:
6177            ;
6178        }
6179
6180L180:
6181        ;
6182    }
6183
6184L200:
6185    return 0;
6186} /* elmhes_ */
6187
6188/* Subroutine */ int eltran_(integer *nm, integer *n, integer *low, integer *
6189        igh, doublereal *a, integer *int_, doublereal *z)
6190{
6191    /* System generated locals */
6192    integer a_dim1, a_offset, z_dim1, z_offset, i_1, i_2;
6193
6194    /* Local variables */
6195    static integer i, j, kl, mm, mp, mp1;
6196
6197
6198
6199/*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ELMTRANS,
6200*/
6201/*     NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON. */
6202/*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971). */
6203
6204/*     THIS SUBROUTINE ACCUMULATES THE STABILIZED ELEMENTARY */
6205/*     SIMILARITY TRANSFORMATIONS USED IN THE REDUCTION OF A */
6206/*     REAL GENERAL MATRIX TO UPPER HESSENBERG FORM BY  ELMHES. */
6207
6208/*     ON INPUT */
6209
6210/*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
6211/*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
6212/*          DIMENSION STATEMENT. */
6213
6214/*        N IS THE ORDER OF THE MATRIX. */
6215
6216/*        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
6217/*          SUBROUTINE  BALANC.  IF  BALANC  HAS NOT BEEN USED, */
6218/*          SET LOW=1, IGH=N. */
6219
6220