Doxygen Source Code Documentation
eis_comlr.c File Reference
#include "f2c.h"
Go to the source code of this file.
Functions | |
int | comlr_ (integer *nm, integer *n, integer *low, integer *igh, doublereal *hr, doublereal *hi, doublereal *wr, doublereal *wi, integer *ierr) |
Function Documentation
|
Definition at line 8 of file eis_comlr.c. References abs, cdiv_(), csroot_(), and l.
00011 { 00012 /* System generated locals */ 00013 integer hr_dim1, hr_offset, hi_dim1, hi_offset, i__1, i__2; 00014 doublereal d__1, d__2, d__3, d__4; 00015 00016 /* Local variables */ 00017 extern /* Subroutine */ int cdiv_(doublereal *, doublereal *, doublereal * 00018 , doublereal *, doublereal *, doublereal *); 00019 static integer i__, j, l, m, en, ll, mm; 00020 static doublereal si, ti, xi, yi, sr, tr, xr, yr; 00021 static integer im1; 00022 extern /* Subroutine */ int csroot_(doublereal *, doublereal *, 00023 doublereal *, doublereal *); 00024 static integer mp1, itn, its; 00025 static doublereal zzi, zzr; 00026 static integer enm1; 00027 static doublereal tst1, tst2; 00028 00029 00030 00031 /* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE COMLR, */ 00032 /* NUM. MATH. 12, 369-376(1968) BY MARTIN AND WILKINSON. */ 00033 /* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 396-403(1971). */ 00034 00035 /* THIS SUBROUTINE FINDS THE EIGENVALUES OF A COMPLEX */ 00036 /* UPPER HESSENBERG MATRIX BY THE MODIFIED LR METHOD. */ 00037 00038 /* ON INPUT */ 00039 00040 /* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */ 00041 /* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */ 00042 /* DIMENSION STATEMENT. */ 00043 00044 /* N IS THE ORDER OF THE MATRIX. */ 00045 00046 /* LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */ 00047 /* SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED, */ 00048 /* SET LOW=1, IGH=N. */ 00049 00050 /* HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS, */ 00051 /* RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX. */ 00052 /* THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN THE */ 00053 /* MULTIPLIERS WHICH WERE USED IN THE REDUCTION BY COMHES, */ 00054 /* IF PERFORMED. */ 00055 00056 /* ON OUTPUT */ 00057 00058 /* THE UPPER HESSENBERG PORTIONS OF HR AND HI HAVE BEEN */ 00059 /* DESTROYED. THEREFORE, THEY MUST BE SAVED BEFORE */ 00060 /* CALLING COMLR IF SUBSEQUENT CALCULATION OF */ 00061 /* EIGENVECTORS IS TO BE PERFORMED. */ 00062 00063 /* WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, */ 00064 /* RESPECTIVELY, OF THE EIGENVALUES. IF AN ERROR */ 00065 /* 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 /* CALLS CDIV FOR COMPLEX DIVISION. */ 00074 /* CALLS CSROOT FOR COMPLEX SQUARE ROOT. */ 00075 00076 /* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */ 00077 /* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 00078 */ 00079 00080 /* THIS VERSION DATED AUGUST 1983. */ 00081 00082 /* ------------------------------------------------------------------ 00083 */ 00084 00085 /* Parameter adjustments */ 00086 --wi; 00087 --wr; 00088 hi_dim1 = *nm; 00089 hi_offset = hi_dim1 + 1; 00090 hi -= hi_offset; 00091 hr_dim1 = *nm; 00092 hr_offset = hr_dim1 + 1; 00093 hr -= hr_offset; 00094 00095 /* Function Body */ 00096 *ierr = 0; 00097 /* .......... STORE ROOTS ISOLATED BY CBAL .......... */ 00098 i__1 = *n; 00099 for (i__ = 1; i__ <= i__1; ++i__) { 00100 if (i__ >= *low && i__ <= *igh) { 00101 goto L200; 00102 } 00103 wr[i__] = hr[i__ + i__ * hr_dim1]; 00104 wi[i__] = hi[i__ + i__ * hi_dim1]; 00105 L200: 00106 ; 00107 } 00108 00109 en = *igh; 00110 tr = 0.; 00111 ti = 0.; 00112 itn = *n * 30; 00113 /* .......... SEARCH FOR NEXT EIGENVALUE .......... */ 00114 L220: 00115 if (en < *low) { 00116 goto L1001; 00117 } 00118 its = 0; 00119 enm1 = en - 1; 00120 /* .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT */ 00121 /* FOR L=EN STEP -1 UNTIL LOW D0 -- .......... */ 00122 L240: 00123 i__1 = en; 00124 for (ll = *low; ll <= i__1; ++ll) { 00125 l = en + *low - ll; 00126 if (l == *low) { 00127 goto L300; 00128 } 00129 tst1 = (d__1 = hr[l - 1 + (l - 1) * hr_dim1], abs(d__1)) + (d__2 = hi[ 00130 l - 1 + (l - 1) * hi_dim1], abs(d__2)) + (d__3 = hr[l + l * 00131 hr_dim1], abs(d__3)) + (d__4 = hi[l + l * hi_dim1], abs(d__4)) 00132 ; 00133 tst2 = tst1 + (d__1 = hr[l + (l - 1) * hr_dim1], abs(d__1)) + (d__2 = 00134 hi[l + (l - 1) * hi_dim1], abs(d__2)); 00135 if (tst2 == tst1) { 00136 goto L300; 00137 } 00138 /* L260: */ 00139 } 00140 /* .......... FORM SHIFT .......... */ 00141 L300: 00142 if (l == en) { 00143 goto L660; 00144 } 00145 if (itn == 0) { 00146 goto L1000; 00147 } 00148 if (its == 10 || its == 20) { 00149 goto L320; 00150 } 00151 sr = hr[en + en * hr_dim1]; 00152 si = hi[en + en * hi_dim1]; 00153 xr = hr[enm1 + en * hr_dim1] * hr[en + enm1 * hr_dim1] - hi[enm1 + en * 00154 hi_dim1] * hi[en + enm1 * hi_dim1]; 00155 xi = hr[enm1 + en * hr_dim1] * hi[en + enm1 * hi_dim1] + hi[enm1 + en * 00156 hi_dim1] * hr[en + enm1 * hr_dim1]; 00157 if (xr == 0. && xi == 0.) { 00158 goto L340; 00159 } 00160 yr = (hr[enm1 + enm1 * hr_dim1] - sr) / 2.; 00161 yi = (hi[enm1 + enm1 * hi_dim1] - si) / 2.; 00162 /* Computing 2nd power */ 00163 d__2 = yr; 00164 /* Computing 2nd power */ 00165 d__3 = yi; 00166 d__1 = d__2 * d__2 - d__3 * d__3 + xr; 00167 d__4 = yr * 2. * yi + xi; 00168 csroot_(&d__1, &d__4, &zzr, &zzi); 00169 if (yr * zzr + yi * zzi >= 0.) { 00170 goto L310; 00171 } 00172 zzr = -zzr; 00173 zzi = -zzi; 00174 L310: 00175 d__1 = yr + zzr; 00176 d__2 = yi + zzi; 00177 cdiv_(&xr, &xi, &d__1, &d__2, &xr, &xi); 00178 sr -= xr; 00179 si -= xi; 00180 goto L340; 00181 /* .......... FORM EXCEPTIONAL SHIFT .......... */ 00182 L320: 00183 sr = (d__1 = hr[en + enm1 * hr_dim1], abs(d__1)) + (d__2 = hr[enm1 + (en 00184 - 2) * hr_dim1], abs(d__2)); 00185 si = (d__1 = hi[en + enm1 * hi_dim1], abs(d__1)) + (d__2 = hi[enm1 + (en 00186 - 2) * hi_dim1], abs(d__2)); 00187 00188 L340: 00189 i__1 = en; 00190 for (i__ = *low; i__ <= i__1; ++i__) { 00191 hr[i__ + i__ * hr_dim1] -= sr; 00192 hi[i__ + i__ * hi_dim1] -= si; 00193 /* L360: */ 00194 } 00195 00196 tr += sr; 00197 ti += si; 00198 ++its; 00199 --itn; 00200 /* .......... LOOK FOR TWO CONSECUTIVE SMALL */ 00201 /* SUB-DIAGONAL ELEMENTS .......... */ 00202 xr = (d__1 = hr[enm1 + enm1 * hr_dim1], abs(d__1)) + (d__2 = hi[enm1 + 00203 enm1 * hi_dim1], abs(d__2)); 00204 yr = (d__1 = hr[en + enm1 * hr_dim1], abs(d__1)) + (d__2 = hi[en + enm1 * 00205 hi_dim1], abs(d__2)); 00206 zzr = (d__1 = hr[en + en * hr_dim1], abs(d__1)) + (d__2 = hi[en + en * 00207 hi_dim1], abs(d__2)); 00208 /* .......... FOR M=EN-1 STEP -1 UNTIL L DO -- .......... */ 00209 i__1 = enm1; 00210 for (mm = l; mm <= i__1; ++mm) { 00211 m = enm1 + l - mm; 00212 if (m == l) { 00213 goto L420; 00214 } 00215 yi = yr; 00216 yr = (d__1 = hr[m + (m - 1) * hr_dim1], abs(d__1)) + (d__2 = hi[m + ( 00217 m - 1) * hi_dim1], abs(d__2)); 00218 xi = zzr; 00219 zzr = xr; 00220 xr = (d__1 = hr[m - 1 + (m - 1) * hr_dim1], abs(d__1)) + (d__2 = hi[m 00221 - 1 + (m - 1) * hi_dim1], abs(d__2)); 00222 tst1 = zzr / yi * (zzr + xr + xi); 00223 tst2 = tst1 + yr; 00224 if (tst2 == tst1) { 00225 goto L420; 00226 } 00227 /* L380: */ 00228 } 00229 /* .......... TRIANGULAR DECOMPOSITION H=L*R .......... */ 00230 L420: 00231 mp1 = m + 1; 00232 00233 i__1 = en; 00234 for (i__ = mp1; i__ <= i__1; ++i__) { 00235 im1 = i__ - 1; 00236 xr = hr[im1 + im1 * hr_dim1]; 00237 xi = hi[im1 + im1 * hi_dim1]; 00238 yr = hr[i__ + im1 * hr_dim1]; 00239 yi = hi[i__ + im1 * hi_dim1]; 00240 if (abs(xr) + abs(xi) >= abs(yr) + abs(yi)) { 00241 goto L460; 00242 } 00243 /* .......... INTERCHANGE ROWS OF HR AND HI .......... */ 00244 i__2 = en; 00245 for (j = im1; j <= i__2; ++j) { 00246 zzr = hr[im1 + j * hr_dim1]; 00247 hr[im1 + j * hr_dim1] = hr[i__ + j * hr_dim1]; 00248 hr[i__ + j * hr_dim1] = zzr; 00249 zzi = hi[im1 + j * hi_dim1]; 00250 hi[im1 + j * hi_dim1] = hi[i__ + j * hi_dim1]; 00251 hi[i__ + j * hi_dim1] = zzi; 00252 /* L440: */ 00253 } 00254 00255 cdiv_(&xr, &xi, &yr, &yi, &zzr, &zzi); 00256 wr[i__] = 1.; 00257 goto L480; 00258 L460: 00259 cdiv_(&yr, &yi, &xr, &xi, &zzr, &zzi); 00260 wr[i__] = -1.; 00261 L480: 00262 hr[i__ + im1 * hr_dim1] = zzr; 00263 hi[i__ + im1 * hi_dim1] = zzi; 00264 00265 i__2 = en; 00266 for (j = i__; j <= i__2; ++j) { 00267 hr[i__ + j * hr_dim1] = hr[i__ + j * hr_dim1] - zzr * hr[im1 + j * 00268 hr_dim1] + zzi * hi[im1 + j * hi_dim1]; 00269 hi[i__ + j * hi_dim1] = hi[i__ + j * hi_dim1] - zzr * hi[im1 + j * 00270 hi_dim1] - zzi * hr[im1 + j * hr_dim1]; 00271 /* L500: */ 00272 } 00273 00274 /* L520: */ 00275 } 00276 /* .......... COMPOSITION R*L=H .......... */ 00277 i__1 = en; 00278 for (j = mp1; j <= i__1; ++j) { 00279 xr = hr[j + (j - 1) * hr_dim1]; 00280 xi = hi[j + (j - 1) * hi_dim1]; 00281 hr[j + (j - 1) * hr_dim1] = 0.; 00282 hi[j + (j - 1) * hi_dim1] = 0.; 00283 /* .......... INTERCHANGE COLUMNS OF HR AND HI, */ 00284 /* IF NECESSARY .......... */ 00285 if (wr[j] <= 0.) { 00286 goto L580; 00287 } 00288 00289 i__2 = j; 00290 for (i__ = l; i__ <= i__2; ++i__) { 00291 zzr = hr[i__ + (j - 1) * hr_dim1]; 00292 hr[i__ + (j - 1) * hr_dim1] = hr[i__ + j * hr_dim1]; 00293 hr[i__ + j * hr_dim1] = zzr; 00294 zzi = hi[i__ + (j - 1) * hi_dim1]; 00295 hi[i__ + (j - 1) * hi_dim1] = hi[i__ + j * hi_dim1]; 00296 hi[i__ + j * hi_dim1] = zzi; 00297 /* L540: */ 00298 } 00299 00300 L580: 00301 i__2 = j; 00302 for (i__ = l; i__ <= i__2; ++i__) { 00303 hr[i__ + (j - 1) * hr_dim1] = hr[i__ + (j - 1) * hr_dim1] + xr * 00304 hr[i__ + j * hr_dim1] - xi * hi[i__ + j * hi_dim1]; 00305 hi[i__ + (j - 1) * hi_dim1] = hi[i__ + (j - 1) * hi_dim1] + xr * 00306 hi[i__ + j * hi_dim1] + xi * hr[i__ + j * hr_dim1]; 00307 /* L600: */ 00308 } 00309 00310 /* L640: */ 00311 } 00312 00313 goto L240; 00314 /* .......... A ROOT FOUND .......... */ 00315 L660: 00316 wr[en] = hr[en + en * hr_dim1] + tr; 00317 wi[en] = hi[en + en * hi_dim1] + ti; 00318 en = enm1; 00319 goto L220; 00320 /* .......... SET ERROR -- ALL EIGENVALUES HAVE NOT */ 00321 /* CONVERGED AFTER 30*N ITERATIONS .......... */ 00322 L1000: 00323 *ierr = en; 00324 L1001: 00325 return 0; 00326 } /* comlr_ */ |