LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine cppt05 ( character  UPLO,
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 
)

CPPT05

Purpose: 
 CPPT05 tests the error bounds from iterative refinement for the
 computed solution to a system of equations A*X = B, where A is a
 Hermitian 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
          Hermitian 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 COMPLEX array, dimension (N*(N+1)/2)
          The upper or lower triangle of the Hermitian 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 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.
Date
November 2011

Definition at line 159 of file cppt05.f.

159 *
160 * -- LAPACK test routine (version 3.4.0) --
161 * -- LAPACK is a software package provided by Univ. of Tennessee, --
162 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
163 * November 2011
164 *
165 * .. Scalar Arguments ..
166  CHARACTER uplo
167  INTEGER ldb, ldx, ldxact, n, nrhs
168 * ..
169 * .. Array Arguments ..
170  REAL berr( * ), ferr( * ), reslts( * )
171  COMPLEX ap( * ), b( ldb, * ), x( ldx, * ),
172  $ xact( ldxact, * )
173 * ..
174 *
175 * =====================================================================
176 *
177 * .. Parameters ..
178  REAL zero, one
179  parameter ( zero = 0.0e+0, one = 1.0e+0 )
180 * ..
181 * .. Local Scalars ..
182  LOGICAL upper
183  INTEGER i, imax, j, jc, k
184  REAL axbi, diff, eps, errbnd, ovfl, tmp, unfl, xnorm
185  COMPLEX zdum
186 * ..
187 * .. External Functions ..
188  LOGICAL lsame
189  INTEGER icamax
190  REAL slamch
191  EXTERNAL lsame, icamax, slamch
192 * ..
193 * .. Intrinsic Functions ..
194  INTRINSIC abs, aimag, max, min, real
195 * ..
196 * .. Statement Functions ..
197  REAL cabs1
198 * ..
199 * .. Statement Function definitions ..
200  cabs1( zdum ) = abs( REAL( ZDUM ) ) + abs( aimag( zdum ) )
201 * ..
202 * .. Executable Statements ..
203 *
204 * Quick exit if N = 0 or NRHS = 0.
205 *
206  IF( n.LE.0 .OR. nrhs.LE.0 ) THEN
207  reslts( 1 ) = zero
208  reslts( 2 ) = zero
209  RETURN
210  END IF
211 *
212  eps = slamch( 'Epsilon' )
213  unfl = slamch( 'Safe minimum' )
214  ovfl = one / unfl
215  upper = lsame( uplo, 'U' )
216 *
217 * Test 1: Compute the maximum of
218 * norm(X - XACT) / ( norm(X) * FERR )
219 * over all the vectors X and XACT using the infinity-norm.
220 *
221  errbnd = zero
222  DO 30 j = 1, nrhs
223  imax = icamax( n, x( 1, j ), 1 )
224  xnorm = max( cabs1( x( imax, j ) ), unfl )
225  diff = zero
226  DO 10 i = 1, n
227  diff = max( diff, cabs1( x( i, j )-xact( i, j ) ) )
228  10 CONTINUE
229 *
230  IF( xnorm.GT.one ) THEN
231  GO TO 20
232  ELSE IF( diff.LE.ovfl*xnorm ) THEN
233  GO TO 20
234  ELSE
235  errbnd = one / eps
236  GO TO 30
237  END IF
238 *
239  20 CONTINUE
240  IF( diff / xnorm.LE.ferr( j ) ) THEN
241  errbnd = max( errbnd, ( diff / xnorm ) / ferr( j ) )
242  ELSE
243  errbnd = one / eps
244  END IF
245  30 CONTINUE
246  reslts( 1 ) = errbnd
247 *
248 * Test 2: Compute the maximum of BERR / ( (n+1)*EPS + (*) ), where
249 * (*) = (n+1)*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
250 *
251  DO 90 k = 1, nrhs
252  DO 80 i = 1, n
253  tmp = cabs1( b( i, k ) )
254  IF( upper ) THEN
255  jc = ( ( i-1 )*i ) / 2
256  DO 40 j = 1, i - 1
257  tmp = tmp + cabs1( ap( jc+j ) )*cabs1( x( j, k ) )
258  40 CONTINUE
259  tmp = tmp + abs( REAL( AP( JC+I ) ) )*cabs1( x( i, k ) )
260  jc = jc + i + i
261  DO 50 j = i + 1, n
262  tmp = tmp + cabs1( ap( jc ) )*cabs1( x( j, k ) )
263  jc = jc + j
264  50 CONTINUE
265  ELSE
266  jc = i
267  DO 60 j = 1, i - 1
268  tmp = tmp + cabs1( ap( jc ) )*cabs1( x( j, k ) )
269  jc = jc + n - j
270  60 CONTINUE
271  tmp = tmp + abs( REAL( AP( JC ) ) )*cabs1( x( i, k ) )
272  DO 70 j = i + 1, n
273  tmp = tmp + cabs1( ap( jc+j-i ) )*cabs1( x( j, k ) )
274  70 CONTINUE
275  END IF
276  IF( i.EQ.1 ) THEN
277  axbi = tmp
278  ELSE
279  axbi = min( axbi, tmp )
280  END IF
281  80 CONTINUE
282  tmp = berr( k ) / ( ( n+1 )*eps+( n+1 )*unfl /
283  $ max( axbi, ( n+1 )*unfl ) )
284  IF( k.EQ.1 ) THEN
285  reslts( 2 ) = tmp
286  ELSE
287  reslts( 2 ) = max( reslts( 2 ), tmp )
288  END IF
289  90 CONTINUE
290 *
291  RETURN
292 *
293 * End of CPPT05
294 *
icamax
integer function icamax(N, CX, INCX)
ICAMAX
Definition: icamax.f:53
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the caller graph for this function: