1 SUBROUTINE PDGGQRF( N, M, P, A, IA, JA, DESCA, TAUA, B, IB, JB, 2 $ DESCB, TAUB, WORK, LWORK, INFO ) 3* 4* -- ScaLAPACK 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 IA, IB, INFO, JA, JB, LWORK, M, N, P 11* .. 12* .. Array Arguments .. 13 INTEGER DESCA( * ), DESCB( * ) 14 DOUBLE PRECISION A( * ), B( * ), TAUA( * ), TAUB( * ), WORK( * ) 15* .. 16* 17* Purpose 18* ======= 19* 20* PDGGQRF computes a generalized QR factorization of 21* an N-by-M matrix sub( A ) = A(IA:IA+N-1,JA:JA+M-1) and 22* an N-by-P matrix sub( B ) = B(IB:IB+N-1,JB:JB+P-1): 23* 24* sub( A ) = Q*R, sub( B ) = Q*T*Z, 25* 26* where Q is an N-by-N orthogonal matrix, Z is a P-by-P orthogonal 27* matrix, and R and T assume one of the forms: 28* 29* if N >= M, R = ( R11 ) M , or if N < M, R = ( R11 R12 ) N, 30* ( 0 ) N-M N M-N 31* M 32* 33* where R11 is upper triangular, and 34* 35* if N <= P, T = ( 0 T12 ) N, or if N > P, T = ( T11 ) N-P, 36* P-N N ( T21 ) P 37* P 38* 39* where T12 or T21 is upper triangular. 40* 41* In particular, if sub( B ) is square and nonsingular, the GQR 42* factorization of sub( A ) and sub( B ) implicitly gives the QR 43* factorization of inv( sub( B ) )* sub( A ): 44* 45* inv( sub( B ) )*sub( A )= Z'*(inv(T)*R) 46* 47* where inv( sub( B ) ) denotes the inverse of the matrix sub( B ), 48* and Z' denotes the transpose of matrix Z. 49* 50* Notes 51* ===== 52* 53* Each global data object is described by an associated description 54* vector. This vector stores the information required to establish 55* the mapping between an object element and its corresponding process 56* and memory location. 57* 58* Let A be a generic term for any 2D block cyclicly distributed array. 59* Such a global array has an associated description vector DESCA. 60* In the following comments, the character _ should be read as 61* "of the global array". 62* 63* NOTATION STORED IN EXPLANATION 64* --------------- -------------- -------------------------------------- 65* DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case, 66* DTYPE_A = 1. 67* CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating 68* the BLACS process grid A is distribu- 69* ted over. The context itself is glo- 70* bal, but the handle (the integer 71* value) may vary. 72* M_A (global) DESCA( M_ ) The number of rows in the global 73* array A. 74* N_A (global) DESCA( N_ ) The number of columns in the global 75* array A. 76* MB_A (global) DESCA( MB_ ) The blocking factor used to distribute 77* the rows of the array. 78* NB_A (global) DESCA( NB_ ) The blocking factor used to distribute 79* the columns of the array. 80* RSRC_A (global) DESCA( RSRC_ ) The process row over which the first 81* row of the array A is distributed. 82* CSRC_A (global) DESCA( CSRC_ ) The process column over which the 83* first column of the array A is 84* distributed. 85* LLD_A (local) DESCA( LLD_ ) The leading dimension of the local 86* array. LLD_A >= MAX(1,LOCr(M_A)). 87* 88* Let K be the number of rows or columns of a distributed matrix, 89* and assume that its process grid has dimension p x q. 90* LOCr( K ) denotes the number of elements of K that a process 91* would receive if K were distributed over the p processes of its 92* process column. 93* Similarly, LOCc( K ) denotes the number of elements of K that a 94* process would receive if K were distributed over the q processes of 95* its process row. 96* The values of LOCr() and LOCc() may be determined via a call to the 97* ScaLAPACK tool function, NUMROC: 98* LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ), 99* LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ). 100* An upper bound for these quantities may be computed by: 101* LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A 102* LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A 103* 104* Arguments 105* ========= 106* 107* N (global input) INTEGER 108* The number of rows to be operated on i.e the number of rows 109* of the distributed submatrices sub( A ) and sub( B ). N >= 0. 110* 111* M (global input) INTEGER 112* The number of columns to be operated on i.e the number of 113* columns of the distributed submatrix sub( A ). M >= 0. 114* 115* P (global input) INTEGER 116* The number of columns to be operated on i.e the number of 117* columns of the distributed submatrix sub( B ). P >= 0. 118* 119* A (local input/local output) DOUBLE PRECISION pointer into the 120* local memory to an array of dimension (LLD_A, LOCc(JA+M-1)). 121* On entry, the local pieces of the N-by-M distributed matrix 122* sub( A ) which is to be factored. On exit, the elements on 123* and above the diagonal of sub( A ) contain the min(N,M) by M 124* upper trapezoidal matrix R (R is upper triangular if N >= M); 125* the elements below the diagonal, with the array TAUA, 126* represent the orthogonal matrix Q as a product of min(N,M) 127* elementary reflectors (see Further Details). 128* 129* IA (global input) INTEGER 130* The row index in the global array A indicating the first 131* row of sub( A ). 132* 133* JA (global input) INTEGER 134* The column index in the global array A indicating the 135* first column of sub( A ). 136* 137* DESCA (global and local input) INTEGER array of dimension DLEN_. 138* The array descriptor for the distributed matrix A. 139* 140* TAUA (local output) DOUBLE PRECISION array, dimension 141* LOCc(JA+MIN(N,M)-1). This array contains the scalar factors 142* TAUA of the elementary reflectors which represent the 143* orthogonal matrix Q. TAUA is tied to the distributed matrix 144* A. (see Further Details). 145* 146* B (local input/local output) DOUBLE PRECISION pointer into the 147* local memory to an array of dimension (LLD_B, LOCc(JB+P-1)). 148* On entry, the local pieces of the N-by-P distributed matrix 149* sub( B ) which is to be factored. On exit, if N <= P, the 150* upper triangle of B(IB:IB+N-1,JB+P-N:JB+P-1) contains the 151* N by N upper triangular matrix T; if N > P, the elements on 152* and above the (N-P)-th subdiagonal contain the N by P upper 153* trapezoidal matrix T; the remaining elements, with the array 154* TAUB, represent the orthogonal matrix Z as a product of 155* elementary reflectors (see Further Details). 156* 157* IB (global input) INTEGER 158* The row index in the global array B indicating the first 159* row of sub( B ). 160* 161* JB (global input) INTEGER 162* The column index in the global array B indicating the 163* first column of sub( B ). 164* 165* DESCB (global and local input) INTEGER array of dimension DLEN_. 166* The array descriptor for the distributed matrix B. 167* 168* TAUB (local output) DOUBLE PRECISION array, dimension LOCr(IB+N-1) 169* This array contains the scalar factors of the elementary 170* reflectors which represent the orthogonal unitary matrix Z. 171* TAUB is tied to the distributed matrix B (see Further 172* Details). 173* 174* WORK (local workspace/local output) DOUBLE PRECISION array, 175* dimension (LWORK) 176* On exit, WORK(1) returns the minimal and optimal LWORK. 177* 178* LWORK (local or global input) INTEGER 179* The dimension of the array WORK. 180* LWORK is local input and must be at least 181* LWORK >= MAX( NB_A * ( NpA0 + MqA0 + NB_A ), 182* MAX( (NB_A*(NB_A-1))/2, (PqB0 + NpB0)*NB_A ) + 183* NB_A * NB_A, 184* MB_B * ( NpB0 + PqB0 + MB_B ) ), where 185* 186* IROFFA = MOD( IA-1, MB_A ), ICOFFA = MOD( JA-1, NB_A ), 187* IAROW = INDXG2P( IA, MB_A, MYROW, RSRC_A, NPROW ), 188* IACOL = INDXG2P( JA, NB_A, MYCOL, CSRC_A, NPCOL ), 189* NpA0 = NUMROC( N+IROFFA, MB_A, MYROW, IAROW, NPROW ), 190* MqA0 = NUMROC( M+ICOFFA, NB_A, MYCOL, IACOL, NPCOL ), 191* 192* IROFFB = MOD( IB-1, MB_B ), ICOFFB = MOD( JB-1, NB_B ), 193* IBROW = INDXG2P( IB, MB_B, MYROW, RSRC_B, NPROW ), 194* IBCOL = INDXG2P( JB, NB_B, MYCOL, CSRC_B, NPCOL ), 195* NpB0 = NUMROC( N+IROFFB, MB_B, MYROW, IBROW, NPROW ), 196* PqB0 = NUMROC( P+ICOFFB, NB_B, MYCOL, IBCOL, NPCOL ), 197* 198* and NUMROC, INDXG2P are ScaLAPACK tool functions; 199* MYROW, MYCOL, NPROW and NPCOL can be determined by calling 200* the subroutine BLACS_GRIDINFO. 201* 202* If LWORK = -1, then LWORK is global input and a workspace 203* query is assumed; the routine only calculates the minimum 204* and optimal size for all work arrays. Each of these 205* values is returned in the first entry of the corresponding 206* work array, and no error message is issued by PXERBLA. 207* 208* INFO (global output) INTEGER 209* = 0: successful exit 210* < 0: If the i-th argument is an array and the j-entry had 211* an illegal value, then INFO = -(i*100+j), if the i-th 212* argument is a scalar and had an illegal value, then 213* INFO = -i. 214* 215* Further Details 216* =============== 217* 218* The matrix Q is represented as a product of elementary reflectors 219* 220* Q = H(ja) H(ja+1) . . . H(ja+k-1), where k = min(n,m). 221* 222* Each H(i) has the form 223* 224* H(i) = I - taua * v * v' 225* 226* where taua is a real scalar, and v is a real vector with 227* v(1:i-1) = 0 and v(i) = 1; v(i+1:n) is stored on exit in 228* A(ia+i:ia+n-1,ja+i-1), and taua in TAUA(ja+i-1). 229* To form Q explicitly, use ScaLAPACK subroutine PDORGQR. 230* To use Q to update another matrix, use ScaLAPACK subroutine PDORMQR. 231* 232* The matrix Z is represented as a product of elementary reflectors 233* 234* Z = H(ib) H(ib+1) . . . H(ib+k-1), where k = min(n,p). 235* 236* Each H(i) has the form 237* 238* H(i) = I - taub * v * v' 239* 240* where taub is a real scalar, and v is a real vector with 241* v(p-k+i+1:p) = 0 and v(p-k+i) = 1; v(1:p-k+i-1) is stored on exit in 242* B(ib+n-k+i-1,jb:jb+p-k+i-2), and taub in TAUB(ib+n-k+i-1). 243* To form Z explicitly, use ScaLAPACK subroutine PDORGRQ. 244* To use Z to update another matrix, use ScaLAPACK subroutine PDORMRQ. 245* 246* Alignment requirements 247* ====================== 248* 249* The distributed submatrices sub( A ) and sub( B ) must verify some 250* alignment properties, namely the following expression should be true: 251* 252* ( MB_A.EQ.MB_B .AND. IROFFA.EQ.IROFFB .AND. IAROW.EQ.IBROW ) 253* 254* ===================================================================== 255* 256* .. Parameters .. 257 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_, 258 $ LLD_, MB_, M_, NB_, N_, RSRC_ 259 PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1, 260 $ CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, 261 $ RSRC_ = 7, CSRC_ = 8, LLD_ = 9 ) 262* .. 263* .. Local Scalars .. 264 LOGICAL LQUERY 265 INTEGER IACOL, IAROW, IBCOL, IBROW, ICOFFA, ICOFFB, 266 $ ICTXT, IROFFA, IROFFB, LWMIN, MQA0, MYCOL, 267 $ MYROW, NPA0, NPB0, NPCOL, NPROW, PQB0 268* .. 269* .. External Subroutines .. 270 EXTERNAL BLACS_GRIDINFO, CHK1MAT, PCHK2MAT, PDGEQRF, 271 $ PDGERQF, PDORMQR, PXERBLA 272* .. 273* .. Local Arrays .. 274 INTEGER IDUM1( 1 ), IDUM2( 1 ) 275* .. 276* .. External Functions .. 277 INTEGER INDXG2P, NUMROC 278 EXTERNAL INDXG2P, NUMROC 279* .. 280* .. Intrinsic Functions .. 281 INTRINSIC DBLE, INT, MAX, MIN, MOD 282* .. 283* .. Executable Statements .. 284* 285* Get grid parameters 286* 287 ICTXT = DESCA( CTXT_ ) 288 CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL ) 289* 290* Test the input parameters 291* 292 INFO = 0 293 IF( NPROW.EQ.-1 ) THEN 294 INFO = -707 295 ELSE 296 CALL CHK1MAT( N, 1, M, 2, IA, JA, DESCA, 7, INFO ) 297 CALL CHK1MAT( N, 1, P, 3, IB, JB, DESCB, 12, INFO ) 298 IF( INFO.EQ.0 ) THEN 299 IROFFA = MOD( IA-1, DESCA( MB_ ) ) 300 ICOFFA = MOD( JA-1, DESCA( NB_ ) ) 301 IROFFB = MOD( IB-1, DESCB( MB_ ) ) 302 ICOFFB = MOD( JB-1, DESCB( NB_ ) ) 303 IAROW = INDXG2P( IA, DESCA( MB_ ), MYROW, DESCA( RSRC_ ), 304 $ NPROW ) 305 IACOL = INDXG2P( JA, DESCA( NB_ ), MYCOL, DESCA( CSRC_ ), 306 $ NPCOL ) 307 IBROW = INDXG2P( IB, DESCB( MB_ ), MYROW, DESCB( RSRC_ ), 308 $ NPROW ) 309 IBCOL = INDXG2P( JB, DESCB( NB_ ), MYCOL, DESCB( CSRC_ ), 310 $ NPCOL ) 311 NPA0 = NUMROC( N+IROFFA, DESCA( MB_ ), MYROW, IAROW, NPROW ) 312 MQA0 = NUMROC( M+ICOFFA, DESCA( NB_ ), MYCOL, IACOL, NPCOL ) 313 NPB0 = NUMROC( N+IROFFB, DESCB( MB_ ), MYROW, IBROW, NPROW ) 314 PQB0 = NUMROC( P+ICOFFB, DESCB( NB_ ), MYCOL, IBCOL, NPCOL ) 315 LWMIN = MAX( DESCA( NB_ ) * ( NPA0 + MQA0 + DESCA( NB_ ) ), 316 $ MAX( MAX( ( DESCA( NB_ )*( DESCA( NB_ ) - 1 ) ) / 2, 317 $ ( PQB0 + NPB0 ) * DESCA( NB_ ) ) + 318 $ DESCA( NB_ ) * DESCA( NB_ ), 319 $ DESCB( MB_ ) * ( NPB0 + PQB0 + DESCB( MB_ ) ) ) ) 320* 321 WORK( 1 ) = DBLE( LWMIN ) 322 LQUERY = ( LWORK.EQ.-1 ) 323 IF( IAROW.NE.IBROW .OR. IROFFA.NE.IROFFB ) THEN 324 INFO = -10 325 ELSE IF( DESCA( MB_ ).NE.DESCB( MB_ ) ) THEN 326 INFO = -1203 327 ELSE IF( ICTXT.NE.DESCB( CTXT_ ) ) THEN 328 INFO = -1207 329 ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 330 INFO = -15 331 END IF 332 END IF 333 IF( LQUERY ) THEN 334 IDUM1( 1 ) = -1 335 ELSE 336 IDUM1( 1 ) = 1 337 END IF 338 IDUM2( 1 ) = 15 339 CALL PCHK2MAT( N, 1, M, 2, IA, JA, DESCA, 7, N, 1, P, 3, IB, 340 $ JB, DESCB, 12, 1, IDUM1, IDUM2, INFO ) 341 END IF 342* 343 IF( INFO.NE.0 ) THEN 344 CALL PXERBLA( ICTXT, 'PDGGQRF', -INFO ) 345 RETURN 346 ELSE IF( LQUERY ) THEN 347 RETURN 348 END IF 349* 350* QR factorization of N-by-M matrix sub( A ): sub( A ) = Q*R 351* 352 CALL PDGEQRF( N, M, A, IA, JA, DESCA, TAUA, WORK, LWORK, INFO ) 353 LWMIN = INT( WORK( 1 ) ) 354* 355* Update sub( B ) := Q'*sub( B ). 356* 357 CALL PDORMQR( 'Left', 'Transpose', N, P, MIN( N, M ), A, IA, JA, 358 $ DESCA, TAUA, B, IB, JB, DESCB, WORK, LWORK, INFO ) 359 LWMIN = MIN( LWMIN, INT( WORK( 1 ) ) ) 360* 361* RQ factorization of N-by-P matrix sub( B ): sub( B ) = T*Z. 362* 363 CALL PDGERQF( N, P, B, IB, JB, DESCB, TAUB, WORK, LWORK, INFO ) 364 WORK( 1 ) = DBLE( MAX( LWMIN, INT( WORK( 1 ) ) ) ) 365* 366 RETURN 367* 368* End of PDGGQRF 369* 370 END 371