/* 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_ */
syntax highlighted by Code2HTML, v. 0.9.1