1 SUBROUTINE PDLAED1( N, N1, D, ID, Q, IQ, JQ, DESCQ, RHO, WORK, 2 $ IWORK, INFO ) 3* 4* -- ScaLAPACK auxiliary routine (version 1.7) -- 5* University of Tennessee, Knoxville, Oak Ridge National Laboratory, 6* and University of California, Berkeley. 7* December 31, 1998 8* 9* .. Scalar Arguments .. 10 INTEGER ID, INFO, IQ, JQ, N, N1 11 DOUBLE PRECISION RHO 12* .. 13* .. Array Arguments .. 14 INTEGER DESCQ( * ), IWORK( * ) 15 DOUBLE PRECISION D( * ), Q( * ), WORK( * ) 16* .. 17* 18* Purpose 19* ======= 20* 21* PDLAED1 computes the updated eigensystem of a diagonal 22* matrix after modification by a rank-one symmetric matrix, 23* in parallel. 24* 25* T = Q(in) ( D(in) + RHO * Z*Z' ) Q'(in) = Q(out) * D(out) * Q'(out) 26* 27* where Z = Q'u, u is a vector of length N with ones in the 28* N1 and N1 + 1 th elements and zeros elsewhere. 29* 30* The eigenvectors of the original matrix are stored in Q, and the 31* eigenvalues are in D. The algorithm consists of three stages: 32* 33* The first stage consists of deflating the size of the problem 34* when there are multiple eigenvalues or if there is a zero in 35* the Z vector. For each such occurence the dimension of the 36* secular equation problem is reduced by one. This stage is 37* performed by the routine PDLAED2. 38* 39* The second stage consists of calculating the updated 40* eigenvalues. This is done by finding the roots of the secular 41* equation via the routine SLAED4 (as called by PDLAED3). 42* This routine also calculates the eigenvectors of the current 43* problem. 44* 45* The final stage consists of computing the updated eigenvectors 46* directly using the updated eigenvalues. The eigenvectors for 47* the current problem are multiplied with the eigenvectors from 48* the overall problem. 49* 50* Arguments 51* ========= 52* 53* N (global input) INTEGER 54* The order of the tridiagonal matrix T. N >= 0. 55* 56* 57* N1 (input) INTEGER 58* The location of the last eigenvalue in the leading 59* sub-matrix. 60* min(1,N) <= N1 <= N. 61* 62* D (global input/output) DOUBLE PRECISION array, dimension (N) 63* On entry,the eigenvalues of the rank-1-perturbed matrix. 64* On exit, the eigenvalues of the repaired matrix. 65* 66* ID (global input) INTEGER 67* Q's global row/col index, which points to the beginning 68* of the submatrix which is to be operated on. 69* 70* Q (local output) DOUBLE PRECISION array, 71* global dimension (N, N), 72* local dimension ( LLD_Q, LOCc(JQ+N-1)) 73* Q contains the orthonormal eigenvectors of the symmetric 74* tridiagonal matrix. 75* 76* IQ (global input) INTEGER 77* Q's global row index, which points to the beginning of the 78* submatrix which is to be operated on. 79* 80* JQ (global input) INTEGER 81* Q's global column index, which points to the beginning of 82* the submatrix which is to be operated on. 83* 84* DESCQ (global and local input) INTEGER array of dimension DLEN_. 85* The array descriptor for the distributed matrix Z. 86* 87* RHO (input) DOUBLE PRECISION 88* The subdiagonal entry used to create the rank-1 modification. 89* 90* WORK (local workspace/output) DOUBLE PRECISION array, 91* dimension 6*N + 2*NP*NQ 92* 93* IWORK (local workspace/output) INTEGER array, 94* dimension 7*N + 8*NPCOL + 2 95* 96* INFO (global output) INTEGER 97* = 0: successful exit 98* < 0: If the i-th argument is an array and the j-entry had 99* an illegal value, then INFO = -(i*100+j), if the i-th 100* argument is a scalar and had an illegal value, then 101* INFO = -i. 102* > 0: The algorithm failed to compute the ith eigenvalue. 103* 104* ===================================================================== 105* 106* .. Parameters .. 107* 108 INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_, 109 $ MB_, NB_, RSRC_, CSRC_, LLD_ 110 PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1, 111 $ CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, 112 $ RSRC_ = 7, CSRC_ = 8, LLD_ = 9 ) 113 DOUBLE PRECISION ZERO, ONE 114 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 115* .. 116* .. Local Scalars .. 117 INTEGER COL, COLTYP, IBUF, ICTOT, ICTXT, IDLMDA, IIQ, 118 $ INDCOL, INDROW, INDX, INDXC, INDXP, INDXR, INQ, 119 $ IPQ, IPQ2, IPSM, IPU, IPWORK, IQ1, IQ2, IQCOL, 120 $ IQQ, IQROW, IW, IZ, J, JC, JJ2C, JJC, JJQ, JNQ, 121 $ K, LDQ, LDQ2, LDU, MYCOL, MYROW, NB, NN, NN1, 122 $ NN2, NP, NPCOL, NPROW, NQ 123* .. 124* .. Local Arrays .. 125 INTEGER DESCQ2( DLEN_ ), DESCU( DLEN_ ) 126* .. 127* .. External Functions .. 128 INTEGER NUMROC 129 EXTERNAL NUMROC 130* .. 131* .. External Subroutines .. 132 EXTERNAL BLACS_GRIDINFO, DCOPY, DESCINIT, INFOG1L, 133 $ INFOG2L, PDGEMM, PDLAED2, PDLAED3, PDLAEDZ, 134 $ PDLASET, PXERBLA 135* .. 136* .. Intrinsic Functions .. 137 INTRINSIC MAX, MIN 138* .. 139* .. Executable Statements .. 140* 141* This is just to keep ftnchek and toolpack/1 happy 142 IF( BLOCK_CYCLIC_2D*CSRC_*CTXT_*DLEN_*DTYPE_*LLD_*MB_*M_*NB_*N_* 143 $ RSRC_.LT.0 )RETURN 144* 145* 146* Test the input parameters. 147* 148 CALL BLACS_GRIDINFO( DESCQ( CTXT_ ), NPROW, NPCOL, MYROW, MYCOL ) 149 INFO = 0 150 IF( NPROW.EQ.-1 ) THEN 151 INFO = -( 600+CTXT_ ) 152 ELSE IF( N.LT.0 ) THEN 153 INFO = -1 154 ELSE IF( ID.GT.DESCQ( N_ ) ) THEN 155 INFO = -4 156 ELSE IF( N1.GE.N ) THEN 157 INFO = -2 158 END IF 159 IF( INFO.NE.0 ) THEN 160 CALL PXERBLA( DESCQ( CTXT_ ), 'PDLAED1', -INFO ) 161 RETURN 162 END IF 163* 164* Quick return if possible 165* 166 IF( N.EQ.0 ) 167 $ RETURN 168* 169* The following values are integer pointers which indicate 170* the portion of the workspace used by a particular array 171* in PDLAED2 and PDLAED3. 172* 173 ICTXT = DESCQ( CTXT_ ) 174 NB = DESCQ( NB_ ) 175 LDQ = DESCQ( LLD_ ) 176* 177 CALL INFOG2L( IQ-1+ID, JQ-1+ID, DESCQ, NPROW, NPCOL, MYROW, MYCOL, 178 $ IIQ, JJQ, IQROW, IQCOL ) 179* 180 NP = NUMROC( N, DESCQ( MB_ ), MYROW, IQROW, NPROW ) 181 NQ = NUMROC( N, DESCQ( NB_ ), MYCOL, IQCOL, NPCOL ) 182* 183 LDQ2 = MAX( NP, 1 ) 184 LDU = LDQ2 185* 186 IZ = 1 187 IDLMDA = IZ + N 188 IW = IDLMDA + N 189 IPQ2 = IW + N 190 IPU = IPQ2 + LDQ2*NQ 191 IBUF = IPU + LDU*NQ 192* (IBUF est de taille 3*N au maximum) 193* 194 ICTOT = 1 195 IPSM = ICTOT + NPCOL*4 196 INDX = IPSM + NPCOL*4 197 INDXC = INDX + N 198 INDXP = INDXC + N 199 INDCOL = INDXP + N 200 COLTYP = INDCOL + N 201 INDROW = COLTYP + N 202 INDXR = INDROW + N 203* 204 CALL DESCINIT( DESCQ2, N, N, NB, NB, IQROW, IQCOL, ICTXT, LDQ2, 205 $ INFO ) 206 CALL DESCINIT( DESCU, N, N, NB, NB, IQROW, IQCOL, ICTXT, LDU, 207 $ INFO ) 208* 209* Form the z-vector which consists of the last row of Q_1 and the 210* first row of Q_2. 211* 212 IPWORK = IDLMDA 213 CALL PDLAEDZ( N, N1, ID, Q, IQ, JQ, LDQ, DESCQ, WORK( IZ ), 214 $ WORK( IPWORK ) ) 215* 216* Deflate eigenvalues. 217* 218 IPQ = IIQ + ( JJQ-1 )*LDQ 219 CALL PDLAED2( ICTXT, K, N, N1, NB, D, IQROW, IQCOL, Q( IPQ ), LDQ, 220 $ RHO, WORK( IZ ), WORK( IW ), WORK( IDLMDA ), 221 $ WORK( IPQ2 ), LDQ2, WORK( IBUF ), IWORK( ICTOT ), 222 $ IWORK( IPSM ), NPCOL, IWORK( INDX ), IWORK( INDXC ), 223 $ IWORK( INDXP ), IWORK( INDCOL ), IWORK( COLTYP ), 224 $ NN, NN1, NN2, IQ1, IQ2 ) 225* 226* 227* Solve Secular Equation. 228* 229 IF( K.NE.0 ) THEN 230 CALL PDLASET( 'A', N, N, ZERO, ONE, WORK( IPU ), 1, 1, DESCU ) 231 CALL PDLAED3( ICTXT, K, N, NB, D, IQROW, IQCOL, RHO, 232 $ WORK( IDLMDA ), WORK( IW ), WORK( IZ ), 233 $ WORK( IPU ), LDQ2, WORK( IBUF ), IWORK( INDX ), 234 $ IWORK( INDCOL ), IWORK( INDROW ), IWORK( INDXR ), 235 $ IWORK( INDXC ), IWORK( ICTOT ), NPCOL, INFO ) 236* 237* Compute the updated eigenvectors. 238* 239 IQQ = MIN( IQ1, IQ2 ) 240 IF( NN1.GT.0 ) THEN 241 INQ = IQ - 1 + ID 242 JNQ = JQ - 1 + ID + IQQ - 1 243 CALL PDGEMM( 'N', 'N', N1, NN, NN1, ONE, WORK( IPQ2 ), 1, 244 $ IQ1, DESCQ2, WORK( IPU ), IQ1, IQQ, DESCU, 245 $ ZERO, Q, INQ, JNQ, DESCQ ) 246 END IF 247 IF( NN2.GT.0 ) THEN 248 INQ = IQ - 1 + ID + N1 249 JNQ = JQ - 1 + ID + IQQ - 1 250 CALL PDGEMM( 'N', 'N', N-N1, NN, NN2, ONE, WORK( IPQ2 ), 251 $ N1+1, IQ2, DESCQ2, WORK( IPU ), IQ2, IQQ, 252 $ DESCU, ZERO, Q, INQ, JNQ, DESCQ ) 253 END IF 254* 255 DO 10 J = K + 1, N 256 JC = IWORK( INDX+J-1 ) 257 CALL INFOG1L( JQ-1+JC, NB, NPCOL, MYCOL, IQCOL, JJC, COL ) 258 CALL INFOG1L( JC, NB, NPCOL, MYCOL, IQCOL, JJ2C, COL ) 259 IF( MYCOL.EQ.COL ) THEN 260 IQ2 = IPQ2 + ( JJ2C-1 )*LDQ2 261 INQ = IPQ + ( JJC-1 )*LDQ 262 CALL DCOPY( NP, WORK( IQ2 ), 1, Q( INQ ), 1 ) 263 END IF 264 10 CONTINUE 265 END IF 266* 267 20 CONTINUE 268 RETURN 269* 270* End of PDLAED1 271* 272 END 273