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