Routine: PSGERFS()  File: SRC\psgerfs.f

 
 
# lines: 886
  # code: 886
  # comment: 0
  # blank:0
# Variables:101
# Callers:1
# Callings:3
# Words:565
# Keywords:376
 

 

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