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