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
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