LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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.

Definition at line 187 of file ctbt05.f.

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