/* ratqr.f -- translated by f2c (version 19961017).
   You must link the resulting object file with the libraries:
	-lf2c -lm   (in that order)
*/

#include <math.h>
#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 */
    integer jdef;
    doublereal f;
    integer i__, j, k;
    doublereal p, q, r__, s, delta;
    integer k1, ii, jj;
    doublereal ep, qp=0.0;
    extern doublereal epslon_(doublereal *);
    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_ */

