Routine: PDPORFS()  File: SRC\pdporfs.f

 
 
# lines: 861
  # code: 861
  # comment: 0
  # blank:0
# Variables:100
# Callers:1
# Callings:3
# Words:557
# Keywords:366
 

 

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