Routine: PZTRRFS()  File: SRC\pztrrfs.f

 
 
# lines: 806
  # code: 806
  # comment: 0
  # blank:0
# Variables:98
# Callers:0
# Callings:2
# Words:537
# Keywords:352
 

 

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