|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
00001 SUBROUTINE PZUNM2L( SIDE, TRANS, M, N, K, A, IA, JA, DESCA, TAU, 00002 $ C, IC, JC, DESCC, WORK, LWORK, INFO ) 00003 * 00004 * -- ScaLAPACK routine (version 1.7) -- 00005 * University of Tennessee, Knoxville, Oak Ridge National Laboratory, 00006 * and University of California, Berkeley. 00007 * May 25, 2001 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER SIDE, TRANS 00011 INTEGER IA, IC, INFO, JA, JC, K, LWORK, M, N 00012 * .. 00013 * .. Array Arguments .. 00014 INTEGER DESCA( * ), DESCC( * ) 00015 COMPLEX*16 A( * ), C( * ), TAU( * ), WORK( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * PZUNM2L overwrites the general complex M-by-N distributed matrix 00022 * sub( C ) = C(IC:IC+M-1,JC:JC+N-1) with 00023 * 00024 * SIDE = 'L' SIDE = 'R' 00025 * TRANS = 'N': Q * sub( C ) sub( C ) * Q 00026 * TRANS = 'C': Q**H * sub( C ) sub( C ) * Q**H 00027 * 00028 * where Q is a complex unitary distributed matrix defined as the 00029 * product of K elementary reflectors 00030 * 00031 * Q = H(k) . . . H(2) H(1) 00032 * 00033 * as returned by PZGEQLF. Q is of order M if SIDE = 'L' and of order N 00034 * if SIDE = 'R'. 00035 * 00036 * Notes 00037 * ===== 00038 * 00039 * Each global data object is described by an associated description 00040 * vector. This vector stores the information required to establish 00041 * the mapping between an object element and its corresponding process 00042 * and memory location. 00043 * 00044 * Let A be a generic term for any 2D block cyclicly distributed array. 00045 * Such a global array has an associated description vector DESCA. 00046 * In the following comments, the character _ should be read as 00047 * "of the global array". 00048 * 00049 * NOTATION STORED IN EXPLANATION 00050 * --------------- -------------- -------------------------------------- 00051 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case, 00052 * DTYPE_A = 1. 00053 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating 00054 * the BLACS process grid A is distribu- 00055 * ted over. The context itself is glo- 00056 * bal, but the handle (the integer 00057 * value) may vary. 00058 * M_A (global) DESCA( M_ ) The number of rows in the global 00059 * array A. 00060 * N_A (global) DESCA( N_ ) The number of columns in the global 00061 * array A. 00062 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute 00063 * the rows of the array. 00064 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute 00065 * the columns of the array. 00066 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first 00067 * row of the array A is distributed. 00068 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the 00069 * first column of the array A is 00070 * distributed. 00071 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local 00072 * array. LLD_A >= MAX(1,LOCr(M_A)). 00073 * 00074 * Let K be the number of rows or columns of a distributed matrix, 00075 * and assume that its process grid has dimension p x q. 00076 * LOCr( K ) denotes the number of elements of K that a process 00077 * would receive if K were distributed over the p processes of its 00078 * process column. 00079 * Similarly, LOCc( K ) denotes the number of elements of K that a 00080 * process would receive if K were distributed over the q processes of 00081 * its process row. 00082 * The values of LOCr() and LOCc() may be determined via a call to the 00083 * ScaLAPACK tool function, NUMROC: 00084 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ), 00085 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ). 00086 * An upper bound for these quantities may be computed by: 00087 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A 00088 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A 00089 * 00090 * Arguments 00091 * ========= 00092 * 00093 * SIDE (global input) CHARACTER 00094 * = 'L': apply Q or Q**H from the Left; 00095 * = 'R': apply Q or Q**H from the Right. 00096 * 00097 * TRANS (global input) CHARACTER 00098 * = 'N': No transpose, apply Q; 00099 * = 'C': Conjugate transpose, apply Q**H. 00100 * 00101 * M (global input) INTEGER 00102 * The number of rows to be operated on i.e the number of rows 00103 * of the distributed submatrix sub( C ). M >= 0. 00104 * 00105 * N (global input) INTEGER 00106 * The number of columns to be operated on i.e the number of 00107 * columns of the distributed submatrix sub( C ). N >= 0. 00108 * 00109 * K (global input) INTEGER 00110 * The number of elementary reflectors whose product defines the 00111 * matrix Q. If SIDE = 'L', M >= K >= 0, if SIDE = 'R', 00112 * N >= K >= 0. 00113 * 00114 * A (local input) COMPLEX*16 pointer into the local memory 00115 * to an array of dimension (LLD_A,LOCc(JA+K-1)). On entry, the 00116 * j-th column must contain the vector which defines the elemen- 00117 * tary reflector H(j), JA <= j <= JA+K-1, as returned by 00118 * PZGEQLF in the K columns of its distributed matrix 00119 * argument A(IA:*,JA:JA+K-1). A(IA:*,JA:JA+K-1) is modified by 00120 * the routine but restored on exit. 00121 * If SIDE = 'L', LLD_A >= MAX( 1, LOCr(IA+M-1) ), 00122 * if SIDE = 'R', LLD_A >= MAX( 1, LOCr(IA+N-1) ). 00123 * 00124 * IA (global input) INTEGER 00125 * The row index in the global array A indicating the first 00126 * row of sub( A ). 00127 * 00128 * JA (global input) INTEGER 00129 * The column index in the global array A indicating the 00130 * first column of sub( A ). 00131 * 00132 * DESCA (global and local input) INTEGER array of dimension DLEN_. 00133 * The array descriptor for the distributed matrix A. 00134 * 00135 * TAU (local input) COMPLEX*16, array, dimension LOCc(JA+N-1) 00136 * This array contains the scalar factors TAU(j) of the 00137 * elementary reflectors H(j) as returned by PZGEQLF. 00138 * TAU is tied to the distributed matrix A. 00139 * 00140 * C (local input/local output) COMPLEX*16 pointer into the 00141 * local memory to an array of dimension (LLD_C,LOCc(JC+N-1)). 00142 * On entry, the local pieces of the distributed matrix sub(C). 00143 * On exit, sub( C ) is overwritten by Q*sub( C ) or Q'*sub( C ) 00144 * or sub( C )*Q' or sub( C )*Q. 00145 * 00146 * IC (global input) INTEGER 00147 * The row index in the global array C indicating the first 00148 * row of sub( C ). 00149 * 00150 * JC (global input) INTEGER 00151 * The column index in the global array C indicating the 00152 * first column of sub( C ). 00153 * 00154 * DESCC (global and local input) INTEGER array of dimension DLEN_. 00155 * The array descriptor for the distributed matrix C. 00156 * 00157 * WORK (local workspace/local output) COMPLEX*16 array, 00158 * dimension (LWORK) 00159 * On exit, WORK(1) returns the minimal and optimal LWORK. 00160 * 00161 * LWORK (local or global input) INTEGER 00162 * The dimension of the array WORK. 00163 * LWORK is local input and must be at least 00164 * If SIDE = 'L', LWORK >= MpC0 + MAX( 1, NqC0 ); 00165 * if SIDE = 'R', LWORK >= NqC0 + MAX( MAX( 1, MpC0 ), NUMROC( 00166 * NUMROC( N+ICOFFC,NB_A,0,0,NPCOL ),NB_A,0,0,LCMQ ) ); 00167 * 00168 * where LCMQ = LCM / NPCOL with LCM = ICLM( NPROW, NPCOL ), 00169 * 00170 * IROFFC = MOD( IC-1, MB_C ), ICOFFC = MOD( JC-1, NB_C ), 00171 * ICROW = INDXG2P( IC, MB_C, MYROW, RSRC_C, NPROW ), 00172 * ICCOL = INDXG2P( JC, NB_C, MYCOL, CSRC_C, NPCOL ), 00173 * MpC0 = NUMROC( M+IROFFC, MB_C, MYROW, ICROW, NPROW ), 00174 * NqC0 = NUMROC( N+ICOFFC, NB_C, MYCOL, ICCOL, NPCOL ), 00175 * 00176 * ILCM, INDXG2P and NUMROC are ScaLAPACK tool functions; 00177 * MYROW, MYCOL, NPROW and NPCOL can be determined by calling 00178 * the subroutine BLACS_GRIDINFO. 00179 * 00180 * If LWORK = -1, then LWORK is global input and a workspace 00181 * query is assumed; the routine only calculates the minimum 00182 * and optimal size for all work arrays. Each of these 00183 * values is returned in the first entry of the corresponding 00184 * work array, and no error message is issued by PXERBLA. 00185 * 00186 * 00187 * INFO (local output) INTEGER 00188 * = 0: successful exit 00189 * < 0: If the i-th argument is an array and the j-entry had 00190 * an illegal value, then INFO = -(i*100+j), if the i-th 00191 * argument is a scalar and had an illegal value, then 00192 * INFO = -i. 00193 * 00194 * Alignment requirements 00195 * ====================== 00196 * 00197 * The distributed submatrices A(IA:*, JA:*) and C(IC:IC+M-1,JC:JC+N-1) 00198 * must verify some alignment properties, namely the following 00199 * expressions should be true: 00200 * 00201 * If SIDE = 'L', 00202 * ( MB_A.EQ.MB_C .AND. IROFFA.EQ.IROFFC .AND. IAROW.EQ.ICROW ) 00203 * If SIDE = 'R', 00204 * ( MB_A.EQ.NB_C .AND. IROFFA.EQ.ICOFFC ) 00205 * 00206 * ===================================================================== 00207 * 00208 * .. Parameters .. 00209 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_, 00210 $ LLD_, MB_, M_, NB_, N_, RSRC_ 00211 PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1, 00212 $ CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, 00213 $ RSRC_ = 7, CSRC_ = 8, LLD_ = 9 ) 00214 COMPLEX*16 ONE 00215 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) ) 00216 * .. 00217 * .. Local Scalars .. 00218 LOGICAL LEFT, LQUERY, NOTRAN 00219 CHARACTER COLBTOP, ROWBTOP 00220 INTEGER IACOL, IAROW, ICCOL, ICOFFC, ICROW, ICTXT, ICC, 00221 $ II, IROFFA, IROFFC, J, J1, J2, J3, JCC, JJ, 00222 $ LCM, LCMQ, LWMIN, MI, MP, MPC0, MYCOL, MYROW, 00223 $ NI, NPCOL, NPROW, NQ, NQC0 00224 COMPLEX*16 AJJ 00225 * .. 00226 * .. External Subroutines .. 00227 EXTERNAL BLACS_ABORT, BLACS_GRIDINFO, CHK1MAT, INFOG2L, 00228 $ PB_TOPGET, PB_TOPSET, PXERBLA, PZELSET, 00229 $ PZELSET2, PZLARF, PZLARFC, ZGEBR2D, ZGEBS2D, 00230 $ ZGERV2D, ZGESD2D, ZSCAL 00231 * .. 00232 * .. External Functions .. 00233 LOGICAL LSAME 00234 INTEGER ILCM, INDXG2P, NUMROC 00235 EXTERNAL ILCM, INDXG2P, LSAME, NUMROC 00236 * .. 00237 * .. Intrinsic Functions .. 00238 INTRINSIC DBLE, DCMPLX, DCONJG, MAX, MOD 00239 * .. 00240 * .. Executable Statements .. 00241 * 00242 * Get grid parameters 00243 * 00244 ICTXT = DESCA( CTXT_ ) 00245 CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL ) 00246 * 00247 * Test the input parameters 00248 * 00249 INFO = 0 00250 IF( NPROW.EQ.-1 ) THEN 00251 INFO = -(900+CTXT_) 00252 ELSE 00253 LEFT = LSAME( SIDE, 'L' ) 00254 NOTRAN = LSAME( TRANS, 'N' ) 00255 * 00256 * NQ is the order of Q 00257 * 00258 IF( LEFT ) THEN 00259 NQ = M 00260 CALL CHK1MAT( M, 3, K, 5, IA, JA, DESCA, 9, INFO ) 00261 ELSE 00262 NQ = N 00263 CALL CHK1MAT( N, 4, K, 5, IA, JA, DESCA, 9, INFO ) 00264 END IF 00265 CALL CHK1MAT( M, 3, N, 4, IC, JC, DESCC, 14, INFO ) 00266 IF( INFO.EQ.0 ) THEN 00267 IROFFA = MOD( IA-1, DESCA( MB_ ) ) 00268 IROFFC = MOD( IC-1, DESCC( MB_ ) ) 00269 ICOFFC = MOD( JC-1, DESCC( NB_ ) ) 00270 IAROW = INDXG2P( IA, DESCA( MB_ ), MYROW, DESCA( RSRC_ ), 00271 $ NPROW ) 00272 ICROW = INDXG2P( IC, DESCC( MB_ ), MYROW, DESCC( RSRC_ ), 00273 $ NPROW ) 00274 ICCOL = INDXG2P( JC, DESCC( NB_ ), MYCOL, DESCC( CSRC_ ), 00275 $ NPCOL ) 00276 MPC0 = NUMROC( M+IROFFC, DESCC( MB_ ), MYROW, ICROW, NPROW ) 00277 NQC0 = NUMROC( N+ICOFFC, DESCC( NB_ ), MYCOL, ICCOL, NPCOL ) 00278 * 00279 IF( LEFT ) THEN 00280 LWMIN = MPC0 + MAX( 1, NQC0 ) 00281 ELSE 00282 LCM = ILCM( NPROW, NPCOL ) 00283 LCMQ = LCM / NPCOL 00284 LWMIN = NQC0 + MAX( MAX( 1, MPC0 ), NUMROC( NUMROC( 00285 $ N+ICOFFC, DESCA( NB_ ), 0, 0, NPCOL ), 00286 $ DESCA( NB_ ), 0, 0, LCMQ ) ) 00287 END IF 00288 * 00289 WORK( 1 ) = DCMPLX( DBLE( LWMIN ) ) 00290 LQUERY = ( LWORK.EQ.-1 ) 00291 IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN 00292 INFO = -1 00293 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN 00294 INFO = -2 00295 ELSE IF( K.LT.0 .OR. K.GT.NQ ) THEN 00296 INFO = -5 00297 ELSE IF( .NOT.LEFT .AND. DESCA( MB_ ).NE.DESCC( NB_ ) ) THEN 00298 INFO = -(900+NB_) 00299 ELSE IF( LEFT .AND. IROFFA.NE.IROFFC ) THEN 00300 INFO = -12 00301 ELSE IF( LEFT .AND. IAROW.NE.ICROW ) THEN 00302 INFO = -12 00303 ELSE IF( .NOT.LEFT .AND. IROFFA.NE.ICOFFC ) THEN 00304 INFO = -13 00305 ELSE IF( LEFT .AND. DESCA( MB_ ).NE.DESCC( MB_ ) ) THEN 00306 INFO = -(1400+MB_) 00307 ELSE IF( ICTXT.NE.DESCC( CTXT_ ) ) THEN 00308 INFO = -(1400+CTXT_) 00309 ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00310 INFO = -16 00311 END IF 00312 END IF 00313 END IF 00314 IF( INFO.NE.0 ) THEN 00315 CALL PXERBLA( ICTXT, 'PZUNM2L', -INFO ) 00316 CALL BLACS_ABORT( ICTXT, 1 ) 00317 RETURN 00318 ELSE IF( LQUERY ) THEN 00319 RETURN 00320 END IF 00321 * 00322 * Quick return if possible 00323 * 00324 IF( M.EQ.0 .OR. N.EQ.0 .OR. K.EQ.0 ) 00325 $ RETURN 00326 * 00327 IF( DESCA( M_ ).EQ.1 ) THEN 00328 CALL INFOG2L( IA, JA, DESCA, NPROW, NPCOL, MYROW, MYCOL, II, 00329 $ JJ, IAROW, IACOL ) 00330 CALL INFOG2L( IC, JC, DESCC, NPROW, NPCOL, MYROW, MYCOL, ICC, 00331 $ JCC, ICROW, ICCOL ) 00332 IF( LEFT ) THEN 00333 IF( MYROW.EQ.IAROW ) THEN 00334 NQ = NUMROC( JC+N-1, DESCC( NB_ ), MYCOL, DESCC( CSRC_ ), 00335 $ NPCOL ) 00336 IF( MYCOL.EQ.IACOL ) THEN 00337 IF( NOTRAN ) THEN 00338 AJJ = ONE - TAU( JJ ) 00339 ELSE 00340 AJJ = ONE - DCONJG( TAU( JJ ) ) 00341 END IF 00342 CALL ZGEBS2D( ICTXT, 'Rowwise', ' ', 1, 1, AJJ, 1 ) 00343 CALL ZSCAL( NQ-JCC+1, AJJ, 00344 $ C( ICC+(JCC-1)*DESCC( LLD_ ) ), 00345 $ DESCC( LLD_ ) ) 00346 ELSE 00347 CALL ZGEBR2D( ICTXT, 'Rowwise', ' ', 1, 1, AJJ, 1, 00348 $ IAROW, IACOL ) 00349 CALL ZSCAL( NQ-JCC+1, AJJ, 00350 $ C( ICC+(JCC-1)*DESCC( LLD_ ) ), 00351 $ DESCC( LLD_ ) ) 00352 END IF 00353 END IF 00354 ELSE 00355 IF( MYCOL.EQ.IACOL ) THEN 00356 IF( NOTRAN ) THEN 00357 AJJ = ONE - TAU( JJ ) 00358 ELSE 00359 AJJ = ONE - DCONJG( TAU( JJ ) ) 00360 END IF 00361 END IF 00362 * 00363 IF( IACOL.NE.ICCOL ) THEN 00364 IF( MYCOL.EQ.IACOL ) 00365 $ CALL ZGESD2D( ICTXT, 1, 1, AJJ, 1, MYROW, ICCOL ) 00366 IF( MYCOL.EQ.ICCOL ) 00367 $ CALL ZGERV2D( ICTXT, 1, 1, AJJ, 1, MYROW, IACOL ) 00368 END IF 00369 * 00370 IF( MYCOL.EQ.ICCOL ) THEN 00371 MP = NUMROC( IC+M-1, DESCC( MB_ ), MYROW, DESCC( RSRC_ ), 00372 $ NPROW ) 00373 CALL ZSCAL( MP-ICC+1, AJJ, C( ICC+(JCC-1)* 00374 $ DESCC( LLD_ ) ), 1 ) 00375 END IF 00376 * 00377 END IF 00378 * 00379 ELSE 00380 * 00381 CALL PB_TOPGET( ICTXT, 'Broadcast', 'Rowwise', ROWBTOP ) 00382 CALL PB_TOPGET( ICTXT, 'Broadcast', 'Columnwise', COLBTOP ) 00383 * 00384 IF( LEFT .AND. NOTRAN .OR. .NOT.LEFT .AND. .NOT.NOTRAN ) THEN 00385 J1 = JA 00386 J2 = JA+K-1 00387 J3 = 1 00388 ELSE 00389 J1 = JA+K-1 00390 J2 = JA 00391 J3 = -1 00392 END IF 00393 * 00394 IF( LEFT ) THEN 00395 NI = N 00396 IF( NOTRAN ) THEN 00397 CALL PB_TOPSET( ICTXT, 'Broadcast', 'Rowwise', 'I-ring' ) 00398 ELSE 00399 CALL PB_TOPSET( ICTXT, 'Broadcast', 'Rowwise', 'D-ring' ) 00400 END IF 00401 CALL PB_TOPSET( ICTXT, 'Broadcast', 'Columnwise', ' ' ) 00402 ELSE 00403 MI = M 00404 END IF 00405 * 00406 DO 10 J = J1, J2, J3 00407 * 00408 IF( LEFT ) THEN 00409 * 00410 * H(j) or H(j)' is applied to C(ic:ic+m-k+j-ja,jc:jc+n-1) 00411 * 00412 MI = M - K + J - JA + 1 00413 ELSE 00414 * 00415 * H(j) or H(j)' is applied to C(ic:ic+m-1,jc:jc+n-k+j-ja) 00416 * 00417 NI = N - K + J - JA + 1 00418 END IF 00419 * 00420 * Apply H(j) or H(j)' 00421 * 00422 CALL PZELSET2( AJJ, A, IA+NQ-K+J-JA, J, DESCA, ONE ) 00423 IF( NOTRAN ) THEN 00424 CALL PZLARF( SIDE, MI, NI, A, IA, J, DESCA, 1, TAU, C, 00425 $ IC, JC, DESCC, WORK ) 00426 ELSE 00427 CALL PZLARFC( SIDE, MI, NI, A, IA, J, DESCA, 1, TAU, C, 00428 $ IC, JC, DESCC, WORK ) 00429 END IF 00430 CALL PZELSET( A, IA+NQ-K+J-JA, J, DESCA, AJJ ) 00431 * 00432 10 CONTINUE 00433 * 00434 CALL PB_TOPSET( ICTXT, 'Broadcast', 'Rowwise', ROWBTOP ) 00435 CALL PB_TOPSET( ICTXT, 'Broadcast', 'Columnwise', COLBTOP ) 00436 * 00437 END IF 00438 * 00439 WORK( 1 ) = DCMPLX( DBLE( LWMIN ) ) 00440 * 00441 RETURN 00442 * 00443 * End of PZUNM2L 00444 * 00445 END
