LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ ctbt05()

subroutine ctbt05 ( character  UPLO,
character  TRANS,
character  DIAG,
integer  N,
integer  KD,
integer  NRHS,
complex, dimension( ldab, * )  AB,
integer  LDAB,
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 
)

CTBT05

Purpose:
 CTBT05 tests the error bounds from iterative refinement for the
 computed solution to a system of equations A*X = B, where A is a
 triangular band 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 / ( NZ*EPS + (*) ), where
             (*) = NZ*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
             and NZ = max. number of nonzeros in any row of A, plus 1
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]KD
          KD is INTEGER
          The number of super-diagonals of the matrix A if UPLO = 'U',
          or the number of sub-diagonals if UPLO = 'L'.  KD >= 0.
[in]NRHS
          NRHS is INTEGER
          The number of columns of the matrices X, B, and XACT.
          NRHS >= 0.
[in]AB
          AB is COMPLEX array, dimension (LDAB,N)
          The upper or lower triangular band matrix A, stored in the
          first kd+1 rows of the array. The j-th column of A is stored
          in the j-th column of the array AB as follows:
          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
          If DIAG = 'U', the diagonal elements of A are not referenced
          and are assumed to be 1.
[in]LDAB
          LDAB is INTEGER
          The leading dimension of the array AB.  LDAB >= KD+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 / ( NZ*EPS + (*) )
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016

Definition at line 191 of file ctbt05.f.

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