dstegr2.f

Go to the documentation of this file.

      SUBROUTINE dstegr2( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
     $                   M, W, Z, LDZ, NZC, ISUPPZ, WORK, LWORK, IWORK,
     $                   LIWORK, DOL, DOU, ZOFFSET, INFO )
*
*  -- ScaLAPACK auxiliary routine (version 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver
*     July 4, 2010
*
*     .. Scalar Arguments ..
      CHARACTER          JOBZ, RANGE
      INTEGER            DOL, DOU, IL, INFO, IU, 
     $                   ldz, nzc, liwork, lwork, m, n, zoffset
      DOUBLE PRECISION VL, VU
 
*     ..
*     .. Array Arguments ..
      INTEGER            ISUPPZ( * ), IWORK( * )
      DOUBLE PRECISION   D( * ), E( * ), W( * ), WORK( * )
      DOUBLE PRECISION   Z( LDZ, * )
*     ..
*
*  Purpose
*  =======
*
*  DSTEGR2 computes selected eigenvalues and, optionally, eigenvectors
*  of a real symmetric tridiagonal matrix T. It is invoked in the 
*  ScaLAPACK MRRR driver PDSYEVR and the corresponding Hermitian
*  version either when only eigenvalues are to be computed, or when only
*  a single processor is used (the sequential-like case).
*
*  DSTEGR2 has been adapted from LAPACK's DSTEGR. Please note the
*  following crucial changes.
*
*  1. The calling sequence has two additional INTEGER parameters, 
*     DOL and DOU, that should satisfy M>=DOU>=DOL>=1. 
*     DSTEGR2  ONLY computes the eigenpairs
*     corresponding to eigenvalues DOL through DOU in W. (That is,
*     instead of computing the eigenpairs 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 is 
*     M = DOU - DOL + 1. This concerns the case where only eigenvalues
*     are computed, but on more than one processor. Thus, in this case
*     M refers to the number of eigenvalues computed on this processor.
*  
*  3. 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 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.          
*          Other processors will hold reliable information on other
*          parts of the W array. This information is communicated in
*          the ScaLAPACK driver.
*
*  Z       (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M) )
*          If JOBZ = 'V', and if INFO = 0, then the first M columns of Z
*          contain some of the orthonormal eigenvectors of the matrix T
*          corresponding to the selected eigenvalues, with the i-th
*          column of Z holding the eigenvector associated with W(i).
*          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 can be computed with a workspace
*          query by setting NZC = -1, see below.
*
*  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.
*
*  ISUPPZ  (output) INTEGER ARRAY, dimension ( 2*max(1,M) )
*          The support of the eigenvectors in Z, i.e., the indices
*          indicating the nonzero elements in Z. The i-th computed eigenvector
*          is nonzero only in elements ISUPPZ( 2*i-1 ) through
*          ISUPPZ( 2*i ). This is relevant in the case when the matrix 
*          is split. ISUPPZ is only set if N>2.
*
*  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 the eigenvalues W(1:M), only eigenvectors 
*          Z(:,DOL) to Z(:,DOU) are computed.
*          If DOL > 1, then Z(:,DOL-1-ZOFFSET) is used and overwritten.
*          If DOU < M, then Z(:,DOU+1-ZOFFSET) is used and overwritten.
*
*  ZOFFSET (input) INTEGER
*          Offset for storing the eigenpairs when Z is distributed
*          in 1D-cyclic fashion
*
*  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 DLARRE2,
*                if INFO = 20X, internal error in DLARRV.
*                Here, the digit X = ABS( IINFO ) < 10, where IINFO is 
*                the nonzero error code returned by DLARRE2 or 
*                DLARRV, respectively.
*
*  =====================================================================
*
*     .. 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            I, IIL, IINDBL, IINDW, IINDWK, IINFO, IINSPL,
     $                   iiu, inde2, inderr, indgp, indgrs, indwrk,
     $                   itmp, itmp2, j, jj, liwmin, lwmin, nsplit,
     $                   nzcmin
      DOUBLE PRECISION   BIGNUM, EPS, PIVMIN, RMAX, RMIN, RTOL1, RTOL2,
     $                   SAFMIN, SCALE, SMLNUM, THRESH, TMP, TNRM, WL,
     $                   wu
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      DOUBLE PRECISION   DLAMCH, DLANST
      EXTERNAL           lsame, dlamch, dlanst
*     ..
*     .. External Subroutines ..
      EXTERNAL           dcopy, dlae2, dlaev2, dlarrc, dlarre2,
     $                   dlarrv, dlasrt, dscal, dswap
*     ..
*     .. 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 )
 
*     DSTEGR2 needs WORK of size 6*N, IWORK of size 3*N.
*     In addition, DLARRE2 needs WORK of size 6*N, IWORK of size 5*N.
*     Furthermore, DLARRV needs WORK of size 12*N, IWORK of size 7*N.
      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 DLARRE2.
         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 ) ) )
*
      IF( info.EQ.0 ) THEN
         work( 1 ) = lwmin
         iwork( 1 ) = liwmin
*
         IF( wantz .AND. alleig ) THEN
            nzcmin = n
            iil = 1
            iiu = n
         ELSE IF( wantz .AND. valeig ) THEN
            CALL dlarrc( 'T', n, vl, vu, d, e, safmin, 
     $                            nzcmin, itmp, itmp2, info )
            iil = itmp+1
            iiu = itmp2
         ELSE IF( wantz .AND. indeig ) THEN
            nzcmin = iiu-iil+1
         ELSE 
*           WANTZ .EQ. FALSE.   
            nzcmin = 0
         ENDIF  
         IF( zquery .AND. info.EQ.0 ) THEN
            z( 1,1 ) = nzcmin
         ELSE IF( nzc.LT.nzcmin .AND. .NOT.zquery ) THEN
            info = -14
         END IF
      END IF
 
      IF ( wantz ) THEN
         IF ( dol.LT.1 .OR. dol.GT.nzcmin ) THEN 
            info = -20
         ENDIF
         IF ( dou.LT.1 .OR. dou.GT.nzcmin .OR. dou.LT.dol) THEN 
            info = -21
         ENDIF
      ENDIF
 
      IF( info.NE.0 ) THEN
*
C         Disable sequential error handler
C         for parallel case
C         CALL XERBLA( 'DSTEGR2', -INFO )
*
         RETURN
      ELSE IF( lquery .OR. zquery ) THEN
         RETURN
      END IF
*
*     Quick return if possible
*
      m = 0
      IF( n.EQ.0 )
     $   RETURN
*
      IF( n.EQ.1 ) THEN
         IF( alleig .OR. indeig ) THEN
            m = 1
            w( 1 ) = d( 1 )
         ELSE
            IF( wl.LT.d( 1 ) .AND. wu.GE.d( 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
      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 )
      IF( tnrm.GT.zero .AND. tnrm.LT.rmin ) THEN
         scale = rmin / tnrm
      ELSE IF( tnrm.GT.rmax ) THEN
         scale = rmax / tnrm
      END IF
      IF( scale.NE.one ) 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 DLARRE2      
*     A negative THRESH forces the old splitting criterion based on the
*     size of the off-diagonal. A positive THRESH switches to splitting
*     which preserves relative accuracy. 
*
      iinfo = -1
*     Set the splitting criterion
      IF (iinfo.EQ.0) THEN
         thresh = eps
      ELSE
         thresh = -eps
      ENDIF
*
*     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
      IF( .NOT.wantz ) THEN
*        DLARRE2 computes the eigenvalues to full precision.   
         rtol1 = four * eps
         rtol2 = four * eps
      ELSE   
*        DLARRE2 computes the eigenvalues to less than full precision.
*        DLARRV will refine the eigenvalue approximations, and we can
*        need less accurate initial bisection in DLARRE2.
*        Note: these settings do only affect the subset case and DLARRE2
         rtol1 = sqrt(eps)
         rtol2 = max( sqrt(eps)*5.0d-3, four * eps )
      ENDIF
      CALL dlarre2( range, n, wl, wu, iil, iiu, d, e, 
     $             work(inde2), rtol1, rtol2, thresh, nsplit, 
     $             iwork( iinspl ), m, dol, dou,
     $             w, work( inderr ),
     $             work( indgp ), iwork( iindbl ),
     $             iwork( iindw ), work( indgrs ), pivmin,
     $             work( indwrk ), iwork( iindwk ), iinfo )
      IF( iinfo.NE.0 ) THEN
         info = 100 + abs( iinfo )
         RETURN
      END IF
*     Note that if RANGE .NE. 'V', DLARRE2 computes bounds on the desired
*     part of the spectrum. All desired eigenvalues are contained in
*     (WL,WU]
 
 
      IF( wantz ) THEN
*
*        Compute the desired eigenvectors corresponding to the computed
*        eigenvalues
*
         CALL dlarrv( n, wl, wu, d, e,
     $                pivmin, iwork( iinspl ), m, 
     $                dol, dou, minrgp, rtol1, rtol2, 
     $                w, work( inderr ), work( indgp ), iwork( iindbl ),
     $                iwork( iindw ), work( indgrs ), z, ldz,
     $                isuppz, work( indwrk ), iwork( iindwk ), iinfo )
         IF( iinfo.NE.0 ) THEN
            info = 200 + abs( iinfo )
            RETURN
         END IF
      ELSE
*        DLARRE2 computes eigenvalues of the (shifted) root representation
*        DLARRV returns the eigenvalues of the unshifted matrix.
*        However, if the eigenvectors are not desired by the user, we need
*        to apply the corresponding shifts from DLARRE2 to obtain the 
*        eigenvalues of the original matrix. 
         DO 20 j = 1, m
            itmp = iwork( iindbl+j-1 )
            w( j ) = w( j ) + e( iwork( iinspl+itmp-1 ) )
 20      CONTINUE
      END IF
*
 
*
*     If matrix was scaled, then rescale eigenvalues appropriately.
*
      IF( scale.NE.one ) THEN
         CALL dscal( m, one / scale, w, 1 )
      END IF
*
*     Correct M if needed 
*
      IF ( wantz ) THEN
         IF( dol.NE.1 .OR. dou.NE.m ) THEN
            m = dou - dol +1
         ENDIF
      ENDIF
*
*     If eigenvalues are not in increasing order, then sort them, 
*     possibly along with eigenvectors.
*
      IF( nsplit.GT.1 ) THEN
         IF( .NOT. wantz ) THEN
            CALL dlasrt( 'I', dou - dol +1, w(dol), iinfo )
            IF( iinfo.NE.0 ) THEN
               info = 3
               RETURN
            END IF
         ELSE
            DO 60 j = dol, dou - 1
               i = 0
               tmp = w( j )
               DO 50 jj = j + 1, m
                  IF( w( jj ).LT.tmp ) THEN
                     i = jj
                     tmp = w( jj )
                  END IF
 50            CONTINUE
               IF( i.NE.0 ) THEN
                  w( i ) = w( j )
                  w( j ) = tmp
                  IF( wantz ) THEN
                     CALL dswap( n, z( 1, i-zoffset ), 
     $                                 1, z( 1, j-zoffset ), 1 )
                     itmp = isuppz( 2*i-1 )
                     isuppz( 2*i-1 ) = isuppz( 2*j-1 )
                     isuppz( 2*j-1 ) = itmp
                     itmp = isuppz( 2*i )
                     isuppz( 2*i ) = isuppz( 2*j )
                     isuppz( 2*j ) = itmp
                  END IF
               END IF
 60         CONTINUE
         END IF
      ENDIF
*
      work( 1 ) = lwmin
      iwork( 1 ) = liwmin
      RETURN
*
*     End of DSTEGR2
*
      END