001: SUBROUTINE DLAR1V( N, B1, BN, LAMBDA, D, L, LD, LLD,
002: $ PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA,
003: $ R, ISUPPZ, NRMINV, RESID, RQCORR, WORK )
004: *
005: * -- LAPACK auxiliary routine (version 3.2) --
006: * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
007: * November 2006
008: *
009: * .. Scalar Arguments ..
010: LOGICAL WANTNC
011: INTEGER B1, BN, N, NEGCNT, R
012: DOUBLE PRECISION GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID,
013: $ RQCORR, ZTZ
014: * ..
015: * .. Array Arguments ..
016: INTEGER ISUPPZ( * )
017: DOUBLE PRECISION D( * ), L( * ), LD( * ), LLD( * ),
018: $ WORK( * )
019: DOUBLE PRECISION Z( * )
020: * ..
021: *
022: * Purpose
023: * =======
024: *
025: * DLAR1V computes the (scaled) r-th column of the inverse of
026: * the sumbmatrix in rows B1 through BN of the tridiagonal matrix
027: * L D L^T - sigma I. When sigma is close to an eigenvalue, the
028: * computed vector is an accurate eigenvector. Usually, r corresponds
029: * to the index where the eigenvector is largest in magnitude.
030: * The following steps accomplish this computation :
031: * (a) Stationary qd transform, L D L^T - sigma I = L(+) D(+) L(+)^T,
032: * (b) Progressive qd transform, L D L^T - sigma I = U(-) D(-) U(-)^T,
033: * (c) Computation of the diagonal elements of the inverse of
034: * L D L^T - sigma I by combining the above transforms, and choosing
035: * r as the index where the diagonal of the inverse is (one of the)
036: * largest in magnitude.
037: * (d) Computation of the (scaled) r-th column of the inverse using the
038: * twisted factorization obtained by combining the top part of the
039: * the stationary and the bottom part of the progressive transform.
040: *
041: * Arguments
042: * =========
043: *
044: * N (input) INTEGER
045: * The order of the matrix L D L^T.
046: *
047: * B1 (input) INTEGER
048: * First index of the submatrix of L D L^T.
049: *
050: * BN (input) INTEGER
051: * Last index of the submatrix of L D L^T.
052: *
053: * LAMBDA (input) DOUBLE PRECISION
054: * The shift. In order to compute an accurate eigenvector,
055: * LAMBDA should be a good approximation to an eigenvalue
056: * of L D L^T.
057: *
058: * L (input) DOUBLE PRECISION array, dimension (N-1)
059: * The (n-1) subdiagonal elements of the unit bidiagonal matrix
060: * L, in elements 1 to N-1.
061: *
062: * D (input) DOUBLE PRECISION array, dimension (N)
063: * The n diagonal elements of the diagonal matrix D.
064: *
065: * LD (input) DOUBLE PRECISION array, dimension (N-1)
066: * The n-1 elements L(i)*D(i).
067: *
068: * LLD (input) DOUBLE PRECISION array, dimension (N-1)
069: * The n-1 elements L(i)*L(i)*D(i).
070: *
071: * PIVMIN (input) DOUBLE PRECISION
072: * The minimum pivot in the Sturm sequence.
073: *
074: * GAPTOL (input) DOUBLE PRECISION
075: * Tolerance that indicates when eigenvector entries are negligible
076: * w.r.t. their contribution to the residual.
077: *
078: * Z (input/output) DOUBLE PRECISION array, dimension (N)
079: * On input, all entries of Z must be set to 0.
080: * On output, Z contains the (scaled) r-th column of the
081: * inverse. The scaling is such that Z(R) equals 1.
082: *
083: * WANTNC (input) LOGICAL
084: * Specifies whether NEGCNT has to be computed.
085: *
086: * NEGCNT (output) INTEGER
087: * If WANTNC is .TRUE. then NEGCNT = the number of pivots < pivmin
088: * in the matrix factorization L D L^T, and NEGCNT = -1 otherwise.
089: *
090: * ZTZ (output) DOUBLE PRECISION
091: * The square of the 2-norm of Z.
092: *
093: * MINGMA (output) DOUBLE PRECISION
094: * The reciprocal of the largest (in magnitude) diagonal
095: * element of the inverse of L D L^T - sigma I.
096: *
097: * R (input/output) INTEGER
098: * The twist index for the twisted factorization used to
099: * compute Z.
100: * On input, 0 <= R <= N. If R is input as 0, R is set to
101: * the index where (L D L^T - sigma I)^{-1} is largest
102: * in magnitude. If 1 <= R <= N, R is unchanged.
103: * On output, R contains the twist index used to compute Z.
104: * Ideally, R designates the position of the maximum entry in the
105: * eigenvector.
106: *
107: * ISUPPZ (output) INTEGER array, dimension (2)
108: * The support of the vector in Z, i.e., the vector Z is
109: * nonzero only in elements ISUPPZ(1) through ISUPPZ( 2 ).
110: *
111: * NRMINV (output) DOUBLE PRECISION
112: * NRMINV = 1/SQRT( ZTZ )
113: *
114: * RESID (output) DOUBLE PRECISION
115: * The residual of the FP vector.
116: * RESID = ABS( MINGMA )/SQRT( ZTZ )
117: *
118: * RQCORR (output) DOUBLE PRECISION
119: * The Rayleigh Quotient correction to LAMBDA.
120: * RQCORR = MINGMA*TMP
121: *
122: * WORK (workspace) DOUBLE PRECISION array, dimension (4*N)
123: *
124: * Further Details
125: * ===============
126: *
127: * Based on contributions by
128: * Beresford Parlett, University of California, Berkeley, USA
129: * Jim Demmel, University of California, Berkeley, USA
130: * Inderjit Dhillon, University of Texas, Austin, USA
131: * Osni Marques, LBNL/NERSC, USA
132: * Christof Voemel, University of California, Berkeley, USA
133: *
134: * =====================================================================
135: *
136: * .. Parameters ..
137: DOUBLE PRECISION ZERO, ONE
138: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
139:
140: * ..
141: * .. Local Scalars ..
142: LOGICAL SAWNAN1, SAWNAN2
143: INTEGER I, INDLPL, INDP, INDS, INDUMN, NEG1, NEG2, R1,
144: $ R2
145: DOUBLE PRECISION DMINUS, DPLUS, EPS, S, TMP
146: * ..
147: * .. External Functions ..
148: LOGICAL DISNAN
149: DOUBLE PRECISION DLAMCH
150: EXTERNAL DISNAN, DLAMCH
151: * ..
152: * .. Intrinsic Functions ..
153: INTRINSIC ABS
154: * ..
155: * .. Executable Statements ..
156: *
157: EPS = DLAMCH( 'Precision' )
158:
159:
160: IF( R.EQ.0 ) THEN
161: R1 = B1
162: R2 = BN
163: ELSE
164: R1 = R
165: R2 = R
166: END IF
167:
168: * Storage for LPLUS
169: INDLPL = 0
170: * Storage for UMINUS
171: INDUMN = N
172: INDS = 2*N + 1
173: INDP = 3*N + 1
174:
175: IF( B1.EQ.1 ) THEN
176: WORK( INDS ) = ZERO
177: ELSE
178: WORK( INDS+B1-1 ) = LLD( B1-1 )
179: END IF
180:
181: *
182: * Compute the stationary transform (using the differential form)
183: * until the index R2.
184: *
185: SAWNAN1 = .FALSE.
186: NEG1 = 0
187: S = WORK( INDS+B1-1 ) - LAMBDA
188: DO 50 I = B1, R1 - 1
189: DPLUS = D( I ) + S
190: WORK( INDLPL+I ) = LD( I ) / DPLUS
191: IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
192: WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
193: S = WORK( INDS+I ) - LAMBDA
194: 50 CONTINUE
195: SAWNAN1 = DISNAN( S )
196: IF( SAWNAN1 ) GOTO 60
197: DO 51 I = R1, R2 - 1
198: DPLUS = D( I ) + S
199: WORK( INDLPL+I ) = LD( I ) / DPLUS
200: WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
201: S = WORK( INDS+I ) - LAMBDA
202: 51 CONTINUE
203: SAWNAN1 = DISNAN( S )
204: *
205: 60 CONTINUE
206: IF( SAWNAN1 ) THEN
207: * Runs a slower version of the above loop if a NaN is detected
208: NEG1 = 0
209: S = WORK( INDS+B1-1 ) - LAMBDA
210: DO 70 I = B1, R1 - 1
211: DPLUS = D( I ) + S
212: IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
213: WORK( INDLPL+I ) = LD( I ) / DPLUS
214: IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
215: WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
216: IF( WORK( INDLPL+I ).EQ.ZERO )
217: $ WORK( INDS+I ) = LLD( I )
218: S = WORK( INDS+I ) - LAMBDA
219: 70 CONTINUE
220: DO 71 I = R1, R2 - 1
221: DPLUS = D( I ) + S
222: IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
223: WORK( INDLPL+I ) = LD( I ) / DPLUS
224: WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
225: IF( WORK( INDLPL+I ).EQ.ZERO )
226: $ WORK( INDS+I ) = LLD( I )
227: S = WORK( INDS+I ) - LAMBDA
228: 71 CONTINUE
229: END IF
230: *
231: * Compute the progressive transform (using the differential form)
232: * until the index R1
233: *
234: SAWNAN2 = .FALSE.
235: NEG2 = 0
236: WORK( INDP+BN-1 ) = D( BN ) - LAMBDA
237: DO 80 I = BN - 1, R1, -1
238: DMINUS = LLD( I ) + WORK( INDP+I )
239: TMP = D( I ) / DMINUS
240: IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
241: WORK( INDUMN+I ) = L( I )*TMP
242: WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
243: 80 CONTINUE
244: TMP = WORK( INDP+R1-1 )
245: SAWNAN2 = DISNAN( TMP )
246:
247: IF( SAWNAN2 ) THEN
248: * Runs a slower version of the above loop if a NaN is detected
249: NEG2 = 0
250: DO 100 I = BN-1, R1, -1
251: DMINUS = LLD( I ) + WORK( INDP+I )
252: IF(ABS(DMINUS).LT.PIVMIN) DMINUS = -PIVMIN
253: TMP = D( I ) / DMINUS
254: IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
255: WORK( INDUMN+I ) = L( I )*TMP
256: WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
257: IF( TMP.EQ.ZERO )
258: $ WORK( INDP+I-1 ) = D( I ) - LAMBDA
259: 100 CONTINUE
260: END IF
261: *
262: * Find the index (from R1 to R2) of the largest (in magnitude)
263: * diagonal element of the inverse
264: *
265: MINGMA = WORK( INDS+R1-1 ) + WORK( INDP+R1-1 )
266: IF( MINGMA.LT.ZERO ) NEG1 = NEG1 + 1
267: IF( WANTNC ) THEN
268: NEGCNT = NEG1 + NEG2
269: ELSE
270: NEGCNT = -1
271: ENDIF
272: IF( ABS(MINGMA).EQ.ZERO )
273: $ MINGMA = EPS*WORK( INDS+R1-1 )
274: R = R1
275: DO 110 I = R1, R2 - 1
276: TMP = WORK( INDS+I ) + WORK( INDP+I )
277: IF( TMP.EQ.ZERO )
278: $ TMP = EPS*WORK( INDS+I )
279: IF( ABS( TMP ).LE.ABS( MINGMA ) ) THEN
280: MINGMA = TMP
281: R = I + 1
282: END IF
283: 110 CONTINUE
284: *
285: * Compute the FP vector: solve N^T v = e_r
286: *
287: ISUPPZ( 1 ) = B1
288: ISUPPZ( 2 ) = BN
289: Z( R ) = ONE
290: ZTZ = ONE
291: *
292: * Compute the FP vector upwards from R
293: *
294: IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
295: DO 210 I = R-1, B1, -1
296: Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
297: IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
298: $ THEN
299: Z( I ) = ZERO
300: ISUPPZ( 1 ) = I + 1
301: GOTO 220
302: ENDIF
303: ZTZ = ZTZ + Z( I )*Z( I )
304: 210 CONTINUE
305: 220 CONTINUE
306: ELSE
307: * Run slower loop if NaN occurred.
308: DO 230 I = R - 1, B1, -1
309: IF( Z( I+1 ).EQ.ZERO ) THEN
310: Z( I ) = -( LD( I+1 ) / LD( I ) )*Z( I+2 )
311: ELSE
312: Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
313: END IF
314: IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
315: $ THEN
316: Z( I ) = ZERO
317: ISUPPZ( 1 ) = I + 1
318: GO TO 240
319: END IF
320: ZTZ = ZTZ + Z( I )*Z( I )
321: 230 CONTINUE
322: 240 CONTINUE
323: ENDIF
324:
325: * Compute the FP vector downwards from R in blocks of size BLKSIZ
326: IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
327: DO 250 I = R, BN-1
328: Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
329: IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
330: $ THEN
331: Z( I+1 ) = ZERO
332: ISUPPZ( 2 ) = I
333: GO TO 260
334: END IF
335: ZTZ = ZTZ + Z( I+1 )*Z( I+1 )
336: 250 CONTINUE
337: 260 CONTINUE
338: ELSE
339: * Run slower loop if NaN occurred.
340: DO 270 I = R, BN - 1
341: IF( Z( I ).EQ.ZERO ) THEN
342: Z( I+1 ) = -( LD( I-1 ) / LD( I ) )*Z( I-1 )
343: ELSE
344: Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
345: END IF
346: IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
347: $ THEN
348: Z( I+1 ) = ZERO
349: ISUPPZ( 2 ) = I
350: GO TO 280
351: END IF
352: ZTZ = ZTZ + Z( I+1 )*Z( I+1 )
353: 270 CONTINUE
354: 280 CONTINUE
355: END IF
356: *
357: * Compute quantities for convergence test
358: *
359: TMP = ONE / ZTZ
360: NRMINV = SQRT( TMP )
361: RESID = ABS( MINGMA )*NRMINV
362: RQCORR = MINGMA*TMP
363: *
364: *
365: RETURN
366: *
367: * End of DLAR1V
368: *
369: END
370: