/* bqr.f -- translated by f2c (version 19961017).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
#include "f2c.h"
/* Table of constant values */
static doublereal c_b8 = 1.;
/* Subroutine */ int bqr_(integer *nm, integer *n, integer *mb, doublereal *a,
doublereal *t, doublereal *r__, integer *ierr, integer *nv,
doublereal *rv)
{
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2, i__3;
doublereal d__1;
/* Builtin functions */
double d_sign(doublereal *, doublereal *), sqrt(doublereal);
/* Local variables */
static doublereal f, g;
static integer i__, j, k, l, m;
static doublereal q, s, scale;
static integer imult, m1, m2, m3, m4, m21, m31, ii, ik, jk, kj, jm, kk,
km, ll, mk, mn, ni, mz;
extern doublereal pythag_(doublereal *, doublereal *);
static integer kj1, its;
static doublereal tst1, tst2;
/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BQR, */
/* NUM. MATH. 16, 85-92(1970) BY MARTIN, REINSCH, AND WILKINSON. */
/* HANDBOOK FOR AUTO. COMP., VOL II-LINEAR ALGEBRA, 266-272(1971). */
/* THIS SUBROUTINE FINDS THE EIGENVALUE OF SMALLEST (USUALLY) */
/* MAGNITUDE OF A REAL SYMMETRIC BAND MATRIX USING THE */
/* QR ALGORITHM WITH SHIFTS OF ORIGIN. CONSECUTIVE CALLS */
/* CAN BE MADE TO FIND FURTHER EIGENVALUES. */
/* 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. */
/* MB IS THE (HALF) BAND WIDTH OF THE MATRIX, DEFINED AS THE */
/* NUMBER OF ADJACENT DIAGONALS, INCLUDING THE PRINCIPAL */
/* DIAGONAL, REQUIRED TO SPECIFY THE NON-ZERO PORTION OF THE */
/* LOWER TRIANGLE OF THE MATRIX. */
/* A CONTAINS THE LOWER TRIANGLE OF THE SYMMETRIC BAND INPUT */
/* MATRIX STORED AS AN N BY MB ARRAY. ITS LOWEST SUBDIAGONAL */
/* IS STORED IN THE LAST N+1-MB POSITIONS OF THE FIRST COLUMN, */
/* ITS NEXT SUBDIAGONAL IN THE LAST N+2-MB POSITIONS OF THE */
/* SECOND COLUMN, FURTHER SUBDIAGONALS SIMILARLY, AND FINALLY */
/* ITS PRINCIPAL DIAGONAL IN THE N POSITIONS OF THE LAST COLUMN.
*/
/* CONTENTS OF STORAGES NOT PART OF THE MATRIX ARE ARBITRARY. */
/* ON A SUBSEQUENT CALL, ITS OUTPUT CONTENTS FROM THE PREVIOUS */
/* CALL SHOULD BE PASSED. */
/* T SPECIFIES THE SHIFT (OF EIGENVALUES) APPLIED TO THE DIAGONAL
*/
/* OF A IN FORMING THE INPUT MATRIX. WHAT IS ACTUALLY DETERMINED
*/
/* IS THE EIGENVALUE OF A+TI (I IS THE IDENTITY MATRIX) NEAREST
*/
/* TO T. ON A SUBSEQUENT CALL, THE OUTPUT VALUE OF T FROM THE */
/* PREVIOUS CALL SHOULD BE PASSED IF THE NEXT NEAREST EIGENVALUE
*/
/* IS SOUGHT. */
/* R SHOULD BE SPECIFIED AS ZERO ON THE FIRST CALL, AND AS ITS */
/* OUTPUT VALUE FROM THE PREVIOUS CALL ON A SUBSEQUENT CALL. */
/* IT IS USED TO DETERMINE WHEN THE LAST ROW AND COLUMN OF */
/* THE TRANSFORMED BAND MATRIX CAN BE REGARDED AS NEGLIGIBLE. */
/* NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER RV */
/* AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT. */
/* ON OUTPUT */
/* A CONTAINS THE TRANSFORMED BAND MATRIX. THE MATRIX A+TI */
/* DERIVED FROM THE OUTPUT PARAMETERS IS SIMILAR TO THE */
/* INPUT A+TI TO WITHIN ROUNDING ERRORS. ITS LAST ROW AND */
/* COLUMN ARE NULL (IF IERR IS ZERO). */
/* T CONTAINS THE COMPUTED EIGENVALUE OF A+TI (IF IERR IS ZERO). */
/* R CONTAINS THE MAXIMUM OF ITS INPUT VALUE AND THE NORM OF THE */
/* LAST COLUMN OF THE INPUT MATRIX A. */
/* IERR IS SET TO */
/* ZERO FOR NORMAL RETURN, */
/* N IF THE EIGENVALUE HAS NOT BEEN */
/* DETERMINED AFTER 30 ITERATIONS. */
/* RV IS A TEMPORARY STORAGE ARRAY OF DIMENSION AT LEAST */
/* (2*MB**2+4*MB-3). THE FIRST (3*MB-2) LOCATIONS CORRESPOND */
/* TO THE ALGOL ARRAY B, THE NEXT (2*MB-1) LOCATIONS CORRESPOND
*/
/* TO THE ALGOL ARRAY H, AND THE FINAL (2*MB**2-MB) LOCATIONS */
/* CORRESPOND TO THE MB BY (2*MB-1) ALGOL ARRAY U. */
/* NOTE. FOR A SUBSEQUENT CALL, N SHOULD BE REPLACED BY N-1, BUT */
/* MB SHOULD NOT BE ALTERED EVEN WHEN IT EXCEEDS THE CURRENT N. */
/* 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 */
a_dim1 = *nm;
a_offset = a_dim1 + 1;
a -= a_offset;
--rv;
/* Function Body */
*ierr = 0;
m1 = min(*mb,*n);
m = m1 - 1;
m2 = m + m;
m21 = m2 + 1;
m3 = m21 + m;
m31 = m3 + 1;
m4 = m31 + m2;
mn = m + *n;
mz = *mb - m1;
its = 0;
/* .......... TEST FOR CONVERGENCE .......... */
L40:
g = a[*n + *mb * a_dim1];
if (m == 0) {
goto L360;
}
f = 0.;
i__1 = m;
for (k = 1; k <= i__1; ++k) {
mk = k + mz;
f += (d__1 = a[*n + mk * a_dim1], abs(d__1));
/* L50: */
}
if (its == 0 && f > *r__) {
*r__ = f;
}
tst1 = *r__;
tst2 = tst1 + f;
if (tst2 <= tst1) {
goto L360;
}
if (its == 30) {
goto L1000;
}
++its;
/* .......... FORM SHIFT FROM BOTTOM 2 BY 2 MINOR .......... */
if (f > *r__ * .25 && its < 5) {
goto L90;
}
f = a[*n + (*mb - 1) * a_dim1];
if (f == 0.) {
goto L70;
}
q = (a[*n - 1 + *mb * a_dim1] - g) / (f * 2.);
s = pythag_(&q, &c_b8);
g -= f / (q + d_sign(&s, &q));
L70:
*t += g;
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
/* L80: */
a[i__ + *mb * a_dim1] -= g;
}
L90:
i__1 = m4;
for (k = m31; k <= i__1; ++k) {
/* L100: */
rv[k] = 0.;
}
i__1 = mn;
for (ii = 1; ii <= i__1; ++ii) {
i__ = ii - m;
ni = *n - ii;
if (ni < 0) {
goto L230;
}
/* .......... FORM COLUMN OF SHIFTED MATRIX A-G*I .......... */
/* Computing MAX */
i__2 = 1, i__3 = 2 - i__;
l = max(i__2,i__3);
i__2 = m3;
for (k = 1; k <= i__2; ++k) {
/* L110: */
rv[k] = 0.;
}
i__2 = m1;
for (k = l; k <= i__2; ++k) {
km = k + m;
mk = k + mz;
rv[km] = a[ii + mk * a_dim1];
/* L120: */
}
ll = min(m,ni);
if (ll == 0) {
goto L135;
}
i__2 = ll;
for (k = 1; k <= i__2; ++k) {
km = k + m21;
ik = ii + k;
mk = *mb - k;
rv[km] = a[ik + mk * a_dim1];
/* L130: */
}
/* .......... PRE-MULTIPLY WITH HOUSEHOLDER REFLECTIONS ..........
*/
L135:
ll = m2;
imult = 0;
/* .......... MULTIPLICATION PROCEDURE .......... */
L140:
kj = m4 - m1;
i__2 = ll;
for (j = 1; j <= i__2; ++j) {
kj += m1;
jm = j + m3;
if (rv[jm] == 0.) {
goto L170;
}
f = 0.;
i__3 = m1;
for (k = 1; k <= i__3; ++k) {
++kj;
jk = j + k - 1;
f += rv[kj] * rv[jk];
/* L150: */
}
f /= rv[jm];
kj -= m1;
i__3 = m1;
for (k = 1; k <= i__3; ++k) {
++kj;
jk = j + k - 1;
rv[jk] -= rv[kj] * f;
/* L160: */
}
kj -= m1;
L170:
;
}
if (imult != 0) {
goto L280;
}
/* .......... HOUSEHOLDER REFLECTION .......... */
f = rv[m21];
s = 0.;
rv[m4] = 0.;
scale = 0.;
i__2 = m3;
for (k = m21; k <= i__2; ++k) {
/* L180: */
scale += (d__1 = rv[k], abs(d__1));
}
if (scale == 0.) {
goto L210;
}
i__2 = m3;
for (k = m21; k <= i__2; ++k) {
/* L190: */
/* Computing 2nd power */
d__1 = rv[k] / scale;
s += d__1 * d__1;
}
s = scale * scale * s;
d__1 = sqrt(s);
g = -d_sign(&d__1, &f);
rv[m21] = g;
rv[m4] = s - f * g;
kj = m4 + m2 * m1 + 1;
rv[kj] = f - g;
i__2 = m1;
for (k = 2; k <= i__2; ++k) {
++kj;
km = k + m2;
rv[kj] = rv[km];
/* L200: */
}
/* .......... SAVE COLUMN OF TRIANGULAR FACTOR R .......... */
L210:
i__2 = m1;
for (k = l; k <= i__2; ++k) {
km = k + m;
mk = k + mz;
a[ii + mk * a_dim1] = rv[km];
/* L220: */
}
L230:
/* Computing MAX */
i__2 = 1, i__3 = m1 + 1 - i__;
l = max(i__2,i__3);
if (i__ <= 0) {
goto L300;
}
/* .......... PERFORM ADDITIONAL STEPS .......... */
i__2 = m21;
for (k = 1; k <= i__2; ++k) {
/* L240: */
rv[k] = 0.;
}
/* Computing MIN */
i__2 = m1, i__3 = ni + m1;
ll = min(i__2,i__3);
/* .......... GET ROW OF TRIANGULAR FACTOR R .......... */
i__2 = ll;
for (kk = 1; kk <= i__2; ++kk) {
k = kk - 1;
km = k + m1;
ik = i__ + k;
mk = *mb - k;
rv[km] = a[ik + mk * a_dim1];
/* L250: */
}
/* .......... POST-MULTIPLY WITH HOUSEHOLDER REFLECTIONS .........
. */
ll = m1;
imult = 1;
goto L140;
/* .......... STORE COLUMN OF NEW A MATRIX .......... */
L280:
i__2 = m1;
for (k = l; k <= i__2; ++k) {
mk = k + mz;
a[i__ + mk * a_dim1] = rv[k];
/* L290: */
}
/* .......... UPDATE HOUSEHOLDER REFLECTIONS .......... */
L300:
if (l > 1) {
--l;
}
kj1 = m4 + l * m1;
i__2 = m2;
for (j = l; j <= i__2; ++j) {
jm = j + m3;
rv[jm] = rv[jm + 1];
i__3 = m1;
for (k = 1; k <= i__3; ++k) {
++kj1;
kj = kj1 - m1;
rv[kj] = rv[kj1];
/* L320: */
}
}
/* L350: */
}
goto L40;
/* .......... CONVERGENCE .......... */
L360:
*t += g;
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
/* L380: */
a[i__ + *mb * a_dim1] -= g;
}
i__1 = m1;
for (k = 1; k <= i__1; ++k) {
mk = k + mz;
a[*n + mk * a_dim1] = 0.;
/* L400: */
}
goto L1001;
/* .......... SET ERROR -- NO CONVERGENCE TO */
/* EIGENVALUE AFTER 30 ITERATIONS .......... */
L1000:
*ierr = *n;
L1001:
return 0;
} /* bqr_ */
syntax highlighted by Code2HTML, v. 0.9.1