SUBROUTINE PDLAED0( N, D, E, Q, IQ, JQ, DESCQ, WORK, IWORK, INFO ) * * -- ScaLAPACK auxiliary routine (version 1.7) -- * University of Tennessee, Knoxville, Oak Ridge National Laboratory, * and University of California, Berkeley. * December 31, 1998 * * .. Scalar Arguments .. INTEGER INFO, IQ, JQ, N * .. * .. Array Arguments .. INTEGER DESCQ( * ), IWORK( * ) DOUBLE PRECISION D( * ), E( * ), Q( * ), WORK( * ) * .. * * Purpose * ======= * * PDLAED0 computes all eigenvalues and corresponding eigenvectors of a * symmetric tridiagonal matrix using the divide and conquer method. * * * Arguments * ========= * * 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, * global dimension (N, N), * 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 ) DOUBLE PRECISION array, dimension (LWORK) * LWORK = 6*N + 2*NP*NQ, with * NP = NUMROC( N, MB_Q, MYROW, IQROW, NPROW ) * NQ = NUMROC( N, NB_Q, MYCOL, IQCOL, NPCOL ) * IQROW = INDXG2P( IQ, NB_Q, MYROW, RSRC_Q, NPROW ) * IQCOL = INDXG2P( JQ, MB_Q, MYCOL, CSRC_Q, NPCOL ) * * IWORK (local workspace/output) INTEGER array, dimension (LIWORK) * 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). * * ===================================================================== * * .. 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 ) * .. * .. Local Scalars .. INTEGER I, ID, IDCOL, IDROW, IID, IINFO, IIQ, IM1, IM2, $ IPQ, IQCOL, IQROW, J, JJD, JJQ, LDQ, MATSIZ, $ MYCOL, MYROW, N1, NB, NBL, NBL1, NPCOL, NPROW, $ SUBPBS, TSUBPBS * .. * .. External Subroutines .. EXTERNAL BLACS_GRIDINFO, DGEBR2D, DGEBS2D, DGERV2D, $ DGESD2D, DSTEQR, INFOG2L, PDLAED1, PXERBLA * .. * .. External Functions .. * .. * .. Intrinsic Functions .. INTRINSIC ABS, MIN * .. * .. 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 ) INFO = 0 IF( DESCQ( NB_ ).GT.N .OR. N.LT.2 ) $ INFO = -1 IF( INFO.NE.0 ) THEN CALL PXERBLA( DESCQ( CTXT_ ), 'PDLAED0', -INFO ) RETURN END IF * NB = DESCQ( NB_ ) LDQ = DESCQ( LLD_ ) CALL INFOG2L( IQ, JQ, DESCQ, NPROW, NPCOL, MYROW, MYCOL, IIQ, JJQ, $ IQROW, IQCOL ) * * Determine the size and placement of the submatrices, and save in * the leading elements of IWORK. * TSUBPBS = ( N-1 ) / NB + 1 IWORK( 1 ) = TSUBPBS SUBPBS = 1 10 CONTINUE IF( IWORK( SUBPBS ).GT.1 ) THEN DO 20 J = SUBPBS, 1, -1 IWORK( 2*J ) = ( IWORK( J )+1 ) / 2 IWORK( 2*J-1 ) = IWORK( J ) / 2 20 CONTINUE SUBPBS = 2*SUBPBS GO TO 10 END IF DO 30 J = 2, SUBPBS IWORK( J ) = IWORK( J ) + IWORK( J-1 ) 30 CONTINUE * * Divide the matrix into TSUBPBS submatrices of size at most NB * using rank-1 modifications (cuts). * DO 40 I = NB + 1, N, NB IM1 = I - 1 D( IM1 ) = D( IM1 ) - ABS( E( IM1 ) ) D( I ) = D( I ) - ABS( E( IM1 ) ) 40 CONTINUE * * Solve each submatrix eigenproblem at the bottom of the divide and * conquer tree. D is the same on each process. * DO 50 ID = 1, N, NB CALL INFOG2L( IQ-1+ID, JQ-1+ID, DESCQ, NPROW, NPCOL, MYROW, $ MYCOL, IID, JJD, IDROW, IDCOL ) MATSIZ = MIN( NB, N-ID+1 ) IF( MYROW.EQ.IDROW .AND. MYCOL.EQ.IDCOL ) THEN IPQ = IID + ( JJD-1 )*LDQ CALL DSTEQR( 'I', MATSIZ, D( ID ), E( ID ), Q( IPQ ), LDQ, $ WORK, INFO ) IF( INFO.NE.0 ) THEN CALL PXERBLA( DESCQ( CTXT_ ), 'DSTEQR', -INFO ) RETURN END IF IF( MYROW.NE.IQROW .OR. MYCOL.NE.IQCOL ) THEN CALL DGESD2D( DESCQ( CTXT_ ), MATSIZ, 1, D( ID ), MATSIZ, $ IQROW, IQCOL ) END IF ELSE IF( MYROW.EQ.IQROW .AND. MYCOL.EQ.IQCOL ) THEN CALL DGERV2D( DESCQ( CTXT_ ), MATSIZ, 1, D( ID ), MATSIZ, $ IDROW, IDCOL ) END IF 50 CONTINUE * IF( MYROW.EQ.IQROW .AND. MYCOL.EQ.IQCOL ) THEN CALL DGEBS2D( DESCQ( CTXT_ ), 'A', ' ', N, 1, D, N ) ELSE CALL DGEBR2D( DESCQ( CTXT_ ), 'A', ' ', N, 1, D, N, IQROW, $ IQCOL ) END IF * * Successively merge eigensystems of adjacent submatrices * into eigensystem for the corresponding larger matrix. * * while ( SUBPBS > 1 ) * 60 CONTINUE IF( SUBPBS.GT.1 ) THEN IM2 = SUBPBS - 2 DO 80 I = 0, IM2, 2 IF( I.EQ.0 ) THEN NBL = IWORK( 2 ) NBL1 = IWORK( 1 ) IF( NBL1.EQ.0 ) $ GO TO 70 ID = 1 MATSIZ = MIN( N, NBL*NB ) N1 = NBL1*NB ELSE NBL = IWORK( I+2 ) - IWORK( I ) NBL1 = NBL / 2 IF( NBL1.EQ.0 ) $ GO TO 70 ID = IWORK( I )*NB + 1 MATSIZ = MIN( NB*NBL, N-ID+1 ) N1 = NBL1*NB END IF * * Merge lower order eigensystems (of size N1 and MATSIZ - N1) * into an eigensystem of size MATSIZ. * CALL PDLAED1( MATSIZ, N1, D( ID ), ID, Q, IQ, JQ, DESCQ, $ E( ID+N1-1 ), WORK, IWORK( SUBPBS+1 ), IINFO ) IF( IINFO.NE.0 ) THEN INFO = IINFO*( N+1 ) + ID END IF * 70 CONTINUE IWORK( I / 2+1 ) = IWORK( I+2 ) 80 CONTINUE SUBPBS = SUBPBS / 2 * GO TO 60 END IF * * end while * 90 CONTINUE RETURN * * End of PDLAED0 * END