1: *> \brief \b DLAR1V 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 DLAR1V + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlar1v.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlar1v.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlar1v.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLAR1V( 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: * DOUBLE PRECISION Z( * )
36: * ..
37: *
38: *
39: *> \par Purpose:
40: * =============
41: *>
42: *> \verbatim
43: *>
44: *> DLAR1V 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 DOUBLE PRECISION 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 doubleOTHERauxiliary
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 DLAR1V( 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: DOUBLE PRECISION Z( * )
246: * ..
247: *
248: * =====================================================================
249: *
250: * .. Parameters ..
251: DOUBLE PRECISION ZERO, ONE
252: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
253:
254: * ..
255: * .. Local Scalars ..
256: LOGICAL SAWNAN1, SAWNAN2
257: INTEGER I, INDLPL, INDP, INDS, INDUMN, NEG1, NEG2, R1,
258: $ R2
259: DOUBLE PRECISION DMINUS, DPLUS, EPS, S, TMP
260: * ..
261: * .. External Functions ..
262: LOGICAL DISNAN
263: DOUBLE PRECISION DLAMCH
264: EXTERNAL DISNAN, DLAMCH
265: * ..
266: * .. Intrinsic Functions ..
267: INTRINSIC ABS
268: * ..
269: * .. Executable Statements ..
270: *
271: EPS = DLAMCH( 'Precision' )
272:
273:
274: IF( R.EQ.0 ) THEN
275: R1 = B1
276: R2 = BN
277: ELSE
278: R1 = R
279: R2 = R
280: END IF
281:
282: * Storage for LPLUS
283: INDLPL = 0
284: * Storage for UMINUS
285: INDUMN = N
286: INDS = 2*N + 1
287: INDP = 3*N + 1
288:
289: IF( B1.EQ.1 ) THEN
290: WORK( INDS ) = ZERO
291: ELSE
292: WORK( INDS+B1-1 ) = LLD( B1-1 )
293: END IF
294:
295: *
296: * Compute the stationary transform (using the differential form)
297: * until the index R2.
298: *
299: SAWNAN1 = .FALSE.
300: NEG1 = 0
301: S = WORK( INDS+B1-1 ) - LAMBDA
302: DO 50 I = B1, R1 - 1
303: DPLUS = D( I ) + S
304: WORK( INDLPL+I ) = LD( I ) / DPLUS
305: IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
306: WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
307: S = WORK( INDS+I ) - LAMBDA
308: 50 CONTINUE
309: SAWNAN1 = DISNAN( S )
310: IF( SAWNAN1 ) GOTO 60
311: DO 51 I = R1, R2 - 1
312: DPLUS = D( I ) + S
313: WORK( INDLPL+I ) = LD( I ) / DPLUS
314: WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
315: S = WORK( INDS+I ) - LAMBDA
316: 51 CONTINUE
317: SAWNAN1 = DISNAN( S )
318: *
319: 60 CONTINUE
320: IF( SAWNAN1 ) THEN
321: * Runs a slower version of the above loop if a NaN is detected
322: NEG1 = 0
323: S = WORK( INDS+B1-1 ) - LAMBDA
324: DO 70 I = B1, R1 - 1
325: DPLUS = D( I ) + S
326: IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
327: WORK( INDLPL+I ) = LD( I ) / DPLUS
328: IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
329: WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
330: IF( WORK( INDLPL+I ).EQ.ZERO )
331: $ WORK( INDS+I ) = LLD( I )
332: S = WORK( INDS+I ) - LAMBDA
333: 70 CONTINUE
334: DO 71 I = R1, R2 - 1
335: DPLUS = D( I ) + S
336: IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
337: WORK( INDLPL+I ) = LD( I ) / DPLUS
338: WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
339: IF( WORK( INDLPL+I ).EQ.ZERO )
340: $ WORK( INDS+I ) = LLD( I )
341: S = WORK( INDS+I ) - LAMBDA
342: 71 CONTINUE
343: END IF
344: *
345: * Compute the progressive transform (using the differential form)
346: * until the index R1
347: *
348: SAWNAN2 = .FALSE.
349: NEG2 = 0
350: WORK( INDP+BN-1 ) = D( BN ) - LAMBDA
351: DO 80 I = BN - 1, R1, -1
352: DMINUS = LLD( I ) + WORK( INDP+I )
353: TMP = D( I ) / DMINUS
354: IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
355: WORK( INDUMN+I ) = L( I )*TMP
356: WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
357: 80 CONTINUE
358: TMP = WORK( INDP+R1-1 )
359: SAWNAN2 = DISNAN( TMP )
360:
361: IF( SAWNAN2 ) THEN
362: * Runs a slower version of the above loop if a NaN is detected
363: NEG2 = 0
364: DO 100 I = BN-1, R1, -1
365: DMINUS = LLD( I ) + WORK( INDP+I )
366: IF(ABS(DMINUS).LT.PIVMIN) DMINUS = -PIVMIN
367: TMP = D( I ) / DMINUS
368: IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
369: WORK( INDUMN+I ) = L( I )*TMP
370: WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
371: IF( TMP.EQ.ZERO )
372: $ WORK( INDP+I-1 ) = D( I ) - LAMBDA
373: 100 CONTINUE
374: END IF
375: *
376: * Find the index (from R1 to R2) of the largest (in magnitude)
377: * diagonal element of the inverse
378: *
379: MINGMA = WORK( INDS+R1-1 ) + WORK( INDP+R1-1 )
380: IF( MINGMA.LT.ZERO ) NEG1 = NEG1 + 1
381: IF( WANTNC ) THEN
382: NEGCNT = NEG1 + NEG2
383: ELSE
384: NEGCNT = -1
385: ENDIF
386: IF( ABS(MINGMA).EQ.ZERO )
387: $ MINGMA = EPS*WORK( INDS+R1-1 )
388: R = R1
389: DO 110 I = R1, R2 - 1
390: TMP = WORK( INDS+I ) + WORK( INDP+I )
391: IF( TMP.EQ.ZERO )
392: $ TMP = EPS*WORK( INDS+I )
393: IF( ABS( TMP ).LE.ABS( MINGMA ) ) THEN
394: MINGMA = TMP
395: R = I + 1
396: END IF
397: 110 CONTINUE
398: *
399: * Compute the FP vector: solve N^T v = e_r
400: *
401: ISUPPZ( 1 ) = B1
402: ISUPPZ( 2 ) = BN
403: Z( R ) = ONE
404: ZTZ = ONE
405: *
406: * Compute the FP vector upwards from R
407: *
408: IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
409: DO 210 I = R-1, B1, -1
410: Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
411: IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
412: $ THEN
413: Z( I ) = ZERO
414: ISUPPZ( 1 ) = I + 1
415: GOTO 220
416: ENDIF
417: ZTZ = ZTZ + Z( I )*Z( I )
418: 210 CONTINUE
419: 220 CONTINUE
420: ELSE
421: * Run slower loop if NaN occurred.
422: DO 230 I = R - 1, B1, -1
423: IF( Z( I+1 ).EQ.ZERO ) THEN
424: Z( I ) = -( LD( I+1 ) / LD( I ) )*Z( I+2 )
425: ELSE
426: Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
427: END IF
428: IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
429: $ THEN
430: Z( I ) = ZERO
431: ISUPPZ( 1 ) = I + 1
432: GO TO 240
433: END IF
434: ZTZ = ZTZ + Z( I )*Z( I )
435: 230 CONTINUE
436: 240 CONTINUE
437: ENDIF
438:
439: * Compute the FP vector downwards from R in blocks of size BLKSIZ
440: IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
441: DO 250 I = R, BN-1
442: Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
443: IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
444: $ THEN
445: Z( I+1 ) = ZERO
446: ISUPPZ( 2 ) = I
447: GO TO 260
448: END IF
449: ZTZ = ZTZ + Z( I+1 )*Z( I+1 )
450: 250 CONTINUE
451: 260 CONTINUE
452: ELSE
453: * Run slower loop if NaN occurred.
454: DO 270 I = R, BN - 1
455: IF( Z( I ).EQ.ZERO ) THEN
456: Z( I+1 ) = -( LD( I-1 ) / LD( I ) )*Z( I-1 )
457: ELSE
458: Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
459: END IF
460: IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
461: $ THEN
462: Z( I+1 ) = ZERO
463: ISUPPZ( 2 ) = I
464: GO TO 280
465: END IF
466: ZTZ = ZTZ + Z( I+1 )*Z( I+1 )
467: 270 CONTINUE
468: 280 CONTINUE
469: END IF
470: *
471: * Compute quantities for convergence test
472: *
473: TMP = ONE / ZTZ
474: NRMINV = SQRT( TMP )
475: RESID = ABS( MINGMA )*NRMINV
476: RQCORR = MINGMA*TMP
477: *
478: *
479: RETURN
480: *
481: * End of DLAR1V
482: *
483: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlar1v.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack