 LAPACK  3.8.0 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 + (*) )```
Date
December 2016

Definition at line 158 of file dppt05.f.

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