LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.
Date
November 2011

Definition at line 191 of file stbt05.f.

191 *
192 * -- LAPACK test routine (version 3.4.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 * November 2011
196 *
197 * .. Scalar Arguments ..
198  CHARACTER diag, trans, uplo
199  INTEGER kd, ldab, ldb, ldx, ldxact, n, nrhs
200 * ..
201 * .. Array Arguments ..
202  REAL ab( ldab, * ), b( ldb, * ), berr( * ),
203  $ ferr( * ), reslts( * ), 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 * ..
218 * .. External Functions ..
219  LOGICAL lsame
220  INTEGER isamax
221  REAL slamch
222  EXTERNAL lsame, isamax, slamch
223 * ..
224 * .. Intrinsic Functions ..
225  INTRINSIC abs, max, min
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  nz = min( kd, n-1 ) + 1
244 *
245 * Test 1: Compute the maximum of
246 * norm(X - XACT) / ( norm(X) * FERR )
247 * over all the vectors X and XACT using the infinity-norm.
248 *
249  errbnd = zero
250  DO 30 j = 1, nrhs
251  imax = isamax( n, x( 1, j ), 1 )
252  xnorm = max( abs( x( imax, j ) ), unfl )
253  diff = zero
254  DO 10 i = 1, n
255  diff = max( diff, abs( x( i, j )-xact( i, j ) ) )
256  10 CONTINUE
257 *
258  IF( xnorm.GT.one ) THEN
259  GO TO 20
260  ELSE IF( diff.LE.ovfl*xnorm ) THEN
261  GO TO 20
262  ELSE
263  errbnd = one / eps
264  GO TO 30
265  END IF
266 *
267  20 CONTINUE
268  IF( diff / xnorm.LE.ferr( j ) ) THEN
269  errbnd = max( errbnd, ( diff / xnorm ) / ferr( j ) )
270  ELSE
271  errbnd = one / eps
272  END IF
273  30 CONTINUE
274  reslts( 1 ) = errbnd
275 *
276 * Test 2: Compute the maximum of BERR / ( NZ*EPS + (*) ), where
277 * (*) = NZ*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
278 *
279  ifu = 0
280  IF( unit )
281  $ ifu = 1
282  DO 90 k = 1, nrhs
283  DO 80 i = 1, n
284  tmp = abs( b( i, k ) )
285  IF( upper ) THEN
286  IF( .NOT.notran ) THEN
287  DO 40 j = max( i-kd, 1 ), i - ifu
288  tmp = tmp + abs( ab( kd+1-i+j, i ) )*
289  $ abs( x( j, k ) )
290  40 CONTINUE
291  IF( unit )
292  $ tmp = tmp + abs( x( i, k ) )
293  ELSE
294  IF( unit )
295  $ tmp = tmp + abs( x( i, k ) )
296  DO 50 j = i + ifu, min( i+kd, n )
297  tmp = tmp + abs( ab( kd+1+i-j, j ) )*
298  $ abs( x( j, k ) )
299  50 CONTINUE
300  END IF
301  ELSE
302  IF( notran ) THEN
303  DO 60 j = max( i-kd, 1 ), i - ifu
304  tmp = tmp + abs( ab( 1+i-j, j ) )*abs( x( j, k ) )
305  60 CONTINUE
306  IF( unit )
307  $ tmp = tmp + abs( x( i, k ) )
308  ELSE
309  IF( unit )
310  $ tmp = tmp + abs( x( i, k ) )
311  DO 70 j = i + ifu, min( i+kd, n )
312  tmp = tmp + abs( ab( 1+j-i, i ) )*abs( x( j, k ) )
313  70 CONTINUE
314  END IF
315  END IF
316  IF( i.EQ.1 ) THEN
317  axbi = tmp
318  ELSE
319  axbi = min( axbi, tmp )
320  END IF
321  80 CONTINUE
322  tmp = berr( k ) / ( nz*eps+nz*unfl / max( axbi, nz*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 STBT05
333 *
isamax
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.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: