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