clar1v.c

Go to the documentation of this file.

/* 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"
#include "blaswrap.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 */
    integer i__;
    real s;
    integer r1, r2;
    real eps, tmp;
    integer neg1, neg2, indp, inds;
    real dplus;
    extern doublereal slamch_(char *);
    integer indlpl, indumn;
    extern logical sisnan_(real *);
    real dminus;
    logical sawnan1, sawnan2;


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

/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */

/*  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 */

/*  ===================================================================== */

/*     .. Parameters .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. External Functions .. */
/*     .. */
/*     .. Intrinsic Functions .. */
/*     .. */
/*     .. Executable Statements .. */

    /* 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_ */