Routine: PCPORFS()  File: SRC\pcporfs.f

 
 
# lines: 870
  # code: 870
  # comment: 0
  # blank:0
# Variables:103
# Callers:1
# Callings:3
# Words:553
# Keywords:369
 

 

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