00001 /* tred1.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 tred1_(integer *nm, integer *n, doublereal *a, 
00009         doublereal *d__, doublereal *e, doublereal *e2)
00010 {
00011     /* System generated locals */
00012     integer a_dim1, a_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;
00022     static integer ii, jp1;
00023 
00024 
00025 
00026 /*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED1, */
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 */
00031 /*     TO A SYMMETRIC TRIDIAGONAL MATRIX USING */
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 /*        A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANS- */
00048 /*          FORMATIONS USED IN THE REDUCTION IN ITS STRICT LOWER */
00049 /*          TRIANGLE.  THE FULL UPPER TRIANGLE OF A IS UNALTERED. */
00050 
00051 /*        D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX. */
00052 
00053 /*        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL */
00054 /*          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO. */
00055 
00056 /*        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. */
00057 /*          E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED. */
00058 
00059 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00060 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00061 */
00062 
00063 /*     THIS VERSION DATED AUGUST 1983. */
00064 
00065 /*     ------------------------------------------------------------------ 
00066 */
00067 
00068     /* Parameter adjustments */
00069     --e2;
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         d__[i__] = a[*n + i__ * a_dim1];
00080         a[*n + i__ * a_dim1] = a[i__ + i__ * a_dim1];
00081 /* L100: */
00082     }
00083 /*     .......... FOR I=N STEP -1 UNTIL 1 DO -- .......... */
00084     i__1 = *n;
00085     for (ii = 1; ii <= i__1; ++ii) {
00086         i__ = *n + 1 - ii;
00087         l = i__ - 1;
00088         h__ = 0.;
00089         scale = 0.;
00090         if (l < 1) {
00091             goto L130;
00092         }
00093 /*     .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) .......... */
00094         i__2 = l;
00095         for (k = 1; k <= i__2; ++k) {
00096 /* L120: */
00097             scale += (d__1 = d__[k], abs(d__1));
00098         }
00099 
00100         if (scale != 0.) {
00101             goto L140;
00102         }
00103 
00104         i__2 = l;
00105         for (j = 1; j <= i__2; ++j) {
00106             d__[j] = a[l + j * a_dim1];
00107             a[l + j * a_dim1] = a[i__ + j * a_dim1];
00108             a[i__ + j * a_dim1] = 0.;
00109 /* L125: */
00110         }
00111 
00112 L130:
00113         e[i__] = 0.;
00114         e2[i__] = 0.;
00115         goto L300;
00116 
00117 L140:
00118         i__2 = l;
00119         for (k = 1; k <= i__2; ++k) {
00120             d__[k] /= scale;
00121             h__ += d__[k] * d__[k];
00122 /* L150: */
00123         }
00124 
00125         e2[i__] = scale * scale * h__;
00126         f = d__[l];
00127         d__1 = sqrt(h__);
00128         g = -d_sign(&d__1, &f);
00129         e[i__] = scale * g;
00130         h__ -= f * g;
00131         d__[l] = f - g;
00132         if (l == 1) {
00133             goto L285;
00134         }
00135 /*     .......... FORM A*U .......... */
00136         i__2 = l;
00137         for (j = 1; j <= i__2; ++j) {
00138 /* L170: */
00139             e[j] = 0.;
00140         }
00141 
00142         i__2 = l;
00143         for (j = 1; j <= i__2; ++j) {
00144             f = d__[j];
00145             g = e[j] + a[j + j * a_dim1] * f;
00146             jp1 = j + 1;
00147             if (l < jp1) {
00148                 goto L220;
00149             }
00150 
00151             i__3 = l;
00152             for (k = jp1; k <= i__3; ++k) {
00153                 g += a[k + j * a_dim1] * d__[k];
00154                 e[k] += a[k + j * a_dim1] * f;
00155 /* L200: */
00156             }
00157 
00158 L220:
00159             e[j] = g;
00160 /* L240: */
00161         }
00162 /*     .......... FORM P .......... */
00163         f = 0.;
00164 
00165         i__2 = l;
00166         for (j = 1; j <= i__2; ++j) {
00167             e[j] /= h__;
00168             f += e[j] * d__[j];
00169 /* L245: */
00170         }
00171 
00172         h__ = f / (h__ + h__);
00173 /*     .......... FORM Q .......... */
00174         i__2 = l;
00175         for (j = 1; j <= i__2; ++j) {
00176 /* L250: */
00177             e[j] -= h__ * d__[j];
00178         }
00179 /*     .......... FORM REDUCED A .......... */
00180         i__2 = l;
00181         for (j = 1; j <= i__2; ++j) {
00182             f = d__[j];
00183             g = e[j];
00184 
00185             i__3 = l;
00186             for (k = j; k <= i__3; ++k) {
00187 /* L260: */
00188                 a[k + j * a_dim1] = a[k + j * a_dim1] - f * e[k] - g * d__[k];
00189             }
00190 
00191 /* L280: */
00192         }
00193 
00194 L285:
00195         i__2 = l;
00196         for (j = 1; j <= i__2; ++j) {
00197             f = d__[j];
00198             d__[j] = a[l + j * a_dim1];
00199             a[l + j * a_dim1] = a[i__ + j * a_dim1];
00200             a[i__ + j * a_dim1] = f * scale;
00201 /* L290: */
00202         }
00203 
00204 L300:
00205         ;
00206     }
00207 
00208     return 0;
00209 } /* tred1_ */
00210