/* ratqr.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 ratqr_(integer *n, doublereal *eps1, doublereal *d__,
doublereal *e, doublereal *e2, integer *m, doublereal *w, integer *
ind, doublereal *bd, logical *type__, integer *idef, integer *ierr)
{
/* System generated locals */
integer i__1, i__2;
doublereal d__1, d__2, d__3;
/* Local variables */
static integer jdef;
static doublereal f;
static integer i__, j, k;
static doublereal p, q, r__, s, delta;
static integer k1, ii, jj;
static doublereal ep, qp;
extern doublereal epslon_(doublereal *);
static doublereal err, tot;
/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE RATQR, */
/* NUM. MATH. 11, 264-272(1968) BY REINSCH AND BAUER. */
/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 257-265(1971). */
/* THIS SUBROUTINE FINDS THE ALGEBRAICALLY SMALLEST OR LARGEST */
/* EIGENVALUES OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE */
/* RATIONAL QR METHOD WITH NEWTON CORRECTIONS. */
/* ON INPUT */
/* N IS THE ORDER OF THE MATRIX. */
/* EPS1 IS A THEORETICAL ABSOLUTE ERROR TOLERANCE FOR THE */
/* COMPUTED EIGENVALUES. IF THE INPUT EPS1 IS NON-POSITIVE, */
/* OR INDEED SMALLER THAN ITS DEFAULT VALUE, IT IS RESET */
/* AT EACH ITERATION TO THE RESPECTIVE DEFAULT VALUE, */
/* NAMELY, THE PRODUCT OF THE RELATIVE MACHINE PRECISION */
/* AND THE MAGNITUDE OF THE CURRENT EIGENVALUE ITERATE. */
/* THE THEORETICAL ABSOLUTE ERROR IN THE K-TH EIGENVALUE */
/* IS USUALLY NOT GREATER THAN K TIMES EPS1. */
/* D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. */
/* E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX */
/* IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. */
/* E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. */
/* E2(1) IS ARBITRARY. */
/* M IS THE NUMBER OF EIGENVALUES TO BE FOUND. */
/* IDEF SHOULD BE SET TO 1 IF THE INPUT MATRIX IS KNOWN TO BE */
/* POSITIVE DEFINITE, TO -1 IF THE INPUT MATRIX IS KNOWN TO */
/* BE NEGATIVE DEFINITE, AND TO 0 OTHERWISE. */
/* TYPE SHOULD BE SET TO .TRUE. IF THE SMALLEST EIGENVALUES */
/* ARE TO BE FOUND, AND TO .FALSE. IF THE LARGEST EIGENVALUES */
/* ARE TO BE FOUND. */
/* ON OUTPUT */
/* EPS1 IS UNALTERED UNLESS IT HAS BEEN RESET TO ITS */
/* (LAST) DEFAULT VALUE. */
/* D AND E ARE UNALTERED (UNLESS W OVERWRITES D). */
/* ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED */
/* AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE */
/* MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES. */
/* E2(1) IS SET TO 0.0D0 IF THE SMALLEST EIGENVALUES HAVE BEEN */
/* FOUND, AND TO 2.0D0 IF THE LARGEST EIGENVALUES HAVE BEEN */
/* FOUND. E2 IS OTHERWISE UNALTERED (UNLESS OVERWRITTEN BY BD).
*/
/* W CONTAINS THE M ALGEBRAICALLY SMALLEST EIGENVALUES IN */
/* ASCENDING ORDER, OR THE M LARGEST EIGENVALUES IN */
/* DESCENDING ORDER. IF AN ERROR EXIT IS MADE BECAUSE OF */
/* AN INCORRECT SPECIFICATION OF IDEF, NO EIGENVALUES */
/* ARE FOUND. IF THE NEWTON ITERATES FOR A PARTICULAR */
/* EIGENVALUE ARE NOT MONOTONE, THE BEST ESTIMATE OBTAINED */
/* IS RETURNED AND IERR IS SET. W MAY COINCIDE WITH D. */
/* IND CONTAINS IN ITS FIRST M POSITIONS THE SUBMATRIX INDICES */
/* ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W -- */
/* 1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX FROM */
/* THE TOP, 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC..
*/
/* BD CONTAINS REFINED BOUNDS FOR THE THEORETICAL ERRORS OF THE */
/* CORRESPONDING EIGENVALUES IN W. THESE BOUNDS ARE USUALLY */
/* WITHIN THE TOLERANCE SPECIFIED BY EPS1. BD MAY COINCIDE */
/* WITH E2. */
/* IERR IS SET TO */
/* ZERO FOR NORMAL RETURN, */
/* 6*N+1 IF IDEF IS SET TO 1 AND TYPE TO .TRUE. */
/* WHEN THE MATRIX IS NOT POSITIVE DEFINITE, OR */
/* IF IDEF IS SET TO -1 AND TYPE TO .FALSE. */
/* WHEN THE MATRIX IS NOT NEGATIVE DEFINITE, */
/* 5*N+K IF SUCCESSIVE ITERATES TO THE K-TH EIGENVALUE */
/* ARE NOT MONOTONE INCREASING, WHERE K REFERS */
/* TO THE LAST SUCH OCCURRENCE. */
/* NOTE THAT SUBROUTINE TRIDIB IS GENERALLY FASTER AND MORE */
/* ACCURATE THAN RATQR IF THE EIGENVALUES ARE CLUSTERED. */
/* 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 */
--bd;
--ind;
--w;
--e2;
--e;
--d__;
/* Function Body */
*ierr = 0;
jdef = *idef;
/* .......... COPY D ARRAY INTO W .......... */
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
/* L20: */
w[i__] = d__[i__];
}
if (*type__) {
goto L40;
}
j = 1;
goto L400;
L40:
err = 0.;
s = 0.;
/* .......... LOOK FOR SMALL SUB-DIAGONAL ENTRIES AND DEFINE */
/* INITIAL SHIFT FROM LOWER GERSCHGORIN BOUND. */
/* COPY E2 ARRAY INTO BD .......... */
tot = w[1];
q = 0.;
j = 0;
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
p = q;
if (i__ == 1) {
goto L60;
}
d__3 = (d__1 = d__[i__], abs(d__1)) + (d__2 = d__[i__ - 1], abs(d__2))
;
if (p > epslon_(&d__3)) {
goto L80;
}
L60:
e2[i__] = 0.;
L80:
bd[i__] = e2[i__];
/* .......... COUNT ALSO IF ELEMENT OF E2 HAS UNDERFLOWED ........
.. */
if (e2[i__] == 0.) {
++j;
}
ind[i__] = j;
q = 0.;
if (i__ != *n) {
q = (d__1 = e[i__ + 1], abs(d__1));
}
/* Computing MIN */
d__1 = w[i__] - p - q;
tot = min(d__1,tot);
/* L100: */
}
if (jdef == 1 && tot < 0.) {
goto L140;
}
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
/* L110: */
w[i__] -= tot;
}
goto L160;
L140:
tot = 0.;
L160:
i__1 = *m;
for (k = 1; k <= i__1; ++k) {
/* .......... NEXT QR TRANSFORMATION .......... */
L180:
tot += s;
delta = w[*n] - s;
i__ = *n;
f = (d__1 = epslon_(&tot), abs(d__1));
if (*eps1 < f) {
*eps1 = f;
}
if (delta > *eps1) {
goto L190;
}
if (delta < -(*eps1)) {
goto L1000;
}
goto L300;
/* .......... REPLACE SMALL SUB-DIAGONAL SQUARES BY ZERO */
/* TO REDUCE THE INCIDENCE OF UNDERFLOWS .......... */
L190:
if (k == *n) {
goto L210;
}
k1 = k + 1;
i__2 = *n;
for (j = k1; j <= i__2; ++j) {
d__2 = w[j] + w[j - 1];
/* Computing 2nd power */
d__1 = epslon_(&d__2);
if (bd[j] <= d__1 * d__1) {
bd[j] = 0.;
}
/* L200: */
}
L210:
f = bd[*n] / delta;
qp = delta + f;
p = 1.;
if (k == *n) {
goto L260;
}
k1 = *n - k;
/* .......... FOR I=N-1 STEP -1 UNTIL K DO -- .......... */
i__2 = k1;
for (ii = 1; ii <= i__2; ++ii) {
i__ = *n - ii;
q = w[i__] - s - f;
r__ = q / qp;
p = p * r__ + 1.;
ep = f * r__;
w[i__ + 1] = qp + ep;
delta = q - ep;
if (delta > *eps1) {
goto L220;
}
if (delta < -(*eps1)) {
goto L1000;
}
goto L300;
L220:
f = bd[i__] / q;
qp = delta + f;
bd[i__ + 1] = qp * ep;
/* L240: */
}
L260:
w[k] = qp;
s = qp / p;
if (tot + s > tot) {
goto L180;
}
/* .......... SET ERROR -- IRREGULAR END OF ITERATION. */
/* DEFLATE MINIMUM DIAGONAL ELEMENT .......... */
*ierr = *n * 5 + k;
s = 0.;
delta = qp;
i__2 = *n;
for (j = k; j <= i__2; ++j) {
if (w[j] > delta) {
goto L280;
}
i__ = j;
delta = w[j];
L280:
;
}
/* .......... CONVERGENCE .......... */
L300:
if (i__ < *n) {
bd[i__ + 1] = bd[i__] * f / qp;
}
ii = ind[i__];
if (i__ == k) {
goto L340;
}
k1 = i__ - k;
/* .......... FOR J=I-1 STEP -1 UNTIL K DO -- .......... */
i__2 = k1;
for (jj = 1; jj <= i__2; ++jj) {
j = i__ - jj;
w[j + 1] = w[j] - s;
bd[j + 1] = bd[j];
ind[j + 1] = ind[j];
/* L320: */
}
L340:
w[k] = tot;
err += abs(delta);
bd[k] = err;
ind[k] = ii;
/* L360: */
}
if (*type__) {
goto L1001;
}
f = bd[1];
e2[1] = 2.;
bd[1] = f;
j = 2;
/* .......... NEGATE ELEMENTS OF W FOR LARGEST VALUES .......... */
L400:
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
/* L500: */
w[i__] = -w[i__];
}
jdef = -jdef;
switch (j) {
case 1: goto L40;
case 2: goto L1001;
}
/* .......... SET ERROR -- IDEF SPECIFIED INCORRECTLY .......... */
L1000:
*ierr = *n * 6 + 1;
L1001:
return 0;
} /* ratqr_ */
syntax highlighted by Code2HTML, v. 0.9.1