#include "blaswrap.h"
/* clar1v.f -- translated by f2c (version 20061008).
   You must link the resulting object file with libf2c:
	on Microsoft Windows system, link with libf2c.lib;
	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
	or, if you install libf2c.a in a standard place, with -lf2c -lm
	-- in that order, at the end of the command line, as in
		cc *.o -lf2c -lm
	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,

		http://www.netlib.org/f2c/libf2c.zip
*/

#include "f2c.h"

/* Subroutine */ int clar1v_(integer *n, integer *b1, integer *bn, real *
	lambda, real *d__, real *l, real *ld, real *lld, real *pivmin, real *
	gaptol, complex *z__, logical *wantnc, integer *negcnt, real *ztz, 
	real *mingma, integer *r__, integer *isuppz, real *nrminv, real *
	resid, real *rqcorr, real *work)
{
    /* System generated locals */
    integer i__1, i__2, i__3, i__4;
    real r__1;
    complex q__1, q__2;

    /* Builtin functions */
    double c_abs(complex *), sqrt(doublereal);

    /* Local variables */
    static integer i__;
    static real s;
    static integer r1, r2;
    static real eps, tmp;
    static integer neg1, neg2, indp, inds;
    static real dplus;
    extern doublereal slamch_(char *);
    static integer indlpl, indumn;
    extern logical sisnan_(real *);
    static real dminus;
    static logical sawnan1, sawnan2;


/*  -- LAPACK auxiliary routine (version 3.1) --   
       Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..   
       November 2006   


    Purpose   
    =======   

    CLAR1V computes the (scaled) r-th column of the inverse of   
    the sumbmatrix in rows B1 through BN of the tridiagonal matrix   
    L D L^T - sigma I. When sigma is close to an eigenvalue, the   
    computed vector is an accurate eigenvector. Usually, r corresponds   
    to the index where the eigenvector is largest in magnitude.   
    The following steps accomplish this computation :   
    (a) Stationary qd transform,  L D L^T - sigma I = L(+) D(+) L(+)^T,   
    (b) Progressive qd transform, L D L^T - sigma I = U(-) D(-) U(-)^T,   
    (c) Computation of the diagonal elements of the inverse of   
        L D L^T - sigma I by combining the above transforms, and choosing   
        r as the index where the diagonal of the inverse is (one of the)   
        largest in magnitude.   
    (d) Computation of the (scaled) r-th column of the inverse using the   
        twisted factorization obtained by combining the top part of the   
        the stationary and the bottom part of the progressive transform.   

    Arguments   
    =========   

    N        (input) INTEGER   
             The order of the matrix L D L^T.   

    B1       (input) INTEGER   
             First index of the submatrix of L D L^T.   

    BN       (input) INTEGER   
             Last index of the submatrix of L D L^T.   

    LAMBDA    (input) REAL   
             The shift. In order to compute an accurate eigenvector,   
             LAMBDA should be a good approximation to an eigenvalue   
             of L D L^T.   

    L        (input) REAL             array, dimension (N-1)   
             The (n-1) subdiagonal elements of the unit bidiagonal matrix   
             L, in elements 1 to N-1.   

    D        (input) REAL             array, dimension (N)   
             The n diagonal elements of the diagonal matrix D.   

    LD       (input) REAL             array, dimension (N-1)   
             The n-1 elements L(i)*D(i).   

    LLD      (input) REAL             array, dimension (N-1)   
             The n-1 elements L(i)*L(i)*D(i).   

    PIVMIN   (input) REAL   
             The minimum pivot in the Sturm sequence.   

    GAPTOL   (input) REAL   
             Tolerance that indicates when eigenvector entries are negligible   
             w.r.t. their contribution to the residual.   

    Z        (input/output) COMPLEX          array, dimension (N)   
             On input, all entries of Z must be set to 0.   
             On output, Z contains the (scaled) r-th column of the   
             inverse. The scaling is such that Z(R) equals 1.   

    WANTNC   (input) LOGICAL   
             Specifies whether NEGCNT has to be computed.   

    NEGCNT   (output) INTEGER   
             If WANTNC is .TRUE. then NEGCNT = the number of pivots < pivmin   
             in the  matrix factorization L D L^T, and NEGCNT = -1 otherwise.   

    ZTZ      (output) REAL   
             The square of the 2-norm of Z.   

    MINGMA   (output) REAL   
             The reciprocal of the largest (in magnitude) diagonal   
             element of the inverse of L D L^T - sigma I.   

    R        (input/output) INTEGER   
             The twist index for the twisted factorization used to   
             compute Z.   
             On input, 0 <= R <= N. If R is input as 0, R is set to   
             the index where (L D L^T - sigma I)^{-1} is largest   
             in magnitude. If 1 <= R <= N, R is unchanged.   
             On output, R contains the twist index used to compute Z.   
             Ideally, R designates the position of the maximum entry in the   
             eigenvector.   

    ISUPPZ   (output) INTEGER array, dimension (2)   
             The support of the vector in Z, i.e., the vector Z is   
             nonzero only in elements ISUPPZ(1) through ISUPPZ( 2 ).   

    NRMINV   (output) REAL   
             NRMINV = 1/SQRT( ZTZ )   

    RESID    (output) REAL   
             The residual of the FP vector.   
             RESID = ABS( MINGMA )/SQRT( ZTZ )   

    RQCORR   (output) REAL   
             The Rayleigh Quotient correction to LAMBDA.   
             RQCORR = MINGMA*TMP   

    WORK     (workspace) REAL             array, dimension (4*N)   

    Further Details   
    ===============   

    Based on contributions by   
       Beresford Parlett, University of California, Berkeley, USA   
       Jim Demmel, University of California, Berkeley, USA   
       Inderjit Dhillon, University of Texas, Austin, USA   
       Osni Marques, LBNL/NERSC, USA   
       Christof Voemel, University of California, Berkeley, USA   

    =====================================================================   


       Parameter adjustments */
    --work;
    --isuppz;
    --z__;
    --lld;
    --ld;
    --l;
    --d__;

    /* Function Body */
    eps = slamch_("Precision");
    if (*r__ == 0) {
	r1 = *b1;
	r2 = *bn;
    } else {
	r1 = *r__;
	r2 = *r__;
    }
/*     Storage for LPLUS */
    indlpl = 0;
/*     Storage for UMINUS */
    indumn = *n;
    inds = (*n << 1) + 1;
    indp = *n * 3 + 1;
    if (*b1 == 1) {
	work[inds] = 0.f;
    } else {
	work[inds + *b1 - 1] = lld[*b1 - 1];
    }

/*     Compute the stationary transform (using the differential form)   
       until the index R2. */

    sawnan1 = FALSE_;
    neg1 = 0;
    s = work[inds + *b1 - 1] - *lambda;
    i__1 = r1 - 1;
    for (i__ = *b1; i__ <= i__1; ++i__) {
	dplus = d__[i__] + s;
	work[indlpl + i__] = ld[i__] / dplus;
	if (dplus < 0.f) {
	    ++neg1;
	}
	work[inds + i__] = s * work[indlpl + i__] * l[i__];
	s = work[inds + i__] - *lambda;
/* L50: */
    }
    sawnan1 = sisnan_(&s);
    if (sawnan1) {
	goto L60;
    }
    i__1 = r2 - 1;
    for (i__ = r1; i__ <= i__1; ++i__) {
	dplus = d__[i__] + s;
	work[indlpl + i__] = ld[i__] / dplus;
	work[inds + i__] = s * work[indlpl + i__] * l[i__];
	s = work[inds + i__] - *lambda;
/* L51: */
    }
    sawnan1 = sisnan_(&s);

L60:
    if (sawnan1) {
/*        Runs a slower version of the above loop if a NaN is detected */
	neg1 = 0;
	s = work[inds + *b1 - 1] - *lambda;
	i__1 = r1 - 1;
	for (i__ = *b1; i__ <= i__1; ++i__) {
	    dplus = d__[i__] + s;
	    if (dabs(dplus) < *pivmin) {
		dplus = -(*pivmin);
	    }
	    work[indlpl + i__] = ld[i__] / dplus;
	    if (dplus < 0.f) {
		++neg1;
	    }
	    work[inds + i__] = s * work[indlpl + i__] * l[i__];
	    if (work[indlpl + i__] == 0.f) {
		work[inds + i__] = lld[i__];
	    }
	    s = work[inds + i__] - *lambda;
/* L70: */
	}
	i__1 = r2 - 1;
	for (i__ = r1; i__ <= i__1; ++i__) {
	    dplus = d__[i__] + s;
	    if (dabs(dplus) < *pivmin) {
		dplus = -(*pivmin);
	    }
	    work[indlpl + i__] = ld[i__] / dplus;
	    work[inds + i__] = s * work[indlpl + i__] * l[i__];
	    if (work[indlpl + i__] == 0.f) {
		work[inds + i__] = lld[i__];
	    }
	    s = work[inds + i__] - *lambda;
/* L71: */
	}
    }

/*     Compute the progressive transform (using the differential form)   
       until the index R1 */

    sawnan2 = FALSE_;
    neg2 = 0;
    work[indp + *bn - 1] = d__[*bn] - *lambda;
    i__1 = r1;
    for (i__ = *bn - 1; i__ >= i__1; --i__) {
	dminus = lld[i__] + work[indp + i__];
	tmp = d__[i__] / dminus;
	if (dminus < 0.f) {
	    ++neg2;
	}
	work[indumn + i__] = l[i__] * tmp;
	work[indp + i__ - 1] = work[indp + i__] * tmp - *lambda;
/* L80: */
    }
    tmp = work[indp + r1 - 1];
    sawnan2 = sisnan_(&tmp);
    if (sawnan2) {
/*        Runs a slower version of the above loop if a NaN is detected */
	neg2 = 0;
	i__1 = r1;
	for (i__ = *bn - 1; i__ >= i__1; --i__) {
	    dminus = lld[i__] + work[indp + i__];
	    if (dabs(dminus) < *pivmin) {
		dminus = -(*pivmin);
	    }
	    tmp = d__[i__] / dminus;
	    if (dminus < 0.f) {
		++neg2;
	    }
	    work[indumn + i__] = l[i__] * tmp;
	    work[indp + i__ - 1] = work[indp + i__] * tmp - *lambda;
	    if (tmp == 0.f) {
		work[indp + i__ - 1] = d__[i__] - *lambda;
	    }
/* L100: */
	}
    }

/*     Find the index (from R1 to R2) of the largest (in magnitude)   
       diagonal element of the inverse */

    *mingma = work[inds + r1 - 1] + work[indp + r1 - 1];
    if (*mingma < 0.f) {
	++neg1;
    }
    if (*wantnc) {
	*negcnt = neg1 + neg2;
    } else {
	*negcnt = -1;
    }
    if (dabs(*mingma) == 0.f) {
	*mingma = eps * work[inds + r1 - 1];
    }
    *r__ = r1;
    i__1 = r2 - 1;
    for (i__ = r1; i__ <= i__1; ++i__) {
	tmp = work[inds + i__] + work[indp + i__];
	if (tmp == 0.f) {
	    tmp = eps * work[inds + i__];
	}
	if (dabs(tmp) <= dabs(*mingma)) {
	    *mingma = tmp;
	    *r__ = i__ + 1;
	}
/* L110: */
    }

/*     Compute the FP vector: solve N^T v = e_r */

    isuppz[1] = *b1;
    isuppz[2] = *bn;
    i__1 = *r__;
    z__[i__1].r = 1.f, z__[i__1].i = 0.f;
    *ztz = 1.f;

/*     Compute the FP vector upwards from R */

    if (! sawnan1 && ! sawnan2) {
	i__1 = *b1;
	for (i__ = *r__ - 1; i__ >= i__1; --i__) {
	    i__2 = i__;
	    i__3 = indlpl + i__;
	    i__4 = i__ + 1;
	    q__2.r = work[i__3] * z__[i__4].r, q__2.i = work[i__3] * z__[i__4]
		    .i;
	    q__1.r = -q__2.r, q__1.i = -q__2.i;
	    z__[i__2].r = q__1.r, z__[i__2].i = q__1.i;
	    if ((c_abs(&z__[i__]) + c_abs(&z__[i__ + 1])) * (r__1 = ld[i__], 
		    dabs(r__1)) < *gaptol) {
		i__2 = i__;
		z__[i__2].r = 0.f, z__[i__2].i = 0.f;
		isuppz[1] = i__ + 1;
		goto L220;
	    }
	    i__2 = i__;
	    i__3 = i__;
	    q__1.r = z__[i__2].r * z__[i__3].r - z__[i__2].i * z__[i__3].i, 
		    q__1.i = z__[i__2].r * z__[i__3].i + z__[i__2].i * z__[
		    i__3].r;
	    *ztz += q__1.r;
/* L210: */
	}
L220:
	;
    } else {
/*        Run slower loop if NaN occurred. */
	i__1 = *b1;
	for (i__ = *r__ - 1; i__ >= i__1; --i__) {
	    i__2 = i__ + 1;
	    if (z__[i__2].r == 0.f && z__[i__2].i == 0.f) {
		i__2 = i__;
		r__1 = -(ld[i__ + 1] / ld[i__]);
		i__3 = i__ + 2;
		q__1.r = r__1 * z__[i__3].r, q__1.i = r__1 * z__[i__3].i;
		z__[i__2].r = q__1.r, z__[i__2].i = q__1.i;
	    } else {
		i__2 = i__;
		i__3 = indlpl + i__;
		i__4 = i__ + 1;
		q__2.r = work[i__3] * z__[i__4].r, q__2.i = work[i__3] * z__[
			i__4].i;
		q__1.r = -q__2.r, q__1.i = -q__2.i;
		z__[i__2].r = q__1.r, z__[i__2].i = q__1.i;
	    }
	    if ((c_abs(&z__[i__]) + c_abs(&z__[i__ + 1])) * (r__1 = ld[i__], 
		    dabs(r__1)) < *gaptol) {
		i__2 = i__;
		z__[i__2].r = 0.f, z__[i__2].i = 0.f;
		isuppz[1] = i__ + 1;
		goto L240;
	    }
	    i__2 = i__;
	    i__3 = i__;
	    q__1.r = z__[i__2].r * z__[i__3].r - z__[i__2].i * z__[i__3].i, 
		    q__1.i = z__[i__2].r * z__[i__3].i + z__[i__2].i * z__[
		    i__3].r;
	    *ztz += q__1.r;
/* L230: */
	}
L240:
	;
    }
/*     Compute the FP vector downwards from R in blocks of size BLKSIZ */
    if (! sawnan1 && ! sawnan2) {
	i__1 = *bn - 1;
	for (i__ = *r__; i__ <= i__1; ++i__) {
	    i__2 = i__ + 1;
	    i__3 = indumn + i__;
	    i__4 = i__;
	    q__2.r = work[i__3] * z__[i__4].r, q__2.i = work[i__3] * z__[i__4]
		    .i;
	    q__1.r = -q__2.r, q__1.i = -q__2.i;
	    z__[i__2].r = q__1.r, z__[i__2].i = q__1.i;
	    if ((c_abs(&z__[i__]) + c_abs(&z__[i__ + 1])) * (r__1 = ld[i__], 
		    dabs(r__1)) < *gaptol) {
		i__2 = i__ + 1;
		z__[i__2].r = 0.f, z__[i__2].i = 0.f;
		isuppz[2] = i__;
		goto L260;
	    }
	    i__2 = i__ + 1;
	    i__3 = i__ + 1;
	    q__1.r = z__[i__2].r * z__[i__3].r - z__[i__2].i * z__[i__3].i, 
		    q__1.i = z__[i__2].r * z__[i__3].i + z__[i__2].i * z__[
		    i__3].r;
	    *ztz += q__1.r;
/* L250: */
	}
L260:
	;
    } else {
/*        Run slower loop if NaN occurred. */
	i__1 = *bn - 1;
	for (i__ = *r__; i__ <= i__1; ++i__) {
	    i__2 = i__;
	    if (z__[i__2].r == 0.f && z__[i__2].i == 0.f) {
		i__2 = i__ + 1;
		r__1 = -(ld[i__ - 1] / ld[i__]);
		i__3 = i__ - 1;
		q__1.r = r__1 * z__[i__3].r, q__1.i = r__1 * z__[i__3].i;
		z__[i__2].r = q__1.r, z__[i__2].i = q__1.i;
	    } else {
		i__2 = i__ + 1;
		i__3 = indumn + i__;
		i__4 = i__;
		q__2.r = work[i__3] * z__[i__4].r, q__2.i = work[i__3] * z__[
			i__4].i;
		q__1.r = -q__2.r, q__1.i = -q__2.i;
		z__[i__2].r = q__1.r, z__[i__2].i = q__1.i;
	    }
	    if ((c_abs(&z__[i__]) + c_abs(&z__[i__ + 1])) * (r__1 = ld[i__], 
		    dabs(r__1)) < *gaptol) {
		i__2 = i__ + 1;
		z__[i__2].r = 0.f, z__[i__2].i = 0.f;
		isuppz[2] = i__;
		goto L280;
	    }
	    i__2 = i__ + 1;
	    i__3 = i__ + 1;
	    q__1.r = z__[i__2].r * z__[i__3].r - z__[i__2].i * z__[i__3].i, 
		    q__1.i = z__[i__2].r * z__[i__3].i + z__[i__2].i * z__[
		    i__3].r;
	    *ztz += q__1.r;
/* L270: */
	}
L280:
	;
    }

/*     Compute quantities for convergence test */

    tmp = 1.f / *ztz;
    *nrminv = sqrt(tmp);
    *resid = dabs(*mingma) * *nrminv;
    *rqcorr = *mingma * tmp;


    return 0;

/*     End of CLAR1V */

} /* clar1v_ */