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