LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dppt05()

subroutine dppt05 ( character  uplo,
integer  n,
integer  nrhs,
double precision, dimension( * )  ap,
double precision, dimension( ldb, * )  b,
integer  ldb,
double precision, dimension( ldx, * )  x,
integer  ldx,
double precision, dimension( ldxact, * )  xact,
integer  ldxact,
double precision, dimension( * )  ferr,
double precision, dimension( * )  berr,
double precision, dimension( * )  reslts 
)

DPPT05

Purpose:
 DPPT05 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 dppt05.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 DOUBLE PRECISION AP( * ), B( LDB, * ), BERR( * ), FERR( * ),
167 $ RESLTS( * ), X( LDX, * ), XACT( LDXACT, * )
168* ..
169*
170* =====================================================================
171*
172* .. Parameters ..
173 DOUBLE PRECISION ZERO, ONE
174 parameter( zero = 0.0d+0, one = 1.0d+0 )
175* ..
176* .. Local Scalars ..
177 LOGICAL UPPER
178 INTEGER I, IMAX, J, JC, K
179 DOUBLE PRECISION AXBI, DIFF, EPS, ERRBND, OVFL, TMP, UNFL, XNORM
180* ..
181* .. External Functions ..
182 LOGICAL LSAME
183 INTEGER IDAMAX
184 DOUBLE PRECISION DLAMCH
185 EXTERNAL lsame, idamax, dlamch
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 = dlamch( 'Epsilon' )
201 unfl = dlamch( '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 = idamax( 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 DPPT05
280*
idamax
integer function idamax(n, dx, incx)
IDAMAX
Definition idamax.f:71
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
Here is the caller graph for this function: