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