SUBROUTINE DSTEGR2B( 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 DOUBLE PRECISION PIVMIN, SCALE, WL, WU LOGICAL VSTART, FINISH * .. * .. Array Arguments .. INTEGER ISUPPZ( * ), IWORK( * ) DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ) DOUBLE PRECISION Z( LDZ, * ) * .. * * Purpose * ======= * * DSTEGR2B should only be called after a call to DSTEGR2A. * From eigenvalues and initial representations computed by DSTEGR2A, * DSTEGR2B 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'. * DSTEGR2B 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. DSTEGR2B needs more workspace in Z than the sequential DSTEGR. * 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) 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. * * M (input) INTEGER * The total number of eigenvalues found * in DSTEGR2A. 0 <= M <= N. * * W (input) DOUBLE PRECISION 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) DOUBLE PRECISION 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) 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. * * NEEDIL (input/output) INTEGER * NEEDIU (input/output) INTEGER * Describes which are the left and right outermost eigenvalues * still to be computed. Initially computed by DLARRE2A, * modified in the course of the algorithm. * * INDWLC (output) DOUBLE PRECISION * 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) DOUBLE PRECISION * The minimum pivot in the sturm sequence for T. * * SCALE (input) DOUBLE PRECISION * The scaling factor for T. Used for unscaling the eigenvalues * at the very end of the algorithm. * * WL (input) DOUBLE PRECISION * WU (input) DOUBLE PRECISION * 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 DLARRV2. * Here, the digit X = ABS( IINFO ) < 10, where IINFO is * the nonzero error code returned by DLARRV2. * * .. Parameters .. DOUBLE PRECISION ONE, FOUR, MINRGP PARAMETER ( ONE = 1.0D0, $ FOUR = 4.0D0, $ MINRGP = 1.0D-3 ) * .. * .. Local Scalars .. LOGICAL LQUERY, WANTZ, ZQUERY INTEGER IINDBL, IINDW, IINDWK, IINFO, IINSPL, INDERR, $ INDGP, INDGRS, INDSDM, INDWRK, ITMP, J, LIWMIN, $ LWMIN DOUBLE PRECISION EPS, RTOL1, RTOL2 * .. * .. External Functions .. LOGICAL LSAME DOUBLE PRECISION DLAMCH, DLANST EXTERNAL LSAME, DLAMCH, DLANST * .. * .. External Subroutines .. EXTERNAL DLARRV2, DSCAL * .. * .. Intrinsic Functions .. INTRINSIC DBLE, MAX, MIN, SQRT * .. * .. Executable Statements .. * * Test the input parameters. * WANTZ = LSAME( JOBZ, 'V' ) * LQUERY = ( ( LWORK.EQ.-1 ).OR.( LIWORK.EQ.-1 ) ) ZQUERY = ( NZC.EQ.-1 ) * DSTEGR2B needs WORK of size 6*N, IWORK of size 3*N. * In addition, DLARRE2A needed WORK of size 6*N, IWORK of size 5*N. * Workspace is kept consistent even though DLARRE2A is not called here. * Furthermore, DLARRV2 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 = DLAMCH( '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.0D-3, FOUR * EPS ) IF( WANTZ ) THEN * * Compute the desired eigenvectors corresponding to the computed * eigenvalues * CALL DLARRV2( 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 * DLARRE2A computed eigenvalues of the (shifted) root representation * DLARRV2 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 DLARRE2A 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 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 * * No sorting of eigenpairs is done here, done later in the * calling subroutine * WORK( 1 ) = LWMIN IWORK( 1 ) = LIWMIN ENDIF RETURN * * End of DSTEGR2B * END