1 SUBROUTINE PCLAPV2( DIREC, ROWCOL, M, N, A, IA, JA, DESCA, IPIV, 2 $ IP, JP, DESCIP ) 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 DIREC, ROWCOL 11 INTEGER IA, IP, JA, JP, M, N 12* .. 13* .. Array Arguments .. 14 INTEGER DESCA( * ), DESCIP( * ), IPIV( * ) 15 COMPLEX A( * ) 16* .. 17* 18* Purpose 19* ======= 20* 21* PCLAPV2 applies either P (permutation matrix indicated by IPIV) 22* or inv( P ) to a M-by-N distributed matrix sub( A ) denoting 23* A(IA:IA+M-1,JA:JA+N-1), resulting in row or column pivoting. The 24* pivot vector should be aligned with the distributed matrix A. For 25* pivoting the rows of sub( A ), IPIV should be distributed along a 26* process column and replicated over all process rows. Similarly, 27* IPIV should be distributed along a process row and replicated over 28* all process columns for column pivoting. 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* DIREC (global input) CHARACTER 88* Specifies in which order the permutation is applied: 89* = 'F' (Forward) Applies pivots Forward from top of matrix. 90* Computes P * sub( A ); 91* = 'B' (Backward) Applies pivots Backward from bottom of 92* matrix. Computes inv( P ) * sub( A ). 93* 94* ROWCOL (global input) CHARACTER 95* Specifies if the rows or columns are to be permuted: 96* = 'R' Rows will be permuted, 97* = 'C' Columns will be permuted. 98* 99* M (global input) INTEGER 100* The number of rows to be operated on, i.e. the number of rows 101* of the distributed submatrix sub( A ). M >= 0. 102* 103* N (global input) INTEGER 104* The number of columns to be operated on, i.e. the number of 105* columns of the distributed submatrix sub( A ). N >= 0. 106* 107* A (local input/local output) COMPLEX pointer into the 108* local memory to an array of dimension (LLD_A, LOCc(JA+N-1)). 109* On entry, this local array contains the local pieces of the 110* distributed matrix sub( A ) to which the row or columns 111* interchanges will be applied. On exit, this array contains 112* the local pieces of the permuted distributed matrix. 113* 114* IA (global input) INTEGER 115* The row index in the global array A indicating the first 116* row of sub( A ). 117* 118* JA (global input) INTEGER 119* The column index in the global array A indicating the 120* first column of sub( A ). 121* 122* DESCA (global and local input) INTEGER array of dimension DLEN_. 123* The array descriptor for the distributed matrix A. 124* 125* IPIV (input) INTEGER array, dimension >= LOCr(M_A)+MB_A if 126* ROWCOL = 'R', LOCc(N_A)+NB_A otherwise. It contains 127* the pivoting information. IPIV(i) is the global row (column), 128* local row (column) i was swapped with. The last piece of the 129* array of size MB_A (resp. NB_A) is used as workspace. IPIV is 130* tied to the distributed matrix A. 131* 132* IP (global input) INTEGER 133* IPIV's global row index, which points to the beginning of the 134* submatrix which is to be operated on. 135* 136* JP (global input) INTEGER 137* IPIV's global column index, which points to the beginning of 138* the submatrix which is to be operated on. 139* 140* DESCIP (global and local input) INTEGER array of dimension 8 141* The array descriptor for the distributed matrix IPIV. 142* 143* ===================================================================== 144* 145* .. Parameters .. 146 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_, 147 $ LLD_, MB_, M_, NB_, N_, RSRC_ 148 PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1, 149 $ CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, 150 $ RSRC_ = 7, CSRC_ = 8, LLD_ = 9 ) 151* .. 152* .. Local Scalars .. 153 LOGICAL FORWRD, ROWPVT 154 INTEGER I, IB, ICTXT, ICURCOL, ICURROW, IIP, IP1, ITMP, 155 $ IPVWRK, J, JB, JJP, JP1, K, MA, MBA, MYCOL, 156 $ MYROW, NBA, NPCOL, NPROW 157* .. 158* .. External Subroutines .. 159 EXTERNAL BLACS_GRIDINFO, IGEBS2D, IGEBR2D, INFOG2L, 160 $ PCSWAP 161* .. 162* .. External Functions .. 163 LOGICAL LSAME 164 INTEGER ICEIL, NUMROC 165 EXTERNAL ICEIL, LSAME, NUMROC 166* .. 167* .. Intrinsic Functions .. 168 INTRINSIC MIN, MOD 169* .. 170* .. Executable Statements .. 171* 172 ROWPVT = LSAME( ROWCOL, 'R' ) 173 IF( ROWPVT ) THEN 174 IF( M.LE.1 .OR. N.LT.1 ) 175 $ RETURN 176 ELSE 177 IF( M.LT.1 .OR. N.LE.1 ) 178 $ RETURN 179 END IF 180 FORWRD = LSAME( DIREC, 'F' ) 181* 182* 183* Get grid and matrix parameters 184* 185 MA = DESCA( M_ ) 186 MBA = DESCA( MB_ ) 187 NBA = DESCA( NB_ ) 188 ICTXT = DESCA( CTXT_ ) 189 CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL ) 190* 191* If I'm applying pivots from beginning to end (e.g., repeating 192* pivoting done earlier). Thus this section computes P * sub( A ). 193* 194 IF( FORWRD ) THEN 195 CALL INFOG2L( IP, JP, DESCIP, NPROW, NPCOL, MYROW, MYCOL, 196 $ IIP, JJP, ICURROW, ICURCOL ) 197* 198* If I'm pivoting the rows of sub( A ) 199* 200 IF( ROWPVT ) THEN 201 IPVWRK = NUMROC( DESCIP( M_ ), DESCIP( MB_ ), MYROW, 202 $ DESCIP( RSRC_ ), NPROW ) + 1 - 203 $ DESCIP( MB_ ) 204* 205* Loop over rows of sub( A ) 206* 207 I = IA 208 IB = MIN( M, ICEIL( IA, MBA ) * MBA - IA + 1 ) 209 10 CONTINUE 210* 211* Find local pointer into IPIV, and broadcast this block's 212* pivot information to everyone in process column 213* 214 IF( MYROW.EQ.ICURROW ) THEN 215 CALL IGEBS2D( ICTXT, 'Columnwise', ' ', IB, 1, 216 $ IPIV( IIP ), IB ) 217 ITMP = IIP 218 IIP = IIP + IB 219 ELSE 220 ITMP = IPVWRK 221 CALL IGEBR2D( ICTXT, 'Columnwise', ' ', IB, 1, 222 $ IPIV( ITMP ), IB, ICURROW, MYCOL ) 223 END IF 224* 225* Pivot the block of rows 226* 227 DO 20 K = I, I+IB-1 228 IP1 = IPIV( ITMP ) - IP + IA 229 IF( IP1.NE.K ) 230 $ CALL PCSWAP( N, A, K, JA, DESCA, MA, A, IP1, JA, 231 $ DESCA, MA ) 232 ITMP = ITMP + 1 233 20 CONTINUE 234* 235* Go on to next row of processes, increment row counter, 236* and figure number of rows to pivot next 237* 238 ICURROW = MOD( ICURROW+1, NPROW ) 239 I = I + IB 240 IB = MIN( MBA, M-I+IA ) 241 IF( IB .GT. 0 ) GOTO 10 242* 243* If I am pivoting the columns of sub( A ) 244* 245 ELSE 246 IPVWRK = NUMROC( DESCIP( N_ ), DESCIP( NB_ ), MYCOL, 247 $ DESCIP( CSRC_ ), NPCOL ) + 1 - 248 $ DESCIP( NB_ ) 249* 250* Loop over columns of sub( A ) 251* 252 J = JA 253 JB = MIN( N, ICEIL( JA, NBA ) * NBA - JA + 1 ) 254 30 CONTINUE 255* 256* Find local pointer into IPIV, and broadcast this block's 257* pivot information to everyone in process row 258* 259 IF( MYCOL.EQ.ICURCOL ) THEN 260 CALL IGEBS2D( ICTXT, 'Rowwise', ' ', JB, 1, 261 $ IPIV( JJP ), JB ) 262 ITMP = JJP 263 JJP = JJP + JB 264 ELSE 265 ITMP = IPVWRK 266 CALL IGEBR2D( ICTXT, 'Rowwise', ' ', JB, 1, 267 $ IPIV( ITMP ), JB, MYROW, ICURCOL ) 268 END IF 269* 270* Pivot the block of columns 271* 272 DO 40 K = J, J+JB-1 273 JP1 = IPIV( ITMP ) - JP + JA 274 IF( JP1.NE.K ) 275 $ CALL PCSWAP( M, A, IA, K, DESCA, 1, A, IA, JP1, 276 $ DESCA, 1 ) 277 ITMP = ITMP + 1 278 40 CONTINUE 279* 280* Go on to next column of processes, increment column 281* counter, and figure number of columns to pivot next 282* 283 ICURCOL = MOD( ICURCOL+1, NPCOL ) 284 J = J + JB 285 JB = MIN( NBA, N-J+JA ) 286 IF( JB .GT. 0 ) GOTO 30 287 END IF 288* 289* If I want to apply pivots in reverse order, i.e. reversing 290* pivoting done earlier. Thus this section computes 291* inv( P ) * sub( A ). 292* 293 ELSE 294* 295* If I'm pivoting the rows of sub( A ) 296* 297 IF( ROWPVT ) THEN 298 CALL INFOG2L( IP+M-1, JP, DESCIP, NPROW, NPCOL, MYROW, 299 $ MYCOL, IIP, JJP, ICURROW, ICURCOL ) 300* 301 IPVWRK = NUMROC( DESCIP( M_ ), DESCIP( MB_ ), MYROW, 302 $ DESCIP( RSRC_ ), NPROW ) + 1 - 303 $ DESCIP( MB_ ) 304* 305* If I'm not in the current process row, my IIP points out 306* past end of pivot vector (since I don't own a piece of the 307* last row). Adjust IIP so it points at last pivot entry. 308* 309 IF( MYROW.NE.ICURROW ) IIP = IIP - 1 310* 311* Loop over rows in reverse order, starting at last row 312* 313 I = IA + M - 1 314 IB = MOD( I, MBA ) 315 IF( IB .EQ. 0 ) IB = MBA 316 IB = MIN( IB, M ) 317 50 CONTINUE 318* 319* Find local pointer into IPIV, and broadcast this block's 320* pivot information to everyone in process column 321* 322 IF( MYROW.EQ.ICURROW ) THEN 323 ITMP = IIP 324 IIP = IIP - IB 325 CALL IGEBS2D( ICTXT, 'Columnwise', ' ', IB, 1, 326 $ IPIV( IIP+1 ), IB ) 327 ELSE 328 CALL IGEBR2D( ICTXT, 'Columnwise', ' ', IB, 1, 329 $ IPIV( IPVWRK ), IB, ICURROW, MYCOL ) 330 ITMP = IPVWRK + IB - 1 331 END IF 332* 333* Pivot the block of rows 334* 335 DO 60 K = I, I-IB+1, -1 336 IP1 = IPIV( ITMP ) - IP + IA 337 IF( IP1.NE.K ) 338 $ CALL PCSWAP( N, A, K, JA, DESCA, MA, A, IP1, JA, 339 $ DESCA, MA ) 340 ITMP = ITMP - 1 341 60 CONTINUE 342* 343* Go to previous row of processes, decrement row counter, 344* and figure number of rows to be pivoted next 345* 346 ICURROW = MOD( NPROW+ICURROW-1, NPROW ) 347 I = I - IB 348 IB = MIN( MBA, I-IA+1 ) 349 IF( IB .GT. 0 ) GOTO 50 350* 351* Otherwise, I'm pivoting the columns of sub( A ) 352* 353 ELSE 354 CALL INFOG2L( IP, JP+N-1, DESCIP, NPROW, NPCOL, MYROW, 355 $ MYCOL, IIP, JJP, ICURROW, ICURCOL ) 356 IPVWRK = NUMROC( DESCIP( N_ ), DESCIP( NB_ ), MYCOL, 357 $ DESCIP( CSRC_ ), NPCOL ) + 1 - 358 $ DESCIP( NB_ ) 359* 360* If I'm not in the current process column, my JJP points out 361* past end of pivot vector (since I don't own a piece of the 362* last column). Adjust JJP so it points at last pivot entry. 363* 364 IF( MYCOL.NE.ICURCOL ) JJP = JJP - 1 365* 366* Loop over columns in reverse order starting at last column 367* 368 J = JA + N - 1 369 JB = MOD( J, NBA ) 370 IF( JB .EQ. 0 ) JB = NBA 371 JB = MIN( JB, N ) 372 70 CONTINUE 373* 374* Find local pointer into IPIV, and broadcast this block's 375* pivot information to everyone in process row 376* 377 IF( MYCOL.EQ.ICURCOL ) THEN 378 ITMP = JJP 379 JJP = JJP - JB 380 CALL IGEBS2D( ICTXT, 'Rowwise', ' ', JB, 1, 381 $ IPIV( JJP+1 ), JB ) 382 ELSE 383 CALL IGEBR2D( ICTXT, 'Rowwise', ' ', JB, 1, 384 $ IPIV( IPVWRK ), JB, MYROW, ICURCOL ) 385 ITMP = IPVWRK + JB - 1 386 END IF 387* 388* Pivot a block of columns 389* 390 DO 80 K = J, J-JB+1, -1 391 JP1 = IPIV( ITMP ) - JP + JA 392 IF( JP1.NE.K ) 393 $ CALL PCSWAP( M, A, IA, K, DESCA, 1, A, IA, JP1, 394 $ DESCA, 1 ) 395 ITMP = ITMP - 1 396 80 CONTINUE 397* 398* Go to previous row of processes, decrement row counter, 399* and figure number of rows to be pivoted next 400* 401 ICURCOL = MOD( NPCOL+ICURCOL-1, NPCOL ) 402 J = J - JB 403 JB = MIN( NBA, J-JA+1 ) 404 IF( JB .GT. 0 ) GOTO 70 405 END IF 406* 407 END IF 408* 409 RETURN 410* 411* End PCLAPV2 412* 413 END 414