00001 /* tred3.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 tred3_(integer *n, integer *nv, doublereal *a, 
00009         doublereal *d__, doublereal *e, doublereal *e2)
00010 {
00011     /* System generated locals */
00012     integer i__1, i__2, i__3;
00013     doublereal d__1;
00014 
00015     /* Builtin functions */
00016     double sqrt(doublereal), 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, jk, iz, jm1;
00023 
00024 
00025 
00026 /*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED3, */
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, STORED AS */
00031 /*     A ONE-DIMENSIONAL ARRAY, TO A SYMMETRIC TRIDIAGONAL MATRIX */
00032 /*     USING ORTHOGONAL SIMILARITY TRANSFORMATIONS. */
00033 
00034 /*     ON INPUT */
00035 
00036 /*        N IS THE ORDER OF THE MATRIX. */
00037 
00038 /*        NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A */
00039 /*          AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT. */
00040 
00041 /*        A CONTAINS THE LOWER TRIANGLE OF THE REAL SYMMETRIC */
00042 /*          INPUT MATRIX, STORED ROW-WISE AS A ONE-DIMENSIONAL */
00043 /*          ARRAY, IN ITS FIRST N*(N+1)/2 POSITIONS. */
00044 
00045 /*     ON OUTPUT */
00046 
00047 /*        A CONTAINS INFORMATION ABOUT THE ORTHOGONAL */
00048 /*          TRANSFORMATIONS USED IN THE REDUCTION. */
00049 
00050 /*        D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX. */
00051 
00052 /*        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL */
00053 /*          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO. */
00054 
00055 /*        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. */
00056 /*          E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED. */
00057 
00058 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00059 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00060 */
00061 
00062 /*     THIS VERSION DATED AUGUST 1983. */
00063 
00064 /*     ------------------------------------------------------------------ 
00065 */
00066 
00067 /*     .......... FOR I=N STEP -1 UNTIL 1 DO -- .......... */
00068     /* Parameter adjustments */
00069     --e2;
00070     --e;
00071     --d__;
00072     --a;
00073 
00074     /* Function Body */
00075     i__1 = *n;
00076     for (ii = 1; ii <= i__1; ++ii) {
00077         i__ = *n + 1 - ii;
00078         l = i__ - 1;
00079         iz = i__ * l / 2;
00080         h__ = 0.;
00081         scale = 0.;
00082         if (l < 1) {
00083             goto L130;
00084         }
00085 /*     .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) .......... */
00086         i__2 = l;
00087         for (k = 1; k <= i__2; ++k) {
00088             ++iz;
00089             d__[k] = a[iz];
00090             scale += (d__1 = d__[k], abs(d__1));
00091 /* L120: */
00092         }
00093 
00094         if (scale != 0.) {
00095             goto L140;
00096         }
00097 L130:
00098         e[i__] = 0.;
00099         e2[i__] = 0.;
00100         goto L290;
00101 
00102 L140:
00103         i__2 = l;
00104         for (k = 1; k <= i__2; ++k) {
00105             d__[k] /= scale;
00106             h__ += d__[k] * d__[k];
00107 /* L150: */
00108         }
00109 
00110         e2[i__] = scale * scale * h__;
00111         f = d__[l];
00112         d__1 = sqrt(h__);
00113         g = -d_sign(&d__1, &f);
00114         e[i__] = scale * g;
00115         h__ -= f * g;
00116         d__[l] = f - g;
00117         a[iz] = scale * d__[l];
00118         if (l == 1) {
00119             goto L290;
00120         }
00121         jk = 1;
00122 
00123         i__2 = l;
00124         for (j = 1; j <= i__2; ++j) {
00125             f = d__[j];
00126             g = 0.;
00127             jm1 = j - 1;
00128             if (jm1 < 1) {
00129                 goto L220;
00130             }
00131 
00132             i__3 = jm1;
00133             for (k = 1; k <= i__3; ++k) {
00134                 g += a[jk] * d__[k];
00135                 e[k] += a[jk] * f;
00136                 ++jk;
00137 /* L200: */
00138             }
00139 
00140 L220:
00141             e[j] = g + a[jk] * f;
00142             ++jk;
00143 /* L240: */
00144         }
00145 /*     .......... FORM P .......... */
00146         f = 0.;
00147 
00148         i__2 = l;
00149         for (j = 1; j <= i__2; ++j) {
00150             e[j] /= h__;
00151             f += e[j] * d__[j];
00152 /* L245: */
00153         }
00154 
00155         hh = f / (h__ + h__);
00156 /*     .......... FORM Q .......... */
00157         i__2 = l;
00158         for (j = 1; j <= i__2; ++j) {
00159 /* L250: */
00160             e[j] -= hh * d__[j];
00161         }
00162 
00163         jk = 1;
00164 /*     .......... FORM REDUCED A .......... */
00165         i__2 = l;
00166         for (j = 1; j <= i__2; ++j) {
00167             f = d__[j];
00168             g = e[j];
00169 
00170             i__3 = j;
00171             for (k = 1; k <= i__3; ++k) {
00172                 a[jk] = a[jk] - f * e[k] - g * d__[k];
00173                 ++jk;
00174 /* L260: */
00175             }
00176 
00177 /* L280: */
00178         }
00179 
00180 L290:
00181         d__[i__] = a[iz + 1];
00182         a[iz + 1] = scale * sqrt(h__);
00183 /* L300: */
00184     }
00185 
00186     return 0;
00187 } /* tred3_ */
00188