|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
00001 SUBROUTINE PDGELQ2( M, N, A, IA, JA, DESCA, TAU, WORK, LWORK, 00002 $ 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 INTEGER IA, INFO, JA, LWORK, M, N 00011 * .. 00012 * .. Array Arguments .. 00013 INTEGER DESCA( * ) 00014 DOUBLE PRECISION A( * ), TAU( * ), WORK( * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * PDGELQ2 computes a LQ factorization of a real distributed M-by-N 00021 * matrix sub( A ) = A(IA:IA+M-1,JA:JA+N-1) = L * Q. 00022 * 00023 * Notes 00024 * ===== 00025 * 00026 * Each global data object is described by an associated description 00027 * vector. This vector stores the information required to establish 00028 * the mapping between an object element and its corresponding process 00029 * and memory location. 00030 * 00031 * Let A be a generic term for any 2D block cyclicly distributed array. 00032 * Such a global array has an associated description vector DESCA. 00033 * In the following comments, the character _ should be read as 00034 * "of the global array". 00035 * 00036 * NOTATION STORED IN EXPLANATION 00037 * --------------- -------------- -------------------------------------- 00038 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case, 00039 * DTYPE_A = 1. 00040 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating 00041 * the BLACS process grid A is distribu- 00042 * ted over. The context itself is glo- 00043 * bal, but the handle (the integer 00044 * value) may vary. 00045 * M_A (global) DESCA( M_ ) The number of rows in the global 00046 * array A. 00047 * N_A (global) DESCA( N_ ) The number of columns in the global 00048 * array A. 00049 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute 00050 * the rows of the array. 00051 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute 00052 * the columns of the array. 00053 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first 00054 * row of the array A is distributed. 00055 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the 00056 * first column of the array A is 00057 * distributed. 00058 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local 00059 * array. LLD_A >= MAX(1,LOCr(M_A)). 00060 * 00061 * Let K be the number of rows or columns of a distributed matrix, 00062 * and assume that its process grid has dimension p x q. 00063 * LOCr( K ) denotes the number of elements of K that a process 00064 * would receive if K were distributed over the p processes of its 00065 * process column. 00066 * Similarly, LOCc( K ) denotes the number of elements of K that a 00067 * process would receive if K were distributed over the q processes of 00068 * its process row. 00069 * The values of LOCr() and LOCc() may be determined via a call to the 00070 * ScaLAPACK tool function, NUMROC: 00071 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ), 00072 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ). 00073 * An upper bound for these quantities may be computed by: 00074 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A 00075 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A 00076 * 00077 * Arguments 00078 * ========= 00079 * 00080 * M (global input) INTEGER 00081 * The number of rows to be operated on, i.e. the number of rows 00082 * of the distributed submatrix sub( A ). M >= 0. 00083 * 00084 * N (global input) INTEGER 00085 * The number of columns to be operated on, i.e. the number of 00086 * columns of the distributed submatrix sub( A ). N >= 0. 00087 * 00088 * A (local input/local output) DOUBLE PRECISION pointer into the 00089 * local memory to an array of dimension (LLD_A, LOCc(JA+N-1)). 00090 * On entry, the local pieces of the M-by-N distributed matrix 00091 * sub( A ) which is to be factored. On exit, the elements on 00092 * and below the diagonal of sub( A ) contain the M by min(M,N) 00093 * lower trapezoidal matrix L (L is lower triangular if M <= N); 00094 * the elements above the diagonal, with the array TAU, repre- 00095 * sent the orthogonal matrix Q as a product of elementary 00096 * reflectors (see Further Details). 00097 * 00098 * IA (global input) INTEGER 00099 * The row index in the global array A indicating the first 00100 * row of sub( A ). 00101 * 00102 * JA (global input) INTEGER 00103 * The column index in the global array A indicating the 00104 * first column of sub( A ). 00105 * 00106 * DESCA (global and local input) INTEGER array of dimension DLEN_. 00107 * The array descriptor for the distributed matrix A. 00108 * 00109 * TAU (local output) DOUBLE PRECISION array, dimension 00110 * LOCr(IA+MIN(M,N)-1). This array contains the scalar factors 00111 * of the elementary reflectors. TAU is tied to the distributed 00112 * matrix A. 00113 * 00114 * WORK (local workspace/local output) DOUBLE PRECISION array, 00115 * dimension (LWORK) 00116 * On exit, WORK(1) returns the minimal and optimal LWORK. 00117 * 00118 * LWORK (local or global input) INTEGER 00119 * The dimension of the array WORK. 00120 * LWORK is local input and must be at least 00121 * LWORK >= Nq0 + MAX( 1, Mp0 ), where 00122 * 00123 * IROFF = MOD( IA-1, MB_A ), ICOFF = MOD( JA-1, NB_A ), 00124 * IAROW = INDXG2P( IA, MB_A, MYROW, RSRC_A, NPROW ), 00125 * IACOL = INDXG2P( JA, NB_A, MYCOL, CSRC_A, NPCOL ), 00126 * Mp0 = NUMROC( M+IROFF, MB_A, MYROW, IAROW, NPROW ), 00127 * Nq0 = NUMROC( N+ICOFF, NB_A, MYCOL, IACOL, NPCOL ), 00128 * 00129 * and NUMROC, INDXG2P are ScaLAPACK tool functions; 00130 * MYROW, MYCOL, NPROW and NPCOL can be determined by calling 00131 * the subroutine BLACS_GRIDINFO. 00132 * 00133 * If LWORK = -1, then LWORK is global input and a workspace 00134 * query is assumed; the routine only calculates the minimum 00135 * and optimal size for all work arrays. Each of these 00136 * values is returned in the first entry of the corresponding 00137 * work array, and no error message is issued by PXERBLA. 00138 * 00139 * INFO (local output) INTEGER 00140 * = 0: successful exit 00141 * < 0: If the i-th argument is an array and the j-entry had 00142 * an illegal value, then INFO = -(i*100+j), if the i-th 00143 * argument is a scalar and had an illegal value, then 00144 * INFO = -i. 00145 * 00146 * Further Details 00147 * =============== 00148 * 00149 * The matrix Q is represented as a product of elementary reflectors 00150 * 00151 * Q = H(ia+k-1) H(ia+k-2) . . . H(ia), where k = min(m,n). 00152 * 00153 * Each H(i) has the form 00154 * 00155 * H(i) = I - tau * v * v' 00156 * 00157 * where tau is a real scalar, and v is a real vector with v(1:i-1)=0 00158 * and v(i) = 1; v(i+1:n) is stored on exit in A(ia+i-1,ja+i:ja+n-1), 00159 * and tau in TAU(ia+i-1). 00160 * 00161 * ===================================================================== 00162 * 00163 * .. Parameters .. 00164 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_, 00165 $ LLD_, MB_, M_, NB_, N_, RSRC_ 00166 PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1, 00167 $ CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, 00168 $ RSRC_ = 7, CSRC_ = 8, LLD_ = 9 ) 00169 DOUBLE PRECISION ONE 00170 PARAMETER ( ONE = 1.0D+0 ) 00171 * .. 00172 * .. Local Scalars .. 00173 LOGICAL LQUERY 00174 CHARACTER COLBTOP, ROWBTOP 00175 INTEGER IACOL, IAROW, I, ICTXT, J, K, LWMIN, MP, MYCOL, 00176 $ MYROW, NPCOL, NPROW, NQ 00177 DOUBLE PRECISION AII 00178 * .. 00179 * .. External Subroutines .. 00180 EXTERNAL BLACS_ABORT, BLACS_GRIDINFO, CHK1MAT, PDELSET, 00181 $ PDLARF, PDLARFG, PB_TOPGET, PB_TOPSET, PXERBLA 00182 * .. 00183 * .. External Functions .. 00184 INTEGER INDXG2P, NUMROC 00185 EXTERNAL INDXG2P, NUMROC 00186 * .. 00187 * .. Intrinsic Functions .. 00188 INTRINSIC DBLE, MAX, MIN, MOD 00189 * .. 00190 * .. Executable Statements .. 00191 * 00192 * Get grid parameters 00193 * 00194 ICTXT = DESCA( CTXT_ ) 00195 CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL ) 00196 * 00197 * Test the input parameters 00198 * 00199 INFO = 0 00200 IF( NPROW.EQ.-1 ) THEN 00201 INFO = -(600+CTXT_) 00202 ELSE 00203 CALL CHK1MAT( M, 1, N, 2, IA, JA, DESCA, 6, INFO ) 00204 IF( INFO.EQ.0 ) THEN 00205 IAROW = INDXG2P( IA, DESCA( MB_ ), MYROW, DESCA( RSRC_ ), 00206 $ NPROW ) 00207 IACOL = INDXG2P( JA, DESCA( NB_ ), MYCOL, DESCA( CSRC_ ), 00208 $ NPCOL ) 00209 MP = NUMROC( M+MOD( IA-1, DESCA( MB_ ) ), DESCA( MB_ ), 00210 $ MYROW, IAROW, NPROW ) 00211 NQ = NUMROC( N+MOD( JA-1, DESCA( NB_ ) ), DESCA( NB_ ), 00212 $ MYCOL, IACOL, NPCOL ) 00213 LWMIN = NQ + MAX( 1, MP ) 00214 * 00215 WORK( 1 ) = DBLE( LWMIN ) 00216 LQUERY = ( LWORK.EQ.-1 ) 00217 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) 00218 $ INFO = -9 00219 END IF 00220 END IF 00221 * 00222 IF( INFO.NE.0 ) THEN 00223 CALL PXERBLA( ICTXT, 'PDGELQ2', -INFO ) 00224 CALL BLACS_ABORT( ICTXT, 1 ) 00225 RETURN 00226 ELSE IF( LQUERY ) THEN 00227 RETURN 00228 END IF 00229 * 00230 * Quick return if possible 00231 * 00232 IF( M.EQ.0 .OR. N.EQ.0 ) 00233 $ RETURN 00234 * 00235 CALL PB_TOPGET( ICTXT, 'Broadcast', 'Rowwise', ROWBTOP ) 00236 CALL PB_TOPGET( ICTXT, 'Broadcast', 'Columnwise', COLBTOP ) 00237 CALL PB_TOPSET( ICTXT, 'Broadcast', 'Rowwise', ' ' ) 00238 CALL PB_TOPSET( ICTXT, 'Broadcast', 'Columnwise', 'I-ring' ) 00239 * 00240 K = MIN( M, N ) 00241 DO 10 I = IA, IA+K-1 00242 J = JA + I - IA 00243 * 00244 * Generate elementary reflector H(i) to annihilate 00245 * A(i,j+1:ja+n-1) 00246 * 00247 CALL PDLARFG( N-J+JA, AII, I, J, A, I, MIN( J+1, JA+N-1 ), 00248 $ DESCA, DESCA( M_ ), TAU ) 00249 * 00250 IF( I.LT.IA+M-1 ) THEN 00251 * 00252 * Apply H(i) to A(i+1:ia+m-1,j:ja+n-1) from the right 00253 * 00254 CALL PDELSET( A, I, J, DESCA, ONE ) 00255 CALL PDLARF( 'Right', M-I+IA-1, N-J+JA, A, I, J, DESCA, 00256 $ DESCA( M_ ), TAU, A, I+1, J, DESCA, WORK ) 00257 END IF 00258 CALL PDELSET( A, I, J, DESCA, AII ) 00259 * 00260 10 CONTINUE 00261 * 00262 CALL PB_TOPSET( ICTXT, 'Broadcast', 'Rowwise', ROWBTOP ) 00263 CALL PB_TOPSET( ICTXT, 'Broadcast', 'Columnwise', COLBTOP ) 00264 * 00265 WORK( 1 ) = DBLE( LWMIN ) 00266 * 00267 RETURN 00268 * 00269 * End of PDGELQ2 00270 * 00271 END
