LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dppt05.f
Go to the documentation of this file.
1*> \brief \b DPPT05
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8* Definition:
9* ===========
10*
11* SUBROUTINE DPPT05( UPLO, N, NRHS, AP, B, LDB, X, LDX, XACT,
12* LDXACT, FERR, BERR, RESLTS )
13*
14* .. Scalar Arguments ..
15* CHARACTER UPLO
16* INTEGER LDB, LDX, LDXACT, N, NRHS
17* ..
18* .. Array Arguments ..
19* DOUBLE PRECISION AP( * ), B( LDB, * ), BERR( * ), FERR( * ),
20* $ RESLTS( * ), X( LDX, * ), XACT( LDXACT, * )
21* ..
22*
23*
24*> \par Purpose:
25* =============
26*>
27*> \verbatim
28*>
29*> DPPT05 tests the error bounds from iterative refinement for the
30*> computed solution to a system of equations A*X = B, where A is a
31*> symmetric matrix in packed storage format.
32*>
33*> RESLTS(1) = test of the error bound
34*> = norm(X - XACT) / ( norm(X) * FERR )
35*>
36*> A large value is returned if this ratio is not less than one.
37*>
38*> RESLTS(2) = residual from the iterative refinement routine
39*> = the maximum of BERR / ( (n+1)*EPS + (*) ), where
40*> (*) = (n+1)*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
41*> \endverbatim
42*
43* Arguments:
44* ==========
45*
46*> \param[in] UPLO
47*> \verbatim
48*> UPLO is CHARACTER*1
49*> Specifies whether the upper or lower triangular part of the
50*> symmetric matrix A is stored.
51*> = 'U': Upper triangular
52*> = 'L': Lower triangular
53*> \endverbatim
54*>
55*> \param[in] N
56*> \verbatim
57*> N is INTEGER
58*> The number of rows of the matrices X, B, and XACT, and the
59*> order of the matrix A. N >= 0.
60*> \endverbatim
61*>
62*> \param[in] NRHS
63*> \verbatim
64*> NRHS is INTEGER
65*> The number of columns of the matrices X, B, and XACT.
66*> NRHS >= 0.
67*> \endverbatim
68*>
69*> \param[in] AP
70*> \verbatim
71*> AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
72*> The upper or lower triangle of the symmetric matrix A, packed
73*> columnwise in a linear array. The j-th column of A is stored
74*> in the array AP as follows:
75*> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
76*> if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
77*> \endverbatim
78*>
79*> \param[in] B
80*> \verbatim
81*> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
82*> The right hand side vectors for the system of linear
83*> equations.
84*> \endverbatim
85*>
86*> \param[in] LDB
87*> \verbatim
88*> LDB is INTEGER
89*> The leading dimension of the array B. LDB >= max(1,N).
90*> \endverbatim
91*>
92*> \param[in] X
93*> \verbatim
94*> X is DOUBLE PRECISION array, dimension (LDX,NRHS)
95*> The computed solution vectors. Each vector is stored as a
96*> column of the matrix X.
97*> \endverbatim
98*>
99*> \param[in] LDX
100*> \verbatim
101*> LDX is INTEGER
102*> The leading dimension of the array X. LDX >= max(1,N).
103*> \endverbatim
104*>
105*> \param[in] XACT
106*> \verbatim
107*> XACT is DOUBLE PRECISION array, dimension (LDX,NRHS)
108*> The exact solution vectors. Each vector is stored as a
109*> column of the matrix XACT.
110*> \endverbatim
111*>
112*> \param[in] LDXACT
113*> \verbatim
114*> LDXACT is INTEGER
115*> The leading dimension of the array XACT. LDXACT >= max(1,N).
116*> \endverbatim
117*>
118*> \param[in] FERR
119*> \verbatim
120*> FERR is DOUBLE PRECISION array, dimension (NRHS)
121*> The estimated forward error bounds for each solution vector
122*> X. If XTRUE is the true solution, FERR bounds the magnitude
123*> of the largest entry in (X - XTRUE) divided by the magnitude
124*> of the largest entry in X.
125*> \endverbatim
126*>
127*> \param[in] BERR
128*> \verbatim
129*> BERR is DOUBLE PRECISION array, dimension (NRHS)
130*> The componentwise relative backward error of each solution
131*> vector (i.e., the smallest relative change in any entry of A
132*> or B that makes X an exact solution).
133*> \endverbatim
134*>
135*> \param[out] RESLTS
136*> \verbatim
137*> RESLTS is DOUBLE PRECISION array, dimension (2)
138*> The maximum over the NRHS solution vectors of the ratios:
139*> RESLTS(1) = norm(X - XACT) / ( norm(X) * FERR )
140*> RESLTS(2) = BERR / ( (n+1)*EPS + (*) )
141*> \endverbatim
142*
143* Authors:
144* ========
145*
146*> \author Univ. of Tennessee
147*> \author Univ. of California Berkeley
148*> \author Univ. of Colorado Denver
149*> \author NAG Ltd.
150*
151*> \ingroup double_lin
152*
153* =====================================================================
154 SUBROUTINE dppt05( UPLO, N, NRHS, AP, B, LDB, X, LDX, XACT,
155 $ LDXACT, FERR, BERR, RESLTS )
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*
281 END
subroutine dppt05(uplo, n, nrhs, ap, b, ldb, x, ldx, xact, ldxact, ferr, berr, reslts)
DPPT05
Definition dppt05.f:156