Routine: PSGELS()  File: SRC\psgels.f

 
 
# lines: 586
  # code: 586
  # comment: 0
  # blank:0
# Variables:68
# Callers:0
# Callings:9
# Words:286
# Keywords:171
 

 

..
     .. Array Arguments ..
     ..
  Purpose
  =======
  PSGELS solves overdetermined or underdetermined real linear
  systems involving an M-by-N matrix sub( A ) = A(IA:IA+M-1,JA:JA+N-1),
  or its transpose, using a QR or LQ factorization of sub( A ).  It is
  assumed that sub( A ) has full rank.
  The following options are provided:
  1. If TRANS = 'N' and m >= n:  find the least squares solution of
     an overdetermined system, i.e., solve the least squares problem
                  minimize || sub( B ) - sub( A )*X ||.
  2. If TRANS = 'N' and m < n:  find the minimum norm solution of
     an underdetermined system sub( A ) * X = sub( B ).
  3. If TRANS = 'T' and m >= n:  find the minimum norm solution of
     an undetermined system sub( A )**T * X = sub( B ).
  4. If TRANS = 'T' and m < n:  find the least squares solution of
     an overdetermined system, i.e., solve the least squares problem
                  minimize || sub( B ) - sub( A )**T * X ||.
  where sub( B ) denotes B( IB:IB+M-1, JB:JB+NRHS-1 ) when TRANS = 'N'
  and B( IB:IB+N-1, JB:JB+NRHS-1 ) otherwise. Several right hand side
  vectors b and solution vectors x can be handled in a single call;
  When TRANS = 'N', the solution vectors are stored as the columns of
  the N-by-NRHS right hand side matrix sub( B ) and the M-by-NRHS
  right hand side matrix sub( B ) otherwise.
  Notes
  =====
  Each global data object is described by an associated description
  vector.  This vector stores the information required to establish
  the mapping between an object element and its corresponding process
  and memory location.
  Let A be a generic term for any 2D block cyclicly distributed array.
  Such a global array has an associated description vector DESCA.
  In the following comments, the character _ should be read as
  "of the global array".
  NOTATION        STORED IN      EXPLANATION
  --------------- -------------- --------------------------------------
  DTYPE_A(global) DESCA( DTYPE_ )The descriptor type.  In this case,
                                 DTYPE_A = 1.
  CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
                                 the BLACS process grid A is distribu-
                                 ted over. The context itself is glo-
                                 bal, but the handle (the integer
                                 value) may vary.
  M_A    (global) DESCA( M_ )    The number of rows in the global
                                 array A.
  N_A    (global) DESCA( N_ )    The number of columns in the global
                                 array A.
  MB_A   (global) DESCA( MB_ )   The blocking factor used to distribute
                                 the rows of the array.
  NB_A   (global) DESCA( NB_ )   The blocking factor used to distribute
                                 the columns of the array.
  RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
                                 row of the array A is distributed.
  CSRC_A (global) DESCA( CSRC_ ) The process column over which the
                                 first column of the array A is
                                 distributed.
  LLD_A  (local)  DESCA( LLD_ )  The leading dimension of the local
                                 array.  LLD_A >= MAX(1,LOCr(M_A)).
  Let K be the number of rows or columns of a distributed matrix,
  and assume that its process grid has dimension p x q.
  LOCr( K ) denotes the number of elements of K that a process
  would receive if K were distributed over the p processes of its
  process column.
  Similarly, LOCc( K ) denotes the number of elements of K that a
  process would receive if K were distributed over the q processes of
  its process row.
  The values of LOCr() and LOCc() may be determined via a call to the
  ScaLAPACK tool function, NUMROC:
          LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
          LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
  An upper bound for these quantities may be computed by:
          LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
          LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
  Arguments
  =========
  TRANS   (global input) CHARACTER
          = 'N': the linear system involves sub( A );
          = 'T': the linear system involves sub( A )**T.
  M       (global input) INTEGER
          The number of rows to be operated on, i.e. the number of
          rows of the distributed submatrix sub( A ). M >= 0.
  N       (global input) INTEGER
          The number of columns to be operated on, i.e. the number of
          columns of the distributed submatrix sub( A ). N >= 0.
  NRHS    (global input) INTEGER
          The number of right hand sides, i.e. the number of columns
          of the distributed submatrices sub( B ) and X.  NRHS >= 0.
  A       (local input/local output) REAL pointer into the
          local memory to an array of local dimension
          ( LLD_A, LOCc(JA+N-1) ).  On entry, the M-by-N matrix A.
          if M >= N, sub( A ) is overwritten by details of its QR
            factorization as returned by PSGEQRF;
          if M <  N, sub( A ) is overwritten by details of its LQ
            factorization as returned by PSGELQF.
  IA      (global input) INTEGER
          The row index in the global array A indicating the first
          row of sub( A ).
  JA      (global input) INTEGER
          The column index in the global array A indicating the
          first column of sub( A ).
  DESCA   (global and local input) INTEGER array of dimension DLEN_.
          The array descriptor for the distributed matrix A.
  B       (local input/local output) REAL pointer into the
          local memory to an array of local dimension
          (LLD_B, LOCc(JB+NRHS-1)).  On entry, this array contains the
          local pieces of the distributed matrix B of right hand side
          vectors, stored columnwise;
          sub( B ) is M-by-NRHS if TRANS='N', and N-by-NRHS otherwise.
          On exit, sub( B ) is overwritten by the solution vectors,
          stored columnwise:  if TRANS = 'N' and M >= N, rows 1 to N
          of sub( B ) contain the least squares solution vectors; the
          residual sum of squares for the solution in each column is
          given by the sum of squares of elements N+1 to M in that
          column; if TRANS = 'N' and M < N, rows 1 to N of sub( B )
          contain the minimum norm solution vectors; if TRANS = 'T'
          and M >= N, rows 1 to M of sub( B ) contain the minimum norm
          solution vectors; if TRANS = 'T' and M < N, rows 1 to M of
          sub( B ) contain the least squares solution vectors; the
          residual sum of squares for the solution in each column is
          given by the sum of squares of elements M+1 to N in that
          column.
  IB      (global input) INTEGER
          The row index in the global array B indicating the first
          row of sub( B ).
  JB      (global input) INTEGER
          The column index in the global array B indicating the
          first column of sub( B ).
  DESCB   (global and local input) INTEGER array of dimension DLEN_.
          The array descriptor for the distributed matrix B.
  WORK    (local workspace/local output) REAL array,
                                                  dimension (LWORK)
          On exit, WORK(1) returns the minimal and optimal LWORK.
  LWORK   (local or global input) INTEGER
          The dimension of the array WORK.
          LWORK is local input and must be at least
          LWORK >= LTAU + MAX( LWF, LWS ) where
          If M >= N, then
            LTAU = NUMROC( JA+MIN(M,N)-1, NB_A, MYCOL, CSRC_A, NPCOL ),
            LWF  = NB_A * ( MpA0 + NqA0 + NB_A )
            LWS  = MAX( (NB_A*(NB_A-1))/2, (NRHSqB0 + MpB0)*NB_A ) +
                   NB_A * NB_A
          Else
            LTAU = NUMROC( IA+MIN(M,N)-1, MB_A, MYROW, RSRC_A, NPROW ),
            LWF  = MB_A * ( MpA0 + NqA0 + MB_A )
            LWS  = MAX( (MB_A*(MB_A-1))/2, ( NpB0 + MAX( NqA0 +
                   NUMROC( NUMROC( N+IROFFB, MB_A, 0, 0, NPROW ),
                   MB_A, 0, 0, LCMP ), NRHSqB0 ) )*MB_A ) +
                   MB_A * MB_A
          End if
          where LCMP = LCM / NPROW with LCM = ILCM( NPROW, NPCOL ),
          IROFFA = MOD( IA-1, MB_A ), ICOFFA = MOD( JA-1, NB_A ),
          IAROW = INDXG2P( IA, MB_A, MYROW, RSRC_A, NPROW ),
          IACOL = INDXG2P( JA, NB_A, MYCOL, CSRC_A, NPCOL ),
          MpA0 = NUMROC( M+IROFFA, MB_A, MYROW, IAROW, NPROW ),
          NqA0 = NUMROC( N+ICOFFA, NB_A, MYCOL, IACOL, NPCOL ),
          IROFFB = MOD( IB-1, MB_B ), ICOFFB = MOD( JB-1, NB_B ),
          IBROW = INDXG2P( IB, MB_B, MYROW, RSRC_B, NPROW ),
          IBCOL = INDXG2P( JB, NB_B, MYCOL, CSRC_B, NPCOL ),
          MpB0 = NUMROC( M+IROFFB, MB_B, MYROW, IBROW, NPROW ),
          NpB0 = NUMROC( N+IROFFB, MB_B, MYROW, IBROW, NPROW ),
          NRHSqB0 = NUMROC( NRHS+ICOFFB, NB_B, MYCOL, IBCOL, NPCOL ),
          ILCM, INDXG2P and NUMROC are ScaLAPACK tool functions;
          MYROW, MYCOL, NPROW and NPCOL can be determined by calling
          the subroutine BLACS_GRIDINFO.
          If LWORK = -1, then LWORK is global input and a workspace
          query is assumed; the routine only calculates the minimum
          and optimal size for all work arrays. Each of these
          values is returned in the first entry of the corresponding
          work array, and no error message is issued by PXERBLA.
  INFO    (global output) INTEGER
          = 0:  successful exit
          < 0:  If the i-th argument is an array and the j-entry had
                an illegal value, then INFO = -(i*100+j), if the i-th
                argument is a scalar and had an illegal value, then
                INFO = -i.
  =====================================================================
     .. Parameters ..

 
Display dynamic version Find AutoScroll Reload FontSize: - + Hide Comments Hide Blanks Frame FullScreen MailPrint

 
001        SUBROUTINE PSGELS( TRANS , M , N , NRHS , A , IA , JA , DESCA , B , IB , JB ,
002       $DESCB , WORK , LWORK , INFO )
003  
004  *     -- ScaLAPACK routine(version 1.7) --
005  *     University of Tennessee , Knoxville , Oak Ridge National Laboratory ,
006  *     and University of California , Berkeley.
007  *     May 1 , 1997
008  
009  *     .. Scalar Arguments ..
010        CHARACTER TRANS
011        INTEGER IA , IB , INFO , JA , JB , LWORK , M , N , NRHS
012        INTEGER BLOCK_CYCLIC_2D , CSRC_ , CTXT_ , DLEN_ , DTYPE_ ,
013       $LLD_ , MB_ , M_ , NB_ , N_ , RSRC_
014        PARAMETER( BLOCK_CYCLIC_2D = 1 , DLEN_ = 9 , DTYPE_ = 1 ,
015       $CTXT_ = 2 , M_ = 3 , N_ = 4 , MB_ = 5 , NB_ = 6 ,
016       $RSRC_ = 7 , CSRC_ = 8 , LLD_ = 9 )
017        REAL ZERO , ONE
018        PARAMETER( ZERO = 0.0E + 0 , ONE = 1.0E + 0 )
019  *     ..
020  *     .. Local Scalars ..
021        LOGICAL LQUERY , TPSD
022        INTEGER BROW , IACOL , IAROW , IASCL , IBCOL , IBROW , IBSCL ,
023       $ICOFFA , ICOFFB , ICTXT , IPW , IROFFA , IROFFB ,
024       $LCM , LCMP , LTAU , LWF , LWMIN , LWS , MPA0 , MPB0 ,
025       $MYCOL , MYROW , NPB0 , NPCOL , NPROW , NQA0 ,
026       $NRHSQB0 , SCLLEN
027        REAL ANRM , BIGNUM , BNRM , SMLNUM
028  *     ..
029  *     .. Local Arrays ..
030        INTEGER IDUM1( 2 ) , IDUM2( 2 )
031        REAL RWORK( 1 )
032  *     ..
033  *     .. External Functions ..
034        LOGICAL LSAME
035        INTEGER ILCM
036        INTEGER INDXG2P , NUMROC
037        REAL PSLANGE , PSLAMCH
038        EXTERNAL ILCM , INDXG2P , LSAME , NUMROC , PSLANGE ,
039       $PSLAMCH
040  *     ..
041  *     .. External Subroutines ..
042        EXTERNAL BLACS_GRIDINFO , CHK1MAT , PCHK2MAT , PSGELQF ,
043       $PSGEQRF , PSLABAD , PSLASCL , PSLASET ,
044       $PSORMLQ , PSORMQR , PSTRSM , PXERBLA
045  *     ..
046  *     .. Intrinsic Functions ..
047        INTRINSIC ICHAR , MAX , MIN , MOD , REAL
048  *     ..
049  *     .. Executable Statements ..
050  
051  *     Get grid parameters
052  
053        ICTXT = DESCA( CTXT_ )
054        CALL BLACS_GRIDINFO( ICTXT , NPROW , NPCOL , MYROW , MYCOL )
055  
056  *     Test the input parameters
057  
058        INFO = 0
059        IF( NPROW.EQ. - 1 ) THEN
060            INFO = - ( 800 + CTXT_ )
061        ELSE
062            CALL CHK1MAT( M , 2 , N , 3 , IA , JA , DESCA , 8 , INFO )
063            CALL CHK1MAT( N , 3 , NRHS , 4 , IB , JB , DESCB , 12 , INFO )
064            IF( INFO.EQ.0 ) THEN
065                IROFFA = MOD( IA - 1 , DESCA( MB_ ) )
066                ICOFFA = MOD( JA - 1 , DESCA( NB_ ) )
067                IROFFB = MOD( IB - 1 , DESCB( MB_ ) )
068                ICOFFB = MOD( JB - 1 , DESCB( NB_ ) )
069                IAROW = INDXG2P( IA , DESCA( MB_ ) , MYROW , DESCA( RSRC_ ) ,
070       $        NPROW )
071                IACOL = INDXG2P( IA , DESCA( NB_ ) , MYCOL , DESCA( CSRC_ ) ,
072       $        NPCOL )
073                MPA0 = NUMROC( M + IROFFA , DESCA( MB_ ) , MYROW , IAROW , NPROW )
074                NQA0 = NUMROC( N + ICOFFA , DESCA( NB_ ) , MYCOL , IACOL , NPCOL )
075  
076                IBROW = INDXG2P( IB , DESCB( MB_ ) , MYROW , DESCB( RSRC_ ) ,
077       $        NPROW )
078                IBCOL = INDXG2P( IB , DESCB( NB_ ) , MYCOL , DESCB( CSRC_ ) ,
079       $        NPCOL )
080                NRHSQB0 = NUMROC( NRHS + ICOFFB , DESCB( NB_ ) , MYCOL , IBCOL ,
081       $        NPCOL )
082                IF( M.GE.N ) THEN
083                    MPB0 = NUMROC( M + IROFFB , DESCB( MB_ ) , MYROW , IBROW ,
084       $            NPROW )
085                    LTAU = NUMROC( JA + MIN(M , N) - 1 , DESCA( NB_ ) , MYCOL ,
086       $            DESCA( CSRC_ ) , NPCOL )
087                    LWF = DESCA( NB_ ) * ( MPA0 + NQA0 + DESCA( NB_ ) )
088                    LWS = MAX(( DESCA( NB_ )*( DESCA( NB_ ) - 1 ) ) / 2 ,
089       $( MPB0 + NRHSQB0 ) * DESCA( NB_ ) ) +
090       $            DESCA( NB_ )*DESCA( NB_ )
091                ELSE
092                    LCM = ILCM( NPROW , NPCOL )
093                    LCMP = LCM / NPROW
094                    NPB0 = NUMROC( N + IROFFB , DESCB( MB_ ) , MYROW , IBROW ,
095       $            NPROW )
096                    LTAU = NUMROC( IA + MIN(M , N) - 1 , DESCA( MB_ ) , MYROW ,
097       $            DESCA( RSRC_ ) , NPROW )
098                    LWF = DESCA( MB_ ) * ( MPA0 + NQA0 + DESCA( MB_ ) )
099                    LWS = MAX(( DESCA( MB_ )*( DESCA( MB_ ) - 1 ) ) / 2 ,
100       $( NPB0 + MAX( NQA0 + NUMROC( NUMROC( N + IROFFB ,
101       $            DESCA( MB_ ) , 0 , 0 , NPROW ) , DESCA( MB_ ) , 0 , 0 ,
102       $            LCMP ) , NRHSQB0 ) )*DESCA( MB_ ) ) +
103       $            DESCA( MB_ ) * DESCA( MB_ )
104                END IF
105                LWMIN = LTAU + MAX( LWF , LWS )
106                WORK( 1 ) = REAL( LWMIN )
107                LQUERY =( LWORK.EQ. - 1 )
108  
109                TPSD = .TRUE.
110                IF( LSAME( TRANS , 'N' ) )
111       $            TPSD = .FALSE.
112  
113                    IF( .NOT.( LSAME( TRANS , 'N' ) .OR.
114       $            LSAME( TRANS , 'T' ) ) ) THEN
115                    INFO = - 1
116                ELSE IF( M.LT.0 ) THEN
117                    INFO = - 2
118                ELSE IF( N.LT.0 ) THEN
119                    INFO = - 3
120                ELSE IF( NRHS.LT.0 ) THEN
121                    INFO = - 4
122                ELSE IF( M.GE.N .AND. IROFFA.NE.IROFFB ) THEN
123                    INFO = - 10
124                ELSE IF( M.GE.N .AND. IAROW.NE.IBROW ) THEN
125                    INFO = - 10
126                ELSE IF( M.LT.N .AND. ICOFFA.NE.IROFFB ) THEN
127                    INFO = - 10
128                ELSE IF( M.GE.N .AND. DESCA( MB_ ).NE.DESCB( MB_ ) ) THEN
129                    INFO = - ( 1200 + MB_ )
130                ELSE IF( M.LT.N .AND. DESCA( NB_ ).NE.DESCB( MB_ ) ) THEN
131                    INFO = - ( 1200 + MB_ )
132                ELSE IF( ICTXT.NE.DESCB( CTXT_ ) ) THEN
133                    INFO = - ( 1200 + CTXT_ )
134                ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
135                    INFO = - 14
136                END IF
137            END IF
138  
139            IF( .NOT.TPSD ) THEN
140                IDUM1( 1 ) = ICHAR( 'N' )
141            ELSE
142                IDUM1( 1 ) = ICHAR( 'T' )
143            END IF
144            IDUM2( 1 ) = 1
145            IF( LWORK.EQ. - 1 ) THEN
146                IDUM1( 2 ) = - 1
147            ELSE
148                IDUM1( 2 ) = 1
149            END IF
150            IDUM2( 2 ) = 14
151            CALL PCHK2MAT( M , 2 , N , 3 , IA , JA , DESCA , 8 , N , 3 , NRHS , 4 ,
152       $    IB , JB , DESCB , 12 , 2 , IDUM1 , IDUM2 , INFO )
153        END IF
154  
155        IF( INFO.NE.0 ) THEN
156            CALL PXERBLA( ICTXT , 'PSGELS' , - INFO )
157            RETURN
158        ELSE IF( LQUERY ) THEN
159            RETURN
160        END IF
161  
162  *     Quick return if possible
163  
164        IF( MIN( M , N , NRHS ).EQ.0 ) THEN
165            CALL PSLASET ( 'Full' , MAX( M , N ) , NRHS , ZERO , ZERO , B ,
166       $    IB , JB , DESCB )
167            RETURN
168        END IF
169  
170  *     Get machine parameters
171  
172        SMLNUM = PSLAMCH( ICTXT , 'S' )
173        SMLNUM = SMLNUM / PSLAMCH( ICTXT , 'P' )
174        BIGNUM = ONE / SMLNUM
175        CALL PSLABAD ( ICTXT , SMLNUM , BIGNUM )
176  
177  *     Scale A , B if max entry outside range[SMLNUM , BIGNUM]
178  
179        ANRM = PSLANGE( 'M' , M , N , A , IA , JA , DESCA , RWORK )
180        IASCL = 0
181        IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
182  
183  *         Scale matrix norm up to SMLNUM
184  
185            CALL PSLASCL ( 'G' , ANRM , SMLNUM , M , N , A , IA , JA , DESCA ,
186       $    INFO )
187            IASCL = 1
188        ELSE IF( ANRM.GT.BIGNUM ) THEN
189  
190  *         Scale matrix norm down to BIGNUM
191  
192            CALL PSLASCL ( 'G' , ANRM , BIGNUM , M , N , A , IA , JA , DESCA ,
193       $    INFO )
194            IASCL = 2
195        ELSE IF( ANRM.EQ.ZERO ) THEN
196  
197  *         Matrix all zero. Return zero solution.
198  
199            CALL PSLASET ( 'F' , MAX( M , N ) , NRHS , ZERO , ZERO , B , IB , JB ,
200       $    DESCB )
201            GO TO 10
202        END IF
203  
204        BROW = M
205        IF( TPSD )
206       $    BROW = N
207  
208            BNRM = PSLANGE( 'M' , BROW , NRHS , B , IB , JB , DESCB , RWORK )
209  
210            IBSCL = 0
211            IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
212  
213  *             Scale matrix norm up to SMLNUM
214  
215                CALL PSLASCL ( 'G' , BNRM , SMLNUM , BROW , NRHS , B , IB , JB ,
216       $        DESCB , INFO )
217                IBSCL = 1
218            ELSE IF( BNRM.GT.BIGNUM ) THEN
219  
220  *             Scale matrix norm down to BIGNUM
221  
222                CALL PSLASCL ( 'G' , BNRM , BIGNUM , BROW , NRHS , B , IB , JB ,
223       $        DESCB , INFO )
224                IBSCL = 2
225            END IF
226  
227            IPW = LTAU + 1
228  
229            IF( M.GE.N ) THEN
230  
231  *             compute QR factorization of A
232  
233                CALL PSGEQRF ( M , N , A , IA , JA , DESCA , WORK , WORK( IPW ) ,
234       $        LWORK - LTAU , INFO )
235  
236  *             workspace at least N , optimally N*NB
237  
238                IF( .NOT.TPSD ) THEN
239  
240  *                 Least - Squares Problem min || A * X - B ||
241  
242  *                 B(IB : IB + M - 1 , JB : JB + NRHS - 1) := Q' * B(IB : IB + M - 1 , JB : JB + NRHS - 1)
243  
244                    CALL PSORMQR ( 'Left' , 'Transpose' , M , NRHS , N , A , IA , JA ,
245       $            DESCA , WORK , B , IB , JB , DESCB , WORK( IPW ) ,
246       $            LWORK - LTAU , INFO )
247  
248  *                 workspace at least NRHS , optimally NRHS*NB
249  
250  *                 B(IB : IB + N - 1 , JB : JB + NRHS - 1) := inv(R) *
251  *                 B(IB : IB + N - 1 , JB : JB + NRHS - 1)
252  
253                    CALL PSTRSM( 'Left' , 'Upper' , 'No transpose' , 'Non - unit' , N ,
254       $            NRHS , ONE , A , IA , JA , DESCA , B , IB , JB , DESCB )
255  
256                    SCLLEN = N
257  
258                ELSE
259  
260  *                 Overdetermined system of equations sub( A )' * X = sub( B )
261  
262  *                 sub( B ) := inv(R') * sub( B )
263  
264                    CALL PSTRSM( 'Left' , 'Upper' , 'Transpose' , 'Non - unit' , N ,
265       $            NRHS , ONE , A , IA , JA , DESCA , B , IB , JB , DESCB )
266  
267  *                 B(IB + N : IB + M - 1 , JB : JB + NRHS - 1) = ZERO
268  
269                    CALL PSLASET ( 'All' , M - N , NRHS , ZERO , ZERO , B , IB + N , JB ,
270       $            DESCB )
271  
272  *                 B(IB : IB + M - 1 , JB : JB + NRHS - 1) := Q(1 : N ,: ) *
273  *                 B(IB : IB + N - 1 , JB : JB + NRHS - 1)
274  
275                    CALL PSORMQR ( 'Left' , 'No transpose' , M , NRHS , N , A , IA , JA ,
276       $            DESCA , WORK , B , IB , JB , DESCB , WORK( IPW ) ,
277       $            LWORK - LTAU , INFO )
278  
279  *                 workspace at least NRHS , optimally NRHS*NB
280  
281                    SCLLEN = M
282  
283                END IF
284  
285            ELSE
286  
287  *             Compute LQ factorization of sub( A )
288  
289                CALL PSGELQF ( M , N , A , IA , JA , DESCA , WORK , WORK( IPW ) ,
290       $        LWORK - LTAU , INFO )
291  
292  *             workspace at least M , optimally M*NB.
293  
294                IF( .NOT.TPSD ) THEN
295  
296  *                 underdetermined system of equations sub( A ) * X = sub( B )
297  
298  *                 B(IB : IB + M - 1 , JB : JB + NRHS - 1) := inv(L) *
299  *                 B(IB : IB + M - 1 , JB : JB + NRHS - 1)
300  
301                    CALL PSTRSM( 'Left' , 'Lower' , 'No transpose' , 'Non - unit' , M ,
302       $            NRHS , ONE , A , IA , JA , DESCA , B , IB , JB , DESCB )
303  
304  *                 B(IB + M : IB + N - 1 , JB : JB + NRHS - 1) = 0
305  
306                    CALL PSLASET ( 'All' , N - M , NRHS , ZERO , ZERO , B , IB + M , JB ,
307       $            DESCB )
308  
309  *                 B(IB : IB + N - 1 , JB : JB + NRHS - 1) := Q(1 : N ,: )' *
310  *                 B(IB : IB + M - 1 , JB : JB + NRHS - 1)
311  
312                    CALL PSORMLQ ( 'Left' , 'Transpose' , N , NRHS , M , A , IA , JA ,
313       $            DESCA , WORK , B , IB , JB , DESCB , WORK( IPW ) ,
314       $            LWORK - LTAU , INFO )
315  
316  *                 workspace at least NRHS , optimally NRHS*NB
317  
318                    SCLLEN = N
319  
320                ELSE
321  
322  *                 overdetermined system min || A' * X - B ||
323  
324  *                 B(IB : IB + N - 1 , JB : JB + NRHS - 1) := Q * B(IB : IB + N - 1 , JB : JB + NRHS - 1)
325  
326                    CALL PSORMLQ ( 'Left' , 'No transpose' , N , NRHS , M , A , IA , JA ,
327       $            DESCA , WORK , B , IB , JB , DESCB , WORK( IPW ) ,
328       $            LWORK - LTAU , INFO )
329  
330  *                 workspace at least NRHS , optimally NRHS*NB
331  
332  *                 B(IB : IB + M - 1 , JB : JB + NRHS - 1) := inv(L') *
333  *                 B(IB : IB + M - 1 , JB : JB + NRHS - 1)
334  
335                    CALL PSTRSM( 'Left' , 'Lower' , 'Transpose' , 'Non - unit' , M ,
336       $            NRHS , ONE , A , IA , JA , DESCA , B , IB , JB , DESCB )
337  
338                    SCLLEN = M
339  
340                END IF
341  
342            END IF
343  
344  *         Undo scaling
345  
346            IF( IASCL.EQ.1 ) THEN
347                CALL PSLASCL ( 'G' , ANRM , SMLNUM , SCLLEN , NRHS , B , IB , JB ,
348       $        DESCB , INFO )
349            ELSE IF( IASCL.EQ.2 ) THEN
350                CALL PSLASCL ( 'G' , ANRM , BIGNUM , SCLLEN , NRHS , B , IB , JB ,
351       $        DESCB , INFO )
352            END IF
353            IF( IBSCL.EQ.1 ) THEN
354                CALL PSLASCL ( 'G' , SMLNUM , BNRM , SCLLEN , NRHS , B , IB , JB ,
355       $        DESCB , INFO )
356            ELSE IF( IBSCL.EQ.2 ) THEN
357                CALL PSLASCL ( 'G' , BIGNUM , BNRM , SCLLEN , NRHS , B , IB , JB ,
358       $        DESCB , INFO )
359            END IF
360  
361     10 CONTINUE
362  
363        WORK( 1 ) = REAL( LWMIN )
364  
365        RETURN
366  
367  *     End of PSGELS
368  
369        END