LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ stpt05()

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.

Definition at line 172 of file stpt05.f.

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