LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ spot05()

subroutine spot05 ( character  UPLO,
integer  N,
integer  NRHS,
real, dimension( lda, * )  A,
integer  LDA,
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 
)

SPOT05

Purpose:
 SPOT05 tests the error bounds from iterative refinement for the
 computed solution to a system of equations A*X = B, where A is a
 symmetric n by n 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 / ( (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 upper or lower triangular part of the
          symmetric matrix A is stored.
          = 'U':  Upper triangular
          = 'L':  Lower 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]A
          A is REAL array, dimension (LDA,N)
          The symmetric matrix A.  If UPLO = 'U', the leading n by n
          upper triangular part of A contains the upper triangular part
          of the matrix A, and the strictly lower triangular part of A
          is not referenced.  If UPLO = 'L', the leading n by n lower
          triangular part of A contains the lower triangular part of
          the matrix A, and the strictly upper triangular part of A is
          not referenced.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[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 162 of file spot05.f.

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