1 SUBROUTINE PCGELS( 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 COMPLEX A( * ), B( * ), WORK( * ) 16* .. 17* 18* Purpose 19* ======= 20* 21* PCGELS solves overdetermined or underdetermined complex linear 22* systems involving an M-by-N matrix sub( A ) = A(IA:IA+M-1,JA:JA+N-1), 23* or its conjugate-transpose, using a QR or LQ factorization of 24* sub( A ). It is 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 = 'C' and m >= n: find the minimum norm solution of 36* an undetermined system sub( A )**H * X = sub( B ). 37* 38* 4. If TRANS = 'C' 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 )**H * 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* = 'C': the linear system involves sub( A )**H. 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) COMPLEX 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 PCGEQRF; 127* if M < N, sub( A ) is overwritten by details of its LQ 128* factorization as returned by PCGELQF. 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) COMPLEX 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 = 'C' 154* and M >= N, rows 1 to M of sub( B ) contain the minimum norm 155* solution vectors; if TRANS = 'C' 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) COMPLEX 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 REAL ZERO, ONE 235 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 236 COMPLEX CZERO, CONE 237 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 238 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 239* .. 240* .. Local Scalars .. 241 LOGICAL LQUERY, TPSD 242 INTEGER BROW, IACOL, IAROW, IASCL, IBCOL, IBROW, IBSCL, 243 $ ICOFFA, ICOFFB, ICTXT, IPW, IROFFA, IROFFB, 244 $ LCM, LCMP, LTAU, LWF, LWMIN, LWS, MPA0, MPB0, 245 $ MYCOL, MYROW, NPB0, NPCOL, NPROW, NQA0, 246 $ NRHSQB0, SCLLEN 247 REAL ANRM, BIGNUM, BNRM, SMLNUM 248* .. 249* .. Local Arrays .. 250 INTEGER IDUM1( 2 ), IDUM2( 2 ) 251 REAL RWORK( 1 ) 252* .. 253* .. External Functions .. 254 LOGICAL LSAME 255 INTEGER ILCM 256 INTEGER INDXG2P, NUMROC 257 REAL PCLANGE, PSLAMCH 258 EXTERNAL ILCM, INDXG2P, LSAME, NUMROC, PCLANGE, 259 $ PSLAMCH 260* .. 261* .. External Subroutines .. 262 EXTERNAL BLACS_GRIDINFO, CHK1MAT, PCHK2MAT, PCGELQF, 263 $ PCGEQRF, PSLABAD, PCLASCL, PCLASET, 264 $ PCTRSM, PCUNMLQ, PCUNMQR, PXERBLA 265* .. 266* .. Intrinsic Functions .. 267 INTRINSIC CMPLX, ICHAR, MAX, MIN, MOD, REAL 268* .. 269* .. Executable Statements .. 270* 271* Get grid parameters 272* 273 ICTXT = DESCA( CTXT_ ) 274 CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL ) 275* 276* Test the input parameters 277* 278 INFO = 0 279 IF( NPROW.EQ.-1 ) THEN 280 INFO = -( 800 + CTXT_ ) 281 ELSE 282 CALL CHK1MAT( M, 2, N, 3, IA, JA, DESCA, 8, INFO ) 283 IF ( M .GE. N ) THEN 284 CALL CHK1MAT( M, 2, NRHS, 4, IB, JB, DESCB, 12, INFO ) 285 ELSE 286 CALL CHK1MAT( N, 3, NRHS, 4, IB, JB, DESCB, 12, INFO ) 287 ENDIF 288 IF( INFO.EQ.0 ) THEN 289 IROFFA = MOD( IA-1, DESCA( MB_ ) ) 290 ICOFFA = MOD( JA-1, DESCA( NB_ ) ) 291 IROFFB = MOD( IB-1, DESCB( MB_ ) ) 292 ICOFFB = MOD( JB-1, DESCB( NB_ ) ) 293 IAROW = INDXG2P( IA, DESCA( MB_ ), MYROW, DESCA( RSRC_ ), 294 $ NPROW ) 295 IACOL = INDXG2P( IA, DESCA( NB_ ), MYCOL, DESCA( CSRC_ ), 296 $ NPCOL ) 297 MPA0 = NUMROC( M+IROFFA, DESCA( MB_ ), MYROW, IAROW, NPROW ) 298 NQA0 = NUMROC( N+ICOFFA, DESCA( NB_ ), MYCOL, IACOL, NPCOL ) 299* 300 IBROW = INDXG2P( IB, DESCB( MB_ ), MYROW, DESCB( RSRC_ ), 301 $ NPROW ) 302 IBCOL = INDXG2P( IB, DESCB( NB_ ), MYCOL, DESCB( CSRC_ ), 303 $ NPCOL ) 304 NRHSQB0 = NUMROC( NRHS+ICOFFB, DESCB( NB_ ), MYCOL, IBCOL, 305 $ NPCOL ) 306 IF( M.GE.N ) THEN 307 MPB0 = NUMROC( M+IROFFB, DESCB( MB_ ), MYROW, IBROW, 308 $ NPROW ) 309 LTAU = NUMROC( JA+MIN(M,N)-1, DESCA( NB_ ), MYCOL, 310 $ DESCA( CSRC_ ), NPCOL ) 311 LWF = DESCA( NB_ ) * ( MPA0 + NQA0 + DESCA( NB_ ) ) 312 LWS = MAX( ( DESCA( NB_ )*( DESCA( NB_ ) - 1 ) ) / 2, 313 $ ( MPB0 + NRHSQB0 ) * DESCA( NB_ ) ) + 314 $ DESCA( NB_ )*DESCA( NB_ ) 315 ELSE 316 LCM = ILCM( NPROW, NPCOL ) 317 LCMP = LCM / NPROW 318 NPB0 = NUMROC( N+IROFFB, DESCB( MB_ ), MYROW, IBROW, 319 $ NPROW ) 320 LTAU = NUMROC( IA+MIN(M,N)-1, DESCA( MB_ ), MYROW, 321 $ DESCA( RSRC_ ), NPROW ) 322 LWF = DESCA( MB_ ) * ( MPA0 + NQA0 + DESCA( MB_ ) ) 323 LWS = MAX( ( DESCA( MB_ )*( DESCA( MB_ ) - 1 ) ) / 2, 324 $ ( NPB0 + MAX( NQA0 + NUMROC( NUMROC( N+IROFFB, 325 $ DESCA( MB_ ), 0, 0, NPROW ), DESCA( MB_ ), 0, 0, 326 $ LCMP ), NRHSQB0 ) )*DESCA( MB_ ) ) + 327 $ DESCA( MB_ ) * DESCA( MB_ ) 328 END IF 329 LWMIN = LTAU + MAX( LWF, LWS ) 330 WORK( 1 ) = CMPLX( REAL( LWMIN ) ) 331 LQUERY = ( LWORK.EQ.-1 ) 332* 333 TPSD = .TRUE. 334 IF( LSAME( TRANS, 'N' ) ) 335 $ TPSD = .FALSE. 336* 337 IF( .NOT.( LSAME( TRANS, 'N' ) .OR. 338 $ LSAME( TRANS, 'C' ) ) ) THEN 339 INFO = -1 340 ELSE IF( M.LT.0 ) THEN 341 INFO = -2 342 ELSE IF( N.LT.0 ) THEN 343 INFO = -3 344 ELSE IF( NRHS.LT.0 ) THEN 345 INFO = -4 346 ELSE IF( M.GE.N .AND. IROFFA.NE.IROFFB ) THEN 347 INFO = -10 348 ELSE IF( M.GE.N .AND. IAROW.NE.IBROW ) THEN 349 INFO = -10 350 ELSE IF( M.LT.N .AND. ICOFFA.NE.IROFFB ) THEN 351 INFO = -10 352 ELSE IF( M.GE.N .AND. DESCA( MB_ ).NE.DESCB( MB_ ) ) THEN 353 INFO = -( 1200 + MB_ ) 354 ELSE IF( M.LT.N .AND. DESCA( NB_ ).NE.DESCB( MB_ ) ) THEN 355 INFO = -( 1200 + MB_ ) 356 ELSE IF( ICTXT.NE.DESCB( CTXT_ ) ) THEN 357 INFO = -( 1200 + CTXT_ ) 358 ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 359 INFO = -14 360 END IF 361 END IF 362* 363 IF( .NOT.TPSD ) THEN 364 IDUM1( 1 ) = ICHAR( 'N' ) 365 ELSE 366 IDUM1( 1 ) = ICHAR( 'C' ) 367 END IF 368 IDUM2( 1 ) = 1 369 IF( LWORK.EQ.-1 ) THEN 370 IDUM1( 2 ) = -1 371 ELSE 372 IDUM1( 2 ) = 1 373 END IF 374 IDUM2( 2 ) = 14 375 CALL PCHK2MAT( M, 2, N, 3, IA, JA, DESCA, 8, N, 3, NRHS, 4, 376 $ IB, JB, DESCB, 12, 2, IDUM1, IDUM2, INFO ) 377 END IF 378* 379 IF( INFO.NE.0 ) THEN 380 CALL PXERBLA( ICTXT, 'PCGELS', -INFO ) 381 RETURN 382 ELSE IF( LQUERY ) THEN 383 RETURN 384 END IF 385* 386* Quick return if possible 387* 388 IF( MIN( M, N, NRHS ).EQ.0 ) THEN 389 CALL PCLASET( 'Full', MAX( M, N ), NRHS, CZERO, CZERO, B, 390 $ IB, JB, DESCB ) 391 RETURN 392 END IF 393* 394* Get machine parameters 395* 396 SMLNUM = PSLAMCH( ICTXT, 'S' ) 397 SMLNUM = SMLNUM / PSLAMCH( ICTXT, 'P' ) 398 BIGNUM = ONE / SMLNUM 399 CALL PSLABAD( ICTXT, SMLNUM, BIGNUM ) 400* 401* Scale A, B if max entry outside range [SMLNUM,BIGNUM] 402* 403 ANRM = PCLANGE( 'M', M, N, A, IA, JA, DESCA, RWORK ) 404 IASCL = 0 405 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 406* 407* Scale matrix norm up to SMLNUM 408* 409 CALL PCLASCL( 'G', ANRM, SMLNUM, M, N, A, IA, JA, DESCA, 410 $ INFO ) 411 IASCL = 1 412 ELSE IF( ANRM.GT.BIGNUM ) THEN 413* 414* Scale matrix norm down to BIGNUM 415* 416 CALL PCLASCL( 'G', ANRM, BIGNUM, M, N, A, IA, JA, DESCA, 417 $ INFO ) 418 IASCL = 2 419 ELSE IF( ANRM.EQ.ZERO ) THEN 420* 421* Matrix all zero. Return zero solution. 422* 423 CALL PCLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, IB, 424 $ JB, DESCB ) 425 GO TO 10 426 END IF 427* 428 BROW = M 429 IF( TPSD ) 430 $ BROW = N 431* 432 BNRM = PCLANGE( 'M', BROW, NRHS, B, IB, JB, DESCB, RWORK ) 433* 434 IBSCL = 0 435 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 436* 437* Scale matrix norm up to SMLNUM 438* 439 CALL PCLASCL( 'G', BNRM, SMLNUM, BROW, NRHS, B, IB, JB, 440 $ DESCB, INFO ) 441 IBSCL = 1 442 ELSE IF( BNRM.GT.BIGNUM ) THEN 443* 444* Scale matrix norm down to BIGNUM 445* 446 CALL PCLASCL( 'G', BNRM, BIGNUM, BROW, NRHS, B, IB, JB, 447 $ DESCB, INFO ) 448 IBSCL = 2 449 END IF 450* 451 IPW = LTAU + 1 452* 453 IF( M.GE.N ) THEN 454* 455* compute QR factorization of A 456* 457 CALL PCGEQRF( M, N, A, IA, JA, DESCA, WORK, WORK( IPW ), 458 $ LWORK-LTAU, INFO ) 459* 460* workspace at least N, optimally N*NB 461* 462 IF( .NOT.TPSD ) THEN 463* 464* Least-Squares Problem min || A * X - B || 465* 466* B(IB:IB+M-1,JB:JB+NRHS-1) := Q' * B(IB:IB+M-1,JB:JB+NRHS-1) 467* 468 CALL PCUNMQR( 'Left', 'Conjugate transpose', M, NRHS, N, A, 469 $ IA, JA, DESCA, WORK, B, IB, JB, DESCB, 470 $ WORK( IPW ), LWORK-LTAU, INFO ) 471* 472* workspace at least NRHS, optimally NRHS*NB 473* 474* B(IB:IB+N-1,JB:JB+NRHS-1) := inv(R) * 475* B(IB:IB+N-1,JB:JB+NRHS-1) 476* 477 CALL PCTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', N, 478 $ NRHS, CONE, A, IA, JA, DESCA, B, IB, JB, 479 $ DESCB ) 480* 481 SCLLEN = N 482* 483 ELSE 484* 485* Overdetermined system of equations sub( A )' * X = sub( B ) 486* 487* sub( B ) := inv(R') * sub( B ) 488* 489 CALL PCTRSM( 'Left', 'Upper', 'Conjugate transpose', 490 $ 'Non-unit', N, NRHS, CONE, A, IA, JA, DESCA, 491 $ B, IB, JB, DESCB ) 492* 493* B(IB+N:IB+M-1,JB:JB+NRHS-1) = ZERO 494* 495 CALL PCLASET( 'All', M-N, NRHS, CZERO, CZERO, B, IB+N, JB, 496 $ DESCB ) 497* 498* B(IB:IB+M-1,JB:JB+NRHS-1) := Q(1:N,:) * 499* B(IB:IB+N-1,JB:JB+NRHS-1) 500* 501 CALL PCUNMQR( 'Left', 'No transpose', M, NRHS, N, A, IA, JA, 502 $ DESCA, WORK, B, IB, JB, DESCB, WORK( IPW ), 503 $ LWORK-LTAU, INFO ) 504* 505* workspace at least NRHS, optimally NRHS*NB 506* 507 SCLLEN = M 508* 509 END IF 510* 511 ELSE 512* 513* Compute LQ factorization of sub( A ) 514* 515 CALL PCGELQF( M, N, A, IA, JA, DESCA, WORK, WORK( IPW ), 516 $ LWORK-LTAU, INFO ) 517* 518* workspace at least M, optimally M*NB. 519* 520 IF( .NOT.TPSD ) THEN 521* 522* underdetermined system of equations sub( A ) * X = sub( B ) 523* 524* B(IB:IB+M-1,JB:JB+NRHS-1) := inv(L) * 525* B(IB:IB+M-1,JB:JB+NRHS-1) 526* 527 CALL PCTRSM( 'Left', 'Lower', 'No transpose', 'Non-unit', M, 528 $ NRHS, CONE, A, IA, JA, DESCA, B, IB, JB, 529 $ DESCB ) 530* 531* B(IB+M:IB+N-1,JB:JB+NRHS-1) = 0 532* 533 CALL PCLASET( 'All', N-M, NRHS, CZERO, CZERO, B, IB+M, JB, 534 $ DESCB ) 535* 536* B(IB:IB+N-1,JB:JB+NRHS-1) := Q(1:N,:)' * 537* B(IB:IB+M-1,JB:JB+NRHS-1) 538* 539 CALL PCUNMLQ( 'Left', 'Conjugate transpose', N, NRHS, M, A, 540 $ IA, JA, DESCA, WORK, B, IB, JB, DESCB, 541 $ WORK( IPW ), LWORK-LTAU, INFO ) 542* 543* workspace at least NRHS, optimally NRHS*NB 544* 545 SCLLEN = N 546* 547 ELSE 548* 549* overdetermined system min || A' * X - B || 550* 551* B(IB:IB+N-1,JB:JB+NRHS-1) := Q * B(IB:IB+N-1,JB:JB+NRHS-1) 552* 553 CALL PCUNMLQ( 'Left', 'No transpose', N, NRHS, M, A, IA, JA, 554 $ DESCA, WORK, B, IB, JB, DESCB, WORK( IPW ), 555 $ LWORK-LTAU, INFO ) 556* 557* workspace at least NRHS, optimally NRHS*NB 558* 559* B(IB:IB+M-1,JB:JB+NRHS-1) := inv(L') * 560* B(IB:IB+M-1,JB:JB+NRHS-1) 561* 562 CALL PCTRSM( 'Left', 'Lower', 'Conjugate transpose', 563 $ 'Non-unit', M, NRHS, CONE, A, IA, JA, DESCA, 564 $ B, IB, JB, DESCB ) 565* 566 SCLLEN = M 567* 568 END IF 569* 570 END IF 571* 572* Undo scaling 573* 574 IF( IASCL.EQ.1 ) THEN 575 CALL PCLASCL( 'G', ANRM, SMLNUM, SCLLEN, NRHS, B, IB, JB, 576 $ DESCB, INFO ) 577 ELSE IF( IASCL.EQ.2 ) THEN 578 CALL PCLASCL( 'G', ANRM, BIGNUM, SCLLEN, NRHS, B, IB, JB, 579 $ DESCB, INFO ) 580 END IF 581 IF( IBSCL.EQ.1 ) THEN 582 CALL PCLASCL( 'G', SMLNUM, BNRM, SCLLEN, NRHS, B, IB, JB, 583 $ DESCB, INFO ) 584 ELSE IF( IBSCL.EQ.2 ) THEN 585 CALL PCLASCL( 'G', BIGNUM, BNRM, SCLLEN, NRHS, B, IB, JB, 586 $ DESCB, INFO ) 587 END IF 588* 589 10 CONTINUE 590* 591 WORK( 1 ) = CMPLX( REAL( LWMIN ) ) 592* 593 RETURN 594* 595* End of PCGELS 596* 597 END 598