|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
00001 SUBROUTINE PCGELS( TRANS, M, N, NRHS, A, IA, JA, DESCA, B, IB, JB, 00002 $ DESCB, 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 1, 1997 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER TRANS 00011 INTEGER IA, IB, INFO, JA, JB, LWORK, M, N, NRHS 00012 * .. 00013 * .. Array Arguments .. 00014 INTEGER DESCA( * ), DESCB( * ) 00015 COMPLEX A( * ), B( * ), WORK( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * PCGELS solves overdetermined or underdetermined complex linear 00022 * systems involving an M-by-N matrix sub( A ) = A(IA:IA+M-1,JA:JA+N-1), 00023 * or its conjugate-transpose, using a QR or LQ factorization of 00024 * sub( A ). It is assumed that sub( A ) has full rank. 00025 * 00026 * The following options are provided: 00027 * 00028 * 1. If TRANS = 'N' and m >= n: find the least squares solution of 00029 * an overdetermined system, i.e., solve the least squares problem 00030 * minimize || sub( B ) - sub( A )*X ||. 00031 * 00032 * 2. If TRANS = 'N' and m < n: find the minimum norm solution of 00033 * an underdetermined system sub( A ) * X = sub( B ). 00034 * 00035 * 3. If TRANS = 'C' and m >= n: find the minimum norm solution of 00036 * an undetermined system sub( A )**H * X = sub( B ). 00037 * 00038 * 4. If TRANS = 'C' and m < n: find the least squares solution of 00039 * an overdetermined system, i.e., solve the least squares problem 00040 * minimize || sub( B ) - sub( A )**H * X ||. 00041 * 00042 * where sub( B ) denotes B( IB:IB+M-1, JB:JB+NRHS-1 ) when TRANS = 'N' 00043 * and B( IB:IB+N-1, JB:JB+NRHS-1 ) otherwise. Several right hand side 00044 * vectors b and solution vectors x can be handled in a single call; 00045 * When TRANS = 'N', the solution vectors are stored as the columns of 00046 * the N-by-NRHS right hand side matrix sub( B ) and the M-by-NRHS 00047 * right hand side matrix sub( B ) otherwise. 00048 * 00049 * Notes 00050 * ===== 00051 * 00052 * Each global data object is described by an associated description 00053 * vector. This vector stores the information required to establish 00054 * the mapping between an object element and its corresponding process 00055 * and memory location. 00056 * 00057 * Let A be a generic term for any 2D block cyclicly distributed array. 00058 * Such a global array has an associated description vector DESCA. 00059 * In the following comments, the character _ should be read as 00060 * "of the global array". 00061 * 00062 * NOTATION STORED IN EXPLANATION 00063 * --------------- -------------- -------------------------------------- 00064 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case, 00065 * DTYPE_A = 1. 00066 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating 00067 * the BLACS process grid A is distribu- 00068 * ted over. The context itself is glo- 00069 * bal, but the handle (the integer 00070 * value) may vary. 00071 * M_A (global) DESCA( M_ ) The number of rows in the global 00072 * array A. 00073 * N_A (global) DESCA( N_ ) The number of columns in the global 00074 * array A. 00075 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute 00076 * the rows of the array. 00077 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute 00078 * the columns of the array. 00079 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first 00080 * row of the array A is distributed. 00081 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the 00082 * first column of the array A is 00083 * distributed. 00084 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local 00085 * array. LLD_A >= MAX(1,LOCr(M_A)). 00086 * 00087 * Let K be the number of rows or columns of a distributed matrix, 00088 * and assume that its process grid has dimension p x q. 00089 * LOCr( K ) denotes the number of elements of K that a process 00090 * would receive if K were distributed over the p processes of its 00091 * process column. 00092 * Similarly, LOCc( K ) denotes the number of elements of K that a 00093 * process would receive if K were distributed over the q processes of 00094 * its process row. 00095 * The values of LOCr() and LOCc() may be determined via a call to the 00096 * ScaLAPACK tool function, NUMROC: 00097 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ), 00098 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ). 00099 * An upper bound for these quantities may be computed by: 00100 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A 00101 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A 00102 * 00103 * Arguments 00104 * ========= 00105 * 00106 * TRANS (global input) CHARACTER 00107 * = 'N': the linear system involves sub( A ); 00108 * = 'C': the linear system involves sub( A )**H. 00109 * 00110 * M (global input) INTEGER 00111 * The number of rows to be operated on, i.e. the number of 00112 * rows of the distributed submatrix sub( A ). M >= 0. 00113 * 00114 * N (global input) INTEGER 00115 * The number of columns to be operated on, i.e. the number of 00116 * columns of the distributed submatrix sub( A ). N >= 0. 00117 * 00118 * NRHS (global input) INTEGER 00119 * The number of right hand sides, i.e. the number of columns 00120 * of the distributed submatrices sub( B ) and X. NRHS >= 0. 00121 * 00122 * A (local input/local output) COMPLEX pointer into the 00123 * local memory to an array of local dimension 00124 * ( LLD_A, LOCc(JA+N-1) ). On entry, the M-by-N matrix A. 00125 * if M >= N, sub( A ) is overwritten by details of its QR 00126 * factorization as returned by PCGEQRF; 00127 * if M < N, sub( A ) is overwritten by details of its LQ 00128 * factorization as returned by PCGELQF. 00129 * 00130 * IA (global input) INTEGER 00131 * The row index in the global array A indicating the first 00132 * row of sub( A ). 00133 * 00134 * JA (global input) INTEGER 00135 * The column index in the global array A indicating the 00136 * first column of sub( A ). 00137 * 00138 * DESCA (global and local input) INTEGER array of dimension DLEN_. 00139 * The array descriptor for the distributed matrix A. 00140 * 00141 * B (local input/local output) COMPLEX pointer into the 00142 * local memory to an array of local dimension 00143 * (LLD_B, LOCc(JB+NRHS-1)). On entry, this array contains the 00144 * local pieces of the distributed matrix B of right hand side 00145 * vectors, stored columnwise; 00146 * sub( B ) is M-by-NRHS if TRANS='N', and N-by-NRHS otherwise. 00147 * On exit, sub( B ) is overwritten by the solution vectors, 00148 * stored columnwise: if TRANS = 'N' and M >= N, rows 1 to N 00149 * of sub( B ) contain the least squares solution vectors; the 00150 * residual sum of squares for the solution in each column is 00151 * given by the sum of squares of elements N+1 to M in that 00152 * column; if TRANS = 'N' and M < N, rows 1 to N of sub( B ) 00153 * contain the minimum norm solution vectors; if TRANS = 'C' 00154 * and M >= N, rows 1 to M of sub( B ) contain the minimum norm 00155 * solution vectors; if TRANS = 'C' and M < N, rows 1 to M of 00156 * sub( B ) contain the least squares solution vectors; the 00157 * residual sum of squares for the solution in each column is 00158 * given by the sum of squares of elements M+1 to N in that 00159 * column. 00160 * 00161 * IB (global input) INTEGER 00162 * The row index in the global array B indicating the first 00163 * row of sub( B ). 00164 * 00165 * JB (global input) INTEGER 00166 * The column index in the global array B indicating the 00167 * first column of sub( B ). 00168 * 00169 * DESCB (global and local input) INTEGER array of dimension DLEN_. 00170 * The array descriptor for the distributed matrix B. 00171 * 00172 * WORK (local workspace/local output) COMPLEX array, 00173 * dimension (LWORK) 00174 * On exit, WORK(1) returns the minimal and optimal LWORK. 00175 * 00176 * LWORK (local or global input) INTEGER 00177 * The dimension of the array WORK. 00178 * LWORK is local input and must be at least 00179 * LWORK >= LTAU + MAX( LWF, LWS ) where 00180 * If M >= N, then 00181 * LTAU = NUMROC( JA+MIN(M,N)-1, NB_A, MYCOL, CSRC_A, NPCOL ), 00182 * LWF = NB_A * ( MpA0 + NqA0 + NB_A ) 00183 * LWS = MAX( (NB_A*(NB_A-1))/2, (NRHSqB0 + MpB0)*NB_A ) + 00184 * NB_A * NB_A 00185 * Else 00186 * LTAU = NUMROC( IA+MIN(M,N)-1, MB_A, MYROW, RSRC_A, NPROW ), 00187 * LWF = MB_A * ( MpA0 + NqA0 + MB_A ) 00188 * LWS = MAX( (MB_A*(MB_A-1))/2, ( NpB0 + MAX( NqA0 + 00189 * NUMROC( NUMROC( N+IROFFB, MB_A, 0, 0, NPROW ), 00190 * MB_A, 0, 0, LCMP ), NRHSqB0 ) )*MB_A ) + 00191 * MB_A * MB_A 00192 * End if 00193 * 00194 * where LCMP = LCM / NPROW with LCM = ILCM( NPROW, NPCOL ), 00195 * 00196 * IROFFA = MOD( IA-1, MB_A ), ICOFFA = MOD( JA-1, NB_A ), 00197 * IAROW = INDXG2P( IA, MB_A, MYROW, RSRC_A, NPROW ), 00198 * IACOL = INDXG2P( JA, NB_A, MYCOL, CSRC_A, NPCOL ), 00199 * MpA0 = NUMROC( M+IROFFA, MB_A, MYROW, IAROW, NPROW ), 00200 * NqA0 = NUMROC( N+ICOFFA, NB_A, MYCOL, IACOL, NPCOL ), 00201 * 00202 * IROFFB = MOD( IB-1, MB_B ), ICOFFB = MOD( JB-1, NB_B ), 00203 * IBROW = INDXG2P( IB, MB_B, MYROW, RSRC_B, NPROW ), 00204 * IBCOL = INDXG2P( JB, NB_B, MYCOL, CSRC_B, NPCOL ), 00205 * MpB0 = NUMROC( M+IROFFB, MB_B, MYROW, IBROW, NPROW ), 00206 * NpB0 = NUMROC( N+IROFFB, MB_B, MYROW, IBROW, NPROW ), 00207 * NRHSqB0 = NUMROC( NRHS+ICOFFB, NB_B, MYCOL, IBCOL, NPCOL ), 00208 * 00209 * ILCM, INDXG2P and NUMROC are ScaLAPACK tool functions; 00210 * MYROW, MYCOL, NPROW and NPCOL can be determined by calling 00211 * the subroutine BLACS_GRIDINFO. 00212 * 00213 * If LWORK = -1, then LWORK is global input and a workspace 00214 * query is assumed; the routine only calculates the minimum 00215 * and optimal size for all work arrays. Each of these 00216 * values is returned in the first entry of the corresponding 00217 * work array, and no error message is issued by PXERBLA. 00218 * 00219 * INFO (global output) INTEGER 00220 * = 0: successful exit 00221 * < 0: If the i-th argument is an array and the j-entry had 00222 * an illegal value, then INFO = -(i*100+j), if the i-th 00223 * argument is a scalar and had an illegal value, then 00224 * INFO = -i. 00225 * 00226 * ===================================================================== 00227 * 00228 * .. Parameters .. 00229 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_, 00230 $ LLD_, MB_, M_, NB_, N_, RSRC_ 00231 PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1, 00232 $ CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, 00233 $ RSRC_ = 7, CSRC_ = 8, LLD_ = 9 ) 00234 REAL ZERO, ONE 00235 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00236 COMPLEX CZERO, CONE 00237 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 00238 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 00239 * .. 00240 * .. Local Scalars .. 00241 LOGICAL LQUERY, TPSD 00242 INTEGER BROW, IACOL, IAROW, IASCL, IBCOL, IBROW, IBSCL, 00243 $ ICOFFA, ICOFFB, ICTXT, IPW, IROFFA, IROFFB, 00244 $ LCM, LCMP, LTAU, LWF, LWMIN, LWS, MPA0, MPB0, 00245 $ MYCOL, MYROW, NPB0, NPCOL, NPROW, NQA0, 00246 $ NRHSQB0, SCLLEN 00247 REAL ANRM, BIGNUM, BNRM, SMLNUM 00248 * .. 00249 * .. Local Arrays .. 00250 INTEGER IDUM1( 2 ), IDUM2( 2 ) 00251 REAL RWORK( 1 ) 00252 * .. 00253 * .. External Functions .. 00254 LOGICAL LSAME 00255 INTEGER ILCM 00256 INTEGER INDXG2P, NUMROC 00257 REAL PCLANGE, PSLAMCH 00258 EXTERNAL ILCM, INDXG2P, LSAME, NUMROC, PCLANGE, 00259 $ PSLAMCH 00260 * .. 00261 * .. External Subroutines .. 00262 EXTERNAL BLACS_GRIDINFO, CHK1MAT, PCHK2MAT, PCGELQF, 00263 $ PCGEQRF, PSLABAD, PCLASCL, PCLASET, 00264 $ PCTRSM, PCUNMLQ, PCUNMQR, PXERBLA 00265 * .. 00266 * .. Intrinsic Functions .. 00267 INTRINSIC CMPLX, ICHAR, MAX, MIN, MOD, REAL 00268 * .. 00269 * .. Executable Statements .. 00270 * 00271 * Get grid parameters 00272 * 00273 ICTXT = DESCA( CTXT_ ) 00274 CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL ) 00275 * 00276 * Test the input parameters 00277 * 00278 INFO = 0 00279 IF( NPROW.EQ.-1 ) THEN 00280 INFO = -( 800 + CTXT_ ) 00281 ELSE 00282 CALL CHK1MAT( M, 2, N, 3, IA, JA, DESCA, 8, INFO ) 00283 IF ( M .GE. N ) THEN 00284 CALL CHK1MAT( M, 2, NRHS, 4, IB, JB, DESCB, 12, INFO ) 00285 ELSE 00286 CALL CHK1MAT( N, 3, NRHS, 4, IB, JB, DESCB, 12, INFO ) 00287 ENDIF 00288 IF( INFO.EQ.0 ) THEN 00289 IROFFA = MOD( IA-1, DESCA( MB_ ) ) 00290 ICOFFA = MOD( JA-1, DESCA( NB_ ) ) 00291 IROFFB = MOD( IB-1, DESCB( MB_ ) ) 00292 ICOFFB = MOD( JB-1, DESCB( NB_ ) ) 00293 IAROW = INDXG2P( IA, DESCA( MB_ ), MYROW, DESCA( RSRC_ ), 00294 $ NPROW ) 00295 IACOL = INDXG2P( IA, DESCA( NB_ ), MYCOL, DESCA( CSRC_ ), 00296 $ NPCOL ) 00297 MPA0 = NUMROC( M+IROFFA, DESCA( MB_ ), MYROW, IAROW, NPROW ) 00298 NQA0 = NUMROC( N+ICOFFA, DESCA( NB_ ), MYCOL, IACOL, NPCOL ) 00299 * 00300 IBROW = INDXG2P( IB, DESCB( MB_ ), MYROW, DESCB( RSRC_ ), 00301 $ NPROW ) 00302 IBCOL = INDXG2P( IB, DESCB( NB_ ), MYCOL, DESCB( CSRC_ ), 00303 $ NPCOL ) 00304 NRHSQB0 = NUMROC( NRHS+ICOFFB, DESCB( NB_ ), MYCOL, IBCOL, 00305 $ NPCOL ) 00306 IF( M.GE.N ) THEN 00307 MPB0 = NUMROC( M+IROFFB, DESCB( MB_ ), MYROW, IBROW, 00308 $ NPROW ) 00309 LTAU = NUMROC( JA+MIN(M,N)-1, DESCA( NB_ ), MYCOL, 00310 $ DESCA( CSRC_ ), NPCOL ) 00311 LWF = DESCA( NB_ ) * ( MPA0 + NQA0 + DESCA( NB_ ) ) 00312 LWS = MAX( ( DESCA( NB_ )*( DESCA( NB_ ) - 1 ) ) / 2, 00313 $ ( MPB0 + NRHSQB0 ) * DESCA( NB_ ) ) + 00314 $ DESCA( NB_ )*DESCA( NB_ ) 00315 ELSE 00316 LCM = ILCM( NPROW, NPCOL ) 00317 LCMP = LCM / NPROW 00318 NPB0 = NUMROC( N+IROFFB, DESCB( MB_ ), MYROW, IBROW, 00319 $ NPROW ) 00320 LTAU = NUMROC( IA+MIN(M,N)-1, DESCA( MB_ ), MYROW, 00321 $ DESCA( RSRC_ ), NPROW ) 00322 LWF = DESCA( MB_ ) * ( MPA0 + NQA0 + DESCA( MB_ ) ) 00323 LWS = MAX( ( DESCA( MB_ )*( DESCA( MB_ ) - 1 ) ) / 2, 00324 $ ( NPB0 + MAX( NQA0 + NUMROC( NUMROC( N+IROFFB, 00325 $ DESCA( MB_ ), 0, 0, NPROW ), DESCA( MB_ ), 0, 0, 00326 $ LCMP ), NRHSQB0 ) )*DESCA( MB_ ) ) + 00327 $ DESCA( MB_ ) * DESCA( MB_ ) 00328 END IF 00329 LWMIN = LTAU + MAX( LWF, LWS ) 00330 WORK( 1 ) = CMPLX( REAL( LWMIN ) ) 00331 LQUERY = ( LWORK.EQ.-1 ) 00332 * 00333 TPSD = .TRUE. 00334 IF( LSAME( TRANS, 'N' ) ) 00335 $ TPSD = .FALSE. 00336 * 00337 IF( .NOT.( LSAME( TRANS, 'N' ) .OR. 00338 $ LSAME( TRANS, 'C' ) ) ) THEN 00339 INFO = -1 00340 ELSE IF( M.LT.0 ) THEN 00341 INFO = -2 00342 ELSE IF( N.LT.0 ) THEN 00343 INFO = -3 00344 ELSE IF( NRHS.LT.0 ) THEN 00345 INFO = -4 00346 ELSE IF( M.GE.N .AND. IROFFA.NE.IROFFB ) THEN 00347 INFO = -10 00348 ELSE IF( M.GE.N .AND. IAROW.NE.IBROW ) THEN 00349 INFO = -10 00350 ELSE IF( M.LT.N .AND. ICOFFA.NE.IROFFB ) THEN 00351 INFO = -10 00352 ELSE IF( M.GE.N .AND. DESCA( MB_ ).NE.DESCB( MB_ ) ) THEN 00353 INFO = -( 1200 + MB_ ) 00354 ELSE IF( M.LT.N .AND. DESCA( NB_ ).NE.DESCB( MB_ ) ) THEN 00355 INFO = -( 1200 + MB_ ) 00356 ELSE IF( ICTXT.NE.DESCB( CTXT_ ) ) THEN 00357 INFO = -( 1200 + CTXT_ ) 00358 ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00359 INFO = -14 00360 END IF 00361 END IF 00362 * 00363 IF( .NOT.TPSD ) THEN 00364 IDUM1( 1 ) = ICHAR( 'N' ) 00365 ELSE 00366 IDUM1( 1 ) = ICHAR( 'C' ) 00367 END IF 00368 IDUM2( 1 ) = 1 00369 IF( LWORK.EQ.-1 ) THEN 00370 IDUM1( 2 ) = -1 00371 ELSE 00372 IDUM1( 2 ) = 1 00373 END IF 00374 IDUM2( 2 ) = 14 00375 CALL PCHK2MAT( M, 2, N, 3, IA, JA, DESCA, 8, N, 3, NRHS, 4, 00376 $ IB, JB, DESCB, 12, 2, IDUM1, IDUM2, INFO ) 00377 END IF 00378 * 00379 IF( INFO.NE.0 ) THEN 00380 CALL PXERBLA( ICTXT, 'PCGELS', -INFO ) 00381 RETURN 00382 ELSE IF( LQUERY ) THEN 00383 RETURN 00384 END IF 00385 * 00386 * Quick return if possible 00387 * 00388 IF( MIN( M, N, NRHS ).EQ.0 ) THEN 00389 CALL PCLASET( 'Full', MAX( M, N ), NRHS, CZERO, CZERO, B, 00390 $ IB, JB, DESCB ) 00391 RETURN 00392 END IF 00393 * 00394 * Get machine parameters 00395 * 00396 SMLNUM = PSLAMCH( ICTXT, 'S' ) 00397 SMLNUM = SMLNUM / PSLAMCH( ICTXT, 'P' ) 00398 BIGNUM = ONE / SMLNUM 00399 CALL PSLABAD( ICTXT, SMLNUM, BIGNUM ) 00400 * 00401 * Scale A, B if max entry outside range [SMLNUM,BIGNUM] 00402 * 00403 ANRM = PCLANGE( 'M', M, N, A, IA, JA, DESCA, RWORK ) 00404 IASCL = 0 00405 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00406 * 00407 * Scale matrix norm up to SMLNUM 00408 * 00409 CALL PCLASCL( 'G', ANRM, SMLNUM, M, N, A, IA, JA, DESCA, 00410 $ INFO ) 00411 IASCL = 1 00412 ELSE IF( ANRM.GT.BIGNUM ) THEN 00413 * 00414 * Scale matrix norm down to BIGNUM 00415 * 00416 CALL PCLASCL( 'G', ANRM, BIGNUM, M, N, A, IA, JA, DESCA, 00417 $ INFO ) 00418 IASCL = 2 00419 ELSE IF( ANRM.EQ.ZERO ) THEN 00420 * 00421 * Matrix all zero. Return zero solution. 00422 * 00423 CALL PCLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, IB, 00424 $ JB, DESCB ) 00425 GO TO 10 00426 END IF 00427 * 00428 BROW = M 00429 IF( TPSD ) 00430 $ BROW = N 00431 * 00432 BNRM = PCLANGE( 'M', BROW, NRHS, B, IB, JB, DESCB, RWORK ) 00433 * 00434 IBSCL = 0 00435 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 00436 * 00437 * Scale matrix norm up to SMLNUM 00438 * 00439 CALL PCLASCL( 'G', BNRM, SMLNUM, BROW, NRHS, B, IB, JB, 00440 $ DESCB, INFO ) 00441 IBSCL = 1 00442 ELSE IF( BNRM.GT.BIGNUM ) THEN 00443 * 00444 * Scale matrix norm down to BIGNUM 00445 * 00446 CALL PCLASCL( 'G', BNRM, BIGNUM, BROW, NRHS, B, IB, JB, 00447 $ DESCB, INFO ) 00448 IBSCL = 2 00449 END IF 00450 * 00451 IPW = LTAU + 1 00452 * 00453 IF( M.GE.N ) THEN 00454 * 00455 * compute QR factorization of A 00456 * 00457 CALL PCGEQRF( M, N, A, IA, JA, DESCA, WORK, WORK( IPW ), 00458 $ LWORK-LTAU, INFO ) 00459 * 00460 * workspace at least N, optimally N*NB 00461 * 00462 IF( .NOT.TPSD ) THEN 00463 * 00464 * Least-Squares Problem min || A * X - B || 00465 * 00466 * B(IB:IB+M-1,JB:JB+NRHS-1) := Q' * B(IB:IB+M-1,JB:JB+NRHS-1) 00467 * 00468 CALL PCUNMQR( 'Left', 'Conjugate transpose', M, NRHS, N, A, 00469 $ IA, JA, DESCA, WORK, B, IB, JB, DESCB, 00470 $ WORK( IPW ), LWORK-LTAU, INFO ) 00471 * 00472 * workspace at least NRHS, optimally NRHS*NB 00473 * 00474 * B(IB:IB+N-1,JB:JB+NRHS-1) := inv(R) * 00475 * B(IB:IB+N-1,JB:JB+NRHS-1) 00476 * 00477 CALL PCTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', N, 00478 $ NRHS, CONE, A, IA, JA, DESCA, B, IB, JB, 00479 $ DESCB ) 00480 * 00481 SCLLEN = N 00482 * 00483 ELSE 00484 * 00485 * Overdetermined system of equations sub( A )' * X = sub( B ) 00486 * 00487 * sub( B ) := inv(R') * sub( B ) 00488 * 00489 CALL PCTRSM( 'Left', 'Upper', 'Conjugate transpose', 00490 $ 'Non-unit', N, NRHS, CONE, A, IA, JA, DESCA, 00491 $ B, IB, JB, DESCB ) 00492 * 00493 * B(IB+N:IB+M-1,JB:JB+NRHS-1) = ZERO 00494 * 00495 CALL PCLASET( 'All', M-N, NRHS, CZERO, CZERO, B, IB+N, JB, 00496 $ DESCB ) 00497 * 00498 * B(IB:IB+M-1,JB:JB+NRHS-1) := Q(1:N,:) * 00499 * B(IB:IB+N-1,JB:JB+NRHS-1) 00500 * 00501 CALL PCUNMQR( 'Left', 'No transpose', M, NRHS, N, A, IA, JA, 00502 $ DESCA, WORK, B, IB, JB, DESCB, WORK( IPW ), 00503 $ LWORK-LTAU, INFO ) 00504 * 00505 * workspace at least NRHS, optimally NRHS*NB 00506 * 00507 SCLLEN = M 00508 * 00509 END IF 00510 * 00511 ELSE 00512 * 00513 * Compute LQ factorization of sub( A ) 00514 * 00515 CALL PCGELQF( M, N, A, IA, JA, DESCA, WORK, WORK( IPW ), 00516 $ LWORK-LTAU, INFO ) 00517 * 00518 * workspace at least M, optimally M*NB. 00519 * 00520 IF( .NOT.TPSD ) THEN 00521 * 00522 * underdetermined system of equations sub( A ) * X = sub( B ) 00523 * 00524 * B(IB:IB+M-1,JB:JB+NRHS-1) := inv(L) * 00525 * B(IB:IB+M-1,JB:JB+NRHS-1) 00526 * 00527 CALL PCTRSM( 'Left', 'Lower', 'No transpose', 'Non-unit', M, 00528 $ NRHS, CONE, A, IA, JA, DESCA, B, IB, JB, 00529 $ DESCB ) 00530 * 00531 * B(IB+M:IB+N-1,JB:JB+NRHS-1) = 0 00532 * 00533 CALL PCLASET( 'All', N-M, NRHS, CZERO, CZERO, B, IB+M, JB, 00534 $ DESCB ) 00535 * 00536 * B(IB:IB+N-1,JB:JB+NRHS-1) := Q(1:N,:)' * 00537 * B(IB:IB+M-1,JB:JB+NRHS-1) 00538 * 00539 CALL PCUNMLQ( 'Left', 'Conjugate transpose', N, NRHS, M, A, 00540 $ IA, JA, DESCA, WORK, B, IB, JB, DESCB, 00541 $ WORK( IPW ), LWORK-LTAU, INFO ) 00542 * 00543 * workspace at least NRHS, optimally NRHS*NB 00544 * 00545 SCLLEN = N 00546 * 00547 ELSE 00548 * 00549 * overdetermined system min || A' * X - B || 00550 * 00551 * B(IB:IB+N-1,JB:JB+NRHS-1) := Q * B(IB:IB+N-1,JB:JB+NRHS-1) 00552 * 00553 CALL PCUNMLQ( 'Left', 'No transpose', N, NRHS, M, A, IA, JA, 00554 $ DESCA, WORK, B, IB, JB, DESCB, WORK( IPW ), 00555 $ LWORK-LTAU, INFO ) 00556 * 00557 * workspace at least NRHS, optimally NRHS*NB 00558 * 00559 * B(IB:IB+M-1,JB:JB+NRHS-1) := inv(L') * 00560 * B(IB:IB+M-1,JB:JB+NRHS-1) 00561 * 00562 CALL PCTRSM( 'Left', 'Lower', 'Conjugate transpose', 00563 $ 'Non-unit', M, NRHS, CONE, A, IA, JA, DESCA, 00564 $ B, IB, JB, DESCB ) 00565 * 00566 SCLLEN = M 00567 * 00568 END IF 00569 * 00570 END IF 00571 * 00572 * Undo scaling 00573 * 00574 IF( IASCL.EQ.1 ) THEN 00575 CALL PCLASCL( 'G', ANRM, SMLNUM, SCLLEN, NRHS, B, IB, JB, 00576 $ DESCB, INFO ) 00577 ELSE IF( IASCL.EQ.2 ) THEN 00578 CALL PCLASCL( 'G', ANRM, BIGNUM, SCLLEN, NRHS, B, IB, JB, 00579 $ DESCB, INFO ) 00580 END IF 00581 IF( IBSCL.EQ.1 ) THEN 00582 CALL PCLASCL( 'G', SMLNUM, BNRM, SCLLEN, NRHS, B, IB, JB, 00583 $ DESCB, INFO ) 00584 ELSE IF( IBSCL.EQ.2 ) THEN 00585 CALL PCLASCL( 'G', BIGNUM, BNRM, SCLLEN, NRHS, B, IB, JB, 00586 $ DESCB, INFO ) 00587 END IF 00588 * 00589 10 CONTINUE 00590 * 00591 WORK( 1 ) = CMPLX( REAL( LWMIN ) ) 00592 * 00593 RETURN 00594 * 00595 * End of PCGELS 00596 * 00597 END
