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