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