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

◆ stbt05()

subroutine stbt05 ( character  uplo,
character  trans,
character  diag,
integer  n,
integer  kd,
integer  nrhs,
real, dimension( ldab, * )  ab,
integer  ldab,
real, dimension( ldb, * )  b,
integer  ldb,
real, dimension( ldx, * )  x,
integer  ldx,
real, dimension( ldxact, * )  xact,
integer  ldxact,
real, dimension( * )  ferr,
real, dimension( * )  berr,
real, dimension( * )  reslts 
)

STBT05

Purpose:
 STBT05 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 REAL 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 REAL 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 REAL 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 REAL 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 stbt05.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 AB( LDAB, * ), B( LDB, * ), BERR( * ),
200 $ FERR( * ), RESLTS( * ), 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* ..
215* .. External Functions ..
216 LOGICAL LSAME
217 INTEGER ISAMAX
218 REAL SLAMCH
219 EXTERNAL lsame, isamax, slamch
220* ..
221* .. Intrinsic Functions ..
222 INTRINSIC abs, max, min
223* ..
224* .. Executable Statements ..
225*
226* Quick exit if N = 0 or NRHS = 0.
227*
228 IF( n.LE.0 .OR. nrhs.LE.0 ) THEN
229 reslts( 1 ) = zero
230 reslts( 2 ) = zero
231 RETURN
232 END IF
233*
234 eps = slamch( 'Epsilon' )
235 unfl = slamch( 'Safe minimum' )
236 ovfl = one / unfl
237 upper = lsame( uplo, 'U' )
238 notran = lsame( trans, 'N' )
239 unit = lsame( diag, 'U' )
240 nz = min( kd, n-1 ) + 1
241*
242* Test 1: Compute the maximum of
243* norm(X - XACT) / ( norm(X) * FERR )
244* over all the vectors X and XACT using the infinity-norm.
245*
246 errbnd = zero
247 DO 30 j = 1, nrhs
248 imax = isamax( n, x( 1, j ), 1 )
249 xnorm = max( abs( x( imax, j ) ), unfl )
250 diff = zero
251 DO 10 i = 1, n
252 diff = max( diff, abs( x( i, j )-xact( i, j ) ) )
253 10 CONTINUE
254*
255 IF( xnorm.GT.one ) THEN
256 GO TO 20
257 ELSE IF( diff.LE.ovfl*xnorm ) THEN
258 GO TO 20
259 ELSE
260 errbnd = one / eps
261 GO TO 30
262 END IF
263*
264 20 CONTINUE
265 IF( diff / xnorm.LE.ferr( j ) ) THEN
266 errbnd = max( errbnd, ( diff / xnorm ) / ferr( j ) )
267 ELSE
268 errbnd = one / eps
269 END IF
270 30 CONTINUE
271 reslts( 1 ) = errbnd
272*
273* Test 2: Compute the maximum of BERR / ( NZ*EPS + (*) ), where
274* (*) = NZ*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
275*
276 ifu = 0
277 IF( unit )
278 $ ifu = 1
279 DO 90 k = 1, nrhs
280 DO 80 i = 1, n
281 tmp = abs( b( i, k ) )
282 IF( upper ) THEN
283 IF( .NOT.notran ) THEN
284 DO 40 j = max( i-kd, 1 ), i - ifu
285 tmp = tmp + abs( ab( kd+1-i+j, i ) )*
286 $ abs( x( j, k ) )
287 40 CONTINUE
288 IF( unit )
289 $ tmp = tmp + abs( x( i, k ) )
290 ELSE
291 IF( unit )
292 $ tmp = tmp + abs( x( i, k ) )
293 DO 50 j = i + ifu, min( i+kd, n )
294 tmp = tmp + abs( ab( kd+1+i-j, j ) )*
295 $ abs( x( j, k ) )
296 50 CONTINUE
297 END IF
298 ELSE
299 IF( notran ) THEN
300 DO 60 j = max( i-kd, 1 ), i - ifu
301 tmp = tmp + abs( ab( 1+i-j, j ) )*abs( x( j, k ) )
302 60 CONTINUE
303 IF( unit )
304 $ tmp = tmp + abs( x( i, k ) )
305 ELSE
306 IF( unit )
307 $ tmp = tmp + abs( x( i, k ) )
308 DO 70 j = i + ifu, min( i+kd, n )
309 tmp = tmp + abs( ab( 1+j-i, i ) )*abs( 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 ) / ( nz*eps+nz*unfl / max( axbi, nz*unfl ) )
320 IF( k.EQ.1 ) THEN
321 reslts( 2 ) = tmp
322 ELSE
323 reslts( 2 ) = max( reslts( 2 ), tmp )
324 END IF
325 90 CONTINUE
326*
327 RETURN
328*
329* End of STBT05
330*
isamax
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
Here is the caller graph for this function: