LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ dsyevx_2stage()

subroutine dsyevx_2stage ( 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_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices

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

Purpose:
 DSYEVX_2STAGE computes selected eigenvalues and, optionally, eigenvectors
 of a real symmetric matrix A using the 2stage technique for
 the reduction to tridiagonal.  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.
                  Not available in this release.
[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  
          If JOBZ = 'N' and N > 1, LWORK must be queried.
                                   LWORK = MAX(1, 8*N, dimension) where
                                   dimension = max(stage1,stage2) + (KD+1)*N + 3*N
                                             = N*KD + N*max(KD+1,FACTOPTNB) 
                                               + max(2*KD*KD, KD*NTHREADS) 
                                               + (KD+1)*N + 3*N
                                   where KD is the blocking size of the reduction,
                                   FACTOPTNB is the blocking used by the QR or LQ
                                   algorithm, usually FACTOPTNB=128 is a good choice
                                   NTHREADS is the number of threads used when
                                   openMP compilation is enabled, otherwise =1.
          If JOBZ = 'V' and N > 1, LWORK must be queried. Not yet available

          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
Further Details:
  All details about the 2stage techniques are available in:

  Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
  Parallel reduction to condensed forms for symmetric eigenvalue problems
  using aggregated fine-grained and memory-aware kernels. In Proceedings
  of 2011 International Conference for High Performance Computing,
  Networking, Storage and Analysis (SC '11), New York, NY, USA,
  Article 8 , 11 pages.
  http://doi.acm.org/10.1145/2063384.2063394

  A. Haidar, J. Kurzak, P. Luszczek, 2013.
  An improved parallel singular value algorithm and its implementation 
  for multicore hardware, In Proceedings of 2013 International Conference
  for High Performance Computing, Networking, Storage and Analysis (SC '13).
  Denver, Colorado, USA, 2013.
  Article 90, 12 pages.
  http://doi.acm.org/10.1145/2503210.2503292

  A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
  A novel hybrid CPU-GPU generalized eigensolver for electronic structure 
  calculations based on fine-grained memory aware tasks.
  International Journal of High Performance Computing Applications.
  Volume 28 Issue 2, Pages 196-209, May 2014.
  http://hpc.sagepub.com/content/28/2/196 

Definition at line 302 of file dsyevx_2stage.f.

302 *
303  IMPLICIT NONE
304 *
305 * -- LAPACK driver routine (version 3.8.0) --
306 * -- LAPACK is a software package provided by Univ. of Tennessee, --
307 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
308 * June 2016
309 *
310 * .. Scalar Arguments ..
311  CHARACTER JOBZ, RANGE, UPLO
312  INTEGER IL, INFO, IU, LDA, LDZ, LWORK, M, N
313  DOUBLE PRECISION ABSTOL, VL, VU
314 * ..
315 * .. Array Arguments ..
316  INTEGER IFAIL( * ), IWORK( * )
317  DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * ), Z( LDZ, * )
318 * ..
319 *
320 * =====================================================================
321 *
322 * .. Parameters ..
323  DOUBLE PRECISION ZERO, ONE
324  parameter( zero = 0.0d+0, one = 1.0d+0 )
325 * ..
326 * .. Local Scalars ..
327  LOGICAL ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
328  $ WANTZ
329  CHARACTER ORDER
330  INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
331  $ INDISP, INDIWO, INDTAU, INDWKN, INDWRK, ISCALE,
332  $ ITMP1, J, JJ, LLWORK, LLWRKN,
333  $ NSPLIT, LWMIN, LHTRD, LWTRD, KD, IB, INDHOUS
334  DOUBLE PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
335  $ SIGMA, SMLNUM, TMP1, VLL, VUU
336 * ..
337 * .. External Functions ..
338  LOGICAL LSAME
339  INTEGER ILAENV2STAGE
340  DOUBLE PRECISION DLAMCH, DLANSY
341  EXTERNAL lsame, dlamch, dlansy, ilaenv2stage
342 * ..
343 * .. External Subroutines ..
344  EXTERNAL dcopy, dlacpy, dorgtr, dormtr, dscal, dstebz,
346  $ dsytrd_2stage
347 * ..
348 * .. Intrinsic Functions ..
349  INTRINSIC max, min, sqrt
350 * ..
351 * .. Executable Statements ..
352 *
353 * Test the input parameters.
354 *
355  lower = lsame( uplo, 'L' )
356  wantz = lsame( jobz, 'V' )
357  alleig = lsame( range, 'A' )
358  valeig = lsame( range, 'V' )
359  indeig = lsame( range, 'I' )
360  lquery = ( lwork.EQ.-1 )
361 *
362  info = 0
363  IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
364  info = -1
365  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
366  info = -2
367  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
368  info = -3
369  ELSE IF( n.LT.0 ) THEN
370  info = -4
371  ELSE IF( lda.LT.max( 1, n ) ) THEN
372  info = -6
373  ELSE
374  IF( valeig ) THEN
375  IF( n.GT.0 .AND. vu.LE.vl )
376  $ info = -8
377  ELSE IF( indeig ) THEN
378  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
379  info = -9
380  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
381  info = -10
382  END IF
383  END IF
384  END IF
385  IF( info.EQ.0 ) THEN
386  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
387  info = -15
388  END IF
389  END IF
390 *
391  IF( info.EQ.0 ) THEN
392  IF( n.LE.1 ) THEN
393  lwmin = 1
394  work( 1 ) = lwmin
395  ELSE
396  kd = ilaenv2stage( 1, 'DSYTRD_2STAGE', jobz,
397  $ n, -1, -1, -1 )
398  ib = ilaenv2stage( 2, 'DSYTRD_2STAGE', jobz,
399  $ n, kd, -1, -1 )
400  lhtrd = ilaenv2stage( 3, 'DSYTRD_2STAGE', jobz,
401  $ n, kd, ib, -1 )
402  lwtrd = ilaenv2stage( 4, 'DSYTRD_2STAGE', jobz,
403  $ n, kd, ib, -1 )
404  lwmin = max( 8*n, 3*n + lhtrd + lwtrd )
405  work( 1 ) = lwmin
406  END IF
407 *
408  IF( lwork.LT.lwmin .AND. .NOT.lquery )
409  $ info = -17
410  END IF
411 *
412  IF( info.NE.0 ) THEN
413  CALL xerbla( 'DSYEVX_2STAGE', -info )
414  RETURN
415  ELSE IF( lquery ) THEN
416  RETURN
417  END IF
418 *
419 * Quick return if possible
420 *
421  m = 0
422  IF( n.EQ.0 ) THEN
423  RETURN
424  END IF
425 *
426  IF( n.EQ.1 ) THEN
427  IF( alleig .OR. indeig ) THEN
428  m = 1
429  w( 1 ) = a( 1, 1 )
430  ELSE
431  IF( vl.LT.a( 1, 1 ) .AND. vu.GE.a( 1, 1 ) ) THEN
432  m = 1
433  w( 1 ) = a( 1, 1 )
434  END IF
435  END IF
436  IF( wantz )
437  $ z( 1, 1 ) = one
438  RETURN
439  END IF
440 *
441 * Get machine constants.
442 *
443  safmin = dlamch( 'Safe minimum' )
444  eps = dlamch( 'Precision' )
445  smlnum = safmin / eps
446  bignum = one / smlnum
447  rmin = sqrt( smlnum )
448  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
449 *
450 * Scale matrix to allowable range, if necessary.
451 *
452  iscale = 0
453  abstll = abstol
454  IF( valeig ) THEN
455  vll = vl
456  vuu = vu
457  END IF
458  anrm = dlansy( 'M', uplo, n, a, lda, work )
459  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
460  iscale = 1
461  sigma = rmin / anrm
462  ELSE IF( anrm.GT.rmax ) THEN
463  iscale = 1
464  sigma = rmax / anrm
465  END IF
466  IF( iscale.EQ.1 ) THEN
467  IF( lower ) THEN
468  DO 10 j = 1, n
469  CALL dscal( n-j+1, sigma, a( j, j ), 1 )
470  10 CONTINUE
471  ELSE
472  DO 20 j = 1, n
473  CALL dscal( j, sigma, a( 1, j ), 1 )
474  20 CONTINUE
475  END IF
476  IF( abstol.GT.0 )
477  $ abstll = abstol*sigma
478  IF( valeig ) THEN
479  vll = vl*sigma
480  vuu = vu*sigma
481  END IF
482  END IF
483 *
484 * Call DSYTRD_2STAGE to reduce symmetric matrix to tridiagonal form.
485 *
486  indtau = 1
487  inde = indtau + n
488  indd = inde + n
489  indhous = indd + n
490  indwrk = indhous + lhtrd
491  llwork = lwork - indwrk + 1
492 *
493  CALL dsytrd_2stage( jobz, uplo, n, a, lda, work( indd ),
494  $ work( inde ), work( indtau ), work( indhous ),
495  $ lhtrd, work( indwrk ), llwork, iinfo )
496 *
497 * If all eigenvalues are desired and ABSTOL is less than or equal to
498 * zero, then call DSTERF or DORGTR and SSTEQR. If this fails for
499 * some eigenvalue, then try DSTEBZ.
500 *
501  test = .false.
502  IF( indeig ) THEN
503  IF( il.EQ.1 .AND. iu.EQ.n ) THEN
504  test = .true.
505  END IF
506  END IF
507  IF( ( alleig .OR. test ) .AND. ( abstol.LE.zero ) ) THEN
508  CALL dcopy( n, work( indd ), 1, w, 1 )
509  indee = indwrk + 2*n
510  IF( .NOT.wantz ) THEN
511  CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
512  CALL dsterf( n, w, work( indee ), info )
513  ELSE
514  CALL dlacpy( 'A', n, n, a, lda, z, ldz )
515  CALL dorgtr( uplo, n, z, ldz, work( indtau ),
516  $ work( indwrk ), llwork, iinfo )
517  CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
518  CALL dsteqr( jobz, n, w, work( indee ), z, ldz,
519  $ work( indwrk ), info )
520  IF( info.EQ.0 ) THEN
521  DO 30 i = 1, n
522  ifail( i ) = 0
523  30 CONTINUE
524  END IF
525  END IF
526  IF( info.EQ.0 ) THEN
527  m = n
528  GO TO 40
529  END IF
530  info = 0
531  END IF
532 *
533 * Otherwise, call DSTEBZ and, if eigenvectors are desired, SSTEIN.
534 *
535  IF( wantz ) THEN
536  order = 'B'
537  ELSE
538  order = 'E'
539  END IF
540  indibl = 1
541  indisp = indibl + n
542  indiwo = indisp + n
543  CALL dstebz( range, order, n, vll, vuu, il, iu, abstll,
544  $ work( indd ), work( inde ), m, nsplit, w,
545  $ iwork( indibl ), iwork( indisp ), work( indwrk ),
546  $ iwork( indiwo ), info )
547 *
548  IF( wantz ) THEN
549  CALL dstein( n, work( indd ), work( inde ), m, w,
550  $ iwork( indibl ), iwork( indisp ), z, ldz,
551  $ work( indwrk ), iwork( indiwo ), ifail, info )
552 *
553 * Apply orthogonal matrix used in reduction to tridiagonal
554 * form to eigenvectors returned by DSTEIN.
555 *
556  indwkn = inde
557  llwrkn = lwork - indwkn + 1
558  CALL dormtr( 'L', uplo, 'N', n, m, a, lda, work( indtau ), z,
559  $ ldz, work( indwkn ), llwrkn, iinfo )
560  END IF
561 *
562 * If matrix was scaled, then rescale eigenvalues appropriately.
563 *
564  40 CONTINUE
565  IF( iscale.EQ.1 ) THEN
566  IF( info.EQ.0 ) THEN
567  imax = m
568  ELSE
569  imax = info - 1
570  END IF
571  CALL dscal( imax, one / sigma, w, 1 )
572  END IF
573 *
574 * If eigenvalues are not in order, then sort them, along with
575 * eigenvectors.
576 *
577  IF( wantz ) THEN
578  DO 60 j = 1, m - 1
579  i = 0
580  tmp1 = w( j )
581  DO 50 jj = j + 1, m
582  IF( w( jj ).LT.tmp1 ) THEN
583  i = jj
584  tmp1 = w( jj )
585  END IF
586  50 CONTINUE
587 *
588  IF( i.NE.0 ) THEN
589  itmp1 = iwork( indibl+i-1 )
590  w( i ) = w( j )
591  iwork( indibl+i-1 ) = iwork( indibl+j-1 )
592  w( j ) = tmp1
593  iwork( indibl+j-1 ) = itmp1
594  CALL dswap( n, z( 1, i ), 1, z( 1, j ), 1 )
595  IF( info.NE.0 ) THEN
596  itmp1 = ifail( i )
597  ifail( i ) = ifail( j )
598  ifail( j ) = itmp1
599  END IF
600  END IF
601  60 CONTINUE
602  END IF
603 *
604 * Set WORK(1) to optimal workspace size.
605 *
606  work( 1 ) = lwmin
607 *
608  RETURN
609 *
610 * End of DSYEVX_2STAGE
611 *
Here is the call graph for this function:
Here is the caller graph for this function:
dstein
subroutine dstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
DSTEIN
Definition: dstein.f:176
dstebz
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
dlansy
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:124
dsterf
subroutine dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:88
dsytrd_2stage
subroutine dsytrd_2stage(VECT, UPLO, N, A, LDA, D, E, TAU, HOUS2, LHOUS2, WORK, LWORK, INFO)
DSYTRD_2STAGE
Definition: dsytrd_2stage.f:226
dcopy
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:84
dsteqr
subroutine dsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
DSTEQR
Definition: dsteqr.f:133
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
dswap
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:84
dormtr
subroutine dormtr(SIDE, UPLO, TRANS, M, N, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMTR
Definition: dormtr.f:173
ilaenv2stage
integer function ilaenv2stage(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV2STAGE
Definition: ilaenv2stage.f:151
dscal
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:81
dorgtr
subroutine dorgtr(UPLO, N, A, LDA, TAU, WORK, LWORK, INFO)
DORGTR
Definition: dorgtr.f:125
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:70
dlacpy
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