00011 {
00012     /* System generated locals */
00013     integer ar_dim1, ar_offset, ai_dim1, ai_offset, i__1, i__2, i__3;
00014     doublereal d__1, d__2;
00015 
00016     /* Builtin functions */
00017     double sqrt(doublereal);
00018 
00019     /* Local variables */
00020     static doublereal f, g, h__;
00021     static integer i__, j, m;
00022     static doublereal scale;
00023     static integer la;
00024     static doublereal fi;
00025     static integer ii, jj;
00026     static doublereal fr;
00027     static integer mp;
00028     extern doublereal pythag_(doublereal *, doublereal *);
00029     static integer kp1;
00030 
00031 
00032 
00033 /*     THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF */
00034 /*     THE ALGOL PROCEDURE ORTHES, NUM. MATH. 12, 349-368(1968) */
00035 /*     BY MARTIN AND WILKINSON. */
00036 /*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). */
00037 
00038 /*     GIVEN A COMPLEX GENERAL MATRIX, THIS SUBROUTINE */
00039 /*     REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS */
00040 /*     LOW THROUGH IGH TO UPPER HESSENBERG FORM BY */
00041 /*     UNITARY SIMILARITY TRANSFORMATIONS. */
00042 
00043 /*     ON INPUT */
00044 
00045 /*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
00046 /*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
00047 /*          DIMENSION STATEMENT. */
00048 
00049 /*        N IS THE ORDER OF THE MATRIX. */
00050 
00051 /*        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
00052 /*          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED, */
00053 /*          SET LOW=1, IGH=N. */
00054 
00055 /*        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */
00056 /*          RESPECTIVELY, OF THE COMPLEX INPUT MATRIX. */
00057 
00058 /*     ON OUTPUT */
00059 
00060 /*        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */
00061 /*          RESPECTIVELY, OF THE HESSENBERG MATRIX.  INFORMATION */
00062 /*          ABOUT THE UNITARY TRANSFORMATIONS USED IN THE REDUCTION */
00063 /*          IS STORED IN THE REMAINING TRIANGLES UNDER THE */
00064 /*          HESSENBERG MATRIX. */
00065 
00066 /*        ORTR AND ORTI CONTAIN FURTHER INFORMATION ABOUT THE */
00067 /*          TRANSFORMATIONS.  ONLY ELEMENTS LOW THROUGH IGH ARE USED. */
00068 
00069 /*     CALLS PYTHAG FOR  DSQRT(A*A + B*B) . */
00070 
00071 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00072 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00073 */
00074 
00075 /*     THIS VERSION DATED AUGUST 1983. */
00076 
00077 /*     ------------------------------------------------------------------ 
00078 */
00079 
00080     /* Parameter adjustments */
00081     ai_dim1 = *nm;
00082     ai_offset = ai_dim1 + 1;
00083     ai -= ai_offset;
00084     ar_dim1 = *nm;
00085     ar_offset = ar_dim1 + 1;
00086     ar -= ar_offset;
00087     --orti;
00088     --ortr;
00089 
00090     /* Function Body */
00091     la = *igh - 1;
00092     kp1 = *low + 1;
00093     if (la < kp1) {
00094         goto L200;
00095     }
00096 
00097     i__1 = la;
00098     for (m = kp1; m <= i__1; ++m) {
00099         h__ = 0.;
00100         ortr[m] = 0.;
00101         orti[m] = 0.;
00102         scale = 0.;
00103 /*     .......... SCALE COLUMN (ALGOL TOL THEN NOT NEEDED) .......... 
00104 */
00105         i__2 = *igh;
00106         for (i__ = m; i__ <= i__2; ++i__) {
00107 /* L90: */
00108             scale = scale + (d__1 = ar[i__ + (m - 1) * ar_dim1], abs(d__1)) + 
00109                     (d__2 = ai[i__ + (m - 1) * ai_dim1], abs(d__2));
00110         }
00111 
00112         if (scale == 0.) {
00113             goto L180;
00114         }
00115         mp = m + *igh;
00116 /*     .......... FOR I=IGH STEP -1 UNTIL M DO -- .......... */
00117         i__2 = *igh;
00118         for (ii = m; ii <= i__2; ++ii) {
00119             i__ = mp - ii;
00120             ortr[i__] = ar[i__ + (m - 1) * ar_dim1] / scale;
00121             orti[i__] = ai[i__ + (m - 1) * ai_dim1] / scale;
00122             h__ = h__ + ortr[i__] * ortr[i__] + orti[i__] * orti[i__];
00123 /* L100: */
00124         }
00125 
00126         g = sqrt(h__);
00127         f = pythag_(&ortr[m], &orti[m]);
00128         if (f == 0.) {
00129             goto L103;
00130         }
00131         h__ += f * g;
00132         g /= f;
00133         ortr[m] = (g + 1.) * ortr[m];
00134         orti[m] = (g + 1.) * orti[m];
00135         goto L105;
00136 
00137 L103:
00138         ortr[m] = g;
00139         ar[m + (m - 1) * ar_dim1] = scale;
00140 /*     .......... FORM (I-(U*UT)/H) * A .......... */
00141 L105:
00142         i__2 = *n;
00143         for (j = m; j <= i__2; ++j) {
00144             fr = 0.;
00145             fi = 0.;
00146 /*     .......... FOR I=IGH STEP -1 UNTIL M DO -- .......... */
00147             i__3 = *igh;
00148             for (ii = m; ii <= i__3; ++ii) {
00149                 i__ = mp - ii;
00150                 fr = fr + ortr[i__] * ar[i__ + j * ar_dim1] + orti[i__] * ai[
00151                         i__ + j * ai_dim1];
00152                 fi = fi + ortr[i__] * ai[i__ + j * ai_dim1] - orti[i__] * ar[
00153                         i__ + j * ar_dim1];
00154 /* L110: */
00155             }
00156 
00157             fr /= h__;
00158             fi /= h__;
00159 
00160             i__3 = *igh;
00161             for (i__ = m; i__ <= i__3; ++i__) {
00162                 ar[i__ + j * ar_dim1] = ar[i__ + j * ar_dim1] - fr * ortr[i__]
00163                          + fi * orti[i__];
00164                 ai[i__ + j * ai_dim1] = ai[i__ + j * ai_dim1] - fr * orti[i__]
00165                          - fi * ortr[i__];
00166 /* L120: */
00167             }
00168 
00169 /* L130: */
00170         }
00171 /*     .......... FORM (I-(U*UT)/H)*A*(I-(U*UT)/H) .......... */
00172         i__2 = *igh;
00173         for (i__ = 1; i__ <= i__2; ++i__) {
00174             fr = 0.;
00175             fi = 0.;
00176 /*     .......... FOR J=IGH STEP -1 UNTIL M DO -- .......... */
00177             i__3 = *igh;
00178             for (jj = m; jj <= i__3; ++jj) {
00179                 j = mp - jj;
00180                 fr = fr + ortr[j] * ar[i__ + j * ar_dim1] - orti[j] * ai[i__ 
00181                         + j * ai_dim1];
00182                 fi = fi + ortr[j] * ai[i__ + j * ai_dim1] + orti[j] * ar[i__ 
00183                         + j * ar_dim1];
00184 /* L140: */
00185             }
00186 
00187             fr /= h__;
00188             fi /= h__;
00189 
00190             i__3 = *igh;
00191             for (j = m; j <= i__3; ++j) {
00192                 ar[i__ + j * ar_dim1] = ar[i__ + j * ar_dim1] - fr * ortr[j] 
00193                         - fi * orti[j];
00194                 ai[i__ + j * ai_dim1] = ai[i__ + j * ai_dim1] + fr * orti[j] 
00195                         - fi * ortr[j];
00196 /* L150: */
00197             }
00198 
00199 /* L160: */
00200         }
00201 
00202         ortr[m] = scale * ortr[m];
00203         orti[m] = scale * orti[m];
00204         ar[m + (m - 1) * ar_dim1] = -g * ar[m + (m - 1) * ar_dim1];
00205         ai[m + (m - 1) * ai_dim1] = -g * ai[m + (m - 1) * ai_dim1];
00206 L180:
00207         ;
00208     }
00209 
00210 L200:
00211     return 0;
00212 } /* corth_ */