LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ dsyevx()

subroutine dsyevx ( character  JOBZ,
character  RANGE,
character  UPLO,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision  VL,
double precision  VU,
integer  IL,
integer  IU,
double precision  ABSTOL,
integer  M,
double precision, dimension( * )  W,
double precision, dimension( ldz, * )  Z,
integer  LDZ,
double precision, dimension( * )  WORK,
integer  LWORK,
integer, dimension( * )  IWORK,
integer, dimension( * )  IFAIL,
integer  INFO 
)

DSYEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices

Download DSYEVX + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 DSYEVX computes selected eigenvalues and, optionally, eigenvectors
 of a real symmetric matrix A.  Eigenvalues and eigenvectors can be
 selected by specifying either a range of values or a range of indices
 for the desired eigenvalues.
Parameters
[in]JOBZ
          JOBZ is CHARACTER*1
          = 'N':  Compute eigenvalues only;
          = 'V':  Compute eigenvalues and eigenvectors.
[in]RANGE
          RANGE is CHARACTER*1
          = 'A': all eigenvalues will be found.
          = 'V': all eigenvalues in the half-open interval (VL,VU]
                 will be found.
          = 'I': the IL-th through IU-th eigenvalues will be found.
[in]UPLO
          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in,out]A
          A is DOUBLE PRECISION array, dimension (LDA, N)
          On entry, the symmetric matrix A.  If UPLO = 'U', the
          leading N-by-N upper triangular part of A contains the
          upper triangular part of the matrix A.  If UPLO = 'L',
          the leading N-by-N lower triangular part of A contains
          the lower triangular part of the matrix A.
          On exit, the lower triangle (if UPLO='L') or the upper
          triangle (if UPLO='U') of A, including the diagonal, is
          destroyed.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[in]VL
          VL is DOUBLE PRECISION
          If RANGE='V', the lower bound of the interval to
          be searched for eigenvalues. VL < VU.
          Not referenced if RANGE = 'A' or 'I'.
[in]VU
          VU is DOUBLE PRECISION
          If RANGE='V', the upper bound of the interval to
          be searched for eigenvalues. VL < VU.
          Not referenced if RANGE = 'A' or 'I'.
[in]IL
          IL is INTEGER
          If RANGE='I', the index of the
          smallest eigenvalue to be returned.
          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
          Not referenced if RANGE = 'A' or 'V'.
[in]IU
          IU is INTEGER
          If RANGE='I', the index of the
          largest eigenvalue to be returned.
          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
          Not referenced if RANGE = 'A' or 'V'.
[in]ABSTOL
          ABSTOL is DOUBLE PRECISION
          The absolute error tolerance for the eigenvalues.
          An approximate eigenvalue is accepted as converged
          when it is determined to lie in an interval [a,b]
          of width less than or equal to

                  ABSTOL + EPS *   max( |a|,|b| ) ,

          where EPS is the machine precision.  If ABSTOL is less than
          or equal to zero, then  EPS*|T|  will be used in its place,
          where |T| is the 1-norm of the tridiagonal matrix obtained
          by reducing A to tridiagonal form.

          Eigenvalues will be computed most accurately when ABSTOL is
          set to twice the underflow threshold 2*DLAMCH('S'), not zero.
          If this routine returns with INFO>0, indicating that some
          eigenvectors did not converge, try setting ABSTOL to
          2*DLAMCH('S').

          See "Computing Small Singular Values of Bidiagonal Matrices
          with Guaranteed High Relative Accuracy," by Demmel and
          Kahan, LAPACK Working Note #3.
[out]M
          M is INTEGER
          The total number of eigenvalues found.  0 <= M <= N.
          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
[out]W
          W is DOUBLE PRECISION array, dimension (N)
          On normal exit, the first M elements contain the selected
          eigenvalues in ascending order.
[out]Z
          Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M))
          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
          contain the orthonormal eigenvectors of the matrix A
          corresponding to the selected eigenvalues, with the i-th
          column of Z holding the eigenvector associated with W(i).
          If an eigenvector fails to converge, then that column of Z
          contains the latest approximation to the eigenvector, and the
          index of the eigenvector is returned in IFAIL.
          If JOBZ = 'N', then Z is not referenced.
          Note: the user must ensure that at least max(1,M) columns are
          supplied in the array Z; if RANGE = 'V', the exact value of M
          is not known in advance and an upper bound must be used.
[in]LDZ
          LDZ is INTEGER
          The leading dimension of the array Z.  LDZ >= 1, and if
          JOBZ = 'V', LDZ >= max(1,N).
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The length of the array WORK.  LWORK >= 1, when N <= 1;
          otherwise 8*N.
          For optimal efficiency, LWORK >= (NB+3)*N,
          where NB is the max of the blocksize for DSYTRD and DORMTR
          returned by ILAENV.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]IWORK
          IWORK is INTEGER array, dimension (5*N)
[out]IFAIL
          IFAIL is INTEGER array, dimension (N)
          If JOBZ = 'V', then if INFO = 0, the first M elements of
          IFAIL are zero.  If INFO > 0, then IFAIL contains the
          indices of the eigenvectors that failed to converge.
          If JOBZ = 'N', then IFAIL is not referenced.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i, then i eigenvectors failed to converge.
                Their indices are stored in array IFAIL.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 250 of file dsyevx.f.

253 *
254 * -- LAPACK driver routine --
255 * -- LAPACK is a software package provided by Univ. of Tennessee, --
256 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
257 *
258 * .. Scalar Arguments ..
259  CHARACTER JOBZ, RANGE, UPLO
260  INTEGER IL, INFO, IU, LDA, LDZ, LWORK, M, N
261  DOUBLE PRECISION ABSTOL, VL, VU
262 * ..
263 * .. Array Arguments ..
264  INTEGER IFAIL( * ), IWORK( * )
265  DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * ), Z( LDZ, * )
266 * ..
267 *
268 * =====================================================================
269 *
270 * .. Parameters ..
271  DOUBLE PRECISION ZERO, ONE
272  parameter( zero = 0.0d+0, one = 1.0d+0 )
273 * ..
274 * .. Local Scalars ..
275  LOGICAL ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
276  $ WANTZ
277  CHARACTER ORDER
278  INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
279  $ INDISP, INDIWO, INDTAU, INDWKN, INDWRK, ISCALE,
280  $ ITMP1, J, JJ, LLWORK, LLWRKN, LWKMIN,
281  $ LWKOPT, NB, NSPLIT
282  DOUBLE PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
283  $ SIGMA, SMLNUM, TMP1, VLL, VUU
284 * ..
285 * .. External Functions ..
286  LOGICAL LSAME
287  INTEGER ILAENV
288  DOUBLE PRECISION DLAMCH, DLANSY
289  EXTERNAL lsame, ilaenv, dlamch, dlansy
290 * ..
291 * .. External Subroutines ..
292  EXTERNAL dcopy, dlacpy, dorgtr, dormtr, dscal, dstebz,
294 * ..
295 * .. Intrinsic Functions ..
296  INTRINSIC max, min, sqrt
297 * ..
298 * .. Executable Statements ..
299 *
300 * Test the input parameters.
301 *
302  lower = lsame( uplo, 'L' )
303  wantz = lsame( jobz, 'V' )
304  alleig = lsame( range, 'A' )
305  valeig = lsame( range, 'V' )
306  indeig = lsame( range, 'I' )
307  lquery = ( lwork.EQ.-1 )
308 *
309  info = 0
310  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
311  info = -1
312  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
313  info = -2
314  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
315  info = -3
316  ELSE IF( n.LT.0 ) THEN
317  info = -4
318  ELSE IF( lda.LT.max( 1, n ) ) THEN
319  info = -6
320  ELSE
321  IF( valeig ) THEN
322  IF( n.GT.0 .AND. vu.LE.vl )
323  $ info = -8
324  ELSE IF( indeig ) THEN
325  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
326  info = -9
327  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
328  info = -10
329  END IF
330  END IF
331  END IF
332  IF( info.EQ.0 ) THEN
333  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
334  info = -15
335  END IF
336  END IF
337 *
338  IF( info.EQ.0 ) THEN
339  IF( n.LE.1 ) THEN
340  lwkmin = 1
341  work( 1 ) = lwkmin
342  ELSE
343  lwkmin = 8*n
344  nb = ilaenv( 1, 'DSYTRD', uplo, n, -1, -1, -1 )
345  nb = max( nb, ilaenv( 1, 'DORMTR', uplo, n, -1, -1, -1 ) )
346  lwkopt = max( lwkmin, ( nb + 3 )*n )
347  work( 1 ) = lwkopt
348  END IF
349 *
350  IF( lwork.LT.lwkmin .AND. .NOT.lquery )
351  $ info = -17
352  END IF
353 *
354  IF( info.NE.0 ) THEN
355  CALL xerbla( 'DSYEVX', -info )
356  RETURN
357  ELSE IF( lquery ) THEN
358  RETURN
359  END IF
360 *
361 * Quick return if possible
362 *
363  m = 0
364  IF( n.EQ.0 ) THEN
365  RETURN
366  END IF
367 *
368  IF( n.EQ.1 ) THEN
369  IF( alleig .OR. indeig ) THEN
370  m = 1
371  w( 1 ) = a( 1, 1 )
372  ELSE
373  IF( vl.LT.a( 1, 1 ) .AND. vu.GE.a( 1, 1 ) ) THEN
374  m = 1
375  w( 1 ) = a( 1, 1 )
376  END IF
377  END IF
378  IF( wantz )
379  $ z( 1, 1 ) = one
380  RETURN
381  END IF
382 *
383 * Get machine constants.
384 *
385  safmin = dlamch( 'Safe minimum' )
386  eps = dlamch( 'Precision' )
387  smlnum = safmin / eps
388  bignum = one / smlnum
389  rmin = sqrt( smlnum )
390  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
391 *
392 * Scale matrix to allowable range, if necessary.
393 *
394  iscale = 0
395  abstll = abstol
396  IF( valeig ) THEN
397  vll = vl
398  vuu = vu
399  END IF
400  anrm = dlansy( 'M', uplo, n, a, lda, work )
401  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
402  iscale = 1
403  sigma = rmin / anrm
404  ELSE IF( anrm.GT.rmax ) THEN
405  iscale = 1
406  sigma = rmax / anrm
407  END IF
408  IF( iscale.EQ.1 ) THEN
409  IF( lower ) THEN
410  DO 10 j = 1, n
411  CALL dscal( n-j+1, sigma, a( j, j ), 1 )
412  10 CONTINUE
413  ELSE
414  DO 20 j = 1, n
415  CALL dscal( j, sigma, a( 1, j ), 1 )
416  20 CONTINUE
417  END IF
418  IF( abstol.GT.0 )
419  $ abstll = abstol*sigma
420  IF( valeig ) THEN
421  vll = vl*sigma
422  vuu = vu*sigma
423  END IF
424  END IF
425 *
426 * Call DSYTRD to reduce symmetric matrix to tridiagonal form.
427 *
428  indtau = 1
429  inde = indtau + n
430  indd = inde + n
431  indwrk = indd + n
432  llwork = lwork - indwrk + 1
433  CALL dsytrd( uplo, n, a, lda, work( indd ), work( inde ),
434  $ work( indtau ), work( indwrk ), llwork, iinfo )
435 *
436 * If all eigenvalues are desired and ABSTOL is less than or equal to
437 * zero, then call DSTERF or DORGTR and SSTEQR. If this fails for
438 * some eigenvalue, then try DSTEBZ.
439 *
440  test = .false.
441  IF( indeig ) THEN
442  IF( il.EQ.1 .AND. iu.EQ.n ) THEN
443  test = .true.
444  END IF
445  END IF
446  IF( ( alleig .OR. test ) .AND. ( abstol.LE.zero ) ) THEN
447  CALL dcopy( n, work( indd ), 1, w, 1 )
448  indee = indwrk + 2*n
449  IF( .NOT.wantz ) THEN
450  CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
451  CALL dsterf( n, w, work( indee ), info )
452  ELSE
453  CALL dlacpy( 'A', n, n, a, lda, z, ldz )
454  CALL dorgtr( uplo, n, z, ldz, work( indtau ),
455  $ work( indwrk ), llwork, iinfo )
456  CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
457  CALL dsteqr( jobz, n, w, work( indee ), z, ldz,
458  $ work( indwrk ), info )
459  IF( info.EQ.0 ) THEN
460  DO 30 i = 1, n
461  ifail( i ) = 0
462  30 CONTINUE
463  END IF
464  END IF
465  IF( info.EQ.0 ) THEN
466  m = n
467  GO TO 40
468  END IF
469  info = 0
470  END IF
471 *
472 * Otherwise, call DSTEBZ and, if eigenvectors are desired, SSTEIN.
473 *
474  IF( wantz ) THEN
475  order = 'B'
476  ELSE
477  order = 'E'
478  END IF
479  indibl = 1
480  indisp = indibl + n
481  indiwo = indisp + n
482  CALL dstebz( range, order, n, vll, vuu, il, iu, abstll,
483  $ work( indd ), work( inde ), m, nsplit, w,
484  $ iwork( indibl ), iwork( indisp ), work( indwrk ),
485  $ iwork( indiwo ), info )
486 *
487  IF( wantz ) THEN
488  CALL dstein( n, work( indd ), work( inde ), m, w,
489  $ iwork( indibl ), iwork( indisp ), z, ldz,
490  $ work( indwrk ), iwork( indiwo ), ifail, info )
491 *
492 * Apply orthogonal matrix used in reduction to tridiagonal
493 * form to eigenvectors returned by DSTEIN.
494 *
495  indwkn = inde
496  llwrkn = lwork - indwkn + 1
497  CALL dormtr( 'L', uplo, 'N', n, m, a, lda, work( indtau ), z,
498  $ ldz, work( indwkn ), llwrkn, iinfo )
499  END IF
500 *
501 * If matrix was scaled, then rescale eigenvalues appropriately.
502 *
503  40 CONTINUE
504  IF( iscale.EQ.1 ) THEN
505  IF( info.EQ.0 ) THEN
506  imax = m
507  ELSE
508  imax = info - 1
509  END IF
510  CALL dscal( imax, one / sigma, w, 1 )
511  END IF
512 *
513 * If eigenvalues are not in order, then sort them, along with
514 * eigenvectors.
515 *
516  IF( wantz ) THEN
517  DO 60 j = 1, m - 1
518  i = 0
519  tmp1 = w( j )
520  DO 50 jj = j + 1, m
521  IF( w( jj ).LT.tmp1 ) THEN
522  i = jj
523  tmp1 = w( jj )
524  END IF
525  50 CONTINUE
526 *
527  IF( i.NE.0 ) THEN
528  itmp1 = iwork( indibl+i-1 )
529  w( i ) = w( j )
530  iwork( indibl+i-1 ) = iwork( indibl+j-1 )
531  w( j ) = tmp1
532  iwork( indibl+j-1 ) = itmp1
533  CALL dswap( n, z( 1, i ), 1, z( 1, j ), 1 )
534  IF( info.NE.0 ) THEN
535  itmp1 = ifail( i )
536  ifail( i ) = ifail( j )
537  ifail( j ) = itmp1
538  END IF
539  END IF
540  60 CONTINUE
541  END IF
542 *
543 * Set WORK(1) to optimal workspace size.
544 *
545  work( 1 ) = lwkopt
546 *
547  RETURN
548 *
549 * End of DSYEVX
550 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.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:273
subroutine dsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
DSTEQR
Definition: dsteqr.f:131
subroutine dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:86
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:82
subroutine dstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
DSTEIN
Definition: dstein.f:174
subroutine dormtr(SIDE, UPLO, TRANS, M, N, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMTR
Definition: dormtr.f:171
subroutine dorgtr(UPLO, N, A, LDA, TAU, WORK, LWORK, INFO)
DORGTR
Definition: dorgtr.f:123
double precision function dlansy(NORM, UPLO, N, A, LDA, WORK)
DLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition: dlansy.f:122
subroutine dsytrd(UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO)
DSYTRD
Definition: dsytrd.f:192
Here is the call graph for this function:
Here is the caller graph for this function: