LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ sppt05()

subroutine sppt05 ( character  UPLO,
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 
)

SPPT05

Purpose:
 SPPT05 tests the error bounds from iterative refinement for the
 computed solution to a system of equations A*X = B, where A is a
 symmetric 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 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]AP
          AP is REAL array, dimension (N*(N+1)/2)
          The upper or lower triangle of the symmetric 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.
[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 154 of file sppt05.f.

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