Annotation of rpl/lapack/lapack/zlar1v.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b ZLAR1V
! 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 November 2011
! 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: * =====================================================================
1.1 bertrand 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: *
1.9 ! bertrand 233: * -- LAPACK auxiliary routine (version 3.4.0) --
1.1 bertrand 234: * -- LAPACK is a software package provided by Univ. of Tennessee, --
235: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 236: * November 2011
1.1 bertrand 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