|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
00001 SUBROUTINE PSPOCON( UPLO, N, A, IA, JA, DESCA, ANORM, RCOND, WORK, 00002 $ LWORK, IWORK, LIWORK, INFO ) 00003 * 00004 * -- ScaLAPACK routine (version 1.7) -- 00005 * University of Tennessee, Knoxville, Oak Ridge National Laboratory, 00006 * and University of California, Berkeley. 00007 * May 25, 2001 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER UPLO 00011 INTEGER IA, INFO, JA, LIWORK, LWORK, N 00012 REAL ANORM, RCOND 00013 * .. 00014 * .. Array Arguments .. 00015 INTEGER DESCA( * ), IWORK( * ) 00016 REAL A( * ), WORK( * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * PSPOCON estimates the reciprocal of the condition number (in the 00023 * 1-norm) of a real symmetric positive definite distributed matrix 00024 * using the Cholesky factorization A = U**T*U or A = L*L**T computed by 00025 * PSPOTRF. 00026 * 00027 * An estimate is obtained for norm(inv(A(IA:IA+N-1,JA:JA+N-1))), and 00028 * the reciprocal of the condition number is computed as 00029 * RCOND = 1 / ( norm( A(IA:IA+N-1,JA:JA+N-1) ) * 00030 * norm( inv(A(IA:IA+N-1,JA:JA+N-1)) ) ). 00031 * 00032 * Notes 00033 * ===== 00034 * 00035 * Each global data object is described by an associated description 00036 * vector. This vector stores the information required to establish 00037 * the mapping between an object element and its corresponding process 00038 * and memory location. 00039 * 00040 * Let A be a generic term for any 2D block cyclicly distributed array. 00041 * Such a global array has an associated description vector DESCA. 00042 * In the following comments, the character _ should be read as 00043 * "of the global array". 00044 * 00045 * NOTATION STORED IN EXPLANATION 00046 * --------------- -------------- -------------------------------------- 00047 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case, 00048 * DTYPE_A = 1. 00049 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating 00050 * the BLACS process grid A is distribu- 00051 * ted over. The context itself is glo- 00052 * bal, but the handle (the integer 00053 * value) may vary. 00054 * M_A (global) DESCA( M_ ) The number of rows in the global 00055 * array A. 00056 * N_A (global) DESCA( N_ ) The number of columns in the global 00057 * array A. 00058 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute 00059 * the rows of the array. 00060 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute 00061 * the columns of the array. 00062 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first 00063 * row of the array A is distributed. 00064 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the 00065 * first column of the array A is 00066 * distributed. 00067 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local 00068 * array. LLD_A >= MAX(1,LOCr(M_A)). 00069 * 00070 * Let K be the number of rows or columns of a distributed matrix, 00071 * and assume that its process grid has dimension p x q. 00072 * LOCr( K ) denotes the number of elements of K that a process 00073 * would receive if K were distributed over the p processes of its 00074 * process column. 00075 * Similarly, LOCc( K ) denotes the number of elements of K that a 00076 * process would receive if K were distributed over the q processes of 00077 * its process row. 00078 * The values of LOCr() and LOCc() may be determined via a call to the 00079 * ScaLAPACK tool function, NUMROC: 00080 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ), 00081 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ). 00082 * An upper bound for these quantities may be computed by: 00083 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A 00084 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A 00085 * 00086 * Arguments 00087 * ========= 00088 * 00089 * UPLO (global input) CHARACTER 00090 * Specifies whether the factor stored in 00091 * A(IA:IA+N-1,JA:JA+N-1) is upper or lower triangular. 00092 * = 'U': Upper triangular 00093 * = 'L': Lower triangular 00094 * 00095 * N (global input) INTEGER 00096 * The order of the distributed matrix A(IA:IA+N-1,JA:JA+N-1). 00097 * N >= 0. 00098 * 00099 * A (local input) REAL pointer into the local memory to 00100 * an array of dimension ( LLD_A, LOCc(JA+N-1) ). On entry, this 00101 * array contains the local pieces of the factors L or U from 00102 * the Cholesky factorization A(IA:IA+N-1,JA:JA+N-1) = U'*U or 00103 * L*L', as computed by PSPOTRF. 00104 * 00105 * IA (global input) INTEGER 00106 * The row index in the global array A indicating the first 00107 * row of sub( A ). 00108 * 00109 * JA (global input) INTEGER 00110 * The column index in the global array A indicating the 00111 * first column of sub( A ). 00112 * 00113 * DESCA (global and local input) INTEGER array of dimension DLEN_. 00114 * The array descriptor for the distributed matrix A. 00115 * 00116 * ANORM (global input) REAL 00117 * The 1-norm (or infinity-norm) of the symmetric distributed 00118 * matrix A(IA:IA+N-1,JA:JA+N-1). 00119 * 00120 * RCOND (global output) REAL 00121 * The reciprocal of the condition number of the distributed 00122 * matrix A(IA:IA+N-1,JA:JA+N-1), computed as 00123 * RCOND = 1 / ( norm( A(IA:IA+N-1,JA:JA+N-1) ) * 00124 * norm( inv(A(IA:IA+N-1,JA:JA+N-1)) ) ). 00125 * 00126 * WORK (local workspace/local output) REAL array, 00127 * dimension (LWORK) 00128 * On exit, WORK(1) returns the minimal and optimal LWORK. 00129 * 00130 * LWORK (local or global input) INTEGER 00131 * The dimension of the array WORK. 00132 * LWORK is local input and must be at least 00133 * LWORK >= 2*LOCr(N+MOD(IA-1,MB_A)) + 2*LOCc(N+MOD(JA-1,NB_A))+ 00134 * MAX( 2, MAX(NB_A*CEIL(NPROW-1,NPCOL),LOCc(N+MOD(JA-1,NB_A)) + 00135 * NB_A*CEIL(NPCOL-1,NPROW)) ). 00136 * 00137 * If LWORK = -1, then LWORK is global input and a workspace 00138 * query is assumed; the routine only calculates the minimum 00139 * and optimal size for all work arrays. Each of these 00140 * values is returned in the first entry of the corresponding 00141 * work array, and no error message is issued by PXERBLA. 00142 * 00143 * IWORK (local workspace/local output) INTEGER array, 00144 * dimension (LIWORK) 00145 * On exit, IWORK(1) returns the minimal and optimal LIWORK. 00146 * 00147 * LIWORK (local or global input) INTEGER 00148 * The dimension of the array IWORK. 00149 * LIWORK is local input and must be at least 00150 * LIWORK >= LOCr(N+MOD(IA-1,MB_A)). 00151 * 00152 * If LIWORK = -1, then LIWORK is global input and a workspace 00153 * query is assumed; the routine only calculates the minimum 00154 * and optimal size for all work arrays. Each of these 00155 * values is returned in the first entry of the corresponding 00156 * work array, and no error message is issued by PXERBLA. 00157 * 00158 * 00159 * INFO (global output) INTEGER 00160 * = 0: successful exit 00161 * < 0: If the i-th argument is an array and the j-entry had 00162 * an illegal value, then INFO = -(i*100+j), if the i-th 00163 * argument is a scalar and had an illegal value, then 00164 * INFO = -i. 00165 * 00166 * ===================================================================== 00167 * 00168 * .. Parameters .. 00169 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_, 00170 $ LLD_, MB_, M_, NB_, N_, RSRC_ 00171 PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1, 00172 $ CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, 00173 $ RSRC_ = 7, CSRC_ = 8, LLD_ = 9 ) 00174 REAL ONE, ZERO 00175 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00176 * .. 00177 * .. Local Scalars .. 00178 LOGICAL LQUERY, UPPER 00179 CHARACTER CBTOP, COLCTOP, NORMIN, ROWCTOP 00180 INTEGER IACOL, IAROW, ICOFF, ICTXT, IIA, IPNL, IPNU, 00181 $ IPV, IPW, IPX, IROFF, IV, IX, IXX, JJA, JV, 00182 $ JX, KASE, LIWMIN, LWMIN, MYCOL, MYROW, NP, 00183 $ NPCOL, NPROW, NPMOD, NQ, NQMOD 00184 REAL AINVNM, SCALE, SL, SU, SMLNUM 00185 REAL WMAX 00186 * .. 00187 * .. Local Arrays .. 00188 INTEGER DESCV( DLEN_ ), DESCX( DLEN_ ), IDUM1( 3 ), 00189 $ IDUM2( 3 ) 00190 * .. 00191 * .. External Subroutines .. 00192 EXTERNAL BLACS_GRIDINFO, CHK1MAT, DESCSET, INFOG2L, 00193 $ PCHK1MAT, PSAMAX, PSLATRS, PSLACON, 00194 $ PSRSCL, PB_TOPGET, PB_TOPSET, PXERBLA, SGEBR2D, 00195 $ SGEBS2D 00196 * .. 00197 * .. External Functions .. 00198 LOGICAL LSAME 00199 INTEGER ICEIL, INDXG2P, NUMROC 00200 REAL PSLAMCH 00201 EXTERNAL ICEIL, INDXG2P, LSAME, NUMROC, PSLAMCH 00202 * .. 00203 * .. Intrinsic Functions .. 00204 INTRINSIC ABS, ICHAR, MAX, MOD, REAL 00205 * .. 00206 * .. Executable Statements .. 00207 * 00208 * Get grid parameters 00209 * 00210 ICTXT = DESCA( CTXT_ ) 00211 CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL ) 00212 * 00213 * Test the input parameters 00214 * 00215 INFO = 0 00216 IF( NPROW.EQ.-1 ) THEN 00217 INFO = -(600+CTXT_) 00218 ELSE 00219 CALL CHK1MAT( N, 2, N, 2, IA, JA, DESCA, 6, INFO ) 00220 IF( INFO.EQ.0 ) THEN 00221 UPPER = LSAME( UPLO, 'U' ) 00222 IAROW = INDXG2P( IA, DESCA( MB_ ), MYROW, DESCA( RSRC_ ), 00223 $ NPROW ) 00224 IACOL = INDXG2P( JA, DESCA( NB_ ), MYCOL, DESCA( CSRC_ ), 00225 $ NPCOL ) 00226 NPMOD = NUMROC( N + MOD( IA-1, DESCA( MB_ ) ), DESCA( MB_ ), 00227 $ MYROW, IAROW, NPROW ) 00228 NQMOD = NUMROC( N + MOD( JA-1, DESCA( NB_ ) ), DESCA( NB_ ), 00229 $ MYCOL, IACOL, NPCOL ) 00230 LWMIN = 2*NPMOD + 2*NQMOD + 00231 $ MAX( 2, MAX( DESCA( NB_ )* 00232 $ MAX( 1, ICEIL( NPROW-1, NPCOL ) ), NQMOD + 00233 $ DESCA( NB_ )* 00234 $ MAX( 1, ICEIL( NPCOL-1, NPROW ) ) ) ) 00235 WORK( 1 ) = REAL( LWMIN ) 00236 LIWMIN = NPMOD 00237 IWORK( 1 ) = LIWMIN 00238 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 00239 * 00240 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00241 INFO = -1 00242 ELSE IF( ANORM.LT.ZERO ) THEN 00243 INFO = -7 00244 ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00245 INFO = -10 00246 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 00247 IWORK( 1 ) = LIWMIN 00248 INFO = -12 00249 END IF 00250 END IF 00251 * 00252 IF( UPPER ) THEN 00253 IDUM1( 1 ) = ICHAR( 'U' ) 00254 ELSE 00255 IDUM1( 1 ) = ICHAR( 'L' ) 00256 END IF 00257 IDUM2( 1 ) = 1 00258 IF( LWORK.EQ.-1 ) THEN 00259 IDUM1( 2 ) = -1 00260 ELSE 00261 IDUM1( 2 ) = 1 00262 END IF 00263 IDUM2( 2 ) = 10 00264 IF( LIWORK.EQ.-1 ) THEN 00265 IDUM1( 3 ) = -1 00266 ELSE 00267 IDUM1( 3 ) = 1 00268 END IF 00269 IDUM2( 3 ) = 12 00270 CALL PCHK1MAT( N, 2, N, 2, IA, JA, DESCA, 6, 3, IDUM1, IDUM2, 00271 $ INFO ) 00272 END IF 00273 * 00274 IF( INFO.NE.0 ) THEN 00275 CALL PXERBLA( ICTXT, 'PSPOCON', -INFO ) 00276 RETURN 00277 ELSE IF( LQUERY ) THEN 00278 RETURN 00279 END IF 00280 * 00281 * Quick return if possible 00282 * 00283 RCOND = ZERO 00284 IF( N.EQ.0 ) THEN 00285 RCOND = ONE 00286 RETURN 00287 ELSE IF( ANORM.EQ.ZERO ) THEN 00288 RETURN 00289 ELSE IF( N.EQ.1 ) THEN 00290 RCOND = ONE 00291 RETURN 00292 END IF 00293 * 00294 CALL PB_TOPGET( ICTXT, 'Combine', 'Columnwise', COLCTOP ) 00295 CALL PB_TOPGET( ICTXT, 'Combine', 'Rowwise', ROWCTOP ) 00296 CALL PB_TOPSET( ICTXT, 'Combine', 'Columnwise', '1-tree' ) 00297 CALL PB_TOPSET( ICTXT, 'Combine', 'Rowwise', '1-tree' ) 00298 * 00299 SMLNUM = PSLAMCH( ICTXT, 'Safe minimum' ) 00300 IROFF = MOD( IA-1, DESCA( MB_ ) ) 00301 ICOFF = MOD( JA-1, DESCA( NB_ ) ) 00302 CALL INFOG2L( IA, JA, DESCA, NPROW, NPCOL, MYROW, MYCOL, IIA, JJA, 00303 $ IAROW, IACOL ) 00304 NP = NUMROC( N+IROFF, DESCA( MB_ ), MYROW, IAROW, NPROW ) 00305 NQ = NUMROC( N+ICOFF, DESCA( NB_ ), MYCOL, IACOL, NPCOL ) 00306 IV = IROFF + 1 00307 IX = IV 00308 JV = ICOFF + 1 00309 JX = JV 00310 * 00311 IPX = 1 00312 IPV = IPX + NP 00313 IPNL = IPV + NP 00314 IPNU = IPNL + NQ 00315 IPW = IPNU + NQ 00316 * 00317 CALL DESCSET( DESCV, N+IROFF, 1, DESCA( MB_ ), 1, IAROW, MYCOL, 00318 $ ICTXT, MAX( 1, NP ) ) 00319 CALL DESCSET( DESCX, N+IROFF, 1, DESCA( MB_ ), 1, IAROW, MYCOL, 00320 $ ICTXT, MAX( 1, NP ) ) 00321 * 00322 * Estimate the 1-norm (or I-norm) of inv(A). 00323 * 00324 AINVNM = ZERO 00325 KASE = 0 00326 NORMIN = 'N' 00327 * 00328 10 CONTINUE 00329 CALL PSLACON( N, WORK( IPV ), IV, JV, DESCV, WORK( IPX ), IX, JX, 00330 $ DESCX, IWORK, AINVNM, KASE ) 00331 IF( KASE.NE.0 ) THEN 00332 IF( UPPER ) THEN 00333 * 00334 * Multiply by inv(U'). 00335 * 00336 DESCX( CSRC_ ) = IACOL 00337 CALL PSLATRS( 'Upper', 'Transpose', 'Non-unit', NORMIN, 00338 $ N, A, IA, JA, DESCA, WORK( IPX ), IX, JX, 00339 $ DESCX, SL, WORK( IPNL ), WORK( IPW ) ) 00340 DESCX( CSRC_ ) = MYCOL 00341 NORMIN = 'Y' 00342 * 00343 * Multiply by inv(U). 00344 * 00345 DESCX( CSRC_ ) = IACOL 00346 CALL PSLATRS( 'Upper', 'No transpose', 'Non-unit', NORMIN, 00347 $ N, A, IA, JA, DESCA, WORK( IPX ), IX, JX, 00348 $ DESCX, SU, WORK( IPNU ), WORK( IPW ) ) 00349 DESCX( CSRC_ ) = MYCOL 00350 ELSE 00351 * 00352 * Multiply by inv(L). 00353 * 00354 DESCX( CSRC_ ) = IACOL 00355 CALL PSLATRS( 'Lower', 'No transpose', 'Non-unit', NORMIN, 00356 $ N, A, IA, JA, DESCA, WORK( IPX ), IX, JX, 00357 $ DESCX, SL, WORK( IPNL ), WORK( IPW ) ) 00358 DESCX( CSRC_ ) = MYCOL 00359 NORMIN = 'Y' 00360 * 00361 * Multiply by inv(L'). 00362 * 00363 DESCX( CSRC_ ) = IACOL 00364 CALL PSLATRS( 'Lower', 'Transpose', 'Non-unit', NORMIN, 00365 $ N, A, IA, JA, DESCA, WORK( IPX ), IX, JX, 00366 $ DESCX, SU, WORK( IPNU ), WORK( IPW ) ) 00367 DESCX( CSRC_ ) = MYCOL 00368 END IF 00369 * 00370 * Multiply by 1/SCALE if doing so will not cause overflow. 00371 * 00372 SCALE = SL*SU 00373 IF( SCALE.NE.ONE ) THEN 00374 CALL PSAMAX( N, WMAX, IXX, WORK( IPX ), IX, JX, DESCX, 1 ) 00375 IF( DESCX( M_ ).EQ.1 .AND. N.EQ.1 ) THEN 00376 CALL PB_TOPGET( ICTXT, 'Broadcast', 'Columnwise', CBTOP ) 00377 IF( MYROW.EQ.IAROW ) THEN 00378 CALL SGEBS2D( ICTXT, 'Column', CBTOP, 1, 1, WMAX, 1 ) 00379 ELSE 00380 CALL SGEBR2D( ICTXT, 'Column', CBTOP, 1, 1, WMAX, 1, 00381 $ IAROW, MYCOL ) 00382 END IF 00383 END IF 00384 IF( SCALE.LT.ABS( WMAX )*SMLNUM .OR. SCALE.EQ.ZERO ) 00385 $ GO TO 20 00386 CALL PSRSCL( N, SCALE, WORK( IPX ), IX, JX, DESCX, 1 ) 00387 END IF 00388 GO TO 10 00389 END IF 00390 * 00391 * Compute the estimate of the reciprocal condition number. 00392 * 00393 IF( AINVNM.NE.ZERO ) 00394 $ RCOND = ( ONE / AINVNM ) / ANORM 00395 * 00396 20 CONTINUE 00397 * 00398 CALL PB_TOPSET( ICTXT, 'Combine', 'Columnwise', COLCTOP ) 00399 CALL PB_TOPSET( ICTXT, 'Combine', 'Rowwise', ROWCTOP ) 00400 * 00401 RETURN 00402 * 00403 * End of PSPOCON 00404 * 00405 END
