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: *> \ingroup complex16OTHERcomputational
178: *
179: * =====================================================================
180: SUBROUTINE ZSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
181: $ IWORK, IFAIL, INFO )
182: *
183: * -- LAPACK computational routine --
184: * -- LAPACK is a software package provided by Univ. of Tennessee, --
185: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
186: *
187: * .. Scalar Arguments ..
188: INTEGER INFO, LDZ, M, N
189: * ..
190: * .. Array Arguments ..
191: INTEGER IBLOCK( * ), IFAIL( * ), ISPLIT( * ),
192: $ IWORK( * )
193: DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
194: COMPLEX*16 Z( LDZ, * )
195: * ..
196: *
197: * =====================================================================
198: *
199: * .. Parameters ..
200: COMPLEX*16 CZERO, CONE
201: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
202: $ CONE = ( 1.0D+0, 0.0D+0 ) )
203: DOUBLE PRECISION ZERO, ONE, TEN, ODM3, ODM1
204: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TEN = 1.0D+1,
205: $ ODM3 = 1.0D-3, ODM1 = 1.0D-1 )
206: INTEGER MAXITS, EXTRA
207: PARAMETER ( MAXITS = 5, EXTRA = 2 )
208: * ..
209: * .. Local Scalars ..
210: INTEGER B1, BLKSIZ, BN, GPIND, I, IINFO, INDRV1,
211: $ INDRV2, INDRV3, INDRV4, INDRV5, ITS, J, J1,
212: $ JBLK, JMAX, JR, NBLK, NRMCHK
213: DOUBLE PRECISION DTPCRT, EPS, EPS1, NRM, ONENRM, ORTOL, PERTOL,
214: $ SCL, SEP, TOL, XJ, XJM, ZTR
215: * ..
216: * .. Local Arrays ..
217: INTEGER ISEED( 4 )
218: * ..
219: * .. External Functions ..
220: INTEGER IDAMAX
221: DOUBLE PRECISION DLAMCH, DNRM2
222: EXTERNAL IDAMAX, DLAMCH, DNRM2
223: * ..
224: * .. External Subroutines ..
225: EXTERNAL DCOPY, DLAGTF, DLAGTS, DLARNV, DSCAL, XERBLA
226: * ..
227: * .. Intrinsic Functions ..
228: INTRINSIC ABS, DBLE, DCMPLX, MAX, SQRT
229: * ..
230: * .. Executable Statements ..
231: *
232: * Test the input parameters.
233: *
234: INFO = 0
235: DO 10 I = 1, M
236: IFAIL( I ) = 0
237: 10 CONTINUE
238: *
239: IF( N.LT.0 ) THEN
240: INFO = -1
241: ELSE IF( M.LT.0 .OR. M.GT.N ) THEN
242: INFO = -4
243: ELSE IF( LDZ.LT.MAX( 1, N ) ) THEN
244: INFO = -9
245: ELSE
246: DO 20 J = 2, M
247: IF( IBLOCK( J ).LT.IBLOCK( J-1 ) ) THEN
248: INFO = -6
249: GO TO 30
250: END IF
251: IF( IBLOCK( J ).EQ.IBLOCK( J-1 ) .AND. W( J ).LT.W( J-1 ) )
252: $ THEN
253: INFO = -5
254: GO TO 30
255: END IF
256: 20 CONTINUE
257: 30 CONTINUE
258: END IF
259: *
260: IF( INFO.NE.0 ) THEN
261: CALL XERBLA( 'ZSTEIN', -INFO )
262: RETURN
263: END IF
264: *
265: * Quick return if possible
266: *
267: IF( N.EQ.0 .OR. M.EQ.0 ) THEN
268: RETURN
269: ELSE IF( N.EQ.1 ) THEN
270: Z( 1, 1 ) = CONE
271: RETURN
272: END IF
273: *
274: * Get machine constants.
275: *
276: EPS = DLAMCH( 'Precision' )
277: *
278: * Initialize seed for random number generator DLARNV.
279: *
280: DO 40 I = 1, 4
281: ISEED( I ) = 1
282: 40 CONTINUE
283: *
284: * Initialize pointers.
285: *
286: INDRV1 = 0
287: INDRV2 = INDRV1 + N
288: INDRV3 = INDRV2 + N
289: INDRV4 = INDRV3 + N
290: INDRV5 = INDRV4 + N
291: *
292: * Compute eigenvectors of matrix blocks.
293: *
294: J1 = 1
295: DO 180 NBLK = 1, IBLOCK( M )
296: *
297: * Find starting and ending indices of block nblk.
298: *
299: IF( NBLK.EQ.1 ) THEN
300: B1 = 1
301: ELSE
302: B1 = ISPLIT( NBLK-1 ) + 1
303: END IF
304: BN = ISPLIT( NBLK )
305: BLKSIZ = BN - B1 + 1
306: IF( BLKSIZ.EQ.1 )
307: $ GO TO 60
308: GPIND = J1
309: *
310: * Compute reorthogonalization criterion and stopping criterion.
311: *
312: ONENRM = ABS( D( B1 ) ) + ABS( E( B1 ) )
313: ONENRM = MAX( ONENRM, ABS( D( BN ) )+ABS( E( BN-1 ) ) )
314: DO 50 I = B1 + 1, BN - 1
315: ONENRM = MAX( ONENRM, ABS( D( I ) )+ABS( E( I-1 ) )+
316: $ ABS( E( I ) ) )
317: 50 CONTINUE
318: ORTOL = ODM3*ONENRM
319: *
320: DTPCRT = SQRT( ODM1 / BLKSIZ )
321: *
322: * Loop through eigenvalues of block nblk.
323: *
324: 60 CONTINUE
325: JBLK = 0
326: DO 170 J = J1, M
327: IF( IBLOCK( J ).NE.NBLK ) THEN
328: J1 = J
329: GO TO 180
330: END IF
331: JBLK = JBLK + 1
332: XJ = W( J )
333: *
334: * Skip all the work if the block size is one.
335: *
336: IF( BLKSIZ.EQ.1 ) THEN
337: WORK( INDRV1+1 ) = ONE
338: GO TO 140
339: END IF
340: *
341: * If eigenvalues j and j-1 are too close, add a relatively
342: * small perturbation.
343: *
344: IF( JBLK.GT.1 ) THEN
345: EPS1 = ABS( EPS*XJ )
346: PERTOL = TEN*EPS1
347: SEP = XJ - XJM
348: IF( SEP.LT.PERTOL )
349: $ XJ = XJM + PERTOL
350: END IF
351: *
352: ITS = 0
353: NRMCHK = 0
354: *
355: * Get random starting vector.
356: *
357: CALL DLARNV( 2, ISEED, BLKSIZ, WORK( INDRV1+1 ) )
358: *
359: * Copy the matrix T so it won't be destroyed in factorization.
360: *
361: CALL DCOPY( BLKSIZ, D( B1 ), 1, WORK( INDRV4+1 ), 1 )
362: CALL DCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV2+2 ), 1 )
363: CALL DCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV3+1 ), 1 )
364: *
365: * Compute LU factors with partial pivoting ( PT = LU )
366: *
367: TOL = ZERO
368: CALL DLAGTF( BLKSIZ, WORK( INDRV4+1 ), XJ, WORK( INDRV2+2 ),
369: $ WORK( INDRV3+1 ), TOL, WORK( INDRV5+1 ), IWORK,
370: $ IINFO )
371: *
372: * Update iteration count.
373: *
374: 70 CONTINUE
375: ITS = ITS + 1
376: IF( ITS.GT.MAXITS )
377: $ GO TO 120
378: *
379: * Normalize and scale the righthand side vector Pb.
380: *
381: JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 )
382: SCL = BLKSIZ*ONENRM*MAX( EPS,
383: $ ABS( WORK( INDRV4+BLKSIZ ) ) ) /
384: $ ABS( WORK( INDRV1+JMAX ) )
385: CALL DSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 )
386: *
387: * Solve the system LU = Pb.
388: *
389: CALL DLAGTS( -1, BLKSIZ, WORK( INDRV4+1 ), WORK( INDRV2+2 ),
390: $ WORK( INDRV3+1 ), WORK( INDRV5+1 ), IWORK,
391: $ WORK( INDRV1+1 ), TOL, IINFO )
392: *
393: * Reorthogonalize by modified Gram-Schmidt if eigenvalues are
394: * close enough.
395: *
396: IF( JBLK.EQ.1 )
397: $ GO TO 110
398: IF( ABS( XJ-XJM ).GT.ORTOL )
399: $ GPIND = J
400: IF( GPIND.NE.J ) THEN
401: DO 100 I = GPIND, J - 1
402: ZTR = ZERO
403: DO 80 JR = 1, BLKSIZ
404: ZTR = ZTR + WORK( INDRV1+JR )*
405: $ DBLE( Z( B1-1+JR, I ) )
406: 80 CONTINUE
407: DO 90 JR = 1, BLKSIZ
408: WORK( INDRV1+JR ) = WORK( INDRV1+JR ) -
409: $ ZTR*DBLE( Z( B1-1+JR, I ) )
410: 90 CONTINUE
411: 100 CONTINUE
412: END IF
413: *
414: * Check the infinity norm of the iterate.
415: *
416: 110 CONTINUE
417: JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 )
418: NRM = ABS( WORK( INDRV1+JMAX ) )
419: *
420: * Continue for additional iterations after norm reaches
421: * stopping criterion.
422: *
423: IF( NRM.LT.DTPCRT )
424: $ GO TO 70
425: NRMCHK = NRMCHK + 1
426: IF( NRMCHK.LT.EXTRA+1 )
427: $ GO TO 70
428: *
429: GO TO 130
430: *
431: * If stopping criterion was not satisfied, update info and
432: * store eigenvector number in array ifail.
433: *
434: 120 CONTINUE
435: INFO = INFO + 1
436: IFAIL( INFO ) = J
437: *
438: * Accept iterate as jth eigenvector.
439: *
440: 130 CONTINUE
441: SCL = ONE / DNRM2( BLKSIZ, WORK( INDRV1+1 ), 1 )
442: JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 )
443: IF( WORK( INDRV1+JMAX ).LT.ZERO )
444: $ SCL = -SCL
445: CALL DSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 )
446: 140 CONTINUE
447: DO 150 I = 1, N
448: Z( I, J ) = CZERO
449: 150 CONTINUE
450: DO 160 I = 1, BLKSIZ
451: Z( B1+I-1, J ) = DCMPLX( WORK( INDRV1+I ), ZERO )
452: 160 CONTINUE
453: *
454: * Save the shift to check eigenvalue spacing at next
455: * iteration.
456: *
457: XJM = XJ
458: *
459: 170 CONTINUE
460: 180 CONTINUE
461: *
462: RETURN
463: *
464: * End of ZSTEIN
465: *
466: 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