LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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