1 DOUBLE PRECISION FUNCTION PDLANGE( NORM, M, N, A, IA, JA, DESCA, 2 $ WORK ) 3 IMPLICIT NONE 4* 5* -- ScaLAPACK auxiliary routine (version 1.7) -- 6* University of Tennessee, Knoxville, Oak Ridge National Laboratory, 7* and University of California, Berkeley. 8* May 1, 1997 9* 10* .. Scalar Arguments .. 11 CHARACTER NORM 12 INTEGER IA, JA, M, N 13* .. 14* .. Array Arguments .. 15 INTEGER DESCA( * ) 16 DOUBLE PRECISION A( * ), WORK( * ) 17* .. 18* 19* Purpose 20* ======= 21* 22* PDLANGE returns the value of the one norm, or the Frobenius norm, 23* or the infinity norm, or the element of largest absolute value of a 24* distributed matrix sub( A ) = A(IA:IA+M-1, JA:JA+N-1). 25* 26* PDLANGE returns the value 27* 28* ( max(abs(A(i,j))), NORM = 'M' or 'm' with IA <= i <= IA+M-1, 29* ( and JA <= j <= JA+N-1, 30* ( 31* ( norm1( sub( A ) ), NORM = '1', 'O' or 'o' 32* ( 33* ( normI( sub( A ) ), NORM = 'I' or 'i' 34* ( 35* ( normF( sub( A ) ), NORM = 'F', 'f', 'E' or 'e' 36* 37* where norm1 denotes the one norm of a matrix (maximum column sum), 38* normI denotes the infinity norm of a matrix (maximum row sum) and 39* normF denotes the Frobenius norm of a matrix (square root of sum of 40* squares). Note that max(abs(A(i,j))) is not a matrix norm. 41* 42* Notes 43* ===== 44* 45* Each global data object is described by an associated description 46* vector. This vector stores the information required to establish 47* the mapping between an object element and its corresponding process 48* and memory location. 49* 50* Let A be a generic term for any 2D block cyclicly distributed array. 51* Such a global array has an associated description vector DESCA. 52* In the following comments, the character _ should be read as 53* "of the global array". 54* 55* NOTATION STORED IN EXPLANATION 56* --------------- -------------- -------------------------------------- 57* DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case, 58* DTYPE_A = 1. 59* CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating 60* the BLACS process grid A is distribu- 61* ted over. The context itself is glo- 62* bal, but the handle (the integer 63* value) may vary. 64* M_A (global) DESCA( M_ ) The number of rows in the global 65* array A. 66* N_A (global) DESCA( N_ ) The number of columns in the global 67* array A. 68* MB_A (global) DESCA( MB_ ) The blocking factor used to distribute 69* the rows of the array. 70* NB_A (global) DESCA( NB_ ) The blocking factor used to distribute 71* the columns of the array. 72* RSRC_A (global) DESCA( RSRC_ ) The process row over which the first 73* row of the array A is distributed. 74* CSRC_A (global) DESCA( CSRC_ ) The process column over which the 75* first column of the array A is 76* distributed. 77* LLD_A (local) DESCA( LLD_ ) The leading dimension of the local 78* array. LLD_A >= MAX(1,LOCr(M_A)). 79* 80* Let K be the number of rows or columns of a distributed matrix, 81* and assume that its process grid has dimension p x q. 82* LOCr( K ) denotes the number of elements of K that a process 83* would receive if K were distributed over the p processes of its 84* process column. 85* Similarly, LOCc( K ) denotes the number of elements of K that a 86* process would receive if K were distributed over the q processes of 87* its process row. 88* The values of LOCr() and LOCc() may be determined via a call to the 89* ScaLAPACK tool function, NUMROC: 90* LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ), 91* LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ). 92* An upper bound for these quantities may be computed by: 93* LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A 94* LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A 95* 96* Arguments 97* ========= 98* 99* NORM (global input) CHARACTER 100* Specifies the value to be returned in PDLANGE as described 101* above. 102* 103* M (global input) INTEGER 104* The number of rows to be operated on i.e the number of rows 105* of the distributed submatrix sub( A ). When M = 0, PDLANGE 106* is set to zero. M >= 0. 107* 108* N (global input) INTEGER 109* The number of columns to be operated on i.e the number of 110* columns of the distributed submatrix sub( A ). When N = 0, 111* PDLANGE is set to zero. N >= 0. 112* 113* A (local input) DOUBLE PRECISION pointer into the local memory 114* to an array of dimension (LLD_A, LOCc(JA+N-1)) containing the 115* local pieces of the distributed matrix sub( A ). 116* 117* IA (global input) INTEGER 118* The row index in the global array A indicating the first 119* row of sub( A ). 120* 121* JA (global input) INTEGER 122* The column index in the global array A indicating the 123* first column of sub( A ). 124* 125* DESCA (global and local input) INTEGER array of dimension DLEN_. 126* The array descriptor for the distributed matrix A. 127* 128* WORK (local workspace) DOUBLE PRECISION array dimension (LWORK) 129* LWORK >= 0 if NORM = 'M' or 'm' (not referenced), 130* Nq0 if NORM = '1', 'O' or 'o', 131* Mp0 if NORM = 'I' or 'i', 132* 0 if NORM = 'F', 'f', 'E' or 'e' (not referenced), 133* where 134* 135* IROFFA = MOD( IA-1, MB_A ), ICOFFA = MOD( JA-1, NB_A ), 136* IAROW = INDXG2P( IA, MB_A, MYROW, RSRC_A, NPROW ), 137* IACOL = INDXG2P( JA, NB_A, MYCOL, CSRC_A, NPCOL ), 138* Mp0 = NUMROC( M+IROFFA, MB_A, MYROW, IAROW, NPROW ), 139* Nq0 = NUMROC( N+ICOFFA, NB_A, MYCOL, IACOL, NPCOL ), 140* 141* INDXG2P and NUMROC are ScaLAPACK tool functions; MYROW, 142* MYCOL, NPROW and NPCOL can be determined by calling the 143* subroutine BLACS_GRIDINFO. 144* 145* ===================================================================== 146* 147* .. Parameters .. 148 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_, 149 $ LLD_, MB_, M_, NB_, N_, RSRC_ 150 PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1, 151 $ CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, 152 $ RSRC_ = 7, CSRC_ = 8, LLD_ = 9 ) 153 DOUBLE PRECISION ONE, ZERO 154 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 155* .. 156* .. Local Scalars .. 157 INTEGER I, IACOL, IAROW, ICTXT, II, ICOFF, IOFFA, 158 $ IROFF, J, JJ, LDA, MP, MYCOL, MYROW, NPCOL, 159 $ NPROW, NQ 160 DOUBLE PRECISION SUM, VALUE 161* .. 162* .. Local Arrays .. 163 DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 ) 164* .. 165* .. External Subroutines .. 166 EXTERNAL BLACS_GRIDINFO, DCOMBSSQ, DGEBR2D, 167 $ DGEBS2D, DGAMX2D, DGSUM2D, DLASSQ, 168 $ INFOG2L, PDTREECOMB 169* .. 170* .. External Functions .. 171 LOGICAL LSAME 172 INTEGER IDAMAX, NUMROC 173 EXTERNAL LSAME, IDAMAX, NUMROC 174* .. 175* .. Intrinsic Functions .. 176 INTRINSIC ABS, MAX, MIN, MOD, SQRT 177* .. 178* .. Executable Statements .. 179* 180* Get grid parameters. 181* 182 ICTXT = DESCA( CTXT_ ) 183 CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL ) 184* 185 CALL INFOG2L( IA, JA, DESCA, NPROW, NPCOL, MYROW, MYCOL, II, JJ, 186 $ IAROW, IACOL ) 187 IROFF = MOD( IA-1, DESCA( MB_ ) ) 188 ICOFF = MOD( JA-1, DESCA( NB_ ) ) 189 MP = NUMROC( M+IROFF, DESCA( MB_ ), MYROW, IAROW, NPROW ) 190 NQ = NUMROC( N+ICOFF, DESCA( NB_ ), MYCOL, IACOL, NPCOL ) 191 IF( MYROW.EQ.IAROW ) 192 $ MP = MP - IROFF 193 IF( MYCOL.EQ.IACOL ) 194 $ NQ = NQ - ICOFF 195 LDA = DESCA( LLD_ ) 196* 197 IF( MIN( M, N ).EQ.0 ) THEN 198* 199 VALUE = ZERO 200* 201************************************************************************ 202* max norm 203* 204 ELSE IF( LSAME( NORM, 'M' ) ) THEN 205* 206* Find max(abs(A(i,j))). 207* 208 VALUE = ZERO 209 IF( NQ.GT.0 .AND. MP.GT.0 ) THEN 210 IOFFA = (JJ-1)*LDA 211 DO 20 J = JJ, JJ+NQ-1 212 DO 10 I = II, MP+II-1 213 VALUE = MAX( VALUE, ABS( A( IOFFA+I ) ) ) 214 10 CONTINUE 215 IOFFA = IOFFA + LDA 216 20 CONTINUE 217 END IF 218 CALL DGAMX2D( ICTXT, 'All', ' ', 1, 1, VALUE, 1, I, J, -1, 219 $ 0, 0 ) 220* 221************************************************************************ 222* one norm 223* 224 ELSE IF( LSAME( NORM, 'O' ) .OR. NORM.EQ.'1' ) THEN 225* 226* Find norm1( sub( A ) ). 227* 228 IF( NQ.GT.0 ) THEN 229 IOFFA = ( JJ - 1 ) * LDA 230 DO 40 J = JJ, JJ+NQ-1 231 SUM = ZERO 232 IF( MP.GT.0 ) THEN 233 DO 30 I = II, MP+II-1 234 SUM = SUM + ABS( A( IOFFA+I ) ) 235 30 CONTINUE 236 END IF 237 IOFFA = IOFFA + LDA 238 WORK( J-JJ+1 ) = SUM 239 40 CONTINUE 240 END IF 241* 242* Find sum of global matrix columns and store on row 0 of 243* process grid 244* 245 CALL DGSUM2D( ICTXT, 'Columnwise', ' ', 1, NQ, WORK, 1, 246 $ 0, MYCOL ) 247* 248* Find maximum sum of columns for 1-norm 249* 250 IF( MYROW.EQ.0 ) THEN 251 IF( NQ.GT.0 ) THEN 252 VALUE = WORK( IDAMAX( NQ, WORK, 1 ) ) 253 ELSE 254 VALUE = ZERO 255 END IF 256 CALL DGAMX2D( ICTXT, 'Rowwise', ' ', 1, 1, VALUE, 1, I, J, 257 $ -1, 0, 0 ) 258 END IF 259* 260************************************************************************ 261* inf norm 262* 263 ELSE IF( LSAME( NORM, 'I' ) ) THEN 264* 265* Find normI( sub( A ) ). 266* 267 IF( MP.GT.0 ) THEN 268 IOFFA = II + ( JJ - 1 ) * LDA 269 DO 60 I = II, II+MP-1 270 SUM = ZERO 271 IF( NQ.GT.0 ) THEN 272 DO 50 J = IOFFA, IOFFA + NQ*LDA - 1, LDA 273 SUM = SUM + ABS( A( J ) ) 274 50 CONTINUE 275 END IF 276 WORK( I-II+1 ) = SUM 277 IOFFA = IOFFA + 1 278 60 CONTINUE 279 END IF 280* 281* Find sum of global matrix rows and store on column 0 of 282* process grid 283* 284 CALL DGSUM2D( ICTXT, 'Rowwise', ' ', MP, 1, WORK, MAX( 1, MP ), 285 $ MYROW, 0 ) 286* 287* Find maximum sum of rows for supnorm 288* 289 IF( MYCOL.EQ.0 ) THEN 290 IF( MP.GT.0 ) THEN 291 VALUE = WORK( IDAMAX( MP, WORK, 1 ) ) 292 ELSE 293 VALUE = ZERO 294 END IF 295 CALL DGAMX2D( ICTXT, 'Columnwise', ' ', 1, 1, VALUE, 1, I, 296 $ J, -1, 0, 0 ) 297 END IF 298* 299************************************************************************ 300* Frobenius norm 301* SSQ(1) is scale 302* SSQ(2) is sum-of-squares 303* 304 ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN 305* 306* Find normF( sub( A ) ). 307* 308 SSQ(1) = ZERO 309 SSQ(2) = ONE 310 IOFFA = II + ( JJ - 1 ) * LDA 311 IF( NQ.GT.0 ) THEN 312 DO 70 J = IOFFA, IOFFA + NQ*LDA - 1, LDA 313 COLSSQ(1) = ZERO 314 COLSSQ(2) = ONE 315 CALL DLASSQ( MP, A( J ), 1, COLSSQ(1), COLSSQ(2) ) 316 CALL DCOMBSSQ( SSQ, COLSSQ ) 317 70 CONTINUE 318 END IF 319* 320* Perform the global scaled sum 321* 322 CALL PDTREECOMB( ICTXT, 'All', 2, SSQ, 0, 0, DCOMBSSQ ) 323 VALUE = SSQ( 1 ) * SQRT( SSQ( 2 ) ) 324* 325 END IF 326* 327 IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN 328 CALL DGEBS2D( ICTXT, 'All', ' ', 1, 1, VALUE, 1 ) 329 ELSE 330 CALL DGEBR2D( ICTXT, 'All', ' ', 1, 1, VALUE, 1, 0, 0 ) 331 END IF 332* 333 PDLANGE = VALUE 334* 335 RETURN 336* 337* End of PDLANGE 338* 339 END 340