LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dtpt05.f
Go to the documentation of this file.
1*> \brief \b DTPT05
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 DTPT05( UPLO, TRANS, DIAG, N, NRHS, AP, B, LDB, X, LDX,
12* XACT, LDXACT, FERR, BERR, RESLTS )
13*
14* .. Scalar Arguments ..
15* CHARACTER DIAG, TRANS, 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*> DTPT05 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*> triangular 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 matrix A is upper or lower triangular.
50*> = 'U': Upper triangular
51*> = 'L': Lower triangular
52*> \endverbatim
53*>
54*> \param[in] TRANS
55*> \verbatim
56*> TRANS is CHARACTER*1
57*> Specifies the form of the system of equations.
58*> = 'N': A * X = B (No transpose)
59*> = 'T': A'* X = B (Transpose)
60*> = 'C': A'* X = B (Conjugate transpose = Transpose)
61*> \endverbatim
62*>
63*> \param[in] DIAG
64*> \verbatim
65*> DIAG is CHARACTER*1
66*> Specifies whether or not the matrix A is unit triangular.
67*> = 'N': Non-unit triangular
68*> = 'U': Unit triangular
69*> \endverbatim
70*>
71*> \param[in] N
72*> \verbatim
73*> N is INTEGER
74*> The number of rows of the matrices X, B, and XACT, and the
75*> order of the matrix A. N >= 0.
76*> \endverbatim
77*>
78*> \param[in] NRHS
79*> \verbatim
80*> NRHS is INTEGER
81*> The number of columns of the matrices X, B, and XACT.
82*> NRHS >= 0.
83*> \endverbatim
84*>
85*> \param[in] AP
86*> \verbatim
87*> AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
88*> The upper or lower triangular matrix A, packed columnwise in
89*> a linear array. The j-th column of A is stored in the array
90*> AP as follows:
91*> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
92*> if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
93*> If DIAG = 'U', the diagonal elements of A are not referenced
94*> and are assumed to be 1.
95*> \endverbatim
96*>
97*> \param[in] B
98*> \verbatim
99*> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
100*> The right hand side vectors for the system of linear
101*> equations.
102*> \endverbatim
103*>
104*> \param[in] LDB
105*> \verbatim
106*> LDB is INTEGER
107*> The leading dimension of the array B. LDB >= max(1,N).
108*> \endverbatim
109*>
110*> \param[in] X
111*> \verbatim
112*> X is DOUBLE PRECISION array, dimension (LDX,NRHS)
113*> The computed solution vectors. Each vector is stored as a
114*> column of the matrix X.
115*> \endverbatim
116*>
117*> \param[in] LDX
118*> \verbatim
119*> LDX is INTEGER
120*> The leading dimension of the array X. LDX >= max(1,N).
121*> \endverbatim
122*>
123*> \param[in] XACT
124*> \verbatim
125*> XACT is DOUBLE PRECISION array, dimension (LDX,NRHS)
126*> The exact solution vectors. Each vector is stored as a
127*> column of the matrix XACT.
128*> \endverbatim
129*>
130*> \param[in] LDXACT
131*> \verbatim
132*> LDXACT is INTEGER
133*> The leading dimension of the array XACT. LDXACT >= max(1,N).
134*> \endverbatim
135*>
136*> \param[in] FERR
137*> \verbatim
138*> FERR is DOUBLE PRECISION array, dimension (NRHS)
139*> The estimated forward error bounds for each solution vector
140*> X. If XTRUE is the true solution, FERR bounds the magnitude
141*> of the largest entry in (X - XTRUE) divided by the magnitude
142*> of the largest entry in X.
143*> \endverbatim
144*>
145*> \param[in] BERR
146*> \verbatim
147*> BERR is DOUBLE PRECISION array, dimension (NRHS)
148*> The componentwise relative backward error of each solution
149*> vector (i.e., the smallest relative change in any entry of A
150*> or B that makes X an exact solution).
151*> \endverbatim
152*>
153*> \param[out] RESLTS
154*> \verbatim
155*> RESLTS is DOUBLE PRECISION array, dimension (2)
156*> The maximum over the NRHS solution vectors of the ratios:
157*> RESLTS(1) = norm(X - XACT) / ( norm(X) * FERR )
158*> RESLTS(2) = BERR / ( (n+1)*EPS + (*) )
159*> \endverbatim
160*
161* Authors:
162* ========
163*
164*> \author Univ. of Tennessee
165*> \author Univ. of California Berkeley
166*> \author Univ. of Colorado Denver
167*> \author NAG Ltd.
168*
169*> \ingroup double_lin
170*
171* =====================================================================
172 SUBROUTINE dtpt05( UPLO, TRANS, DIAG, N, NRHS, AP, B, LDB, X, LDX,
173 $ XACT, LDXACT, FERR, BERR, RESLTS )
174*
175* -- LAPACK test routine --
176* -- LAPACK is a software package provided by Univ. of Tennessee, --
177* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
178*
179* .. Scalar Arguments ..
180 CHARACTER DIAG, TRANS, UPLO
181 INTEGER LDB, LDX, LDXACT, N, NRHS
182* ..
183* .. Array Arguments ..
184 DOUBLE PRECISION AP( * ), B( LDB, * ), BERR( * ), FERR( * ),
185 $ reslts( * ), x( ldx, * ), xact( ldxact, * )
186* ..
187*
188* =====================================================================
189*
190* .. Parameters ..
191 DOUBLE PRECISION ZERO, ONE
192 parameter( zero = 0.0d+0, one = 1.0d+0 )
193* ..
194* .. Local Scalars ..
195 LOGICAL NOTRAN, UNIT, UPPER
196 INTEGER I, IFU, IMAX, J, JC, K
197 DOUBLE PRECISION AXBI, DIFF, EPS, ERRBND, OVFL, TMP, UNFL, XNORM
198* ..
199* .. External Functions ..
200 LOGICAL LSAME
201 INTEGER IDAMAX
202 DOUBLE PRECISION DLAMCH
203 EXTERNAL lsame, idamax, dlamch
204* ..
205* .. Intrinsic Functions ..
206 INTRINSIC abs, max, min
207* ..
208* .. Executable Statements ..
209*
210* Quick exit if N = 0 or NRHS = 0.
211*
212 IF( n.LE.0 .OR. nrhs.LE.0 ) THEN
213 reslts( 1 ) = zero
214 reslts( 2 ) = zero
215 RETURN
216 END IF
217*
218 eps = dlamch( 'Epsilon' )
219 unfl = dlamch( 'Safe minimum' )
220 ovfl = one / unfl
221 upper = lsame( uplo, 'U' )
222 notran = lsame( trans, 'N' )
223 unit = lsame( diag, 'U' )
224*
225* Test 1: Compute the maximum of
226* norm(X - XACT) / ( norm(X) * FERR )
227* over all the vectors X and XACT using the infinity-norm.
228*
229 errbnd = zero
230 DO 30 j = 1, nrhs
231 imax = idamax( n, x( 1, j ), 1 )
232 xnorm = max( abs( x( imax, j ) ), unfl )
233 diff = zero
234 DO 10 i = 1, n
235 diff = max( diff, abs( x( i, j )-xact( i, j ) ) )
236 10 CONTINUE
237*
238 IF( xnorm.GT.one ) THEN
239 GO TO 20
240 ELSE IF( diff.LE.ovfl*xnorm ) THEN
241 GO TO 20
242 ELSE
243 errbnd = one / eps
244 GO TO 30
245 END IF
246*
247 20 CONTINUE
248 IF( diff / xnorm.LE.ferr( j ) ) THEN
249 errbnd = max( errbnd, ( diff / xnorm ) / ferr( j ) )
250 ELSE
251 errbnd = one / eps
252 END IF
253 30 CONTINUE
254 reslts( 1 ) = errbnd
255*
256* Test 2: Compute the maximum of BERR / ( (n+1)*EPS + (*) ), where
257* (*) = (n+1)*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
258*
259 ifu = 0
260 IF( unit )
261 $ ifu = 1
262 DO 90 k = 1, nrhs
263 DO 80 i = 1, n
264 tmp = abs( b( i, k ) )
265 IF( upper ) THEN
266 jc = ( ( i-1 )*i ) / 2
267 IF( .NOT.notran ) THEN
268 DO 40 j = 1, i - ifu
269 tmp = tmp + abs( ap( jc+j ) )*abs( x( j, k ) )
270 40 CONTINUE
271 IF( unit )
272 $ tmp = tmp + abs( x( i, k ) )
273 ELSE
274 jc = jc + i
275 IF( unit ) THEN
276 tmp = tmp + abs( x( i, k ) )
277 jc = jc + i
278 END IF
279 DO 50 j = i + ifu, n
280 tmp = tmp + abs( ap( jc ) )*abs( x( j, k ) )
281 jc = jc + j
282 50 CONTINUE
283 END IF
284 ELSE
285 IF( notran ) THEN
286 jc = i
287 DO 60 j = 1, i - ifu
288 tmp = tmp + abs( ap( jc ) )*abs( x( j, k ) )
289 jc = jc + n - j
290 60 CONTINUE
291 IF( unit )
292 $ tmp = tmp + abs( x( i, k ) )
293 ELSE
294 jc = ( i-1 )*( n-i ) + ( i*( i+1 ) ) / 2
295 IF( unit )
296 $ tmp = tmp + abs( x( i, k ) )
297 DO 70 j = i + ifu, n
298 tmp = tmp + abs( ap( jc+j-i ) )*abs( x( j, k ) )
299 70 CONTINUE
300 END IF
301 END IF
302 IF( i.EQ.1 ) THEN
303 axbi = tmp
304 ELSE
305 axbi = min( axbi, tmp )
306 END IF
307 80 CONTINUE
308 tmp = berr( k ) / ( ( n+1 )*eps+( n+1 )*unfl /
309 $ max( axbi, ( n+1 )*unfl ) )
310 IF( k.EQ.1 ) THEN
311 reslts( 2 ) = tmp
312 ELSE
313 reslts( 2 ) = max( reslts( 2 ), tmp )
314 END IF
315 90 CONTINUE
316*
317 RETURN
318*
319* End of DTPT05
320*
321 END
dtpt05
subroutine dtpt05(uplo, trans, diag, n, nrhs, ap, b, ldb, x, ldx, xact, ldxact, ferr, berr, reslts)
DTPT05
Definition dtpt05.f:174