1 SUBROUTINE PCLACON( N, V, IV, JV, DESCV, X, IX, JX, DESCX, EST, 2 $ KASE ) 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 INTEGER IV, IX, JV, JX, KASE, N 11 REAL EST 12* .. 13* .. Array Arguments .. 14 INTEGER DESCV( * ), DESCX( * ) 15 COMPLEX V( * ), X( * ) 16* .. 17* 18* Purpose 19* ======= 20* 21* PCLACON estimates the 1-norm of a square, complex distributed matrix 22* A. Reverse communication is used for evaluating matrix-vector 23* products. X and V are aligned with the distributed matrix A, this 24* information is implicitly contained within IV, IX, DESCV, and DESCX. 25* 26* Notes 27* ===== 28* 29* Each global data object is described by an associated description 30* vector. This vector stores the information required to establish 31* the mapping between an object element and its corresponding process 32* and memory location. 33* 34* Let A be a generic term for any 2D block cyclicly distributed array. 35* Such a global array has an associated description vector DESCA. 36* In the following comments, the character _ should be read as 37* "of the global array". 38* 39* NOTATION STORED IN EXPLANATION 40* --------------- -------------- -------------------------------------- 41* DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case, 42* DTYPE_A = 1. 43* CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating 44* the BLACS process grid A is distribu- 45* ted over. The context itself is glo- 46* bal, but the handle (the integer 47* value) may vary. 48* M_A (global) DESCA( M_ ) The number of rows in the global 49* array A. 50* N_A (global) DESCA( N_ ) The number of columns in the global 51* array A. 52* MB_A (global) DESCA( MB_ ) The blocking factor used to distribute 53* the rows of the array. 54* NB_A (global) DESCA( NB_ ) The blocking factor used to distribute 55* the columns of the array. 56* RSRC_A (global) DESCA( RSRC_ ) The process row over which the first 57* row of the array A is distributed. 58* CSRC_A (global) DESCA( CSRC_ ) The process column over which the 59* first column of the array A is 60* distributed. 61* LLD_A (local) DESCA( LLD_ ) The leading dimension of the local 62* array. LLD_A >= MAX(1,LOCr(M_A)). 63* 64* Let K be the number of rows or columns of a distributed matrix, 65* and assume that its process grid has dimension p x q. 66* LOCr( K ) denotes the number of elements of K that a process 67* would receive if K were distributed over the p processes of its 68* process column. 69* Similarly, LOCc( K ) denotes the number of elements of K that a 70* process would receive if K were distributed over the q processes of 71* its process row. 72* The values of LOCr() and LOCc() may be determined via a call to the 73* ScaLAPACK tool function, NUMROC: 74* LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ), 75* LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ). 76* An upper bound for these quantities may be computed by: 77* LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A 78* LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A 79* 80* Arguments 81* ========= 82* 83* N (global input) INTEGER 84* The length of the distributed vectors V and X. N >= 0. 85* 86* V (local workspace) COMPLEX pointer into the local 87* memory to an array of dimension LOCr(N+MOD(IV-1,MB_V)). On 88* the final return, V = A*W, where EST = norm(V)/norm(W) 89* (W is not returned). 90* 91* IV (global input) INTEGER 92* The row index in the global array V indicating the first 93* row of sub( V ). 94* 95* JV (global input) INTEGER 96* The column index in the global array V indicating the 97* first column of sub( V ). 98* 99* DESCV (global and local input) INTEGER array of dimension DLEN_. 100* The array descriptor for the distributed matrix V. 101* 102* X (local input/local output) COMPLEX pointer into the 103* local memory to an array of dimension 104* LOCr(N+MOD(IX-1,MB_X)). On an intermediate return, X 105* should be overwritten by 106* A * X, if KASE=1, 107* A' * X, if KASE=2, 108* where A' is the conjugate transpose of A, and PCLACON must 109* be re-called with all the other parameters unchanged. 110* 111* IX (global input) INTEGER 112* The row index in the global array X indicating the first 113* row of sub( X ). 114* 115* JX (global input) INTEGER 116* The column index in the global array X indicating the 117* first column of sub( X ). 118* 119* DESCX (global and local input) INTEGER array of dimension DLEN_. 120* The array descriptor for the distributed matrix X. 121* 122* 123* EST (global output) REAL 124* An estimate (a lower bound) for norm(A). 125* 126* KASE (local input/local output) INTEGER 127* On the initial call to PCLACON, KASE should be 0. 128* On an intermediate return, KASE will be 1 or 2, indicating 129* whether X should be overwritten by A * X or A' * X. 130* On the final return from PCLACON, KASE will again be 0. 131* 132* Further Details 133* =============== 134* 135* The serial version CLACON has been contributed by Nick Higham, 136* University of Manchester. It was originally named SONEST, dated 137* March 16, 1988. 138* 139* Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of 140* a real or complex matrix, with applications to condition estimation", 141* ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988. 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 INTEGER ITMAX 152 PARAMETER ( ITMAX = 5 ) 153 REAL ONE, TWO 154 PARAMETER ( ONE = 1.0E+0, TWO = 2.0E+0 ) 155 COMPLEX CZERO, CONE 156 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 157 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 158* .. 159* .. Local Scalars .. 160 INTEGER I, ICTXT, IIVX, IMAXROW, IOFFVX, IROFF, ITER, 161 $ IVXCOL, IVXROW, J, JLAST, JJVX, JUMP, K, 162 $ MYCOL, MYROW, NP, NPCOL, NPROW 163 REAL ALTSGN, ESTOLD, SAFMIN, TEMP 164 COMPLEX JLMAX, XMAX 165* .. 166* .. Local Arrays .. 167 COMPLEX WORK( 2 ) 168* .. 169* .. External Subroutines .. 170 EXTERNAL BLACS_GRIDINFO, CCOPY, CGEBR2D, CGEBS2D, 171 $ INFOG2L, PCELGET, PCMAX1, 172 $ PSCSUM1, SGEBR2D, SGEBS2D 173* .. 174* .. External Functions .. 175 INTEGER INDXG2L, INDXG2P, INDXL2G, NUMROC 176 REAL PSLAMCH 177 EXTERNAL INDXG2L, INDXG2P, INDXL2G, NUMROC, PSLAMCH 178* .. 179* .. Intrinsic Functions .. 180 INTRINSIC ABS, CMPLX, REAL 181* .. 182* .. Save statement .. 183 SAVE 184* .. 185* .. Executable Statements .. 186* 187* Get grid parameters. 188* 189 ICTXT = DESCX( CTXT_ ) 190 CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL ) 191* 192 CALL INFOG2L( IX, JX, DESCX, NPROW, NPCOL, MYROW, MYCOL, 193 $ IIVX, JJVX, IVXROW, IVXCOL ) 194 IF( MYCOL.NE.IVXCOL ) 195 $ RETURN 196 IROFF = MOD( IX-1, DESCX( MB_ ) ) 197 NP = NUMROC( N+IROFF, DESCX( MB_ ), MYROW, IVXROW, NPROW ) 198 IF( MYROW.EQ.IVXROW ) 199 $ NP = NP - IROFF 200 IOFFVX = IIVX + (JJVX-1)*DESCX( LLD_ ) 201* 202 SAFMIN = PSLAMCH( ICTXT, 'Safe minimum' ) 203 IF( KASE.EQ.0 ) THEN 204 DO 10 I = IOFFVX, IOFFVX+NP-1 205 X( I ) = CMPLX( ONE / REAL( N ) ) 206 10 CONTINUE 207 KASE = 1 208 JUMP = 1 209 RETURN 210 END IF 211* 212 GO TO ( 20, 40, 70, 90, 120 )JUMP 213* 214* ................ ENTRY (JUMP = 1) 215* FIRST ITERATION. X HAS BEEN OVERWRITTEN BY A*X 216* 217 20 CONTINUE 218 IF( N.EQ.1 ) THEN 219 IF( MYROW.EQ.IVXROW ) THEN 220 V( IOFFVX ) = X( IOFFVX ) 221 EST = ABS( V( IOFFVX ) ) 222 CALL SGEBS2D( ICTXT, 'Columnwise', ' ', 1, 1, EST, 1 ) 223 ELSE 224 CALL SGEBR2D( ICTXT, 'Columnwise', ' ', 1, 1, EST, 1, 225 $ IVXROW, MYCOL ) 226 END IF 227* ... QUIT 228 GO TO 130 229 END IF 230 CALL PSCSUM1( N, EST, X, IX, JX, DESCX, 1 ) 231 IF( DESCX( M_ ).EQ.1 .AND. N.EQ.1 ) THEN 232 IF( MYROW.EQ.IVXROW ) THEN 233 CALL SGEBS2D( ICTXT, 'Columnwise', ' ', 1, 1, EST, 1 ) 234 ELSE 235 CALL SGEBR2D( ICTXT, 'Columnwise', ' ', 1, 1, EST, 1, 236 $ IVXROW, MYCOL ) 237 END IF 238 END IF 239* 240 DO 30 I = IOFFVX, IOFFVX+NP-1 241 IF( ABS( X( I ) ).GT.SAFMIN ) THEN 242 X( I ) = X( I ) / CMPLX( ABS( X( I ) ) ) 243 ELSE 244 X( I ) = CONE 245 END IF 246 30 CONTINUE 247 KASE = 2 248 JUMP = 2 249 RETURN 250* 251* ................ ENTRY (JUMP = 2) 252* FIRST ITERATION. X HAS BEEN OVERWRITTEN BY CTRANS(A)*X 253* 254 40 CONTINUE 255 CALL PCMAX1( N, XMAX, J, X, IX, JX, DESCX, 1 ) 256 IF( DESCX( M_ ).EQ.1 .AND. N.EQ.1 ) THEN 257 IF( MYROW.EQ.IVXROW ) THEN 258 WORK( 1 ) = XMAX 259 WORK( 2 ) = CMPLX( REAL( J ) ) 260 CALL CGEBS2D( ICTXT, 'Columnwise', ' ', 2, 1, WORK, 2 ) 261 ELSE 262 CALL CGEBR2D( ICTXT, 'Columnwise', ' ', 2, 1, WORK, 2, 263 $ IVXROW, MYCOL ) 264 XMAX = WORK( 1 ) 265 J = NINT( REAL( WORK( 2 ) ) ) 266 END IF 267 END IF 268 ITER = 2 269* 270* MAIN LOOP - ITERATIONS 2, 3,...,ITMAX 271* 272 50 CONTINUE 273 DO 60 I = IOFFVX, IOFFVX+NP-1 274 X( I ) = CZERO 275 60 CONTINUE 276 IMAXROW = INDXG2P( J, DESCX( MB_ ), MYROW, DESCX( RSRC_ ), NPROW ) 277 IF( MYROW.EQ.IMAXROW ) THEN 278 I = INDXG2L( J, DESCX( MB_ ), MYROW, DESCX( RSRC_ ), NPROW ) 279 X( I ) = CONE 280 END IF 281 KASE = 1 282 JUMP = 3 283 RETURN 284* 285* ................ ENTRY (JUMP = 3) 286* X HAS BEEN OVERWRITTEN BY A*X 287* 288 70 CONTINUE 289 CALL CCOPY( NP, X( IOFFVX ), 1, V( IOFFVX ), 1 ) 290 ESTOLD = EST 291 CALL PSCSUM1( N, EST, V, IV, JV, DESCV, 1 ) 292 IF( DESCV( M_ ).EQ.1 .AND. N.EQ.1 ) THEN 293 IF( MYROW.EQ.IVXROW ) THEN 294 CALL SGEBS2D( ICTXT, 'Columnwise', ' ', 1, 1, EST, 1 ) 295 ELSE 296 CALL SGEBR2D( ICTXT, 'Columnwise', ' ', 1, 1, EST, 1, 297 $ IVXROW, MYCOL ) 298 END IF 299 END IF 300* 301* TEST FOR CYCLING 302 IF( EST.LE.ESTOLD ) 303 $ GO TO 100 304* 305 DO 80 I = IOFFVX, IOFFVX+NP-1 306 IF( ABS( X( I ) ).GT.SAFMIN ) THEN 307 X( I ) = X( I ) / CMPLX( ABS( X( I ) ) ) 308 ELSE 309 X( I ) = CONE 310 END IF 311 80 CONTINUE 312 KASE = 2 313 JUMP = 4 314 RETURN 315* 316* ................ ENTRY (JUMP = 4) 317* X HAS BEEN OVERWRITTEN BY CTRANS(A)*X 318* 319 90 CONTINUE 320 JLAST = J 321 CALL PCMAX1( N, XMAX, J, X, IX, JX, DESCX, 1 ) 322 IF( DESCX( M_ ).EQ.1 .AND. N.EQ.1 ) THEN 323 IF( MYROW.EQ.IVXROW ) THEN 324 WORK( 1 ) = XMAX 325 WORK( 2 ) = CMPLX( REAL( J ) ) 326 CALL CGEBS2D( ICTXT, 'Columnwise', ' ', 2, 1, WORK, 2 ) 327 ELSE 328 CALL CGEBR2D( ICTXT, 'Columnwise', ' ', 2, 1, WORK, 2, 329 $ IVXROW, MYCOL ) 330 XMAX = WORK( 1 ) 331 J = NINT( REAL( WORK( 2 ) ) ) 332 END IF 333 END IF 334 CALL PCELGET( 'Columnwise', ' ', JLMAX, X, JLAST, JX, DESCX ) 335 IF( ( REAL( JLMAX ).NE.ABS( REAL( XMAX ) ) ).AND. 336 $ ( ITER.LT.ITMAX ) ) THEN 337 ITER = ITER + 1 338 GO TO 50 339 END IF 340* 341* ITERATION COMPLETE. FINAL STAGE. 342* 343 100 CONTINUE 344 DO 110 I = IOFFVX, IOFFVX+NP-1 345 K = INDXL2G( I-IOFFVX+IIVX, DESCX( MB_ ), MYROW, 346 $ DESCX( RSRC_ ), NPROW )-IX+1 347 IF( MOD( K, 2 ).EQ.0 ) THEN 348 ALTSGN = -ONE 349 ELSE 350 ALTSGN = ONE 351 END IF 352 X( I ) = CMPLX( ALTSGN*( ONE+REAL( K-1 ) / REAL( N-1 ) ) ) 353 110 CONTINUE 354 KASE = 1 355 JUMP = 5 356 RETURN 357* 358* ................ ENTRY (JUMP = 5) 359* X HAS BEEN OVERWRITTEN BY A*X 360* 361 120 CONTINUE 362 CALL PSCSUM1( N, TEMP, X, IX, JX, DESCX, 1 ) 363 IF( DESCX( M_ ).EQ.1 .AND. N.EQ.1 ) THEN 364 IF( MYROW.EQ.IVXROW ) THEN 365 CALL SGEBS2D( ICTXT, 'Columnwise', ' ', 1, 1, TEMP, 1 ) 366 ELSE 367 CALL SGEBR2D( ICTXT, 'Columnwise', ' ', 1, 1, TEMP, 1, 368 $ IVXROW, MYCOL ) 369 END IF 370 END IF 371 TEMP = TWO*( TEMP / REAL( 3*N ) ) 372 IF( TEMP.GT.EST ) THEN 373 CALL CCOPY( NP, X( IOFFVX ), 1, V( IOFFVX ), 1 ) 374 EST = TEMP 375 END IF 376* 377 130 CONTINUE 378 KASE = 0 379* 380 RETURN 381* 382* End of PCLACON 383* 384 END 385