LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ 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 *
integer function icamax(N, CX, INCX)
ICAMAX
Definition: icamax.f:71
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the caller graph for this function: