00001 /* bqr.f -- translated by f2c (version 19961017).
00002    You must link the resulting object file with the libraries:
00003         -lf2c -lm   (in that order)
00004 */
00005 
00006 #include "f2c.h"
00007 
00008 /* Table of constant values */
00009 
00010 static doublereal c_b8 = 1.;
00011 
00012 /* Subroutine */ int bqr_(integer *nm, integer *n, integer *mb, doublereal *a,
00013          doublereal *t, doublereal *r__, integer *ierr, integer *nv, 
00014         doublereal *rv)
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_ */
00420