1 SUBROUTINE PDGELS( TRANS, M, N, NRHS, A, IA, JA, DESCA, B, IB, JB, 2 $ DESCB, 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 CHARACTER TRANS 11 INTEGER IA, IB, INFO, JA, JB, LWORK, M, N, NRHS 12* .. 13* .. Array Arguments .. 14 INTEGER DESCA( * ), DESCB( * ) 15 DOUBLE PRECISION A( * ), B( * ), WORK( * ) 16* .. 17* 18* Purpose 19* ======= 20* 21* PDGELS solves overdetermined or underdetermined real linear 22* systems involving an M-by-N matrix sub( A ) = A(IA:IA+M-1,JA:JA+N-1), 23* or its transpose, using a QR or LQ factorization of sub( A ). It is 24* assumed that sub( A ) has full rank. 25* 26* The following options are provided: 27* 28* 1. If TRANS = 'N' and m >= n: find the least squares solution of 29* an overdetermined system, i.e., solve the least squares problem 30* minimize || sub( B ) - sub( A )*X ||. 31* 32* 2. If TRANS = 'N' and m < n: find the minimum norm solution of 33* an underdetermined system sub( A ) * X = sub( B ). 34* 35* 3. If TRANS = 'T' and m >= n: find the minimum norm solution of 36* an undetermined system sub( A )**T * X = sub( B ). 37* 38* 4. If TRANS = 'T' and m < n: find the least squares solution of 39* an overdetermined system, i.e., solve the least squares problem 40* minimize || sub( B ) - sub( A )**T * X ||. 41* 42* where sub( B ) denotes B( IB:IB+M-1, JB:JB+NRHS-1 ) when TRANS = 'N' 43* and B( IB:IB+N-1, JB:JB+NRHS-1 ) otherwise. Several right hand side 44* vectors b and solution vectors x can be handled in a single call; 45* When TRANS = 'N', the solution vectors are stored as the columns of 46* the N-by-NRHS right hand side matrix sub( B ) and the M-by-NRHS 47* right hand side matrix sub( B ) otherwise. 48* 49* Notes 50* ===== 51* 52* Each global data object is described by an associated description 53* vector. This vector stores the information required to establish 54* the mapping between an object element and its corresponding process 55* and memory location. 56* 57* Let A be a generic term for any 2D block cyclicly distributed array. 58* Such a global array has an associated description vector DESCA. 59* In the following comments, the character _ should be read as 60* "of the global array". 61* 62* NOTATION STORED IN EXPLANATION 63* --------------- -------------- -------------------------------------- 64* DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case, 65* DTYPE_A = 1. 66* CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating 67* the BLACS process grid A is distribu- 68* ted over. The context itself is glo- 69* bal, but the handle (the integer 70* value) may vary. 71* M_A (global) DESCA( M_ ) The number of rows in the global 72* array A. 73* N_A (global) DESCA( N_ ) The number of columns in the global 74* array A. 75* MB_A (global) DESCA( MB_ ) The blocking factor used to distribute 76* the rows of the array. 77* NB_A (global) DESCA( NB_ ) The blocking factor used to distribute 78* the columns of the array. 79* RSRC_A (global) DESCA( RSRC_ ) The process row over which the first 80* row of the array A is distributed. 81* CSRC_A (global) DESCA( CSRC_ ) The process column over which the 82* first column of the array A is 83* distributed. 84* LLD_A (local) DESCA( LLD_ ) The leading dimension of the local 85* array. LLD_A >= MAX(1,LOCr(M_A)). 86* 87* Let K be the number of rows or columns of a distributed matrix, 88* and assume that its process grid has dimension p x q. 89* LOCr( K ) denotes the number of elements of K that a process 90* would receive if K were distributed over the p processes of its 91* process column. 92* Similarly, LOCc( K ) denotes the number of elements of K that a 93* process would receive if K were distributed over the q processes of 94* its process row. 95* The values of LOCr() and LOCc() may be determined via a call to the 96* ScaLAPACK tool function, NUMROC: 97* LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ), 98* LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ). 99* An upper bound for these quantities may be computed by: 100* LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A 101* LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A 102* 103* Arguments 104* ========= 105* 106* TRANS (global input) CHARACTER 107* = 'N': the linear system involves sub( A ); 108* = 'T': the linear system involves sub( A )**T. 109* 110* M (global input) INTEGER 111* The number of rows to be operated on, i.e. the number of 112* rows of the distributed submatrix sub( A ). M >= 0. 113* 114* N (global input) INTEGER 115* The number of columns to be operated on, i.e. the number of 116* columns of the distributed submatrix sub( A ). N >= 0. 117* 118* NRHS (global input) INTEGER 119* The number of right hand sides, i.e. the number of columns 120* of the distributed submatrices sub( B ) and X. NRHS >= 0. 121* 122* A (local input/local output) DOUBLE PRECISION pointer into the 123* local memory to an array of local dimension 124* ( LLD_A, LOCc(JA+N-1) ). On entry, the M-by-N matrix A. 125* if M >= N, sub( A ) is overwritten by details of its QR 126* factorization as returned by PDGEQRF; 127* if M < N, sub( A ) is overwritten by details of its LQ 128* factorization as returned by PDGELQF. 129* 130* IA (global input) INTEGER 131* The row index in the global array A indicating the first 132* row of sub( A ). 133* 134* JA (global input) INTEGER 135* The column index in the global array A indicating the 136* first column of sub( A ). 137* 138* DESCA (global and local input) INTEGER array of dimension DLEN_. 139* The array descriptor for the distributed matrix A. 140* 141* B (local input/local output) DOUBLE PRECISION pointer into the 142* local memory to an array of local dimension 143* (LLD_B, LOCc(JB+NRHS-1)). On entry, this array contains the 144* local pieces of the distributed matrix B of right hand side 145* vectors, stored columnwise; 146* sub( B ) is M-by-NRHS if TRANS='N', and N-by-NRHS otherwise. 147* On exit, sub( B ) is overwritten by the solution vectors, 148* stored columnwise: if TRANS = 'N' and M >= N, rows 1 to N 149* of sub( B ) contain the least squares solution vectors; the 150* residual sum of squares for the solution in each column is 151* given by the sum of squares of elements N+1 to M in that 152* column; if TRANS = 'N' and M < N, rows 1 to N of sub( B ) 153* contain the minimum norm solution vectors; if TRANS = 'T' 154* and M >= N, rows 1 to M of sub( B ) contain the minimum norm 155* solution vectors; if TRANS = 'T' and M < N, rows 1 to M of 156* sub( B ) contain the least squares solution vectors; the 157* residual sum of squares for the solution in each column is 158* given by the sum of squares of elements M+1 to N in that 159* column. 160* 161* IB (global input) INTEGER 162* The row index in the global array B indicating the first 163* row of sub( B ). 164* 165* JB (global input) INTEGER 166* The column index in the global array B indicating the 167* first column of sub( B ). 168* 169* DESCB (global and local input) INTEGER array of dimension DLEN_. 170* The array descriptor for the distributed matrix B. 171* 172* WORK (local workspace/local output) DOUBLE PRECISION array, 173* dimension (LWORK) 174* On exit, WORK(1) returns the minimal and optimal LWORK. 175* 176* LWORK (local or global input) INTEGER 177* The dimension of the array WORK. 178* LWORK is local input and must be at least 179* LWORK >= LTAU + MAX( LWF, LWS ) where 180* If M >= N, then 181* LTAU = NUMROC( JA+MIN(M,N)-1, NB_A, MYCOL, CSRC_A, NPCOL ), 182* LWF = NB_A * ( MpA0 + NqA0 + NB_A ) 183* LWS = MAX( (NB_A*(NB_A-1))/2, (NRHSqB0 + MpB0)*NB_A ) + 184* NB_A * NB_A 185* Else 186* LTAU = NUMROC( IA+MIN(M,N)-1, MB_A, MYROW, RSRC_A, NPROW ), 187* LWF = MB_A * ( MpA0 + NqA0 + MB_A ) 188* LWS = MAX( (MB_A*(MB_A-1))/2, ( NpB0 + MAX( NqA0 + 189* NUMROC( NUMROC( N+IROFFB, MB_A, 0, 0, NPROW ), 190* MB_A, 0, 0, LCMP ), NRHSqB0 ) )*MB_A ) + 191* MB_A * MB_A 192* End if 193* 194* where LCMP = LCM / NPROW with LCM = ILCM( NPROW, NPCOL ), 195* 196* IROFFA = MOD( IA-1, MB_A ), ICOFFA = MOD( JA-1, NB_A ), 197* IAROW = INDXG2P( IA, MB_A, MYROW, RSRC_A, NPROW ), 198* IACOL = INDXG2P( JA, NB_A, MYCOL, CSRC_A, NPCOL ), 199* MpA0 = NUMROC( M+IROFFA, MB_A, MYROW, IAROW, NPROW ), 200* NqA0 = NUMROC( N+ICOFFA, NB_A, MYCOL, IACOL, NPCOL ), 201* 202* IROFFB = MOD( IB-1, MB_B ), ICOFFB = MOD( JB-1, NB_B ), 203* IBROW = INDXG2P( IB, MB_B, MYROW, RSRC_B, NPROW ), 204* IBCOL = INDXG2P( JB, NB_B, MYCOL, CSRC_B, NPCOL ), 205* MpB0 = NUMROC( M+IROFFB, MB_B, MYROW, IBROW, NPROW ), 206* NpB0 = NUMROC( N+IROFFB, MB_B, MYROW, IBROW, NPROW ), 207* NRHSqB0 = NUMROC( NRHS+ICOFFB, NB_B, MYCOL, IBCOL, NPCOL ), 208* 209* ILCM, INDXG2P and NUMROC are ScaLAPACK tool functions; 210* MYROW, MYCOL, NPROW and NPCOL can be determined by calling 211* the subroutine BLACS_GRIDINFO. 212* 213* If LWORK = -1, then LWORK is global input and a workspace 214* query is assumed; the routine only calculates the minimum 215* and optimal size for all work arrays. Each of these 216* values is returned in the first entry of the corresponding 217* work array, and no error message is issued by PXERBLA. 218* 219* INFO (global output) INTEGER 220* = 0: successful exit 221* < 0: If the i-th argument is an array and the j-entry had 222* an illegal value, then INFO = -(i*100+j), if the i-th 223* argument is a scalar and had an illegal value, then 224* INFO = -i. 225* 226* ===================================================================== 227* 228* .. Parameters .. 229 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_, 230 $ LLD_, MB_, M_, NB_, N_, RSRC_ 231 PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1, 232 $ CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, 233 $ RSRC_ = 7, CSRC_ = 8, LLD_ = 9 ) 234 DOUBLE PRECISION ZERO, ONE 235 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 236* .. 237* .. Local Scalars .. 238 LOGICAL LQUERY, TPSD 239 INTEGER BROW, IACOL, IAROW, IASCL, IBCOL, IBROW, IBSCL, 240 $ ICOFFA, ICOFFB, ICTXT, IPW, IROFFA, IROFFB, 241 $ LCM, LCMP, LTAU, LWF, LWMIN, LWS, MPA0, MPB0, 242 $ MYCOL, MYROW, NPB0, NPCOL, NPROW, NQA0, 243 $ NRHSQB0, SCLLEN 244 DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM 245* .. 246* .. Local Arrays .. 247 INTEGER IDUM1( 2 ), IDUM2( 2 ) 248 DOUBLE PRECISION RWORK( 1 ) 249* .. 250* .. External Functions .. 251 LOGICAL LSAME 252 INTEGER ILCM 253 INTEGER INDXG2P, NUMROC 254 DOUBLE PRECISION PDLAMCH, PDLANGE 255 EXTERNAL ILCM, INDXG2P, LSAME, NUMROC, PDLAMCH, 256 $ PDLANGE 257* .. 258* .. External Subroutines .. 259 EXTERNAL BLACS_GRIDINFO, CHK1MAT, PCHK2MAT, PDGELQF, 260 $ PDGEQRF, PDLABAD, PDLASCL, PDLASET, 261 $ PDORMLQ, PDORMQR, PDTRSM, PXERBLA 262* .. 263* .. Intrinsic Functions .. 264 INTRINSIC DBLE, ICHAR, MAX, MIN, MOD 265* .. 266* .. Executable Statements .. 267* 268* Get grid parameters 269* 270 ICTXT = DESCA( CTXT_ ) 271 CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL ) 272* 273* Test the input parameters 274* 275 INFO = 0 276 IF( NPROW.EQ.-1 ) THEN 277 INFO = -( 800 + CTXT_ ) 278 ELSE 279 CALL CHK1MAT( M, 2, N, 3, IA, JA, DESCA, 8, INFO ) 280 IF ( M .GE. N ) THEN 281 CALL CHK1MAT( M, 2, NRHS, 4, IB, JB, DESCB, 12, INFO ) 282 ELSE 283 CALL CHK1MAT( N, 3, NRHS, 4, IB, JB, DESCB, 12, INFO ) 284 ENDIF 285 IF( INFO.EQ.0 ) THEN 286 IROFFA = MOD( IA-1, DESCA( MB_ ) ) 287 ICOFFA = MOD( JA-1, DESCA( NB_ ) ) 288 IROFFB = MOD( IB-1, DESCB( MB_ ) ) 289 ICOFFB = MOD( JB-1, DESCB( NB_ ) ) 290 IAROW = INDXG2P( IA, DESCA( MB_ ), MYROW, DESCA( RSRC_ ), 291 $ NPROW ) 292 IACOL = INDXG2P( IA, DESCA( NB_ ), MYCOL, DESCA( CSRC_ ), 293 $ NPCOL ) 294 MPA0 = NUMROC( M+IROFFA, DESCA( MB_ ), MYROW, IAROW, NPROW ) 295 NQA0 = NUMROC( N+ICOFFA, DESCA( NB_ ), MYCOL, IACOL, NPCOL ) 296* 297 IBROW = INDXG2P( IB, DESCB( MB_ ), MYROW, DESCB( RSRC_ ), 298 $ NPROW ) 299 IBCOL = INDXG2P( IB, DESCB( NB_ ), MYCOL, DESCB( CSRC_ ), 300 $ NPCOL ) 301 NRHSQB0 = NUMROC( NRHS+ICOFFB, DESCB( NB_ ), MYCOL, IBCOL, 302 $ NPCOL ) 303 IF( M.GE.N ) THEN 304 MPB0 = NUMROC( M+IROFFB, DESCB( MB_ ), MYROW, IBROW, 305 $ NPROW ) 306 LTAU = NUMROC( JA+MIN(M,N)-1, DESCA( NB_ ), MYCOL, 307 $ DESCA( CSRC_ ), NPCOL ) 308 LWF = DESCA( NB_ ) * ( MPA0 + NQA0 + DESCA( NB_ ) ) 309 LWS = MAX( ( DESCA( NB_ )*( DESCA( NB_ ) - 1 ) ) / 2, 310 $ ( MPB0 + NRHSQB0 ) * DESCA( NB_ ) ) + 311 $ DESCA( NB_ )*DESCA( NB_ ) 312 ELSE 313 LCM = ILCM( NPROW, NPCOL ) 314 LCMP = LCM / NPROW 315 NPB0 = NUMROC( N+IROFFB, DESCB( MB_ ), MYROW, IBROW, 316 $ NPROW ) 317 LTAU = NUMROC( IA+MIN(M,N)-1, DESCA( MB_ ), MYROW, 318 $ DESCA( RSRC_ ), NPROW ) 319 LWF = DESCA( MB_ ) * ( MPA0 + NQA0 + DESCA( MB_ ) ) 320 LWS = MAX( ( DESCA( MB_ )*( DESCA( MB_ ) - 1 ) ) / 2, 321 $ ( NPB0 + MAX( NQA0 + NUMROC( NUMROC( N+IROFFB, 322 $ DESCA( MB_ ), 0, 0, NPROW ), DESCA( MB_ ), 0, 0, 323 $ LCMP ), NRHSQB0 ) )*DESCA( MB_ ) ) + 324 $ DESCA( MB_ ) * DESCA( MB_ ) 325 END IF 326 LWMIN = LTAU + MAX( LWF, LWS ) 327 WORK( 1 ) = DBLE( LWMIN ) 328 LQUERY = ( LWORK.EQ.-1 ) 329* 330 TPSD = .TRUE. 331 IF( LSAME( TRANS, 'N' ) ) 332 $ TPSD = .FALSE. 333* 334 IF( .NOT.( LSAME( TRANS, 'N' ) .OR. 335 $ LSAME( TRANS, 'T' ) ) ) THEN 336 INFO = -1 337 ELSE IF( M.LT.0 ) THEN 338 INFO = -2 339 ELSE IF( N.LT.0 ) THEN 340 INFO = -3 341 ELSE IF( NRHS.LT.0 ) THEN 342 INFO = -4 343 ELSE IF( M.GE.N .AND. IROFFA.NE.IROFFB ) THEN 344 INFO = -10 345 ELSE IF( M.GE.N .AND. IAROW.NE.IBROW ) THEN 346 INFO = -10 347 ELSE IF( M.LT.N .AND. ICOFFA.NE.IROFFB ) THEN 348 INFO = -10 349 ELSE IF( M.GE.N .AND. DESCA( MB_ ).NE.DESCB( MB_ ) ) THEN 350 INFO = -( 1200 + MB_ ) 351 ELSE IF( M.LT.N .AND. DESCA( NB_ ).NE.DESCB( MB_ ) ) THEN 352 INFO = -( 1200 + MB_ ) 353 ELSE IF( ICTXT.NE.DESCB( CTXT_ ) ) THEN 354 INFO = -( 1200 + CTXT_ ) 355 ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 356 INFO = -14 357 END IF 358 END IF 359* 360 IF( .NOT.TPSD ) THEN 361 IDUM1( 1 ) = ICHAR( 'N' ) 362 ELSE 363 IDUM1( 1 ) = ICHAR( 'T' ) 364 END IF 365 IDUM2( 1 ) = 1 366 IF( LWORK.EQ.-1 ) THEN 367 IDUM1( 2 ) = -1 368 ELSE 369 IDUM1( 2 ) = 1 370 END IF 371 IDUM2( 2 ) = 14 372 CALL PCHK2MAT( M, 2, N, 3, IA, JA, DESCA, 8, N, 3, NRHS, 4, 373 $ IB, JB, DESCB, 12, 2, IDUM1, IDUM2, INFO ) 374 END IF 375* 376 IF( INFO.NE.0 ) THEN 377 CALL PXERBLA( ICTXT, 'PDGELS', -INFO ) 378 RETURN 379 ELSE IF( LQUERY ) THEN 380 RETURN 381 END IF 382* 383* Quick return if possible 384* 385 IF( MIN( M, N, NRHS ).EQ.0 ) THEN 386 CALL PDLASET( 'Full', MAX( M, N ), NRHS, ZERO, ZERO, B, 387 $ IB, JB, DESCB ) 388 RETURN 389 END IF 390* 391* Get machine parameters 392* 393 SMLNUM = PDLAMCH( ICTXT, 'S' ) 394 SMLNUM = SMLNUM / PDLAMCH( ICTXT, 'P' ) 395 BIGNUM = ONE / SMLNUM 396 CALL PDLABAD( ICTXT, SMLNUM, BIGNUM ) 397* 398* Scale A, B if max entry outside range [SMLNUM,BIGNUM] 399* 400 ANRM = PDLANGE( 'M', M, N, A, IA, JA, DESCA, RWORK ) 401 IASCL = 0 402 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 403* 404* Scale matrix norm up to SMLNUM 405* 406 CALL PDLASCL( 'G', ANRM, SMLNUM, M, N, A, IA, JA, DESCA, 407 $ INFO ) 408 IASCL = 1 409 ELSE IF( ANRM.GT.BIGNUM ) THEN 410* 411* Scale matrix norm down to BIGNUM 412* 413 CALL PDLASCL( 'G', ANRM, BIGNUM, M, N, A, IA, JA, DESCA, 414 $ INFO ) 415 IASCL = 2 416 ELSE IF( ANRM.EQ.ZERO ) THEN 417* 418* Matrix all zero. Return zero solution. 419* 420 CALL PDLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, IB, JB, 421 $ DESCB ) 422 GO TO 10 423 END IF 424* 425 BROW = M 426 IF( TPSD ) 427 $ BROW = N 428* 429 BNRM = PDLANGE( 'M', BROW, NRHS, B, IB, JB, DESCB, RWORK ) 430* 431 IBSCL = 0 432 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 433* 434* Scale matrix norm up to SMLNUM 435* 436 CALL PDLASCL( 'G', BNRM, SMLNUM, BROW, NRHS, B, IB, JB, 437 $ DESCB, INFO ) 438 IBSCL = 1 439 ELSE IF( BNRM.GT.BIGNUM ) THEN 440* 441* Scale matrix norm down to BIGNUM 442* 443 CALL PDLASCL( 'G', BNRM, BIGNUM, BROW, NRHS, B, IB, JB, 444 $ DESCB, INFO ) 445 IBSCL = 2 446 END IF 447* 448 IPW = LTAU + 1 449* 450 IF( M.GE.N ) THEN 451* 452* compute QR factorization of A 453* 454 CALL PDGEQRF( M, N, A, IA, JA, DESCA, WORK, WORK( IPW ), 455 $ LWORK-LTAU, INFO ) 456* 457* workspace at least N, optimally N*NB 458* 459 IF( .NOT.TPSD ) THEN 460* 461* Least-Squares Problem min || A * X - B || 462* 463* B(IB:IB+M-1,JB:JB+NRHS-1) := Q' * B(IB:IB+M-1,JB:JB+NRHS-1) 464* 465 CALL PDORMQR( 'Left', 'Transpose', M, NRHS, N, A, IA, JA, 466 $ DESCA, WORK, B, IB, JB, DESCB, WORK( IPW ), 467 $ LWORK-LTAU, INFO ) 468* 469* workspace at least NRHS, optimally NRHS*NB 470* 471* B(IB:IB+N-1,JB:JB+NRHS-1) := inv(R) * 472* B(IB:IB+N-1,JB:JB+NRHS-1) 473* 474 CALL PDTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', N, 475 $ NRHS, ONE, A, IA, JA, DESCA, B, IB, JB, DESCB ) 476* 477 SCLLEN = N 478* 479 ELSE 480* 481* Overdetermined system of equations sub( A )' * X = sub( B ) 482* 483* sub( B ) := inv(R') * sub( B ) 484* 485 CALL PDTRSM( 'Left', 'Upper', 'Transpose', 'Non-unit', N, 486 $ NRHS, ONE, A, IA, JA, DESCA, B, IB, JB, DESCB ) 487* 488* B(IB+N:IB+M-1,JB:JB+NRHS-1) = ZERO 489* 490 CALL PDLASET( 'All', M-N, NRHS, ZERO, ZERO, B, IB+N, JB, 491 $ DESCB ) 492* 493* B(IB:IB+M-1,JB:JB+NRHS-1) := Q(1:N,:) * 494* B(IB:IB+N-1,JB:JB+NRHS-1) 495* 496 CALL PDORMQR( 'Left', 'No transpose', M, NRHS, N, A, IA, JA, 497 $ DESCA, WORK, B, IB, JB, DESCB, WORK( IPW ), 498 $ LWORK-LTAU, INFO ) 499* 500* workspace at least NRHS, optimally NRHS*NB 501* 502 SCLLEN = M 503* 504 END IF 505* 506 ELSE 507* 508* Compute LQ factorization of sub( A ) 509* 510 CALL PDGELQF( M, N, A, IA, JA, DESCA, WORK, WORK( IPW ), 511 $ LWORK-LTAU, INFO ) 512* 513* workspace at least M, optimally M*NB. 514* 515 IF( .NOT.TPSD ) THEN 516* 517* underdetermined system of equations sub( A ) * X = sub( B ) 518* 519* B(IB:IB+M-1,JB:JB+NRHS-1) := inv(L) * 520* B(IB:IB+M-1,JB:JB+NRHS-1) 521* 522 CALL PDTRSM( 'Left', 'Lower', 'No transpose', 'Non-unit', M, 523 $ NRHS, ONE, A, IA, JA, DESCA, B, IB, JB, DESCB ) 524* 525* B(IB+M:IB+N-1,JB:JB+NRHS-1) = 0 526* 527 CALL PDLASET( 'All', N-M, NRHS, ZERO, ZERO, B, IB+M, JB, 528 $ DESCB ) 529* 530* B(IB:IB+N-1,JB:JB+NRHS-1) := Q(1:N,:)' * 531* B(IB:IB+M-1,JB:JB+NRHS-1) 532* 533 CALL PDORMLQ( 'Left', 'Transpose', N, NRHS, M, A, IA, JA, 534 $ DESCA, WORK, B, IB, JB, DESCB, WORK( IPW ), 535 $ LWORK-LTAU, INFO ) 536* 537* workspace at least NRHS, optimally NRHS*NB 538* 539 SCLLEN = N 540* 541 ELSE 542* 543* overdetermined system min || A' * X - B || 544* 545* B(IB:IB+N-1,JB:JB+NRHS-1) := Q * B(IB:IB+N-1,JB:JB+NRHS-1) 546* 547 CALL PDORMLQ( 'Left', 'No transpose', N, NRHS, M, A, IA, JA, 548 $ DESCA, WORK, B, IB, JB, DESCB, WORK( IPW ), 549 $ LWORK-LTAU, INFO ) 550* 551* workspace at least NRHS, optimally NRHS*NB 552* 553* B(IB:IB+M-1,JB:JB+NRHS-1) := inv(L') * 554* B(IB:IB+M-1,JB:JB+NRHS-1) 555* 556 CALL PDTRSM( 'Left', 'Lower', 'Transpose', 'Non-unit', M, 557 $ NRHS, ONE, A, IA, JA, DESCA, B, IB, JB, DESCB ) 558* 559 SCLLEN = M 560* 561 END IF 562* 563 END IF 564* 565* Undo scaling 566* 567 IF( IASCL.EQ.1 ) THEN 568 CALL PDLASCL( 'G', ANRM, SMLNUM, SCLLEN, NRHS, B, IB, JB, 569 $ DESCB, INFO ) 570 ELSE IF( IASCL.EQ.2 ) THEN 571 CALL PDLASCL( 'G', ANRM, BIGNUM, SCLLEN, NRHS, B, IB, JB, 572 $ DESCB, INFO ) 573 END IF 574 IF( IBSCL.EQ.1 ) THEN 575 CALL PDLASCL( 'G', SMLNUM, BNRM, SCLLEN, NRHS, B, IB, JB, 576 $ DESCB, INFO ) 577 ELSE IF( IBSCL.EQ.2 ) THEN 578 CALL PDLASCL( 'G', BIGNUM, BNRM, SCLLEN, NRHS, B, IB, JB, 579 $ DESCB, INFO ) 580 END IF 581* 582 10 CONTINUE 583* 584 WORK( 1 ) = DBLE( LWMIN ) 585* 586 RETURN 587* 588* End of PDGELS 589* 590 END 591