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_ */
|
