LAPACK  3.8.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.
Date
December 2016

Definition at line 158 of file sppt05.f.

158 *
159 * -- LAPACK test routine (version 3.7.0) --
160 * -- LAPACK is a software package provided by Univ. of Tennessee, --
161 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
162 * December 2016
163 *
164 * .. Scalar Arguments ..
165  CHARACTER uplo
166  INTEGER ldb, ldx, ldxact, n, nrhs
167 * ..
168 * .. Array Arguments ..
169  REAL ap( * ), b( ldb, * ), berr( * ), ferr( * ),
170  $ reslts( * ), x( ldx, * ), xact( ldxact, * )
171 * ..
172 *
173 * =====================================================================
174 *
175 * .. Parameters ..
176  REAL zero, one
177  parameter( zero = 0.0e+0, one = 1.0e+0 )
178 * ..
179 * .. Local Scalars ..
180  LOGICAL upper
181  INTEGER i, imax, j, jc, k
182  REAL axbi, diff, eps, errbnd, ovfl, tmp, unfl, xnorm
183 * ..
184 * .. External Functions ..
185  LOGICAL lsame
186  INTEGER isamax
187  REAL slamch
188  EXTERNAL lsame, isamax, slamch
189 * ..
190 * .. Intrinsic Functions ..
191  INTRINSIC abs, max, min
192 * ..
193 * .. Executable Statements ..
194 *
195 * Quick exit if N = 0 or NRHS = 0.
196 *
197  IF( n.LE.0 .OR. nrhs.LE.0 ) THEN
198  reslts( 1 ) = zero
199  reslts( 2 ) = zero
200  RETURN
201  END IF
202 *
203  eps = slamch( 'Epsilon' )
204  unfl = slamch( 'Safe minimum' )
205  ovfl = one / unfl
206  upper = lsame( uplo, 'U' )
207 *
208 * Test 1: Compute the maximum of
209 * norm(X - XACT) / ( norm(X) * FERR )
210 * over all the vectors X and XACT using the infinity-norm.
211 *
212  errbnd = zero
213  DO 30 j = 1, nrhs
214  imax = isamax( n, x( 1, j ), 1 )
215  xnorm = max( abs( x( imax, j ) ), unfl )
216  diff = zero
217  DO 10 i = 1, n
218  diff = max( diff, abs( x( i, j )-xact( i, j ) ) )
219  10 CONTINUE
220 *
221  IF( xnorm.GT.one ) THEN
222  GO TO 20
223  ELSE IF( diff.LE.ovfl*xnorm ) THEN
224  GO TO 20
225  ELSE
226  errbnd = one / eps
227  GO TO 30
228  END IF
229 *
230  20 CONTINUE
231  IF( diff / xnorm.LE.ferr( j ) ) THEN
232  errbnd = max( errbnd, ( diff / xnorm ) / ferr( j ) )
233  ELSE
234  errbnd = one / eps
235  END IF
236  30 CONTINUE
237  reslts( 1 ) = errbnd
238 *
239 * Test 2: Compute the maximum of BERR / ( (n+1)*EPS + (*) ), where
240 * (*) = (n+1)*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
241 *
242  DO 90 k = 1, nrhs
243  DO 80 i = 1, n
244  tmp = abs( b( i, k ) )
245  IF( upper ) THEN
246  jc = ( ( i-1 )*i ) / 2
247  DO 40 j = 1, i
248  tmp = tmp + abs( ap( jc+j ) )*abs( x( j, k ) )
249  40 CONTINUE
250  jc = jc + i
251  DO 50 j = i + 1, n
252  tmp = tmp + abs( ap( jc ) )*abs( x( j, k ) )
253  jc = jc + j
254  50 CONTINUE
255  ELSE
256  jc = i
257  DO 60 j = 1, i - 1
258  tmp = tmp + abs( ap( jc ) )*abs( x( j, k ) )
259  jc = jc + n - j
260  60 CONTINUE
261  DO 70 j = i, n
262  tmp = tmp + abs( ap( jc+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 SPPT05
283 *
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:73
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
Here is the caller graph for this function: