1 SUBROUTINE PDSTEDC( COMPZ, N, D, E, Q, IQ, JQ, DESCQ, WORK, LWORK, 2 $ IWORK, LIWORK, INFO ) 3* 4* -- ScaLAPACK routine (version 1.7) -- 5* University of Tennessee, Knoxville, Oak Ridge National Laboratory, 6* and University of California, Berkeley. 7* March 13, 2000 8* 9* .. Scalar Arguments .. 10 CHARACTER COMPZ 11 INTEGER INFO, IQ, JQ, LIWORK, LWORK, N 12* .. 13* .. Array Arguments .. 14 INTEGER DESCQ( * ), IWORK( * ) 15 DOUBLE PRECISION D( * ), E( * ), Q( * ), WORK( * ) 16* .. 17* 18* Purpose 19* ======= 20* PDSTEDC computes all eigenvalues and eigenvectors of a 21* symmetric tridiagonal matrix in parallel, using the divide and 22* conquer algorithm. 23* 24* This code makes very mild assumptions about floating point 25* arithmetic. It will work on machines with a guard digit in 26* add/subtract, or on those binary machines without guard digits 27* which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. 28* It could conceivably fail on hexadecimal or decimal machines 29* without guard digits, but we know of none. See DLAED3 for details. 30* 31* Arguments 32* ========= 33* 34* COMPZ (input) CHARACTER*1 35* = 'N': Compute eigenvalues only. (NOT IMPLEMENTED YET) 36* = 'I': Compute eigenvectors of tridiagonal matrix also. 37* = 'V': Compute eigenvectors of original dense symmetric 38* matrix also. On entry, Z contains the orthogonal 39* matrix used to reduce the original matrix to 40* tridiagonal form. (NOT IMPLEMENTED YET) 41* 42* N (global input) INTEGER 43* The order of the tridiagonal matrix T. N >= 0. 44* 45* D (global input/output) DOUBLE PRECISION array, dimension (N) 46* On entry, the diagonal elements of the tridiagonal matrix. 47* On exit, if INFO = 0, the eigenvalues in descending order. 48* 49* E (global input/output) DOUBLE PRECISION array, dimension (N-1) 50* On entry, the subdiagonal elements of the tridiagonal matrix. 51* On exit, E has been destroyed. 52* 53* Q (local output) DOUBLE PRECISION array, 54* local dimension ( LLD_Q, LOCc(JQ+N-1)) 55* Q contains the orthonormal eigenvectors of the symmetric 56* tridiagonal matrix. 57* On output, Q is distributed across the P processes in block 58* cyclic format. 59* 60* IQ (global input) INTEGER 61* Q's global row index, which points to the beginning of the 62* submatrix which is to be operated on. 63* 64* JQ (global input) INTEGER 65* Q's global column index, which points to the beginning of 66* the submatrix which is to be operated on. 67* 68* DESCQ (global and local input) INTEGER array of dimension DLEN_. 69* The array descriptor for the distributed matrix Z. 70* 71* 72* WORK (local workspace/output) DOUBLE PRECISION array, 73* dimension (LWORK) 74* On output, WORK(1) returns the workspace needed. 75* 76* LWORK (local input/output) INTEGER, 77* the dimension of the array WORK. 78* LWORK = 6*N + 2*NP*NQ 79* NP = NUMROC( N, NB, MYROW, DESCQ( RSRC_ ), NPROW ) 80* NQ = NUMROC( N, NB, MYCOL, DESCQ( CSRC_ ), NPCOL ) 81* 82* If LWORK = -1, the LWORK is global input and a workspace 83* query is assumed; the routine only calculates the minimum 84* size for the WORK array. The required workspace is returned 85* as the first element of WORK and no error message is issued 86* by PXERBLA. 87* 88* IWORK (local workspace/output) INTEGER array, dimension (LIWORK) 89* On exit, if LIWORK > 0, IWORK(1) returns the optimal LIWORK. 90* 91* LIWORK (input) INTEGER 92* The dimension of the array IWORK. 93* LIWORK = 2 + 7*N + 8*NPCOL 94* 95* INFO (global output) INTEGER 96* = 0: successful exit 97* < 0: If the i-th argument is an array and the j-entry had 98* an illegal value, then INFO = -(i*100+j), if the i-th 99* argument is a scalar and had an illegal value, then 100* INFO = -i. 101* > 0: The algorithm failed to compute the INFO/(N+1) th 102* eigenvalue while working on the submatrix lying in 103* global rows and columns mod(INFO,N+1). 104* 105* Further Details 106* ======= ======= 107* 108* Contributed by Francoise Tisseur, University of Manchester. 109* 110* Reference: F. Tisseur and J. Dongarra, "A Parallel Divide and 111* Conquer Algorithm for the Symmetric Eigenvalue Problem 112* on Distributed Memory Architectures", 113* SIAM J. Sci. Comput., 6:20 (1999), pp. 2223--2236. 114* (see also LAPACK Working Note 132) 115* http://www.netlib.org/lapack/lawns/lawn132.ps 116* 117* ===================================================================== 118* 119* .. Parameters .. 120 INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_, 121 $ MB_, NB_, RSRC_, CSRC_, LLD_ 122 PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1, 123 $ CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, 124 $ RSRC_ = 7, CSRC_ = 8, LLD_ = 9 ) 125 DOUBLE PRECISION ZERO, ONE 126 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 127* .. 128* .. Local Scalars .. 129 LOGICAL LQUERY 130 INTEGER ICOFFQ, IIQ, IPQ, IQCOL, IQROW, IROFFQ, JJQ, 131 $ LDQ, LIWMIN, LWMIN, MYCOL, MYROW, NB, NP, 132 $ NPCOL, NPROW, NQ 133 DOUBLE PRECISION ORGNRM 134* .. 135* .. External Functions .. 136 LOGICAL LSAME 137 INTEGER INDXG2P, NUMROC 138 DOUBLE PRECISION DLANST 139 EXTERNAL INDXG2P, LSAME, NUMROC, DLANST 140* .. 141* .. External Subroutines .. 142 EXTERNAL BLACS_GRIDINFO, CHK1MAT, DLASCL, DSTEDC, 143 $ INFOG2L, PDLAED0, PDLASRT, PXERBLA 144* .. 145* .. Intrinsic Functions .. 146 INTRINSIC DBLE, MOD 147* .. 148* .. Executable Statements .. 149* 150* This is just to keep ftnchek and toolpack/1 happy 151 IF( BLOCK_CYCLIC_2D*CSRC_*CTXT_*DLEN_*DTYPE_*LLD_*MB_*M_*NB_*N_* 152 $ RSRC_.LT.0 )RETURN 153* 154* Test the input parameters. 155* 156 CALL BLACS_GRIDINFO( DESCQ( CTXT_ ), NPROW, NPCOL, MYROW, MYCOL ) 157 LDQ = DESCQ( LLD_ ) 158 NB = DESCQ( NB_ ) 159 NP = NUMROC( N, NB, MYROW, DESCQ( RSRC_ ), NPROW ) 160 NQ = NUMROC( N, NB, MYCOL, DESCQ( CSRC_ ), NPCOL ) 161 INFO = 0 162 IF( NPROW.EQ.-1 ) THEN 163 INFO = -( 600+CTXT_ ) 164 ELSE 165 CALL CHK1MAT( N, 2, N, 2, IQ, JQ, DESCQ, 8, INFO ) 166 IF( INFO.EQ.0 ) THEN 167 NB = DESCQ( NB_ ) 168 IROFFQ = MOD( IQ-1, DESCQ( MB_ ) ) 169 ICOFFQ = MOD( JQ-1, DESCQ( NB_ ) ) 170 IQROW = INDXG2P( IQ, NB, MYROW, DESCQ( RSRC_ ), NPROW ) 171 IQCOL = INDXG2P( JQ, NB, MYCOL, DESCQ( CSRC_ ), NPCOL ) 172 LWMIN = 6*N + 2*NP*NQ 173 LIWMIN = 2 + 7*N + 8*NPCOL 174* 175 WORK( 1 ) = DBLE( LWMIN ) 176 IWORK( 1 ) = LIWMIN 177 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 178 IF( .NOT.LSAME( COMPZ, 'I' ) ) THEN 179 INFO = -1 180 ELSE IF( N.LT.0 ) THEN 181 INFO = -2 182 ELSE IF( IROFFQ.NE.ICOFFQ .OR. ICOFFQ.NE.0 ) THEN 183 INFO = -5 184 ELSE IF( DESCQ( MB_ ).NE.DESCQ( NB_ ) ) THEN 185 INFO = -( 700+NB_ ) 186 ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 187 INFO = -10 188 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 189 INFO = -12 190 END IF 191 END IF 192 END IF 193 IF( INFO.NE.0 ) THEN 194 CALL PXERBLA( DESCQ( CTXT_ ), 'PDSTEDC', -INFO ) 195 RETURN 196 ELSE IF( LQUERY ) THEN 197 RETURN 198 END IF 199* 200* Quick return 201* 202 IF( N.EQ.0 ) 203 $ GO TO 10 204 CALL INFOG2L( IQ, JQ, DESCQ, NPROW, NPCOL, MYROW, MYCOL, IIQ, JJQ, 205 $ IQROW, IQCOL ) 206 IF( N.EQ.1 ) THEN 207 IF( MYROW.EQ.IQROW .AND. MYCOL.EQ.IQCOL ) 208 $ Q( 1 ) = ONE 209 GO TO 10 210 END IF 211* 212* If N is smaller than the minimum divide size NB, then 213* solve the problem with the serial divide and conquer 214* code locally. 215* 216 IF( N.LE.NB ) THEN 217 IF( ( MYROW.EQ.IQROW ) .AND. ( MYCOL.EQ.IQCOL ) ) THEN 218 IPQ = IIQ + ( JJQ-1 )*LDQ 219 CALL DSTEDC( 'I', N, D, E, Q( IPQ ), LDQ, WORK, LWORK, 220 $ IWORK, LIWORK, INFO ) 221 IF( INFO.NE.0 ) THEN 222 INFO = ( N+1 ) + N 223 GO TO 10 224 END IF 225 END IF 226 GO TO 10 227 END IF 228* 229* If P=NPROW*NPCOL=1, solve the problem with DSTEDC. 230* 231 IF( NPCOL*NPROW.EQ.1 ) THEN 232 IPQ = IIQ + ( JJQ-1 )*LDQ 233 CALL DSTEDC( 'I', N, D, E, Q( IPQ ), LDQ, WORK, LWORK, IWORK, 234 $ LIWORK, INFO ) 235 GO TO 10 236 END IF 237* 238* Scale matrix to allowable range, if necessary. 239* 240 ORGNRM = DLANST( 'M', N, D, E ) 241 IF( ORGNRM.NE.ZERO ) THEN 242 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO ) 243 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N-1, 1, E, N-1, INFO ) 244 END IF 245* 246 CALL PDLAED0( N, D, E, Q, IQ, JQ, DESCQ, WORK, IWORK, INFO ) 247* 248* Sort eigenvalues and corresponding eigenvectors 249* 250 CALL PDLASRT( 'I', N, D, Q, IQ, JQ, DESCQ, WORK, LWORK, IWORK, 251 $ LIWORK, INFO ) 252* 253* Scale back. 254* 255 IF( ORGNRM.NE.ZERO ) 256 $ CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO ) 257* 258 10 CONTINUE 259* 260 IF( LWORK.GT.0 ) 261 $ WORK( 1 ) = DBLE( LWMIN ) 262 IF( LIWORK.GT.0 ) 263 $ IWORK( 1 ) = LIWMIN 264 RETURN 265* 266* End of PDSTEDC 267* 268 END 269