1: *> \brief \b ZLAR1V computes the (scaled) r-th column of the inverse of the submatrix in rows b1 through bn of the tridiagonal matrix LDLT - λI.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZLAR1V + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlar1v.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlar1v.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlar1v.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZLAR1V( N, B1, BN, LAMBDA, D, L, LD, LLD,
22: * PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA,
23: * R, ISUPPZ, NRMINV, RESID, RQCORR, WORK )
24: *
25: * .. Scalar Arguments ..
26: * LOGICAL WANTNC
27: * INTEGER B1, BN, N, NEGCNT, R
28: * DOUBLE PRECISION GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID,
29: * $ RQCORR, ZTZ
30: * ..
31: * .. Array Arguments ..
32: * INTEGER ISUPPZ( * )
33: * DOUBLE PRECISION D( * ), L( * ), LD( * ), LLD( * ),
34: * $ WORK( * )
35: * COMPLEX*16 Z( * )
36: * ..
37: *
38: *
39: *> \par Purpose:
40: * =============
41: *>
42: *> \verbatim
43: *>
44: *> ZLAR1V computes the (scaled) r-th column of the inverse of
45: *> the sumbmatrix in rows B1 through BN of the tridiagonal matrix
46: *> L D L**T - sigma I. When sigma is close to an eigenvalue, the
47: *> computed vector is an accurate eigenvector. Usually, r corresponds
48: *> to the index where the eigenvector is largest in magnitude.
49: *> The following steps accomplish this computation :
50: *> (a) Stationary qd transform, L D L**T - sigma I = L(+) D(+) L(+)**T,
51: *> (b) Progressive qd transform, L D L**T - sigma I = U(-) D(-) U(-)**T,
52: *> (c) Computation of the diagonal elements of the inverse of
53: *> L D L**T - sigma I by combining the above transforms, and choosing
54: *> r as the index where the diagonal of the inverse is (one of the)
55: *> largest in magnitude.
56: *> (d) Computation of the (scaled) r-th column of the inverse using the
57: *> twisted factorization obtained by combining the top part of the
58: *> the stationary and the bottom part of the progressive transform.
59: *> \endverbatim
60: *
61: * Arguments:
62: * ==========
63: *
64: *> \param[in] N
65: *> \verbatim
66: *> N is INTEGER
67: *> The order of the matrix L D L**T.
68: *> \endverbatim
69: *>
70: *> \param[in] B1
71: *> \verbatim
72: *> B1 is INTEGER
73: *> First index of the submatrix of L D L**T.
74: *> \endverbatim
75: *>
76: *> \param[in] BN
77: *> \verbatim
78: *> BN is INTEGER
79: *> Last index of the submatrix of L D L**T.
80: *> \endverbatim
81: *>
82: *> \param[in] LAMBDA
83: *> \verbatim
84: *> LAMBDA is DOUBLE PRECISION
85: *> The shift. In order to compute an accurate eigenvector,
86: *> LAMBDA should be a good approximation to an eigenvalue
87: *> of L D L**T.
88: *> \endverbatim
89: *>
90: *> \param[in] L
91: *> \verbatim
92: *> L is DOUBLE PRECISION array, dimension (N-1)
93: *> The (n-1) subdiagonal elements of the unit bidiagonal matrix
94: *> L, in elements 1 to N-1.
95: *> \endverbatim
96: *>
97: *> \param[in] D
98: *> \verbatim
99: *> D is DOUBLE PRECISION array, dimension (N)
100: *> The n diagonal elements of the diagonal matrix D.
101: *> \endverbatim
102: *>
103: *> \param[in] LD
104: *> \verbatim
105: *> LD is DOUBLE PRECISION array, dimension (N-1)
106: *> The n-1 elements L(i)*D(i).
107: *> \endverbatim
108: *>
109: *> \param[in] LLD
110: *> \verbatim
111: *> LLD is DOUBLE PRECISION array, dimension (N-1)
112: *> The n-1 elements L(i)*L(i)*D(i).
113: *> \endverbatim
114: *>
115: *> \param[in] PIVMIN
116: *> \verbatim
117: *> PIVMIN is DOUBLE PRECISION
118: *> The minimum pivot in the Sturm sequence.
119: *> \endverbatim
120: *>
121: *> \param[in] GAPTOL
122: *> \verbatim
123: *> GAPTOL is DOUBLE PRECISION
124: *> Tolerance that indicates when eigenvector entries are negligible
125: *> w.r.t. their contribution to the residual.
126: *> \endverbatim
127: *>
128: *> \param[in,out] Z
129: *> \verbatim
130: *> Z is COMPLEX*16 array, dimension (N)
131: *> On input, all entries of Z must be set to 0.
132: *> On output, Z contains the (scaled) r-th column of the
133: *> inverse. The scaling is such that Z(R) equals 1.
134: *> \endverbatim
135: *>
136: *> \param[in] WANTNC
137: *> \verbatim
138: *> WANTNC is LOGICAL
139: *> Specifies whether NEGCNT has to be computed.
140: *> \endverbatim
141: *>
142: *> \param[out] NEGCNT
143: *> \verbatim
144: *> NEGCNT is INTEGER
145: *> If WANTNC is .TRUE. then NEGCNT = the number of pivots < pivmin
146: *> in the matrix factorization L D L**T, and NEGCNT = -1 otherwise.
147: *> \endverbatim
148: *>
149: *> \param[out] ZTZ
150: *> \verbatim
151: *> ZTZ is DOUBLE PRECISION
152: *> The square of the 2-norm of Z.
153: *> \endverbatim
154: *>
155: *> \param[out] MINGMA
156: *> \verbatim
157: *> MINGMA is DOUBLE PRECISION
158: *> The reciprocal of the largest (in magnitude) diagonal
159: *> element of the inverse of L D L**T - sigma I.
160: *> \endverbatim
161: *>
162: *> \param[in,out] R
163: *> \verbatim
164: *> R is INTEGER
165: *> The twist index for the twisted factorization used to
166: *> compute Z.
167: *> On input, 0 <= R <= N. If R is input as 0, R is set to
168: *> the index where (L D L**T - sigma I)^{-1} is largest
169: *> in magnitude. If 1 <= R <= N, R is unchanged.
170: *> On output, R contains the twist index used to compute Z.
171: *> Ideally, R designates the position of the maximum entry in the
172: *> eigenvector.
173: *> \endverbatim
174: *>
175: *> \param[out] ISUPPZ
176: *> \verbatim
177: *> ISUPPZ is INTEGER array, dimension (2)
178: *> The support of the vector in Z, i.e., the vector Z is
179: *> nonzero only in elements ISUPPZ(1) through ISUPPZ( 2 ).
180: *> \endverbatim
181: *>
182: *> \param[out] NRMINV
183: *> \verbatim
184: *> NRMINV is DOUBLE PRECISION
185: *> NRMINV = 1/SQRT( ZTZ )
186: *> \endverbatim
187: *>
188: *> \param[out] RESID
189: *> \verbatim
190: *> RESID is DOUBLE PRECISION
191: *> The residual of the FP vector.
192: *> RESID = ABS( MINGMA )/SQRT( ZTZ )
193: *> \endverbatim
194: *>
195: *> \param[out] RQCORR
196: *> \verbatim
197: *> RQCORR is DOUBLE PRECISION
198: *> The Rayleigh Quotient correction to LAMBDA.
199: *> RQCORR = MINGMA*TMP
200: *> \endverbatim
201: *>
202: *> \param[out] WORK
203: *> \verbatim
204: *> WORK is DOUBLE PRECISION array, dimension (4*N)
205: *> \endverbatim
206: *
207: * Authors:
208: * ========
209: *
210: *> \author Univ. of Tennessee
211: *> \author Univ. of California Berkeley
212: *> \author Univ. of Colorado Denver
213: *> \author NAG Ltd.
214: *
215: *> \ingroup complex16OTHERauxiliary
216: *
217: *> \par Contributors:
218: * ==================
219: *>
220: *> Beresford Parlett, University of California, Berkeley, USA \n
221: *> Jim Demmel, University of California, Berkeley, USA \n
222: *> Inderjit Dhillon, University of Texas, Austin, USA \n
223: *> Osni Marques, LBNL/NERSC, USA \n
224: *> Christof Voemel, University of California, Berkeley, USA
225: *
226: * =====================================================================
227: SUBROUTINE ZLAR1V( N, B1, BN, LAMBDA, D, L, LD, LLD,
228: $ PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA,
229: $ R, ISUPPZ, NRMINV, RESID, RQCORR, WORK )
230: *
231: * -- LAPACK auxiliary routine --
232: * -- LAPACK is a software package provided by Univ. of Tennessee, --
233: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
234: *
235: * .. Scalar Arguments ..
236: LOGICAL WANTNC
237: INTEGER B1, BN, N, NEGCNT, R
238: DOUBLE PRECISION GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID,
239: $ RQCORR, ZTZ
240: * ..
241: * .. Array Arguments ..
242: INTEGER ISUPPZ( * )
243: DOUBLE PRECISION D( * ), L( * ), LD( * ), LLD( * ),
244: $ WORK( * )
245: COMPLEX*16 Z( * )
246: * ..
247: *
248: * =====================================================================
249: *
250: * .. Parameters ..
251: DOUBLE PRECISION ZERO, ONE
252: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
253: COMPLEX*16 CONE
254: PARAMETER ( CONE = ( 1.0D0, 0.0D0 ) )
255:
256: * ..
257: * .. Local Scalars ..
258: LOGICAL SAWNAN1, SAWNAN2
259: INTEGER I, INDLPL, INDP, INDS, INDUMN, NEG1, NEG2, R1,
260: $ R2
261: DOUBLE PRECISION DMINUS, DPLUS, EPS, S, TMP
262: * ..
263: * .. External Functions ..
264: LOGICAL DISNAN
265: DOUBLE PRECISION DLAMCH
266: EXTERNAL DISNAN, DLAMCH
267: * ..
268: * .. Intrinsic Functions ..
269: INTRINSIC ABS, DBLE
270: * ..
271: * .. Executable Statements ..
272: *
273: EPS = DLAMCH( 'Precision' )
274:
275:
276: IF( R.EQ.0 ) THEN
277: R1 = B1
278: R2 = BN
279: ELSE
280: R1 = R
281: R2 = R
282: END IF
283:
284: * Storage for LPLUS
285: INDLPL = 0
286: * Storage for UMINUS
287: INDUMN = N
288: INDS = 2*N + 1
289: INDP = 3*N + 1
290:
291: IF( B1.EQ.1 ) THEN
292: WORK( INDS ) = ZERO
293: ELSE
294: WORK( INDS+B1-1 ) = LLD( B1-1 )
295: END IF
296:
297: *
298: * Compute the stationary transform (using the differential form)
299: * until the index R2.
300: *
301: SAWNAN1 = .FALSE.
302: NEG1 = 0
303: S = WORK( INDS+B1-1 ) - LAMBDA
304: DO 50 I = B1, R1 - 1
305: DPLUS = D( I ) + S
306: WORK( INDLPL+I ) = LD( I ) / DPLUS
307: IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
308: WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
309: S = WORK( INDS+I ) - LAMBDA
310: 50 CONTINUE
311: SAWNAN1 = DISNAN( S )
312: IF( SAWNAN1 ) GOTO 60
313: DO 51 I = R1, R2 - 1
314: DPLUS = D( I ) + S
315: WORK( INDLPL+I ) = LD( I ) / DPLUS
316: WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
317: S = WORK( INDS+I ) - LAMBDA
318: 51 CONTINUE
319: SAWNAN1 = DISNAN( S )
320: *
321: 60 CONTINUE
322: IF( SAWNAN1 ) THEN
323: * Runs a slower version of the above loop if a NaN is detected
324: NEG1 = 0
325: S = WORK( INDS+B1-1 ) - LAMBDA
326: DO 70 I = B1, R1 - 1
327: DPLUS = D( I ) + S
328: IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
329: WORK( INDLPL+I ) = LD( I ) / DPLUS
330: IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
331: WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
332: IF( WORK( INDLPL+I ).EQ.ZERO )
333: $ WORK( INDS+I ) = LLD( I )
334: S = WORK( INDS+I ) - LAMBDA
335: 70 CONTINUE
336: DO 71 I = R1, R2 - 1
337: DPLUS = D( I ) + S
338: IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
339: WORK( INDLPL+I ) = LD( I ) / DPLUS
340: WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
341: IF( WORK( INDLPL+I ).EQ.ZERO )
342: $ WORK( INDS+I ) = LLD( I )
343: S = WORK( INDS+I ) - LAMBDA
344: 71 CONTINUE
345: END IF
346: *
347: * Compute the progressive transform (using the differential form)
348: * until the index R1
349: *
350: SAWNAN2 = .FALSE.
351: NEG2 = 0
352: WORK( INDP+BN-1 ) = D( BN ) - LAMBDA
353: DO 80 I = BN - 1, R1, -1
354: DMINUS = LLD( I ) + WORK( INDP+I )
355: TMP = D( I ) / DMINUS
356: IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
357: WORK( INDUMN+I ) = L( I )*TMP
358: WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
359: 80 CONTINUE
360: TMP = WORK( INDP+R1-1 )
361: SAWNAN2 = DISNAN( TMP )
362:
363: IF( SAWNAN2 ) THEN
364: * Runs a slower version of the above loop if a NaN is detected
365: NEG2 = 0
366: DO 100 I = BN-1, R1, -1
367: DMINUS = LLD( I ) + WORK( INDP+I )
368: IF(ABS(DMINUS).LT.PIVMIN) DMINUS = -PIVMIN
369: TMP = D( I ) / DMINUS
370: IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
371: WORK( INDUMN+I ) = L( I )*TMP
372: WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
373: IF( TMP.EQ.ZERO )
374: $ WORK( INDP+I-1 ) = D( I ) - LAMBDA
375: 100 CONTINUE
376: END IF
377: *
378: * Find the index (from R1 to R2) of the largest (in magnitude)
379: * diagonal element of the inverse
380: *
381: MINGMA = WORK( INDS+R1-1 ) + WORK( INDP+R1-1 )
382: IF( MINGMA.LT.ZERO ) NEG1 = NEG1 + 1
383: IF( WANTNC ) THEN
384: NEGCNT = NEG1 + NEG2
385: ELSE
386: NEGCNT = -1
387: ENDIF
388: IF( ABS(MINGMA).EQ.ZERO )
389: $ MINGMA = EPS*WORK( INDS+R1-1 )
390: R = R1
391: DO 110 I = R1, R2 - 1
392: TMP = WORK( INDS+I ) + WORK( INDP+I )
393: IF( TMP.EQ.ZERO )
394: $ TMP = EPS*WORK( INDS+I )
395: IF( ABS( TMP ).LE.ABS( MINGMA ) ) THEN
396: MINGMA = TMP
397: R = I + 1
398: END IF
399: 110 CONTINUE
400: *
401: * Compute the FP vector: solve N^T v = e_r
402: *
403: ISUPPZ( 1 ) = B1
404: ISUPPZ( 2 ) = BN
405: Z( R ) = CONE
406: ZTZ = ONE
407: *
408: * Compute the FP vector upwards from R
409: *
410: IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
411: DO 210 I = R-1, B1, -1
412: Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
413: IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
414: $ THEN
415: Z( I ) = ZERO
416: ISUPPZ( 1 ) = I + 1
417: GOTO 220
418: ENDIF
419: ZTZ = ZTZ + DBLE( Z( I )*Z( I ) )
420: 210 CONTINUE
421: 220 CONTINUE
422: ELSE
423: * Run slower loop if NaN occurred.
424: DO 230 I = R - 1, B1, -1
425: IF( Z( I+1 ).EQ.ZERO ) THEN
426: Z( I ) = -( LD( I+1 ) / LD( I ) )*Z( I+2 )
427: ELSE
428: Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
429: END IF
430: IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
431: $ THEN
432: Z( I ) = ZERO
433: ISUPPZ( 1 ) = I + 1
434: GO TO 240
435: END IF
436: ZTZ = ZTZ + DBLE( Z( I )*Z( I ) )
437: 230 CONTINUE
438: 240 CONTINUE
439: ENDIF
440:
441: * Compute the FP vector downwards from R in blocks of size BLKSIZ
442: IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
443: DO 250 I = R, BN-1
444: Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
445: IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
446: $ THEN
447: Z( I+1 ) = ZERO
448: ISUPPZ( 2 ) = I
449: GO TO 260
450: END IF
451: ZTZ = ZTZ + DBLE( Z( I+1 )*Z( I+1 ) )
452: 250 CONTINUE
453: 260 CONTINUE
454: ELSE
455: * Run slower loop if NaN occurred.
456: DO 270 I = R, BN - 1
457: IF( Z( I ).EQ.ZERO ) THEN
458: Z( I+1 ) = -( LD( I-1 ) / LD( I ) )*Z( I-1 )
459: ELSE
460: Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
461: END IF
462: IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
463: $ THEN
464: Z( I+1 ) = ZERO
465: ISUPPZ( 2 ) = I
466: GO TO 280
467: END IF
468: ZTZ = ZTZ + DBLE( Z( I+1 )*Z( I+1 ) )
469: 270 CONTINUE
470: 280 CONTINUE
471: END IF
472: *
473: * Compute quantities for convergence test
474: *
475: TMP = ONE / ZTZ
476: NRMINV = SQRT( TMP )
477: RESID = ABS( MINGMA )*NRMINV
478: RQCORR = MINGMA*TMP
479: *
480: *
481: RETURN
482: *
483: * End of ZLAR1V
484: *
485: 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