SUBROUTINE PDSTEDC( COMPZ, N, D, E, Q, IQ, JQ, DESCQ, WORK, LWORK, $ IWORK, LIWORK, INFO ) * * -- ScaLAPACK routine (version 1.7) -- * University of Tennessee, Knoxville, Oak Ridge National Laboratory, * and University of California, Berkeley. * March 13, 2000 * * .. Scalar Arguments .. CHARACTER COMPZ INTEGER INFO, IQ, JQ, LIWORK, LWORK, N * .. * .. Array Arguments .. INTEGER DESCQ( * ), IWORK( * ) DOUBLE PRECISION D( * ), E( * ), Q( * ), WORK( * ) * .. * * Purpose * ======= * PDSTEDC computes all eigenvalues and eigenvectors of a * symmetric tridiagonal matrix in parallel, using the divide and * conquer algorithm. * * This code makes very mild assumptions about floating point * arithmetic. It will work on machines with a guard digit in * add/subtract, or on those binary machines without guard digits * which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. * It could conceivably fail on hexadecimal or decimal machines * without guard digits, but we know of none. See DLAED3 for details. * * Arguments * ========= * * COMPZ (input) CHARACTER*1 * = 'N': Compute eigenvalues only. (NOT IMPLEMENTED YET) * = 'I': Compute eigenvectors of tridiagonal matrix also. * = 'V': Compute eigenvectors of original dense symmetric * matrix also. On entry, Z contains the orthogonal * matrix used to reduce the original matrix to * tridiagonal form. (NOT IMPLEMENTED YET) * * N (global input) INTEGER * The order of the tridiagonal matrix T. N >= 0. * * D (global input/output) DOUBLE PRECISION array, dimension (N) * On entry, the diagonal elements of the tridiagonal matrix. * On exit, if INFO = 0, the eigenvalues in descending order. * * E (global input/output) DOUBLE PRECISION array, dimension (N-1) * On entry, the subdiagonal elements of the tridiagonal matrix. * On exit, E has been destroyed. * * Q (local output) DOUBLE PRECISION array, * local dimension ( LLD_Q, LOCc(JQ+N-1)) * Q contains the orthonormal eigenvectors of the symmetric * tridiagonal matrix. * On output, Q is distributed across the P processes in block * cyclic format. * * IQ (global input) INTEGER * Q's global row index, which points to the beginning of the * submatrix which is to be operated on. * * JQ (global input) INTEGER * Q's global column index, which points to the beginning of * the submatrix which is to be operated on. * * DESCQ (global and local input) INTEGER array of dimension DLEN_. * The array descriptor for the distributed matrix Z. * * * WORK (local workspace/output) DOUBLE PRECISION array, * dimension (LWORK) * On output, WORK(1) returns the workspace needed. * * LWORK (local input/output) INTEGER, * the dimension of the array WORK. * LWORK = 6*N + 2*NP*NQ * NP = NUMROC( N, NB, MYROW, DESCQ( RSRC_ ), NPROW ) * NQ = NUMROC( N, NB, MYCOL, DESCQ( CSRC_ ), NPCOL ) * * If LWORK = -1, the LWORK is global input and a workspace * query is assumed; the routine only calculates the minimum * size for the WORK array. The required workspace is returned * as the first element of WORK and no error message is issued * by PXERBLA. * * IWORK (local workspace/output) INTEGER array, dimension (LIWORK) * On exit, if LIWORK > 0, IWORK(1) returns the optimal LIWORK. * * LIWORK (input) INTEGER * The dimension of the array IWORK. * LIWORK = 2 + 7*N + 8*NPCOL * * INFO (global output) INTEGER * = 0: successful exit * < 0: If the i-th argument is an array and the j-entry had * an illegal value, then INFO = -(i*100+j), if the i-th * argument is a scalar and had an illegal value, then * INFO = -i. * > 0: The algorithm failed to compute the INFO/(N+1) th * eigenvalue while working on the submatrix lying in * global rows and columns mod(INFO,N+1). * * Further Details * ======= ======= * * Contributed by Francoise Tisseur, University of Manchester. * * Reference: F. Tisseur and J. Dongarra, "A Parallel Divide and * Conquer Algorithm for the Symmetric Eigenvalue Problem * on Distributed Memory Architectures", * SIAM J. Sci. Comput., 6:20 (1999), pp. 2223--2236. * (see also LAPACK Working Note 132) * http://www.netlib.org/lapack/lawns/lawn132.ps * * ===================================================================== * * .. Parameters .. INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_, $ MB_, NB_, RSRC_, CSRC_, LLD_ PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1, $ CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, $ RSRC_ = 7, CSRC_ = 8, LLD_ = 9 ) DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) * .. * .. Local Scalars .. LOGICAL LQUERY INTEGER ICOFFQ, IIQ, IPQ, IQCOL, IQROW, IROFFQ, JJQ, $ LDQ, LIWMIN, LWMIN, MYCOL, MYROW, NB, NP, $ NPCOL, NPROW, NQ DOUBLE PRECISION ORGNRM * .. * .. External Functions .. LOGICAL LSAME INTEGER INDXG2P, NUMROC DOUBLE PRECISION DLANST EXTERNAL INDXG2P, LSAME, NUMROC, DLANST * .. * .. External Subroutines .. EXTERNAL BLACS_GRIDINFO, CHK1MAT, DLASCL, DSTEDC, $ INFOG2L, PDLAED0, PDLASRT, PXERBLA * .. * .. Intrinsic Functions .. INTRINSIC DBLE, MOD * .. * .. Executable Statements .. * * This is just to keep ftnchek and toolpack/1 happy IF( BLOCK_CYCLIC_2D*CSRC_*CTXT_*DLEN_*DTYPE_*LLD_*MB_*M_*NB_*N_* $ RSRC_.LT.0 )RETURN * * Test the input parameters. * CALL BLACS_GRIDINFO( DESCQ( CTXT_ ), NPROW, NPCOL, MYROW, MYCOL ) LDQ = DESCQ( LLD_ ) NB = DESCQ( NB_ ) NP = NUMROC( N, NB, MYROW, DESCQ( RSRC_ ), NPROW ) NQ = NUMROC( N, NB, MYCOL, DESCQ( CSRC_ ), NPCOL ) INFO = 0 IF( NPROW.EQ.-1 ) THEN INFO = -( 600+CTXT_ ) ELSE CALL CHK1MAT( N, 2, N, 2, IQ, JQ, DESCQ, 8, INFO ) IF( INFO.EQ.0 ) THEN NB = DESCQ( NB_ ) IROFFQ = MOD( IQ-1, DESCQ( MB_ ) ) ICOFFQ = MOD( JQ-1, DESCQ( NB_ ) ) IQROW = INDXG2P( IQ, NB, MYROW, DESCQ( RSRC_ ), NPROW ) IQCOL = INDXG2P( JQ, NB, MYCOL, DESCQ( CSRC_ ), NPCOL ) LWMIN = 6*N + 2*NP*NQ LIWMIN = 2 + 7*N + 8*NPCOL * WORK( 1 ) = DBLE( LWMIN ) IWORK( 1 ) = LIWMIN LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) IF( .NOT.LSAME( COMPZ, 'I' ) ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( IROFFQ.NE.ICOFFQ .OR. ICOFFQ.NE.0 ) THEN INFO = -5 ELSE IF( DESCQ( MB_ ).NE.DESCQ( NB_ ) ) THEN INFO = -( 700+NB_ ) ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN INFO = -10 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN INFO = -12 END IF END IF END IF IF( INFO.NE.0 ) THEN CALL PXERBLA( DESCQ( CTXT_ ), 'PDSTEDC', -INFO ) RETURN ELSE IF( LQUERY ) THEN RETURN END IF * * Quick return * IF( N.EQ.0 ) $ GO TO 10 CALL INFOG2L( IQ, JQ, DESCQ, NPROW, NPCOL, MYROW, MYCOL, IIQ, JJQ, $ IQROW, IQCOL ) IF( N.EQ.1 ) THEN IF( MYROW.EQ.IQROW .AND. MYCOL.EQ.IQCOL ) $ Q( 1 ) = ONE GO TO 10 END IF * * If N is smaller than the minimum divide size NB, then * solve the problem with the serial divide and conquer * code locally. * IF( N.LE.NB ) THEN IF( ( MYROW.EQ.IQROW ) .AND. ( MYCOL.EQ.IQCOL ) ) THEN IPQ = IIQ + ( JJQ-1 )*LDQ CALL DSTEDC( 'I', N, D, E, Q( IPQ ), LDQ, WORK, LWORK, $ IWORK, LIWORK, INFO ) IF( INFO.NE.0 ) THEN INFO = ( N+1 ) + N GO TO 10 END IF END IF GO TO 10 END IF * * If P=NPROW*NPCOL=1, solve the problem with DSTEDC. * IF( NPCOL*NPROW.EQ.1 ) THEN IPQ = IIQ + ( JJQ-1 )*LDQ CALL DSTEDC( 'I', N, D, E, Q( IPQ ), LDQ, WORK, LWORK, IWORK, $ LIWORK, INFO ) GO TO 10 END IF * * Scale matrix to allowable range, if necessary. * ORGNRM = DLANST( 'M', N, D, E ) IF( ORGNRM.NE.ZERO ) THEN CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO ) CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N-1, 1, E, N-1, INFO ) END IF * CALL PDLAED0( N, D, E, Q, IQ, JQ, DESCQ, WORK, IWORK, INFO ) * * Sort eigenvalues and corresponding eigenvectors * CALL PDLASRT( 'I', N, D, Q, IQ, JQ, DESCQ, WORK, LWORK, IWORK, $ LIWORK, INFO ) * * Scale back. * IF( ORGNRM.NE.ZERO ) $ CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO ) * 10 CONTINUE * IF( LWORK.GT.0 ) $ WORK( 1 ) = DBLE( LWMIN ) IF( LIWORK.GT.0 ) $ IWORK( 1 ) = LIWMIN RETURN * * End of PDSTEDC * END