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

◆ dpot05()

subroutine dpot05 ( character uplo,
integer n,
integer nrhs,
double precision, dimension( lda, * ) a,
integer lda,
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 )

DPOT05

Purpose:
!>
!> DPOT05 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 DOUBLE PRECISION 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 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 162 of file dpot05.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 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), BERR( * ), FERR( * ),
175 $ RESLTS( * ), X( LDX, * ), XACT( LDXACT, * )
176* ..
177*
178* =====================================================================
179*
180* .. Parameters ..
181 DOUBLE PRECISION ZERO, ONE
182 parameter( zero = 0.0d+0, one = 1.0d+0 )
183* ..
184* .. Local Scalars ..
185 LOGICAL UPPER
186 INTEGER I, IMAX, J, K
187 DOUBLE PRECISION AXBI, DIFF, EPS, ERRBND, OVFL, TMP, UNFL, XNORM
188* ..
189* .. External Functions ..
190 LOGICAL LSAME
191 INTEGER IDAMAX
192 DOUBLE PRECISION DLAMCH
193 EXTERNAL lsame, idamax, dlamch
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 = dlamch( 'Epsilon' )
209 unfl = dlamch( '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 = idamax( 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 DPOT05
283*
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: