1: SUBROUTINE DSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
2: $ IWORK, IFAIL, INFO )
3: *
4: * -- LAPACK routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
8: *
9: * .. Scalar Arguments ..
10: INTEGER INFO, LDZ, M, N
11: * ..
12: * .. Array Arguments ..
13: INTEGER IBLOCK( * ), IFAIL( * ), ISPLIT( * ),
14: $ IWORK( * )
15: DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
16: * ..
17: *
18: * Purpose
19: * =======
20: *
21: * DSTEIN computes the eigenvectors of a real symmetric tridiagonal
22: * matrix T corresponding to specified eigenvalues, using inverse
23: * iteration.
24: *
25: * The maximum number of iterations allowed for each eigenvector is
26: * specified by an internal parameter MAXITS (currently set to 5).
27: *
28: * Arguments
29: * =========
30: *
31: * N (input) INTEGER
32: * The order of the matrix. N >= 0.
33: *
34: * D (input) DOUBLE PRECISION array, dimension (N)
35: * The n diagonal elements of the tridiagonal matrix T.
36: *
37: * E (input) DOUBLE PRECISION array, dimension (N-1)
38: * The (n-1) subdiagonal elements of the tridiagonal matrix
39: * T, in elements 1 to N-1.
40: *
41: * M (input) INTEGER
42: * The number of eigenvectors to be found. 0 <= M <= N.
43: *
44: * W (input) DOUBLE PRECISION array, dimension (N)
45: * The first M elements of W contain the eigenvalues for
46: * which eigenvectors are to be computed. The eigenvalues
47: * should be grouped by split-off block and ordered from
48: * smallest to largest within the block. ( The output array
49: * W from DSTEBZ with ORDER = 'B' is expected here. )
50: *
51: * IBLOCK (input) INTEGER array, dimension (N)
52: * The submatrix indices associated with the corresponding
53: * eigenvalues in W; IBLOCK(i)=1 if eigenvalue W(i) belongs to
54: * the first submatrix from the top, =2 if W(i) belongs to
55: * the second submatrix, etc. ( The output array IBLOCK
56: * from DSTEBZ is expected here. )
57: *
58: * ISPLIT (input) INTEGER array, dimension (N)
59: * The splitting points, at which T breaks up into submatrices.
60: * The first submatrix consists of rows/columns 1 to
61: * ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
62: * through ISPLIT( 2 ), etc.
63: * ( The output array ISPLIT from DSTEBZ is expected here. )
64: *
65: * Z (output) DOUBLE PRECISION array, dimension (LDZ, M)
66: * The computed eigenvectors. The eigenvector associated
67: * with the eigenvalue W(i) is stored in the i-th column of
68: * Z. Any vector which fails to converge is set to its current
69: * iterate after MAXITS iterations.
70: *
71: * LDZ (input) INTEGER
72: * The leading dimension of the array Z. LDZ >= max(1,N).
73: *
74: * WORK (workspace) DOUBLE PRECISION array, dimension (5*N)
75: *
76: * IWORK (workspace) INTEGER array, dimension (N)
77: *
78: * IFAIL (output) INTEGER array, dimension (M)
79: * On normal exit, all elements of IFAIL are zero.
80: * If one or more eigenvectors fail to converge after
81: * MAXITS iterations, then their indices are stored in
82: * array IFAIL.
83: *
84: * INFO (output) INTEGER
85: * = 0: successful exit.
86: * < 0: if INFO = -i, the i-th argument had an illegal value
87: * > 0: if INFO = i, then i eigenvectors failed to converge
88: * in MAXITS iterations. Their indices are stored in
89: * array IFAIL.
90: *
91: * Internal Parameters
92: * ===================
93: *
94: * MAXITS INTEGER, default = 5
95: * The maximum number of iterations performed.
96: *
97: * EXTRA INTEGER, default = 2
98: * The number of iterations performed after norm growth
99: * criterion is satisfied, should be at least 1.
100: *
101: * =====================================================================
102: *
103: * .. Parameters ..
104: DOUBLE PRECISION ZERO, ONE, TEN, ODM3, ODM1
105: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TEN = 1.0D+1,
106: $ ODM3 = 1.0D-3, ODM1 = 1.0D-1 )
107: INTEGER MAXITS, EXTRA
108: PARAMETER ( MAXITS = 5, EXTRA = 2 )
109: * ..
110: * .. Local Scalars ..
111: INTEGER B1, BLKSIZ, BN, GPIND, I, IINFO, INDRV1,
112: $ INDRV2, INDRV3, INDRV4, INDRV5, ITS, J, J1,
113: $ JBLK, JMAX, NBLK, NRMCHK
114: DOUBLE PRECISION DTPCRT, EPS, EPS1, NRM, ONENRM, ORTOL, PERTOL,
115: $ SCL, SEP, TOL, XJ, XJM, ZTR
116: * ..
117: * .. Local Arrays ..
118: INTEGER ISEED( 4 )
119: * ..
120: * .. External Functions ..
121: INTEGER IDAMAX
122: DOUBLE PRECISION DASUM, DDOT, DLAMCH, DNRM2
123: EXTERNAL IDAMAX, DASUM, DDOT, DLAMCH, DNRM2
124: * ..
125: * .. External Subroutines ..
126: EXTERNAL DAXPY, DCOPY, DLAGTF, DLAGTS, DLARNV, DSCAL,
127: $ XERBLA
128: * ..
129: * .. Intrinsic Functions ..
130: INTRINSIC ABS, MAX, SQRT
131: * ..
132: * .. Executable Statements ..
133: *
134: * Test the input parameters.
135: *
136: INFO = 0
137: DO 10 I = 1, M
138: IFAIL( I ) = 0
139: 10 CONTINUE
140: *
141: IF( N.LT.0 ) THEN
142: INFO = -1
143: ELSE IF( M.LT.0 .OR. M.GT.N ) THEN
144: INFO = -4
145: ELSE IF( LDZ.LT.MAX( 1, N ) ) THEN
146: INFO = -9
147: ELSE
148: DO 20 J = 2, M
149: IF( IBLOCK( J ).LT.IBLOCK( J-1 ) ) THEN
150: INFO = -6
151: GO TO 30
152: END IF
153: IF( IBLOCK( J ).EQ.IBLOCK( J-1 ) .AND. W( J ).LT.W( J-1 ) )
154: $ THEN
155: INFO = -5
156: GO TO 30
157: END IF
158: 20 CONTINUE
159: 30 CONTINUE
160: END IF
161: *
162: IF( INFO.NE.0 ) THEN
163: CALL XERBLA( 'DSTEIN', -INFO )
164: RETURN
165: END IF
166: *
167: * Quick return if possible
168: *
169: IF( N.EQ.0 .OR. M.EQ.0 ) THEN
170: RETURN
171: ELSE IF( N.EQ.1 ) THEN
172: Z( 1, 1 ) = ONE
173: RETURN
174: END IF
175: *
176: * Get machine constants.
177: *
178: EPS = DLAMCH( 'Precision' )
179: *
180: * Initialize seed for random number generator DLARNV.
181: *
182: DO 40 I = 1, 4
183: ISEED( I ) = 1
184: 40 CONTINUE
185: *
186: * Initialize pointers.
187: *
188: INDRV1 = 0
189: INDRV2 = INDRV1 + N
190: INDRV3 = INDRV2 + N
191: INDRV4 = INDRV3 + N
192: INDRV5 = INDRV4 + N
193: *
194: * Compute eigenvectors of matrix blocks.
195: *
196: J1 = 1
197: DO 160 NBLK = 1, IBLOCK( M )
198: *
199: * Find starting and ending indices of block nblk.
200: *
201: IF( NBLK.EQ.1 ) THEN
202: B1 = 1
203: ELSE
204: B1 = ISPLIT( NBLK-1 ) + 1
205: END IF
206: BN = ISPLIT( NBLK )
207: BLKSIZ = BN - B1 + 1
208: IF( BLKSIZ.EQ.1 )
209: $ GO TO 60
210: GPIND = B1
211: *
212: * Compute reorthogonalization criterion and stopping criterion.
213: *
214: ONENRM = ABS( D( B1 ) ) + ABS( E( B1 ) )
215: ONENRM = MAX( ONENRM, ABS( D( BN ) )+ABS( E( BN-1 ) ) )
216: DO 50 I = B1 + 1, BN - 1
217: ONENRM = MAX( ONENRM, ABS( D( I ) )+ABS( E( I-1 ) )+
218: $ ABS( E( I ) ) )
219: 50 CONTINUE
220: ORTOL = ODM3*ONENRM
221: *
222: DTPCRT = SQRT( ODM1 / BLKSIZ )
223: *
224: * Loop through eigenvalues of block nblk.
225: *
226: 60 CONTINUE
227: JBLK = 0
228: DO 150 J = J1, M
229: IF( IBLOCK( J ).NE.NBLK ) THEN
230: J1 = J
231: GO TO 160
232: END IF
233: JBLK = JBLK + 1
234: XJ = W( J )
235: *
236: * Skip all the work if the block size is one.
237: *
238: IF( BLKSIZ.EQ.1 ) THEN
239: WORK( INDRV1+1 ) = ONE
240: GO TO 120
241: END IF
242: *
243: * If eigenvalues j and j-1 are too close, add a relatively
244: * small perturbation.
245: *
246: IF( JBLK.GT.1 ) THEN
247: EPS1 = ABS( EPS*XJ )
248: PERTOL = TEN*EPS1
249: SEP = XJ - XJM
250: IF( SEP.LT.PERTOL )
251: $ XJ = XJM + PERTOL
252: END IF
253: *
254: ITS = 0
255: NRMCHK = 0
256: *
257: * Get random starting vector.
258: *
259: CALL DLARNV( 2, ISEED, BLKSIZ, WORK( INDRV1+1 ) )
260: *
261: * Copy the matrix T so it won't be destroyed in factorization.
262: *
263: CALL DCOPY( BLKSIZ, D( B1 ), 1, WORK( INDRV4+1 ), 1 )
264: CALL DCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV2+2 ), 1 )
265: CALL DCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV3+1 ), 1 )
266: *
267: * Compute LU factors with partial pivoting ( PT = LU )
268: *
269: TOL = ZERO
270: CALL DLAGTF( BLKSIZ, WORK( INDRV4+1 ), XJ, WORK( INDRV2+2 ),
271: $ WORK( INDRV3+1 ), TOL, WORK( INDRV5+1 ), IWORK,
272: $ IINFO )
273: *
274: * Update iteration count.
275: *
276: 70 CONTINUE
277: ITS = ITS + 1
278: IF( ITS.GT.MAXITS )
279: $ GO TO 100
280: *
281: * Normalize and scale the righthand side vector Pb.
282: *
283: SCL = BLKSIZ*ONENRM*MAX( EPS,
284: $ ABS( WORK( INDRV4+BLKSIZ ) ) ) /
285: $ DASUM( BLKSIZ, WORK( INDRV1+1 ), 1 )
286: CALL DSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 )
287: *
288: * Solve the system LU = Pb.
289: *
290: CALL DLAGTS( -1, BLKSIZ, WORK( INDRV4+1 ), WORK( INDRV2+2 ),
291: $ WORK( INDRV3+1 ), WORK( INDRV5+1 ), IWORK,
292: $ WORK( INDRV1+1 ), TOL, IINFO )
293: *
294: * Reorthogonalize by modified Gram-Schmidt if eigenvalues are
295: * close enough.
296: *
297: IF( JBLK.EQ.1 )
298: $ GO TO 90
299: IF( ABS( XJ-XJM ).GT.ORTOL )
300: $ GPIND = J
301: IF( GPIND.NE.J ) THEN
302: DO 80 I = GPIND, J - 1
303: ZTR = -DDOT( BLKSIZ, WORK( INDRV1+1 ), 1, Z( B1, I ),
304: $ 1 )
305: CALL DAXPY( BLKSIZ, ZTR, Z( B1, I ), 1,
306: $ WORK( INDRV1+1 ), 1 )
307: 80 CONTINUE
308: END IF
309: *
310: * Check the infinity norm of the iterate.
311: *
312: 90 CONTINUE
313: JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 )
314: NRM = ABS( WORK( INDRV1+JMAX ) )
315: *
316: * Continue for additional iterations after norm reaches
317: * stopping criterion.
318: *
319: IF( NRM.LT.DTPCRT )
320: $ GO TO 70
321: NRMCHK = NRMCHK + 1
322: IF( NRMCHK.LT.EXTRA+1 )
323: $ GO TO 70
324: *
325: GO TO 110
326: *
327: * If stopping criterion was not satisfied, update info and
328: * store eigenvector number in array ifail.
329: *
330: 100 CONTINUE
331: INFO = INFO + 1
332: IFAIL( INFO ) = J
333: *
334: * Accept iterate as jth eigenvector.
335: *
336: 110 CONTINUE
337: SCL = ONE / DNRM2( BLKSIZ, WORK( INDRV1+1 ), 1 )
338: JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 )
339: IF( WORK( INDRV1+JMAX ).LT.ZERO )
340: $ SCL = -SCL
341: CALL DSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 )
342: 120 CONTINUE
343: DO 130 I = 1, N
344: Z( I, J ) = ZERO
345: 130 CONTINUE
346: DO 140 I = 1, BLKSIZ
347: Z( B1+I-1, J ) = WORK( INDRV1+I )
348: 140 CONTINUE
349: *
350: * Save the shift to check eigenvalue spacing at next
351: * iteration.
352: *
353: XJM = XJ
354: *
355: 150 CONTINUE
356: 160 CONTINUE
357: *
358: RETURN
359: *
360: * End of DSTEIN
361: *
362: 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