Annotation of rpl/lapack/lapack/zstein.f, revision 1.8
1.8 ! bertrand 1: *> \brief \b ZSTEIN
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download ZSTEIN + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zstein.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zstein.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zstein.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
! 22: * IWORK, IFAIL, INFO )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * INTEGER INFO, LDZ, M, N
! 26: * ..
! 27: * .. Array Arguments ..
! 28: * INTEGER IBLOCK( * ), IFAIL( * ), ISPLIT( * ),
! 29: * $ IWORK( * )
! 30: * DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
! 31: * COMPLEX*16 Z( LDZ, * )
! 32: * ..
! 33: *
! 34: *
! 35: *> \par Purpose:
! 36: * =============
! 37: *>
! 38: *> \verbatim
! 39: *>
! 40: *> ZSTEIN computes the eigenvectors of a real symmetric tridiagonal
! 41: *> matrix T corresponding to specified eigenvalues, using inverse
! 42: *> iteration.
! 43: *>
! 44: *> The maximum number of iterations allowed for each eigenvector is
! 45: *> specified by an internal parameter MAXITS (currently set to 5).
! 46: *>
! 47: *> Although the eigenvectors are real, they are stored in a complex
! 48: *> array, which may be passed to ZUNMTR or ZUPMTR for back
! 49: *> transformation to the eigenvectors of a complex Hermitian matrix
! 50: *> which was reduced to tridiagonal form.
! 51: *>
! 52: *> \endverbatim
! 53: *
! 54: * Arguments:
! 55: * ==========
! 56: *
! 57: *> \param[in] N
! 58: *> \verbatim
! 59: *> N is INTEGER
! 60: *> The order of the matrix. N >= 0.
! 61: *> \endverbatim
! 62: *>
! 63: *> \param[in] D
! 64: *> \verbatim
! 65: *> D is DOUBLE PRECISION array, dimension (N)
! 66: *> The n diagonal elements of the tridiagonal matrix T.
! 67: *> \endverbatim
! 68: *>
! 69: *> \param[in] E
! 70: *> \verbatim
! 71: *> E is DOUBLE PRECISION array, dimension (N-1)
! 72: *> The (n-1) subdiagonal elements of the tridiagonal matrix
! 73: *> T, stored in elements 1 to N-1.
! 74: *> \endverbatim
! 75: *>
! 76: *> \param[in] M
! 77: *> \verbatim
! 78: *> M is INTEGER
! 79: *> The number of eigenvectors to be found. 0 <= M <= N.
! 80: *> \endverbatim
! 81: *>
! 82: *> \param[in] W
! 83: *> \verbatim
! 84: *> W is DOUBLE PRECISION array, dimension (N)
! 85: *> The first M elements of W contain the eigenvalues for
! 86: *> which eigenvectors are to be computed. The eigenvalues
! 87: *> should be grouped by split-off block and ordered from
! 88: *> smallest to largest within the block. ( The output array
! 89: *> W from DSTEBZ with ORDER = 'B' is expected here. )
! 90: *> \endverbatim
! 91: *>
! 92: *> \param[in] IBLOCK
! 93: *> \verbatim
! 94: *> IBLOCK is INTEGER array, dimension (N)
! 95: *> The submatrix indices associated with the corresponding
! 96: *> eigenvalues in W; IBLOCK(i)=1 if eigenvalue W(i) belongs to
! 97: *> the first submatrix from the top, =2 if W(i) belongs to
! 98: *> the second submatrix, etc. ( The output array IBLOCK
! 99: *> from DSTEBZ is expected here. )
! 100: *> \endverbatim
! 101: *>
! 102: *> \param[in] ISPLIT
! 103: *> \verbatim
! 104: *> ISPLIT is INTEGER array, dimension (N)
! 105: *> The splitting points, at which T breaks up into submatrices.
! 106: *> The first submatrix consists of rows/columns 1 to
! 107: *> ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
! 108: *> through ISPLIT( 2 ), etc.
! 109: *> ( The output array ISPLIT from DSTEBZ is expected here. )
! 110: *> \endverbatim
! 111: *>
! 112: *> \param[out] Z
! 113: *> \verbatim
! 114: *> Z is COMPLEX*16 array, dimension (LDZ, M)
! 115: *> The computed eigenvectors. The eigenvector associated
! 116: *> with the eigenvalue W(i) is stored in the i-th column of
! 117: *> Z. Any vector which fails to converge is set to its current
! 118: *> iterate after MAXITS iterations.
! 119: *> The imaginary parts of the eigenvectors are set to zero.
! 120: *> \endverbatim
! 121: *>
! 122: *> \param[in] LDZ
! 123: *> \verbatim
! 124: *> LDZ is INTEGER
! 125: *> The leading dimension of the array Z. LDZ >= max(1,N).
! 126: *> \endverbatim
! 127: *>
! 128: *> \param[out] WORK
! 129: *> \verbatim
! 130: *> WORK is DOUBLE PRECISION array, dimension (5*N)
! 131: *> \endverbatim
! 132: *>
! 133: *> \param[out] IWORK
! 134: *> \verbatim
! 135: *> IWORK is INTEGER array, dimension (N)
! 136: *> \endverbatim
! 137: *>
! 138: *> \param[out] IFAIL
! 139: *> \verbatim
! 140: *> IFAIL is INTEGER array, dimension (M)
! 141: *> On normal exit, all elements of IFAIL are zero.
! 142: *> If one or more eigenvectors fail to converge after
! 143: *> MAXITS iterations, then their indices are stored in
! 144: *> array IFAIL.
! 145: *> \endverbatim
! 146: *>
! 147: *> \param[out] INFO
! 148: *> \verbatim
! 149: *> INFO is INTEGER
! 150: *> = 0: successful exit
! 151: *> < 0: if INFO = -i, the i-th argument had an illegal value
! 152: *> > 0: if INFO = i, then i eigenvectors failed to converge
! 153: *> in MAXITS iterations. Their indices are stored in
! 154: *> array IFAIL.
! 155: *> \endverbatim
! 156: *
! 157: *> \par Internal Parameters:
! 158: * =========================
! 159: *>
! 160: *> \verbatim
! 161: *> MAXITS INTEGER, default = 5
! 162: *> The maximum number of iterations performed.
! 163: *>
! 164: *> EXTRA INTEGER, default = 2
! 165: *> The number of iterations performed after norm growth
! 166: *> criterion is satisfied, should be at least 1.
! 167: *> \endverbatim
! 168: *
! 169: * Authors:
! 170: * ========
! 171: *
! 172: *> \author Univ. of Tennessee
! 173: *> \author Univ. of California Berkeley
! 174: *> \author Univ. of Colorado Denver
! 175: *> \author NAG Ltd.
! 176: *
! 177: *> \date November 2011
! 178: *
! 179: *> \ingroup complex16OTHERcomputational
! 180: *
! 181: * =====================================================================
1.1 bertrand 182: SUBROUTINE ZSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
183: $ IWORK, IFAIL, INFO )
184: *
1.8 ! bertrand 185: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 186: * -- LAPACK is a software package provided by Univ. of Tennessee, --
187: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8 ! bertrand 188: * November 2011
1.1 bertrand 189: *
190: * .. Scalar Arguments ..
191: INTEGER INFO, LDZ, M, N
192: * ..
193: * .. Array Arguments ..
194: INTEGER IBLOCK( * ), IFAIL( * ), ISPLIT( * ),
195: $ IWORK( * )
196: DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
197: COMPLEX*16 Z( LDZ, * )
198: * ..
199: *
200: * =====================================================================
201: *
202: * .. Parameters ..
203: COMPLEX*16 CZERO, CONE
204: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
205: $ CONE = ( 1.0D+0, 0.0D+0 ) )
206: DOUBLE PRECISION ZERO, ONE, TEN, ODM3, ODM1
207: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TEN = 1.0D+1,
208: $ ODM3 = 1.0D-3, ODM1 = 1.0D-1 )
209: INTEGER MAXITS, EXTRA
210: PARAMETER ( MAXITS = 5, EXTRA = 2 )
211: * ..
212: * .. Local Scalars ..
213: INTEGER B1, BLKSIZ, BN, GPIND, I, IINFO, INDRV1,
214: $ INDRV2, INDRV3, INDRV4, INDRV5, ITS, J, J1,
215: $ JBLK, JMAX, JR, NBLK, NRMCHK
216: DOUBLE PRECISION DTPCRT, EPS, EPS1, NRM, ONENRM, ORTOL, PERTOL,
217: $ SCL, SEP, TOL, XJ, XJM, ZTR
218: * ..
219: * .. Local Arrays ..
220: INTEGER ISEED( 4 )
221: * ..
222: * .. External Functions ..
223: INTEGER IDAMAX
224: DOUBLE PRECISION DASUM, DLAMCH, DNRM2
225: EXTERNAL IDAMAX, DASUM, DLAMCH, DNRM2
226: * ..
227: * .. External Subroutines ..
228: EXTERNAL DCOPY, DLAGTF, DLAGTS, DLARNV, DSCAL, XERBLA
229: * ..
230: * .. Intrinsic Functions ..
231: INTRINSIC ABS, DBLE, DCMPLX, MAX, SQRT
232: * ..
233: * .. Executable Statements ..
234: *
235: * Test the input parameters.
236: *
237: INFO = 0
238: DO 10 I = 1, M
239: IFAIL( I ) = 0
240: 10 CONTINUE
241: *
242: IF( N.LT.0 ) THEN
243: INFO = -1
244: ELSE IF( M.LT.0 .OR. M.GT.N ) THEN
245: INFO = -4
246: ELSE IF( LDZ.LT.MAX( 1, N ) ) THEN
247: INFO = -9
248: ELSE
249: DO 20 J = 2, M
250: IF( IBLOCK( J ).LT.IBLOCK( J-1 ) ) THEN
251: INFO = -6
252: GO TO 30
253: END IF
254: IF( IBLOCK( J ).EQ.IBLOCK( J-1 ) .AND. W( J ).LT.W( J-1 ) )
255: $ THEN
256: INFO = -5
257: GO TO 30
258: END IF
259: 20 CONTINUE
260: 30 CONTINUE
261: END IF
262: *
263: IF( INFO.NE.0 ) THEN
264: CALL XERBLA( 'ZSTEIN', -INFO )
265: RETURN
266: END IF
267: *
268: * Quick return if possible
269: *
270: IF( N.EQ.0 .OR. M.EQ.0 ) THEN
271: RETURN
272: ELSE IF( N.EQ.1 ) THEN
273: Z( 1, 1 ) = CONE
274: RETURN
275: END IF
276: *
277: * Get machine constants.
278: *
279: EPS = DLAMCH( 'Precision' )
280: *
281: * Initialize seed for random number generator DLARNV.
282: *
283: DO 40 I = 1, 4
284: ISEED( I ) = 1
285: 40 CONTINUE
286: *
287: * Initialize pointers.
288: *
289: INDRV1 = 0
290: INDRV2 = INDRV1 + N
291: INDRV3 = INDRV2 + N
292: INDRV4 = INDRV3 + N
293: INDRV5 = INDRV4 + N
294: *
295: * Compute eigenvectors of matrix blocks.
296: *
297: J1 = 1
298: DO 180 NBLK = 1, IBLOCK( M )
299: *
300: * Find starting and ending indices of block nblk.
301: *
302: IF( NBLK.EQ.1 ) THEN
303: B1 = 1
304: ELSE
305: B1 = ISPLIT( NBLK-1 ) + 1
306: END IF
307: BN = ISPLIT( NBLK )
308: BLKSIZ = BN - B1 + 1
309: IF( BLKSIZ.EQ.1 )
310: $ GO TO 60
311: GPIND = B1
312: *
313: * Compute reorthogonalization criterion and stopping criterion.
314: *
315: ONENRM = ABS( D( B1 ) ) + ABS( E( B1 ) )
316: ONENRM = MAX( ONENRM, ABS( D( BN ) )+ABS( E( BN-1 ) ) )
317: DO 50 I = B1 + 1, BN - 1
318: ONENRM = MAX( ONENRM, ABS( D( I ) )+ABS( E( I-1 ) )+
319: $ ABS( E( I ) ) )
320: 50 CONTINUE
321: ORTOL = ODM3*ONENRM
322: *
323: DTPCRT = SQRT( ODM1 / BLKSIZ )
324: *
325: * Loop through eigenvalues of block nblk.
326: *
327: 60 CONTINUE
328: JBLK = 0
329: DO 170 J = J1, M
330: IF( IBLOCK( J ).NE.NBLK ) THEN
331: J1 = J
332: GO TO 180
333: END IF
334: JBLK = JBLK + 1
335: XJ = W( J )
336: *
337: * Skip all the work if the block size is one.
338: *
339: IF( BLKSIZ.EQ.1 ) THEN
340: WORK( INDRV1+1 ) = ONE
341: GO TO 140
342: END IF
343: *
344: * If eigenvalues j and j-1 are too close, add a relatively
345: * small perturbation.
346: *
347: IF( JBLK.GT.1 ) THEN
348: EPS1 = ABS( EPS*XJ )
349: PERTOL = TEN*EPS1
350: SEP = XJ - XJM
351: IF( SEP.LT.PERTOL )
352: $ XJ = XJM + PERTOL
353: END IF
354: *
355: ITS = 0
356: NRMCHK = 0
357: *
358: * Get random starting vector.
359: *
360: CALL DLARNV( 2, ISEED, BLKSIZ, WORK( INDRV1+1 ) )
361: *
362: * Copy the matrix T so it won't be destroyed in factorization.
363: *
364: CALL DCOPY( BLKSIZ, D( B1 ), 1, WORK( INDRV4+1 ), 1 )
365: CALL DCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV2+2 ), 1 )
366: CALL DCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV3+1 ), 1 )
367: *
368: * Compute LU factors with partial pivoting ( PT = LU )
369: *
370: TOL = ZERO
371: CALL DLAGTF( BLKSIZ, WORK( INDRV4+1 ), XJ, WORK( INDRV2+2 ),
372: $ WORK( INDRV3+1 ), TOL, WORK( INDRV5+1 ), IWORK,
373: $ IINFO )
374: *
375: * Update iteration count.
376: *
377: 70 CONTINUE
378: ITS = ITS + 1
379: IF( ITS.GT.MAXITS )
380: $ GO TO 120
381: *
382: * Normalize and scale the righthand side vector Pb.
383: *
384: SCL = BLKSIZ*ONENRM*MAX( EPS,
385: $ ABS( WORK( INDRV4+BLKSIZ ) ) ) /
386: $ DASUM( BLKSIZ, WORK( INDRV1+1 ), 1 )
387: CALL DSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 )
388: *
389: * Solve the system LU = Pb.
390: *
391: CALL DLAGTS( -1, BLKSIZ, WORK( INDRV4+1 ), WORK( INDRV2+2 ),
392: $ WORK( INDRV3+1 ), WORK( INDRV5+1 ), IWORK,
393: $ WORK( INDRV1+1 ), TOL, IINFO )
394: *
395: * Reorthogonalize by modified Gram-Schmidt if eigenvalues are
396: * close enough.
397: *
398: IF( JBLK.EQ.1 )
399: $ GO TO 110
400: IF( ABS( XJ-XJM ).GT.ORTOL )
401: $ GPIND = J
402: IF( GPIND.NE.J ) THEN
403: DO 100 I = GPIND, J - 1
404: ZTR = ZERO
405: DO 80 JR = 1, BLKSIZ
406: ZTR = ZTR + WORK( INDRV1+JR )*
407: $ DBLE( Z( B1-1+JR, I ) )
408: 80 CONTINUE
409: DO 90 JR = 1, BLKSIZ
410: WORK( INDRV1+JR ) = WORK( INDRV1+JR ) -
411: $ ZTR*DBLE( Z( B1-1+JR, I ) )
412: 90 CONTINUE
413: 100 CONTINUE
414: END IF
415: *
416: * Check the infinity norm of the iterate.
417: *
418: 110 CONTINUE
419: JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 )
420: NRM = ABS( WORK( INDRV1+JMAX ) )
421: *
422: * Continue for additional iterations after norm reaches
423: * stopping criterion.
424: *
425: IF( NRM.LT.DTPCRT )
426: $ GO TO 70
427: NRMCHK = NRMCHK + 1
428: IF( NRMCHK.LT.EXTRA+1 )
429: $ GO TO 70
430: *
431: GO TO 130
432: *
433: * If stopping criterion was not satisfied, update info and
434: * store eigenvector number in array ifail.
435: *
436: 120 CONTINUE
437: INFO = INFO + 1
438: IFAIL( INFO ) = J
439: *
440: * Accept iterate as jth eigenvector.
441: *
442: 130 CONTINUE
443: SCL = ONE / DNRM2( BLKSIZ, WORK( INDRV1+1 ), 1 )
444: JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 )
445: IF( WORK( INDRV1+JMAX ).LT.ZERO )
446: $ SCL = -SCL
447: CALL DSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 )
448: 140 CONTINUE
449: DO 150 I = 1, N
450: Z( I, J ) = CZERO
451: 150 CONTINUE
452: DO 160 I = 1, BLKSIZ
453: Z( B1+I-1, J ) = DCMPLX( WORK( INDRV1+I ), ZERO )
454: 160 CONTINUE
455: *
456: * Save the shift to check eigenvalue spacing at next
457: * iteration.
458: *
459: XJM = XJ
460: *
461: 170 CONTINUE
462: 180 CONTINUE
463: *
464: RETURN
465: *
466: * End of ZSTEIN
467: *
468: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zstein.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack