Doxygen Source Code Documentation
eis_hqr.c File Reference
#include "f2c.h"
Go to the source code of this file.
Functions | |
int | hqr_ (integer *nm, integer *n, integer *low, integer *igh, doublereal *h__, doublereal *wr, doublereal *wi, integer *ierr) |
Function Documentation
|
Definition at line 8 of file eis_hqr.c. References abs, d_sign(), l, min, p, and q. Referenced by rg_().
00010 { 00011 /* System generated locals */ 00012 integer h_dim1, h_offset, i__1, i__2, i__3; 00013 doublereal d__1, d__2; 00014 00015 /* Builtin functions */ 00016 double sqrt(doublereal), d_sign(doublereal *, doublereal *); 00017 00018 /* Local variables */ 00019 static doublereal norm; 00020 static integer i__, j, k, l, m; 00021 static doublereal p, q, r__, s, t, w, x, y; 00022 static integer na, en, ll, mm; 00023 static doublereal zz; 00024 static logical notlas; 00025 static integer mp2, itn, its, enm2; 00026 static doublereal tst1, tst2; 00027 00028 00029 00030 /* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE HQR, */ 00031 /* NUM. MATH. 14, 219-231(1970) BY MARTIN, PETERS, AND WILKINSON. */ 00032 /* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 359-371(1971). */ 00033 00034 /* THIS SUBROUTINE FINDS THE EIGENVALUES OF A REAL */ 00035 /* UPPER HESSENBERG MATRIX BY THE QR METHOD. */ 00036 00037 /* ON INPUT */ 00038 00039 /* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */ 00040 /* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */ 00041 /* DIMENSION STATEMENT. */ 00042 00043 /* N IS THE ORDER OF THE MATRIX. */ 00044 00045 /* LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */ 00046 /* SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, */ 00047 /* SET LOW=1, IGH=N. */ 00048 00049 /* H CONTAINS THE UPPER HESSENBERG MATRIX. INFORMATION ABOUT */ 00050 /* THE TRANSFORMATIONS USED IN THE REDUCTION TO HESSENBERG */ 00051 /* FORM BY ELMHES OR ORTHES, IF PERFORMED, IS STORED */ 00052 /* IN THE REMAINING TRIANGLE UNDER THE HESSENBERG MATRIX. */ 00053 00054 /* ON OUTPUT */ 00055 00056 /* H HAS BEEN DESTROYED. THEREFORE, IT MUST BE SAVED */ 00057 /* BEFORE CALLING HQR IF SUBSEQUENT CALCULATION AND */ 00058 /* BACK TRANSFORMATION OF EIGENVECTORS IS TO BE PERFORMED. */ 00059 00060 /* WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, */ 00061 /* RESPECTIVELY, OF THE EIGENVALUES. THE EIGENVALUES */ 00062 /* ARE UNORDERED EXCEPT THAT COMPLEX CONJUGATE PAIRS */ 00063 /* OF VALUES APPEAR CONSECUTIVELY WITH THE EIGENVALUE */ 00064 /* HAVING THE POSITIVE IMAGINARY PART FIRST. IF AN */ 00065 /* ERROR EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT */ 00066 /* FOR INDICES IERR+1,...,N. */ 00067 00068 /* IERR IS SET TO */ 00069 /* ZERO FOR NORMAL RETURN, */ 00070 /* J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED */ 00071 /* WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. */ 00072 00073 /* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */ 00074 /* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 00075 */ 00076 00077 /* THIS VERSION DATED AUGUST 1983. */ 00078 00079 /* ------------------------------------------------------------------ 00080 */ 00081 00082 /* Parameter adjustments */ 00083 --wi; 00084 --wr; 00085 h_dim1 = *nm; 00086 h_offset = h_dim1 + 1; 00087 h__ -= h_offset; 00088 00089 /* Function Body */ 00090 *ierr = 0; 00091 norm = 0.; 00092 k = 1; 00093 /* .......... STORE ROOTS ISOLATED BY BALANC */ 00094 /* AND COMPUTE MATRIX NORM .......... */ 00095 i__1 = *n; 00096 for (i__ = 1; i__ <= i__1; ++i__) { 00097 00098 i__2 = *n; 00099 for (j = k; j <= i__2; ++j) { 00100 /* L40: */ 00101 norm += (d__1 = h__[i__ + j * h_dim1], abs(d__1)); 00102 } 00103 00104 k = i__; 00105 if (i__ >= *low && i__ <= *igh) { 00106 goto L50; 00107 } 00108 wr[i__] = h__[i__ + i__ * h_dim1]; 00109 wi[i__] = 0.; 00110 L50: 00111 ; 00112 } 00113 00114 en = *igh; 00115 t = 0.; 00116 itn = *n * 30; 00117 /* .......... SEARCH FOR NEXT EIGENVALUES .......... */ 00118 L60: 00119 if (en < *low) { 00120 goto L1001; 00121 } 00122 its = 0; 00123 na = en - 1; 00124 enm2 = na - 1; 00125 /* .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT */ 00126 /* FOR L=EN STEP -1 UNTIL LOW DO -- .......... */ 00127 L70: 00128 i__1 = en; 00129 for (ll = *low; ll <= i__1; ++ll) { 00130 l = en + *low - ll; 00131 if (l == *low) { 00132 goto L100; 00133 } 00134 s = (d__1 = h__[l - 1 + (l - 1) * h_dim1], abs(d__1)) + (d__2 = h__[l 00135 + l * h_dim1], abs(d__2)); 00136 if (s == 0.) { 00137 s = norm; 00138 } 00139 tst1 = s; 00140 tst2 = tst1 + (d__1 = h__[l + (l - 1) * h_dim1], abs(d__1)); 00141 if (tst2 == tst1) { 00142 goto L100; 00143 } 00144 /* L80: */ 00145 } 00146 /* .......... FORM SHIFT .......... */ 00147 L100: 00148 x = h__[en + en * h_dim1]; 00149 if (l == en) { 00150 goto L270; 00151 } 00152 y = h__[na + na * h_dim1]; 00153 w = h__[en + na * h_dim1] * h__[na + en * h_dim1]; 00154 if (l == na) { 00155 goto L280; 00156 } 00157 if (itn == 0) { 00158 goto L1000; 00159 } 00160 if (its != 10 && its != 20) { 00161 goto L130; 00162 } 00163 /* .......... FORM EXCEPTIONAL SHIFT .......... */ 00164 t += x; 00165 00166 i__1 = en; 00167 for (i__ = *low; i__ <= i__1; ++i__) { 00168 /* L120: */ 00169 h__[i__ + i__ * h_dim1] -= x; 00170 } 00171 00172 s = (d__1 = h__[en + na * h_dim1], abs(d__1)) + (d__2 = h__[na + enm2 * 00173 h_dim1], abs(d__2)); 00174 x = s * .75; 00175 y = x; 00176 w = s * -.4375 * s; 00177 L130: 00178 ++its; 00179 --itn; 00180 /* .......... LOOK FOR TWO CONSECUTIVE SMALL */ 00181 /* SUB-DIAGONAL ELEMENTS. */ 00182 /* FOR M=EN-2 STEP -1 UNTIL L DO -- .......... */ 00183 i__1 = enm2; 00184 for (mm = l; mm <= i__1; ++mm) { 00185 m = enm2 + l - mm; 00186 zz = h__[m + m * h_dim1]; 00187 r__ = x - zz; 00188 s = y - zz; 00189 p = (r__ * s - w) / h__[m + 1 + m * h_dim1] + h__[m + (m + 1) * 00190 h_dim1]; 00191 q = h__[m + 1 + (m + 1) * h_dim1] - zz - r__ - s; 00192 r__ = h__[m + 2 + (m + 1) * h_dim1]; 00193 s = abs(p) + abs(q) + abs(r__); 00194 p /= s; 00195 q /= s; 00196 r__ /= s; 00197 if (m == l) { 00198 goto L150; 00199 } 00200 tst1 = abs(p) * ((d__1 = h__[m - 1 + (m - 1) * h_dim1], abs(d__1)) + 00201 abs(zz) + (d__2 = h__[m + 1 + (m + 1) * h_dim1], abs(d__2))); 00202 tst2 = tst1 + (d__1 = h__[m + (m - 1) * h_dim1], abs(d__1)) * (abs(q) 00203 + abs(r__)); 00204 if (tst2 == tst1) { 00205 goto L150; 00206 } 00207 /* L140: */ 00208 } 00209 00210 L150: 00211 mp2 = m + 2; 00212 00213 i__1 = en; 00214 for (i__ = mp2; i__ <= i__1; ++i__) { 00215 h__[i__ + (i__ - 2) * h_dim1] = 0.; 00216 if (i__ == mp2) { 00217 goto L160; 00218 } 00219 h__[i__ + (i__ - 3) * h_dim1] = 0.; 00220 L160: 00221 ; 00222 } 00223 /* .......... DOUBLE QR STEP INVOLVING ROWS L TO EN AND */ 00224 /* COLUMNS M TO EN .......... */ 00225 i__1 = na; 00226 for (k = m; k <= i__1; ++k) { 00227 notlas = k != na; 00228 if (k == m) { 00229 goto L170; 00230 } 00231 p = h__[k + (k - 1) * h_dim1]; 00232 q = h__[k + 1 + (k - 1) * h_dim1]; 00233 r__ = 0.; 00234 if (notlas) { 00235 r__ = h__[k + 2 + (k - 1) * h_dim1]; 00236 } 00237 x = abs(p) + abs(q) + abs(r__); 00238 if (x == 0.) { 00239 goto L260; 00240 } 00241 p /= x; 00242 q /= x; 00243 r__ /= x; 00244 L170: 00245 d__1 = sqrt(p * p + q * q + r__ * r__); 00246 s = d_sign(&d__1, &p); 00247 if (k == m) { 00248 goto L180; 00249 } 00250 h__[k + (k - 1) * h_dim1] = -s * x; 00251 goto L190; 00252 L180: 00253 if (l != m) { 00254 h__[k + (k - 1) * h_dim1] = -h__[k + (k - 1) * h_dim1]; 00255 } 00256 L190: 00257 p += s; 00258 x = p / s; 00259 y = q / s; 00260 zz = r__ / s; 00261 q /= p; 00262 r__ /= p; 00263 if (notlas) { 00264 goto L225; 00265 } 00266 /* .......... ROW MODIFICATION .......... */ 00267 i__2 = *n; 00268 for (j = k; j <= i__2; ++j) { 00269 p = h__[k + j * h_dim1] + q * h__[k + 1 + j * h_dim1]; 00270 h__[k + j * h_dim1] -= p * x; 00271 h__[k + 1 + j * h_dim1] -= p * y; 00272 /* L200: */ 00273 } 00274 00275 /* Computing MIN */ 00276 i__2 = en, i__3 = k + 3; 00277 j = min(i__2,i__3); 00278 /* .......... COLUMN MODIFICATION .......... */ 00279 i__2 = j; 00280 for (i__ = 1; i__ <= i__2; ++i__) { 00281 p = x * h__[i__ + k * h_dim1] + y * h__[i__ + (k + 1) * h_dim1]; 00282 h__[i__ + k * h_dim1] -= p; 00283 h__[i__ + (k + 1) * h_dim1] -= p * q; 00284 /* L210: */ 00285 } 00286 goto L255; 00287 L225: 00288 /* .......... ROW MODIFICATION .......... */ 00289 i__2 = *n; 00290 for (j = k; j <= i__2; ++j) { 00291 p = h__[k + j * h_dim1] + q * h__[k + 1 + j * h_dim1] + r__ * h__[ 00292 k + 2 + j * h_dim1]; 00293 h__[k + j * h_dim1] -= p * x; 00294 h__[k + 1 + j * h_dim1] -= p * y; 00295 h__[k + 2 + j * h_dim1] -= p * zz; 00296 /* L230: */ 00297 } 00298 00299 /* Computing MIN */ 00300 i__2 = en, i__3 = k + 3; 00301 j = min(i__2,i__3); 00302 /* .......... COLUMN MODIFICATION .......... */ 00303 i__2 = j; 00304 for (i__ = 1; i__ <= i__2; ++i__) { 00305 p = x * h__[i__ + k * h_dim1] + y * h__[i__ + (k + 1) * h_dim1] + 00306 zz * h__[i__ + (k + 2) * h_dim1]; 00307 h__[i__ + k * h_dim1] -= p; 00308 h__[i__ + (k + 1) * h_dim1] -= p * q; 00309 h__[i__ + (k + 2) * h_dim1] -= p * r__; 00310 /* L240: */ 00311 } 00312 L255: 00313 00314 L260: 00315 ; 00316 } 00317 00318 goto L70; 00319 /* .......... ONE ROOT FOUND .......... */ 00320 L270: 00321 wr[en] = x + t; 00322 wi[en] = 0.; 00323 en = na; 00324 goto L60; 00325 /* .......... TWO ROOTS FOUND .......... */ 00326 L280: 00327 p = (y - x) / 2.; 00328 q = p * p + w; 00329 zz = sqrt((abs(q))); 00330 x += t; 00331 if (q < 0.) { 00332 goto L320; 00333 } 00334 /* .......... REAL PAIR .......... */ 00335 zz = p + d_sign(&zz, &p); 00336 wr[na] = x + zz; 00337 wr[en] = wr[na]; 00338 if (zz != 0.) { 00339 wr[en] = x - w / zz; 00340 } 00341 wi[na] = 0.; 00342 wi[en] = 0.; 00343 goto L330; 00344 /* .......... COMPLEX PAIR .......... */ 00345 L320: 00346 wr[na] = x + p; 00347 wr[en] = x + p; 00348 wi[na] = zz; 00349 wi[en] = -zz; 00350 L330: 00351 en = enm2; 00352 goto L60; 00353 /* .......... SET ERROR -- ALL EIGENVALUES HAVE NOT */ 00354 /* CONVERGED AFTER 30*N ITERATIONS .......... */ 00355 L1000: 00356 *ierr = en; 00357 L1001: 00358 return 0; 00359 } /* hqr_ */ |