/* htrid3.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 htrid3_(integer *nm, integer *n, doublereal *a,
doublereal *d__, doublereal *e, doublereal *e2, doublereal *tau)
{
/* System generated locals */
integer a_dim1, a_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 jm1, jp1;
/* THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF */
/* THE ALGOL PROCEDURE TRED3, 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, STORED AS */
/* A SINGLE SQUARE ARRAY, 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. */
/* A CONTAINS THE LOWER TRIANGLE OF THE COMPLEX HERMITIAN INPUT */
/* MATRIX. THE REAL PARTS OF THE MATRIX ELEMENTS ARE STORED */
/* IN THE FULL LOWER TRIANGLE OF A, AND THE IMAGINARY PARTS */
/* ARE STORED IN THE TRANSPOSED POSITIONS OF THE STRICT UPPER */
/* TRIANGLE OF A. NO STORAGE IS REQUIRED FOR THE ZERO */
/* IMAGINARY PARTS OF THE DIAGONAL ELEMENTS. */
/* ON OUTPUT */
/* A CONTAINS INFORMATION ABOUT THE UNITARY TRANSFORMATIONS */
/* USED IN THE REDUCTION. */
/* 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__;
a_dim1 = *nm;
a_offset = a_dim1 + 1;
a -= a_offset;
/* Function Body */
tau[(*n << 1) + 1] = 1.;
tau[(*n << 1) + 2] = 0.;
/* .......... 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 = a[i__ + k * a_dim1], abs(d__1)) + (d__2 =
a[k + i__ * a_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) {
a[i__ + k * a_dim1] /= scale;
a[k + i__ * a_dim1] /= scale;
h__ = h__ + a[i__ + k * a_dim1] * a[i__ + k * a_dim1] + a[k + i__
* a_dim1] * a[k + i__ * a_dim1];
/* L150: */
}
e2[i__] = scale * scale * h__;
g = sqrt(h__);
e[i__] = scale * g;
f = pythag_(&a[i__ + l * a_dim1], &a[l + i__ * a_dim1]);
/* .......... FORM NEXT DIAGONAL ELEMENT OF MATRIX T .......... */
if (f == 0.) {
goto L160;
}
tau[(l << 1) + 1] = (a[l + i__ * a_dim1] * tau[(i__ << 1) + 2] - a[
i__ + l * a_dim1] * tau[(i__ << 1) + 1]) / f;
si = (a[i__ + l * a_dim1] * tau[(i__ << 1) + 2] + a[l + i__ * a_dim1]
* tau[(i__ << 1) + 1]) / f;
h__ += f * g;
g = g / f + 1.;
a[i__ + l * a_dim1] = g * a[i__ + l * a_dim1];
a[l + i__ * a_dim1] = g * a[l + i__ * a_dim1];
if (l == 1) {
goto L270;
}
goto L170;
L160:
tau[(l << 1) + 1] = -tau[(i__ << 1) + 1];
si = tau[(i__ << 1) + 2];
a[i__ + l * a_dim1] = g;
L170:
f = 0.;
i__2 = l;
for (j = 1; j <= i__2; ++j) {
g = 0.;
gi = 0.;
if (j == 1) {
goto L190;
}
jm1 = j - 1;
/* .......... FORM ELEMENT OF A*U .......... */
i__3 = jm1;
for (k = 1; k <= i__3; ++k) {
g = g + a[j + k * a_dim1] * a[i__ + k * a_dim1] + a[k + j *
a_dim1] * a[k + i__ * a_dim1];
gi = gi - a[j + k * a_dim1] * a[k + i__ * a_dim1] + a[k + j *
a_dim1] * a[i__ + k * a_dim1];
/* L180: */
}
L190:
g += a[j + j * a_dim1] * a[i__ + j * a_dim1];
gi -= a[j + j * a_dim1] * a[j + i__ * a_dim1];
jp1 = j + 1;
if (l < jp1) {
goto L220;
}
i__3 = l;
for (k = jp1; k <= i__3; ++k) {
g = g + a[k + j * a_dim1] * a[i__ + k * a_dim1] - a[j + k *
a_dim1] * a[k + i__ * a_dim1];
gi = gi - a[k + j * a_dim1] * a[k + i__ * a_dim1] - a[j + k *
a_dim1] * a[i__ + k * a_dim1];
/* L200: */
}
/* .......... FORM ELEMENT OF P .......... */
L220:
e[j] = g / h__;
tau[(j << 1) + 2] = gi / h__;
f = f + e[j] * a[i__ + j * a_dim1] - tau[(j << 1) + 2] * a[j +
i__ * a_dim1];
/* L240: */
}
hh = f / (h__ + h__);
/* .......... FORM REDUCED A .......... */
i__2 = l;
for (j = 1; j <= i__2; ++j) {
f = a[i__ + j * a_dim1];
g = e[j] - hh * f;
e[j] = g;
fi = -a[j + i__ * a_dim1];
gi = tau[(j << 1) + 2] - hh * fi;
tau[(j << 1) + 2] = -gi;
a[j + j * a_dim1] -= (f * g + fi * gi) * 2.;
if (j == 1) {
goto L260;
}
jm1 = j - 1;
i__3 = jm1;
for (k = 1; k <= i__3; ++k) {
a[j + k * a_dim1] = a[j + k * a_dim1] - f * e[k] - g * a[i__
+ k * a_dim1] + fi * tau[(k << 1) + 2] + gi * a[k +
i__ * a_dim1];
a[k + j * a_dim1] = a[k + j * a_dim1] - f * tau[(k << 1) + 2]
- g * a[k + i__ * a_dim1] - fi * e[k] - gi * a[i__ +
k * a_dim1];
/* L250: */
}
L260:
;
}
L270:
i__2 = l;
for (k = 1; k <= i__2; ++k) {
a[i__ + k * a_dim1] = scale * a[i__ + k * a_dim1];
a[k + i__ * a_dim1] = scale * a[k + i__ * a_dim1];
/* L280: */
}
tau[(l << 1) + 2] = -si;
L290:
d__[i__] = a[i__ + i__ * a_dim1];
a[i__ + i__ * a_dim1] = scale * sqrt(h__);
/* L300: */
}
return 0;
} /* htrid3_ */
syntax highlighted by Code2HTML, v. 0.9.1