00001 /* combak.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 combak_(integer *nm, integer *low, integer *igh, 
00009         doublereal *ar, doublereal *ai, integer *int__, integer *m, 
00010         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 integer i__, j, la, mm, mp;
00018     static doublereal xi, xr;
00019     static integer kp1, mp1;
00020 
00021 
00022 
00023 /*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE COMBAK, */
00024 /*     NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON. */
00025 /*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). */
00026 
00027 /*     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX GENERAL */
00028 /*     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING */
00029 /*     UPPER HESSENBERG MATRIX DETERMINED BY  COMHES. */
00030 
00031 /*     ON INPUT */
00032 
00033 /*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
00034 /*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
00035 /*          DIMENSION STATEMENT. */
00036 
00037 /*        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
00038 /*          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED, */
00039 /*          SET LOW=1 AND IGH EQUAL TO THE ORDER OF THE MATRIX. */
00040 
00041 /*        AR AND AI CONTAIN THE MULTIPLIERS WHICH WERE USED IN THE */
00042 /*          REDUCTION BY  COMHES  IN THEIR LOWER TRIANGLES */
00043 /*          BELOW THE SUBDIAGONAL. */
00044 
00045 /*        INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS */
00046 /*          INTERCHANGED IN THE REDUCTION BY  COMHES. */
00047 /*          ONLY ELEMENTS LOW THROUGH IGH ARE USED. */
00048 
00049 /*        M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED. */
00050 
00051 /*        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, */
00052 /*          RESPECTIVELY, OF THE EIGENVECTORS TO BE */
00053 /*          BACK TRANSFORMED IN THEIR FIRST M COLUMNS. */
00054 
00055 /*     ON OUTPUT */
00056 
00057 /*        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, */
00058 /*          RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS */
00059 /*          IN THEIR FIRST M COLUMNS. */
00060 
00061 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00062 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00063 */
00064 
00065 /*     THIS VERSION DATED AUGUST 1983. */
00066 
00067 /*     ------------------------------------------------------------------ 
00068 */
00069 
00070     /* Parameter adjustments */
00071     --int__;
00072     ai_dim1 = *nm;
00073     ai_offset = ai_dim1 + 1;
00074     ai -= ai_offset;
00075     ar_dim1 = *nm;
00076     ar_offset = ar_dim1 + 1;
00077     ar -= ar_offset;
00078     zi_dim1 = *nm;
00079     zi_offset = zi_dim1 + 1;
00080     zi -= zi_offset;
00081     zr_dim1 = *nm;
00082     zr_offset = zr_dim1 + 1;
00083     zr -= zr_offset;
00084 
00085     /* Function Body */
00086     if (*m == 0) {
00087         goto L200;
00088     }
00089     la = *igh - 1;
00090     kp1 = *low + 1;
00091     if (la < kp1) {
00092         goto L200;
00093     }
00094 /*     .......... FOR MP=IGH-1 STEP -1 UNTIL LOW+1 DO -- .......... */
00095     i__1 = la;
00096     for (mm = kp1; mm <= i__1; ++mm) {
00097         mp = *low + *igh - mm;
00098         mp1 = mp + 1;
00099 
00100         i__2 = *igh;
00101         for (i__ = mp1; i__ <= i__2; ++i__) {
00102             xr = ar[i__ + (mp - 1) * ar_dim1];
00103             xi = ai[i__ + (mp - 1) * ai_dim1];
00104             if (xr == 0. && xi == 0.) {
00105                 goto L110;
00106             }
00107 
00108             i__3 = *m;
00109             for (j = 1; j <= i__3; ++j) {
00110                 zr[i__ + j * zr_dim1] = zr[i__ + j * zr_dim1] + xr * zr[mp + 
00111                         j * zr_dim1] - xi * zi[mp + j * zi_dim1];
00112                 zi[i__ + j * zi_dim1] = zi[i__ + j * zi_dim1] + xr * zi[mp + 
00113                         j * zi_dim1] + xi * zr[mp + j * zr_dim1];
00114 /* L100: */
00115             }
00116 
00117 L110:
00118             ;
00119         }
00120 
00121         i__ = int__[mp];
00122         if (i__ == mp) {
00123             goto L140;
00124         }
00125 
00126         i__2 = *m;
00127         for (j = 1; j <= i__2; ++j) {
00128             xr = zr[i__ + j * zr_dim1];
00129             zr[i__ + j * zr_dim1] = zr[mp + j * zr_dim1];
00130             zr[mp + j * zr_dim1] = xr;
00131             xi = zi[i__ + j * zi_dim1];
00132             zi[i__ + j * zi_dim1] = zi[mp + j * zi_dim1];
00133             zi[mp + j * zi_dim1] = xi;
00134 /* L130: */
00135         }
00136 
00137 L140:
00138         ;
00139     }
00140 
00141 L200:
00142     return 0;
00143 } /* combak_ */
00144