00001 /* cortb.f -- translated by f2c (version 19961017).
00002    You must link the resulting object file with the libraries:
00003         -lf2c -lm   (in that order)
00004 */
00005 
00006 #include "f2c.h"
00007 
00008 /* Subroutine */ int cortb_(integer *nm, integer *low, integer *igh, 
00009         doublereal *ar, doublereal *ai, doublereal *ortr, doublereal *orti, 
00010         integer *m, doublereal *zr, doublereal *zi)
00011 {
00012     /* System generated locals */
00013     integer ar_dim1, ar_offset, ai_dim1, ai_offset, zr_dim1, zr_offset, 
00014             zi_dim1, zi_offset, i__1, i__2, i__3;
00015 
00016     /* Local variables */
00017     static doublereal h__;
00018     static integer i__, j, la;
00019     static doublereal gi, gr;
00020     static integer mm, mp, kp1, mp1;
00021 
00022 
00023 
00024 /*     THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF */
00025 /*     THE ALGOL PROCEDURE ORTBAK, NUM. MATH. 12, 349-368(1968) */
00026 /*     BY MARTIN AND WILKINSON. */
00027 /*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). */
00028 
00029 /*     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX GENERAL */
00030 /*     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING */
00031 /*     UPPER HESSENBERG MATRIX DETERMINED BY  CORTH. */
00032 
00033 /*     ON INPUT */
00034 
00035 /*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
00036 /*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
00037 /*          DIMENSION STATEMENT. */
00038 
00039 /*        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
00040 /*          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED, */
00041 /*          SET LOW=1 AND IGH EQUAL TO THE ORDER OF THE MATRIX. */
00042 
00043 /*        AR AND AI CONTAIN INFORMATION ABOUT THE UNITARY */
00044 /*          TRANSFORMATIONS USED IN THE REDUCTION BY  CORTH */
00045 /*          IN THEIR STRICT LOWER TRIANGLES. */
00046 
00047 /*        ORTR AND ORTI CONTAIN FURTHER INFORMATION ABOUT THE */
00048 /*          TRANSFORMATIONS USED IN THE REDUCTION BY  CORTH. */
00049 /*          ONLY ELEMENTS LOW THROUGH IGH ARE USED. */
00050 
00051 /*        M IS THE NUMBER OF COLUMNS OF ZR AND ZI TO BE BACK TRANSFORMED. 
00052 */
00053 
00054 /*        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, */
00055 /*          RESPECTIVELY, OF THE EIGENVECTORS TO BE */
00056 /*          BACK TRANSFORMED IN THEIR FIRST M COLUMNS. */
00057 
00058 /*     ON OUTPUT */
00059 
00060 /*        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, */
00061 /*          RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS */
00062 /*          IN THEIR FIRST M COLUMNS. */
00063 
00064 /*        ORTR AND ORTI HAVE BEEN ALTERED. */
00065 
00066 /*     NOTE THAT CORTB PRESERVES VECTOR EUCLIDEAN NORMS. */
00067 
00068 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00069 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00070 */
00071 
00072 /*     THIS VERSION DATED AUGUST 1983. */
00073 
00074 /*     ------------------------------------------------------------------ 
00075 */
00076 
00077     /* Parameter adjustments */
00078     --orti;
00079     --ortr;
00080     ai_dim1 = *nm;
00081     ai_offset = ai_dim1 + 1;
00082     ai -= ai_offset;
00083     ar_dim1 = *nm;
00084     ar_offset = ar_dim1 + 1;
00085     ar -= ar_offset;
00086     zi_dim1 = *nm;
00087     zi_offset = zi_dim1 + 1;
00088     zi -= zi_offset;
00089     zr_dim1 = *nm;
00090     zr_offset = zr_dim1 + 1;
00091     zr -= zr_offset;
00092 
00093     /* Function Body */
00094     if (*m == 0) {
00095         goto L200;
00096     }
00097     la = *igh - 1;
00098     kp1 = *low + 1;
00099     if (la < kp1) {
00100         goto L200;
00101     }
00102 /*     .......... FOR MP=IGH-1 STEP -1 UNTIL LOW+1 DO -- .......... */
00103     i__1 = la;
00104     for (mm = kp1; mm <= i__1; ++mm) {
00105         mp = *low + *igh - mm;
00106         if (ar[mp + (mp - 1) * ar_dim1] == 0. && ai[mp + (mp - 1) * ai_dim1] 
00107                 == 0.) {
00108             goto L140;
00109         }
00110 /*     .......... H BELOW IS NEGATIVE OF H FORMED IN CORTH .......... 
00111 */
00112         h__ = ar[mp + (mp - 1) * ar_dim1] * ortr[mp] + ai[mp + (mp - 1) * 
00113                 ai_dim1] * orti[mp];
00114         mp1 = mp + 1;
00115 
00116         i__2 = *igh;
00117         for (i__ = mp1; i__ <= i__2; ++i__) {
00118             ortr[i__] = ar[i__ + (mp - 1) * ar_dim1];
00119             orti[i__] = ai[i__ + (mp - 1) * ai_dim1];
00120 /* L100: */
00121         }
00122 
00123         i__2 = *m;
00124         for (j = 1; j <= i__2; ++j) {
00125             gr = 0.;
00126             gi = 0.;
00127 
00128             i__3 = *igh;
00129             for (i__ = mp; i__ <= i__3; ++i__) {
00130                 gr = gr + ortr[i__] * zr[i__ + j * zr_dim1] + orti[i__] * zi[
00131                         i__ + j * zi_dim1];
00132                 gi = gi + ortr[i__] * zi[i__ + j * zi_dim1] - orti[i__] * zr[
00133                         i__ + j * zr_dim1];
00134 /* L110: */
00135             }
00136 
00137             gr /= h__;
00138             gi /= h__;
00139 
00140             i__3 = *igh;
00141             for (i__ = mp; i__ <= i__3; ++i__) {
00142                 zr[i__ + j * zr_dim1] = zr[i__ + j * zr_dim1] + gr * ortr[i__]
00143                          - gi * orti[i__];
00144                 zi[i__ + j * zi_dim1] = zi[i__ + j * zi_dim1] + gr * orti[i__]
00145                          + gi * ortr[i__];
00146 /* L120: */
00147             }
00148 
00149 /* L130: */
00150         }
00151 
00152 L140:
00153         ;
00154     }
00155 
00156 L200:
00157     return 0;
00158 } /* cortb_ */
00159