LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ 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 *
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:71
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the caller graph for this function: