1 SUBROUTINE DSTEVR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, ABSTOL,
2 $ M, W, Z, LDZ, ISUPPZ, WORK, LWORK, IWORK,
3 $ LIWORK, INFO )
4 *
5 * -- LAPACK driver routine (version 3.2) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * November 2006
9 *
10 * .. Scalar Arguments ..
11 CHARACTER JOBZ, RANGE
12 INTEGER IL, INFO, IU, LDZ, LIWORK, LWORK, M, N
13 DOUBLE PRECISION ABSTOL, VL, VU
14 * ..
15 * .. Array Arguments ..
16 INTEGER ISUPPZ( * ), IWORK( * )
17 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * DSTEVR computes selected eigenvalues and, optionally, eigenvectors
24 * of a real symmetric tridiagonal matrix T. Eigenvalues and
25 * eigenvectors can be selected by specifying either a range of values
26 * or a range of indices for the desired eigenvalues.
27 *
28 * Whenever possible, DSTEVR calls DSTEMR to compute the
29 * eigenspectrum using Relatively Robust Representations. DSTEMR
30 * computes eigenvalues by the dqds algorithm, while orthogonal
31 * eigenvectors are computed from various "good" L D L^T representations
32 * (also known as Relatively Robust Representations). Gram-Schmidt
33 * orthogonalization is avoided as far as possible. More specifically,
34 * the various steps of the algorithm are as follows. For the i-th
35 * unreduced block of T,
36 * (a) Compute T - sigma_i = L_i D_i L_i^T, such that L_i D_i L_i^T
37 * is a relatively robust representation,
38 * (b) Compute the eigenvalues, lambda_j, of L_i D_i L_i^T to high
39 * relative accuracy by the dqds algorithm,
40 * (c) If there is a cluster of close eigenvalues, "choose" sigma_i
41 * close to the cluster, and go to step (a),
42 * (d) Given the approximate eigenvalue lambda_j of L_i D_i L_i^T,
43 * compute the corresponding eigenvector by forming a
44 * rank-revealing twisted factorization.
45 * The desired accuracy of the output can be specified by the input
46 * parameter ABSTOL.
47 *
48 * For more details, see "A new O(n^2) algorithm for the symmetric
49 * tridiagonal eigenvalue/eigenvector problem", by Inderjit Dhillon,
50 * Computer Science Division Technical Report No. UCB//CSD-97-971,
51 * UC Berkeley, May 1997.
52 *
53 *
54 * Note 1 : DSTEVR calls DSTEMR when the full spectrum is requested
55 * on machines which conform to the ieee-754 floating point standard.
56 * DSTEVR calls DSTEBZ and DSTEIN on non-ieee machines and
57 * when partial spectrum requests are made.
58 *
59 * Normal execution of DSTEMR may create NaNs and infinities and
60 * hence may abort due to a floating point exception in environments
61 * which do not handle NaNs and infinities in the ieee standard default
62 * manner.
63 *
64 * Arguments
65 * =========
66 *
67 * JOBZ (input) CHARACTER*1
68 * = 'N': Compute eigenvalues only;
69 * = 'V': Compute eigenvalues and eigenvectors.
70 *
71 * RANGE (input) CHARACTER*1
72 * = 'A': all eigenvalues will be found.
73 * = 'V': all eigenvalues in the half-open interval (VL,VU]
74 * will be found.
75 * = 'I': the IL-th through IU-th eigenvalues will be found.
76 ********** For RANGE = 'V' or 'I' and IU - IL < N - 1, DSTEBZ and
77 ********** DSTEIN are called
78 *
79 * N (input) INTEGER
80 * The order of the matrix. N >= 0.
81 *
82 * D (input/output) DOUBLE PRECISION array, dimension (N)
83 * On entry, the n diagonal elements of the tridiagonal matrix
84 * A.
85 * On exit, D may be multiplied by a constant factor chosen
86 * to avoid over/underflow in computing the eigenvalues.
87 *
88 * E (input/output) DOUBLE PRECISION array, dimension (max(1,N-1))
89 * On entry, the (n-1) subdiagonal elements of the tridiagonal
90 * matrix A in elements 1 to N-1 of E.
91 * On exit, E may be multiplied by a constant factor chosen
92 * to avoid over/underflow in computing the eigenvalues.
93 *
94 * VL (input) DOUBLE PRECISION
95 * VU (input) DOUBLE PRECISION
96 * If RANGE='V', the lower and upper bounds of the interval to
97 * be searched for eigenvalues. VL < VU.
98 * Not referenced if RANGE = 'A' or 'I'.
99 *
100 * IL (input) INTEGER
101 * IU (input) INTEGER
102 * If RANGE='I', the indices (in ascending order) of the
103 * smallest and largest eigenvalues to be returned.
104 * 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
105 * Not referenced if RANGE = 'A' or 'V'.
106 *
107 * ABSTOL (input) DOUBLE PRECISION
108 * The absolute error tolerance for the eigenvalues.
109 * An approximate eigenvalue is accepted as converged
110 * when it is determined to lie in an interval [a,b]
111 * of width less than or equal to
112 *
113 * ABSTOL + EPS * max( |a|,|b| ) ,
114 *
115 * where EPS is the machine precision. If ABSTOL is less than
116 * or equal to zero, then EPS*|T| will be used in its place,
117 * where |T| is the 1-norm of the tridiagonal matrix obtained
118 * by reducing A to tridiagonal form.
119 *
120 * See "Computing Small Singular Values of Bidiagonal Matrices
121 * with Guaranteed High Relative Accuracy," by Demmel and
122 * Kahan, LAPACK Working Note #3.
123 *
124 * If high relative accuracy is important, set ABSTOL to
125 * DLAMCH( 'Safe minimum' ). Doing so will guarantee that
126 * eigenvalues are computed to high relative accuracy when
127 * possible in future releases. The current code does not
128 * make any guarantees about high relative accuracy, but
129 * future releases will. See J. Barlow and J. Demmel,
130 * "Computing Accurate Eigensystems of Scaled Diagonally
131 * Dominant Matrices", LAPACK Working Note #7, for a discussion
132 * of which matrices define their eigenvalues to high relative
133 * accuracy.
134 *
135 * M (output) INTEGER
136 * The total number of eigenvalues found. 0 <= M <= N.
137 * If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
138 *
139 * W (output) DOUBLE PRECISION array, dimension (N)
140 * The first M elements contain the selected eigenvalues in
141 * ascending order.
142 *
143 * Z (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M) )
144 * If JOBZ = 'V', then if INFO = 0, the first M columns of Z
145 * contain the orthonormal eigenvectors of the matrix A
146 * corresponding to the selected eigenvalues, with the i-th
147 * column of Z holding the eigenvector associated with W(i).
148 * Note: the user must ensure that at least max(1,M) columns are
149 * supplied in the array Z; if RANGE = 'V', the exact value of M
150 * is not known in advance and an upper bound must be used.
151 *
152 * LDZ (input) INTEGER
153 * The leading dimension of the array Z. LDZ >= 1, and if
154 * JOBZ = 'V', LDZ >= max(1,N).
155 *
156 * ISUPPZ (output) INTEGER array, dimension ( 2*max(1,M) )
157 * The support of the eigenvectors in Z, i.e., the indices
158 * indicating the nonzero elements in Z. The i-th eigenvector
159 * is nonzero only in elements ISUPPZ( 2*i-1 ) through
160 * ISUPPZ( 2*i ).
161 ********** Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
162 *
163 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
164 * On exit, if INFO = 0, WORK(1) returns the optimal (and
165 * minimal) LWORK.
166 *
167 * LWORK (input) INTEGER
168 * The dimension of the array WORK. LWORK >= max(1,20*N).
169 *
170 * If LWORK = -1, then a workspace query is assumed; the routine
171 * only calculates the optimal sizes of the WORK and IWORK
172 * arrays, returns these values as the first entries of the WORK
173 * and IWORK arrays, and no error message related to LWORK or
174 * LIWORK is issued by XERBLA.
175 *
176 * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
177 * On exit, if INFO = 0, IWORK(1) returns the optimal (and
178 * minimal) LIWORK.
179 *
180 * LIWORK (input) INTEGER
181 * The dimension of the array IWORK. LIWORK >= max(1,10*N).
182 *
183 * If LIWORK = -1, then a workspace query is assumed; the
184 * routine only calculates the optimal sizes of the WORK and
185 * IWORK arrays, returns these values as the first entries of
186 * the WORK and IWORK arrays, and no error message related to
187 * LWORK or LIWORK is issued by XERBLA.
188 *
189 * INFO (output) INTEGER
190 * = 0: successful exit
191 * < 0: if INFO = -i, the i-th argument had an illegal value
192 * > 0: Internal error
193 *
194 * Further Details
195 * ===============
196 *
197 * Based on contributions by
198 * Inderjit Dhillon, IBM Almaden, USA
199 * Osni Marques, LBNL/NERSC, USA
200 * Ken Stanley, Computer Science Division, University of
201 * California at Berkeley, USA
202 *
203 * =====================================================================
204 *
205 * .. Parameters ..
206 DOUBLE PRECISION ZERO, ONE, TWO
207 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
208 * ..
209 * .. Local Scalars ..
210 LOGICAL ALLEIG, INDEIG, TEST, LQUERY, VALEIG, WANTZ,
211 $ TRYRAC
212 CHARACTER ORDER
213 INTEGER I, IEEEOK, IMAX, INDIBL, INDIFL, INDISP,
214 $ INDIWO, ISCALE, ITMP1, J, JJ, LIWMIN, LWMIN,
215 $ NSPLIT
216 DOUBLE PRECISION BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM,
217 $ TMP1, TNRM, VLL, VUU
218 * ..
219 * .. External Functions ..
220 LOGICAL LSAME
221 INTEGER ILAENV
222 DOUBLE PRECISION DLAMCH, DLANST
223 EXTERNAL LSAME, ILAENV, DLAMCH, DLANST
224 * ..
225 * .. External Subroutines ..
226 EXTERNAL DCOPY, DSCAL, DSTEBZ, DSTEMR, DSTEIN, DSTERF,
227 $ DSWAP, XERBLA
228 * ..
229 * .. Intrinsic Functions ..
230 INTRINSIC MAX, MIN, SQRT
231 * ..
232 * .. Executable Statements ..
233 *
234 *
235 * Test the input parameters.
236 *
237 IEEEOK = ILAENV( 10, 'DSTEVR', 'N', 1, 2, 3, 4 )
238 *
239 WANTZ = LSAME( JOBZ, 'V' )
240 ALLEIG = LSAME( RANGE, 'A' )
241 VALEIG = LSAME( RANGE, 'V' )
242 INDEIG = LSAME( RANGE, 'I' )
243 *
244 LQUERY = ( ( LWORK.EQ.-1 ) .OR. ( LIWORK.EQ.-1 ) )
245 LWMIN = MAX( 1, 20*N )
246 LIWMIN = MAX( 1, 10*N )
247 *
248 *
249 INFO = 0
250 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
251 INFO = -1
252 ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN
253 INFO = -2
254 ELSE IF( N.LT.0 ) THEN
255 INFO = -3
256 ELSE
257 IF( VALEIG ) THEN
258 IF( N.GT.0 .AND. VU.LE.VL )
259 $ INFO = -7
260 ELSE IF( INDEIG ) THEN
261 IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN
262 INFO = -8
263 ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN
264 INFO = -9
265 END IF
266 END IF
267 END IF
268 IF( INFO.EQ.0 ) THEN
269 IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
270 INFO = -14
271 END IF
272 END IF
273 *
274 IF( INFO.EQ.0 ) THEN
275 WORK( 1 ) = LWMIN
276 IWORK( 1 ) = LIWMIN
277 *
278 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
279 INFO = -17
280 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
281 INFO = -19
282 END IF
283 END IF
284 *
285 IF( INFO.NE.0 ) THEN
286 CALL XERBLA( 'DSTEVR', -INFO )
287 RETURN
288 ELSE IF( LQUERY ) THEN
289 RETURN
290 END IF
291 *
292 * Quick return if possible
293 *
294 M = 0
295 IF( N.EQ.0 )
296 $ RETURN
297 *
298 IF( N.EQ.1 ) THEN
299 IF( ALLEIG .OR. INDEIG ) THEN
300 M = 1
301 W( 1 ) = D( 1 )
302 ELSE
303 IF( VL.LT.D( 1 ) .AND. VU.GE.D( 1 ) ) THEN
304 M = 1
305 W( 1 ) = D( 1 )
306 END IF
307 END IF
308 IF( WANTZ )
309 $ Z( 1, 1 ) = ONE
310 RETURN
311 END IF
312 *
313 * Get machine constants.
314 *
315 SAFMIN = DLAMCH( 'Safe minimum' )
316 EPS = DLAMCH( 'Precision' )
317 SMLNUM = SAFMIN / EPS
318 BIGNUM = ONE / SMLNUM
319 RMIN = SQRT( SMLNUM )
320 RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) )
321 *
322 *
323 * Scale matrix to allowable range, if necessary.
324 *
325 ISCALE = 0
326 VLL = VL
327 VUU = VU
328 *
329 TNRM = DLANST( 'M', N, D, E )
330 IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN
331 ISCALE = 1
332 SIGMA = RMIN / TNRM
333 ELSE IF( TNRM.GT.RMAX ) THEN
334 ISCALE = 1
335 SIGMA = RMAX / TNRM
336 END IF
337 IF( ISCALE.EQ.1 ) THEN
338 CALL DSCAL( N, SIGMA, D, 1 )
339 CALL DSCAL( N-1, SIGMA, E( 1 ), 1 )
340 IF( VALEIG ) THEN
341 VLL = VL*SIGMA
342 VUU = VU*SIGMA
343 END IF
344 END IF
345
346 * Initialize indices into workspaces. Note: These indices are used only
347 * if DSTERF or DSTEMR fail.
348
349 * IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in DSTEBZ and
350 * stores the block indices of each of the M<=N eigenvalues.
351 INDIBL = 1
352 * IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in DSTEBZ and
353 * stores the starting and finishing indices of each block.
354 INDISP = INDIBL + N
355 * IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors
356 * that corresponding to eigenvectors that fail to converge in
357 * DSTEIN. This information is discarded; if any fail, the driver
358 * returns INFO > 0.
359 INDIFL = INDISP + N
360 * INDIWO is the offset of the remaining integer workspace.
361 INDIWO = INDISP + N
362 *
363 * If all eigenvalues are desired, then
364 * call DSTERF or DSTEMR. If this fails for some eigenvalue, then
365 * try DSTEBZ.
366 *
367 *
368 TEST = .FALSE.
369 IF( INDEIG ) THEN
370 IF( IL.EQ.1 .AND. IU.EQ.N ) THEN
371 TEST = .TRUE.
372 END IF
373 END IF
374 IF( ( ALLEIG .OR. TEST ) .AND. IEEEOK.EQ.1 ) THEN
375 CALL DCOPY( N-1, E( 1 ), 1, WORK( 1 ), 1 )
376 IF( .NOT.WANTZ ) THEN
377 CALL DCOPY( N, D, 1, W, 1 )
378 CALL DSTERF( N, W, WORK, INFO )
379 ELSE
380 CALL DCOPY( N, D, 1, WORK( N+1 ), 1 )
381 IF (ABSTOL .LE. TWO*N*EPS) THEN
382 TRYRAC = .TRUE.
383 ELSE
384 TRYRAC = .FALSE.
385 END IF
386 CALL DSTEMR( JOBZ, 'A', N, WORK( N+1 ), WORK, VL, VU, IL,
387 $ IU, M, W, Z, LDZ, N, ISUPPZ, TRYRAC,
388 $ WORK( 2*N+1 ), LWORK-2*N, IWORK, LIWORK, INFO )
389 *
390 END IF
391 IF( INFO.EQ.0 ) THEN
392 M = N
393 GO TO 10
394 END IF
395 INFO = 0
396 END IF
397 *
398 * Otherwise, call DSTEBZ and, if eigenvectors are desired, DSTEIN.
399 *
400 IF( WANTZ ) THEN
401 ORDER = 'B'
402 ELSE
403 ORDER = 'E'
404 END IF
405
406 CALL DSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTOL, D, E, M,
407 $ NSPLIT, W, IWORK( INDIBL ), IWORK( INDISP ), WORK,
408 $ IWORK( INDIWO ), INFO )
409 *
410 IF( WANTZ ) THEN
411 CALL DSTEIN( N, D, E, M, W, IWORK( INDIBL ), IWORK( INDISP ),
412 $ Z, LDZ, WORK, IWORK( INDIWO ), IWORK( INDIFL ),
413 $ INFO )
414 END IF
415 *
416 * If matrix was scaled, then rescale eigenvalues appropriately.
417 *
418 10 CONTINUE
419 IF( ISCALE.EQ.1 ) THEN
420 IF( INFO.EQ.0 ) THEN
421 IMAX = M
422 ELSE
423 IMAX = INFO - 1
424 END IF
425 CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
426 END IF
427 *
428 * If eigenvalues are not in order, then sort them, along with
429 * eigenvectors.
430 *
431 IF( WANTZ ) THEN
432 DO 30 J = 1, M - 1
433 I = 0
434 TMP1 = W( J )
435 DO 20 JJ = J + 1, M
436 IF( W( JJ ).LT.TMP1 ) THEN
437 I = JJ
438 TMP1 = W( JJ )
439 END IF
440 20 CONTINUE
441 *
442 IF( I.NE.0 ) THEN
443 ITMP1 = IWORK( I )
444 W( I ) = W( J )
445 IWORK( I ) = IWORK( J )
446 W( J ) = TMP1
447 IWORK( J ) = ITMP1
448 CALL DSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 )
449 END IF
450 30 CONTINUE
451 END IF
452 *
453 * Causes problems with tests 19 & 20:
454 * IF (wantz .and. INDEIG ) Z( 1,1) = Z(1,1) / 1.002 + .002
455 *
456 *
457 WORK( 1 ) = LWMIN
458 IWORK( 1 ) = LIWMIN
459 RETURN
460 *
461 * End of DSTEVR
462 *
463 END
2 $ M, W, Z, LDZ, ISUPPZ, WORK, LWORK, IWORK,
3 $ LIWORK, INFO )
4 *
5 * -- LAPACK driver routine (version 3.2) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * November 2006
9 *
10 * .. Scalar Arguments ..
11 CHARACTER JOBZ, RANGE
12 INTEGER IL, INFO, IU, LDZ, LIWORK, LWORK, M, N
13 DOUBLE PRECISION ABSTOL, VL, VU
14 * ..
15 * .. Array Arguments ..
16 INTEGER ISUPPZ( * ), IWORK( * )
17 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * DSTEVR computes selected eigenvalues and, optionally, eigenvectors
24 * of a real symmetric tridiagonal matrix T. Eigenvalues and
25 * eigenvectors can be selected by specifying either a range of values
26 * or a range of indices for the desired eigenvalues.
27 *
28 * Whenever possible, DSTEVR calls DSTEMR to compute the
29 * eigenspectrum using Relatively Robust Representations. DSTEMR
30 * computes eigenvalues by the dqds algorithm, while orthogonal
31 * eigenvectors are computed from various "good" L D L^T representations
32 * (also known as Relatively Robust Representations). Gram-Schmidt
33 * orthogonalization is avoided as far as possible. More specifically,
34 * the various steps of the algorithm are as follows. For the i-th
35 * unreduced block of T,
36 * (a) Compute T - sigma_i = L_i D_i L_i^T, such that L_i D_i L_i^T
37 * is a relatively robust representation,
38 * (b) Compute the eigenvalues, lambda_j, of L_i D_i L_i^T to high
39 * relative accuracy by the dqds algorithm,
40 * (c) If there is a cluster of close eigenvalues, "choose" sigma_i
41 * close to the cluster, and go to step (a),
42 * (d) Given the approximate eigenvalue lambda_j of L_i D_i L_i^T,
43 * compute the corresponding eigenvector by forming a
44 * rank-revealing twisted factorization.
45 * The desired accuracy of the output can be specified by the input
46 * parameter ABSTOL.
47 *
48 * For more details, see "A new O(n^2) algorithm for the symmetric
49 * tridiagonal eigenvalue/eigenvector problem", by Inderjit Dhillon,
50 * Computer Science Division Technical Report No. UCB//CSD-97-971,
51 * UC Berkeley, May 1997.
52 *
53 *
54 * Note 1 : DSTEVR calls DSTEMR when the full spectrum is requested
55 * on machines which conform to the ieee-754 floating point standard.
56 * DSTEVR calls DSTEBZ and DSTEIN on non-ieee machines and
57 * when partial spectrum requests are made.
58 *
59 * Normal execution of DSTEMR may create NaNs and infinities and
60 * hence may abort due to a floating point exception in environments
61 * which do not handle NaNs and infinities in the ieee standard default
62 * manner.
63 *
64 * Arguments
65 * =========
66 *
67 * JOBZ (input) CHARACTER*1
68 * = 'N': Compute eigenvalues only;
69 * = 'V': Compute eigenvalues and eigenvectors.
70 *
71 * RANGE (input) CHARACTER*1
72 * = 'A': all eigenvalues will be found.
73 * = 'V': all eigenvalues in the half-open interval (VL,VU]
74 * will be found.
75 * = 'I': the IL-th through IU-th eigenvalues will be found.
76 ********** For RANGE = 'V' or 'I' and IU - IL < N - 1, DSTEBZ and
77 ********** DSTEIN are called
78 *
79 * N (input) INTEGER
80 * The order of the matrix. N >= 0.
81 *
82 * D (input/output) DOUBLE PRECISION array, dimension (N)
83 * On entry, the n diagonal elements of the tridiagonal matrix
84 * A.
85 * On exit, D may be multiplied by a constant factor chosen
86 * to avoid over/underflow in computing the eigenvalues.
87 *
88 * E (input/output) DOUBLE PRECISION array, dimension (max(1,N-1))
89 * On entry, the (n-1) subdiagonal elements of the tridiagonal
90 * matrix A in elements 1 to N-1 of E.
91 * On exit, E may be multiplied by a constant factor chosen
92 * to avoid over/underflow in computing the eigenvalues.
93 *
94 * VL (input) DOUBLE PRECISION
95 * VU (input) DOUBLE PRECISION
96 * If RANGE='V', the lower and upper bounds of the interval to
97 * be searched for eigenvalues. VL < VU.
98 * Not referenced if RANGE = 'A' or 'I'.
99 *
100 * IL (input) INTEGER
101 * IU (input) INTEGER
102 * If RANGE='I', the indices (in ascending order) of the
103 * smallest and largest eigenvalues to be returned.
104 * 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
105 * Not referenced if RANGE = 'A' or 'V'.
106 *
107 * ABSTOL (input) DOUBLE PRECISION
108 * The absolute error tolerance for the eigenvalues.
109 * An approximate eigenvalue is accepted as converged
110 * when it is determined to lie in an interval [a,b]
111 * of width less than or equal to
112 *
113 * ABSTOL + EPS * max( |a|,|b| ) ,
114 *
115 * where EPS is the machine precision. If ABSTOL is less than
116 * or equal to zero, then EPS*|T| will be used in its place,
117 * where |T| is the 1-norm of the tridiagonal matrix obtained
118 * by reducing A to tridiagonal form.
119 *
120 * See "Computing Small Singular Values of Bidiagonal Matrices
121 * with Guaranteed High Relative Accuracy," by Demmel and
122 * Kahan, LAPACK Working Note #3.
123 *
124 * If high relative accuracy is important, set ABSTOL to
125 * DLAMCH( 'Safe minimum' ). Doing so will guarantee that
126 * eigenvalues are computed to high relative accuracy when
127 * possible in future releases. The current code does not
128 * make any guarantees about high relative accuracy, but
129 * future releases will. See J. Barlow and J. Demmel,
130 * "Computing Accurate Eigensystems of Scaled Diagonally
131 * Dominant Matrices", LAPACK Working Note #7, for a discussion
132 * of which matrices define their eigenvalues to high relative
133 * accuracy.
134 *
135 * M (output) INTEGER
136 * The total number of eigenvalues found. 0 <= M <= N.
137 * If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
138 *
139 * W (output) DOUBLE PRECISION array, dimension (N)
140 * The first M elements contain the selected eigenvalues in
141 * ascending order.
142 *
143 * Z (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M) )
144 * If JOBZ = 'V', then if INFO = 0, the first M columns of Z
145 * contain the orthonormal eigenvectors of the matrix A
146 * corresponding to the selected eigenvalues, with the i-th
147 * column of Z holding the eigenvector associated with W(i).
148 * Note: the user must ensure that at least max(1,M) columns are
149 * supplied in the array Z; if RANGE = 'V', the exact value of M
150 * is not known in advance and an upper bound must be used.
151 *
152 * LDZ (input) INTEGER
153 * The leading dimension of the array Z. LDZ >= 1, and if
154 * JOBZ = 'V', LDZ >= max(1,N).
155 *
156 * ISUPPZ (output) INTEGER array, dimension ( 2*max(1,M) )
157 * The support of the eigenvectors in Z, i.e., the indices
158 * indicating the nonzero elements in Z. The i-th eigenvector
159 * is nonzero only in elements ISUPPZ( 2*i-1 ) through
160 * ISUPPZ( 2*i ).
161 ********** Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
162 *
163 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
164 * On exit, if INFO = 0, WORK(1) returns the optimal (and
165 * minimal) LWORK.
166 *
167 * LWORK (input) INTEGER
168 * The dimension of the array WORK. LWORK >= max(1,20*N).
169 *
170 * If LWORK = -1, then a workspace query is assumed; the routine
171 * only calculates the optimal sizes of the WORK and IWORK
172 * arrays, returns these values as the first entries of the WORK
173 * and IWORK arrays, and no error message related to LWORK or
174 * LIWORK is issued by XERBLA.
175 *
176 * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
177 * On exit, if INFO = 0, IWORK(1) returns the optimal (and
178 * minimal) LIWORK.
179 *
180 * LIWORK (input) INTEGER
181 * The dimension of the array IWORK. LIWORK >= max(1,10*N).
182 *
183 * If LIWORK = -1, then a workspace query is assumed; the
184 * routine only calculates the optimal sizes of the WORK and
185 * IWORK arrays, returns these values as the first entries of
186 * the WORK and IWORK arrays, and no error message related to
187 * LWORK or LIWORK is issued by XERBLA.
188 *
189 * INFO (output) INTEGER
190 * = 0: successful exit
191 * < 0: if INFO = -i, the i-th argument had an illegal value
192 * > 0: Internal error
193 *
194 * Further Details
195 * ===============
196 *
197 * Based on contributions by
198 * Inderjit Dhillon, IBM Almaden, USA
199 * Osni Marques, LBNL/NERSC, USA
200 * Ken Stanley, Computer Science Division, University of
201 * California at Berkeley, USA
202 *
203 * =====================================================================
204 *
205 * .. Parameters ..
206 DOUBLE PRECISION ZERO, ONE, TWO
207 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
208 * ..
209 * .. Local Scalars ..
210 LOGICAL ALLEIG, INDEIG, TEST, LQUERY, VALEIG, WANTZ,
211 $ TRYRAC
212 CHARACTER ORDER
213 INTEGER I, IEEEOK, IMAX, INDIBL, INDIFL, INDISP,
214 $ INDIWO, ISCALE, ITMP1, J, JJ, LIWMIN, LWMIN,
215 $ NSPLIT
216 DOUBLE PRECISION BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM,
217 $ TMP1, TNRM, VLL, VUU
218 * ..
219 * .. External Functions ..
220 LOGICAL LSAME
221 INTEGER ILAENV
222 DOUBLE PRECISION DLAMCH, DLANST
223 EXTERNAL LSAME, ILAENV, DLAMCH, DLANST
224 * ..
225 * .. External Subroutines ..
226 EXTERNAL DCOPY, DSCAL, DSTEBZ, DSTEMR, DSTEIN, DSTERF,
227 $ DSWAP, XERBLA
228 * ..
229 * .. Intrinsic Functions ..
230 INTRINSIC MAX, MIN, SQRT
231 * ..
232 * .. Executable Statements ..
233 *
234 *
235 * Test the input parameters.
236 *
237 IEEEOK = ILAENV( 10, 'DSTEVR', 'N', 1, 2, 3, 4 )
238 *
239 WANTZ = LSAME( JOBZ, 'V' )
240 ALLEIG = LSAME( RANGE, 'A' )
241 VALEIG = LSAME( RANGE, 'V' )
242 INDEIG = LSAME( RANGE, 'I' )
243 *
244 LQUERY = ( ( LWORK.EQ.-1 ) .OR. ( LIWORK.EQ.-1 ) )
245 LWMIN = MAX( 1, 20*N )
246 LIWMIN = MAX( 1, 10*N )
247 *
248 *
249 INFO = 0
250 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
251 INFO = -1
252 ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN
253 INFO = -2
254 ELSE IF( N.LT.0 ) THEN
255 INFO = -3
256 ELSE
257 IF( VALEIG ) THEN
258 IF( N.GT.0 .AND. VU.LE.VL )
259 $ INFO = -7
260 ELSE IF( INDEIG ) THEN
261 IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN
262 INFO = -8
263 ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN
264 INFO = -9
265 END IF
266 END IF
267 END IF
268 IF( INFO.EQ.0 ) THEN
269 IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
270 INFO = -14
271 END IF
272 END IF
273 *
274 IF( INFO.EQ.0 ) THEN
275 WORK( 1 ) = LWMIN
276 IWORK( 1 ) = LIWMIN
277 *
278 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
279 INFO = -17
280 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
281 INFO = -19
282 END IF
283 END IF
284 *
285 IF( INFO.NE.0 ) THEN
286 CALL XERBLA( 'DSTEVR', -INFO )
287 RETURN
288 ELSE IF( LQUERY ) THEN
289 RETURN
290 END IF
291 *
292 * Quick return if possible
293 *
294 M = 0
295 IF( N.EQ.0 )
296 $ RETURN
297 *
298 IF( N.EQ.1 ) THEN
299 IF( ALLEIG .OR. INDEIG ) THEN
300 M = 1
301 W( 1 ) = D( 1 )
302 ELSE
303 IF( VL.LT.D( 1 ) .AND. VU.GE.D( 1 ) ) THEN
304 M = 1
305 W( 1 ) = D( 1 )
306 END IF
307 END IF
308 IF( WANTZ )
309 $ Z( 1, 1 ) = ONE
310 RETURN
311 END IF
312 *
313 * Get machine constants.
314 *
315 SAFMIN = DLAMCH( 'Safe minimum' )
316 EPS = DLAMCH( 'Precision' )
317 SMLNUM = SAFMIN / EPS
318 BIGNUM = ONE / SMLNUM
319 RMIN = SQRT( SMLNUM )
320 RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) )
321 *
322 *
323 * Scale matrix to allowable range, if necessary.
324 *
325 ISCALE = 0
326 VLL = VL
327 VUU = VU
328 *
329 TNRM = DLANST( 'M', N, D, E )
330 IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN
331 ISCALE = 1
332 SIGMA = RMIN / TNRM
333 ELSE IF( TNRM.GT.RMAX ) THEN
334 ISCALE = 1
335 SIGMA = RMAX / TNRM
336 END IF
337 IF( ISCALE.EQ.1 ) THEN
338 CALL DSCAL( N, SIGMA, D, 1 )
339 CALL DSCAL( N-1, SIGMA, E( 1 ), 1 )
340 IF( VALEIG ) THEN
341 VLL = VL*SIGMA
342 VUU = VU*SIGMA
343 END IF
344 END IF
345
346 * Initialize indices into workspaces. Note: These indices are used only
347 * if DSTERF or DSTEMR fail.
348
349 * IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in DSTEBZ and
350 * stores the block indices of each of the M<=N eigenvalues.
351 INDIBL = 1
352 * IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in DSTEBZ and
353 * stores the starting and finishing indices of each block.
354 INDISP = INDIBL + N
355 * IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors
356 * that corresponding to eigenvectors that fail to converge in
357 * DSTEIN. This information is discarded; if any fail, the driver
358 * returns INFO > 0.
359 INDIFL = INDISP + N
360 * INDIWO is the offset of the remaining integer workspace.
361 INDIWO = INDISP + N
362 *
363 * If all eigenvalues are desired, then
364 * call DSTERF or DSTEMR. If this fails for some eigenvalue, then
365 * try DSTEBZ.
366 *
367 *
368 TEST = .FALSE.
369 IF( INDEIG ) THEN
370 IF( IL.EQ.1 .AND. IU.EQ.N ) THEN
371 TEST = .TRUE.
372 END IF
373 END IF
374 IF( ( ALLEIG .OR. TEST ) .AND. IEEEOK.EQ.1 ) THEN
375 CALL DCOPY( N-1, E( 1 ), 1, WORK( 1 ), 1 )
376 IF( .NOT.WANTZ ) THEN
377 CALL DCOPY( N, D, 1, W, 1 )
378 CALL DSTERF( N, W, WORK, INFO )
379 ELSE
380 CALL DCOPY( N, D, 1, WORK( N+1 ), 1 )
381 IF (ABSTOL .LE. TWO*N*EPS) THEN
382 TRYRAC = .TRUE.
383 ELSE
384 TRYRAC = .FALSE.
385 END IF
386 CALL DSTEMR( JOBZ, 'A', N, WORK( N+1 ), WORK, VL, VU, IL,
387 $ IU, M, W, Z, LDZ, N, ISUPPZ, TRYRAC,
388 $ WORK( 2*N+1 ), LWORK-2*N, IWORK, LIWORK, INFO )
389 *
390 END IF
391 IF( INFO.EQ.0 ) THEN
392 M = N
393 GO TO 10
394 END IF
395 INFO = 0
396 END IF
397 *
398 * Otherwise, call DSTEBZ and, if eigenvectors are desired, DSTEIN.
399 *
400 IF( WANTZ ) THEN
401 ORDER = 'B'
402 ELSE
403 ORDER = 'E'
404 END IF
405
406 CALL DSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTOL, D, E, M,
407 $ NSPLIT, W, IWORK( INDIBL ), IWORK( INDISP ), WORK,
408 $ IWORK( INDIWO ), INFO )
409 *
410 IF( WANTZ ) THEN
411 CALL DSTEIN( N, D, E, M, W, IWORK( INDIBL ), IWORK( INDISP ),
412 $ Z, LDZ, WORK, IWORK( INDIWO ), IWORK( INDIFL ),
413 $ INFO )
414 END IF
415 *
416 * If matrix was scaled, then rescale eigenvalues appropriately.
417 *
418 10 CONTINUE
419 IF( ISCALE.EQ.1 ) THEN
420 IF( INFO.EQ.0 ) THEN
421 IMAX = M
422 ELSE
423 IMAX = INFO - 1
424 END IF
425 CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
426 END IF
427 *
428 * If eigenvalues are not in order, then sort them, along with
429 * eigenvectors.
430 *
431 IF( WANTZ ) THEN
432 DO 30 J = 1, M - 1
433 I = 0
434 TMP1 = W( J )
435 DO 20 JJ = J + 1, M
436 IF( W( JJ ).LT.TMP1 ) THEN
437 I = JJ
438 TMP1 = W( JJ )
439 END IF
440 20 CONTINUE
441 *
442 IF( I.NE.0 ) THEN
443 ITMP1 = IWORK( I )
444 W( I ) = W( J )
445 IWORK( I ) = IWORK( J )
446 W( J ) = TMP1
447 IWORK( J ) = ITMP1
448 CALL DSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 )
449 END IF
450 30 CONTINUE
451 END IF
452 *
453 * Causes problems with tests 19 & 20:
454 * IF (wantz .and. INDEIG ) Z( 1,1) = Z(1,1) / 1.002 + .002
455 *
456 *
457 WORK( 1 ) = LWMIN
458 IWORK( 1 ) = LIWMIN
459 RETURN
460 *
461 * End of DSTEVR
462 *
463 END