Doxygen Source Code Documentation
eis_bqr.c File Reference
#include "f2c.h"Go to the source code of this file.
Functions | |
| int | bqr_ (integer *nm, integer *n, integer *mb, doublereal *a, doublereal *t, doublereal *r__, integer *ierr, integer *nv, doublereal *rv) |
Variables | |
| doublereal | c_b8 = 1. |
Function Documentation
|
||||||||||||||||||||||||||||||||||||||||
|
Definition at line 12 of file eis_bqr.c. References a, abs, c_b8, d_sign(), l, m1, m2, max, min, pythag_(), q, and scale.
00015 {
00016 /* System generated locals */
00017 integer a_dim1, a_offset, i__1, i__2, i__3;
00018 doublereal d__1;
00019
00020 /* Builtin functions */
00021 double d_sign(doublereal *, doublereal *), sqrt(doublereal);
00022
00023 /* Local variables */
00024 static doublereal f, g;
00025 static integer i__, j, k, l, m;
00026 static doublereal q, s, scale;
00027 static integer imult, m1, m2, m3, m4, m21, m31, ii, ik, jk, kj, jm, kk,
00028 km, ll, mk, mn, ni, mz;
00029 extern doublereal pythag_(doublereal *, doublereal *);
00030 static integer kj1, its;
00031 static doublereal tst1, tst2;
00032
00033
00034
00035 /* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BQR, */
00036 /* NUM. MATH. 16, 85-92(1970) BY MARTIN, REINSCH, AND WILKINSON. */
00037 /* HANDBOOK FOR AUTO. COMP., VOL II-LINEAR ALGEBRA, 266-272(1971). */
00038
00039 /* THIS SUBROUTINE FINDS THE EIGENVALUE OF SMALLEST (USUALLY) */
00040 /* MAGNITUDE OF A REAL SYMMETRIC BAND MATRIX USING THE */
00041 /* QR ALGORITHM WITH SHIFTS OF ORIGIN. CONSECUTIVE CALLS */
00042 /* CAN BE MADE TO FIND FURTHER EIGENVALUES. */
00043
00044 /* ON INPUT */
00045
00046 /* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
00047 /* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
00048 /* DIMENSION STATEMENT. */
00049
00050 /* N IS THE ORDER OF THE MATRIX. */
00051
00052 /* MB IS THE (HALF) BAND WIDTH OF THE MATRIX, DEFINED AS THE */
00053 /* NUMBER OF ADJACENT DIAGONALS, INCLUDING THE PRINCIPAL */
00054 /* DIAGONAL, REQUIRED TO SPECIFY THE NON-ZERO PORTION OF THE */
00055 /* LOWER TRIANGLE OF THE MATRIX. */
00056
00057 /* A CONTAINS THE LOWER TRIANGLE OF THE SYMMETRIC BAND INPUT */
00058 /* MATRIX STORED AS AN N BY MB ARRAY. ITS LOWEST SUBDIAGONAL */
00059 /* IS STORED IN THE LAST N+1-MB POSITIONS OF THE FIRST COLUMN, */
00060 /* ITS NEXT SUBDIAGONAL IN THE LAST N+2-MB POSITIONS OF THE */
00061 /* SECOND COLUMN, FURTHER SUBDIAGONALS SIMILARLY, AND FINALLY */
00062 /* ITS PRINCIPAL DIAGONAL IN THE N POSITIONS OF THE LAST COLUMN.
00063 */
00064 /* CONTENTS OF STORAGES NOT PART OF THE MATRIX ARE ARBITRARY. */
00065 /* ON A SUBSEQUENT CALL, ITS OUTPUT CONTENTS FROM THE PREVIOUS */
00066 /* CALL SHOULD BE PASSED. */
00067
00068 /* T SPECIFIES THE SHIFT (OF EIGENVALUES) APPLIED TO THE DIAGONAL
00069 */
00070 /* OF A IN FORMING THE INPUT MATRIX. WHAT IS ACTUALLY DETERMINED
00071 */
00072 /* IS THE EIGENVALUE OF A+TI (I IS THE IDENTITY MATRIX) NEAREST
00073 */
00074 /* TO T. ON A SUBSEQUENT CALL, THE OUTPUT VALUE OF T FROM THE */
00075 /* PREVIOUS CALL SHOULD BE PASSED IF THE NEXT NEAREST EIGENVALUE
00076 */
00077 /* IS SOUGHT. */
00078
00079 /* R SHOULD BE SPECIFIED AS ZERO ON THE FIRST CALL, AND AS ITS */
00080 /* OUTPUT VALUE FROM THE PREVIOUS CALL ON A SUBSEQUENT CALL. */
00081 /* IT IS USED TO DETERMINE WHEN THE LAST ROW AND COLUMN OF */
00082 /* THE TRANSFORMED BAND MATRIX CAN BE REGARDED AS NEGLIGIBLE. */
00083
00084 /* NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER RV */
00085 /* AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT. */
00086
00087 /* ON OUTPUT */
00088
00089 /* A CONTAINS THE TRANSFORMED BAND MATRIX. THE MATRIX A+TI */
00090 /* DERIVED FROM THE OUTPUT PARAMETERS IS SIMILAR TO THE */
00091 /* INPUT A+TI TO WITHIN ROUNDING ERRORS. ITS LAST ROW AND */
00092 /* COLUMN ARE NULL (IF IERR IS ZERO). */
00093
00094 /* T CONTAINS THE COMPUTED EIGENVALUE OF A+TI (IF IERR IS ZERO). */
00095
00096 /* R CONTAINS THE MAXIMUM OF ITS INPUT VALUE AND THE NORM OF THE */
00097 /* LAST COLUMN OF THE INPUT MATRIX A. */
00098
00099 /* IERR IS SET TO */
00100 /* ZERO FOR NORMAL RETURN, */
00101 /* N IF THE EIGENVALUE HAS NOT BEEN */
00102 /* DETERMINED AFTER 30 ITERATIONS. */
00103
00104 /* RV IS A TEMPORARY STORAGE ARRAY OF DIMENSION AT LEAST */
00105 /* (2*MB**2+4*MB-3). THE FIRST (3*MB-2) LOCATIONS CORRESPOND */
00106 /* TO THE ALGOL ARRAY B, THE NEXT (2*MB-1) LOCATIONS CORRESPOND
00107 */
00108 /* TO THE ALGOL ARRAY H, AND THE FINAL (2*MB**2-MB) LOCATIONS */
00109 /* CORRESPOND TO THE MB BY (2*MB-1) ALGOL ARRAY U. */
00110
00111 /* NOTE. FOR A SUBSEQUENT CALL, N SHOULD BE REPLACED BY N-1, BUT */
00112 /* MB SHOULD NOT BE ALTERED EVEN WHEN IT EXCEEDS THE CURRENT N. */
00113
00114 /* CALLS PYTHAG FOR DSQRT(A*A + B*B) . */
00115
00116 /* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00117 /* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
00118 */
00119
00120 /* THIS VERSION DATED AUGUST 1983. */
00121
00122 /* ------------------------------------------------------------------
00123 */
00124
00125 /* Parameter adjustments */
00126 a_dim1 = *nm;
00127 a_offset = a_dim1 + 1;
00128 a -= a_offset;
00129 --rv;
00130
00131 /* Function Body */
00132 *ierr = 0;
00133 m1 = min(*mb,*n);
00134 m = m1 - 1;
00135 m2 = m + m;
00136 m21 = m2 + 1;
00137 m3 = m21 + m;
00138 m31 = m3 + 1;
00139 m4 = m31 + m2;
00140 mn = m + *n;
00141 mz = *mb - m1;
00142 its = 0;
00143 /* .......... TEST FOR CONVERGENCE .......... */
00144 L40:
00145 g = a[*n + *mb * a_dim1];
00146 if (m == 0) {
00147 goto L360;
00148 }
00149 f = 0.;
00150
00151 i__1 = m;
00152 for (k = 1; k <= i__1; ++k) {
00153 mk = k + mz;
00154 f += (d__1 = a[*n + mk * a_dim1], abs(d__1));
00155 /* L50: */
00156 }
00157
00158 if (its == 0 && f > *r__) {
00159 *r__ = f;
00160 }
00161 tst1 = *r__;
00162 tst2 = tst1 + f;
00163 if (tst2 <= tst1) {
00164 goto L360;
00165 }
00166 if (its == 30) {
00167 goto L1000;
00168 }
00169 ++its;
00170 /* .......... FORM SHIFT FROM BOTTOM 2 BY 2 MINOR .......... */
00171 if (f > *r__ * .25 && its < 5) {
00172 goto L90;
00173 }
00174 f = a[*n + (*mb - 1) * a_dim1];
00175 if (f == 0.) {
00176 goto L70;
00177 }
00178 q = (a[*n - 1 + *mb * a_dim1] - g) / (f * 2.);
00179 s = pythag_(&q, &c_b8);
00180 g -= f / (q + d_sign(&s, &q));
00181 L70:
00182 *t += g;
00183
00184 i__1 = *n;
00185 for (i__ = 1; i__ <= i__1; ++i__) {
00186 /* L80: */
00187 a[i__ + *mb * a_dim1] -= g;
00188 }
00189
00190 L90:
00191 i__1 = m4;
00192 for (k = m31; k <= i__1; ++k) {
00193 /* L100: */
00194 rv[k] = 0.;
00195 }
00196
00197 i__1 = mn;
00198 for (ii = 1; ii <= i__1; ++ii) {
00199 i__ = ii - m;
00200 ni = *n - ii;
00201 if (ni < 0) {
00202 goto L230;
00203 }
00204 /* .......... FORM COLUMN OF SHIFTED MATRIX A-G*I .......... */
00205 /* Computing MAX */
00206 i__2 = 1, i__3 = 2 - i__;
00207 l = max(i__2,i__3);
00208
00209 i__2 = m3;
00210 for (k = 1; k <= i__2; ++k) {
00211 /* L110: */
00212 rv[k] = 0.;
00213 }
00214
00215 i__2 = m1;
00216 for (k = l; k <= i__2; ++k) {
00217 km = k + m;
00218 mk = k + mz;
00219 rv[km] = a[ii + mk * a_dim1];
00220 /* L120: */
00221 }
00222
00223 ll = min(m,ni);
00224 if (ll == 0) {
00225 goto L135;
00226 }
00227
00228 i__2 = ll;
00229 for (k = 1; k <= i__2; ++k) {
00230 km = k + m21;
00231 ik = ii + k;
00232 mk = *mb - k;
00233 rv[km] = a[ik + mk * a_dim1];
00234 /* L130: */
00235 }
00236 /* .......... PRE-MULTIPLY WITH HOUSEHOLDER REFLECTIONS ..........
00237 */
00238 L135:
00239 ll = m2;
00240 imult = 0;
00241 /* .......... MULTIPLICATION PROCEDURE .......... */
00242 L140:
00243 kj = m4 - m1;
00244
00245 i__2 = ll;
00246 for (j = 1; j <= i__2; ++j) {
00247 kj += m1;
00248 jm = j + m3;
00249 if (rv[jm] == 0.) {
00250 goto L170;
00251 }
00252 f = 0.;
00253
00254 i__3 = m1;
00255 for (k = 1; k <= i__3; ++k) {
00256 ++kj;
00257 jk = j + k - 1;
00258 f += rv[kj] * rv[jk];
00259 /* L150: */
00260 }
00261
00262 f /= rv[jm];
00263 kj -= m1;
00264
00265 i__3 = m1;
00266 for (k = 1; k <= i__3; ++k) {
00267 ++kj;
00268 jk = j + k - 1;
00269 rv[jk] -= rv[kj] * f;
00270 /* L160: */
00271 }
00272
00273 kj -= m1;
00274 L170:
00275 ;
00276 }
00277
00278 if (imult != 0) {
00279 goto L280;
00280 }
00281 /* .......... HOUSEHOLDER REFLECTION .......... */
00282 f = rv[m21];
00283 s = 0.;
00284 rv[m4] = 0.;
00285 scale = 0.;
00286
00287 i__2 = m3;
00288 for (k = m21; k <= i__2; ++k) {
00289 /* L180: */
00290 scale += (d__1 = rv[k], abs(d__1));
00291 }
00292
00293 if (scale == 0.) {
00294 goto L210;
00295 }
00296
00297 i__2 = m3;
00298 for (k = m21; k <= i__2; ++k) {
00299 /* L190: */
00300 /* Computing 2nd power */
00301 d__1 = rv[k] / scale;
00302 s += d__1 * d__1;
00303 }
00304
00305 s = scale * scale * s;
00306 d__1 = sqrt(s);
00307 g = -d_sign(&d__1, &f);
00308 rv[m21] = g;
00309 rv[m4] = s - f * g;
00310 kj = m4 + m2 * m1 + 1;
00311 rv[kj] = f - g;
00312
00313 i__2 = m1;
00314 for (k = 2; k <= i__2; ++k) {
00315 ++kj;
00316 km = k + m2;
00317 rv[kj] = rv[km];
00318 /* L200: */
00319 }
00320 /* .......... SAVE COLUMN OF TRIANGULAR FACTOR R .......... */
00321 L210:
00322 i__2 = m1;
00323 for (k = l; k <= i__2; ++k) {
00324 km = k + m;
00325 mk = k + mz;
00326 a[ii + mk * a_dim1] = rv[km];
00327 /* L220: */
00328 }
00329
00330 L230:
00331 /* Computing MAX */
00332 i__2 = 1, i__3 = m1 + 1 - i__;
00333 l = max(i__2,i__3);
00334 if (i__ <= 0) {
00335 goto L300;
00336 }
00337 /* .......... PERFORM ADDITIONAL STEPS .......... */
00338 i__2 = m21;
00339 for (k = 1; k <= i__2; ++k) {
00340 /* L240: */
00341 rv[k] = 0.;
00342 }
00343
00344 /* Computing MIN */
00345 i__2 = m1, i__3 = ni + m1;
00346 ll = min(i__2,i__3);
00347 /* .......... GET ROW OF TRIANGULAR FACTOR R .......... */
00348 i__2 = ll;
00349 for (kk = 1; kk <= i__2; ++kk) {
00350 k = kk - 1;
00351 km = k + m1;
00352 ik = i__ + k;
00353 mk = *mb - k;
00354 rv[km] = a[ik + mk * a_dim1];
00355 /* L250: */
00356 }
00357 /* .......... POST-MULTIPLY WITH HOUSEHOLDER REFLECTIONS .........
00358 . */
00359 ll = m1;
00360 imult = 1;
00361 goto L140;
00362 /* .......... STORE COLUMN OF NEW A MATRIX .......... */
00363 L280:
00364 i__2 = m1;
00365 for (k = l; k <= i__2; ++k) {
00366 mk = k + mz;
00367 a[i__ + mk * a_dim1] = rv[k];
00368 /* L290: */
00369 }
00370 /* .......... UPDATE HOUSEHOLDER REFLECTIONS .......... */
00371 L300:
00372 if (l > 1) {
00373 --l;
00374 }
00375 kj1 = m4 + l * m1;
00376
00377 i__2 = m2;
00378 for (j = l; j <= i__2; ++j) {
00379 jm = j + m3;
00380 rv[jm] = rv[jm + 1];
00381
00382 i__3 = m1;
00383 for (k = 1; k <= i__3; ++k) {
00384 ++kj1;
00385 kj = kj1 - m1;
00386 rv[kj] = rv[kj1];
00387 /* L320: */
00388 }
00389 }
00390
00391 /* L350: */
00392 }
00393
00394 goto L40;
00395 /* .......... CONVERGENCE .......... */
00396 L360:
00397 *t += g;
00398
00399 i__1 = *n;
00400 for (i__ = 1; i__ <= i__1; ++i__) {
00401 /* L380: */
00402 a[i__ + *mb * a_dim1] -= g;
00403 }
00404
00405 i__1 = m1;
00406 for (k = 1; k <= i__1; ++k) {
00407 mk = k + mz;
00408 a[*n + mk * a_dim1] = 0.;
00409 /* L400: */
00410 }
00411
00412 goto L1001;
00413 /* .......... SET ERROR -- NO CONVERGENCE TO */
00414 /* EIGENVALUE AFTER 30 ITERATIONS .......... */
00415 L1000:
00416 *ierr = *n;
00417 L1001:
00418 return 0;
00419 } /* bqr_ */
|
Variable Documentation
|
|
Definition at line 10 of file eis_bqr.c. Referenced by bqr_(). |
