1 SUBROUTINE DSBEVX( JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q, LDQ, VL,
2 $ VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, IWORK,
3 $ IFAIL, 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, UPLO
12 INTEGER IL, INFO, IU, KD, LDAB, LDQ, LDZ, M, N
13 DOUBLE PRECISION ABSTOL, VL, VU
14 * ..
15 * .. Array Arguments ..
16 INTEGER IFAIL( * ), IWORK( * )
17 DOUBLE PRECISION AB( LDAB, * ), Q( LDQ, * ), W( * ), WORK( * ),
18 $ Z( LDZ, * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * DSBEVX computes selected eigenvalues and, optionally, eigenvectors
25 * of a real symmetric band matrix A. Eigenvalues and eigenvectors can
26 * be selected by specifying either a range of values or a range of
27 * indices for the desired eigenvalues.
28 *
29 * Arguments
30 * =========
31 *
32 * JOBZ (input) CHARACTER*1
33 * = 'N': Compute eigenvalues only;
34 * = 'V': Compute eigenvalues and eigenvectors.
35 *
36 * RANGE (input) CHARACTER*1
37 * = 'A': all eigenvalues will be found;
38 * = 'V': all eigenvalues in the half-open interval (VL,VU]
39 * will be found;
40 * = 'I': the IL-th through IU-th eigenvalues will be found.
41 *
42 * UPLO (input) CHARACTER*1
43 * = 'U': Upper triangle of A is stored;
44 * = 'L': Lower triangle of A is stored.
45 *
46 * N (input) INTEGER
47 * The order of the matrix A. N >= 0.
48 *
49 * KD (input) INTEGER
50 * The number of superdiagonals of the matrix A if UPLO = 'U',
51 * or the number of subdiagonals if UPLO = 'L'. KD >= 0.
52 *
53 * AB (input/output) DOUBLE PRECISION array, dimension (LDAB, N)
54 * On entry, the upper or lower triangle of the symmetric band
55 * matrix A, stored in the first KD+1 rows of the array. The
56 * j-th column of A is stored in the j-th column of the array AB
57 * as follows:
58 * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
59 * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
60 *
61 * On exit, AB is overwritten by values generated during the
62 * reduction to tridiagonal form. If UPLO = 'U', the first
63 * superdiagonal and the diagonal of the tridiagonal matrix T
64 * are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
65 * the diagonal and first subdiagonal of T are returned in the
66 * first two rows of AB.
67 *
68 * LDAB (input) INTEGER
69 * The leading dimension of the array AB. LDAB >= KD + 1.
70 *
71 * Q (output) DOUBLE PRECISION array, dimension (LDQ, N)
72 * If JOBZ = 'V', the N-by-N orthogonal matrix used in the
73 * reduction to tridiagonal form.
74 * If JOBZ = 'N', the array Q is not referenced.
75 *
76 * LDQ (input) INTEGER
77 * The leading dimension of the array Q. If JOBZ = 'V', then
78 * LDQ >= max(1,N).
79 *
80 * VL (input) DOUBLE PRECISION
81 * VU (input) DOUBLE PRECISION
82 * If RANGE='V', the lower and upper bounds of the interval to
83 * be searched for eigenvalues. VL < VU.
84 * Not referenced if RANGE = 'A' or 'I'.
85 *
86 * IL (input) INTEGER
87 * IU (input) INTEGER
88 * If RANGE='I', the indices (in ascending order) of the
89 * smallest and largest eigenvalues to be returned.
90 * 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
91 * Not referenced if RANGE = 'A' or 'V'.
92 *
93 * ABSTOL (input) DOUBLE PRECISION
94 * The absolute error tolerance for the eigenvalues.
95 * An approximate eigenvalue is accepted as converged
96 * when it is determined to lie in an interval [a,b]
97 * of width less than or equal to
98 *
99 * ABSTOL + EPS * max( |a|,|b| ) ,
100 *
101 * where EPS is the machine precision. If ABSTOL is less than
102 * or equal to zero, then EPS*|T| will be used in its place,
103 * where |T| is the 1-norm of the tridiagonal matrix obtained
104 * by reducing AB to tridiagonal form.
105 *
106 * Eigenvalues will be computed most accurately when ABSTOL is
107 * set to twice the underflow threshold 2*DLAMCH('S'), not zero.
108 * If this routine returns with INFO>0, indicating that some
109 * eigenvectors did not converge, try setting ABSTOL to
110 * 2*DLAMCH('S').
111 *
112 * See "Computing Small Singular Values of Bidiagonal Matrices
113 * with Guaranteed High Relative Accuracy," by Demmel and
114 * Kahan, LAPACK Working Note #3.
115 *
116 * M (output) INTEGER
117 * The total number of eigenvalues found. 0 <= M <= N.
118 * If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
119 *
120 * W (output) DOUBLE PRECISION array, dimension (N)
121 * The first M elements contain the selected eigenvalues in
122 * ascending order.
123 *
124 * Z (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M))
125 * If JOBZ = 'V', then if INFO = 0, the first M columns of Z
126 * contain the orthonormal eigenvectors of the matrix A
127 * corresponding to the selected eigenvalues, with the i-th
128 * column of Z holding the eigenvector associated with W(i).
129 * If an eigenvector fails to converge, then that column of Z
130 * contains the latest approximation to the eigenvector, and the
131 * index of the eigenvector is returned in IFAIL.
132 * If JOBZ = 'N', then Z is not referenced.
133 * Note: the user must ensure that at least max(1,M) columns are
134 * supplied in the array Z; if RANGE = 'V', the exact value of M
135 * is not known in advance and an upper bound must be used.
136 *
137 * LDZ (input) INTEGER
138 * The leading dimension of the array Z. LDZ >= 1, and if
139 * JOBZ = 'V', LDZ >= max(1,N).
140 *
141 * WORK (workspace) DOUBLE PRECISION array, dimension (7*N)
142 *
143 * IWORK (workspace) INTEGER array, dimension (5*N)
144 *
145 * IFAIL (output) INTEGER array, dimension (N)
146 * If JOBZ = 'V', then if INFO = 0, the first M elements of
147 * IFAIL are zero. If INFO > 0, then IFAIL contains the
148 * indices of the eigenvectors that failed to converge.
149 * If JOBZ = 'N', then IFAIL is not referenced.
150 *
151 * INFO (output) INTEGER
152 * = 0: successful exit.
153 * < 0: if INFO = -i, the i-th argument had an illegal value.
154 * > 0: if INFO = i, then i eigenvectors failed to converge.
155 * Their indices are stored in array IFAIL.
156 *
157 * =====================================================================
158 *
159 * .. Parameters ..
160 DOUBLE PRECISION ZERO, ONE
161 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
162 * ..
163 * .. Local Scalars ..
164 LOGICAL ALLEIG, INDEIG, LOWER, TEST, VALEIG, WANTZ
165 CHARACTER ORDER
166 INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
167 $ INDISP, INDIWO, INDWRK, ISCALE, ITMP1, J, JJ,
168 $ NSPLIT
169 DOUBLE PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
170 $ SIGMA, SMLNUM, TMP1, VLL, VUU
171 * ..
172 * .. External Functions ..
173 LOGICAL LSAME
174 DOUBLE PRECISION DLAMCH, DLANSB
175 EXTERNAL LSAME, DLAMCH, DLANSB
176 * ..
177 * .. External Subroutines ..
178 EXTERNAL DCOPY, DGEMV, DLACPY, DLASCL, DSBTRD, DSCAL,
179 $ DSTEBZ, DSTEIN, DSTEQR, DSTERF, DSWAP, XERBLA
180 * ..
181 * .. Intrinsic Functions ..
182 INTRINSIC MAX, MIN, SQRT
183 * ..
184 * .. Executable Statements ..
185 *
186 * Test the input parameters.
187 *
188 WANTZ = LSAME( JOBZ, 'V' )
189 ALLEIG = LSAME( RANGE, 'A' )
190 VALEIG = LSAME( RANGE, 'V' )
191 INDEIG = LSAME( RANGE, 'I' )
192 LOWER = LSAME( UPLO, 'L' )
193 *
194 INFO = 0
195 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
196 INFO = -1
197 ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN
198 INFO = -2
199 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
200 INFO = -3
201 ELSE IF( N.LT.0 ) THEN
202 INFO = -4
203 ELSE IF( KD.LT.0 ) THEN
204 INFO = -5
205 ELSE IF( LDAB.LT.KD+1 ) THEN
206 INFO = -7
207 ELSE IF( WANTZ .AND. LDQ.LT.MAX( 1, N ) ) THEN
208 INFO = -9
209 ELSE
210 IF( VALEIG ) THEN
211 IF( N.GT.0 .AND. VU.LE.VL )
212 $ INFO = -11
213 ELSE IF( INDEIG ) THEN
214 IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN
215 INFO = -12
216 ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN
217 INFO = -13
218 END IF
219 END IF
220 END IF
221 IF( INFO.EQ.0 ) THEN
222 IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) )
223 $ INFO = -18
224 END IF
225 *
226 IF( INFO.NE.0 ) THEN
227 CALL XERBLA( 'DSBEVX', -INFO )
228 RETURN
229 END IF
230 *
231 * Quick return if possible
232 *
233 M = 0
234 IF( N.EQ.0 )
235 $ RETURN
236 *
237 IF( N.EQ.1 ) THEN
238 M = 1
239 IF( LOWER ) THEN
240 TMP1 = AB( 1, 1 )
241 ELSE
242 TMP1 = AB( KD+1, 1 )
243 END IF
244 IF( VALEIG ) THEN
245 IF( .NOT.( VL.LT.TMP1 .AND. VU.GE.TMP1 ) )
246 $ M = 0
247 END IF
248 IF( M.EQ.1 ) THEN
249 W( 1 ) = TMP1
250 IF( WANTZ )
251 $ Z( 1, 1 ) = ONE
252 END IF
253 RETURN
254 END IF
255 *
256 * Get machine constants.
257 *
258 SAFMIN = DLAMCH( 'Safe minimum' )
259 EPS = DLAMCH( 'Precision' )
260 SMLNUM = SAFMIN / EPS
261 BIGNUM = ONE / SMLNUM
262 RMIN = SQRT( SMLNUM )
263 RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) )
264 *
265 * Scale matrix to allowable range, if necessary.
266 *
267 ISCALE = 0
268 ABSTLL = ABSTOL
269 IF( VALEIG ) THEN
270 VLL = VL
271 VUU = VU
272 ELSE
273 VLL = ZERO
274 VUU = ZERO
275 END IF
276 ANRM = DLANSB( 'M', UPLO, N, KD, AB, LDAB, WORK )
277 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
278 ISCALE = 1
279 SIGMA = RMIN / ANRM
280 ELSE IF( ANRM.GT.RMAX ) THEN
281 ISCALE = 1
282 SIGMA = RMAX / ANRM
283 END IF
284 IF( ISCALE.EQ.1 ) THEN
285 IF( LOWER ) THEN
286 CALL DLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
287 ELSE
288 CALL DLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
289 END IF
290 IF( ABSTOL.GT.0 )
291 $ ABSTLL = ABSTOL*SIGMA
292 IF( VALEIG ) THEN
293 VLL = VL*SIGMA
294 VUU = VU*SIGMA
295 END IF
296 END IF
297 *
298 * Call DSBTRD to reduce symmetric band matrix to tridiagonal form.
299 *
300 INDD = 1
301 INDE = INDD + N
302 INDWRK = INDE + N
303 CALL DSBTRD( JOBZ, UPLO, N, KD, AB, LDAB, WORK( INDD ),
304 $ WORK( INDE ), Q, LDQ, WORK( INDWRK ), IINFO )
305 *
306 * If all eigenvalues are desired and ABSTOL is less than or equal
307 * to zero, then call DSTERF or SSTEQR. If this fails for some
308 * eigenvalue, then try DSTEBZ.
309 *
310 TEST = .FALSE.
311 IF (INDEIG) THEN
312 IF (IL.EQ.1 .AND. IU.EQ.N) THEN
313 TEST = .TRUE.
314 END IF
315 END IF
316 IF ((ALLEIG .OR. TEST) .AND. (ABSTOL.LE.ZERO)) THEN
317 CALL DCOPY( N, WORK( INDD ), 1, W, 1 )
318 INDEE = INDWRK + 2*N
319 IF( .NOT.WANTZ ) THEN
320 CALL DCOPY( N-1, WORK( INDE ), 1, WORK( INDEE ), 1 )
321 CALL DSTERF( N, W, WORK( INDEE ), INFO )
322 ELSE
323 CALL DLACPY( 'A', N, N, Q, LDQ, Z, LDZ )
324 CALL DCOPY( N-1, WORK( INDE ), 1, WORK( INDEE ), 1 )
325 CALL DSTEQR( JOBZ, N, W, WORK( INDEE ), Z, LDZ,
326 $ WORK( INDWRK ), INFO )
327 IF( INFO.EQ.0 ) THEN
328 DO 10 I = 1, N
329 IFAIL( I ) = 0
330 10 CONTINUE
331 END IF
332 END IF
333 IF( INFO.EQ.0 ) THEN
334 M = N
335 GO TO 30
336 END IF
337 INFO = 0
338 END IF
339 *
340 * Otherwise, call DSTEBZ and, if eigenvectors are desired, SSTEIN.
341 *
342 IF( WANTZ ) THEN
343 ORDER = 'B'
344 ELSE
345 ORDER = 'E'
346 END IF
347 INDIBL = 1
348 INDISP = INDIBL + N
349 INDIWO = INDISP + N
350 CALL DSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTLL,
351 $ WORK( INDD ), WORK( INDE ), M, NSPLIT, W,
352 $ IWORK( INDIBL ), IWORK( INDISP ), WORK( INDWRK ),
353 $ IWORK( INDIWO ), INFO )
354 *
355 IF( WANTZ ) THEN
356 CALL DSTEIN( N, WORK( INDD ), WORK( INDE ), M, W,
357 $ IWORK( INDIBL ), IWORK( INDISP ), Z, LDZ,
358 $ WORK( INDWRK ), IWORK( INDIWO ), IFAIL, INFO )
359 *
360 * Apply orthogonal matrix used in reduction to tridiagonal
361 * form to eigenvectors returned by DSTEIN.
362 *
363 DO 20 J = 1, M
364 CALL DCOPY( N, Z( 1, J ), 1, WORK( 1 ), 1 )
365 CALL DGEMV( 'N', N, N, ONE, Q, LDQ, WORK, 1, ZERO,
366 $ Z( 1, J ), 1 )
367 20 CONTINUE
368 END IF
369 *
370 * If matrix was scaled, then rescale eigenvalues appropriately.
371 *
372 30 CONTINUE
373 IF( ISCALE.EQ.1 ) THEN
374 IF( INFO.EQ.0 ) THEN
375 IMAX = M
376 ELSE
377 IMAX = INFO - 1
378 END IF
379 CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
380 END IF
381 *
382 * If eigenvalues are not in order, then sort them, along with
383 * eigenvectors.
384 *
385 IF( WANTZ ) THEN
386 DO 50 J = 1, M - 1
387 I = 0
388 TMP1 = W( J )
389 DO 40 JJ = J + 1, M
390 IF( W( JJ ).LT.TMP1 ) THEN
391 I = JJ
392 TMP1 = W( JJ )
393 END IF
394 40 CONTINUE
395 *
396 IF( I.NE.0 ) THEN
397 ITMP1 = IWORK( INDIBL+I-1 )
398 W( I ) = W( J )
399 IWORK( INDIBL+I-1 ) = IWORK( INDIBL+J-1 )
400 W( J ) = TMP1
401 IWORK( INDIBL+J-1 ) = ITMP1
402 CALL DSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 )
403 IF( INFO.NE.0 ) THEN
404 ITMP1 = IFAIL( I )
405 IFAIL( I ) = IFAIL( J )
406 IFAIL( J ) = ITMP1
407 END IF
408 END IF
409 50 CONTINUE
410 END IF
411 *
412 RETURN
413 *
414 * End of DSBEVX
415 *
416 END
2 $ VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, IWORK,
3 $ IFAIL, 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, UPLO
12 INTEGER IL, INFO, IU, KD, LDAB, LDQ, LDZ, M, N
13 DOUBLE PRECISION ABSTOL, VL, VU
14 * ..
15 * .. Array Arguments ..
16 INTEGER IFAIL( * ), IWORK( * )
17 DOUBLE PRECISION AB( LDAB, * ), Q( LDQ, * ), W( * ), WORK( * ),
18 $ Z( LDZ, * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * DSBEVX computes selected eigenvalues and, optionally, eigenvectors
25 * of a real symmetric band matrix A. Eigenvalues and eigenvectors can
26 * be selected by specifying either a range of values or a range of
27 * indices for the desired eigenvalues.
28 *
29 * Arguments
30 * =========
31 *
32 * JOBZ (input) CHARACTER*1
33 * = 'N': Compute eigenvalues only;
34 * = 'V': Compute eigenvalues and eigenvectors.
35 *
36 * RANGE (input) CHARACTER*1
37 * = 'A': all eigenvalues will be found;
38 * = 'V': all eigenvalues in the half-open interval (VL,VU]
39 * will be found;
40 * = 'I': the IL-th through IU-th eigenvalues will be found.
41 *
42 * UPLO (input) CHARACTER*1
43 * = 'U': Upper triangle of A is stored;
44 * = 'L': Lower triangle of A is stored.
45 *
46 * N (input) INTEGER
47 * The order of the matrix A. N >= 0.
48 *
49 * KD (input) INTEGER
50 * The number of superdiagonals of the matrix A if UPLO = 'U',
51 * or the number of subdiagonals if UPLO = 'L'. KD >= 0.
52 *
53 * AB (input/output) DOUBLE PRECISION array, dimension (LDAB, N)
54 * On entry, the upper or lower triangle of the symmetric band
55 * matrix A, stored in the first KD+1 rows of the array. The
56 * j-th column of A is stored in the j-th column of the array AB
57 * as follows:
58 * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
59 * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
60 *
61 * On exit, AB is overwritten by values generated during the
62 * reduction to tridiagonal form. If UPLO = 'U', the first
63 * superdiagonal and the diagonal of the tridiagonal matrix T
64 * are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
65 * the diagonal and first subdiagonal of T are returned in the
66 * first two rows of AB.
67 *
68 * LDAB (input) INTEGER
69 * The leading dimension of the array AB. LDAB >= KD + 1.
70 *
71 * Q (output) DOUBLE PRECISION array, dimension (LDQ, N)
72 * If JOBZ = 'V', the N-by-N orthogonal matrix used in the
73 * reduction to tridiagonal form.
74 * If JOBZ = 'N', the array Q is not referenced.
75 *
76 * LDQ (input) INTEGER
77 * The leading dimension of the array Q. If JOBZ = 'V', then
78 * LDQ >= max(1,N).
79 *
80 * VL (input) DOUBLE PRECISION
81 * VU (input) DOUBLE PRECISION
82 * If RANGE='V', the lower and upper bounds of the interval to
83 * be searched for eigenvalues. VL < VU.
84 * Not referenced if RANGE = 'A' or 'I'.
85 *
86 * IL (input) INTEGER
87 * IU (input) INTEGER
88 * If RANGE='I', the indices (in ascending order) of the
89 * smallest and largest eigenvalues to be returned.
90 * 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
91 * Not referenced if RANGE = 'A' or 'V'.
92 *
93 * ABSTOL (input) DOUBLE PRECISION
94 * The absolute error tolerance for the eigenvalues.
95 * An approximate eigenvalue is accepted as converged
96 * when it is determined to lie in an interval [a,b]
97 * of width less than or equal to
98 *
99 * ABSTOL + EPS * max( |a|,|b| ) ,
100 *
101 * where EPS is the machine precision. If ABSTOL is less than
102 * or equal to zero, then EPS*|T| will be used in its place,
103 * where |T| is the 1-norm of the tridiagonal matrix obtained
104 * by reducing AB to tridiagonal form.
105 *
106 * Eigenvalues will be computed most accurately when ABSTOL is
107 * set to twice the underflow threshold 2*DLAMCH('S'), not zero.
108 * If this routine returns with INFO>0, indicating that some
109 * eigenvectors did not converge, try setting ABSTOL to
110 * 2*DLAMCH('S').
111 *
112 * See "Computing Small Singular Values of Bidiagonal Matrices
113 * with Guaranteed High Relative Accuracy," by Demmel and
114 * Kahan, LAPACK Working Note #3.
115 *
116 * M (output) INTEGER
117 * The total number of eigenvalues found. 0 <= M <= N.
118 * If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
119 *
120 * W (output) DOUBLE PRECISION array, dimension (N)
121 * The first M elements contain the selected eigenvalues in
122 * ascending order.
123 *
124 * Z (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M))
125 * If JOBZ = 'V', then if INFO = 0, the first M columns of Z
126 * contain the orthonormal eigenvectors of the matrix A
127 * corresponding to the selected eigenvalues, with the i-th
128 * column of Z holding the eigenvector associated with W(i).
129 * If an eigenvector fails to converge, then that column of Z
130 * contains the latest approximation to the eigenvector, and the
131 * index of the eigenvector is returned in IFAIL.
132 * If JOBZ = 'N', then Z is not referenced.
133 * Note: the user must ensure that at least max(1,M) columns are
134 * supplied in the array Z; if RANGE = 'V', the exact value of M
135 * is not known in advance and an upper bound must be used.
136 *
137 * LDZ (input) INTEGER
138 * The leading dimension of the array Z. LDZ >= 1, and if
139 * JOBZ = 'V', LDZ >= max(1,N).
140 *
141 * WORK (workspace) DOUBLE PRECISION array, dimension (7*N)
142 *
143 * IWORK (workspace) INTEGER array, dimension (5*N)
144 *
145 * IFAIL (output) INTEGER array, dimension (N)
146 * If JOBZ = 'V', then if INFO = 0, the first M elements of
147 * IFAIL are zero. If INFO > 0, then IFAIL contains the
148 * indices of the eigenvectors that failed to converge.
149 * If JOBZ = 'N', then IFAIL is not referenced.
150 *
151 * INFO (output) INTEGER
152 * = 0: successful exit.
153 * < 0: if INFO = -i, the i-th argument had an illegal value.
154 * > 0: if INFO = i, then i eigenvectors failed to converge.
155 * Their indices are stored in array IFAIL.
156 *
157 * =====================================================================
158 *
159 * .. Parameters ..
160 DOUBLE PRECISION ZERO, ONE
161 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
162 * ..
163 * .. Local Scalars ..
164 LOGICAL ALLEIG, INDEIG, LOWER, TEST, VALEIG, WANTZ
165 CHARACTER ORDER
166 INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
167 $ INDISP, INDIWO, INDWRK, ISCALE, ITMP1, J, JJ,
168 $ NSPLIT
169 DOUBLE PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
170 $ SIGMA, SMLNUM, TMP1, VLL, VUU
171 * ..
172 * .. External Functions ..
173 LOGICAL LSAME
174 DOUBLE PRECISION DLAMCH, DLANSB
175 EXTERNAL LSAME, DLAMCH, DLANSB
176 * ..
177 * .. External Subroutines ..
178 EXTERNAL DCOPY, DGEMV, DLACPY, DLASCL, DSBTRD, DSCAL,
179 $ DSTEBZ, DSTEIN, DSTEQR, DSTERF, DSWAP, XERBLA
180 * ..
181 * .. Intrinsic Functions ..
182 INTRINSIC MAX, MIN, SQRT
183 * ..
184 * .. Executable Statements ..
185 *
186 * Test the input parameters.
187 *
188 WANTZ = LSAME( JOBZ, 'V' )
189 ALLEIG = LSAME( RANGE, 'A' )
190 VALEIG = LSAME( RANGE, 'V' )
191 INDEIG = LSAME( RANGE, 'I' )
192 LOWER = LSAME( UPLO, 'L' )
193 *
194 INFO = 0
195 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
196 INFO = -1
197 ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN
198 INFO = -2
199 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
200 INFO = -3
201 ELSE IF( N.LT.0 ) THEN
202 INFO = -4
203 ELSE IF( KD.LT.0 ) THEN
204 INFO = -5
205 ELSE IF( LDAB.LT.KD+1 ) THEN
206 INFO = -7
207 ELSE IF( WANTZ .AND. LDQ.LT.MAX( 1, N ) ) THEN
208 INFO = -9
209 ELSE
210 IF( VALEIG ) THEN
211 IF( N.GT.0 .AND. VU.LE.VL )
212 $ INFO = -11
213 ELSE IF( INDEIG ) THEN
214 IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN
215 INFO = -12
216 ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN
217 INFO = -13
218 END IF
219 END IF
220 END IF
221 IF( INFO.EQ.0 ) THEN
222 IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) )
223 $ INFO = -18
224 END IF
225 *
226 IF( INFO.NE.0 ) THEN
227 CALL XERBLA( 'DSBEVX', -INFO )
228 RETURN
229 END IF
230 *
231 * Quick return if possible
232 *
233 M = 0
234 IF( N.EQ.0 )
235 $ RETURN
236 *
237 IF( N.EQ.1 ) THEN
238 M = 1
239 IF( LOWER ) THEN
240 TMP1 = AB( 1, 1 )
241 ELSE
242 TMP1 = AB( KD+1, 1 )
243 END IF
244 IF( VALEIG ) THEN
245 IF( .NOT.( VL.LT.TMP1 .AND. VU.GE.TMP1 ) )
246 $ M = 0
247 END IF
248 IF( M.EQ.1 ) THEN
249 W( 1 ) = TMP1
250 IF( WANTZ )
251 $ Z( 1, 1 ) = ONE
252 END IF
253 RETURN
254 END IF
255 *
256 * Get machine constants.
257 *
258 SAFMIN = DLAMCH( 'Safe minimum' )
259 EPS = DLAMCH( 'Precision' )
260 SMLNUM = SAFMIN / EPS
261 BIGNUM = ONE / SMLNUM
262 RMIN = SQRT( SMLNUM )
263 RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) )
264 *
265 * Scale matrix to allowable range, if necessary.
266 *
267 ISCALE = 0
268 ABSTLL = ABSTOL
269 IF( VALEIG ) THEN
270 VLL = VL
271 VUU = VU
272 ELSE
273 VLL = ZERO
274 VUU = ZERO
275 END IF
276 ANRM = DLANSB( 'M', UPLO, N, KD, AB, LDAB, WORK )
277 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
278 ISCALE = 1
279 SIGMA = RMIN / ANRM
280 ELSE IF( ANRM.GT.RMAX ) THEN
281 ISCALE = 1
282 SIGMA = RMAX / ANRM
283 END IF
284 IF( ISCALE.EQ.1 ) THEN
285 IF( LOWER ) THEN
286 CALL DLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
287 ELSE
288 CALL DLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
289 END IF
290 IF( ABSTOL.GT.0 )
291 $ ABSTLL = ABSTOL*SIGMA
292 IF( VALEIG ) THEN
293 VLL = VL*SIGMA
294 VUU = VU*SIGMA
295 END IF
296 END IF
297 *
298 * Call DSBTRD to reduce symmetric band matrix to tridiagonal form.
299 *
300 INDD = 1
301 INDE = INDD + N
302 INDWRK = INDE + N
303 CALL DSBTRD( JOBZ, UPLO, N, KD, AB, LDAB, WORK( INDD ),
304 $ WORK( INDE ), Q, LDQ, WORK( INDWRK ), IINFO )
305 *
306 * If all eigenvalues are desired and ABSTOL is less than or equal
307 * to zero, then call DSTERF or SSTEQR. If this fails for some
308 * eigenvalue, then try DSTEBZ.
309 *
310 TEST = .FALSE.
311 IF (INDEIG) THEN
312 IF (IL.EQ.1 .AND. IU.EQ.N) THEN
313 TEST = .TRUE.
314 END IF
315 END IF
316 IF ((ALLEIG .OR. TEST) .AND. (ABSTOL.LE.ZERO)) THEN
317 CALL DCOPY( N, WORK( INDD ), 1, W, 1 )
318 INDEE = INDWRK + 2*N
319 IF( .NOT.WANTZ ) THEN
320 CALL DCOPY( N-1, WORK( INDE ), 1, WORK( INDEE ), 1 )
321 CALL DSTERF( N, W, WORK( INDEE ), INFO )
322 ELSE
323 CALL DLACPY( 'A', N, N, Q, LDQ, Z, LDZ )
324 CALL DCOPY( N-1, WORK( INDE ), 1, WORK( INDEE ), 1 )
325 CALL DSTEQR( JOBZ, N, W, WORK( INDEE ), Z, LDZ,
326 $ WORK( INDWRK ), INFO )
327 IF( INFO.EQ.0 ) THEN
328 DO 10 I = 1, N
329 IFAIL( I ) = 0
330 10 CONTINUE
331 END IF
332 END IF
333 IF( INFO.EQ.0 ) THEN
334 M = N
335 GO TO 30
336 END IF
337 INFO = 0
338 END IF
339 *
340 * Otherwise, call DSTEBZ and, if eigenvectors are desired, SSTEIN.
341 *
342 IF( WANTZ ) THEN
343 ORDER = 'B'
344 ELSE
345 ORDER = 'E'
346 END IF
347 INDIBL = 1
348 INDISP = INDIBL + N
349 INDIWO = INDISP + N
350 CALL DSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTLL,
351 $ WORK( INDD ), WORK( INDE ), M, NSPLIT, W,
352 $ IWORK( INDIBL ), IWORK( INDISP ), WORK( INDWRK ),
353 $ IWORK( INDIWO ), INFO )
354 *
355 IF( WANTZ ) THEN
356 CALL DSTEIN( N, WORK( INDD ), WORK( INDE ), M, W,
357 $ IWORK( INDIBL ), IWORK( INDISP ), Z, LDZ,
358 $ WORK( INDWRK ), IWORK( INDIWO ), IFAIL, INFO )
359 *
360 * Apply orthogonal matrix used in reduction to tridiagonal
361 * form to eigenvectors returned by DSTEIN.
362 *
363 DO 20 J = 1, M
364 CALL DCOPY( N, Z( 1, J ), 1, WORK( 1 ), 1 )
365 CALL DGEMV( 'N', N, N, ONE, Q, LDQ, WORK, 1, ZERO,
366 $ Z( 1, J ), 1 )
367 20 CONTINUE
368 END IF
369 *
370 * If matrix was scaled, then rescale eigenvalues appropriately.
371 *
372 30 CONTINUE
373 IF( ISCALE.EQ.1 ) THEN
374 IF( INFO.EQ.0 ) THEN
375 IMAX = M
376 ELSE
377 IMAX = INFO - 1
378 END IF
379 CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
380 END IF
381 *
382 * If eigenvalues are not in order, then sort them, along with
383 * eigenvectors.
384 *
385 IF( WANTZ ) THEN
386 DO 50 J = 1, M - 1
387 I = 0
388 TMP1 = W( J )
389 DO 40 JJ = J + 1, M
390 IF( W( JJ ).LT.TMP1 ) THEN
391 I = JJ
392 TMP1 = W( JJ )
393 END IF
394 40 CONTINUE
395 *
396 IF( I.NE.0 ) THEN
397 ITMP1 = IWORK( INDIBL+I-1 )
398 W( I ) = W( J )
399 IWORK( INDIBL+I-1 ) = IWORK( INDIBL+J-1 )
400 W( J ) = TMP1
401 IWORK( INDIBL+J-1 ) = ITMP1
402 CALL DSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 )
403 IF( INFO.NE.0 ) THEN
404 ITMP1 = IFAIL( I )
405 IFAIL( I ) = IFAIL( J )
406 IFAIL( J ) = ITMP1
407 END IF
408 END IF
409 50 CONTINUE
410 END IF
411 *
412 RETURN
413 *
414 * End of DSBEVX
415 *
416 END