00001 /* cbal.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 cbal_(integer *nm, integer *n, doublereal *ar, 
00009         doublereal *ai, integer *low, integer *igh, doublereal *scale)
00010 {
00011     /* System generated locals */
00012     integer ar_dim1, ar_offset, ai_dim1, ai_offset, i__1, i__2;
00013     doublereal d__1, d__2;
00014 
00015     /* Local variables */
00016     static integer iexc;
00017     static doublereal c__, f, g;
00018     static integer i__, j, k, l, m;
00019     static doublereal r__, s, radix, b2;
00020     static integer jj;
00021     static logical noconv;
00022 
00023 
00024 
00025 /*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE */
00026 /*     CBALANCE, WHICH IS A COMPLEX VERSION OF BALANCE, */
00027 /*     NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH. */
00028 /*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971). */
00029 
00030 /*     THIS SUBROUTINE BALANCES A COMPLEX MATRIX AND ISOLATES */
00031 /*     EIGENVALUES WHENEVER POSSIBLE. */
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 /*        N IS THE ORDER OF THE MATRIX. */
00040 
00041 /*        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */
00042 /*          RESPECTIVELY, OF THE COMPLEX MATRIX TO BE BALANCED. */
00043 
00044 /*     ON OUTPUT */
00045 
00046 /*        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */
00047 /*          RESPECTIVELY, OF THE BALANCED MATRIX. */
00048 
00049 /*        LOW AND IGH ARE TWO INTEGERS SUCH THAT AR(I,J) AND AI(I,J) */
00050 /*          ARE EQUAL TO ZERO IF */
00051 /*           (1) I IS GREATER THAN J AND */
00052 /*           (2) J=1,...,LOW-1 OR I=IGH+1,...,N. */
00053 
00054 /*        SCALE CONTAINS INFORMATION DETERMINING THE */
00055 /*           PERMUTATIONS AND SCALING FACTORS USED. */
00056 
00057 /*     SUPPOSE THAT THE PRINCIPAL SUBMATRIX IN ROWS LOW THROUGH IGH */
00058 /*     HAS BEEN BALANCED, THAT P(J) DENOTES THE INDEX INTERCHANGED */
00059 /*     WITH J DURING THE PERMUTATION STEP, AND THAT THE ELEMENTS */
00060 /*     OF THE DIAGONAL MATRIX USED ARE DENOTED BY D(I,J).  THEN */
00061 /*        SCALE(J) = P(J),    FOR J = 1,...,LOW-1 */
00062 /*                 = D(J,J)       J = LOW,...,IGH */
00063 /*                 = P(J)         J = IGH+1,...,N. */
00064 /*     THE ORDER IN WHICH THE INTERCHANGES ARE MADE IS N TO IGH+1, */
00065 /*     THEN 1 TO LOW-1. */
00066 
00067 /*     NOTE THAT 1 IS RETURNED FOR IGH IF IGH IS ZERO FORMALLY. */
00068 
00069 /*     THE ALGOL PROCEDURE EXC CONTAINED IN CBALANCE APPEARS IN */
00070 /*     CBAL  IN LINE.  (NOTE THAT THE ALGOL ROLES OF IDENTIFIERS */
00071 /*     K,L HAVE BEEN REVERSED.) */
00072 
00073 /*     ARITHMETIC IS REAL THROUGHOUT. */
00074 
00075 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00076 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00077 */
00078 
00079 /*     THIS VERSION DATED AUGUST 1983. */
00080 
00081 /*     ------------------------------------------------------------------ 
00082 */
00083 
00084     /* Parameter adjustments */
00085     --scale;
00086     ai_dim1 = *nm;
00087     ai_offset = ai_dim1 + 1;
00088     ai -= ai_offset;
00089     ar_dim1 = *nm;
00090     ar_offset = ar_dim1 + 1;
00091     ar -= ar_offset;
00092 
00093     /* Function Body */
00094     radix = 16.;
00095 
00096     b2 = radix * radix;
00097     k = 1;
00098     l = *n;
00099     goto L100;
00100 /*     .......... IN-LINE PROCEDURE FOR ROW AND */
00101 /*                COLUMN EXCHANGE .......... */
00102 L20:
00103     scale[m] = (doublereal) j;
00104     if (j == m) {
00105         goto L50;
00106     }
00107 
00108     i__1 = l;
00109     for (i__ = 1; i__ <= i__1; ++i__) {
00110         f = ar[i__ + j * ar_dim1];
00111         ar[i__ + j * ar_dim1] = ar[i__ + m * ar_dim1];
00112         ar[i__ + m * ar_dim1] = f;
00113         f = ai[i__ + j * ai_dim1];
00114         ai[i__ + j * ai_dim1] = ai[i__ + m * ai_dim1];
00115         ai[i__ + m * ai_dim1] = f;
00116 /* L30: */
00117     }
00118 
00119     i__1 = *n;
00120     for (i__ = k; i__ <= i__1; ++i__) {
00121         f = ar[j + i__ * ar_dim1];
00122         ar[j + i__ * ar_dim1] = ar[m + i__ * ar_dim1];
00123         ar[m + i__ * ar_dim1] = f;
00124         f = ai[j + i__ * ai_dim1];
00125         ai[j + i__ * ai_dim1] = ai[m + i__ * ai_dim1];
00126         ai[m + i__ * ai_dim1] = f;
00127 /* L40: */
00128     }
00129 
00130 L50:
00131     switch (iexc) {
00132         case 1:  goto L80;
00133         case 2:  goto L130;
00134     }
00135 /*     .......... SEARCH FOR ROWS ISOLATING AN EIGENVALUE */
00136 /*                AND PUSH THEM DOWN .......... */
00137 L80:
00138     if (l == 1) {
00139         goto L280;
00140     }
00141     --l;
00142 /*     .......... FOR J=L STEP -1 UNTIL 1 DO -- .......... */
00143 L100:
00144     i__1 = l;
00145     for (jj = 1; jj <= i__1; ++jj) {
00146         j = l + 1 - jj;
00147 
00148         i__2 = l;
00149         for (i__ = 1; i__ <= i__2; ++i__) {
00150             if (i__ == j) {
00151                 goto L110;
00152             }
00153             if (ar[j + i__ * ar_dim1] != 0. || ai[j + i__ * ai_dim1] != 0.) {
00154                 goto L120;
00155             }
00156 L110:
00157             ;
00158         }
00159 
00160         m = l;
00161         iexc = 1;
00162         goto L20;
00163 L120:
00164         ;
00165     }
00166 
00167     goto L140;
00168 /*     .......... SEARCH FOR COLUMNS ISOLATING AN EIGENVALUE */
00169 /*                AND PUSH THEM LEFT .......... */
00170 L130:
00171     ++k;
00172 
00173 L140:
00174     i__1 = l;
00175     for (j = k; j <= i__1; ++j) {
00176 
00177         i__2 = l;
00178         for (i__ = k; i__ <= i__2; ++i__) {
00179             if (i__ == j) {
00180                 goto L150;
00181             }
00182             if (ar[i__ + j * ar_dim1] != 0. || ai[i__ + j * ai_dim1] != 0.) {
00183                 goto L170;
00184             }
00185 L150:
00186             ;
00187         }
00188 
00189         m = k;
00190         iexc = 2;
00191         goto L20;
00192 L170:
00193         ;
00194     }
00195 /*     .......... NOW BALANCE THE SUBMATRIX IN ROWS K TO L .......... */
00196     i__1 = l;
00197     for (i__ = k; i__ <= i__1; ++i__) {
00198 /* L180: */
00199         scale[i__] = 1.;
00200     }
00201 /*     .......... ITERATIVE LOOP FOR NORM REDUCTION .......... */
00202 L190:
00203     noconv = FALSE_;
00204 
00205     i__1 = l;
00206     for (i__ = k; i__ <= i__1; ++i__) {
00207         c__ = 0.;
00208         r__ = 0.;
00209 
00210         i__2 = l;
00211         for (j = k; j <= i__2; ++j) {
00212             if (j == i__) {
00213                 goto L200;
00214             }
00215             c__ = c__ + (d__1 = ar[j + i__ * ar_dim1], abs(d__1)) + (d__2 = 
00216                     ai[j + i__ * ai_dim1], abs(d__2));
00217             r__ = r__ + (d__1 = ar[i__ + j * ar_dim1], abs(d__1)) + (d__2 = 
00218                     ai[i__ + j * ai_dim1], abs(d__2));
00219 L200:
00220             ;
00221         }
00222 /*     .......... GUARD AGAINST ZERO C OR R DUE TO UNDERFLOW .........
00223 . */
00224         if (c__ == 0. || r__ == 0.) {
00225             goto L270;
00226         }
00227         g = r__ / radix;
00228         f = 1.;
00229         s = c__ + r__;
00230 L210:
00231         if (c__ >= g) {
00232             goto L220;
00233         }
00234         f *= radix;
00235         c__ *= b2;
00236         goto L210;
00237 L220:
00238         g = r__ * radix;
00239 L230:
00240         if (c__ < g) {
00241             goto L240;
00242         }
00243         f /= radix;
00244         c__ /= b2;
00245         goto L230;
00246 /*     .......... NOW BALANCE .......... */
00247 L240:
00248         if ((c__ + r__) / f >= s * .95) {
00249             goto L270;
00250         }
00251         g = 1. / f;
00252         scale[i__] *= f;
00253         noconv = TRUE_;
00254 
00255         i__2 = *n;
00256         for (j = k; j <= i__2; ++j) {
00257             ar[i__ + j * ar_dim1] *= g;
00258             ai[i__ + j * ai_dim1] *= g;
00259 /* L250: */
00260         }
00261 
00262         i__2 = l;
00263         for (j = 1; j <= i__2; ++j) {
00264             ar[j + i__ * ar_dim1] *= f;
00265             ai[j + i__ * ai_dim1] *= f;
00266 /* L260: */
00267         }
00268 
00269 L270:
00270         ;
00271     }
00272 
00273     if (noconv) {
00274         goto L190;
00275     }
00276 
00277 L280:
00278     *low = k;
00279     *igh = l;
00280     return 0;
00281 } /* cbal_ */
00282