dstegr2a.f

Go to the documentation of this file.

      SUBROUTINE dstegr2a( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
     $                   M, W, Z, LDZ, NZC, WORK, LWORK, IWORK,
     $                   LIWORK, DOL, DOU, NEEDIL, NEEDIU,
     $                   INDERR, NSPLIT, PIVMIN, SCALE, WL, WU,
     $                   INFO )
*
*  -- ScaLAPACK auxiliary routine (version 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver
*     July 4, 2010
*
      IMPLICIT NONE
*
*     .. Scalar Arguments ..
      CHARACTER          JOBZ, RANGE
      INTEGER            DOL, DOU, IL, INDERR, INFO, IU, LDZ, LIWORK,
     $                   LWORK, M, N, NEEDIL, NEEDIU, NSPLIT, NZC
      DOUBLE PRECISION PIVMIN, SCALE, VL, VU, WL, WU
 
*     ..
*     .. Array Arguments ..
      INTEGER            IWORK( * )
      DOUBLE PRECISION   D( * ), E( * ), W( * ), WORK( * )
      DOUBLE PRECISION   Z( LDZ, * )
*     ..
*
*  Purpose
*  =======
*
*  DSTEGR2A computes selected eigenvalues and initial representations.
*  needed for eigenvector computations in DSTEGR2B. It is invoked in the 
*  ScaLAPACK MRRR driver PDSYEVR and the corresponding Hermitian
*  version when both eigenvalues and eigenvectors are computed in parallel.
*  on multiple processors. For this case, DSTEGR2A implements the FIRST 
*  part of the MRRR algorithm, parallel eigenvalue computation and finding
*  the root RRR. At the end of DSTEGR2A,
*  other processors might have a part of the spectrum that is needed to
*  continue the computation locally. Once this eigenvalue information has
*  been received by the processor, the computation can then proceed by calling 
*  the SECOND part of the parallel MRRR algorithm, DSTEGR2B.
*
*  Please note:
*  1. The calling sequence has two additional INTEGER parameters, 
*     (compared to LAPACK's DSTEGR), these are
*     DOL and DOU and should satisfy M>=DOU>=DOL>=1. 
*     These parameters are only relevant for the case JOBZ = 'V'.
*
*     Globally invoked over all processors, DSTEGR2A computes 
*     ALL the eigenVALUES specified by RANGE. 
*     RANGE= 'A': all eigenvalues will be found.
*          = 'V': all eigenvalues in (VL,VU] will be found.
*          = 'I': the IL-th through IU-th eigenvalues will be found.
*
*     DSTEGR2A LOCALLY only computes the eigenvalues 
*     corresponding to eigenvalues DOL through DOU in W. (That is,
*     instead of computing the eigenvectors belonging to W(1) 
*     through W(M), only the eigenvectors belonging to eigenvalues
*     W(DOL) through W(DOU) are computed. In this case, only the
*     eigenvalues DOL:DOU are guaranteed to be fully accurate.
*
*  2. M is NOT the number of eigenvalues specified by RANGE, but it is 
*     M = DOU - DOL + 1. Instead, M refers to the number of eigenvalues computed on 
*     this processor.
*
*  3. While no eigenvectors are computed in DSTEGR2A itself (this is
*     done later in DSTEGR2B), the interface
*     If JOBZ = 'V' then, depending on RANGE and DOL, DOU, DSTEGR2A 
*     might need more workspace in Z then the original DSTEGR. 
*     In particular, the arrays W and Z might not contain all the wanted eigenpairs
*     locally, instead this information is distributed over other 
*     processors.
*  
*  Arguments
*  =========
*
*  JOBZ    (input) CHARACTER*1
*          = 'N':  Compute eigenvalues only;
*          = 'V':  Compute eigenvalues and eigenvectors.
*
*  RANGE   (input) 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.
*
*  N       (input) INTEGER
*          The order of the matrix.  N >= 0.
*
*  D       (input/output) DOUBLE PRECISION array, dimension (N)
*          On entry, the N diagonal elements of the tridiagonal matrix
*          T. On exit, D is overwritten.
*
*  E       (input/output) DOUBLE PRECISION array, dimension (N)
*          On entry, the (N-1) subdiagonal elements of the tridiagonal
*          matrix T in elements 1 to N-1 of E. E(N) need not be set on
*          input, but is used internally as workspace.
*          On exit, E is overwritten.
*
*  VL      (input) DOUBLE PRECISION
*  VU      (input) DOUBLE PRECISION
*          If RANGE='V', the lower and upper bounds of the interval to
*          be searched for eigenvalues. VL < VU.
*          Not referenced if RANGE = 'A' or 'I'.
*
*  IL      (input) INTEGER
*  IU      (input) INTEGER
*          If RANGE='I', the indices (in ascending order) of the
*          smallest and largest eigenvalues to be returned.
*          1 <= IL <= IU <= N, if N > 0.
*          Not referenced if RANGE = 'A' or 'V'.
*
*  M       (output) INTEGER
*          Globally summed over all processors, M equals 
*          the total number of eigenvalues found.  0 <= M <= N.
*          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
*          The local output equals M = DOU - DOL + 1.
*
*  W       (output) DOUBLE PRECISION array, dimension (N)
*          The first M elements contain approximations to the selected 
*          eigenvalues in ascending order. Note that immediately after 
*          exiting this routine, only the eigenvalues from
*          position DOL:DOU are to reliable on this processor
*          because the eigenvalue computation is done in parallel.          
*          The other entries outside DOL:DOU are very crude preliminary
*          approximations. Other processors hold reliable information on 
*          these other parts of the W array. 
*          This information is communicated in the ScaLAPACK driver.
*
*  Z       (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M) )
*          DSTEGR2A does not compute eigenvectors, this is done 
*          in DSTEGR2B. The argument Z as well as all related
*          other arguments only appear to keep the interface consistent
*          and to signal to the user that this subroutine is meant to 
*          be used when eigenvectors are computed.
*
*  LDZ     (input) INTEGER
*          The leading dimension of the array Z.  LDZ >= 1, and if
*          JOBZ = 'V', then LDZ >= max(1,N).
*
*  NZC     (input) INTEGER
*          The number of eigenvectors to be held in the array Z.  
*          If RANGE = 'A', then NZC >= max(1,N).
*          If RANGE = 'V', then NZC >= the number of eigenvalues in (VL,VU].
*          If RANGE = 'I', then NZC >= IU-IL+1.
*          If NZC = -1, then a workspace query is assumed; the
*          routine calculates the number of columns of the array Z that
*          are needed to hold the eigenvectors. 
*          This value is returned as the first entry of the Z array, and
*          no error message related to NZC is issued.
*
*  WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
*          On exit, if INFO = 0, WORK(1) returns the optimal
*          (and minimal) LWORK.
*
*  LWORK   (input) INTEGER
*          The dimension of the array WORK. LWORK >= max(1,18*N)
*          if JOBZ = 'V', and LWORK >= max(1,12*N) if JOBZ = 'N'.
*          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.
*
*  IWORK   (workspace/output) INTEGER array, dimension (LIWORK)
*          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
*
*  LIWORK  (input) INTEGER
*          The dimension of the array IWORK.  LIWORK >= max(1,10*N)
*          if the eigenvectors are desired, and LIWORK >= max(1,8*N)
*          if only the eigenvalues are to be computed.
*          If LIWORK = -1, then a workspace query is assumed; the
*          routine only calculates the optimal size of the IWORK array,
*          returns this value as the first entry of the IWORK array, and
*          no error message related to LIWORK is issued.
*
*  DOL     (input) INTEGER
*  DOU     (input) INTEGER
*          From all the eigenvalues W(1:M), only eigenvalues
*          W(DOL:DOU) are computed.
*
*  NEEDIL  (output) INTEGER
*  NEEDIU  (output) INTEGER
*          The indices of the leftmost and rightmost eigenvalues
*          needed to accurately compute the relevant part of the 
*          representation tree. This information can be used to 
*          find out which processors have the relevant eigenvalue
*          information needed so that it can be communicated.
*
*  INDERR  (output) INTEGER
*          INDERR points to the place in the work space where 
*          the eigenvalue uncertainties (errors) are stored.
*
*  NSPLIT  (output) INTEGER
*          The number of blocks T splits into. 1 <= NSPLIT <= N.
*
*  PIVMIN  (output) DOUBLE PRECISION
*          The minimum pivot in the sturm sequence for T.
*
*  SCALE   (output) DOUBLE PRECISION 
*          The scaling factor for the tridiagonal T.
*
*  WL      (output) DOUBLE PRECISION
*  WU      (output) DOUBLE PRECISION
*          The interval (WL, WU] contains all the wanted eigenvalues.         
*          It is either given by the user or computed in DLARRE2A.
*
*  INFO    (output) INTEGER
*          On exit, INFO
*          = 0:  successful exit
*          other:if INFO = -i, the i-th argument had an illegal value
*                if INFO = 10X, internal error in DLARRE2A,
*                Here, the digit X = ABS( IINFO ) < 10, where IINFO is 
*                the nonzero error code returned by DLARRE2A.
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE, FOUR, MINRGP
      PARAMETER          ( ZERO = 0.0d0, one = 1.0d0,
     $                     four = 4.0d0,
     $                     minrgp = 1.0d-3 )
*     ..
*     .. Local Scalars ..
      LOGICAL            ALLEIG, INDEIG, LQUERY, VALEIG, WANTZ, ZQUERY
      INTEGER            IIL, IINDBL, IINDW, IINDWK, IINFO, IINSPL, IIU,
     $                   INDE2, INDGP, INDGRS, INDSDM, INDWRK, ITMP,
     $                   ITMP2, J, LIWMIN, LWMIN, NZCMIN
      DOUBLE PRECISION   BIGNUM, EPS, RMAX, RMIN, RTOL1, RTOL2, SAFMIN,
     $                   smlnum, thresh, tnrm
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      DOUBLE PRECISION   DLAMCH, DLANST
      EXTERNAL           LSAME, DLAMCH, DLANST
*     ..
*     .. External Subroutines ..
      EXTERNAL           dlarrc, dlarre2a, dscal
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          dble, max, min, sqrt
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      wantz = lsame( jobz, 'V' )
      alleig = lsame( range, 'A' )
      valeig = lsame( range, 'V' )
      indeig = lsame( range, 'I' )
*
      lquery = ( ( lwork.EQ.-1 ).OR.( liwork.EQ.-1 ) )
      zquery = ( nzc.EQ.-1 )
 
*     DSTEGR2A needs WORK of size 6*N, IWORK of size 3*N.
*     In addition, DLARRE2A needs WORK of size 6*N, IWORK of size 5*N.
*     Furthermore, DLARRV2 needs WORK of size 12*N, IWORK of size 7*N.
*     Workspace is kept consistent with DSTEGR2B even though 
*     DLARRV2 is not called here.
      IF( wantz ) THEN
         lwmin = 18*n
         liwmin = 10*n
      ELSE
*        need less workspace if only the eigenvalues are wanted         
         lwmin = 12*n
         liwmin = 8*n
      ENDIF
 
      wl = zero
      wu = zero
      iil = 0
      iiu = 0
 
      IF( valeig ) THEN
*        We do not reference VL, VU in the cases RANGE = 'I','A'
*        The interval (WL, WU] contains all the wanted eigenvalues.         
*        It is either given by the user or computed in DLARRE2A.
         wl = vl
         wu = vu
      ELSEIF( indeig ) THEN
*        We do not reference IL, IU in the cases RANGE = 'V','A'
         iil = il
         iiu = iu
      ENDIF  
*
      info = 0
      IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
         info = -1
      ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
         info = -2
      ELSE IF( n.LT.0 ) THEN
         info = -3
      ELSE IF( valeig .AND. n.GT.0 .AND. wu.LE.wl ) THEN
         info = -7
      ELSE IF( indeig .AND. ( iil.LT.1 .OR. iil.GT.n ) ) THEN
         info = -8
      ELSE IF( indeig .AND. ( iiu.LT.iil .OR. iiu.GT.n ) ) THEN
         info = -9
      ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
         info = -13
      ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
         info = -17
      ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
         info = -19
      END IF
*
*     Get machine constants.
*
      safmin = dlamch( 'safe minimum' )
      EPS = DLAMCH( 'precision' )
      SMLNUM = SAFMIN / EPS
      BIGNUM = ONE / SMLNUM
      RMIN = SQRT( SMLNUM )
      RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) )
*
.EQ.      IF( INFO0 ) THEN
         WORK( 1 ) = LWMIN
         IWORK( 1 ) = LIWMIN
*
.AND.         IF( WANTZ  ALLEIG ) THEN
            NZCMIN = N
            IIL = 1
            IIU = N
.AND.         ELSE IF( WANTZ  VALEIG ) THEN
            CALL DLARRC( 't', N, VL, VU, D, E, SAFMIN, 
     $                            NZCMIN, ITMP, ITMP2, INFO )
            IIL = ITMP+1
            IIU = ITMP2
.AND.         ELSE IF( WANTZ  INDEIG ) THEN
            NZCMIN = IIU-IIL+1
         ELSE 
*           WANTZ .EQ. FALSE.   
            NZCMIN = 0
         ENDIF  
.AND..EQ.         IF( ZQUERY  INFO0 ) THEN
            Z( 1,1 ) = NZCMIN
.LT..AND..NOT.         ELSE IF( NZCNZCMIN  ZQUERY ) THEN
            INFO = -14
         END IF
      END IF
 
      IF ( WANTZ ) THEN
.LT..OR..GT.         IF ( DOL1  DOLNZCMIN ) THEN 
            INFO = -20
         ENDIF
.LT..OR..GT..OR..LT.         IF ( DOU1  DOUNZCMIN  DOUDOL) THEN 
            INFO = -21
         ENDIF
      ENDIF
 
.NE.      IF( INFO0 ) THEN
*
C         Disable sequential error handler
C         for parallel case
C         CALL XERBLA( 'DSTEGR2A', -INFO )
*
         RETURN
.OR.      ELSE IF( LQUERY  ZQUERY ) THEN
         RETURN
      END IF
 
*     Initialize NEEDIL and NEEDIU, these values are changed in DLARRE2A
      NEEDIL = DOU
      NEEDIU = DOL
*
*     Quick return if possible
*
      M = 0
.EQ.      IF( N0 )
     $   RETURN
*
.EQ.      IF( N1 ) THEN
.OR.         IF( ALLEIG  INDEIG ) THEN
            M = 1
            W( 1 ) = D( 1 )
         ELSE
.LT..AND..GE.            IF( WLD( 1 )  WUD( 1 ) ) THEN
               M = 1
               W( 1 ) = D( 1 )
            END IF
         END IF
         IF( WANTZ )
     $      Z( 1, 1 ) = ONE
         RETURN
      END IF
*
      INDGRS = 1
      INDERR = 2*N + 1
      INDGP = 3*N + 1
      INDSDM = 4*N + 1
      INDE2 = 5*N + 1
      INDWRK = 6*N + 1
*
      IINSPL = 1
      IINDBL = N + 1
      IINDW = 2*N + 1
      IINDWK = 3*N + 1
*
*     Scale matrix to allowable range, if necessary.
*
      SCALE = ONE
      TNRM = DLANST( 'm', N, D, E )
.GT..AND..LT.      IF( TNRMZERO  TNRMRMIN ) THEN
         SCALE = RMIN / TNRM
.GT.      ELSE IF( TNRMRMAX ) THEN
         SCALE = RMAX / TNRM
      END IF
.NE.      IF( SCALEONE ) THEN
         CALL DSCAL( N, SCALE, D, 1 )
         CALL DSCAL( N-1, SCALE, E, 1 )
         TNRM = TNRM*SCALE
         IF( VALEIG ) THEN
*           If eigenvalues in interval have to be found, 
*           scale (WL, WU] accordingly
            WL = WL*SCALE
            WU = WU*SCALE
         ENDIF
      END IF
*
*     Compute the desired eigenvalues of the tridiagonal after splitting
*     into smaller subblocks if the corresponding off-diagonal elements
*     are small
*     THRESH is the splitting parameter for DLARRA in DLARRE2A      
*     A negative THRESH forces the old splitting criterion based on the
*     size of the off-diagonal.
      THRESH = -EPS
      IINFO = 0
 
*     Store the squares of the offdiagonal values of T
      DO 5 J = 1, N-1
         WORK( INDE2+J-1 ) = E(J)**2
 5    CONTINUE
 
*     Set the tolerance parameters for bisection
.NOT.      IF( WANTZ ) THEN
*        DLARRE2A computes the eigenvalues to full precision.   
         RTOL1 = FOUR * EPS
         RTOL2 = FOUR * EPS
      ELSE   
*        DLARRE2A computes the eigenvalues to less than full precision.
*        DLARRV2 will refine the eigenvalue approximations, and we can
*        need less accurate initial bisection in DLARRE2A.
         RTOL1 = FOUR*SQRT(EPS)
         RTOL2 = MAX( SQRT(EPS)*5.0D-3, FOUR * EPS )
      ENDIF
      CALL DLARRE2A( RANGE, N, WL, WU, IIL, IIU, D, E, 
     $             WORK(INDE2), RTOL1, RTOL2, THRESH, NSPLIT, 
     $             IWORK( IINSPL ), M, DOL, DOU, NEEDIL, NEEDIU,
     $             W, WORK( INDERR ),
     $             WORK( INDGP ), IWORK( IINDBL ),
     $             IWORK( IINDW ), WORK( INDGRS ), 
     $             WORK( INDSDM ), PIVMIN,
     $             WORK( INDWRK ), IWORK( IINDWK ), 
     $             MINRGP, IINFO )
.NE.      IF( IINFO0 ) THEN
         INFO = 100 + ABS( IINFO )
         RETURN
      END IF
*     Note that if RANGE .NE. 'V', DLARRE2A computes bounds on the desired
*     part of the spectrum. All desired eigenvalues are contained in
*     (WL,WU]
 
 
      RETURN
*
*     End of DSTEGR2A
*
      END