#include "blaswrap.h"
/* -- translated by f2c (version 19990503).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
#include "f2c.h"
/* Subroutine */ int orthes_(integer *nm, integer *n, integer *low, integer *
igh, real *a, real *ort)
{
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2, i__3;
real r__1;
/* Builtin functions */
double sqrt(doublereal), r_sign(real *, real *);
/* Local variables */
static real f, g, h__;
static integer i__, j, m;
static real scale;
static integer la, ii, jj, mp, kp1;
#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ORTHES,
NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON.
HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
GIVEN A REAL GENERAL MATRIX, THIS SUBROUTINE
REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS
LOW THROUGH IGH TO UPPER HESSENBERG FORM BY
ORTHOGONAL 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.
LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED,
SET LOW=1, IGH=N.
A CONTAINS THE INPUT MATRIX.
ON OUTPUT
A CONTAINS THE HESSENBERG MATRIX. INFORMATION ABOUT
THE ORTHOGONAL TRANSFORMATIONS USED IN THE REDUCTION
IS STORED IN THE REMAINING TRIANGLE UNDER THE
HESSENBERG MATRIX.
ORT CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS.
ONLY ELEMENTS LOW THROUGH IGH ARE USED.
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 */
a_dim1 = *nm;
a_offset = 1 + a_dim1 * 1;
a -= a_offset;
--ort;
/* Function Body */
la = *igh - 1;
kp1 = *low + 1;
if (la < kp1) {
goto L200;
}
i__1 = la;
for (m = kp1; m <= i__1; ++m) {
h__ = 0.f;
ort[m] = 0.f;
scale = 0.f;
/* .......... SCALE COLUMN (ALGOL TOL THEN NOT NEEDED) .......... */
i__2 = *igh;
for (i__ = m; i__ <= i__2; ++i__) {
/* L90: */
scale += (r__1 = a_ref(i__, m - 1), dabs(r__1));
}
if (scale == 0.f) {
goto L180;
}
mp = m + *igh;
/* .......... FOR I=IGH STEP -1 UNTIL M DO -- .......... */
i__2 = *igh;
for (ii = m; ii <= i__2; ++ii) {
i__ = mp - ii;
ort[i__] = a_ref(i__, m - 1) / scale;
h__ += ort[i__] * ort[i__];
/* L100: */
}
r__1 = sqrt(h__);
g = -r_sign(&r__1, &ort[m]);
h__ -= ort[m] * g;
ort[m] -= g;
/* .......... FORM (I-(U*UT)/H) * A .......... */
i__2 = *n;
for (j = m; j <= i__2; ++j) {
f = 0.f;
/* .......... FOR I=IGH STEP -1 UNTIL M DO -- .......... */
i__3 = *igh;
for (ii = m; ii <= i__3; ++ii) {
i__ = mp - ii;
f += ort[i__] * a_ref(i__, j);
/* L110: */
}
f /= h__;
i__3 = *igh;
for (i__ = m; i__ <= i__3; ++i__) {
/* L120: */
a_ref(i__, j) = a_ref(i__, j) - f * ort[i__];
}
/* L130: */
}
/* .......... FORM (I-(U*UT)/H)*A*(I-(U*UT)/H) .......... */
i__2 = *igh;
for (i__ = 1; i__ <= i__2; ++i__) {
f = 0.f;
/* .......... FOR J=IGH STEP -1 UNTIL M DO -- .......... */
i__3 = *igh;
for (jj = m; jj <= i__3; ++jj) {
j = mp - jj;
f += ort[j] * a_ref(i__, j);
/* L140: */
}
f /= h__;
i__3 = *igh;
for (j = m; j <= i__3; ++j) {
/* L150: */
a_ref(i__, j) = a_ref(i__, j) - f * ort[j];
}
/* L160: */
}
ort[m] = scale * ort[m];
a_ref(m, m - 1) = scale * g;
L180:
;
}
L200:
return 0;
} /* orthes_ */
#undef a_ref
/* Subroutine */ int tred1_(integer *nm, integer *n, real *a, real *d__, real
*e, real *e2)
{
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2, i__3;
real r__1;
/* Builtin functions */
double sqrt(doublereal), r_sign(real *, real *);
/* Local variables */
static real f, g, h__;
static integer i__, j, k, l;
static real scale;
static integer ii, jp1;
#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
/* THIS SUBROUTINE IS A TRANSLATION 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 REAL SYMMETRIC MATRIX
TO A SYMMETRIC TRIDIAGONAL MATRIX USING
ORTHOGONAL 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 REAL SYMMETRIC INPUT MATRIX. ONLY THE
LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED.
ON OUTPUT
A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANS-
FORMATIONS USED IN THE REDUCTION IN ITS STRICT LOWER
TRIANGLE. THE FULL UPPER TRIANGLE OF A IS UNALTERED.
D CONTAINS THE DIAGONAL ELEMENTS OF 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.
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 */
--e2;
--e;
--d__;
a_dim1 = *nm;
a_offset = 1 + a_dim1 * 1;
a -= a_offset;
/* Function Body */
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
d__[i__] = a_ref(*n, i__);
a_ref(*n, i__) = a_ref(i__, i__);
/* L100: */
}
/* .......... 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.f;
scale = 0.f;
if (l < 1) {
goto L130;
}
/* .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) .......... */
i__2 = l;
for (k = 1; k <= i__2; ++k) {
/* L120: */
scale += (r__1 = d__[k], dabs(r__1));
}
if (scale != 0.f) {
goto L140;
}
i__2 = l;
for (j = 1; j <= i__2; ++j) {
d__[j] = a_ref(l, j);
a_ref(l, j) = a_ref(i__, j);
a_ref(i__, j) = 0.f;
/* L125: */
}
L130:
e[i__] = 0.f;
e2[i__] = 0.f;
goto L300;
L140:
i__2 = l;
for (k = 1; k <= i__2; ++k) {
d__[k] /= scale;
h__ += d__[k] * d__[k];
/* L150: */
}
e2[i__] = scale * scale * h__;
f = d__[l];
r__1 = sqrt(h__);
g = -r_sign(&r__1, &f);
e[i__] = scale * g;
h__ -= f * g;
d__[l] = f - g;
if (l == 1) {
goto L285;
}
/* .......... FORM A*U .......... */
i__2 = l;
for (j = 1; j <= i__2; ++j) {
/* L170: */
e[j] = 0.f;
}
i__2 = l;
for (j = 1; j <= i__2; ++j) {
f = d__[j];
g = e[j] + a_ref(j, j) * f;
jp1 = j + 1;
if (l < jp1) {
goto L220;
}
i__3 = l;
for (k = jp1; k <= i__3; ++k) {
g += a_ref(k, j) * d__[k];
e[k] += a_ref(k, j) * f;
/* L200: */
}
L220:
e[j] = g;
/* L240: */
}
/* .......... FORM P .......... */
f = 0.f;
i__2 = l;
for (j = 1; j <= i__2; ++j) {
e[j] /= h__;
f += e[j] * d__[j];
/* L245: */
}
h__ = f / (h__ + h__);
/* .......... FORM Q .......... */
i__2 = l;
for (j = 1; j <= i__2; ++j) {
/* L250: */
e[j] -= h__ * d__[j];
}
/* .......... FORM REDUCED A .......... */
i__2 = l;
for (j = 1; j <= i__2; ++j) {
f = d__[j];
g = e[j];
i__3 = l;
for (k = j; k <= i__3; ++k) {
/* L260: */
a_ref(k, j) = a_ref(k, j) - f * e[k] - g * d__[k];
}
/* L280: */
}
L285:
i__2 = l;
for (j = 1; j <= i__2; ++j) {
f = d__[j];
d__[j] = a_ref(l, j);
a_ref(l, j) = a_ref(i__, j);
a_ref(i__, j) = f * scale;
/* L290: */
}
L300:
;
}
return 0;
} /* tred1_ */
#undef a_ref