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