LAPACK  3.6.1 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
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.6.1) --
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 dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:88
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
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
subroutine dorgtr(UPLO, N, A, LDA, TAU, WORK, LWORK, INFO)
DORGTR
Definition: dorgtr.f:125
subroutine dsytrd(UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO)
DSYTRD
Definition: dsytrd.f:194
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 xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:53
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 dormtr(SIDE, UPLO, TRANS, M, N, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMTR
Definition: dormtr.f:173
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55
subroutine dstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
DSTEIN
Definition: dstein.f:176