Routine: PZGELS()  File: SRC\pzgels.f

 
 
# lines: 593
  # code: 593
  # comment: 0
  # blank:0
# Variables:70
# Callers:0
# Callings:9
# Words:295
# Keywords:171
 

 

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