SUBROUTINE CLAR1V( N, B1, BN, LAMBDA, D, L, LD, LLD, \$ PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA, \$ R, ISUPPZ, NRMINV, RESID, RQCORR, WORK ) * * -- LAPACK auxiliary routine (version 3.2) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * November 2006 * * .. Scalar Arguments .. LOGICAL WANTNC INTEGER B1, BN, N, NEGCNT, R REAL GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID, \$ RQCORR, ZTZ * .. * .. Array Arguments .. INTEGER ISUPPZ( * ) REAL D( * ), L( * ), LD( * ), LLD( * ), \$ WORK( * ) COMPLEX Z( * ) * .. * * 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 .. REAL ZERO, ONE PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) COMPLEX CONE PARAMETER ( CONE = ( 1.0E0, 0.0E0 ) ) * .. * .. Local Scalars .. LOGICAL SAWNAN1, SAWNAN2 INTEGER I, INDLPL, INDP, INDS, INDUMN, NEG1, NEG2, R1, \$ R2 REAL DMINUS, DPLUS, EPS, S, TMP * .. * .. External Functions .. LOGICAL SISNAN REAL SLAMCH EXTERNAL SISNAN, SLAMCH * .. * .. Intrinsic Functions .. INTRINSIC ABS, REAL * .. * .. Executable Statements .. * EPS = SLAMCH( 'Precision' ) IF( R.EQ.0 ) THEN R1 = B1 R2 = BN ELSE R1 = R R2 = R END IF * Storage for LPLUS INDLPL = 0 * Storage for UMINUS INDUMN = N INDS = 2*N + 1 INDP = 3*N + 1 IF( B1.EQ.1 ) THEN WORK( INDS ) = ZERO ELSE WORK( INDS+B1-1 ) = LLD( B1-1 ) END IF * * Compute the stationary transform (using the differential form) * until the index R2. * SAWNAN1 = .FALSE. NEG1 = 0 S = WORK( INDS+B1-1 ) - LAMBDA DO 50 I = B1, R1 - 1 DPLUS = D( I ) + S WORK( INDLPL+I ) = LD( I ) / DPLUS IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1 WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I ) S = WORK( INDS+I ) - LAMBDA 50 CONTINUE SAWNAN1 = SISNAN( S ) IF( SAWNAN1 ) GOTO 60 DO 51 I = R1, R2 - 1 DPLUS = D( I ) + S WORK( INDLPL+I ) = LD( I ) / DPLUS WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I ) S = WORK( INDS+I ) - LAMBDA 51 CONTINUE SAWNAN1 = SISNAN( S ) * 60 CONTINUE IF( SAWNAN1 ) THEN * Runs a slower version of the above loop if a NaN is detected NEG1 = 0 S = WORK( INDS+B1-1 ) - LAMBDA DO 70 I = B1, R1 - 1 DPLUS = D( I ) + S IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN WORK( INDLPL+I ) = LD( I ) / DPLUS IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1 WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I ) IF( WORK( INDLPL+I ).EQ.ZERO ) \$ WORK( INDS+I ) = LLD( I ) S = WORK( INDS+I ) - LAMBDA 70 CONTINUE DO 71 I = R1, R2 - 1 DPLUS = D( I ) + S IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN WORK( INDLPL+I ) = LD( I ) / DPLUS WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I ) IF( WORK( INDLPL+I ).EQ.ZERO ) \$ WORK( INDS+I ) = LLD( I ) S = WORK( INDS+I ) - LAMBDA 71 CONTINUE END IF * * Compute the progressive transform (using the differential form) * until the index R1 * SAWNAN2 = .FALSE. NEG2 = 0 WORK( INDP+BN-1 ) = D( BN ) - LAMBDA DO 80 I = BN - 1, R1, -1 DMINUS = LLD( I ) + WORK( INDP+I ) TMP = D( I ) / DMINUS IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1 WORK( INDUMN+I ) = L( I )*TMP WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA 80 CONTINUE TMP = WORK( INDP+R1-1 ) SAWNAN2 = SISNAN( TMP ) IF( SAWNAN2 ) THEN * Runs a slower version of the above loop if a NaN is detected NEG2 = 0 DO 100 I = BN-1, R1, -1 DMINUS = LLD( I ) + WORK( INDP+I ) IF(ABS(DMINUS).LT.PIVMIN) DMINUS = -PIVMIN TMP = D( I ) / DMINUS IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1 WORK( INDUMN+I ) = L( I )*TMP WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA IF( TMP.EQ.ZERO ) \$ WORK( INDP+I-1 ) = D( I ) - LAMBDA 100 CONTINUE END IF * * 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.LT.ZERO ) NEG1 = NEG1 + 1 IF( WANTNC ) THEN NEGCNT = NEG1 + NEG2 ELSE NEGCNT = -1 ENDIF IF( ABS(MINGMA).EQ.ZERO ) \$ MINGMA = EPS*WORK( INDS+R1-1 ) R = R1 DO 110 I = R1, R2 - 1 TMP = WORK( INDS+I ) + WORK( INDP+I ) IF( TMP.EQ.ZERO ) \$ TMP = EPS*WORK( INDS+I ) IF( ABS( TMP ).LE.ABS( MINGMA ) ) THEN MINGMA = TMP R = I + 1 END IF 110 CONTINUE * * Compute the FP vector: solve N^T v = e_r * ISUPPZ( 1 ) = B1 ISUPPZ( 2 ) = BN Z( R ) = CONE ZTZ = ONE * * Compute the FP vector upwards from R * IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN DO 210 I = R-1, B1, -1 Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) ) IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL ) \$ THEN Z( I ) = ZERO ISUPPZ( 1 ) = I + 1 GOTO 220 ENDIF ZTZ = ZTZ + REAL( Z( I )*Z( I ) ) 210 CONTINUE 220 CONTINUE ELSE * Run slower loop if NaN occurred. DO 230 I = R - 1, B1, -1 IF( Z( I+1 ).EQ.ZERO ) THEN Z( I ) = -( LD( I+1 ) / LD( I ) )*Z( I+2 ) ELSE Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) ) END IF IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL ) \$ THEN Z( I ) = ZERO ISUPPZ( 1 ) = I + 1 GO TO 240 END IF ZTZ = ZTZ + REAL( Z( I )*Z( I ) ) 230 CONTINUE 240 CONTINUE ENDIF * Compute the FP vector downwards from R in blocks of size BLKSIZ IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN DO 250 I = R, BN-1 Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) ) IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL ) \$ THEN Z( I+1 ) = ZERO ISUPPZ( 2 ) = I GO TO 260 END IF ZTZ = ZTZ + REAL( Z( I+1 )*Z( I+1 ) ) 250 CONTINUE 260 CONTINUE ELSE * Run slower loop if NaN occurred. DO 270 I = R, BN - 1 IF( Z( I ).EQ.ZERO ) THEN Z( I+1 ) = -( LD( I-1 ) / LD( I ) )*Z( I-1 ) ELSE Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) ) END IF IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL ) \$ THEN Z( I+1 ) = ZERO ISUPPZ( 2 ) = I GO TO 280 END IF ZTZ = ZTZ + REAL( Z( I+1 )*Z( I+1 ) ) 270 CONTINUE 280 CONTINUE END IF * * Compute quantities for convergence test * TMP = ONE / ZTZ NRMINV = SQRT( TMP ) RESID = ABS( MINGMA )*NRMINV RQCORR = MINGMA*TMP * * RETURN * * End of CLAR1V * END