Doxygen Source Code Documentation
eis_cbal.c File Reference
#include "f2c.h"
Go to the source code of this file.
Functions | |
int | cbal_ (integer *nm, integer *n, doublereal *ar, doublereal *ai, integer *low, integer *igh, doublereal *scale) |
Function Documentation
|
Definition at line 8 of file eis_cbal.c. Referenced by cg_().
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_ */ |