00010 {
00011     /* System generated locals */
00012     integer a_dim1, a_offset, z_dim1, z_offset, i__1, i__2, i__3;
00013     doublereal d__1;
00014 
00015     /* Builtin functions */
00016     double d_sign(doublereal *, doublereal *);
00017 
00018     /* Local variables */
00019     static doublereal f, g, h__;
00020     static integer i__, j, k, l;
00021     static doublereal scale, hh;
00022     static integer ii, jp1;
00023 
00024 
00025 
00026 /*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED2, */
00027 /*     NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON. */
00028 /*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). */
00029 
00030 /*     THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX TO A */
00031 /*     SYMMETRIC TRIDIAGONAL MATRIX USING AND ACCUMULATING */
00032 /*     ORTHOGONAL SIMILARITY TRANSFORMATIONS. */
00033 
00034 /*     ON INPUT */
00035 
00036 /*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
00037 /*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
00038 /*          DIMENSION STATEMENT. */
00039 
00040 /*        N IS THE ORDER OF THE MATRIX. */
00041 
00042 /*        A CONTAINS THE REAL SYMMETRIC INPUT MATRIX.  ONLY THE */
00043 /*          LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED. */
00044 
00045 /*     ON OUTPUT */
00046 
00047 /*        D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX. */
00048 
00049 /*        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL */
00050 /*          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO. */
00051 
00052 /*        Z CONTAINS THE ORTHOGONAL TRANSFORMATION MATRIX */
00053 /*          PRODUCED IN THE REDUCTION. */
00054 
00055 /*        A AND Z MAY COINCIDE.  IF DISTINCT, A IS UNALTERED. */
00056 
00057 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00058 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00059 */
00060 
00061 /*     THIS VERSION DATED AUGUST 1983. */
00062 
00063 /*     ------------------------------------------------------------------ 
00064 */
00065 
00066     /* Parameter adjustments */
00067     z_dim1 = *nm;
00068     z_offset = z_dim1 + 1;
00069     z__ -= z_offset;
00070     --e;
00071     --d__;
00072     a_dim1 = *nm;
00073     a_offset = a_dim1 + 1;
00074     a -= a_offset;
00075 
00076     /* Function Body */
00077     i__1 = *n;
00078     for (i__ = 1; i__ <= i__1; ++i__) {
00079 
00080         i__2 = *n;
00081         for (j = i__; j <= i__2; ++j) {
00082 /* L80: */
00083             z__[j + i__ * z_dim1] = a[j + i__ * a_dim1];
00084         }
00085 
00086         d__[i__] = a[*n + i__ * a_dim1];
00087 /* L100: */
00088     }
00089 
00090     if (*n == 1) {
00091         goto L510;
00092     }
00093 /*     .......... FOR I=N STEP -1 UNTIL 2 DO -- .......... */
00094     i__1 = *n;
00095     for (ii = 2; ii <= i__1; ++ii) {
00096         i__ = *n + 2 - ii;
00097         l = i__ - 1;
00098         h__ = 0.;
00099         scale = 0.;
00100         if (l < 2) {
00101             goto L130;
00102         }
00103 /*     .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) .......... */
00104         i__2 = l;
00105         for (k = 1; k <= i__2; ++k) {
00106 /* L120: */
00107             scale += (d__1 = d__[k], abs(d__1));
00108         }
00109 
00110         if (scale != 0.) {
00111             goto L140;
00112         }
00113 L130:
00114         e[i__] = d__[l];
00115 
00116         i__2 = l;
00117         for (j = 1; j <= i__2; ++j) {
00118             d__[j] = z__[l + j * z_dim1];
00119             z__[i__ + j * z_dim1] = 0.;
00120             z__[j + i__ * z_dim1] = 0.;
00121 /* L135: */
00122         }
00123 
00124         goto L290;
00125 
00126 L140:
00127         i__2 = l;
00128         for (k = 1; k <= i__2; ++k) {
00129             d__[k] /= scale;
00130             h__ += d__[k] * d__[k];
00131 /* L150: */
00132         }
00133 
00134         f = d__[l];
00135         d__1 = sqrt(h__);
00136         g = -d_sign(&d__1, &f);
00137         e[i__] = scale * g;
00138         h__ -= f * g;
00139         d__[l] = f - g;
00140 /*     .......... FORM A*U .......... */
00141         i__2 = l;
00142         for (j = 1; j <= i__2; ++j) {
00143 /* L170: */
00144             e[j] = 0.;
00145         }
00146 
00147         i__2 = l;
00148         for (j = 1; j <= i__2; ++j) {
00149             f = d__[j];
00150             z__[j + i__ * z_dim1] = f;
00151             g = e[j] + z__[j + j * z_dim1] * f;
00152             jp1 = j + 1;
00153             if (l < jp1) {
00154                 goto L220;
00155             }
00156 
00157             i__3 = l;
00158             for (k = jp1; k <= i__3; ++k) {
00159                 g += z__[k + j * z_dim1] * d__[k];
00160                 e[k] += z__[k + j * z_dim1] * f;
00161 /* L200: */
00162             }
00163 
00164 L220:
00165             e[j] = g;
00166 /* L240: */
00167         }
00168 /*     .......... FORM P .......... */
00169         f = 0.;
00170 
00171         i__2 = l;
00172         for (j = 1; j <= i__2; ++j) {
00173             e[j] /= h__;
00174             f += e[j] * d__[j];
00175 /* L245: */
00176         }
00177 
00178         hh = f / (h__ + h__);
00179 /*     .......... FORM Q .......... */
00180         i__2 = l;
00181         for (j = 1; j <= i__2; ++j) {
00182 /* L250: */
00183             e[j] -= hh * d__[j];
00184         }
00185 /*     .......... FORM REDUCED A .......... */
00186         i__2 = l;
00187         for (j = 1; j <= i__2; ++j) {
00188             f = d__[j];
00189             g = e[j];
00190 
00191             i__3 = l;
00192             for (k = j; k <= i__3; ++k) {
00193 /* L260: */
00194                 z__[k + j * z_dim1] = z__[k + j * z_dim1] - f * e[k] - g * 
00195                         d__[k];
00196             }
00197 
00198             d__[j] = z__[l + j * z_dim1];
00199             z__[i__ + j * z_dim1] = 0.;
00200 /* L280: */
00201         }
00202 
00203 L290:
00204         d__[i__] = h__;
00205 /* L300: */
00206     }
00207 /*     .......... ACCUMULATION OF TRANSFORMATION MATRICES .......... */
00208     i__1 = *n;
00209     for (i__ = 2; i__ <= i__1; ++i__) {
00210         l = i__ - 1;
00211         z__[*n + l * z_dim1] = z__[l + l * z_dim1];
00212         z__[l + l * z_dim1] = 1.;
00213         h__ = d__[i__];
00214         if (h__ == 0.) {
00215             goto L380;
00216         }
00217 
00218         i__2 = l;
00219         for (k = 1; k <= i__2; ++k) {
00220 /* L330: */
00221             d__[k] = z__[k + i__ * z_dim1] / h__;
00222         }
00223 
00224         i__2 = l;
00225         for (j = 1; j <= i__2; ++j) {
00226             g = 0.;
00227 
00228             i__3 = l;
00229             for (k = 1; k <= i__3; ++k) {
00230 /* L340: */
00231                 g += z__[k + i__ * z_dim1] * z__[k + j * z_dim1];
00232             }
00233 
00234             i__3 = l;
00235             for (k = 1; k <= i__3; ++k) {
00236                 z__[k + j * z_dim1] -= g * d__[k];
00237 /* L360: */
00238             }
00239         }
00240 
00241 L380:
00242         i__3 = l;
00243         for (k = 1; k <= i__3; ++k) {
00244 /* L400: */
00245             z__[k + i__ * z_dim1] = 0.;
00246         }
00247 
00248 /* L500: */
00249     }
00250 
00251 L510:
00252     i__1 = *n;
00253     for (i__ = 1; i__ <= i__1; ++i__) {
00254         d__[i__] = z__[*n + i__ * z_dim1];
00255         z__[*n + i__ * z_dim1] = 0.;
00256 /* L520: */
00257     }
00258 
00259     z__[*n + *n * z_dim1] = 1.;
00260     e[1] = 0.;
00261     return 0;
00262 } /* tred2_ */