LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine ctrt05 ( character  UPLO,
character  TRANS,
character  DIAG,
integer  N,
integer  NRHS,
complex, dimension( lda, * )  A,
integer  LDA,
complex, dimension( ldb, * )  B,
integer  LDB,
complex, dimension( ldx, * )  X,
integer  LDX,
complex, dimension( ldxact, * )  XACT,
integer  LDXACT,
real, dimension( * )  FERR,
real, dimension( * )  BERR,
real, dimension( * )  RESLTS 
)

CTRT05

Purpose: 
 CTRT05 tests the error bounds from iterative refinement for the
 computed solution to a system of equations A*X = B, where A is a
 triangular n by n matrix.

 RESLTS(1) = test of the error bound
           = norm(X - XACT) / ( norm(X) * FERR )

 A large value is returned if this ratio is not less than one.

 RESLTS(2) = residual from the iterative refinement routine
           = the maximum of BERR / ( (n+1)*EPS + (*) ), where
             (*) = (n+1)*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          Specifies whether the matrix A is upper or lower triangular.
          = 'U':  Upper triangular
          = 'L':  Lower triangular
[in]TRANS
          TRANS is CHARACTER*1
          Specifies the form of the system of equations.
          = 'N':  A * X = B  (No transpose)
          = 'T':  A'* X = B  (Transpose)
          = 'C':  A'* X = B  (Conjugate transpose = Transpose)
[in]DIAG
          DIAG is CHARACTER*1
          Specifies whether or not the matrix A is unit triangular.
          = 'N':  Non-unit triangular
          = 'U':  Unit triangular
[in]N
          N is INTEGER
          The number of rows of the matrices X, B, and XACT, and the
          order of the matrix A.  N >= 0.
[in]NRHS
          NRHS is INTEGER
          The number of columns of the matrices X, B, and XACT.
          NRHS >= 0.
[in]A
          A is COMPLEX array, dimension (LDA,N)
          The triangular matrix A.  If UPLO = 'U', the leading n by n
          upper triangular part of the array A contains the upper
          triangular matrix, and the strictly lower triangular part of
          A is not referenced.  If UPLO = 'L', the leading n by n lower
          triangular part of the array A contains the lower triangular
          matrix, and the strictly upper triangular part of A is not
          referenced.  If DIAG = 'U', the diagonal elements of A are
          also not referenced and are assumed to be 1.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[in]B
          B is COMPLEX array, dimension (LDB,NRHS)
          The right hand side vectors for the system of linear
          equations.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B.  LDB >= max(1,N).
[in]X
          X is COMPLEX array, dimension (LDX,NRHS)
          The computed solution vectors.  Each vector is stored as a
          column of the matrix X.
[in]LDX
          LDX is INTEGER
          The leading dimension of the array X.  LDX >= max(1,N).
[in]XACT
          XACT is COMPLEX array, dimension (LDX,NRHS)
          The exact solution vectors.  Each vector is stored as a
          column of the matrix XACT.
[in]LDXACT
          LDXACT is INTEGER
          The leading dimension of the array XACT.  LDXACT >= max(1,N).
[in]FERR
          FERR is REAL array, dimension (NRHS)
          The estimated forward error bounds for each solution vector
          X.  If XTRUE is the true solution, FERR bounds the magnitude
          of the largest entry in (X - XTRUE) divided by the magnitude
          of the largest entry in X.
[in]BERR
          BERR is REAL array, dimension (NRHS)
          The componentwise relative backward error of each solution
          vector (i.e., the smallest relative change in any entry of A
          or B that makes X an exact solution).
[out]RESLTS
          RESLTS is REAL array, dimension (2)
          The maximum over the NRHS solution vectors of the ratios:
          RESLTS(1) = norm(X - XACT) / ( norm(X) * FERR )
          RESLTS(2) = BERR / ( (n+1)*EPS + (*) )
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 184 of file ctrt05.f.

184 *
185 * -- LAPACK test routine (version 3.4.0) --
186 * -- LAPACK is a software package provided by Univ. of Tennessee, --
187 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
188 * November 2011
189 *
190 * .. Scalar Arguments ..
191  CHARACTER diag, trans, uplo
192  INTEGER lda, ldb, ldx, ldxact, n, nrhs
193 * ..
194 * .. Array Arguments ..
195  REAL berr( * ), ferr( * ), reslts( * )
196  COMPLEX a( lda, * ), b( ldb, * ), x( ldx, * ),
197  $ xact( ldxact, * )
198 * ..
199 *
200 * =====================================================================
201 *
202 * .. Parameters ..
203  REAL zero, one
204  parameter ( zero = 0.0e+0, one = 1.0e+0 )
205 * ..
206 * .. Local Scalars ..
207  LOGICAL notran, unit, upper
208  INTEGER i, ifu, imax, j, k
209  REAL axbi, diff, eps, errbnd, ovfl, tmp, unfl, xnorm
210  COMPLEX zdum
211 * ..
212 * .. External Functions ..
213  LOGICAL lsame
214  INTEGER icamax
215  REAL slamch
216  EXTERNAL lsame, icamax, slamch
217 * ..
218 * .. Intrinsic Functions ..
219  INTRINSIC abs, aimag, max, min, real
220 * ..
221 * .. Statement Functions ..
222  REAL cabs1
223 * ..
224 * .. Statement Function definitions ..
225  cabs1( zdum ) = abs( REAL( ZDUM ) ) + abs( aimag( zdum ) )
226 * ..
227 * .. Executable Statements ..
228 *
229 * Quick exit if N = 0 or NRHS = 0.
230 *
231  IF( n.LE.0 .OR. nrhs.LE.0 ) THEN
232  reslts( 1 ) = zero
233  reslts( 2 ) = zero
234  RETURN
235  END IF
236 *
237  eps = slamch( 'Epsilon' )
238  unfl = slamch( 'Safe minimum' )
239  ovfl = one / unfl
240  upper = lsame( uplo, 'U' )
241  notran = lsame( trans, 'N' )
242  unit = lsame( diag, 'U' )
243 *
244 * Test 1: Compute the maximum of
245 * norm(X - XACT) / ( norm(X) * FERR )
246 * over all the vectors X and XACT using the infinity-norm.
247 *
248  errbnd = zero
249  DO 30 j = 1, nrhs
250  imax = icamax( n, x( 1, j ), 1 )
251  xnorm = max( cabs1( x( imax, j ) ), unfl )
252  diff = zero
253  DO 10 i = 1, n
254  diff = max( diff, cabs1( x( i, j )-xact( i, j ) ) )
255  10 CONTINUE
256 *
257  IF( xnorm.GT.one ) THEN
258  GO TO 20
259  ELSE IF( diff.LE.ovfl*xnorm ) THEN
260  GO TO 20
261  ELSE
262  errbnd = one / eps
263  GO TO 30
264  END IF
265 *
266  20 CONTINUE
267  IF( diff / xnorm.LE.ferr( j ) ) THEN
268  errbnd = max( errbnd, ( diff / xnorm ) / ferr( j ) )
269  ELSE
270  errbnd = one / eps
271  END IF
272  30 CONTINUE
273  reslts( 1 ) = errbnd
274 *
275 * Test 2: Compute the maximum of BERR / ( (n+1)*EPS + (*) ), where
276 * (*) = (n+1)*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
277 *
278  ifu = 0
279  IF( unit )
280  $ ifu = 1
281  DO 90 k = 1, nrhs
282  DO 80 i = 1, n
283  tmp = cabs1( b( i, k ) )
284  IF( upper ) THEN
285  IF( .NOT.notran ) THEN
286  DO 40 j = 1, i - ifu
287  tmp = tmp + cabs1( a( j, i ) )*cabs1( x( j, k ) )
288  40 CONTINUE
289  IF( unit )
290  $ tmp = tmp + cabs1( x( i, k ) )
291  ELSE
292  IF( unit )
293  $ tmp = tmp + cabs1( x( i, k ) )
294  DO 50 j = i + ifu, n
295  tmp = tmp + cabs1( a( i, j ) )*cabs1( x( j, k ) )
296  50 CONTINUE
297  END IF
298  ELSE
299  IF( notran ) THEN
300  DO 60 j = 1, i - ifu
301  tmp = tmp + cabs1( a( i, j ) )*cabs1( x( j, k ) )
302  60 CONTINUE
303  IF( unit )
304  $ tmp = tmp + cabs1( x( i, k ) )
305  ELSE
306  IF( unit )
307  $ tmp = tmp + cabs1( x( i, k ) )
308  DO 70 j = i + ifu, n
309  tmp = tmp + cabs1( a( j, i ) )*cabs1( x( j, k ) )
310  70 CONTINUE
311  END IF
312  END IF
313  IF( i.EQ.1 ) THEN
314  axbi = tmp
315  ELSE
316  axbi = min( axbi, tmp )
317  END IF
318  80 CONTINUE
319  tmp = berr( k ) / ( ( n+1 )*eps+( n+1 )*unfl /
320  $ max( axbi, ( n+1 )*unfl ) )
321  IF( k.EQ.1 ) THEN
322  reslts( 2 ) = tmp
323  ELSE
324  reslts( 2 ) = max( reslts( 2 ), tmp )
325  END IF
326  90 CONTINUE
327 *
328  RETURN
329 *
330 * End of CTRT05
331 *
icamax
integer function icamax(N, CX, INCX)
ICAMAX
Definition: icamax.f:53
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the caller graph for this function: