|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
00001 SUBROUTINE PDGELS( 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 DOUBLE PRECISION A( * ), B( * ), WORK( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * PDGELS solves overdetermined or underdetermined real linear 00022 * systems involving an M-by-N matrix sub( A ) = A(IA:IA+M-1,JA:JA+N-1), 00023 * or its transpose, using a QR or LQ factorization of sub( A ). It is 00024 * 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 = 'T' and m >= n: find the minimum norm solution of 00036 * an undetermined system sub( A )**T * X = sub( B ). 00037 * 00038 * 4. If TRANS = 'T' 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 )**T * 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 * = 'T': the linear system involves sub( A )**T. 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) DOUBLE PRECISION 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 PDGEQRF; 00127 * if M < N, sub( A ) is overwritten by details of its LQ 00128 * factorization as returned by PDGELQF. 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) DOUBLE PRECISION 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 = 'T' 00154 * and M >= N, rows 1 to M of sub( B ) contain the minimum norm 00155 * solution vectors; if TRANS = 'T' 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) DOUBLE PRECISION 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 DOUBLE PRECISION ZERO, ONE 00235 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00236 * .. 00237 * .. Local Scalars .. 00238 LOGICAL LQUERY, TPSD 00239 INTEGER BROW, IACOL, IAROW, IASCL, IBCOL, IBROW, IBSCL, 00240 $ ICOFFA, ICOFFB, ICTXT, IPW, IROFFA, IROFFB, 00241 $ LCM, LCMP, LTAU, LWF, LWMIN, LWS, MPA0, MPB0, 00242 $ MYCOL, MYROW, NPB0, NPCOL, NPROW, NQA0, 00243 $ NRHSQB0, SCLLEN 00244 DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM 00245 * .. 00246 * .. Local Arrays .. 00247 INTEGER IDUM1( 2 ), IDUM2( 2 ) 00248 DOUBLE PRECISION RWORK( 1 ) 00249 * .. 00250 * .. External Functions .. 00251 LOGICAL LSAME 00252 INTEGER ILCM 00253 INTEGER INDXG2P, NUMROC 00254 DOUBLE PRECISION PDLAMCH, PDLANGE 00255 EXTERNAL ILCM, INDXG2P, LSAME, NUMROC, PDLAMCH, 00256 $ PDLANGE 00257 * .. 00258 * .. External Subroutines .. 00259 EXTERNAL BLACS_GRIDINFO, CHK1MAT, PCHK2MAT, PDGELQF, 00260 $ PDGEQRF, PDLABAD, PDLASCL, PDLASET, 00261 $ PDORMLQ, PDORMQR, PDTRSM, PXERBLA 00262 * .. 00263 * .. Intrinsic Functions .. 00264 INTRINSIC DBLE, ICHAR, MAX, MIN, MOD 00265 * .. 00266 * .. Executable Statements .. 00267 * 00268 * Get grid parameters 00269 * 00270 ICTXT = DESCA( CTXT_ ) 00271 CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL ) 00272 * 00273 * Test the input parameters 00274 * 00275 INFO = 0 00276 IF( NPROW.EQ.-1 ) THEN 00277 INFO = -( 800 + CTXT_ ) 00278 ELSE 00279 CALL CHK1MAT( M, 2, N, 3, IA, JA, DESCA, 8, INFO ) 00280 IF ( M .GE. N ) THEN 00281 CALL CHK1MAT( M, 2, NRHS, 4, IB, JB, DESCB, 12, INFO ) 00282 ELSE 00283 CALL CHK1MAT( N, 3, NRHS, 4, IB, JB, DESCB, 12, INFO ) 00284 ENDIF 00285 IF( INFO.EQ.0 ) THEN 00286 IROFFA = MOD( IA-1, DESCA( MB_ ) ) 00287 ICOFFA = MOD( JA-1, DESCA( NB_ ) ) 00288 IROFFB = MOD( IB-1, DESCB( MB_ ) ) 00289 ICOFFB = MOD( JB-1, DESCB( NB_ ) ) 00290 IAROW = INDXG2P( IA, DESCA( MB_ ), MYROW, DESCA( RSRC_ ), 00291 $ NPROW ) 00292 IACOL = INDXG2P( IA, DESCA( NB_ ), MYCOL, DESCA( CSRC_ ), 00293 $ NPCOL ) 00294 MPA0 = NUMROC( M+IROFFA, DESCA( MB_ ), MYROW, IAROW, NPROW ) 00295 NQA0 = NUMROC( N+ICOFFA, DESCA( NB_ ), MYCOL, IACOL, NPCOL ) 00296 * 00297 IBROW = INDXG2P( IB, DESCB( MB_ ), MYROW, DESCB( RSRC_ ), 00298 $ NPROW ) 00299 IBCOL = INDXG2P( IB, DESCB( NB_ ), MYCOL, DESCB( CSRC_ ), 00300 $ NPCOL ) 00301 NRHSQB0 = NUMROC( NRHS+ICOFFB, DESCB( NB_ ), MYCOL, IBCOL, 00302 $ NPCOL ) 00303 IF( M.GE.N ) THEN 00304 MPB0 = NUMROC( M+IROFFB, DESCB( MB_ ), MYROW, IBROW, 00305 $ NPROW ) 00306 LTAU = NUMROC( JA+MIN(M,N)-1, DESCA( NB_ ), MYCOL, 00307 $ DESCA( CSRC_ ), NPCOL ) 00308 LWF = DESCA( NB_ ) * ( MPA0 + NQA0 + DESCA( NB_ ) ) 00309 LWS = MAX( ( DESCA( NB_ )*( DESCA( NB_ ) - 1 ) ) / 2, 00310 $ ( MPB0 + NRHSQB0 ) * DESCA( NB_ ) ) + 00311 $ DESCA( NB_ )*DESCA( NB_ ) 00312 ELSE 00313 LCM = ILCM( NPROW, NPCOL ) 00314 LCMP = LCM / NPROW 00315 NPB0 = NUMROC( N+IROFFB, DESCB( MB_ ), MYROW, IBROW, 00316 $ NPROW ) 00317 LTAU = NUMROC( IA+MIN(M,N)-1, DESCA( MB_ ), MYROW, 00318 $ DESCA( RSRC_ ), NPROW ) 00319 LWF = DESCA( MB_ ) * ( MPA0 + NQA0 + DESCA( MB_ ) ) 00320 LWS = MAX( ( DESCA( MB_ )*( DESCA( MB_ ) - 1 ) ) / 2, 00321 $ ( NPB0 + MAX( NQA0 + NUMROC( NUMROC( N+IROFFB, 00322 $ DESCA( MB_ ), 0, 0, NPROW ), DESCA( MB_ ), 0, 0, 00323 $ LCMP ), NRHSQB0 ) )*DESCA( MB_ ) ) + 00324 $ DESCA( MB_ ) * DESCA( MB_ ) 00325 END IF 00326 LWMIN = LTAU + MAX( LWF, LWS ) 00327 WORK( 1 ) = DBLE( LWMIN ) 00328 LQUERY = ( LWORK.EQ.-1 ) 00329 * 00330 TPSD = .TRUE. 00331 IF( LSAME( TRANS, 'N' ) ) 00332 $ TPSD = .FALSE. 00333 * 00334 IF( .NOT.( LSAME( TRANS, 'N' ) .OR. 00335 $ LSAME( TRANS, 'T' ) ) ) THEN 00336 INFO = -1 00337 ELSE IF( M.LT.0 ) THEN 00338 INFO = -2 00339 ELSE IF( N.LT.0 ) THEN 00340 INFO = -3 00341 ELSE IF( NRHS.LT.0 ) THEN 00342 INFO = -4 00343 ELSE IF( M.GE.N .AND. IROFFA.NE.IROFFB ) THEN 00344 INFO = -10 00345 ELSE IF( M.GE.N .AND. IAROW.NE.IBROW ) THEN 00346 INFO = -10 00347 ELSE IF( M.LT.N .AND. ICOFFA.NE.IROFFB ) THEN 00348 INFO = -10 00349 ELSE IF( M.GE.N .AND. DESCA( MB_ ).NE.DESCB( MB_ ) ) THEN 00350 INFO = -( 1200 + MB_ ) 00351 ELSE IF( M.LT.N .AND. DESCA( NB_ ).NE.DESCB( MB_ ) ) THEN 00352 INFO = -( 1200 + MB_ ) 00353 ELSE IF( ICTXT.NE.DESCB( CTXT_ ) ) THEN 00354 INFO = -( 1200 + CTXT_ ) 00355 ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00356 INFO = -14 00357 END IF 00358 END IF 00359 * 00360 IF( .NOT.TPSD ) THEN 00361 IDUM1( 1 ) = ICHAR( 'N' ) 00362 ELSE 00363 IDUM1( 1 ) = ICHAR( 'T' ) 00364 END IF 00365 IDUM2( 1 ) = 1 00366 IF( LWORK.EQ.-1 ) THEN 00367 IDUM1( 2 ) = -1 00368 ELSE 00369 IDUM1( 2 ) = 1 00370 END IF 00371 IDUM2( 2 ) = 14 00372 CALL PCHK2MAT( M, 2, N, 3, IA, JA, DESCA, 8, N, 3, NRHS, 4, 00373 $ IB, JB, DESCB, 12, 2, IDUM1, IDUM2, INFO ) 00374 END IF 00375 * 00376 IF( INFO.NE.0 ) THEN 00377 CALL PXERBLA( ICTXT, 'PDGELS', -INFO ) 00378 RETURN 00379 ELSE IF( LQUERY ) THEN 00380 RETURN 00381 END IF 00382 * 00383 * Quick return if possible 00384 * 00385 IF( MIN( M, N, NRHS ).EQ.0 ) THEN 00386 CALL PDLASET( 'Full', MAX( M, N ), NRHS, ZERO, ZERO, B, 00387 $ IB, JB, DESCB ) 00388 RETURN 00389 END IF 00390 * 00391 * Get machine parameters 00392 * 00393 SMLNUM = PDLAMCH( ICTXT, 'S' ) 00394 SMLNUM = SMLNUM / PDLAMCH( ICTXT, 'P' ) 00395 BIGNUM = ONE / SMLNUM 00396 CALL PDLABAD( ICTXT, SMLNUM, BIGNUM ) 00397 * 00398 * Scale A, B if max entry outside range [SMLNUM,BIGNUM] 00399 * 00400 ANRM = PDLANGE( 'M', M, N, A, IA, JA, DESCA, RWORK ) 00401 IASCL = 0 00402 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00403 * 00404 * Scale matrix norm up to SMLNUM 00405 * 00406 CALL PDLASCL( 'G', ANRM, SMLNUM, M, N, A, IA, JA, DESCA, 00407 $ INFO ) 00408 IASCL = 1 00409 ELSE IF( ANRM.GT.BIGNUM ) THEN 00410 * 00411 * Scale matrix norm down to BIGNUM 00412 * 00413 CALL PDLASCL( 'G', ANRM, BIGNUM, M, N, A, IA, JA, DESCA, 00414 $ INFO ) 00415 IASCL = 2 00416 ELSE IF( ANRM.EQ.ZERO ) THEN 00417 * 00418 * Matrix all zero. Return zero solution. 00419 * 00420 CALL PDLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, IB, JB, 00421 $ DESCB ) 00422 GO TO 10 00423 END IF 00424 * 00425 BROW = M 00426 IF( TPSD ) 00427 $ BROW = N 00428 * 00429 BNRM = PDLANGE( 'M', BROW, NRHS, B, IB, JB, DESCB, RWORK ) 00430 * 00431 IBSCL = 0 00432 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 00433 * 00434 * Scale matrix norm up to SMLNUM 00435 * 00436 CALL PDLASCL( 'G', BNRM, SMLNUM, BROW, NRHS, B, IB, JB, 00437 $ DESCB, INFO ) 00438 IBSCL = 1 00439 ELSE IF( BNRM.GT.BIGNUM ) THEN 00440 * 00441 * Scale matrix norm down to BIGNUM 00442 * 00443 CALL PDLASCL( 'G', BNRM, BIGNUM, BROW, NRHS, B, IB, JB, 00444 $ DESCB, INFO ) 00445 IBSCL = 2 00446 END IF 00447 * 00448 IPW = LTAU + 1 00449 * 00450 IF( M.GE.N ) THEN 00451 * 00452 * compute QR factorization of A 00453 * 00454 CALL PDGEQRF( M, N, A, IA, JA, DESCA, WORK, WORK( IPW ), 00455 $ LWORK-LTAU, INFO ) 00456 * 00457 * workspace at least N, optimally N*NB 00458 * 00459 IF( .NOT.TPSD ) THEN 00460 * 00461 * Least-Squares Problem min || A * X - B || 00462 * 00463 * B(IB:IB+M-1,JB:JB+NRHS-1) := Q' * B(IB:IB+M-1,JB:JB+NRHS-1) 00464 * 00465 CALL PDORMQR( 'Left', 'Transpose', M, NRHS, N, A, IA, JA, 00466 $ DESCA, WORK, B, IB, JB, DESCB, WORK( IPW ), 00467 $ LWORK-LTAU, INFO ) 00468 * 00469 * workspace at least NRHS, optimally NRHS*NB 00470 * 00471 * B(IB:IB+N-1,JB:JB+NRHS-1) := inv(R) * 00472 * B(IB:IB+N-1,JB:JB+NRHS-1) 00473 * 00474 CALL PDTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', N, 00475 $ NRHS, ONE, A, IA, JA, DESCA, B, IB, JB, DESCB ) 00476 * 00477 SCLLEN = N 00478 * 00479 ELSE 00480 * 00481 * Overdetermined system of equations sub( A )' * X = sub( B ) 00482 * 00483 * sub( B ) := inv(R') * sub( B ) 00484 * 00485 CALL PDTRSM( 'Left', 'Upper', 'Transpose', 'Non-unit', N, 00486 $ NRHS, ONE, A, IA, JA, DESCA, B, IB, JB, DESCB ) 00487 * 00488 * B(IB+N:IB+M-1,JB:JB+NRHS-1) = ZERO 00489 * 00490 CALL PDLASET( 'All', M-N, NRHS, ZERO, ZERO, B, IB+N, JB, 00491 $ DESCB ) 00492 * 00493 * B(IB:IB+M-1,JB:JB+NRHS-1) := Q(1:N,:) * 00494 * B(IB:IB+N-1,JB:JB+NRHS-1) 00495 * 00496 CALL PDORMQR( 'Left', 'No transpose', M, NRHS, N, A, IA, JA, 00497 $ DESCA, WORK, B, IB, JB, DESCB, WORK( IPW ), 00498 $ LWORK-LTAU, INFO ) 00499 * 00500 * workspace at least NRHS, optimally NRHS*NB 00501 * 00502 SCLLEN = M 00503 * 00504 END IF 00505 * 00506 ELSE 00507 * 00508 * Compute LQ factorization of sub( A ) 00509 * 00510 CALL PDGELQF( M, N, A, IA, JA, DESCA, WORK, WORK( IPW ), 00511 $ LWORK-LTAU, INFO ) 00512 * 00513 * workspace at least M, optimally M*NB. 00514 * 00515 IF( .NOT.TPSD ) THEN 00516 * 00517 * underdetermined system of equations sub( A ) * X = sub( B ) 00518 * 00519 * B(IB:IB+M-1,JB:JB+NRHS-1) := inv(L) * 00520 * B(IB:IB+M-1,JB:JB+NRHS-1) 00521 * 00522 CALL PDTRSM( 'Left', 'Lower', 'No transpose', 'Non-unit', M, 00523 $ NRHS, ONE, A, IA, JA, DESCA, B, IB, JB, DESCB ) 00524 * 00525 * B(IB+M:IB+N-1,JB:JB+NRHS-1) = 0 00526 * 00527 CALL PDLASET( 'All', N-M, NRHS, ZERO, ZERO, B, IB+M, JB, 00528 $ DESCB ) 00529 * 00530 * B(IB:IB+N-1,JB:JB+NRHS-1) := Q(1:N,:)' * 00531 * B(IB:IB+M-1,JB:JB+NRHS-1) 00532 * 00533 CALL PDORMLQ( 'Left', 'Transpose', N, NRHS, M, A, IA, JA, 00534 $ DESCA, WORK, B, IB, JB, DESCB, WORK( IPW ), 00535 $ LWORK-LTAU, INFO ) 00536 * 00537 * workspace at least NRHS, optimally NRHS*NB 00538 * 00539 SCLLEN = N 00540 * 00541 ELSE 00542 * 00543 * overdetermined system min || A' * X - B || 00544 * 00545 * B(IB:IB+N-1,JB:JB+NRHS-1) := Q * B(IB:IB+N-1,JB:JB+NRHS-1) 00546 * 00547 CALL PDORMLQ( 'Left', 'No transpose', N, NRHS, M, A, IA, JA, 00548 $ DESCA, WORK, B, IB, JB, DESCB, WORK( IPW ), 00549 $ LWORK-LTAU, INFO ) 00550 * 00551 * workspace at least NRHS, optimally NRHS*NB 00552 * 00553 * B(IB:IB+M-1,JB:JB+NRHS-1) := inv(L') * 00554 * B(IB:IB+M-1,JB:JB+NRHS-1) 00555 * 00556 CALL PDTRSM( 'Left', 'Lower', 'Transpose', 'Non-unit', M, 00557 $ NRHS, ONE, A, IA, JA, DESCA, B, IB, JB, DESCB ) 00558 * 00559 SCLLEN = M 00560 * 00561 END IF 00562 * 00563 END IF 00564 * 00565 * Undo scaling 00566 * 00567 IF( IASCL.EQ.1 ) THEN 00568 CALL PDLASCL( 'G', ANRM, SMLNUM, SCLLEN, NRHS, B, IB, JB, 00569 $ DESCB, INFO ) 00570 ELSE IF( IASCL.EQ.2 ) THEN 00571 CALL PDLASCL( 'G', ANRM, BIGNUM, SCLLEN, NRHS, B, IB, JB, 00572 $ DESCB, INFO ) 00573 END IF 00574 IF( IBSCL.EQ.1 ) THEN 00575 CALL PDLASCL( 'G', SMLNUM, BNRM, SCLLEN, NRHS, B, IB, JB, 00576 $ DESCB, INFO ) 00577 ELSE IF( IBSCL.EQ.2 ) THEN 00578 CALL PDLASCL( 'G', BIGNUM, BNRM, SCLLEN, NRHS, B, IB, JB, 00579 $ DESCB, INFO ) 00580 END IF 00581 * 00582 10 CONTINUE 00583 * 00584 WORK( 1 ) = DBLE( LWMIN ) 00585 * 00586 RETURN 00587 * 00588 * End of PDGELS 00589 * 00590 END
