/* htridi.f -- translated by f2c (version 19961017). You must link the resulting object file with the libraries: -lf2c -lm (in that order) */ #include "f2c.h" /* Subroutine */ int htridi_(integer *nm, integer *n, doublereal *ar, doublereal *ai, doublereal *d__, doublereal *e, doublereal *e2, doublereal *tau) { /* System generated locals */ integer ar_dim1, ar_offset, ai_dim1, ai_offset, i__1, i__2, i__3; doublereal d__1, d__2; /* Builtin functions */ double sqrt(doublereal); /* Local variables */ static doublereal f, g, h__; static integer i__, j, k, l; static doublereal scale, fi, gi, hh; static integer ii; static doublereal si; extern doublereal pythag_(doublereal *, doublereal *); static integer jp1; /* THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF */ /* THE ALGOL PROCEDURE TRED1, NUM. MATH. 11, 181-195(1968) */ /* BY MARTIN, REINSCH, AND WILKINSON. */ /* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). */ /* THIS SUBROUTINE REDUCES A COMPLEX HERMITIAN MATRIX */ /* TO A REAL SYMMETRIC TRIDIAGONAL MATRIX USING */ /* UNITARY SIMILARITY TRANSFORMATIONS. */ /* ON INPUT */ /* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */ /* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */ /* DIMENSION STATEMENT. */ /* N IS THE ORDER OF THE MATRIX. */ /* AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */ /* RESPECTIVELY, OF THE COMPLEX HERMITIAN INPUT MATRIX. */ /* ONLY THE LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED. */ /* ON OUTPUT */ /* AR AND AI CONTAIN INFORMATION ABOUT THE UNITARY TRANS- */ /* FORMATIONS USED IN THE REDUCTION IN THEIR FULL LOWER */ /* TRIANGLES. THEIR STRICT UPPER TRIANGLES AND THE */ /* DIAGONAL OF AR ARE UNALTERED. */ /* D CONTAINS THE DIAGONAL ELEMENTS OF THE THE TRIDIAGONAL MATRIX. */ /* E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL */ /* MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO. */ /* E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. */ /* E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED. */ /* TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS. */ /* CALLS PYTHAG FOR DSQRT(A*A + B*B) . */ /* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */ /* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY */ /* THIS VERSION DATED AUGUST 1983. */ /* ------------------------------------------------------------------ */ /* Parameter adjustments */ tau -= 3; --e2; --e; --d__; ai_dim1 = *nm; ai_offset = ai_dim1 + 1; ai -= ai_offset; ar_dim1 = *nm; ar_offset = ar_dim1 + 1; ar -= ar_offset; /* Function Body */ tau[(*n << 1) + 1] = 1.; tau[(*n << 1) + 2] = 0.; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { /* L100: */ d__[i__] = ar[i__ + i__ * ar_dim1]; } /* .......... FOR I=N STEP -1 UNTIL 1 DO -- .......... */ i__1 = *n; for (ii = 1; ii <= i__1; ++ii) { i__ = *n + 1 - ii; l = i__ - 1; h__ = 0.; scale = 0.; if (l < 1) { goto L130; } /* .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) .......... */ i__2 = l; for (k = 1; k <= i__2; ++k) { /* L120: */ scale = scale + (d__1 = ar[i__ + k * ar_dim1], abs(d__1)) + (d__2 = ai[i__ + k * ai_dim1], abs(d__2)); } if (scale != 0.) { goto L140; } tau[(l << 1) + 1] = 1.; tau[(l << 1) + 2] = 0.; L130: e[i__] = 0.; e2[i__] = 0.; goto L290; L140: i__2 = l; for (k = 1; k <= i__2; ++k) { ar[i__ + k * ar_dim1] /= scale; ai[i__ + k * ai_dim1] /= scale; h__ = h__ + ar[i__ + k * ar_dim1] * ar[i__ + k * ar_dim1] + ai[ i__ + k * ai_dim1] * ai[i__ + k * ai_dim1]; /* L150: */ } e2[i__] = scale * scale * h__; g = sqrt(h__); e[i__] = scale * g; f = pythag_(&ar[i__ + l * ar_dim1], &ai[i__ + l * ai_dim1]); /* .......... FORM NEXT DIAGONAL ELEMENT OF MATRIX T .......... */ if (f == 0.) { goto L160; } tau[(l << 1) + 1] = (ai[i__ + l * ai_dim1] * tau[(i__ << 1) + 2] - ar[ i__ + l * ar_dim1] * tau[(i__ << 1) + 1]) / f; si = (ar[i__ + l * ar_dim1] * tau[(i__ << 1) + 2] + ai[i__ + l * ai_dim1] * tau[(i__ << 1) + 1]) / f; h__ += f * g; g = g / f + 1.; ar[i__ + l * ar_dim1] = g * ar[i__ + l * ar_dim1]; ai[i__ + l * ai_dim1] = g * ai[i__ + l * ai_dim1]; if (l == 1) { goto L270; } goto L170; L160: tau[(l << 1) + 1] = -tau[(i__ << 1) + 1]; si = tau[(i__ << 1) + 2]; ar[i__ + l * ar_dim1] = g; L170: f = 0.; i__2 = l; for (j = 1; j <= i__2; ++j) { g = 0.; gi = 0.; /* .......... FORM ELEMENT OF A*U .......... */ i__3 = j; for (k = 1; k <= i__3; ++k) { g = g + ar[j + k * ar_dim1] * ar[i__ + k * ar_dim1] + ai[j + k * ai_dim1] * ai[i__ + k * ai_dim1]; gi = gi - ar[j + k * ar_dim1] * ai[i__ + k * ai_dim1] + ai[j + k * ai_dim1] * ar[i__ + k * ar_dim1]; /* L180: */ } jp1 = j + 1; if (l < jp1) { goto L220; } i__3 = l; for (k = jp1; k <= i__3; ++k) { g = g + ar[k + j * ar_dim1] * ar[i__ + k * ar_dim1] - ai[k + j * ai_dim1] * ai[i__ + k * ai_dim1]; gi = gi - ar[k + j * ar_dim1] * ai[i__ + k * ai_dim1] - ai[k + j * ai_dim1] * ar[i__ + k * ar_dim1]; /* L200: */ } /* .......... FORM ELEMENT OF P .......... */ L220: e[j] = g / h__; tau[(j << 1) + 2] = gi / h__; f = f + e[j] * ar[i__ + j * ar_dim1] - tau[(j << 1) + 2] * ai[i__ + j * ai_dim1]; /* L240: */ } hh = f / (h__ + h__); /* .......... FORM REDUCED A .......... */ i__2 = l; for (j = 1; j <= i__2; ++j) { f = ar[i__ + j * ar_dim1]; g = e[j] - hh * f; e[j] = g; fi = -ai[i__ + j * ai_dim1]; gi = tau[(j << 1) + 2] - hh * fi; tau[(j << 1) + 2] = -gi; i__3 = j; for (k = 1; k <= i__3; ++k) { ar[j + k * ar_dim1] = ar[j + k * ar_dim1] - f * e[k] - g * ar[ i__ + k * ar_dim1] + fi * tau[(k << 1) + 2] + gi * ai[ i__ + k * ai_dim1]; ai[j + k * ai_dim1] = ai[j + k * ai_dim1] - f * tau[(k << 1) + 2] - g * ai[i__ + k * ai_dim1] - fi * e[k] - gi * ar[i__ + k * ar_dim1]; /* L260: */ } } L270: i__3 = l; for (k = 1; k <= i__3; ++k) { ar[i__ + k * ar_dim1] = scale * ar[i__ + k * ar_dim1]; ai[i__ + k * ai_dim1] = scale * ai[i__ + k * ai_dim1]; /* L280: */ } tau[(l << 1) + 2] = -si; L290: hh = d__[i__]; d__[i__] = ar[i__ + i__ * ar_dim1]; ar[i__ + i__ * ar_dim1] = hh; ai[i__ + i__ * ai_dim1] = scale * sqrt(h__); /* L300: */ } return 0; } /* htridi_ */