00001 /* bandr.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 /* Subroutine */ int bandr_(integer *nm, integer *n, integer *mb, doublereal *
00009         a, doublereal *d__, doublereal *e, doublereal *e2, logical *matz, 
00010         doublereal *z__)
00011 {
00012     /* System generated locals */
00013     integer a_dim1, a_offset, z_dim1, z_offset, i__1, i__2, i__3, i__4, i__5, 
00014             i__6;
00015     doublereal d__1;
00016 
00017     /* Builtin functions */
00018     double sqrt(doublereal);
00019 
00020     /* Local variables */
00021     static doublereal dmin__;
00022     static integer maxl, maxr;
00023     static doublereal g;
00024     static integer j, k, l, r__;
00025     static doublereal u, b1, b2, c2, f1, f2;
00026     static integer i1, i2, j1, j2, m1, n2, r1;
00027     static doublereal s2;
00028     static integer kr, mr;
00029     static doublereal dminrt;
00030     static integer ugl;
00031 
00032 
00033 
00034 /*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BANDRD, */
00035 /*     NUM. MATH. 12, 231-241(1968) BY SCHWARZ. */
00036 /*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 273-283(1971). */
00037 
00038 /*     THIS SUBROUTINE REDUCES A REAL SYMMETRIC BAND MATRIX */
00039 /*     TO A SYMMETRIC TRIDIAGONAL MATRIX USING AND OPTIONALLY */
00040 /*     ACCUMULATING ORTHOGONAL SIMILARITY TRANSFORMATIONS. */
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 /*        MB IS THE (HALF) BAND WIDTH OF THE MATRIX, DEFINED AS THE */
00051 /*          NUMBER OF ADJACENT DIAGONALS, INCLUDING THE PRINCIPAL */
00052 /*          DIAGONAL, REQUIRED TO SPECIFY THE NON-ZERO PORTION OF THE */
00053 /*          LOWER TRIANGLE OF THE MATRIX. */
00054 
00055 /*        A CONTAINS THE LOWER TRIANGLE OF THE SYMMETRIC BAND INPUT */
00056 /*          MATRIX STORED AS AN N BY MB ARRAY.  ITS LOWEST SUBDIAGONAL */
00057 /*          IS STORED IN THE LAST N+1-MB POSITIONS OF THE FIRST COLUMN, */
00058 /*          ITS NEXT SUBDIAGONAL IN THE LAST N+2-MB POSITIONS OF THE */
00059 /*          SECOND COLUMN, FURTHER SUBDIAGONALS SIMILARLY, AND FINALLY */
00060 /*          ITS PRINCIPAL DIAGONAL IN THE N POSITIONS OF THE LAST COLUMN. 
00061 */
00062 /*          CONTENTS OF STORAGES NOT PART OF THE MATRIX ARE ARBITRARY. */
00063 
00064 /*        MATZ SHOULD BE SET TO .TRUE. IF THE TRANSFORMATION MATRIX IS */
00065 /*          TO BE ACCUMULATED, AND TO .FALSE. OTHERWISE. */
00066 
00067 /*     ON OUTPUT */
00068 
00069 /*        A HAS BEEN DESTROYED, EXCEPT FOR ITS LAST TWO COLUMNS WHICH */
00070 /*          CONTAIN A COPY OF THE TRIDIAGONAL MATRIX. */
00071 
00072 /*        D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX. */
00073 
00074 /*        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL */
00075 /*          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO. */
00076 
00077 /*        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. */
00078 /*          E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED. */
00079 
00080 /*        Z CONTAINS THE ORTHOGONAL TRANSFORMATION MATRIX PRODUCED IN */
00081 /*          THE REDUCTION IF MATZ HAS BEEN SET TO .TRUE.  OTHERWISE, Z */
00082 /*          IS NOT REFERENCED. */
00083 
00084 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00085 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00086 */
00087 
00088 /*     THIS VERSION DATED AUGUST 1983. */
00089 
00090 /*     ------------------------------------------------------------------ 
00091 */
00092 
00093     /* Parameter adjustments */
00094     z_dim1 = *nm;
00095     z_offset = z_dim1 + 1;
00096     z__ -= z_offset;
00097     --e2;
00098     --e;
00099     --d__;
00100     a_dim1 = *nm;
00101     a_offset = a_dim1 + 1;
00102     a -= a_offset;
00103 
00104     /* Function Body */
00105     dmin__ = 5.4210108624275222e-20;
00106     dminrt = 2.3283064365386963e-10;
00107 /*     .......... INITIALIZE DIAGONAL SCALING MATRIX .......... */
00108     i__1 = *n;
00109     for (j = 1; j <= i__1; ++j) {
00110 /* L30: */
00111         d__[j] = 1.;
00112     }
00113 
00114     if (! (*matz)) {
00115         goto L60;
00116     }
00117 
00118     i__1 = *n;
00119     for (j = 1; j <= i__1; ++j) {
00120 
00121         i__2 = *n;
00122         for (k = 1; k <= i__2; ++k) {
00123 /* L40: */
00124             z__[j + k * z_dim1] = 0.;
00125         }
00126 
00127         z__[j + j * z_dim1] = 1.;
00128 /* L50: */
00129     }
00130 
00131 L60:
00132     m1 = *mb - 1;
00133     if ((i__1 = m1 - 1) < 0) {
00134         goto L900;
00135     } else if (i__1 == 0) {
00136         goto L800;
00137     } else {
00138         goto L70;
00139     }
00140 L70:
00141     n2 = *n - 2;
00142 
00143     i__1 = n2;
00144     for (k = 1; k <= i__1; ++k) {
00145 /* Computing MIN */
00146         i__2 = m1, i__3 = *n - k;
00147         maxr = min(i__2,i__3);
00148 /*     .......... FOR R=MAXR STEP -1 UNTIL 2 DO -- .......... */
00149         i__2 = maxr;
00150         for (r1 = 2; r1 <= i__2; ++r1) {
00151             r__ = maxr + 2 - r1;
00152             kr = k + r__;
00153             mr = *mb - r__;
00154             g = a[kr + mr * a_dim1];
00155             a[kr - 1 + a_dim1] = a[kr - 1 + (mr + 1) * a_dim1];
00156             ugl = k;
00157 
00158             i__3 = *n;
00159             i__4 = m1;
00160             for (j = kr; i__4 < 0 ? j >= i__3 : j <= i__3; j += i__4) {
00161                 j1 = j - 1;
00162                 j2 = j1 - 1;
00163                 if (g == 0.) {
00164                     goto L600;
00165                 }
00166                 b1 = a[j1 + a_dim1] / g;
00167                 b2 = b1 * d__[j1] / d__[j];
00168                 s2 = 1. / (b1 * b2 + 1.);
00169                 if (s2 >= .5) {
00170                     goto L450;
00171                 }
00172                 b1 = g / a[j1 + a_dim1];
00173                 b2 = b1 * d__[j] / d__[j1];
00174                 c2 = 1. - s2;
00175                 d__[j1] = c2 * d__[j1];
00176                 d__[j] = c2 * d__[j];
00177                 f1 = a[j + m1 * a_dim1] * 2.;
00178                 f2 = b1 * a[j1 + *mb * a_dim1];
00179                 a[j + m1 * a_dim1] = -b2 * (b1 * a[j + m1 * a_dim1] - a[j + *
00180                         mb * a_dim1]) - f2 + a[j + m1 * a_dim1];
00181                 a[j1 + *mb * a_dim1] = b2 * (b2 * a[j + *mb * a_dim1] + f1) + 
00182                         a[j1 + *mb * a_dim1];
00183                 a[j + *mb * a_dim1] = b1 * (f2 - f1) + a[j + *mb * a_dim1];
00184 
00185                 i__5 = j2;
00186                 for (l = ugl; l <= i__5; ++l) {
00187                     i2 = *mb - j + l;
00188                     u = a[j1 + (i2 + 1) * a_dim1] + b2 * a[j + i2 * a_dim1];
00189                     a[j + i2 * a_dim1] = -b1 * a[j1 + (i2 + 1) * a_dim1] + a[
00190                             j + i2 * a_dim1];
00191                     a[j1 + (i2 + 1) * a_dim1] = u;
00192 /* L200: */
00193                 }
00194 
00195                 ugl = j;
00196                 a[j1 + a_dim1] += b2 * g;
00197                 if (j == *n) {
00198                     goto L350;
00199                 }
00200 /* Computing MIN */
00201                 i__5 = m1, i__6 = *n - j1;
00202                 maxl = min(i__5,i__6);
00203 
00204                 i__5 = maxl;
00205                 for (l = 2; l <= i__5; ++l) {
00206                     i1 = j1 + l;
00207                     i2 = *mb - l;
00208                     u = a[i1 + i2 * a_dim1] + b2 * a[i1 + (i2 + 1) * a_dim1];
00209                     a[i1 + (i2 + 1) * a_dim1] = -b1 * a[i1 + i2 * a_dim1] + a[
00210                             i1 + (i2 + 1) * a_dim1];
00211                     a[i1 + i2 * a_dim1] = u;
00212 /* L300: */
00213                 }
00214 
00215                 i1 = j + m1;
00216                 if (i1 > *n) {
00217                     goto L350;
00218                 }
00219                 g = b2 * a[i1 + a_dim1];
00220 L350:
00221                 if (! (*matz)) {
00222                     goto L500;
00223                 }
00224 
00225                 i__5 = *n;
00226                 for (l = 1; l <= i__5; ++l) {
00227                     u = z__[l + j1 * z_dim1] + b2 * z__[l + j * z_dim1];
00228                     z__[l + j * z_dim1] = -b1 * z__[l + j1 * z_dim1] + z__[l 
00229                             + j * z_dim1];
00230                     z__[l + j1 * z_dim1] = u;
00231 /* L400: */
00232                 }
00233 
00234                 goto L500;
00235 
00236 L450:
00237                 u = d__[j1];
00238                 d__[j1] = s2 * d__[j];
00239                 d__[j] = s2 * u;
00240                 f1 = a[j + m1 * a_dim1] * 2.;
00241                 f2 = b1 * a[j + *mb * a_dim1];
00242                 u = b1 * (f2 - f1) + a[j1 + *mb * a_dim1];
00243                 a[j + m1 * a_dim1] = b2 * (b1 * a[j + m1 * a_dim1] - a[j1 + *
00244                         mb * a_dim1]) + f2 - a[j + m1 * a_dim1];
00245                 a[j1 + *mb * a_dim1] = b2 * (b2 * a[j1 + *mb * a_dim1] + f1) 
00246                         + a[j + *mb * a_dim1];
00247                 a[j + *mb * a_dim1] = u;
00248 
00249                 i__5 = j2;
00250                 for (l = ugl; l <= i__5; ++l) {
00251                     i2 = *mb - j + l;
00252                     u = b2 * a[j1 + (i2 + 1) * a_dim1] + a[j + i2 * a_dim1];
00253                     a[j + i2 * a_dim1] = -a[j1 + (i2 + 1) * a_dim1] + b1 * a[
00254                             j + i2 * a_dim1];
00255                     a[j1 + (i2 + 1) * a_dim1] = u;
00256 /* L460: */
00257                 }
00258 
00259                 ugl = j;
00260                 a[j1 + a_dim1] = b2 * a[j1 + a_dim1] + g;
00261                 if (j == *n) {
00262                     goto L480;
00263                 }
00264 /* Computing MIN */
00265                 i__5 = m1, i__6 = *n - j1;
00266                 maxl = min(i__5,i__6);
00267 
00268                 i__5 = maxl;
00269                 for (l = 2; l <= i__5; ++l) {
00270                     i1 = j1 + l;
00271                     i2 = *mb - l;
00272                     u = b2 * a[i1 + i2 * a_dim1] + a[i1 + (i2 + 1) * a_dim1];
00273                     a[i1 + (i2 + 1) * a_dim1] = -a[i1 + i2 * a_dim1] + b1 * a[
00274                             i1 + (i2 + 1) * a_dim1];
00275                     a[i1 + i2 * a_dim1] = u;
00276 /* L470: */
00277                 }
00278 
00279                 i1 = j + m1;
00280                 if (i1 > *n) {
00281                     goto L480;
00282                 }
00283                 g = a[i1 + a_dim1];
00284                 a[i1 + a_dim1] = b1 * a[i1 + a_dim1];
00285 L480:
00286                 if (! (*matz)) {
00287                     goto L500;
00288                 }
00289 
00290                 i__5 = *n;
00291                 for (l = 1; l <= i__5; ++l) {
00292                     u = b2 * z__[l + j1 * z_dim1] + z__[l + j * z_dim1];
00293                     z__[l + j * z_dim1] = -z__[l + j1 * z_dim1] + b1 * z__[l 
00294                             + j * z_dim1];
00295                     z__[l + j1 * z_dim1] = u;
00296 /* L490: */
00297                 }
00298 
00299 L500:
00300                 ;
00301             }
00302 
00303 L600:
00304             ;
00305         }
00306 
00307         if (k % 64 != 0) {
00308             goto L700;
00309         }
00310 /*     .......... RESCALE TO AVOID UNDERFLOW OR OVERFLOW .......... */
00311         i__2 = *n;
00312         for (j = k; j <= i__2; ++j) {
00313             if (d__[j] >= dmin__) {
00314                 goto L650;
00315             }
00316 /* Computing MAX */
00317             i__4 = 1, i__3 = *mb + 1 - j;
00318             maxl = max(i__4,i__3);
00319 
00320             i__4 = m1;
00321             for (l = maxl; l <= i__4; ++l) {
00322 /* L610: */
00323                 a[j + l * a_dim1] = dminrt * a[j + l * a_dim1];
00324             }
00325 
00326             if (j == *n) {
00327                 goto L630;
00328             }
00329 /* Computing MIN */
00330             i__4 = m1, i__3 = *n - j;
00331             maxl = min(i__4,i__3);
00332 
00333             i__4 = maxl;
00334             for (l = 1; l <= i__4; ++l) {
00335                 i1 = j + l;
00336                 i2 = *mb - l;
00337                 a[i1 + i2 * a_dim1] = dminrt * a[i1 + i2 * a_dim1];
00338 /* L620: */
00339             }
00340 
00341 L630:
00342             if (! (*matz)) {
00343                 goto L645;
00344             }
00345 
00346             i__4 = *n;
00347             for (l = 1; l <= i__4; ++l) {
00348 /* L640: */
00349                 z__[l + j * z_dim1] = dminrt * z__[l + j * z_dim1];
00350             }
00351 
00352 L645:
00353             a[j + *mb * a_dim1] = dmin__ * a[j + *mb * a_dim1];
00354             d__[j] /= dmin__;
00355 L650:
00356             ;
00357         }
00358 
00359 L700:
00360         ;
00361     }
00362 /*     .......... FORM SQUARE ROOT OF SCALING MATRIX .......... */
00363 L800:
00364     i__1 = *n;
00365     for (j = 2; j <= i__1; ++j) {
00366 /* L810: */
00367         e[j] = sqrt(d__[j]);
00368     }
00369 
00370     if (! (*matz)) {
00371         goto L840;
00372     }
00373 
00374     i__1 = *n;
00375     for (j = 1; j <= i__1; ++j) {
00376 
00377         i__2 = *n;
00378         for (k = 2; k <= i__2; ++k) {
00379 /* L820: */
00380             z__[j + k * z_dim1] = e[k] * z__[j + k * z_dim1];
00381         }
00382 
00383 /* L830: */
00384     }
00385 
00386 L840:
00387     u = 1.;
00388 
00389     i__1 = *n;
00390     for (j = 2; j <= i__1; ++j) {
00391         a[j + m1 * a_dim1] = u * e[j] * a[j + m1 * a_dim1];
00392         u = e[j];
00393 /* Computing 2nd power */
00394         d__1 = a[j + m1 * a_dim1];
00395         e2[j] = d__1 * d__1;
00396         a[j + *mb * a_dim1] = d__[j] * a[j + *mb * a_dim1];
00397         d__[j] = a[j + *mb * a_dim1];
00398         e[j] = a[j + m1 * a_dim1];
00399 /* L850: */
00400     }
00401 
00402     d__[1] = a[*mb * a_dim1 + 1];
00403     e[1] = 0.;
00404     e2[1] = 0.;
00405     goto L1001;
00406 
00407 L900:
00408     i__1 = *n;
00409     for (j = 1; j <= i__1; ++j) {
00410         d__[j] = a[j + *mb * a_dim1];
00411         e[j] = 0.;
00412         e2[j] = 0.;
00413 /* L950: */
00414     }
00415 
00416 L1001:
00417     return 0;
00418 } /* bandr_ */
00419