LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ dppt05()

subroutine dppt05 ( character  UPLO,
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 
)

DPPT05

Purpose:
 DPPT05 tests the error bounds from iterative refinement for the
 computed solution to a system of equations A*X = B, where A is a
 symmetric 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
          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]AP
          AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
          The upper or lower triangle of the symmetric 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 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 154 of file dppt05.f.

156 *
157 * -- LAPACK test routine --
158 * -- LAPACK is a software package provided by Univ. of Tennessee, --
159 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
160 *
161 * .. Scalar Arguments ..
162  CHARACTER UPLO
163  INTEGER LDB, LDX, LDXACT, N, NRHS
164 * ..
165 * .. Array Arguments ..
166  DOUBLE PRECISION AP( * ), B( LDB, * ), BERR( * ), FERR( * ),
167  $ RESLTS( * ), X( LDX, * ), XACT( LDXACT, * )
168 * ..
169 *
170 * =====================================================================
171 *
172 * .. Parameters ..
173  DOUBLE PRECISION ZERO, ONE
174  parameter( zero = 0.0d+0, one = 1.0d+0 )
175 * ..
176 * .. Local Scalars ..
177  LOGICAL UPPER
178  INTEGER I, IMAX, J, JC, K
179  DOUBLE PRECISION AXBI, DIFF, EPS, ERRBND, OVFL, TMP, UNFL, XNORM
180 * ..
181 * .. External Functions ..
182  LOGICAL LSAME
183  INTEGER IDAMAX
184  DOUBLE PRECISION DLAMCH
185  EXTERNAL lsame, idamax, dlamch
186 * ..
187 * .. Intrinsic Functions ..
188  INTRINSIC abs, max, min
189 * ..
190 * .. Executable Statements ..
191 *
192 * Quick exit if N = 0 or NRHS = 0.
193 *
194  IF( n.LE.0 .OR. nrhs.LE.0 ) THEN
195  reslts( 1 ) = zero
196  reslts( 2 ) = zero
197  RETURN
198  END IF
199 *
200  eps = dlamch( 'Epsilon' )
201  unfl = dlamch( 'Safe minimum' )
202  ovfl = one / unfl
203  upper = lsame( uplo, 'U' )
204 *
205 * Test 1: Compute the maximum of
206 * norm(X - XACT) / ( norm(X) * FERR )
207 * over all the vectors X and XACT using the infinity-norm.
208 *
209  errbnd = zero
210  DO 30 j = 1, nrhs
211  imax = idamax( n, x( 1, j ), 1 )
212  xnorm = max( abs( x( imax, j ) ), unfl )
213  diff = zero
214  DO 10 i = 1, n
215  diff = max( diff, abs( x( i, j )-xact( i, j ) ) )
216  10 CONTINUE
217 *
218  IF( xnorm.GT.one ) THEN
219  GO TO 20
220  ELSE IF( diff.LE.ovfl*xnorm ) THEN
221  GO TO 20
222  ELSE
223  errbnd = one / eps
224  GO TO 30
225  END IF
226 *
227  20 CONTINUE
228  IF( diff / xnorm.LE.ferr( j ) ) THEN
229  errbnd = max( errbnd, ( diff / xnorm ) / ferr( j ) )
230  ELSE
231  errbnd = one / eps
232  END IF
233  30 CONTINUE
234  reslts( 1 ) = errbnd
235 *
236 * Test 2: Compute the maximum of BERR / ( (n+1)*EPS + (*) ), where
237 * (*) = (n+1)*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
238 *
239  DO 90 k = 1, nrhs
240  DO 80 i = 1, n
241  tmp = abs( b( i, k ) )
242  IF( upper ) THEN
243  jc = ( ( i-1 )*i ) / 2
244  DO 40 j = 1, i
245  tmp = tmp + abs( ap( jc+j ) )*abs( x( j, k ) )
246  40 CONTINUE
247  jc = jc + i
248  DO 50 j = i + 1, n
249  tmp = tmp + abs( ap( jc ) )*abs( x( j, k ) )
250  jc = jc + j
251  50 CONTINUE
252  ELSE
253  jc = i
254  DO 60 j = 1, i - 1
255  tmp = tmp + abs( ap( jc ) )*abs( x( j, k ) )
256  jc = jc + n - j
257  60 CONTINUE
258  DO 70 j = i, n
259  tmp = tmp + abs( ap( jc+j-i ) )*abs( x( j, k ) )
260  70 CONTINUE
261  END IF
262  IF( i.EQ.1 ) THEN
263  axbi = tmp
264  ELSE
265  axbi = min( axbi, tmp )
266  END IF
267  80 CONTINUE
268  tmp = berr( k ) / ( ( n+1 )*eps+( n+1 )*unfl /
269  $ max( axbi, ( n+1 )*unfl ) )
270  IF( k.EQ.1 ) THEN
271  reslts( 2 ) = tmp
272  ELSE
273  reslts( 2 ) = max( reslts( 2 ), tmp )
274  END IF
275  90 CONTINUE
276 *
277  RETURN
278 *
279 * End of DPPT05
280 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:71
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
Here is the caller graph for this function: