LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
dsbevx.f
Go to the documentation of this file.
1 *> \brief <b> DSBEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices</b>
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DSBEVX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsbevx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsbevx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsbevx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DSBEVX( JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q, LDQ, VL,
22 * VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, IWORK,
23 * IFAIL, INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER JOBZ, RANGE, UPLO
27 * INTEGER IL, INFO, IU, KD, LDAB, LDQ, LDZ, M, N
28 * DOUBLE PRECISION ABSTOL, VL, VU
29 * ..
30 * .. Array Arguments ..
31 * INTEGER IFAIL( * ), IWORK( * )
32 * DOUBLE PRECISION AB( LDAB, * ), Q( LDQ, * ), W( * ), WORK( * ),
33 * $ Z( LDZ, * )
34 * ..
35 *
36 *
37 *> \par Purpose:
38 * =============
39 *>
40 *> \verbatim
41 *>
42 *> DSBEVX computes selected eigenvalues and, optionally, eigenvectors
43 *> of a real symmetric band matrix A. Eigenvalues and eigenvectors can
44 *> be selected by specifying either a range of values or a range of
45 *> indices for the desired eigenvalues.
46 *> \endverbatim
47 *
48 * Arguments:
49 * ==========
50 *
51 *> \param[in] JOBZ
52 *> \verbatim
53 *> JOBZ is CHARACTER*1
54 *> = 'N': Compute eigenvalues only;
55 *> = 'V': Compute eigenvalues and eigenvectors.
56 *> \endverbatim
57 *>
58 *> \param[in] RANGE
59 *> \verbatim
60 *> RANGE is CHARACTER*1
61 *> = 'A': all eigenvalues will be found;
62 *> = 'V': all eigenvalues in the half-open interval (VL,VU]
63 *> will be found;
64 *> = 'I': the IL-th through IU-th eigenvalues will be found.
65 *> \endverbatim
66 *>
67 *> \param[in] UPLO
68 *> \verbatim
69 *> UPLO is CHARACTER*1
70 *> = 'U': Upper triangle of A is stored;
71 *> = 'L': Lower triangle of A is stored.
72 *> \endverbatim
73 *>
74 *> \param[in] N
75 *> \verbatim
76 *> N is INTEGER
77 *> The order of the matrix A. N >= 0.
78 *> \endverbatim
79 *>
80 *> \param[in] KD
81 *> \verbatim
82 *> KD is INTEGER
83 *> The number of superdiagonals of the matrix A if UPLO = 'U',
84 *> or the number of subdiagonals if UPLO = 'L'. KD >= 0.
85 *> \endverbatim
86 *>
87 *> \param[in,out] AB
88 *> \verbatim
89 *> AB is DOUBLE PRECISION array, dimension (LDAB, N)
90 *> On entry, the upper or lower triangle of the symmetric band
91 *> matrix A, stored in the first KD+1 rows of the array. The
92 *> j-th column of A is stored in the j-th column of the array AB
93 *> as follows:
94 *> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
95 *> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
96 *>
97 *> On exit, AB is overwritten by values generated during the
98 *> reduction to tridiagonal form. If UPLO = 'U', the first
99 *> superdiagonal and the diagonal of the tridiagonal matrix T
100 *> are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
101 *> the diagonal and first subdiagonal of T are returned in the
102 *> first two rows of AB.
103 *> \endverbatim
104 *>
105 *> \param[in] LDAB
106 *> \verbatim
107 *> LDAB is INTEGER
108 *> The leading dimension of the array AB. LDAB >= KD + 1.
109 *> \endverbatim
110 *>
111 *> \param[out] Q
112 *> \verbatim
113 *> Q is DOUBLE PRECISION array, dimension (LDQ, N)
114 *> If JOBZ = 'V', the N-by-N orthogonal matrix used in the
115 *> reduction to tridiagonal form.
116 *> If JOBZ = 'N', the array Q is not referenced.
117 *> \endverbatim
118 *>
119 *> \param[in] LDQ
120 *> \verbatim
121 *> LDQ is INTEGER
122 *> The leading dimension of the array Q. If JOBZ = 'V', then
123 *> LDQ >= max(1,N).
124 *> \endverbatim
125 *>
126 *> \param[in] VL
127 *> \verbatim
128 *> VL is DOUBLE PRECISION
129 *> \endverbatim
130 *>
131 *> \param[in] VU
132 *> \verbatim
133 *> VU is DOUBLE PRECISION
134 *> If RANGE='V', the lower and upper bounds of the interval to
135 *> be searched for eigenvalues. VL < VU.
136 *> Not referenced if RANGE = 'A' or 'I'.
137 *> \endverbatim
138 *>
139 *> \param[in] IL
140 *> \verbatim
141 *> IL is INTEGER
142 *> \endverbatim
143 *>
144 *> \param[in] IU
145 *> \verbatim
146 *> IU is INTEGER
147 *> If RANGE='I', the indices (in ascending order) of the
148 *> smallest and largest eigenvalues to be returned.
149 *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
150 *> Not referenced if RANGE = 'A' or 'V'.
151 *> \endverbatim
152 *>
153 *> \param[in] ABSTOL
154 *> \verbatim
155 *> ABSTOL is DOUBLE PRECISION
156 *> The absolute error tolerance for the eigenvalues.
157 *> An approximate eigenvalue is accepted as converged
158 *> when it is determined to lie in an interval [a,b]
159 *> of width less than or equal to
160 *>
161 *> ABSTOL + EPS * max( |a|,|b| ) ,
162 *>
163 *> where EPS is the machine precision. If ABSTOL is less than
164 *> or equal to zero, then EPS*|T| will be used in its place,
165 *> where |T| is the 1-norm of the tridiagonal matrix obtained
166 *> by reducing AB to tridiagonal form.
167 *>
168 *> Eigenvalues will be computed most accurately when ABSTOL is
169 *> set to twice the underflow threshold 2*DLAMCH('S'), not zero.
170 *> If this routine returns with INFO>0, indicating that some
171 *> eigenvectors did not converge, try setting ABSTOL to
172 *> 2*DLAMCH('S').
173 *>
174 *> See "Computing Small Singular Values of Bidiagonal Matrices
175 *> with Guaranteed High Relative Accuracy," by Demmel and
176 *> Kahan, LAPACK Working Note #3.
177 *> \endverbatim
178 *>
179 *> \param[out] M
180 *> \verbatim
181 *> M is INTEGER
182 *> The total number of eigenvalues found. 0 <= M <= N.
183 *> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
184 *> \endverbatim
185 *>
186 *> \param[out] W
187 *> \verbatim
188 *> W is DOUBLE PRECISION array, dimension (N)
189 *> The first M elements contain the selected eigenvalues in
190 *> ascending order.
191 *> \endverbatim
192 *>
193 *> \param[out] Z
194 *> \verbatim
195 *> Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M))
196 *> If JOBZ = 'V', then if INFO = 0, the first M columns of Z
197 *> contain the orthonormal eigenvectors of the matrix A
198 *> corresponding to the selected eigenvalues, with the i-th
199 *> column of Z holding the eigenvector associated with W(i).
200 *> If an eigenvector fails to converge, then that column of Z
201 *> contains the latest approximation to the eigenvector, and the
202 *> index of the eigenvector is returned in IFAIL.
203 *> If JOBZ = 'N', then Z is not referenced.
204 *> Note: the user must ensure that at least max(1,M) columns are
205 *> supplied in the array Z; if RANGE = 'V', the exact value of M
206 *> is not known in advance and an upper bound must be used.
207 *> \endverbatim
208 *>
209 *> \param[in] LDZ
210 *> \verbatim
211 *> LDZ is INTEGER
212 *> The leading dimension of the array Z. LDZ >= 1, and if
213 *> JOBZ = 'V', LDZ >= max(1,N).
214 *> \endverbatim
215 *>
216 *> \param[out] WORK
217 *> \verbatim
218 *> WORK is DOUBLE PRECISION array, dimension (7*N)
219 *> \endverbatim
220 *>
221 *> \param[out] IWORK
222 *> \verbatim
223 *> IWORK is INTEGER array, dimension (5*N)
224 *> \endverbatim
225 *>
226 *> \param[out] IFAIL
227 *> \verbatim
228 *> IFAIL is INTEGER array, dimension (N)
229 *> If JOBZ = 'V', then if INFO = 0, the first M elements of
230 *> IFAIL are zero. If INFO > 0, then IFAIL contains the
231 *> indices of the eigenvectors that failed to converge.
232 *> If JOBZ = 'N', then IFAIL is not referenced.
233 *> \endverbatim
234 *>
235 *> \param[out] INFO
236 *> \verbatim
237 *> INFO is INTEGER
238 *> = 0: successful exit.
239 *> < 0: if INFO = -i, the i-th argument had an illegal value.
240 *> > 0: if INFO = i, then i eigenvectors failed to converge.
241 *> Their indices are stored in array IFAIL.
242 *> \endverbatim
243 *
244 * Authors:
245 * ========
246 *
247 *> \author Univ. of Tennessee
248 *> \author Univ. of California Berkeley
249 *> \author Univ. of Colorado Denver
250 *> \author NAG Ltd.
251 *
252 *> \date November 2011
253 *
254 *> \ingroup doubleOTHEReigen
255 *
256 * =====================================================================
257  SUBROUTINE dsbevx( JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q, LDQ, VL,
258  $ vu, il, iu, abstol, m, w, z, ldz, work, iwork,
259  $ ifail, info )
260 *
261 * -- LAPACK driver routine (version 3.4.0) --
262 * -- LAPACK is a software package provided by Univ. of Tennessee, --
263 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
264 * November 2011
265 *
266 * .. Scalar Arguments ..
267  CHARACTER JOBZ, RANGE, UPLO
268  INTEGER IL, INFO, IU, KD, LDAB, LDQ, LDZ, M, N
269  DOUBLE PRECISION ABSTOL, VL, VU
270 * ..
271 * .. Array Arguments ..
272  INTEGER IFAIL( * ), IWORK( * )
273  DOUBLE PRECISION AB( ldab, * ), Q( ldq, * ), W( * ), WORK( * ),
274  $ z( ldz, * )
275 * ..
276 *
277 * =====================================================================
278 *
279 * .. Parameters ..
280  DOUBLE PRECISION ZERO, ONE
281  parameter( zero = 0.0d0, one = 1.0d0 )
282 * ..
283 * .. Local Scalars ..
284  LOGICAL ALLEIG, INDEIG, LOWER, TEST, VALEIG, WANTZ
285  CHARACTER ORDER
286  INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
287  $ indisp, indiwo, indwrk, iscale, itmp1, j, jj,
288  $ nsplit
289  DOUBLE PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
290  $ sigma, smlnum, tmp1, vll, vuu
291 * ..
292 * .. External Functions ..
293  LOGICAL LSAME
294  DOUBLE PRECISION DLAMCH, DLANSB
295  EXTERNAL lsame, dlamch, dlansb
296 * ..
297 * .. External Subroutines ..
298  EXTERNAL dcopy, dgemv, dlacpy, dlascl, dsbtrd, dscal,
300 * ..
301 * .. Intrinsic Functions ..
302  INTRINSIC max, min, sqrt
303 * ..
304 * .. Executable Statements ..
305 *
306 * Test the input parameters.
307 *
308  wantz = lsame( jobz, 'V' )
309  alleig = lsame( range, 'A' )
310  valeig = lsame( range, 'V' )
311  indeig = lsame( range, 'I' )
312  lower = lsame( uplo, 'L' )
313 *
314  info = 0
315  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
316  info = -1
317  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
318  info = -2
319  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
320  info = -3
321  ELSE IF( n.LT.0 ) THEN
322  info = -4
323  ELSE IF( kd.LT.0 ) THEN
324  info = -5
325  ELSE IF( ldab.LT.kd+1 ) THEN
326  info = -7
327  ELSE IF( wantz .AND. ldq.LT.max( 1, n ) ) THEN
328  info = -9
329  ELSE
330  IF( valeig ) THEN
331  IF( n.GT.0 .AND. vu.LE.vl )
332  $ info = -11
333  ELSE IF( indeig ) THEN
334  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
335  info = -12
336  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
337  info = -13
338  END IF
339  END IF
340  END IF
341  IF( info.EQ.0 ) THEN
342  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
343  $ info = -18
344  END IF
345 *
346  IF( info.NE.0 ) THEN
347  CALL xerbla( 'DSBEVX', -info )
348  RETURN
349  END IF
350 *
351 * Quick return if possible
352 *
353  m = 0
354  IF( n.EQ.0 )
355  $ RETURN
356 *
357  IF( n.EQ.1 ) THEN
358  m = 1
359  IF( lower ) THEN
360  tmp1 = ab( 1, 1 )
361  ELSE
362  tmp1 = ab( kd+1, 1 )
363  END IF
364  IF( valeig ) THEN
365  IF( .NOT.( vl.LT.tmp1 .AND. vu.GE.tmp1 ) )
366  $ m = 0
367  END IF
368  IF( m.EQ.1 ) THEN
369  w( 1 ) = tmp1
370  IF( wantz )
371  $ z( 1, 1 ) = one
372  END IF
373  RETURN
374  END IF
375 *
376 * Get machine constants.
377 *
378  safmin = dlamch( 'Safe minimum' )
379  eps = dlamch( 'Precision' )
380  smlnum = safmin / eps
381  bignum = one / smlnum
382  rmin = sqrt( smlnum )
383  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
384 *
385 * Scale matrix to allowable range, if necessary.
386 *
387  iscale = 0
388  abstll = abstol
389  IF( valeig ) THEN
390  vll = vl
391  vuu = vu
392  ELSE
393  vll = zero
394  vuu = zero
395  END IF
396  anrm = dlansb( 'M', uplo, n, kd, ab, ldab, work )
397  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
398  iscale = 1
399  sigma = rmin / anrm
400  ELSE IF( anrm.GT.rmax ) THEN
401  iscale = 1
402  sigma = rmax / anrm
403  END IF
404  IF( iscale.EQ.1 ) THEN
405  IF( lower ) THEN
406  CALL dlascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
407  ELSE
408  CALL dlascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
409  END IF
410  IF( abstol.GT.0 )
411  $ abstll = abstol*sigma
412  IF( valeig ) THEN
413  vll = vl*sigma
414  vuu = vu*sigma
415  END IF
416  END IF
417 *
418 * Call DSBTRD to reduce symmetric band matrix to tridiagonal form.
419 *
420  indd = 1
421  inde = indd + n
422  indwrk = inde + n
423  CALL dsbtrd( jobz, uplo, n, kd, ab, ldab, work( indd ),
424  $ work( inde ), q, ldq, work( indwrk ), iinfo )
425 *
426 * If all eigenvalues are desired and ABSTOL is less than or equal
427 * to zero, then call DSTERF or SSTEQR. If this fails for some
428 * eigenvalue, then try DSTEBZ.
429 *
430  test = .false.
431  IF (indeig) THEN
432  IF (il.EQ.1 .AND. iu.EQ.n) THEN
433  test = .true.
434  END IF
435  END IF
436  IF ((alleig .OR. test) .AND. (abstol.LE.zero)) THEN
437  CALL dcopy( n, work( indd ), 1, w, 1 )
438  indee = indwrk + 2*n
439  IF( .NOT.wantz ) THEN
440  CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
441  CALL dsterf( n, w, work( indee ), info )
442  ELSE
443  CALL dlacpy( 'A', n, n, q, ldq, z, ldz )
444  CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
445  CALL dsteqr( jobz, n, w, work( indee ), z, ldz,
446  $ work( indwrk ), info )
447  IF( info.EQ.0 ) THEN
448  DO 10 i = 1, n
449  ifail( i ) = 0
450  10 CONTINUE
451  END IF
452  END IF
453  IF( info.EQ.0 ) THEN
454  m = n
455  GO TO 30
456  END IF
457  info = 0
458  END IF
459 *
460 * Otherwise, call DSTEBZ and, if eigenvectors are desired, SSTEIN.
461 *
462  IF( wantz ) THEN
463  order = 'B'
464  ELSE
465  order = 'E'
466  END IF
467  indibl = 1
468  indisp = indibl + n
469  indiwo = indisp + n
470  CALL dstebz( range, order, n, vll, vuu, il, iu, abstll,
471  $ work( indd ), work( inde ), m, nsplit, w,
472  $ iwork( indibl ), iwork( indisp ), work( indwrk ),
473  $ iwork( indiwo ), info )
474 *
475  IF( wantz ) THEN
476  CALL dstein( n, work( indd ), work( inde ), m, w,
477  $ iwork( indibl ), iwork( indisp ), z, ldz,
478  $ work( indwrk ), iwork( indiwo ), ifail, info )
479 *
480 * Apply orthogonal matrix used in reduction to tridiagonal
481 * form to eigenvectors returned by DSTEIN.
482 *
483  DO 20 j = 1, m
484  CALL dcopy( n, z( 1, j ), 1, work( 1 ), 1 )
485  CALL dgemv( 'N', n, n, one, q, ldq, work, 1, zero,
486  $ z( 1, j ), 1 )
487  20 CONTINUE
488  END IF
489 *
490 * If matrix was scaled, then rescale eigenvalues appropriately.
491 *
492  30 CONTINUE
493  IF( iscale.EQ.1 ) THEN
494  IF( info.EQ.0 ) THEN
495  imax = m
496  ELSE
497  imax = info - 1
498  END IF
499  CALL dscal( imax, one / sigma, w, 1 )
500  END IF
501 *
502 * If eigenvalues are not in order, then sort them, along with
503 * eigenvectors.
504 *
505  IF( wantz ) THEN
506  DO 50 j = 1, m - 1
507  i = 0
508  tmp1 = w( j )
509  DO 40 jj = j + 1, m
510  IF( w( jj ).LT.tmp1 ) THEN
511  i = jj
512  tmp1 = w( jj )
513  END IF
514  40 CONTINUE
515 *
516  IF( i.NE.0 ) THEN
517  itmp1 = iwork( indibl+i-1 )
518  w( i ) = w( j )
519  iwork( indibl+i-1 ) = iwork( indibl+j-1 )
520  w( j ) = tmp1
521  iwork( indibl+j-1 ) = itmp1
522  CALL dswap( n, z( 1, i ), 1, z( 1, j ), 1 )
523  IF( info.NE.0 ) THEN
524  itmp1 = ifail( i )
525  ifail( i ) = ifail( j )
526  ifail( j ) = itmp1
527  END IF
528  END IF
529  50 CONTINUE
530  END IF
531 *
532  RETURN
533 *
534 * End of DSBEVX
535 *
536  END
subroutine dsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
DSTEQR
Definition: dsteqr.f:133
subroutine dstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
DSTEIN
Definition: dstein.f:176
subroutine dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:88
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
subroutine dsbtrd(VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ, WORK, INFO)
DSBTRD
Definition: dsbtrd.f:165
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: dlascl.f:141
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55
subroutine dsbevx(JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q, LDQ, VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, IWORK, IFAIL, INFO)
DSBEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matric...
Definition: dsbevx.f:260
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:53
subroutine dstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
DSTEBZ
Definition: dstebz.f:265
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:158