Routine: PDTRRFS()  File: SRC\pdtrrfs.f

 
 
# lines: 796
  # code: 796
  # comment: 0
  # blank:0
# Variables:94
# Callers:0
# Callings:2
# Words:535
# Keywords:349
 

 

 
..
     .. Array Arguments ..
     ..
  Purpose
  =======
  PDTRRFS provides error bounds and backward error estimates for the
  solution to a system of linear equations with a triangular
  coefficient matrix.
  The solution matrix X must be computed by PDTRTRS or some other
  means before entering this routine.  PDTRRFS does not do iterative
  refinement because doing so cannot improve the backward error.
  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
          = 'U':  sub( A ) is upper triangular;
          = 'L':  sub( A ) is lower triangular.
  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)
  DIAG    (global input) CHARACTER*1
          = 'N':  sub( A ) is non-unit triangular;
          = 'U':  sub( A ) is unit 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 original triangular
          distributed matrix sub( A ).
          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.
          If DIAG = 'U', the diagonal elements of sub( A ) are also
          not referenced and are assumed to be 1.
  IA      (global input) INTEGER
          The row index in the global array A indicating the first
          row of sub( A ).
  JA      (global input) INTEGER
          The column index in the global array A indicating the
          first column of sub( A ).
  DESCA   (global and local input) INTEGER array of dimension DLEN_.
          The array descriptor for the distributed matrix A.
  B       (local input) 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 ).
  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 bounds for
          each solution vector of sub( X ).  If XTRUE is the true
          solution, FERR bounds the magnitude of the largest entry
          in (sub( X ) - XTRUE) divided by the magnitude of the
          largest entry 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.
  Notes
  =====
  This routine temporarily returns when N <= 1.
  The distributed submatrices sub( X ) and sub( B ) should be
  distributed the same way on the same processes.  These conditions
  ensure that sub( X ) and sub( B ) are "perfectly" aligned.
  Moreover, this routine requires the distributed submatrices sub( A ),
  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( 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 PDTRRFS( UPLO , TRANS , DIAG , N , NRHS , A , IA , JA , DESCA ,
002       $B , IB , JB , DESCB , X , IX , JX , DESCX , FERR ,
003       $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 DIAG , TRANS , UPLO
012        INTEGER INFO , IA , IB , IX , JA , JB , JX , LIWORK , LWORK ,
013       $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        DOUBLE PRECISION ZERO , ONE
020        PARAMETER( ZERO = 0.0D + 0 , ONE = 1.0D + 0 )
021  *     ..
022  *     .. Local Scalars ..
023        LOGICAL LQUERY , NOTRAN , NOUNIT , UPPER
024        CHARACTER TRANST
025        INTEGER IAROW , IXBCOL , IXBROW , IXCOL , IXROW , ICOFFA ,
026       $ICOFFB , ICOFFX , ICTXT , ICURCOL , IDUM , II , IIXB ,
027       $IIW , IOFFXB , IPB , IPR , IPV , IROFFA , IROFFB ,
028       $IROFFX , IW , J , JBRHS , JJ , JJFBE , JJXB , JN , JW ,
029       $K , KASE , LDXB , LIWMIN , LWMIN , MYCOL , MYRHS ,
030       $MYROW , NP , NP0 , NPCOL , NPMOD , NPROW , NZ
031        DOUBLE PRECISION EPS , EST , LSTRES , S , SAFE1 , SAFE2 , SAFMIN
032  *     ..
033  *     .. Local Arrays ..
034        INTEGER DESCW( DLEN_ ) , IDUM1( 5 ) , IDUM2( 5 )
035  *     ..
036  *     .. External Functions ..
037        LOGICAL LSAME
038        INTEGER ICEIL , INDXG2P , NUMROC
039        DOUBLE PRECISION PDLAMCH
040        EXTERNAL ICEIL , INDXG2P , LSAME , NUMROC , PDLAMCH
041  *     ..
042  *     .. External Subroutines ..
043        EXTERNAL BLACS_GRIDINFO , CHK1MAT , DESCSET , DGAMX2D ,
044       $DGEBR2D , DGEBS2D , INFOG2L , PCHK1MAT ,
045       $PCHK2MAT , PDATRMV , PDAXPY , PDCOPY ,
046       $PDLACON , PDTRMV , PDTRSV , PXERBLA
047  *     ..
048  *     .. Intrinsic Functions ..
049        INTRINSIC ABS , DBLE , ICHAR , MAX , MIN , MOD
050  *     ..
051  *     .. Executable Statements ..
052  
053  *     Get grid parameters
054  
055        ICTXT = DESCA( CTXT_ )
056        CALL BLACS_GRIDINFO( ICTXT , NPROW , NPCOL , MYROW , MYCOL )
057  
058  *     Test the input parameters.
059  
060        INFO = 0
061        IF( NPROW.EQ. - 1 ) THEN
062            INFO = - ( 900 + CTXT_ )
063        ELSE
064            CALL CHK1MAT( N , 4 , N , 4 , IA , JA , DESCA , 9 , INFO )
065            CALL CHK1MAT( N , 4 , NRHS , 5 , IB , JB , DESCB , 13 , INFO )
066            CALL CHK1MAT( N , 4 , NRHS , 5 , IX , JX , DESCX , 17 , INFO )
067            IF( INFO.EQ.0 ) THEN
068                UPPER = LSAME( UPLO , 'U' )
069                NOTRAN = LSAME( TRANS , 'N' )
070                NOUNIT = LSAME( DIAG , 'N' )
071                IROFFA = MOD( IA - 1 , DESCA( MB_ ) )
072                ICOFFA = MOD( JA - 1 , DESCA( NB_ ) )
073                IROFFB = MOD( IB - 1 , DESCB( MB_ ) )
074                ICOFFB = MOD( JB - 1 , DESCB( NB_ ) )
075                IROFFX = MOD( IX - 1 , DESCX( MB_ ) )
076                ICOFFX = MOD( JX - 1 , DESCX( NB_ ) )
077                IAROW = INDXG2P( IA , DESCA( MB_ ) , MYROW , DESCA( RSRC_ ) ,
078       $        NPROW )
079                CALL INFOG2L( IB , JB , DESCB , NPROW , NPCOL , MYROW , MYCOL ,
080       $        IIXB , JJXB , IXBROW , IXBCOL )
081                IXROW = INDXG2P( IX , DESCX( MB_ ) , MYROW , DESCX( RSRC_ ) ,
082       $        NPROW )
083                IXCOL = INDXG2P( JX , DESCX( NB_ ) , MYCOL , DESCX( CSRC_ ) ,
084       $        NPCOL )
085                NPMOD = NUMROC( N + IROFFA , DESCA( MB_ ) , MYROW , IAROW ,
086       $        NPROW )
087                LWMIN = 3*NPMOD
088                WORK( 1 ) = DBLE( LWMIN )
089                LIWMIN = NPMOD
090                IWORK( 1 ) = LIWMIN
091                LQUERY =( LWORK.EQ. - 1 .OR. LIWORK.EQ. - 1 )
092  
093                IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO , 'L' ) ) THEN
094                    INFO = - 1
095                ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS , 'T' ) .AND.
096       $            .NOT.LSAME( TRANS , 'C' ) ) THEN
097                    INFO = - 2
098                ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG , 'U' ) ) THEN
099                    INFO = - 3
100                ELSE IF( N.LT.0 ) THEN
101                    INFO = - 4
102                ELSE IF( NRHS.LT.0 ) THEN
103                    INFO = - 5
104                ELSE IF( IROFFA.NE.0 ) THEN
105                    INFO = - 7
106                ELSE IF( ICOFFA.NE.0 ) THEN
107                    INFO = - 8
108                ELSE IF( DESCA( MB_ ).NE.DESCA( NB_ ) ) THEN
109                    INFO = - ( 900 + NB_ )
110                ELSE IF( IROFFA.NE.IROFFB .OR. IAROW.NE.IXBROW ) THEN
111                    INFO = - 11
112                ELSE IF( DESCA( MB_ ).NE.DESCB( MB_ ) ) THEN
113                    INFO = - ( 1300 + MB_ )
114                ELSE IF( ICTXT.NE.DESCB( CTXT_ ) ) THEN
115                    INFO = - ( 1300 + CTXT_ )
116                ELSE IF( IROFFX.NE.0 .OR. IXBROW.NE.IXROW ) THEN
117                    INFO = - 15
118                ELSE IF( ICOFFB.NE.ICOFFX .OR. IXBCOL.NE.IXCOL ) THEN
119                    INFO = - 16
120                ELSE IF( DESCB( MB_ ).NE.DESCX( MB_ ) ) THEN
121                    INFO = - ( 1700 + MB_ )
122                ELSE IF( DESCB( NB_ ).NE.DESCX( NB_ ) ) THEN
123                    INFO = - ( 1700 + NB_ )
124                ELSE IF( ICTXT.NE.DESCX( CTXT_ ) ) THEN
125                    INFO = - ( 1700 + CTXT_ )
126                ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
127                    INFO = - 21
128                ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
129                    INFO = - 23
130                END IF
131            END IF
132  
133            IF( UPPER ) THEN
134                IDUM1( 1 ) = ICHAR( 'U' )
135            ELSE
136                IDUM1( 1 ) = ICHAR( 'L' )
137            END IF
138            IDUM2( 1 ) = 1
139            IF( NOTRAN ) THEN
140                IDUM1( 2 ) = ICHAR( 'N' )
141            ELSE IF( LSAME( TRANS , 'T' ) ) THEN
142                IDUM1( 2 ) = ICHAR( 'T' )
143            ELSE
144                IDUM1( 2 ) = ICHAR( 'C' )
145            END IF
146            IDUM2( 2 ) = 2
147            IF( NOUNIT ) THEN
148                IDUM1( 3 ) = ICHAR( 'N' )
149            ELSE
150                IDUM1( 3 ) = ICHAR( 'U' )
151            END IF
152            IDUM2( 3 ) = 3
153            IF( LWORK.EQ. - 1 ) THEN
154                IDUM1( 4 ) = - 1
155            ELSE
156                IDUM1( 4 ) = 1
157            END IF
158            IDUM2( 4 ) = 21
159            IF( LIWORK.EQ. - 1 ) THEN
160                IDUM1( 5 ) = - 1
161            ELSE
162                IDUM1( 5 ) = 1
163            END IF
164            IDUM2( 5 ) = 23
165            CALL PCHK1MAT( N , 4 , N , 4 , IA , JA , DESCA , 9 , 0 , IDUM1 , IDUM2 ,
166       $    INFO )
167            CALL PCHK2MAT( N , 4 , NRHS , 5 , IB , JB , DESCB , 13 , N , 4 , NRHS , 5 ,
168       $    IX , JX , DESCX , 17 , 5 , IDUM1 , IDUM2 , INFO )
169        END IF
170        IF( INFO.NE.0 ) THEN
171            CALL PXERBLA( ICTXT , 'PDTRRFS' , - INFO )
172            RETURN
173        ELSE IF( LQUERY ) THEN
174            RETURN
175        END IF
176  
177        JJFBE = JJXB
178        MYRHS = NUMROC( JB + NRHS - 1 , DESCB( NB_ ) , MYCOL , DESCB( CSRC_ ) ,
179       $NPCOL )
180  
181  *     Quick return if possible
182  
183        IF( N.LE.1 .OR. NRHS.EQ.0 ) THEN
184            DO 10 JJ = JJFBE , MYRHS
185                FERR( JJ ) = ZERO
186                BERR( JJ ) = ZERO
187     10     CONTINUE
188            RETURN
189        END IF
190  
191        IF( NOTRAN ) THEN
192            TRANST = 'T'
193        ELSE
194            TRANST = 'N'
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 = maximum number of nonzero entries in each row of A , plus 1
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 90 K = 0 , JBRHS - 1
228  
229  *         Compute residual R = B - op(A) * X ,
230  *         where op(A) = A or A' , depending on TRANS.
231  
232            CALL PDCOPY( N , X , IX , JX + K , DESCX , 1 , WORK( IPR ) , IW , JW ,
233       $    DESCW , 1 )
234            CALL PDTRMV( UPLO , TRANS , DIAG , N , A , IA , JA , DESCA ,
235       $    WORK( IPR ) , IW , JW , DESCW , 1 )
236            CALL PDAXPY( N , - ONE , B , IB , JB + K , DESCB , 1 , WORK( IPR ) , IW ,
237       $    JW , DESCW , 1 )
238  
239  *         Compute componentwise relative backward error from formula
240  
241  *         max(i)( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
242  
243  *         where abs(Z) is the componentwise absolute value of the matrix
244  *         or vector Z. If the i - th component of the denominator is less
245  *         than SAFE2 , then SAFE1 is added to the i - th components of the
246  *         numerator and denominator before dividing.
247  
248            IF( MYCOL.EQ.IXBCOL ) THEN
249                IF( NP.GT.0 ) THEN
250                    DO 20 II = IIXB , IIXB + NP - 1
251                        WORK( IIW + II - IIXB ) = ABS( B( II + IOFFXB ) )
252     20             CONTINUE
253                END IF
254            END IF
255  
256            CALL PDATRMV( UPLO , TRANS , DIAG , N , ONE , A , IA , JA , DESCA , X ,
257       $    IX , JX + K , DESCX , 1 , ONE , WORK( IPB ) , IW , JW ,
258       $    DESCW , 1 )
259  
260            S = ZERO
261            IF( MYCOL.EQ.IXBCOL ) THEN
262                IF( NP.GT.0 ) THEN
263                    DO 30 II = IIW - 1 , IIW + NP - 2
264                        IF( WORK( IPB + II ).GT.SAFE2 ) THEN
265                            S = MAX( S , ABS( WORK( IPR + II ) ) /
266       $                    WORK( IPB + II ) )
267                        ELSE
268                            S = MAX( S ,( ABS( WORK( IPR + II ) ) + SAFE1 ) /
269       $( WORK( IPB + II ) + SAFE1 ) )
270                        END IF
271     30             CONTINUE
272                END IF
273            END IF
274  
275            CALL DGAMX2D( ICTXT , 'All' , ' ' , 1 , 1 , S , 1 , IDUM , IDUM , 1 ,
276       $    - 1 , MYCOL )
277            IF( MYCOL.EQ.IXBCOL )
278       $        BERR( JJFBE ) = S
279  
280  *             Bound error from formula
281  
282  *             norm(X - XTRUE) / norm(X) .le. FERR =
283  *             norm( abs(inv(op(A)))*
284  *             ( abs(R) + NZ*EPS*( abs(op(A))*abs(X) + abs(B) ))) / norm(X)
285  
286  *             where
287  *             norm(Z) is the magnitude of the largest component of Z
288  *             inv(op(A)) is the inverse of op(A)
289  *             abs(Z) is the componentwise absolute value of the matrix or
290  *             vector Z
291  *             NZ is the maximum number of nonzeros in any row of A , plus 1
292  *             EPS is machine epsilon
293  
294  *             The i - th component of abs(R) + NZ*EPS*(abs(op(A))*abs(X) + abs(B))
295  *             is incremented by SAFE1 if the i - th component of
296  *             abs(op(A))*abs(X) + abs(B) is less than SAFE2.
297  
298  *             Use PDLACON to estimate the infinity - norm of the matrix
299  *             inv(op(A)) * diag(W) ,
300  *             where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X) + abs(B) )))
301  
302                IF( MYCOL.EQ.IXBCOL ) THEN
303                    IF( NP.GT.0 ) THEN
304                        DO 40 II = IIW - 1 , IIW + NP - 2
305                            IF( WORK( IPB + II ).GT.SAFE2 ) THEN
306                                WORK( IPB + II ) = ABS( WORK( IPR + II ) ) +
307       $                        NZ*EPS*WORK( IPB + II )
308                            ELSE
309                                WORK( IPB + II ) = ABS( WORK( IPR + II ) ) +
310       $                        NZ*EPS*WORK( IPB + II ) + SAFE1
311                            END IF
312     40                 CONTINUE
313                    END IF
314                END IF
315  
316                KASE = 0
317     50 CONTINUE
318        IF( MYCOL.EQ.IXBCOL ) THEN
319            CALL DGEBS2D( ICTXT , 'Rowwise' , ' ' , NP , 1 , WORK( IPR ) ,
320       $    DESCW( LLD_ ) )
321        ELSE
322            CALL DGEBR2D( ICTXT , 'Rowwise' , ' ' , NP , 1 , WORK( IPR ) ,
323       $    DESCW( LLD_ ) , MYROW , IXBCOL )
324        END IF
325        DESCW( CSRC_ ) = MYCOL
326        CALL PDLACON ( N , WORK( IPV ) , IW , JW , DESCW , WORK( IPR ) ,
327       $IW , JW , DESCW , IWORK , EST , KASE )
328        DESCW( CSRC_ ) = IXBCOL
329  
330        IF( KASE.NE.0 ) THEN
331            IF( KASE.EQ.1 ) THEN
332  
333  *             Multiply by diag(W)*inv(op(A)').
334  
335                CALL PDTRSV( UPLO , TRANST , DIAG , N , A , IA , JA , DESCA ,
336       $        WORK( IPR ) , IW , JW , DESCW , 1 )
337                IF( MYCOL.EQ.IXBCOL ) THEN
338                    IF( NP.GT.0 ) THEN
339                        DO 60 II = IIW - 1 , IIW + NP - 2
340                            WORK( IPR + II ) = WORK( IPB + II )*WORK( IPR + II )
341     60                 CONTINUE
342                    END IF
343                END IF
344            ELSE
345  
346  *             Multiply by inv(op(A))*diag(W).
347  
348                IF( MYCOL.EQ.IXBCOL ) THEN
349                    IF( NP.GT.0 ) THEN
350                        DO 70 II = IIW - 1 , IIW + NP - 2
351                            WORK( IPR + II ) = WORK( IPB + II )*WORK( IPR + II )
352     70                 CONTINUE
353                    END IF
354                END IF
355                CALL PDTRSV( UPLO , TRANS , DIAG , N , A , IA , JA , DESCA ,
356       $        WORK( IPR ) , IW , JW , DESCW , 1 )
357            END IF
358            GO TO 50
359        END IF
360  
361  *     Normalize error.
362  
363        LSTRES = ZERO
364        IF( MYCOL.EQ.IXBCOL ) THEN
365            IF( NP.GT.0 ) THEN
366                DO 80 II = IIXB , IIXB + NP - 1
367                    LSTRES = MAX( LSTRES , ABS( X( IOFFXB + II ) ) )
368     80         CONTINUE
369            END IF
370            CALL DGAMX2D( ICTXT , 'Column' , ' ' , 1 , 1 , LSTRES , 1 , IDUM ,
371       $    IDUM , 1 , - 1 , MYCOL )
372            IF( LSTRES.NE.ZERO )
373       $        FERR( JJFBE ) = EST / LSTRES
374  
375                JJXB = JJXB + 1
376                JJFBE = JJFBE + 1
377                IOFFXB = IOFFXB + LDXB
378  
379            END IF
380  
381     90 CONTINUE
382  
383        ICURCOL = MOD( IXBCOL + 1 , NPCOL )
384  
385  *     Do for each right hand side
386  
387        DO 180 J = JN + 1 , JB + NRHS - 1 , DESCB( NB_ )
388            JBRHS = MIN( JB + NRHS - J , DESCB( NB_ ) )
389            DESCW( CSRC_ ) = ICURCOL
390  
391            DO 170 K = 0 , JBRHS - 1
392  
393  *             Compute residual R = B - op(A) * X ,
394  *             where op(A) = A or A' , depending on TRANS.
395  
396                CALL PDCOPY( N , X , IX , J + K , DESCX , 1 , WORK( IPR ) , IW , JW ,
397       $        DESCW , 1 )
398                CALL PDTRMV( UPLO , TRANS , DIAG , N , A , IA , JA , DESCA ,
399       $        WORK( IPR ) , IW , JW , DESCW , 1 )
400                CALL PDAXPY( N , - ONE , B , IB , J + K , DESCB , 1 , WORK( IPR ) ,
401       $        IW , JW , DESCW , 1 )
402  
403  *             Compute componentwise relative backward error from formula
404  
405  *             max(i)( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
406  
407  *             where abs(Z) is the componentwise absolute value of the
408  *             matrix or vector Z. If the i - th component of the
409  *             denominator is less than SAFE2 , then SAFE1 is added to the
410  *             i - th components of the numerator and denominator before
411  *             dividing.
412  
413                IF( MYCOL.EQ.IXBCOL ) THEN
414                    IF( NP.GT.0 ) THEN
415                        DO 100 II = IIXB , IIXB + NP - 1
416                            WORK( IIW + II - IIXB ) = ABS( B( II + IOFFXB ) )
417    100                 CONTINUE
418                    END IF
419                END IF
420  
421                CALL PDATRMV( UPLO , TRANS , DIAG , N , ONE , A , IA , JA , DESCA ,
422       $        X , IX , J + K , DESCX , 1 , ONE , WORK( IPB ) , IW ,
423       $        JW , DESCW , 1 )
424  
425                S = ZERO
426                IF( MYCOL.EQ.IXBCOL ) THEN
427                    IF( NP.GT.0 ) THEN
428                        DO 110 II = IIW - 1 , IIW + NP - 2
429                            IF( WORK( IPB + II ).GT.SAFE2 ) THEN
430                                S = MAX( S , ABS( WORK( IPR + II ) ) /
431       $                        WORK( IPB + II ) )
432                            ELSE
433                                S = MAX( S ,( ABS( WORK( IPR + II ) ) + SAFE1 ) /
434       $( WORK( IPB + II ) + SAFE1 ) )
435                            END IF
436    110                 CONTINUE
437                    END IF
438                END IF
439  
440                CALL DGAMX2D( ICTXT , 'All' , ' ' , 1 , 1 , S , 1 , IDUM , IDUM , 1 ,
441       $        - 1 , MYCOL )
442                IF( MYCOL.EQ.IXBCOL )
443       $            BERR( JJFBE ) = S
444  
445  *                 Bound error from formula
446  
447  *                 norm(X - XTRUE) / norm(X) .le. FERR =
448  *                 norm( abs(inv(op(A)))*
449  *                 ( abs(R) + NZ*EPS*( abs(op(A))*abs(X) + abs(B) ))) / norm(X)
450  
451  *                 where
452  *                 norm(Z) is the magnitude of the largest component of Z
453  *                 inv(op(A)) is the inverse of op(A)
454  *                 abs(Z) is the componentwise absolute value of the matrix
455  *                 or vector Z
456  *                 NZ is the maximum number of nonzeros in any row of A ,
457  *                 plus 1
458  *                 EPS is machine epsilon
459  
460  *                 The i - th component of
461  *                 abs(R) + NZ*EPS*(abs(op(A))*abs(X) + abs(B))
462  *                 is incremented by SAFE1 if the i - th component of
463  *                 abs(op(A))*abs(X) + abs(B) is less than SAFE2.
464  
465  *                 Use PDLACON to estimate the infinity - norm of the matrix
466  *                 inv(op(A)) * diag(W) ,
467  *                 where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X) + abs(B) )))
468  
469                    IF( MYCOL.EQ.IXBCOL ) THEN
470                        IF( NP.GT.0 ) THEN
471                            DO 120 II = IIW - 1 , IIW + NP - 2
472                                IF( WORK( IPB + II ).GT.SAFE2 ) THEN
473                                    WORK( IPB + II ) = ABS( WORK( IPR + II ) ) +
474       $                            NZ*EPS*WORK( IPB + II )
475                                ELSE
476                                    WORK( IPB + II ) = ABS( WORK( IPR + II ) ) +
477       $                            NZ*EPS*WORK( IPB + II ) + SAFE1
478                                END IF
479    120                     CONTINUE
480                        END IF
481                    END IF
482  
483                    KASE = 0
484    130 CONTINUE
485        IF( MYCOL.EQ.IXBCOL ) THEN
486            CALL DGEBS2D( ICTXT , 'Rowwise' , ' ' , NP , 1 , WORK( IPR ) ,
487       $    DESCW( LLD_ ) )
488        ELSE
489            CALL DGEBR2D( ICTXT , 'Rowwise' , ' ' , NP , 1 , WORK( IPR ) ,
490       $    DESCW( LLD_ ) , MYROW , IXBCOL )
491        END IF
492        DESCW( CSRC_ ) = MYCOL
493        CALL PDLACON ( N , WORK( IPV ) , IW , JW , DESCW , WORK( IPR ) ,
494       $IW , JW , DESCW , IWORK , EST , KASE )
495        DESCW( CSRC_ ) = IXBCOL
496  
497        IF( KASE.NE.0 ) THEN
498            IF( KASE.EQ.1 ) THEN
499  
500  *             Multiply by diag(W)*inv(op(A)').
501  
502                CALL PDTRSV( UPLO , TRANST , DIAG , N , A , IA , JA , DESCA ,
503       $        WORK( IPR ) , IW , JW , DESCW , 1 )
504                IF( MYCOL.EQ.IXBCOL ) THEN
505                    IF( NP.GT.0 ) THEN
506                        DO 140 II = IIW - 1 , IIW + NP - 2
507                            WORK( IPR + II ) = WORK( IPB + II )*
508       $                    WORK( IPR + II )
509    140                 CONTINUE
510                    END IF
511                END IF
512            ELSE
513  
514  *             Multiply by inv(op(A))*diag(W).
515  
516                IF( MYCOL.EQ.IXBCOL ) THEN
517                    IF( NP.GT.0 ) THEN
518                        DO 150 II = IIW - 1 , IIW + NP - 2
519                            WORK( IPR + II ) = WORK( IPB + II )*
520       $                    WORK( IPR + II )
521    150                 CONTINUE
522                    END IF
523                END IF
524                CALL PDTRSV( UPLO , TRANS , DIAG , N , A , IA , JA , DESCA ,
525       $        WORK( IPR ) , IW , JW , DESCW , 1 )
526            END IF
527            GO TO 130
528        END IF
529  
530  *     Normalize error.
531  
532        LSTRES = ZERO
533        IF( MYCOL.EQ.IXBCOL ) THEN
534            IF( NP.GT.0 ) THEN
535                DO 160 II = IIXB , IIXB + NP - 1
536                    LSTRES = MAX( LSTRES , ABS( X( IOFFXB + II ) ) )
537    160         CONTINUE
538            END IF
539            CALL DGAMX2D( ICTXT , 'Column' , ' ' , 1 , 1 , LSTRES , 1 ,
540       $    IDUM , IDUM , 1 , - 1 , MYCOL )
541            IF( LSTRES.NE.ZERO )
542       $        FERR( JJFBE ) = EST / LSTRES
543  
544                JJXB = JJXB + 1
545                JJFBE = JJFBE + 1
546                IOFFXB = IOFFXB + LDXB
547  
548            END IF
549  
550    170 CONTINUE
551  
552        ICURCOL = MOD( ICURCOL + 1 , NPCOL )
553  
554    180 CONTINUE
555  
556        WORK( 1 ) = DBLE( LWMIN )
557        IWORK( 1 ) = LIWMIN
558  
559        RETURN
560  
561  *     End of PDTRRFS
562  
563        END