1: SUBROUTINE ZLAR1V( N, B1, BN, LAMBDA, D, L, LD, LLD,
2: $ PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA,
3: $ R, ISUPPZ, NRMINV, RESID, RQCORR, WORK )
4: *
5: * -- LAPACK auxiliary routine (version 3.2) --
6: * -- LAPACK is a software package provided by Univ. of Tennessee, --
7: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8: * November 2006
9: *
10: * .. Scalar Arguments ..
11: LOGICAL WANTNC
12: INTEGER B1, BN, N, NEGCNT, R
13: DOUBLE PRECISION GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID,
14: $ RQCORR, ZTZ
15: * ..
16: * .. Array Arguments ..
17: INTEGER ISUPPZ( * )
18: DOUBLE PRECISION D( * ), L( * ), LD( * ), LLD( * ),
19: $ WORK( * )
20: COMPLEX*16 Z( * )
21: * ..
22: *
23: * Purpose
24: * =======
25: *
26: * ZLAR1V computes the (scaled) r-th column of the inverse of
27: * the sumbmatrix in rows B1 through BN of the tridiagonal matrix
28: * L D L^T - sigma I. When sigma is close to an eigenvalue, the
29: * computed vector is an accurate eigenvector. Usually, r corresponds
30: * to the index where the eigenvector is largest in magnitude.
31: * The following steps accomplish this computation :
32: * (a) Stationary qd transform, L D L^T - sigma I = L(+) D(+) L(+)^T,
33: * (b) Progressive qd transform, L D L^T - sigma I = U(-) D(-) U(-)^T,
34: * (c) Computation of the diagonal elements of the inverse of
35: * L D L^T - sigma I by combining the above transforms, and choosing
36: * r as the index where the diagonal of the inverse is (one of the)
37: * largest in magnitude.
38: * (d) Computation of the (scaled) r-th column of the inverse using the
39: * twisted factorization obtained by combining the top part of the
40: * the stationary and the bottom part of the progressive transform.
41: *
42: * Arguments
43: * =========
44: *
45: * N (input) INTEGER
46: * The order of the matrix L D L^T.
47: *
48: * B1 (input) INTEGER
49: * First index of the submatrix of L D L^T.
50: *
51: * BN (input) INTEGER
52: * Last index of the submatrix of L D L^T.
53: *
54: * LAMBDA (input) DOUBLE PRECISION
55: * The shift. In order to compute an accurate eigenvector,
56: * LAMBDA should be a good approximation to an eigenvalue
57: * of L D L^T.
58: *
59: * L (input) DOUBLE PRECISION array, dimension (N-1)
60: * The (n-1) subdiagonal elements of the unit bidiagonal matrix
61: * L, in elements 1 to N-1.
62: *
63: * D (input) DOUBLE PRECISION array, dimension (N)
64: * The n diagonal elements of the diagonal matrix D.
65: *
66: * LD (input) DOUBLE PRECISION array, dimension (N-1)
67: * The n-1 elements L(i)*D(i).
68: *
69: * LLD (input) DOUBLE PRECISION array, dimension (N-1)
70: * The n-1 elements L(i)*L(i)*D(i).
71: *
72: * PIVMIN (input) DOUBLE PRECISION
73: * The minimum pivot in the Sturm sequence.
74: *
75: * GAPTOL (input) DOUBLE PRECISION
76: * Tolerance that indicates when eigenvector entries are negligible
77: * w.r.t. their contribution to the residual.
78: *
79: * Z (input/output) COMPLEX*16 array, dimension (N)
80: * On input, all entries of Z must be set to 0.
81: * On output, Z contains the (scaled) r-th column of the
82: * inverse. The scaling is such that Z(R) equals 1.
83: *
84: * WANTNC (input) LOGICAL
85: * Specifies whether NEGCNT has to be computed.
86: *
87: * NEGCNT (output) INTEGER
88: * If WANTNC is .TRUE. then NEGCNT = the number of pivots < pivmin
89: * in the matrix factorization L D L^T, and NEGCNT = -1 otherwise.
90: *
91: * ZTZ (output) DOUBLE PRECISION
92: * The square of the 2-norm of Z.
93: *
94: * MINGMA (output) DOUBLE PRECISION
95: * The reciprocal of the largest (in magnitude) diagonal
96: * element of the inverse of L D L^T - sigma I.
97: *
98: * R (input/output) INTEGER
99: * 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) DOUBLE PRECISION
113: * NRMINV = 1/SQRT( ZTZ )
114: *
115: * RESID (output) DOUBLE PRECISION
116: * The residual of the FP vector.
117: * RESID = ABS( MINGMA )/SQRT( ZTZ )
118: *
119: * RQCORR (output) DOUBLE PRECISION
120: * The Rayleigh Quotient correction to LAMBDA.
121: * RQCORR = MINGMA*TMP
122: *
123: * WORK (workspace) DOUBLE PRECISION 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: DOUBLE PRECISION ZERO, ONE
139: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
140: COMPLEX*16 CONE
141: PARAMETER ( CONE = ( 1.0D0, 0.0D0 ) )
142:
143: * ..
144: * .. Local Scalars ..
145: LOGICAL SAWNAN1, SAWNAN2
146: INTEGER I, INDLPL, INDP, INDS, INDUMN, NEG1, NEG2, R1,
147: $ R2
148: DOUBLE PRECISION DMINUS, DPLUS, EPS, S, TMP
149: * ..
150: * .. External Functions ..
151: LOGICAL DISNAN
152: DOUBLE PRECISION DLAMCH
153: EXTERNAL DISNAN, DLAMCH
154: * ..
155: * .. Intrinsic Functions ..
156: INTRINSIC ABS, DBLE
157: * ..
158: * .. Executable Statements ..
159: *
160: EPS = DLAMCH( 'Precision' )
161:
162:
163: IF( R.EQ.0 ) THEN
164: R1 = B1
165: R2 = BN
166: ELSE
167: R1 = R
168: R2 = R
169: END IF
170:
171: * Storage for LPLUS
172: INDLPL = 0
173: * Storage for UMINUS
174: INDUMN = N
175: INDS = 2*N + 1
176: INDP = 3*N + 1
177:
178: IF( B1.EQ.1 ) THEN
179: WORK( INDS ) = ZERO
180: ELSE
181: WORK( INDS+B1-1 ) = LLD( B1-1 )
182: END IF
183:
184: *
185: * Compute the stationary transform (using the differential form)
186: * until the index R2.
187: *
188: SAWNAN1 = .FALSE.
189: NEG1 = 0
190: S = WORK( INDS+B1-1 ) - LAMBDA
191: DO 50 I = B1, R1 - 1
192: DPLUS = D( I ) + S
193: WORK( INDLPL+I ) = LD( I ) / DPLUS
194: IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
195: WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
196: S = WORK( INDS+I ) - LAMBDA
197: 50 CONTINUE
198: SAWNAN1 = DISNAN( S )
199: IF( SAWNAN1 ) GOTO 60
200: DO 51 I = R1, R2 - 1
201: DPLUS = D( I ) + S
202: WORK( INDLPL+I ) = LD( I ) / DPLUS
203: WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
204: S = WORK( INDS+I ) - LAMBDA
205: 51 CONTINUE
206: SAWNAN1 = DISNAN( S )
207: *
208: 60 CONTINUE
209: IF( SAWNAN1 ) THEN
210: * Runs a slower version of the above loop if a NaN is detected
211: NEG1 = 0
212: S = WORK( INDS+B1-1 ) - LAMBDA
213: DO 70 I = B1, R1 - 1
214: DPLUS = D( I ) + S
215: IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
216: WORK( INDLPL+I ) = LD( I ) / DPLUS
217: IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
218: WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
219: IF( WORK( INDLPL+I ).EQ.ZERO )
220: $ WORK( INDS+I ) = LLD( I )
221: S = WORK( INDS+I ) - LAMBDA
222: 70 CONTINUE
223: DO 71 I = R1, R2 - 1
224: DPLUS = D( I ) + S
225: IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
226: WORK( INDLPL+I ) = LD( I ) / DPLUS
227: WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
228: IF( WORK( INDLPL+I ).EQ.ZERO )
229: $ WORK( INDS+I ) = LLD( I )
230: S = WORK( INDS+I ) - LAMBDA
231: 71 CONTINUE
232: END IF
233: *
234: * Compute the progressive transform (using the differential form)
235: * until the index R1
236: *
237: SAWNAN2 = .FALSE.
238: NEG2 = 0
239: WORK( INDP+BN-1 ) = D( BN ) - LAMBDA
240: DO 80 I = BN - 1, R1, -1
241: DMINUS = LLD( I ) + WORK( INDP+I )
242: TMP = D( I ) / DMINUS
243: IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
244: WORK( INDUMN+I ) = L( I )*TMP
245: WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
246: 80 CONTINUE
247: TMP = WORK( INDP+R1-1 )
248: SAWNAN2 = DISNAN( TMP )
249:
250: IF( SAWNAN2 ) THEN
251: * Runs a slower version of the above loop if a NaN is detected
252: NEG2 = 0
253: DO 100 I = BN-1, R1, -1
254: DMINUS = LLD( I ) + WORK( INDP+I )
255: IF(ABS(DMINUS).LT.PIVMIN) DMINUS = -PIVMIN
256: TMP = D( I ) / DMINUS
257: IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
258: WORK( INDUMN+I ) = L( I )*TMP
259: WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
260: IF( TMP.EQ.ZERO )
261: $ WORK( INDP+I-1 ) = D( I ) - LAMBDA
262: 100 CONTINUE
263: END IF
264: *
265: * Find the index (from R1 to R2) of the largest (in magnitude)
266: * diagonal element of the inverse
267: *
268: MINGMA = WORK( INDS+R1-1 ) + WORK( INDP+R1-1 )
269: IF( MINGMA.LT.ZERO ) NEG1 = NEG1 + 1
270: IF( WANTNC ) THEN
271: NEGCNT = NEG1 + NEG2
272: ELSE
273: NEGCNT = -1
274: ENDIF
275: IF( ABS(MINGMA).EQ.ZERO )
276: $ MINGMA = EPS*WORK( INDS+R1-1 )
277: R = R1
278: DO 110 I = R1, R2 - 1
279: TMP = WORK( INDS+I ) + WORK( INDP+I )
280: IF( TMP.EQ.ZERO )
281: $ TMP = EPS*WORK( INDS+I )
282: IF( ABS( TMP ).LE.ABS( MINGMA ) ) THEN
283: MINGMA = TMP
284: R = I + 1
285: END IF
286: 110 CONTINUE
287: *
288: * Compute the FP vector: solve N^T v = e_r
289: *
290: ISUPPZ( 1 ) = B1
291: ISUPPZ( 2 ) = BN
292: Z( R ) = CONE
293: ZTZ = ONE
294: *
295: * Compute the FP vector upwards from R
296: *
297: IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
298: DO 210 I = R-1, B1, -1
299: Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
300: IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
301: $ THEN
302: Z( I ) = ZERO
303: ISUPPZ( 1 ) = I + 1
304: GOTO 220
305: ENDIF
306: ZTZ = ZTZ + DBLE( Z( I )*Z( I ) )
307: 210 CONTINUE
308: 220 CONTINUE
309: ELSE
310: * Run slower loop if NaN occurred.
311: DO 230 I = R - 1, B1, -1
312: IF( Z( I+1 ).EQ.ZERO ) THEN
313: Z( I ) = -( LD( I+1 ) / LD( I ) )*Z( I+2 )
314: ELSE
315: Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
316: END IF
317: IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
318: $ THEN
319: Z( I ) = ZERO
320: ISUPPZ( 1 ) = I + 1
321: GO TO 240
322: END IF
323: ZTZ = ZTZ + DBLE( Z( I )*Z( I ) )
324: 230 CONTINUE
325: 240 CONTINUE
326: ENDIF
327:
328: * Compute the FP vector downwards from R in blocks of size BLKSIZ
329: IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
330: DO 250 I = R, BN-1
331: Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
332: IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
333: $ THEN
334: Z( I+1 ) = ZERO
335: ISUPPZ( 2 ) = I
336: GO TO 260
337: END IF
338: ZTZ = ZTZ + DBLE( Z( I+1 )*Z( I+1 ) )
339: 250 CONTINUE
340: 260 CONTINUE
341: ELSE
342: * Run slower loop if NaN occurred.
343: DO 270 I = R, BN - 1
344: IF( Z( I ).EQ.ZERO ) THEN
345: Z( I+1 ) = -( LD( I-1 ) / LD( I ) )*Z( I-1 )
346: ELSE
347: Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
348: END IF
349: IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
350: $ THEN
351: Z( I+1 ) = ZERO
352: ISUPPZ( 2 ) = I
353: GO TO 280
354: END IF
355: ZTZ = ZTZ + DBLE( Z( I+1 )*Z( I+1 ) )
356: 270 CONTINUE
357: 280 CONTINUE
358: END IF
359: *
360: * Compute quantities for convergence test
361: *
362: TMP = ONE / ZTZ
363: NRMINV = SQRT( TMP )
364: RESID = ABS( MINGMA )*NRMINV
365: RQCORR = MINGMA*TMP
366: *
367: *
368: RETURN
369: *
370: * End of ZLAR1V
371: *
372: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlar1v.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack