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

◆ dtpt05()

subroutine dtpt05 ( character  uplo,
character  trans,
character  diag,
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 
)

DTPT05

Purpose:
 DTPT05 tests the error bounds from iterative refinement for the
 computed solution to a system of equations A*X = B, where A is a
 triangular 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 matrix A is upper or lower triangular.
          = 'U':  Upper triangular
          = 'L':  Lower triangular
[in]TRANS
          TRANS is CHARACTER*1
          Specifies the form of the system of equations.
          = 'N':  A * X = B  (No transpose)
          = 'T':  A'* X = B  (Transpose)
          = 'C':  A'* X = B  (Conjugate transpose = Transpose)
[in]DIAG
          DIAG is CHARACTER*1
          Specifies whether or not the matrix A is unit triangular.
          = 'N':  Non-unit triangular
          = 'U':  Unit 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 triangular 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.
          If DIAG = 'U', the diagonal elements of A are not referenced
          and are assumed to be 1.
[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 172 of file dtpt05.f.

174*
175* -- LAPACK test routine --
176* -- LAPACK is a software package provided by Univ. of Tennessee, --
177* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
178*
179* .. Scalar Arguments ..
180 CHARACTER DIAG, TRANS, UPLO
181 INTEGER LDB, LDX, LDXACT, N, NRHS
182* ..
183* .. Array Arguments ..
184 DOUBLE PRECISION AP( * ), B( LDB, * ), BERR( * ), FERR( * ),
185 $ RESLTS( * ), X( LDX, * ), XACT( LDXACT, * )
186* ..
187*
188* =====================================================================
189*
190* .. Parameters ..
191 DOUBLE PRECISION ZERO, ONE
192 parameter( zero = 0.0d+0, one = 1.0d+0 )
193* ..
194* .. Local Scalars ..
195 LOGICAL NOTRAN, UNIT, UPPER
196 INTEGER I, IFU, IMAX, J, JC, K
197 DOUBLE PRECISION AXBI, DIFF, EPS, ERRBND, OVFL, TMP, UNFL, XNORM
198* ..
199* .. External Functions ..
200 LOGICAL LSAME
201 INTEGER IDAMAX
202 DOUBLE PRECISION DLAMCH
203 EXTERNAL lsame, idamax, dlamch
204* ..
205* .. Intrinsic Functions ..
206 INTRINSIC abs, max, min
207* ..
208* .. Executable Statements ..
209*
210* Quick exit if N = 0 or NRHS = 0.
211*
212 IF( n.LE.0 .OR. nrhs.LE.0 ) THEN
213 reslts( 1 ) = zero
214 reslts( 2 ) = zero
215 RETURN
216 END IF
217*
218 eps = dlamch( 'Epsilon' )
219 unfl = dlamch( 'Safe minimum' )
220 ovfl = one / unfl
221 upper = lsame( uplo, 'U' )
222 notran = lsame( trans, 'N' )
223 unit = lsame( diag, 'U' )
224*
225* Test 1: Compute the maximum of
226* norm(X - XACT) / ( norm(X) * FERR )
227* over all the vectors X and XACT using the infinity-norm.
228*
229 errbnd = zero
230 DO 30 j = 1, nrhs
231 imax = idamax( n, x( 1, j ), 1 )
232 xnorm = max( abs( x( imax, j ) ), unfl )
233 diff = zero
234 DO 10 i = 1, n
235 diff = max( diff, abs( x( i, j )-xact( i, j ) ) )
236 10 CONTINUE
237*
238 IF( xnorm.GT.one ) THEN
239 GO TO 20
240 ELSE IF( diff.LE.ovfl*xnorm ) THEN
241 GO TO 20
242 ELSE
243 errbnd = one / eps
244 GO TO 30
245 END IF
246*
247 20 CONTINUE
248 IF( diff / xnorm.LE.ferr( j ) ) THEN
249 errbnd = max( errbnd, ( diff / xnorm ) / ferr( j ) )
250 ELSE
251 errbnd = one / eps
252 END IF
253 30 CONTINUE
254 reslts( 1 ) = errbnd
255*
256* Test 2: Compute the maximum of BERR / ( (n+1)*EPS + (*) ), where
257* (*) = (n+1)*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
258*
259 ifu = 0
260 IF( unit )
261 $ ifu = 1
262 DO 90 k = 1, nrhs
263 DO 80 i = 1, n
264 tmp = abs( b( i, k ) )
265 IF( upper ) THEN
266 jc = ( ( i-1 )*i ) / 2
267 IF( .NOT.notran ) THEN
268 DO 40 j = 1, i - ifu
269 tmp = tmp + abs( ap( jc+j ) )*abs( x( j, k ) )
270 40 CONTINUE
271 IF( unit )
272 $ tmp = tmp + abs( x( i, k ) )
273 ELSE
274 jc = jc + i
275 IF( unit ) THEN
276 tmp = tmp + abs( x( i, k ) )
277 jc = jc + i
278 END IF
279 DO 50 j = i + ifu, n
280 tmp = tmp + abs( ap( jc ) )*abs( x( j, k ) )
281 jc = jc + j
282 50 CONTINUE
283 END IF
284 ELSE
285 IF( notran ) THEN
286 jc = i
287 DO 60 j = 1, i - ifu
288 tmp = tmp + abs( ap( jc ) )*abs( x( j, k ) )
289 jc = jc + n - j
290 60 CONTINUE
291 IF( unit )
292 $ tmp = tmp + abs( x( i, k ) )
293 ELSE
294 jc = ( i-1 )*( n-i ) + ( i*( i+1 ) ) / 2
295 IF( unit )
296 $ tmp = tmp + abs( x( i, k ) )
297 DO 70 j = i + ifu, n
298 tmp = tmp + abs( ap( jc+j-i ) )*abs( x( j, k ) )
299 70 CONTINUE
300 END IF
301 END IF
302 IF( i.EQ.1 ) THEN
303 axbi = tmp
304 ELSE
305 axbi = min( axbi, tmp )
306 END IF
307 80 CONTINUE
308 tmp = berr( k ) / ( ( n+1 )*eps+( n+1 )*unfl /
309 $ max( axbi, ( n+1 )*unfl ) )
310 IF( k.EQ.1 ) THEN
311 reslts( 2 ) = tmp
312 ELSE
313 reslts( 2 ) = max( reslts( 2 ), tmp )
314 END IF
315 90 CONTINUE
316*
317 RETURN
318*
319* End of DTPT05
320*
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: