LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine stpt05 ( character  UPLO,
character  TRANS,
character  DIAG,
integer  N,
integer  NRHS,
real, dimension( * )  AP,
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 
)

STPT05

Purpose: 
 STPT05 tests the error bounds from iterative refinement for the
 computed solution to a system of equations A*X = B, where A is a
 triangular matrix in packed storage format.

 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 / ( (n+1)*EPS + (*) ), where
             (*) = (n+1)*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
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]NRHS
          NRHS is INTEGER
          The number of columns of the matrices X, B, and XACT.
          NRHS >= 0.
[in]AP
          AP is REAL array, dimension (N*(N+1)/2)
          The upper or lower triangular matrix A, packed columnwise in
          a linear array.  The j-th column of A is stored in the array
          AP as follows:
          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
          If DIAG = 'U', the diagonal elements of A are not referenced
          and are assumed to be 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 / ( (n+1)*EPS + (*) )
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 176 of file stpt05.f.

176 *
177 * -- LAPACK test routine (version 3.4.0) --
178 * -- LAPACK is a software package provided by Univ. of Tennessee, --
179 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
180 * November 2011
181 *
182 * .. Scalar Arguments ..
183  CHARACTER diag, trans, uplo
184  INTEGER ldb, ldx, ldxact, n, nrhs
185 * ..
186 * .. Array Arguments ..
187  REAL ap( * ), b( ldb, * ), berr( * ), ferr( * ),
188  $ reslts( * ), x( ldx, * ), xact( ldxact, * )
189 * ..
190 *
191 * =====================================================================
192 *
193 * .. Parameters ..
194  REAL zero, one
195  parameter ( zero = 0.0e+0, one = 1.0e+0 )
196 * ..
197 * .. Local Scalars ..
198  LOGICAL notran, unit, upper
199  INTEGER i, ifu, imax, j, jc, k
200  REAL axbi, diff, eps, errbnd, ovfl, tmp, unfl, xnorm
201 * ..
202 * .. External Functions ..
203  LOGICAL lsame
204  INTEGER isamax
205  REAL slamch
206  EXTERNAL lsame, isamax, slamch
207 * ..
208 * .. Intrinsic Functions ..
209  INTRINSIC abs, max, min
210 * ..
211 * .. Executable Statements ..
212 *
213 * Quick exit if N = 0 or NRHS = 0.
214 *
215  IF( n.LE.0 .OR. nrhs.LE.0 ) THEN
216  reslts( 1 ) = zero
217  reslts( 2 ) = zero
218  RETURN
219  END IF
220 *
221  eps = slamch( 'Epsilon' )
222  unfl = slamch( 'Safe minimum' )
223  ovfl = one / unfl
224  upper = lsame( uplo, 'U' )
225  notran = lsame( trans, 'N' )
226  unit = lsame( diag, 'U' )
227 *
228 * Test 1: Compute the maximum of
229 * norm(X - XACT) / ( norm(X) * FERR )
230 * over all the vectors X and XACT using the infinity-norm.
231 *
232  errbnd = zero
233  DO 30 j = 1, nrhs
234  imax = isamax( n, x( 1, j ), 1 )
235  xnorm = max( abs( x( imax, j ) ), unfl )
236  diff = zero
237  DO 10 i = 1, n
238  diff = max( diff, abs( x( i, j )-xact( i, j ) ) )
239  10 CONTINUE
240 *
241  IF( xnorm.GT.one ) THEN
242  GO TO 20
243  ELSE IF( diff.LE.ovfl*xnorm ) THEN
244  GO TO 20
245  ELSE
246  errbnd = one / eps
247  GO TO 30
248  END IF
249 *
250  20 CONTINUE
251  IF( diff / xnorm.LE.ferr( j ) ) THEN
252  errbnd = max( errbnd, ( diff / xnorm ) / ferr( j ) )
253  ELSE
254  errbnd = one / eps
255  END IF
256  30 CONTINUE
257  reslts( 1 ) = errbnd
258 *
259 * Test 2: Compute the maximum of BERR / ( (n+1)*EPS + (*) ), where
260 * (*) = (n+1)*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
261 *
262  ifu = 0
263  IF( unit )
264  $ ifu = 1
265  DO 90 k = 1, nrhs
266  DO 80 i = 1, n
267  tmp = abs( b( i, k ) )
268  IF( upper ) THEN
269  jc = ( ( i-1 )*i ) / 2
270  IF( .NOT.notran ) THEN
271  DO 40 j = 1, i - ifu
272  tmp = tmp + abs( ap( jc+j ) )*abs( x( j, k ) )
273  40 CONTINUE
274  IF( unit )
275  $ tmp = tmp + abs( x( i, k ) )
276  ELSE
277  jc = jc + i
278  IF( unit ) THEN
279  tmp = tmp + abs( x( i, k ) )
280  jc = jc + i
281  END IF
282  DO 50 j = i + ifu, n
283  tmp = tmp + abs( ap( jc ) )*abs( x( j, k ) )
284  jc = jc + j
285  50 CONTINUE
286  END IF
287  ELSE
288  IF( notran ) THEN
289  jc = i
290  DO 60 j = 1, i - ifu
291  tmp = tmp + abs( ap( jc ) )*abs( x( j, k ) )
292  jc = jc + n - j
293  60 CONTINUE
294  IF( unit )
295  $ tmp = tmp + abs( x( i, k ) )
296  ELSE
297  jc = ( i-1 )*( n-i ) + ( i*( i+1 ) ) / 2
298  IF( unit )
299  $ tmp = tmp + abs( x( i, k ) )
300  DO 70 j = i + ifu, n
301  tmp = tmp + abs( ap( jc+j-i ) )*abs( x( j, k ) )
302  70 CONTINUE
303  END IF
304  END IF
305  IF( i.EQ.1 ) THEN
306  axbi = tmp
307  ELSE
308  axbi = min( axbi, tmp )
309  END IF
310  80 CONTINUE
311  tmp = berr( k ) / ( ( n+1 )*eps+( n+1 )*unfl /
312  $ max( axbi, ( n+1 )*unfl ) )
313  IF( k.EQ.1 ) THEN
314  reslts( 2 ) = tmp
315  ELSE
316  reslts( 2 ) = max( reslts( 2 ), tmp )
317  END IF
318  90 CONTINUE
319 *
320  RETURN
321 *
322 * End of STPT05
323 *
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: