1 SUBROUTINE PDLAED0( N, D, E, Q, IQ, JQ, DESCQ, WORK, IWORK, INFO ) 2* 3* -- ScaLAPACK auxiliary routine (version 1.7) -- 4* University of Tennessee, Knoxville, Oak Ridge National Laboratory, 5* and University of California, Berkeley. 6* December 31, 1998 7* 8* .. Scalar Arguments .. 9 INTEGER INFO, IQ, JQ, N 10* .. 11* .. Array Arguments .. 12 INTEGER DESCQ( * ), IWORK( * ) 13 DOUBLE PRECISION D( * ), E( * ), Q( * ), WORK( * ) 14* .. 15* 16* Purpose 17* ======= 18* 19* PDLAED0 computes all eigenvalues and corresponding eigenvectors of a 20* symmetric tridiagonal matrix using the divide and conquer method. 21* 22* 23* Arguments 24* ========= 25* 26* N (global input) INTEGER 27* The order of the tridiagonal matrix T. N >= 0. 28* 29* D (global input/output) DOUBLE PRECISION array, dimension (N) 30* On entry, the diagonal elements of the tridiagonal matrix. 31* On exit, if INFO = 0, the eigenvalues in descending order. 32* 33* E (global input/output) DOUBLE PRECISION array, dimension (N-1) 34* On entry, the subdiagonal elements of the tridiagonal matrix. 35* On exit, E has been destroyed. 36* 37* Q (local output) DOUBLE PRECISION array, 38* global dimension (N, N), 39* local dimension ( LLD_Q, LOCc(JQ+N-1)) 40* Q contains the orthonormal eigenvectors of the symmetric 41* tridiagonal matrix. 42* On output, Q is distributed across the P processes in block 43* cyclic format. 44* 45* IQ (global input) INTEGER 46* Q's global row index, which points to the beginning of the 47* submatrix which is to be operated on. 48* 49* JQ (global input) INTEGER 50* Q's global column index, which points to the beginning of 51* the submatrix which is to be operated on. 52* 53* DESCQ (global and local input) INTEGER array of dimension DLEN_. 54* The array descriptor for the distributed matrix Z. 55* 56* 57* WORK (local workspace ) DOUBLE PRECISION array, dimension (LWORK) 58* LWORK = 6*N + 2*NP*NQ, with 59* NP = NUMROC( N, MB_Q, MYROW, IQROW, NPROW ) 60* NQ = NUMROC( N, NB_Q, MYCOL, IQCOL, NPCOL ) 61* IQROW = INDXG2P( IQ, NB_Q, MYROW, RSRC_Q, NPROW ) 62* IQCOL = INDXG2P( JQ, MB_Q, MYCOL, CSRC_Q, NPCOL ) 63* 64* IWORK (local workspace/output) INTEGER array, dimension (LIWORK) 65* LIWORK = 2 + 7*N + 8*NPCOL 66* 67* INFO (global output) INTEGER 68* = 0: successful exit 69* < 0: If the i-th argument is an array and the j-entry had 70* an illegal value, then INFO = -(i*100+j), if the i-th 71* argument is a scalar and had an illegal value, then 72* INFO = -i. 73* > 0: The algorithm failed to compute the INFO/(N+1) th 74* eigenvalue while working on the submatrix lying in 75* global rows and columns mod(INFO,N+1). 76* 77* ===================================================================== 78* 79* .. Parameters .. 80* 81 INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_, 82 $ MB_, NB_, RSRC_, CSRC_, LLD_ 83 PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1, 84 $ CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, 85 $ RSRC_ = 7, CSRC_ = 8, LLD_ = 9 ) 86* .. 87* .. Local Scalars .. 88 INTEGER I, ID, IDCOL, IDROW, IID, IINFO, IIQ, IM1, IM2, 89 $ IPQ, IQCOL, IQROW, J, JJD, JJQ, LDQ, MATSIZ, 90 $ MYCOL, MYROW, N1, NB, NBL, NBL1, NPCOL, NPROW, 91 $ SUBPBS, TSUBPBS 92* .. 93* .. External Subroutines .. 94 EXTERNAL BLACS_GRIDINFO, DGEBR2D, DGEBS2D, DGERV2D, 95 $ DGESD2D, DSTEQR, INFOG2L, PDLAED1, PXERBLA 96* .. 97* .. External Functions .. 98* .. 99* .. Intrinsic Functions .. 100 INTRINSIC ABS, MIN 101* .. 102* .. Executable Statements .. 103* 104* This is just to keep ftnchek and toolpack/1 happy 105 IF( BLOCK_CYCLIC_2D*CSRC_*CTXT_*DLEN_*DTYPE_*LLD_*MB_*M_*NB_*N_* 106 $ RSRC_.LT.0 )RETURN 107* 108* Test the input parameters. 109* 110 CALL BLACS_GRIDINFO( DESCQ( CTXT_ ), NPROW, NPCOL, MYROW, MYCOL ) 111 INFO = 0 112 IF( DESCQ( NB_ ).GT.N .OR. N.LT.2 ) 113 $ INFO = -1 114 IF( INFO.NE.0 ) THEN 115 CALL PXERBLA( DESCQ( CTXT_ ), 'PDLAED0', -INFO ) 116 RETURN 117 END IF 118* 119 NB = DESCQ( NB_ ) 120 LDQ = DESCQ( LLD_ ) 121 CALL INFOG2L( IQ, JQ, DESCQ, NPROW, NPCOL, MYROW, MYCOL, IIQ, JJQ, 122 $ IQROW, IQCOL ) 123* 124* Determine the size and placement of the submatrices, and save in 125* the leading elements of IWORK. 126* 127 TSUBPBS = ( N-1 ) / NB + 1 128 IWORK( 1 ) = TSUBPBS 129 SUBPBS = 1 130 10 CONTINUE 131 IF( IWORK( SUBPBS ).GT.1 ) THEN 132 DO 20 J = SUBPBS, 1, -1 133 IWORK( 2*J ) = ( IWORK( J )+1 ) / 2 134 IWORK( 2*J-1 ) = IWORK( J ) / 2 135 20 CONTINUE 136 SUBPBS = 2*SUBPBS 137 GO TO 10 138 END IF 139 DO 30 J = 2, SUBPBS 140 IWORK( J ) = IWORK( J ) + IWORK( J-1 ) 141 30 CONTINUE 142* 143* Divide the matrix into TSUBPBS submatrices of size at most NB 144* using rank-1 modifications (cuts). 145* 146 DO 40 I = NB + 1, N, NB 147 IM1 = I - 1 148 D( IM1 ) = D( IM1 ) - ABS( E( IM1 ) ) 149 D( I ) = D( I ) - ABS( E( IM1 ) ) 150 40 CONTINUE 151* 152* Solve each submatrix eigenproblem at the bottom of the divide and 153* conquer tree. D is the same on each process. 154* 155 DO 50 ID = 1, N, NB 156 CALL INFOG2L( IQ-1+ID, JQ-1+ID, DESCQ, NPROW, NPCOL, MYROW, 157 $ MYCOL, IID, JJD, IDROW, IDCOL ) 158 MATSIZ = MIN( NB, N-ID+1 ) 159 IF( MYROW.EQ.IDROW .AND. MYCOL.EQ.IDCOL ) THEN 160 IPQ = IID + ( JJD-1 )*LDQ 161 CALL DSTEQR( 'I', MATSIZ, D( ID ), E( ID ), Q( IPQ ), LDQ, 162 $ WORK, INFO ) 163 IF( INFO.NE.0 ) THEN 164 CALL PXERBLA( DESCQ( CTXT_ ), 'DSTEQR', -INFO ) 165 RETURN 166 END IF 167 IF( MYROW.NE.IQROW .OR. MYCOL.NE.IQCOL ) THEN 168 CALL DGESD2D( DESCQ( CTXT_ ), MATSIZ, 1, D( ID ), MATSIZ, 169 $ IQROW, IQCOL ) 170 END IF 171 ELSE IF( MYROW.EQ.IQROW .AND. MYCOL.EQ.IQCOL ) THEN 172 CALL DGERV2D( DESCQ( CTXT_ ), MATSIZ, 1, D( ID ), MATSIZ, 173 $ IDROW, IDCOL ) 174 END IF 175 50 CONTINUE 176* 177 IF( MYROW.EQ.IQROW .AND. MYCOL.EQ.IQCOL ) THEN 178 CALL DGEBS2D( DESCQ( CTXT_ ), 'A', ' ', N, 1, D, N ) 179 ELSE 180 CALL DGEBR2D( DESCQ( CTXT_ ), 'A', ' ', N, 1, D, N, IQROW, 181 $ IQCOL ) 182 END IF 183* 184* Successively merge eigensystems of adjacent submatrices 185* into eigensystem for the corresponding larger matrix. 186* 187* while ( SUBPBS > 1 ) 188* 189 60 CONTINUE 190 IF( SUBPBS.GT.1 ) THEN 191 IM2 = SUBPBS - 2 192 DO 80 I = 0, IM2, 2 193 IF( I.EQ.0 ) THEN 194 NBL = IWORK( 2 ) 195 NBL1 = IWORK( 1 ) 196 IF( NBL1.EQ.0 ) 197 $ GO TO 70 198 ID = 1 199 MATSIZ = MIN( N, NBL*NB ) 200 N1 = NBL1*NB 201 ELSE 202 NBL = IWORK( I+2 ) - IWORK( I ) 203 NBL1 = NBL / 2 204 IF( NBL1.EQ.0 ) 205 $ GO TO 70 206 ID = IWORK( I )*NB + 1 207 MATSIZ = MIN( NB*NBL, N-ID+1 ) 208 N1 = NBL1*NB 209 END IF 210* 211* Merge lower order eigensystems (of size N1 and MATSIZ - N1) 212* into an eigensystem of size MATSIZ. 213* 214 CALL PDLAED1( MATSIZ, N1, D( ID ), ID, Q, IQ, JQ, DESCQ, 215 $ E( ID+N1-1 ), WORK, IWORK( SUBPBS+1 ), IINFO ) 216 IF( IINFO.NE.0 ) THEN 217 INFO = IINFO*( N+1 ) + ID 218 END IF 219* 220 70 CONTINUE 221 IWORK( I / 2+1 ) = IWORK( I+2 ) 222 80 CONTINUE 223 SUBPBS = SUBPBS / 2 224* 225 GO TO 60 226 END IF 227* 228* end while 229* 230 90 CONTINUE 231 RETURN 232* 233* End of PDLAED0 234* 235 END 236