00001 /* htridi.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 htridi_(integer *nm, integer *n, doublereal *ar, 
00009         doublereal *ai, doublereal *d__, doublereal *e, doublereal *e2, 
00010         doublereal *tau)
00011 {
00012     /* System generated locals */
00013     integer ar_dim1, ar_offset, ai_dim1, ai_offset, i__1, i__2, i__3;
00014     doublereal d__1, d__2;
00015 
00016     /* Builtin functions */
00017     double sqrt(doublereal);
00018 
00019     /* Local variables */
00020     static doublereal f, g, h__;
00021     static integer i__, j, k, l;
00022     static doublereal scale, fi, gi, hh;
00023     static integer ii;
00024     static doublereal si;
00025     extern doublereal pythag_(doublereal *, doublereal *);
00026     static integer jp1;
00027 
00028 
00029 
00030 /*     THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF */
00031 /*     THE ALGOL PROCEDURE TRED1, NUM. MATH. 11, 181-195(1968) */
00032 /*     BY MARTIN, REINSCH, AND WILKINSON. */
00033 /*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). */
00034 
00035 /*     THIS SUBROUTINE REDUCES A COMPLEX HERMITIAN MATRIX */
00036 /*     TO A REAL SYMMETRIC TRIDIAGONAL MATRIX USING */
00037 /*     UNITARY SIMILARITY TRANSFORMATIONS. */
00038 
00039 /*     ON INPUT */
00040 
00041 /*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
00042 /*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
00043 /*          DIMENSION STATEMENT. */
00044 
00045 /*        N IS THE ORDER OF THE MATRIX. */
00046 
00047 /*        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */
00048 /*          RESPECTIVELY, OF THE COMPLEX HERMITIAN INPUT MATRIX. */
00049 /*          ONLY THE LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED. */
00050 
00051 /*     ON OUTPUT */
00052 
00053 /*        AR AND AI CONTAIN INFORMATION ABOUT THE UNITARY TRANS- */
00054 /*          FORMATIONS USED IN THE REDUCTION IN THEIR FULL LOWER */
00055 /*          TRIANGLES.  THEIR STRICT UPPER TRIANGLES AND THE */
00056 /*          DIAGONAL OF AR ARE UNALTERED. */
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     ai_dim1 = *nm;
00086     ai_offset = ai_dim1 + 1;
00087     ai -= ai_offset;
00088     ar_dim1 = *nm;
00089     ar_offset = ar_dim1 + 1;
00090     ar -= ar_offset;
00091 
00092     /* Function Body */
00093     tau[(*n << 1) + 1] = 1.;
00094     tau[(*n << 1) + 2] = 0.;
00095 
00096     i__1 = *n;
00097     for (i__ = 1; i__ <= i__1; ++i__) {
00098 /* L100: */
00099         d__[i__] = ar[i__ + i__ * ar_dim1];
00100     }
00101 /*     .......... FOR I=N STEP -1 UNTIL 1 DO -- .......... */
00102     i__1 = *n;
00103     for (ii = 1; ii <= i__1; ++ii) {
00104         i__ = *n + 1 - ii;
00105         l = i__ - 1;
00106         h__ = 0.;
00107         scale = 0.;
00108         if (l < 1) {
00109             goto L130;
00110         }
00111 /*     .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) .......... */
00112         i__2 = l;
00113         for (k = 1; k <= i__2; ++k) {
00114 /* L120: */
00115             scale = scale + (d__1 = ar[i__ + k * ar_dim1], abs(d__1)) + (d__2 
00116                     = ai[i__ + k * ai_dim1], abs(d__2));
00117         }
00118 
00119         if (scale != 0.) {
00120             goto L140;
00121         }
00122         tau[(l << 1) + 1] = 1.;
00123         tau[(l << 1) + 2] = 0.;
00124 L130:
00125         e[i__] = 0.;
00126         e2[i__] = 0.;
00127         goto L290;
00128 
00129 L140:
00130         i__2 = l;
00131         for (k = 1; k <= i__2; ++k) {
00132             ar[i__ + k * ar_dim1] /= scale;
00133             ai[i__ + k * ai_dim1] /= scale;
00134             h__ = h__ + ar[i__ + k * ar_dim1] * ar[i__ + k * ar_dim1] + ai[
00135                     i__ + k * ai_dim1] * ai[i__ + k * ai_dim1];
00136 /* L150: */
00137         }
00138 
00139         e2[i__] = scale * scale * h__;
00140         g = sqrt(h__);
00141         e[i__] = scale * g;
00142         f = pythag_(&ar[i__ + l * ar_dim1], &ai[i__ + l * ai_dim1]);
00143 /*     .......... FORM NEXT DIAGONAL ELEMENT OF MATRIX T .......... */
00144         if (f == 0.) {
00145             goto L160;
00146         }
00147         tau[(l << 1) + 1] = (ai[i__ + l * ai_dim1] * tau[(i__ << 1) + 2] - ar[
00148                 i__ + l * ar_dim1] * tau[(i__ << 1) + 1]) / f;
00149         si = (ar[i__ + l * ar_dim1] * tau[(i__ << 1) + 2] + ai[i__ + l * 
00150                 ai_dim1] * tau[(i__ << 1) + 1]) / f;
00151         h__ += f * g;
00152         g = g / f + 1.;
00153         ar[i__ + l * ar_dim1] = g * ar[i__ + l * ar_dim1];
00154         ai[i__ + l * ai_dim1] = g * ai[i__ + l * ai_dim1];
00155         if (l == 1) {
00156             goto L270;
00157         }
00158         goto L170;
00159 L160:
00160         tau[(l << 1) + 1] = -tau[(i__ << 1) + 1];
00161         si = tau[(i__ << 1) + 2];
00162         ar[i__ + l * ar_dim1] = g;
00163 L170:
00164         f = 0.;
00165 
00166         i__2 = l;
00167         for (j = 1; j <= i__2; ++j) {
00168             g = 0.;
00169             gi = 0.;
00170 /*     .......... FORM ELEMENT OF A*U .......... */
00171             i__3 = j;
00172             for (k = 1; k <= i__3; ++k) {
00173                 g = g + ar[j + k * ar_dim1] * ar[i__ + k * ar_dim1] + ai[j + 
00174                         k * ai_dim1] * ai[i__ + k * ai_dim1];
00175                 gi = gi - ar[j + k * ar_dim1] * ai[i__ + k * ai_dim1] + ai[j 
00176                         + k * ai_dim1] * ar[i__ + k * ar_dim1];
00177 /* L180: */
00178             }
00179 
00180             jp1 = j + 1;
00181             if (l < jp1) {
00182                 goto L220;
00183             }
00184 
00185             i__3 = l;
00186             for (k = jp1; k <= i__3; ++k) {
00187                 g = g + ar[k + j * ar_dim1] * ar[i__ + k * ar_dim1] - ai[k + 
00188                         j * ai_dim1] * ai[i__ + k * ai_dim1];
00189                 gi = gi - ar[k + j * ar_dim1] * ai[i__ + k * ai_dim1] - ai[k 
00190                         + j * ai_dim1] * ar[i__ + k * ar_dim1];
00191 /* L200: */
00192             }
00193 /*     .......... FORM ELEMENT OF P .......... */
00194 L220:
00195             e[j] = g / h__;
00196             tau[(j << 1) + 2] = gi / h__;
00197             f = f + e[j] * ar[i__ + j * ar_dim1] - tau[(j << 1) + 2] * ai[i__ 
00198                     + j * ai_dim1];
00199 /* L240: */
00200         }
00201 
00202         hh = f / (h__ + h__);
00203 /*     .......... FORM REDUCED A .......... */
00204         i__2 = l;
00205         for (j = 1; j <= i__2; ++j) {
00206             f = ar[i__ + j * ar_dim1];
00207             g = e[j] - hh * f;
00208             e[j] = g;
00209             fi = -ai[i__ + j * ai_dim1];
00210             gi = tau[(j << 1) + 2] - hh * fi;
00211             tau[(j << 1) + 2] = -gi;
00212 
00213             i__3 = j;
00214             for (k = 1; k <= i__3; ++k) {
00215                 ar[j + k * ar_dim1] = ar[j + k * ar_dim1] - f * e[k] - g * ar[
00216                         i__ + k * ar_dim1] + fi * tau[(k << 1) + 2] + gi * ai[
00217                         i__ + k * ai_dim1];
00218                 ai[j + k * ai_dim1] = ai[j + k * ai_dim1] - f * tau[(k << 1) 
00219                         + 2] - g * ai[i__ + k * ai_dim1] - fi * e[k] - gi * 
00220                         ar[i__ + k * ar_dim1];
00221 /* L260: */
00222             }
00223         }
00224 
00225 L270:
00226         i__3 = l;
00227         for (k = 1; k <= i__3; ++k) {
00228             ar[i__ + k * ar_dim1] = scale * ar[i__ + k * ar_dim1];
00229             ai[i__ + k * ai_dim1] = scale * ai[i__ + k * ai_dim1];
00230 /* L280: */
00231         }
00232 
00233         tau[(l << 1) + 2] = -si;
00234 L290:
00235         hh = d__[i__];
00236         d__[i__] = ar[i__ + i__ * ar_dim1];
00237         ar[i__ + i__ * ar_dim1] = hh;
00238         ai[i__ + i__ * ai_dim1] = scale * sqrt(h__);
00239 /* L300: */
00240     }
00241 
00242     return 0;
00243 } /* htridi_ */
00244