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