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