1 SUBROUTINE PSTRTRS( UPLO, TRANS, DIAG, N, NRHS, A, IA, JA, DESCA, 2 $ B, IB, JB, DESCB, 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* May 1, 1997 8* 9* .. Scalar Arguments .. 10 CHARACTER DIAG, TRANS, UPLO 11 INTEGER IA, IB, INFO, JA, JB, N, NRHS 12* .. 13* .. Array Arguments .. 14 INTEGER DESCA( * ), DESCB( * ) 15 REAL A( * ), B( * ) 16* .. 17* 18* Purpose 19* ======= 20* 21* PSTRTRS solves a triangular system of the form 22* 23* sub( A ) * X = sub( B ) or sub( A )**T * X = sub( B ), 24* 25* where sub( A ) denotes A(IA:IA+N-1,JA:JA+N-1) and is a triangular 26* distributed matrix of order N, and B(IB:IB+N-1,JB:JB+NRHS-1) is an 27* N-by-NRHS distributed matrix denoted by sub( B ). A check is made 28* to verify that sub( A ) is nonsingular. 29* 30* Notes 31* ===== 32* 33* Each global data object is described by an associated description 34* vector. This vector stores the information required to establish 35* the mapping between an object element and its corresponding process 36* and memory location. 37* 38* Let A be a generic term for any 2D block cyclicly distributed array. 39* Such a global array has an associated description vector DESCA. 40* In the following comments, the character _ should be read as 41* "of the global array". 42* 43* NOTATION STORED IN EXPLANATION 44* --------------- -------------- -------------------------------------- 45* DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case, 46* DTYPE_A = 1. 47* CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating 48* the BLACS process grid A is distribu- 49* ted over. The context itself is glo- 50* bal, but the handle (the integer 51* value) may vary. 52* M_A (global) DESCA( M_ ) The number of rows in the global 53* array A. 54* N_A (global) DESCA( N_ ) The number of columns in the global 55* array A. 56* MB_A (global) DESCA( MB_ ) The blocking factor used to distribute 57* the rows of the array. 58* NB_A (global) DESCA( NB_ ) The blocking factor used to distribute 59* the columns of the array. 60* RSRC_A (global) DESCA( RSRC_ ) The process row over which the first 61* row of the array A is distributed. 62* CSRC_A (global) DESCA( CSRC_ ) The process column over which the 63* first column of the array A is 64* distributed. 65* LLD_A (local) DESCA( LLD_ ) The leading dimension of the local 66* array. LLD_A >= MAX(1,LOCr(M_A)). 67* 68* Let K be the number of rows or columns of a distributed matrix, 69* and assume that its process grid has dimension p x q. 70* LOCr( K ) denotes the number of elements of K that a process 71* would receive if K were distributed over the p processes of its 72* process column. 73* Similarly, LOCc( K ) denotes the number of elements of K that a 74* process would receive if K were distributed over the q processes of 75* its process row. 76* The values of LOCr() and LOCc() may be determined via a call to the 77* ScaLAPACK tool function, NUMROC: 78* LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ), 79* LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ). 80* An upper bound for these quantities may be computed by: 81* LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A 82* LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A 83* 84* Arguments 85* ========= 86* 87* UPLO (global input) CHARACTER 88* = 'U': sub( A ) is upper triangular; 89* = 'L': sub( A ) is lower triangular. 90* 91* TRANS (global input) CHARACTER 92* Specifies the form of the system of equations: 93* = 'N': Solve sub( A ) * X = sub( B ) (No transpose) 94* = 'T': Solve sub( A )**T * X = sub( B ) (Transpose) 95* = 'C': Solve sub( A )**T * X = sub( B ) (Transpose) 96* 97* DIAG (global input) CHARACTER 98* = 'N': sub( A ) is non-unit triangular; 99* = 'U': sub( A ) is unit triangular. 100* 101* N (global input) INTEGER 102* The number of rows and columns to be operated on i.e the 103* order of the distributed submatrix sub( A ). N >= 0. 104* 105* NRHS (global input) INTEGER 106* The number of right hand sides, i.e., the number of columns 107* of the distributed matrix sub( B ). NRHS >= 0. 108* 109* A (local input) REAL pointer into the local memory 110* to an array of dimension (LLD_A,LOCc(JA+N-1) ). This array 111* contains the local pieces of the distributed triangular 112* matrix sub( A ). If UPLO = 'U', the leading N-by-N upper 113* triangular part of sub( A ) contains the upper triangular 114* matrix, and the strictly lower triangular part of sub( A ) 115* is not referenced. If UPLO = 'L', the leading N-by-N lower 116* triangular part of sub( A ) contains the lower triangular 117* matrix, and the strictly upper triangular part of sub( A ) 118* is not referenced. If DIAG = 'U', the diagonal elements of 119* sub( A ) are also not referenced and are assumed to be 1. 120* 121* IA (global input) INTEGER 122* The row index in the global array A indicating the first 123* row of sub( A ). 124* 125* JA (global input) INTEGER 126* The column index in the global array A indicating the 127* first column of sub( A ). 128* 129* DESCA (global and local input) INTEGER array of dimension DLEN_. 130* The array descriptor for the distributed matrix A. 131* 132* B (local input/local output) REAL pointer into the 133* local memory to an array of dimension 134* (LLD_B,LOCc(JB+NRHS-1)). On entry, this array contains the 135* local pieces of the right hand side distributed matrix 136* sub( B ). On exit, if INFO = 0, sub( B ) is overwritten by 137* the solution matrix X. 138* 139* IB (global input) INTEGER 140* The row index in the global array B indicating the first 141* row of sub( B ). 142* 143* JB (global input) INTEGER 144* The column index in the global array B indicating the 145* first column of sub( B ). 146* 147* DESCB (global and local input) INTEGER array of dimension DLEN_. 148* The array descriptor for the distributed matrix B. 149* 150* INFO (output) INTEGER 151* = 0: successful exit 152* < 0: If the i-th argument is an array and the j-entry had 153* an illegal value, then INFO = -(i*100+j), if the i-th 154* argument is a scalar and had an illegal value, then 155* INFO = -i. 156* > 0: If INFO = i, the i-th diagonal element of sub( A ) is 157* zero, indicating that the submatrix is singular and the 158* solutions X have not been computed. 159* 160* ===================================================================== 161* 162* .. Parameters .. 163 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_, 164 $ LLD_, MB_, M_, NB_, N_, RSRC_ 165 PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1, 166 $ CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, 167 $ RSRC_ = 7, CSRC_ = 8, LLD_ = 9 ) 168 REAL ZERO, ONE 169 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 170* .. 171* .. Local Scalars .. 172 LOGICAL NOTRAN, NOUNIT, UPPER 173 INTEGER I, IAROW, IBROW, ICOFFA, ICTXT, ICURCOL, 174 $ ICURROW, IROFFA, IROFFB, IDUM, II, IOFFA, J, 175 $ JBLK, JJ, JN, LDA, LL, MYCOL, MYROW, NPCOL, 176 $ NPROW 177* .. 178* .. Local Arrays .. 179 INTEGER IDUM1( 3 ), IDUM2( 3 ) 180* .. 181* .. External Subroutines .. 182 EXTERNAL BLACS_GRIDINFO, CHK1MAT, IGAMX2D, INFOG2L, 183 $ PCHK2MAT, PSTRSM, PXERBLA 184* .. 185* .. External Functions .. 186 LOGICAL LSAME 187 INTEGER ICEIL, INDXG2P 188 EXTERNAL ICEIL, INDXG2P, LSAME 189* .. 190* .. Intrinsic Functions .. 191 INTRINSIC ICHAR, MIN, MOD 192* .. 193* .. Executable Statements .. 194* 195* Get grid parameters 196* 197 ICTXT = DESCA( CTXT_ ) 198 CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL ) 199* 200* Test input parameters 201* 202 INFO = 0 203 IF( NPROW.EQ.-1 ) THEN 204 INFO = -907 205 ELSE 206 UPPER = LSAME( UPLO, 'U' ) 207 NOUNIT = LSAME( DIAG, 'N' ) 208 NOTRAN = LSAME( TRANS, 'N' ) 209* 210 CALL CHK1MAT( N, 4, N, 4, IA, JA, DESCA, 9, INFO ) 211 CALL CHK1MAT( N, 4, NRHS, 5, IB, JB, DESCB, 13, INFO ) 212 IF( INFO.EQ.0 ) THEN 213 IROFFA = MOD( IA-1, DESCA( MB_ ) ) 214 ICOFFA = MOD( JA-1, DESCA( NB_ ) ) 215 IROFFB = MOD( IB-1, DESCB( MB_ ) ) 216 IAROW = INDXG2P( IA, DESCA( MB_ ), MYROW, DESCA( RSRC_ ), 217 $ NPROW ) 218 IBROW = INDXG2P( IB, DESCB( MB_ ), MYROW, DESCB( RSRC_ ), 219 $ NPROW ) 220 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 221 INFO = -1 222 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. 223 $ .NOT.LSAME( TRANS, 'C' ) ) THEN 224 INFO = -2 225 ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN 226 INFO = -3 227 ELSE IF( IROFFA.NE.ICOFFA .OR. IROFFA.NE.0 ) THEN 228 INFO = -8 229 ELSE IF( IROFFA.NE.IROFFB .OR. IAROW.NE.IBROW ) THEN 230 INFO = -11 231 ELSE IF( DESCA( MB_ ).NE.DESCA( NB_ ) ) THEN 232 INFO = -904 233 ELSE IF( DESCB( MB_ ).NE.DESCA( NB_ ) ) THEN 234 INFO = -1304 235 END IF 236 END IF 237* 238 IF( UPPER ) THEN 239 IDUM1( 1 ) = ICHAR( 'U' ) 240 ELSE 241 IDUM1( 1 ) = ICHAR( 'L' ) 242 END IF 243 IDUM2( 1 ) = 1 244 IF( NOTRAN ) THEN 245 IDUM1( 2 ) = ICHAR( 'N' ) 246 ELSE IF( LSAME( TRANS, 'T' ) ) THEN 247 IDUM1( 2 ) = ICHAR( 'T' ) 248 ELSE IF( LSAME( TRANS, 'C' ) ) THEN 249 IDUM1( 2 ) = ICHAR( 'C' ) 250 END IF 251 IDUM2( 2 ) = 2 252 IF( NOUNIT ) THEN 253 IDUM1( 3 ) = ICHAR( 'N' ) 254 ELSE 255 IDUM1( 3 ) = ICHAR( 'D' ) 256 END IF 257 IDUM2( 3 ) = 3 258 CALL PCHK2MAT( N, 4, N, 4, IA, JA, DESCA, 9, N, 4, NRHS, 5, 259 $ IB, JB, DESCB, 13, 3, IDUM1, IDUM2, INFO ) 260 END IF 261* 262 IF( INFO.NE.0 ) THEN 263 CALL PXERBLA( ICTXT, 'PSTRTRS', -INFO ) 264 RETURN 265 END IF 266* 267* Quick return if possible 268* 269 IF( N.EQ.0 .OR. NRHS.EQ.0 ) 270 $ RETURN 271* 272* Check for singularity if non-unit. 273* 274 IF( NOUNIT ) THEN 275 CALL INFOG2L( IA, JA, DESCA, NPROW, NPCOL, MYROW, MYCOL, 276 $ II, JJ, ICURROW, ICURCOL ) 277 JN = MIN( ICEIL( JA, DESCA( NB_ ) ) * DESCA( NB_ ), JA+N-1 ) 278 LDA = DESCA( LLD_ ) 279 IOFFA = II + ( JJ - 1 ) * LDA 280* 281* Handle first block separately 282* 283 JBLK = JN-JA+1 284 IF( MYROW.EQ.ICURROW .AND. MYCOL.EQ.ICURCOL ) THEN 285 LL = IOFFA 286 DO 10 I = 0, JBLK-1 287 IF( A( LL ).EQ.ZERO .AND. INFO.EQ.0 ) 288 $ INFO = I + 1 289 LL = IOFFA + LDA + 1 290 10 CONTINUE 291 END IF 292 IF( MYROW.EQ.ICURROW ) 293 $ IOFFA = IOFFA + JBLK 294 IF( MYCOL.EQ.ICURCOL ) 295 $ IOFFA = IOFFA + JBLK*LDA 296 ICURROW = MOD( ICURROW+1, NPROW ) 297 ICURCOL = MOD( ICURCOL+1, NPCOL ) 298* 299* Loop over remaining blocks of columns 300* 301 DO 30 J = JN+1, JA+N-1, DESCA( NB_ ) 302 JBLK = MIN( JA+N-J, DESCA( NB_ ) ) 303 IF( MYROW.EQ.ICURROW .AND. MYCOL.EQ.ICURCOL ) THEN 304 LL = IOFFA 305 DO 20 I = 0, JBLK-1 306 IF( A( LL ).EQ.ZERO .AND. INFO.EQ.0 ) 307 $ INFO = J + I - JA + 1 308 LL = IOFFA + LDA + 1 309 20 CONTINUE 310 END IF 311 IF( MYROW.EQ.ICURROW ) 312 $ IOFFA = IOFFA + JBLK 313 IF( MYCOL.EQ.ICURCOL ) 314 $ IOFFA = IOFFA + JBLK*LDA 315 ICURROW = MOD( ICURROW+1, NPROW ) 316 ICURCOL = MOD( ICURCOL+1, NPCOL ) 317 30 CONTINUE 318 CALL IGAMX2D( ICTXT, 'All', ' ', 1, 1, INFO, 1, IDUM, IDUM, 319 $ -1, -1, MYCOL ) 320 IF( INFO.NE.0 ) 321 $ RETURN 322 END IF 323* 324* Solve A * x = b or A' * x = b. 325* 326 CALL PSTRSM( 'Left', UPLO, TRANS, DIAG, N, NRHS, ONE, A, IA, JA, 327 $ DESCA, B, IB, JB, DESCB ) 328* 329 RETURN 330* 331* End of PSTRTRS 332* 333 END 334