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