Routine: PCGERFS()  File: SRC\pcgerfs.f

 
 
# lines: 895
  # code: 895
  # comment: 0
  # blank:0
# Variables:105
# Callers:1
# Callings:3
# Words:565
# Keywords:378
 

 

 
..
     .. Array Arguments ..
     ..
  Purpose
  =======
  PCGERFS improves the computed solution to a system of linear
  equations and provides error bounds and backward error estimates for
  the solutions.
  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
  In the following comments, sub( A ), sub( X ) and sub( B ) denote
  respectively A(IA:IA+N-1,JA:JA+N-1), X(IX:IX+N-1,JX:JX+NRHS-1) and
  B(IB:IB+N-1,JB:JB+NRHS-1).
  Arguments
  =========
  TRANS   (global input) CHARACTER*1
          Specifies the form of the system of equations.
          = 'N': sub( A ) * sub( X ) = sub( B )          (No transpose)
          = 'T': sub( A )**T * sub( X ) = sub( B )          (Transpose)
          = 'C': sub( A )**H * sub( X ) = sub( B )
                                                  (Conjugate transpose)
  N       (global input) INTEGER
          The order of the matrix sub( A ). N >= 0.
  NRHS    (global input) INTEGER
          The number of right hand sides, i.e., the number of columns
          of the matrices sub( B ) and sub( X ).  NRHS >= 0.
  A       (local input) COMPLEX pointer into the local
          memory to an array of local dimension (LLD_A,LOCc(JA+N-1)).
          This array contains the local pieces of the distributed
          matrix sub( A ).
  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.
  AF      (local input) COMPLEX pointer into the local
          memory to an array of local dimension (LLD_AF,LOCc(JA+N-1)).
          This array contains the local pieces of the distributed
          factors of the matrix sub( A ) = P * L * U as computed by
          PCGETRF.
  IAF     (global input) INTEGER
          The row index in the global array AF indicating the first
          row of sub( AF ).
  JAF     (global input) INTEGER
          The column index in the global array AF indicating the
          first column of sub( AF ).
  DESCAF  (global and local input) INTEGER array of dimension DLEN_.
          The array descriptor for the distributed matrix AF.
  IPIV    (local input) INTEGER array of dimension LOCr(M_AF)+MB_AF.
          This array contains the pivoting information as computed
          by PCGETRF. IPIV(i) -> The global row local row i
          was swapped with. This array is tied to the distributed
          matrix A.
  B       (local input) COMPLEX pointer into the local
          memory to an array of local dimension
          (LLD_B,LOCc(JB+NRHS-1)). This array contains the local
          pieces of the distributed matrix of right hand sides
          sub( B ).
  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.
  X       (local input and output) COMPLEX pointer into the
          local memory to an array of local dimension
          (LLD_X,LOCc(JX+NRHS-1)). On entry, this array contains
          the local pieces of the distributed matrix solution
          sub( X ). On exit, the improved solution vectors.
  IX      (global input) INTEGER
          The row index in the global array X indicating the first
          row of sub( X ).
  JX      (global input) INTEGER
          The column index in the global array X indicating the
          first column of sub( X ).
  DESCX   (global and local input) INTEGER array of dimension DLEN_.
          The array descriptor for the distributed matrix X.
  FERR    (local output) REAL array of local dimension
          LOCc(JB+NRHS-1).
          The estimated forward error bound for each solution vector
          of sub( X ).  If XTRUE is the true solution corresponding
          to sub( X ), FERR is an estimated upper bound for the
          magnitude of the largest element in (sub( X ) - XTRUE)
          divided by the magnitude of the largest element in sub( X ).
          The estimate is as reliable as the estimate for RCOND, and
          is almost always a slight overestimate of the true error.
          This array is tied to the distributed matrix X.
  BERR    (local output) REAL array of local dimension
          LOCc(JB+NRHS-1). The componentwise relative backward
          error of each solution vector (i.e., the smallest re-
          lative change in any entry of sub( A ) or sub( B )
          that makes sub( X ) an exact solution).
          This array is tied to the distributed matrix X.
  WORK    (local workspace/local output) COMPLEX 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 >= 2*LOCr( N + MOD(IA-1,MB_A) )
          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.
  RWORK   (local workspace/local output) REAL array,
                                                   dimension (LRWORK)
          On exit, RWORK(1) returns the minimal and optimal LRWORK.
  LRWORK  (local or global input) INTEGER
          The dimension of the array RWORK.
          LRWORK is local input and must be at least
          LRWORK >= LOCr( N + MOD(IB-1,MB_B) ).
          If LRWORK = -1, then LRWORK 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.
  Internal Parameters
  ===================
  ITMAX is the maximum number of steps of iterative refinement.
  Notes
  =====
  This routine temporarily returns when N <= 1.
  The distributed submatrices op( A ) and op( AF ) (respectively
  sub( X ) and sub( B ) ) should be distributed the same way on the
  same processes. These conditions ensure that sub( A ) and sub( AF )
  (resp. sub( X ) and sub( B ) ) are "perfectly" aligned.
  Moreover, this routine requires the distributed submatrices sub( A ),
  sub( AF ), sub( X ), and sub( B ) to be aligned on a block boundary,
  i.e., if f(x,y) = MOD( x-1, y ):
  f( IA, DESCA( MB_ ) ) = f( JA, DESCA( NB_ ) ) = 0,
  f( IAF, DESCAF( MB_ ) ) = f( JAF, DESCAF( NB_ ) ) = 0,
  f( IB, DESCB( MB_ ) ) = f( JB, DESCB( NB_ ) ) = 0, and
  f( IX, DESCX( MB_ ) ) = f( JX, DESCX( NB_ ) ) = 0.
  =====================================================================
     .. Parameters ..
 
Display dynamic version Find AutoScroll Reload FontSize: - + Hide Comments Hide Blanks Frame FullScreen MailPrint

  
001        SUBROUTINE PCGERFS( TRANS , N , NRHS , A , IA , JA , DESCA , AF , IAF ,
002       $JAF , DESCAF , IPIV , B , IB , JB , DESCB , X , IX ,
003       $JX , DESCX , FERR , BERR , WORK , LWORK , RWORK ,
004       $LRWORK , INFO )
005  
006  *     -- ScaLAPACK routine(version 1.7) --
007  *     University of Tennessee , Knoxville , Oak Ridge National Laboratory ,
008  *     and University of California , Berkeley.
009  *     November 15 , 1997
010  
011  *     .. Scalar Arguments ..
012        CHARACTER TRANS
013        INTEGER IA , IAF , IB , IX , INFO , JA , JAF , JB , JX ,
014       $LRWORK , LWORK , N , NRHS
015        INTEGER BLOCK_CYCLIC_2D , CSRC_ , CTXT_ , DLEN_ , DTYPE_ ,
016       $LLD_ , MB_ , M_ , NB_ , N_ , RSRC_
017        PARAMETER( BLOCK_CYCLIC_2D = 1 , DLEN_ = 9 , DTYPE_ = 1 ,
018       $CTXT_ = 2 , M_ = 3 , N_ = 4 , MB_ = 5 , NB_ = 6 ,
019       $RSRC_ = 7 , CSRC_ = 8 , LLD_ = 9 )
020        INTEGER ITMAX
021        PARAMETER( ITMAX = 5 )
022        REAL ZERO , RONE , TWO , THREE
023        PARAMETER( ZERO = 0.0E + 0 , RONE = 1.0E + 0 , TWO = 2.0E + 0 ,
024       $THREE = 3.0E + 0 )
025        COMPLEX ONE
026        PARAMETER( ONE =( 1.0E + 0 , 0.0E + 0 ) )
027  *     ..
028  *     .. Local Scalars ..
029        LOGICAL LQUERY , NOTRAN
030        CHARACTER TRANSN , TRANST
031        INTEGER COUNT , IACOL , IAFCOL , IAFROW , IAROW , IXBCOL ,
032       $IXBROW , IXCOL , IXROW , ICOFFA , ICOFFAF , ICOFFB ,
033       $ICOFFX , ICTXT , ICURCOL , IDUM , II , IIXB , IIW ,
034       $IOFFXB , IPB , IPR , IPV , IROFFA , IROFFAF , IROFFB ,
035       $IROFFX , IW , J , JBRHS , JJ , JJFBE , JJXB , JN , JW ,
036       $K , KASE , LDXB , LRWMIN , LWMIN , MYCOL , MYRHS ,
037       $MYROW , NP , NP0 , NPCOL , NPMOD , NPROW , NZ
038        REAL EPS , EST , LSTRES , S , SAFE1 , SAFE2 , SAFMIN
039        COMPLEX ZDUM
040  *     ..
041  *     .. Local Arrays ..
042        INTEGER DESCW( DLEN_ ) , IDUM1( 5 ) , IDUM2( 5 )
043  *     ..
044  *     .. External Functions ..
045        LOGICAL LSAME
046        INTEGER ICEIL , INDXG2P , NUMROC
047        REAL PSLAMCH
048        EXTERNAL ICEIL , INDXG2P , LSAME , NUMROC , PSLAMCH
049  *     ..
050  *     .. External Subroutines ..
051        EXTERNAL BLACS_GRIDINFO , CGEBR2D , CGEBS2D , CHK1MAT ,
052       $DESCSET , INFOG2L , PCAGEMV , PCAXPY , PCHK2MAT ,
053       $PCCOPY , PCGEMV , PCGETRS , PCLACON ,
054       $PXERBLA , SGAMX2D
055  *     ..
056  *     .. Intrinsic Functions ..
057        INTRINSIC ABS , AIMAG , CMPLX , ICHAR , MAX , MIN , MOD , REAL
058  *     ..
059  *     .. Statement Functions ..
060        REAL CABS1
061  *     ..
062  *     .. Statement Function definitions ..
063        CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
064  *     ..
065  *     .. Executable Statements ..
066  
067  *     Get grid parameters
068  
069        ICTXT = DESCA( CTXT_ )
070        CALL BLACS_GRIDINFO( ICTXT , NPROW , NPCOL , MYROW , MYCOL )
071  
072  *     Test the input parameters.
073  
074        NOTRAN = LSAME( TRANS , 'N' )
075  
076        INFO = 0
077        IF( NPROW.EQ. - 1 ) THEN
078            INFO = - (700 + CTXT_)
079        ELSE
080            CALL CHK1MAT( N , 2 , N , 2 , IA , JA , DESCA , 7 , INFO )
081            CALL CHK1MAT( N , 2 , N , 2 , IAF , JAF , DESCAF , 11 , INFO )
082            CALL CHK1MAT( N , 2 , NRHS , 3 , IB , JB , DESCB , 16 , INFO )
083            CALL CHK1MAT( N , 2 , NRHS , 3 , IX , JX , DESCX , 20 , INFO )
084            IF( INFO.EQ.0 ) THEN
085                IROFFA = MOD( IA - 1 , DESCA( MB_ ) )
086                ICOFFA = MOD( JA - 1 , DESCA( NB_ ) )
087                IROFFAF = MOD( IAF - 1 , DESCAF( MB_ ) )
088                ICOFFAF = MOD( JAF - 1 , DESCAF( NB_ ) )
089                IROFFB = MOD( IB - 1 , DESCB( MB_ ) )
090                ICOFFB = MOD( JB - 1 , DESCB( NB_ ) )
091                IROFFX = MOD( IX - 1 , DESCX( MB_ ) )
092                ICOFFX = MOD( JX - 1 , DESCX( NB_ ) )
093                IAROW = INDXG2P( IA , DESCA( MB_ ) , MYROW , DESCA( RSRC_ ) ,
094       $        NPROW )
095                IAFCOL = INDXG2P( JAF , DESCAF( NB_ ) , MYCOL ,
096       $        DESCAF( CSRC_ ) , NPCOL )
097                IAFROW = INDXG2P( IAF , DESCAF( MB_ ) , MYROW ,
098       $        DESCAF( RSRC_ ) , NPROW )
099                IACOL = INDXG2P( JA , DESCA( NB_ ) , MYCOL , DESCA( CSRC_ ) ,
100       $        NPCOL )
101                CALL INFOG2L( IB , JB , DESCB , NPROW , NPCOL , MYROW , MYCOL ,
102       $        IIXB , JJXB , IXBROW , IXBCOL )
103                IXROW = INDXG2P( IX , DESCX( MB_ ) , MYROW , DESCX( RSRC_ ) ,
104       $        NPROW )
105                IXCOL = INDXG2P( JX , DESCX( NB_ ) , MYCOL , DESCX( CSRC_ ) ,
106       $        NPCOL )
107                NPMOD = NUMROC( N + IROFFA , DESCA( MB_ ) , MYROW , IAROW ,
108       $        NPROW )
109                LWMIN = 2 * NPMOD
110                LRWMIN = NPMOD
111                WORK( 1 ) = CMPLX( REAL( LWMIN ) )
112                RWORK( 1 ) = REAL( LRWMIN )
113                LQUERY =( LWORK.EQ. - 1 .OR. LRWORK.EQ. - 1 )
114  
115                IF(( .NOT.NOTRAN ) .AND.( .NOT.LSAME( TRANS , 'T' ) ) .AND.
116       $( .NOT.LSAME( TRANS , 'C' ) ) ) THEN
117                INFO = - 1
118            ELSE IF( N.LT.0 ) THEN
119                INFO = - 2
120            ELSE IF( NRHS.LT.0 ) THEN
121                INFO = - 3
122            ELSE IF( IROFFA.NE.0 ) THEN
123                INFO = - 5
124            ELSE IF( ICOFFA.NE.0 ) THEN
125                INFO = - 6
126            ELSE IF( DESCA( MB_ ).NE.DESCA( NB_ ) ) THEN
127                INFO = - ( 700 + NB_ )
128            ELSE IF( DESCA( MB_ ).NE.DESCAF( MB_ ) ) THEN
129                INFO = - ( 1100 + MB_ )
130            ELSE IF( IROFFAF.NE.0 .OR. IAROW.NE.IAFROW ) THEN
131                INFO = - 9
132            ELSE IF( DESCA( NB_ ).NE.DESCAF( NB_ ) ) THEN
133                INFO = - ( 1100 + NB_ )
134            ELSE IF( ICOFFAF.NE.0 .OR. IACOL.NE.IAFCOL ) THEN
135                INFO = - 10
136            ELSE IF( ICTXT.NE.DESCAF( CTXT_ ) ) THEN
137                INFO = - ( 1100 + CTXT_ )
138            ELSE IF( IROFFA.NE.IROFFB .OR. IAROW.NE.IXBROW ) THEN
139                INFO = - 14
140            ELSE IF( DESCA( MB_ ).NE.DESCB( MB_ ) ) THEN
141                INFO = - ( 1600 + MB_ )
142            ELSE IF( ICTXT.NE.DESCB( CTXT_ ) ) THEN
143                INFO = - ( 1600 + CTXT_ )
144            ELSE IF( DESCB( MB_ ).NE.DESCX( MB_ ) ) THEN
145                INFO = - ( 2000 + MB_ )
146            ELSE IF( IROFFX.NE.0 .OR. IXBROW.NE.IXROW ) THEN
147                INFO = - 18
148            ELSE IF( DESCB( NB_ ).NE.DESCX( NB_ ) ) THEN
149                INFO = - ( 2000 + NB_ )
150            ELSE IF( ICOFFB.NE.ICOFFX .OR. IXBCOL.NE.IXCOL ) THEN
151                INFO = - 19
152            ELSE IF( ICTXT.NE.DESCX( CTXT_ ) ) THEN
153                INFO = - ( 2000 + CTXT_ )
154            ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
155                INFO = - 24
156            ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
157                INFO = - 26
158            END IF
159        END IF
160  
161        IF( NOTRAN ) THEN
162            IDUM1( 1 ) = ICHAR( 'N' )
163        ELSE IF( LSAME( TRANS , 'T' ) ) THEN
164            IDUM1( 1 ) = ICHAR( 'T' )
165        ELSE
166            IDUM1( 1 ) = ICHAR( 'C' )
167        END IF
168        IDUM2( 1 ) = 1
169        IDUM1( 2 ) = N
170        IDUM2( 2 ) = 2
171        IDUM1( 3 ) = NRHS
172        IDUM2( 3 ) = 3
173        IF( LWORK.EQ. - 1 ) THEN
174            IDUM1( 4 ) = - 1
175        ELSE
176            IDUM1( 4 ) = 1
177        END IF
178        IDUM2( 4 ) = 24
179        IF( LRWORK.EQ. - 1 ) THEN
180            IDUM1( 5 ) = - 1
181        ELSE
182            IDUM1( 5 ) = 1
183        END IF
184        IDUM2( 5 ) = 26
185        CALL PCHK2MAT( N , 2 , N , 2 , IA , JA , DESCA , 7 , N , 2 , N , 2 , IAF ,
186       $JAF , DESCAF , 11 , 5 , IDUM1 , IDUM2 , INFO )
187        CALL PCHK2MAT( N , 2 , NRHS , 3 , IB , JB , DESCB , 16 , N , 2 , NRHS , 3 ,
188       $IX , JX , DESCX , 20 , 5 , IDUM1 , IDUM2 , INFO )
189        END IF
190        IF( INFO.NE.0 ) THEN
191            CALL PXERBLA( ICTXT , 'PCGERFS' , - INFO )
192            RETURN
193        ELSE IF( LQUERY ) THEN
194            RETURN
195        END IF
196  
197        JJFBE = JJXB
198        MYRHS = NUMROC( JB + NRHS - 1 , DESCB( NB_ ) , MYCOL , DESCB( CSRC_ ) ,
199       $NPCOL )
200  
201  *     Quick return if possible
202  
203        IF( N.LE.1 .OR. NRHS.EQ.0 ) THEN
204            DO 10 JJ = JJFBE , MYRHS
205                FERR( JJ ) = ZERO
206                BERR( JJ ) = ZERO
207     10     CONTINUE
208            RETURN
209        END IF
210  
211        IF( NOTRAN ) THEN
212            TRANSN = 'N'
213            TRANST = 'C'
214        ELSE
215            TRANSN = 'C'
216            TRANST = 'N'
217        END IF
218  
219        NP0 = NUMROC( N + IROFFB , DESCB( MB_ ) , MYROW , IXBROW , NPROW )
220        CALL DESCSET( DESCW , N + IROFFB , 1 , DESCA( MB_ ) , 1 , IXBROW , IXBCOL ,
221       $ICTXT , MAX( 1 , NP0 ) )
222        IPB = 1
223        IPR = 1
224        IPV = IPR + NP0
225        IF( MYROW.EQ.IXBROW ) THEN
226            IIW = 1 + IROFFB
227            NP = NP0 - IROFFB
228        ELSE
229            IIW = 1
230            NP = NP0
231        END IF
232        IW = 1 + IROFFB
233        JW = 1
234        LDXB = DESCB( LLD_ )
235        IOFFXB =( JJXB - 1 )*LDXB
236  
237  *     NZ = 1 + maximum number of nonzero entries in each row of sub( A )
238  
239        NZ = N + 1
240        EPS = PSLAMCH( ICTXT , 'Epsilon' )
241        SAFMIN = PSLAMCH( ICTXT , 'Safe minimum' )
242        SAFE1 = NZ*SAFMIN
243        SAFE2 = SAFE1 / EPS
244        JN = MIN( ICEIL( JB , DESCB( NB_ ) ) * DESCB( NB_ ) , JB + NRHS - 1 )
245  
246  *     Handle first block separately
247  
248        JBRHS = JN - JB + 1
249        DO 100 K = 0 , JBRHS - 1
250  
251            COUNT = 1
252            LSTRES = THREE
253     20 CONTINUE
254  
255  *     Loop until stopping criterion is satisfied.
256  
257  *     Compute residual R = sub(B) - op(sub(A)) * sub(X) ,
258  *     where op(sub(A)) = sub(A) , or sub(A)'(A**T or A**H) ,
259  *     depending on TRANS.
260  
261        CALL PCCOPY( N , B , IB , JB + K , DESCB , 1 , WORK( IPR ) , IW , JW ,
262       $DESCW , 1 )
263        CALL PCGEMV( TRANS , N , N , - ONE , A , IA , JA , DESCA , X , IX ,
264       $JX + K , DESCX , 1 , ONE , WORK( IPR ) , IW , JW ,
265       $DESCW , 1 )
266  
267  *     Compute componentwise relative backward error from formula
268  
269  *     max(i)( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
270  
271  *     where abs(Z) is the componentwise absolute value of the
272  *     matrix or vector Z. If the i - th component of the
273  *     denominator is less than SAFE2 , then SAFE1 is added to the
274  *     i - th components of the numerator and denominator before
275  *     dividing.
276  
277        IF( MYCOL.EQ.IXBCOL ) THEN
278            IF( NP.GT.0 ) THEN
279                DO 30 II = IIXB , IIXB + NP - 1
280                    RWORK( IIW + II - IIXB ) = CABS1( B( II + IOFFXB ) )
281     30         CONTINUE
282            END IF
283        END IF
284  
285        CALL PCAGEMV( TRANS , N , N , RONE , A , IA , JA , DESCA , X , IX , JX + K ,
286       $DESCX , 1 , RONE , RWORK( IPB ) , IW , JW , DESCW , 1 )
287  
288        S = ZERO
289        IF( MYCOL.EQ.IXBCOL ) THEN
290            IF( NP.GT.0 ) THEN
291                DO 40 II = IIW - 1 , IIW + NP - 2
292                    IF( RWORK( IPB + II ).GT.SAFE2 ) THEN
293                        S = MAX( S , CABS1( WORK( IPR + II ) ) /
294       $                RWORK( IPB + II ) )
295                    ELSE
296                        S = MAX( S ,( CABS1( WORK( IPR + II ) ) + SAFE1 ) /
297       $( RWORK( IPB + II ) + SAFE1 ) )
298                    END IF
299     40         CONTINUE
300            END IF
301        END IF
302  
303        CALL SGAMX2D( ICTXT , 'All' , ' ' , 1 , 1 , S , 1 , IDUM , IDUM , 1 ,
304       $- 1 , MYCOL )
305        IF( MYCOL.EQ.IXBCOL )
306       $    BERR( JJFBE ) = S
307  
308  *         Test stopping criterion. Continue iterating if
309  *         1) The residual BERR(J + K) is larger than machine epsilon ,
310  *         and
311  *         2) BERR(J + K) decreased by at least a factor of 2 during the
312  *         last iteration , and
313  *         3) At most ITMAX iterations tried.
314  
315            IF( S.GT.EPS .AND. TWO*S.LE.LSTRES .AND. COUNT.LE.ITMAX ) THEN
316  
317  *             Update solution and try again.
318  
319                CALL PCGETRS ( TRANS , N , 1 , AF , IAF , JAF , DESCAF , IPIV ,
320       $        WORK( IPR ) , IW , JW , DESCW , INFO )
321                CALL PCAXPY( N , ONE , WORK( IPR ) , IW , JW , DESCW , 1 , X , IX ,
322       $        JX + K , DESCX , 1 )
323                LSTRES = S
324                COUNT = COUNT + 1
325                GO TO 20
326            END IF
327  
328  *         Bound error from formula
329  
330  *         norm(sub(X) - XTRUE) / norm(sub(X)) .le. FERR =
331  *         norm( abs(inv(op(sub(A))))*
332  *         ( abs(R) + NZ*EPS*(
333  *         abs(op(sub(A)))*abs(sub(X)) + abs(sub(B))))) / norm(sub(X))
334  
335  *         where
336  *         norm(Z) is the magnitude of the largest component of Z
337  *         inv(op(sub(A))) is the inverse of op(sub(A))
338  *         abs(Z) is the componentwise absolute value of the matrix
339  *         or vector Z
340  *         NZ is the maximum number of nonzeros in any row of sub(A) ,
341  *         plus 1
342  *         EPS is machine epsilon
343  
344  *         The i - th component of
345  *         abs(R) + NZ*EPS*(abs(op(sub(A)))*abs(sub(X)) + abs(sub(B)))
346  *         is incremented by SAFE1 if the i - th component of
347  *         abs(op(sub(A)))*abs(sub(X)) + abs(sub(B)) is less than
348  *         SAFE2.
349  
350  *         Use PCLACON to estimate the infinity - norm of the matrix
351  *         inv(op(sub(A))) * diag(W) , where
352  *         W = abs(R) + NZ*EPS*(abs(op(sub(A)))*abs(sub(X)) + abs(sub(B))).
353  
354            IF( MYCOL.EQ.IXBCOL ) THEN
355                IF( NP.GT.0 ) THEN
356                    DO 50 II = IIW - 1 , IIW + NP - 2
357                        IF( RWORK( IPB + II ).GT.SAFE2 ) THEN
358                            RWORK( IPB + II ) = CABS1( WORK( IPR + II ) ) +
359       $                    NZ*EPS*RWORK( IPB + II )
360                        ELSE
361                            RWORK( IPB + II ) = CABS1( WORK( IPR + II ) ) +
362       $                    NZ*EPS*RWORK( IPB + II ) + SAFE1
363                        END IF
364     50             CONTINUE
365                END IF
366            END IF
367  
368            KASE = 0
369     60 CONTINUE
370        IF( MYCOL.EQ.IXBCOL ) THEN
371            CALL CGEBS2D( ICTXT , 'Rowwise' , ' ' , NP , 1 , WORK( IPR ) ,
372       $    DESCW( LLD_ ) )
373        ELSE
374            CALL CGEBR2D( ICTXT , 'Rowwise' , ' ' , NP , 1 , WORK( IPR ) ,
375       $    DESCW( LLD_ ) , MYROW , IXBCOL )
376        END IF
377        DESCW( CSRC_ ) = MYCOL
378        CALL PCLACON ( N , WORK( IPV ) , IW , JW , DESCW , WORK( IPR ) ,
379       $IW , JW , DESCW , EST , KASE )
380        DESCW( CSRC_ ) = IXBCOL
381  
382        IF( KASE.NE.0 ) THEN
383            IF( KASE.EQ.1 ) THEN
384  
385  *             Multiply by diag(W)*inv(op(sub(A))').
386  
387                CALL PCGETRS ( TRANST , N , 1 , AF , IAF , JAF , DESCAF ,
388       $        IPIV , WORK( IPR ) , IW , JW , DESCW , INFO )
389  
390                IF( MYCOL.EQ.IXBCOL ) THEN
391                    IF( NP.GT.0 ) THEN
392                        DO 70 II = IIW - 1 , IIW + NP - 2
393                            WORK( IPR + II ) = RWORK( IPB + II )*WORK( IPR + II )
394     70                 CONTINUE
395                    END IF
396                END IF
397            ELSE
398  
399  *             Multiply by inv(op(sub(A)))*diag(W).
400  
401                IF( MYCOL.EQ.IXBCOL ) THEN
402                    IF( NP.GT.0 ) THEN
403                        DO 80 II = IIW - 1 , IIW + NP - 2
404                            WORK( IPR + II ) = RWORK( IPB + II )*WORK( IPR + II )
405     80                 CONTINUE
406                    END IF
407                END IF
408  
409                CALL PCGETRS ( TRANSN , N , 1 , AF , IAF , JAF , DESCAF ,
410       $        IPIV , WORK( IPR ) , IW , JW , DESCW , INFO )
411            END IF
412            GO TO 60
413        END IF
414  
415  *     Normalize error.
416  
417        LSTRES = ZERO
418        IF( MYCOL.EQ.IXBCOL ) THEN
419            IF( NP.GT.0 ) THEN
420                DO 90 II = IIXB , IIXB + NP - 1
421                    LSTRES = MAX( LSTRES , CABS1( X( IOFFXB + II ) ) )
422     90         CONTINUE
423            END IF
424            CALL SGAMX2D( ICTXT , 'Column' , ' ' , 1 , 1 , LSTRES , 1 , IDUM ,
425       $    IDUM , 1 , - 1 , MYCOL )
426            IF( LSTRES.NE.ZERO )
427       $        FERR( JJFBE ) = EST / LSTRES
428  
429                JJXB = JJXB + 1
430                JJFBE = JJFBE + 1
431                IOFFXB = IOFFXB + LDXB
432  
433            END IF
434  
435    100 CONTINUE
436  
437        ICURCOL = MOD( IXBCOL + 1 , NPCOL )
438  
439  *     Do for each right hand side
440  
441        DO 200 J = JN + 1 , JB + NRHS - 1 , DESCB( NB_ )
442            JBRHS = MIN( JB + NRHS - J , DESCB( NB_ ) )
443            DESCW( CSRC_ ) = ICURCOL
444  
445            DO 190 K = 0 , JBRHS - 1
446  
447                COUNT = 1
448                LSTRES = THREE
449    110 CONTINUE
450  
451  *     Loop until stopping criterion is satisfied.
452  
453  *     Compute residual R = sub(B) - op(sub(A)) * sub(X) ,
454  *     where op(sub(A)) = sub(A) , or sub(A)'(A**T or A**H) ,
455  *     depending on TRANS.
456  
457        CALL PCCOPY( N , B , IB , J + K , DESCB , 1 , WORK( IPR ) , IW , JW ,
458       $DESCW , 1 )
459        CALL PCGEMV( TRANS , N , N , - ONE , A , IA , JA , DESCA , X ,
460       $IX , J + K , DESCX , 1 , ONE , WORK( IPR ) , IW , JW ,
461       $DESCW , 1 )
462  
463  *     Compute componentwise relative backward error from formula
464  
465  *     max(i)(abs(R(i)) / (abs(op(sub(A)))*abs(sub(X)) +
466  *     abs(sub(B)))(i))
467  
468  *     where abs(Z) is the componentwise absolute value of the
469  *     matrix or vector Z. If the i - th component of the
470  *     denominator is less than SAFE2 , then SAFE1 is added to the
471  *     i - th components of the numerator and denominator before
472  *     dividing.
473  
474        IF( MYCOL.EQ.ICURCOL ) THEN
475            IF( NP.GT.0 ) THEN
476                DO 120 II = IIXB , IIXB + NP - 1
477                    RWORK( IIW + II - IIXB ) = CABS1( B( II + IOFFXB ) )
478    120         CONTINUE
479            END IF
480        END IF
481  
482        CALL PCAGEMV( TRANS , N , N , RONE , A , IA , JA , DESCA , X , IX ,
483       $J + K , DESCX , 1 , RONE , RWORK( IPB ) , IW , JW ,
484       $DESCW , 1 )
485  
486        S = ZERO
487        IF( MYCOL.EQ.ICURCOL ) THEN
488            IF( NP.GT.0 )THEN
489                DO 130 II = IIW - 1 , IIW + NP - 2
490                    IF( RWORK( IPB + II ).GT.SAFE2 ) THEN
491                        S = MAX( S , CABS1( WORK( IPR + II ) ) /
492       $                RWORK( IPB + II ) )
493                    ELSE
494                        S = MAX( S ,( CABS1( WORK( IPR + II ) ) + SAFE1 ) /
495       $( RWORK( IPB + II ) + SAFE1 ) )
496                    END IF
497    130         CONTINUE
498            END IF
499        END IF
500  
501        CALL SGAMX2D( ICTXT , 'All' , ' ' , 1 , 1 , S , 1 , IDUM , IDUM , 1 ,
502       $- 1 , MYCOL )
503        IF( MYCOL.EQ.ICURCOL )
504       $    BERR( JJFBE ) = S
505  
506  *         Test stopping criterion. Continue iterating if
507  *         1) The residual BERR(J + K) is larger than machine epsilon ,
508  *         and
509  *         2) BERR(J + K) decreased by at least a factor of 2 during
510  *         the last iteration , and
511  *         3) At most ITMAX iterations tried.
512  
513            IF( S.GT.EPS .AND. TWO*S.LE.LSTRES .AND.
514       $        COUNT.LE.ITMAX ) THEN
515  
516  *             Update solution and try again.
517  
518                CALL PCGETRS ( TRANS , N , 1 , AF , IAF , JAF , DESCAF , IPIV ,
519       $        WORK( IPR ) , IW , JW , DESCW , INFO )
520                CALL PCAXPY( N , ONE , WORK( IPR ) , IW , JW , DESCW , 1 , X ,
521       $        IX , J + K , DESCX , 1 )
522                LSTRES = S
523                COUNT = COUNT + 1
524                GO TO 110
525            END IF
526  
527  *         Bound error from formula
528  
529  *         norm(sub(X) - XTRUE) / norm(sub(X)) .le. FERR =
530  *         norm( abs(inv(op(sub(A))))*
531  *         ( abs(R) + NZ*EPS*(
532  *         abs(op(sub(A)))*abs(sub(X)) + abs(sub(B))))) / norm(sub(X))
533  
534  *         where
535  *         norm(Z) is the magnitude of the largest component of Z
536  *         inv(op(sub(A))) is the inverse of op(sub(A))
537  *         abs(Z) is the componentwise absolute value of the matrix
538  *         or vector Z
539  *         NZ is the maximum number of nonzeros in any row of sub(A) ,
540  *         plus 1
541  *         EPS is machine epsilon
542  
543  *         The i - th component of
544  *         abs(R) + NZ*EPS*(abs(op(sub(A)))*abs(sub(X)) + abs(sub(B)))
545  *         is incremented by SAFE1 if the i - th component of
546  *         abs(op(sub(A)))*abs(sub(X)) + abs(sub(B)) is less than
547  *         SAFE2.
548  
549  *         Use PCLACON to estimate the infinity - norm of the matrix
550  *         inv(op(sub(A))) * diag(W) , where
551  *         W = abs(R) + NZ*EPS*(abs(op(sub(A)))*abs(sub(X)) + abs(sub(B))).
552  
553            IF( MYCOL.EQ.ICURCOL ) THEN
554                IF( NP.GT.0 ) THEN
555                    DO 140 II = IIW - 1 , IIW + NP - 2
556                        IF( RWORK( IPB + II ).GT.SAFE2 ) THEN
557                            RWORK( IPB + II ) = CABS1( WORK( IPR + II ) ) +
558       $                    NZ*EPS*RWORK( IPB + II )
559                        ELSE
560                            RWORK( IPB + II ) = CABS1( WORK( IPR + II ) ) +
561       $                    NZ*EPS*RWORK( IPB + II ) + SAFE1
562                        END IF
563    140             CONTINUE
564                END IF
565            END IF
566  
567            KASE = 0
568    150 CONTINUE
569        IF( MYCOL.EQ.ICURCOL ) THEN
570            CALL CGEBS2D( ICTXT , 'Rowwise' , ' ' , NP , 1 , WORK( IPR ) ,
571       $    DESCW( LLD_ ) )
572        ELSE
573            CALL CGEBR2D( ICTXT , 'Rowwise' , ' ' , NP , 1 , WORK( IPR ) ,
574       $    DESCW( LLD_ ) , MYROW , ICURCOL )
575        END IF
576        DESCW( CSRC_ ) = MYCOL
577        CALL PCLACON ( N , WORK( IPV ) , IW , JW , DESCW , WORK( IPR ) ,
578       $IW , JW , DESCW , EST , KASE )
579        DESCW( CSRC_ ) = ICURCOL
580  
581        IF( KASE.NE.0 ) THEN
582            IF( KASE.EQ.1 ) THEN
583  
584  *             Multiply by diag(W)*inv(op(sub(A))').
585  
586                CALL PCGETRS ( TRANST , N , 1 , AF , IAF , JAF , DESCAF ,
587       $        IPIV , WORK( IPR ) , IW , JW , DESCW , INFO )
588  
589                IF( MYCOL.EQ.ICURCOL ) THEN
590                    IF( NP.GT.0 ) THEN
591                        DO 160 II = IIW - 1 , IIW + NP - 2
592                            WORK( IPR + II ) = RWORK( IPB + II )*
593       $                    WORK( IPR + II )
594    160                 CONTINUE
595                    END IF
596                END IF
597            ELSE
598  
599  *             Multiply by inv(op(sub(A)))*diag(W).
600  
601                IF( MYCOL.EQ.ICURCOL ) THEN
602                    IF( NP.GT.0 ) THEN
603                        DO 170 II = IIW - 1 , IIW + NP - 2
604                            WORK( IPR + II ) = RWORK( IPB + II )*
605       $                    WORK( IPR + II )
606    170                 CONTINUE
607                    END IF
608                END IF
609  
610                CALL PCGETRS ( TRANSN , N , 1 , AF , IAF , JAF , DESCAF ,
611       $        IPIV , WORK( IPR ) , IW , JW , DESCW ,
612       $        INFO )
613            END IF
614            GO TO 150
615        END IF
616  
617  *     Normalize error.
618  
619        LSTRES = ZERO
620        IF( MYCOL.EQ.ICURCOL ) THEN
621            IF( NP.GT.0 ) THEN
622                DO 180 II = IIXB , IIXB + NP - 1
623                    LSTRES = MAX( LSTRES , CABS1( X( IOFFXB + II ) ) )
624    180         CONTINUE
625            END IF
626            CALL SGAMX2D( ICTXT , 'Column' , ' ' , 1 , 1 , LSTRES ,
627       $    1 , IDUM , IDUM , 1 , - 1 , MYCOL )
628            IF( LSTRES.NE.ZERO )
629       $        FERR( JJFBE ) = EST / LSTRES
630  
631                JJXB = JJXB + 1
632                JJFBE = JJFBE + 1
633                IOFFXB = IOFFXB + LDXB
634  
635            END IF
636  
637    190 CONTINUE
638  
639        ICURCOL = MOD( ICURCOL + 1 , NPCOL )
640  
641    200 CONTINUE
642  
643        WORK( 1 ) = CMPLX( REAL( LWMIN ) )
644        RWORK( 1 ) = REAL( LRWMIN )
645  
646        RETURN
647  
648  *     End of PCGERFS
649  
650        END