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