LAPACK  3.8.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.
Date
June 2016

Definition at line 255 of file dsyevx.f.

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 *
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
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine dsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
DSTEQR
Definition: dsteqr.f:133
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, or the element of largest absolute value of a real symmetric matrix.
Definition: dlansy.f:124
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
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:84
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
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 dstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
DSTEBZ
Definition: dstebz.f:275
Here is the call graph for this function:
Here is the caller graph for this function: