LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine ctpt05 ( character  UPLO,
character  TRANS,
character  DIAG,
integer  N,
integer  NRHS,
complex, dimension( * )  AP,
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 
)

CTPT05

Purpose: 
 CTPT05 tests the error bounds from iterative refinement for the
 computed solution to a system of equations A*X = B, where A is a
 triangular matrix in packed storage format.

 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]AP
          AP is COMPLEX array, dimension (N*(N+1)/2)
          The upper or lower triangular matrix A, packed columnwise in
          a linear array.  The j-th column of A is stored in the array
          AP as follows:
          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
          If DIAG = 'U', the diagonal elements of A are not referenced
          and are assumed to be 1.
[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 177 of file ctpt05.f.

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