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: *> \date September 2012
216: *
217: *> \ingroup complex16OTHERauxiliary
218: *
219: *> \par Contributors:
220: * ==================
221: *>
222: *> Beresford Parlett, University of California, Berkeley, USA \n
223: *> Jim Demmel, University of California, Berkeley, USA \n
224: *> Inderjit Dhillon, University of Texas, Austin, USA \n
225: *> Osni Marques, LBNL/NERSC, USA \n
226: *> Christof Voemel, University of California, Berkeley, USA
227: *
228: * =====================================================================
229: SUBROUTINE ZLAR1V( N, B1, BN, LAMBDA, D, L, LD, LLD,
230: $ PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA,
231: $ R, ISUPPZ, NRMINV, RESID, RQCORR, WORK )
232: *
233: * -- LAPACK auxiliary routine (version 3.4.2) --
234: * -- LAPACK is a software package provided by Univ. of Tennessee, --
235: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
236: * September 2012
237: *
238: * .. Scalar Arguments ..
239: LOGICAL WANTNC
240: INTEGER B1, BN, N, NEGCNT, R
241: DOUBLE PRECISION GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID,
242: $ RQCORR, ZTZ
243: * ..
244: * .. Array Arguments ..
245: INTEGER ISUPPZ( * )
246: DOUBLE PRECISION D( * ), L( * ), LD( * ), LLD( * ),
247: $ WORK( * )
248: COMPLEX*16 Z( * )
249: * ..
250: *
251: * =====================================================================
252: *
253: * .. Parameters ..
254: DOUBLE PRECISION ZERO, ONE
255: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
256: COMPLEX*16 CONE
257: PARAMETER ( CONE = ( 1.0D0, 0.0D0 ) )
258:
259: * ..
260: * .. Local Scalars ..
261: LOGICAL SAWNAN1, SAWNAN2
262: INTEGER I, INDLPL, INDP, INDS, INDUMN, NEG1, NEG2, R1,
263: $ R2
264: DOUBLE PRECISION DMINUS, DPLUS, EPS, S, TMP
265: * ..
266: * .. External Functions ..
267: LOGICAL DISNAN
268: DOUBLE PRECISION DLAMCH
269: EXTERNAL DISNAN, DLAMCH
270: * ..
271: * .. Intrinsic Functions ..
272: INTRINSIC ABS, DBLE
273: * ..
274: * .. Executable Statements ..
275: *
276: EPS = DLAMCH( 'Precision' )
277:
278:
279: IF( R.EQ.0 ) THEN
280: R1 = B1
281: R2 = BN
282: ELSE
283: R1 = R
284: R2 = R
285: END IF
286:
287: * Storage for LPLUS
288: INDLPL = 0
289: * Storage for UMINUS
290: INDUMN = N
291: INDS = 2*N + 1
292: INDP = 3*N + 1
293:
294: IF( B1.EQ.1 ) THEN
295: WORK( INDS ) = ZERO
296: ELSE
297: WORK( INDS+B1-1 ) = LLD( B1-1 )
298: END IF
299:
300: *
301: * Compute the stationary transform (using the differential form)
302: * until the index R2.
303: *
304: SAWNAN1 = .FALSE.
305: NEG1 = 0
306: S = WORK( INDS+B1-1 ) - LAMBDA
307: DO 50 I = B1, R1 - 1
308: DPLUS = D( I ) + S
309: WORK( INDLPL+I ) = LD( I ) / DPLUS
310: IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
311: WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
312: S = WORK( INDS+I ) - LAMBDA
313: 50 CONTINUE
314: SAWNAN1 = DISNAN( S )
315: IF( SAWNAN1 ) GOTO 60
316: DO 51 I = R1, R2 - 1
317: DPLUS = D( I ) + S
318: WORK( INDLPL+I ) = LD( I ) / DPLUS
319: WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
320: S = WORK( INDS+I ) - LAMBDA
321: 51 CONTINUE
322: SAWNAN1 = DISNAN( S )
323: *
324: 60 CONTINUE
325: IF( SAWNAN1 ) THEN
326: * Runs a slower version of the above loop if a NaN is detected
327: NEG1 = 0
328: S = WORK( INDS+B1-1 ) - LAMBDA
329: DO 70 I = B1, R1 - 1
330: DPLUS = D( I ) + S
331: IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
332: WORK( INDLPL+I ) = LD( I ) / DPLUS
333: IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
334: WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
335: IF( WORK( INDLPL+I ).EQ.ZERO )
336: $ WORK( INDS+I ) = LLD( I )
337: S = WORK( INDS+I ) - LAMBDA
338: 70 CONTINUE
339: DO 71 I = R1, R2 - 1
340: DPLUS = D( I ) + S
341: IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
342: WORK( INDLPL+I ) = LD( I ) / DPLUS
343: WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
344: IF( WORK( INDLPL+I ).EQ.ZERO )
345: $ WORK( INDS+I ) = LLD( I )
346: S = WORK( INDS+I ) - LAMBDA
347: 71 CONTINUE
348: END IF
349: *
350: * Compute the progressive transform (using the differential form)
351: * until the index R1
352: *
353: SAWNAN2 = .FALSE.
354: NEG2 = 0
355: WORK( INDP+BN-1 ) = D( BN ) - LAMBDA
356: DO 80 I = BN - 1, R1, -1
357: DMINUS = LLD( I ) + WORK( INDP+I )
358: TMP = D( I ) / DMINUS
359: IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
360: WORK( INDUMN+I ) = L( I )*TMP
361: WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
362: 80 CONTINUE
363: TMP = WORK( INDP+R1-1 )
364: SAWNAN2 = DISNAN( TMP )
365:
366: IF( SAWNAN2 ) THEN
367: * Runs a slower version of the above loop if a NaN is detected
368: NEG2 = 0
369: DO 100 I = BN-1, R1, -1
370: DMINUS = LLD( I ) + WORK( INDP+I )
371: IF(ABS(DMINUS).LT.PIVMIN) DMINUS = -PIVMIN
372: TMP = D( I ) / DMINUS
373: IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
374: WORK( INDUMN+I ) = L( I )*TMP
375: WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
376: IF( TMP.EQ.ZERO )
377: $ WORK( INDP+I-1 ) = D( I ) - LAMBDA
378: 100 CONTINUE
379: END IF
380: *
381: * Find the index (from R1 to R2) of the largest (in magnitude)
382: * diagonal element of the inverse
383: *
384: MINGMA = WORK( INDS+R1-1 ) + WORK( INDP+R1-1 )
385: IF( MINGMA.LT.ZERO ) NEG1 = NEG1 + 1
386: IF( WANTNC ) THEN
387: NEGCNT = NEG1 + NEG2
388: ELSE
389: NEGCNT = -1
390: ENDIF
391: IF( ABS(MINGMA).EQ.ZERO )
392: $ MINGMA = EPS*WORK( INDS+R1-1 )
393: R = R1
394: DO 110 I = R1, R2 - 1
395: TMP = WORK( INDS+I ) + WORK( INDP+I )
396: IF( TMP.EQ.ZERO )
397: $ TMP = EPS*WORK( INDS+I )
398: IF( ABS( TMP ).LE.ABS( MINGMA ) ) THEN
399: MINGMA = TMP
400: R = I + 1
401: END IF
402: 110 CONTINUE
403: *
404: * Compute the FP vector: solve N^T v = e_r
405: *
406: ISUPPZ( 1 ) = B1
407: ISUPPZ( 2 ) = BN
408: Z( R ) = CONE
409: ZTZ = ONE
410: *
411: * Compute the FP vector upwards from R
412: *
413: IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
414: DO 210 I = R-1, B1, -1
415: Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
416: IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
417: $ THEN
418: Z( I ) = ZERO
419: ISUPPZ( 1 ) = I + 1
420: GOTO 220
421: ENDIF
422: ZTZ = ZTZ + DBLE( Z( I )*Z( I ) )
423: 210 CONTINUE
424: 220 CONTINUE
425: ELSE
426: * Run slower loop if NaN occurred.
427: DO 230 I = R - 1, B1, -1
428: IF( Z( I+1 ).EQ.ZERO ) THEN
429: Z( I ) = -( LD( I+1 ) / LD( I ) )*Z( I+2 )
430: ELSE
431: Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
432: END IF
433: IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
434: $ THEN
435: Z( I ) = ZERO
436: ISUPPZ( 1 ) = I + 1
437: GO TO 240
438: END IF
439: ZTZ = ZTZ + DBLE( Z( I )*Z( I ) )
440: 230 CONTINUE
441: 240 CONTINUE
442: ENDIF
443:
444: * Compute the FP vector downwards from R in blocks of size BLKSIZ
445: IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
446: DO 250 I = R, BN-1
447: Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
448: IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
449: $ THEN
450: Z( I+1 ) = ZERO
451: ISUPPZ( 2 ) = I
452: GO TO 260
453: END IF
454: ZTZ = ZTZ + DBLE( Z( I+1 )*Z( I+1 ) )
455: 250 CONTINUE
456: 260 CONTINUE
457: ELSE
458: * Run slower loop if NaN occurred.
459: DO 270 I = R, BN - 1
460: IF( Z( I ).EQ.ZERO ) THEN
461: Z( I+1 ) = -( LD( I-1 ) / LD( I ) )*Z( I-1 )
462: ELSE
463: Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
464: END IF
465: IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
466: $ THEN
467: Z( I+1 ) = ZERO
468: ISUPPZ( 2 ) = I
469: GO TO 280
470: END IF
471: ZTZ = ZTZ + DBLE( Z( I+1 )*Z( I+1 ) )
472: 270 CONTINUE
473: 280 CONTINUE
474: END IF
475: *
476: * Compute quantities for convergence test
477: *
478: TMP = ONE / ZTZ
479: NRMINV = SQRT( TMP )
480: RESID = ABS( MINGMA )*NRMINV
481: RQCORR = MINGMA*TMP
482: *
483: *
484: RETURN
485: *
486: * End of ZLAR1V
487: *
488: 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