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: * =====================================================================
182: SUBROUTINE ZSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
183: $ IWORK, IFAIL, INFO )
184: *
185: * -- LAPACK computational routine (version 3.4.0) --
186: * -- LAPACK is a software package provided by Univ. of Tennessee, --
187: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
188: * November 2011
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