sstegr2b.f

Go to the documentation of this file.

      SUBROUTINE SSTEGR2B( JOBZ, N, D, E, 
     $                   M, W, Z, LDZ, NZC, ISUPPZ, WORK, LWORK, IWORK,
     $                   LIWORK, DOL, DOU, NEEDIL, NEEDIU,
     $                   INDWLC, PIVMIN, SCALE, WL, WU,
     $                   VSTART, FINISH, MAXCLS,
     $                   NDEPTH, PARITY, ZOFFSET, 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
      INTEGER            DOL, DOU, INDWLC, INFO, LDZ, LIWORK, LWORK, M,
     $                   MAXCLS, N, NDEPTH, NEEDIL, NEEDIU, NZC, PARITY,
     $                   ZOFFSET

      REAL             PIVMIN, SCALE, WL, WU
      LOGICAL VSTART, FINISH

*     ..
*     .. Array Arguments ..
      INTEGER            ISUPPZ( * ), IWORK( * )
      REAL               D( * ), E( * ), W( * ), WORK( * )
      REAL               Z( LDZ, * )
*     ..
*
*  Purpose
*  =======
*
*  SSTEGR2B should only be called after a call to SSTEGR2A.
*  From eigenvalues and initial representations computed by SSTEGR2A,
*  SSTEGR2B computes the selected eigenvalues and eigenvectors
*  of the real symmetric tridiagonal matrix in parallel 
*  on multiple processors. It is potentially invoked multiple times
*  on a given processor because the locally relevant representation tree 
*  might depend on spectral information that is "owned" by other processors
*  and might need to be communicated. 
* 
*  Please note:
*  1. The calling sequence has two additional INTEGER parameters, 
*     DOL and DOU, that should satisfy M>=DOU>=DOL>=1. 
*     These parameters are only relevant for the case JOBZ = 'V'.
*     SSTEGR2B  ONLY computes the eigenVECTORS 
*     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 accurately refined
*     to all figures by Rayleigh-Quotient iteration.
*
*  2. The additional arguments VSTART, FINISH, NDEPTH, PARITY, ZOFFSET 
*     are included as a thread-safe implementation equivalent to SAVE variables.
*     These variables store details about the local representation tree which is
*     computed layerwise. For scalability reasons, eigenvalues belonging to the 
*     locally relevant representation tree might be computed on other processors.
*     These need to be communicated before the inspection of the RRRs can proceed
*     on any given layer.           
*     Note that only when the variable FINISH is true, the computation has ended
*     All eigenpairs between DOL and DOU have been computed. M is set = DOU - DOL + 1.
*
*  3. SSTEGR2B needs more workspace in Z than the sequential SSTEGR. 
*     It is used to store the conformal embedding of the local representation tree.  
*  
*  Arguments
*  =========
*
*  JOBZ    (input) CHARACTER*1
*          = 'N':  Compute eigenvalues only;
*          = 'V':  Compute eigenvalues and eigenvectors.
*
*  N       (input) INTEGER
*          The order of the matrix.  N >= 0.
*
*  D       (input/output) REAL array, dimension (N)
*          On entry, the N diagonal elements of the tridiagonal matrix
*          T. On exit, D is overwritten.
*
*  E       (input/output) REAL 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.
*
*  M       (input) INTEGER
*          The total number of eigenvalues found
*          in SSTEGR2A.  0 <= M <= N.
*
*  W       (input) REAL array, dimension (N)
*          The first M elements contain approximations to the selected 
*          eigenvalues in ascending order. Note that only the eigenvalues from
*          the locally relevant part of the representation tree, that is
*          all the clusters that include eigenvalues from DOL:DOU, are reliable 
*          on this processor. (It does not need to know about any others anyway.)
*
*  Z       (output) REAL array, dimension (LDZ, max(1,M) )
*          If JOBZ = 'V', and if INFO = 0, then 
*          a subset of the first M columns of Z
*          contain 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).
*          See DOL, DOU for more information.
*
*  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.  
*
*  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) REAL 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.
*
*  NEEDIL  (input/output) INTEGER 
*  NEEDIU  (input/output) INTEGER
*          Describes which are the left and right outermost eigenvalues 
*          still to be computed. Initially computed by SLARRE2A,
*          modified in the course of the algorithm.
*
*  INDWLC  (output) REAL
*          Pointer into the workspace, location where the local
*          eigenvalue representations are stored. ("Local eigenvalues"
*          are those relative to the individual shifts of the RRRs.)
*
*  PIVMIN  (input) REAL
*          The minimum pivot in the sturm sequence for T.
*
*  SCALE   (input) REAL 
*          The scaling factor for T. Used for unscaling the eigenvalues
*          at the very end of the algorithm.
*
*  WL      (input) REAL
*  WU      (input) REAL
*          The interval (WL, WU] contains all the wanted eigenvalues.         
*
*  VSTART  (input/output) LOGICAL 
*          .TRUE. on initialization, set to .FALSE. afterwards.
*
*  FINISH  (input/output) LOGICAL
*          indicates whether all eigenpairs have been computed
*
*  MAXCLS  (input/output) INTEGER
*          The largest cluster worked on by this processor in the
*          representation tree.
*
*  NDEPTH  (input/output) INTEGER
*          The current depth of the representation tree. Set to
*          zero on initial pass, changed when the deeper levels of
*          the representation tree are generated. 
*
*  PARITY  (input/output) INTEGER
*          An internal parameter needed for the storage of the
*          clusters on the current level of the representation tree.
*
*  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 = 20X, internal error in SLARRV2.
*                Here, the digit X = ABS( IINFO ) < 10, where IINFO is 
*                the nonzero error code returned by SLARRV2.
*
*     .. Parameters ..
      REAL               ONE, FOUR, MINRGP
      PARAMETER          ( ONE = 1.0E0,
     $                     FOUR = 4.0E0,
     $                     MINRGP = 3.0E-3 )
*     ..
*     .. Local Scalars ..
      LOGICAL            LQUERY, WANTZ, ZQUERY
      INTEGER            IINDBL, IINDW, IINDWK, IINFO, IINSPL, INDERR,
     $                   INDGP, INDGRS, INDSDM, INDWRK, ITMP, J, LIWMIN,
     $                   LWMIN
      REAL               EPS, RTOL1, RTOL2
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      REAL               SLAMCH, SLANST
      EXTERNAL           LSAME, SLAMCH, SLANST
*     ..
*     .. External Subroutines ..
      EXTERNAL           SLARRV2, SSCAL
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MAX, MIN, REAL, SQRT
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      WANTZ = LSAME( JOBZ, 'V' )
*
      LQUERY = ( ( LWORK.EQ.-1 ).OR.( LIWORK.EQ.-1 ) )
      ZQUERY = ( NZC.EQ.-1 )

*     SSTEGR2B needs WORK of size 6*N, IWORK of size 3*N.
*     In addition, SLARRE2A needed WORK of size 6*N, IWORK of size 5*N.
*     Workspace is kept consistent even though SLARRE2A is not called here.
*     Furthermore, SLARRV2 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
*
      INFO = 0
*
*     Get machine constants.
*
      EPS = SLAMCH( 'Precision' )
*
      IF( (N.EQ.0).OR.(N.EQ.1) ) THEN 
         FINISH = .TRUE.       
         RETURN
      ENDIF

      IF(ZQUERY.OR.LQUERY)
     $   RETURN
*
      INDGRS = 1
      INDERR = 2*N + 1
      INDGP = 3*N + 1
      INDSDM = 4*N + 1
      INDWRK = 6*N + 1
      INDWLC = INDWRK
*
      IINSPL = 1
      IINDBL = N + 1
      IINDW = 2*N + 1
      IINDWK = 3*N + 1

*     Set the tolerance parameters for bisection
      RTOL1 = FOUR*SQRT(EPS)
      RTOL2 = MAX( SQRT(EPS)*5.0E-3, FOUR * EPS )


      IF( WANTZ ) THEN
*
*        Compute the desired eigenvectors corresponding to the computed
*        eigenvalues
*
         CALL SLARRV2( N, WL, WU, D, E,
     $                PIVMIN, IWORK( IINSPL ), M, 
     $                DOL, DOU, NEEDIL, NEEDIU, MINRGP, RTOL1, RTOL2, 
     $                W, WORK( INDERR ), WORK( INDGP ), IWORK( IINDBL ),
     $                IWORK( IINDW ), WORK( INDGRS ), 
     $                WORK( INDSDM ), Z, LDZ,
     $                ISUPPZ, WORK( INDWRK ), IWORK( IINDWK ), 
     $                VSTART, FINISH, 
     $                MAXCLS, NDEPTH, PARITY, ZOFFSET, IINFO )
         IF( IINFO.NE.0 ) THEN
            INFO = 200 + ABS( IINFO )
            RETURN
         END IF
*
      ELSE
*        SLARRE2A computed eigenvalues of the (shifted) root representation
*        SLARRV2 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 SLARRE2A to obtain the 
*        eigenvalues of the original matrix. 
         DO 30 J = 1, M
            ITMP = IWORK( IINDBL+J-1 )
            W( J ) = W( J ) + E( IWORK( IINSPL+ITMP-1 ) )
 30      CONTINUE
*
         FINISH = .TRUE.
*
      END IF
*

      IF(FINISH) THEN        
*        All eigenpairs have been computed       

*
*        If matrix was scaled, then rescale eigenvalues appropriately.
*
         IF( SCALE.NE.ONE ) THEN
            CALL SSCAL( 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
*
*        No sorting of eigenpairs is done here, done later in the
*        calling subroutine
*
         WORK( 1 ) = LWMIN
         IWORK( 1 ) = LIWMIN
      ENDIF

      RETURN
*
*     End of SSTEGR2B
*
      END