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

◆ ctpt05()

subroutine ctpt05 ( character  uplo,
character  trans,
character  diag,
integer  n,
integer  nrhs,
complex, dimension( * )  ap,
complex, dimension( ldb, * )  b,
integer  ldb,
complex, dimension( ldx, * )  x,
integer  ldx,
complex, dimension( ldxact, * )  xact,
integer  ldxact,
real, dimension( * )  ferr,
real, dimension( * )  berr,
real, dimension( * )  reslts 
)

CTPT05

Purpose:
 CTPT05 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 COMPLEX 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 COMPLEX 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 COMPLEX 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 COMPLEX 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 173 of file ctpt05.f.

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