LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ctbt05.f
Go to the documentation of this file.
1*> \brief \b CTBT05
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 CTBT05( 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* REAL BERR( * ), FERR( * ), RESLTS( * )
20* COMPLEX AB( LDAB, * ), B( LDB, * ), X( LDX, * ),
21* $ XACT( LDXACT, * )
22* ..
23*
24*
25*> \par Purpose:
26* =============
27*>
28*> \verbatim
29*>
30*> CTBT05 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 COMPLEX 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 COMPLEX 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 COMPLEX 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 COMPLEX 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 REAL 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 REAL 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 REAL 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 complex_lin
185*
186* =====================================================================
187 SUBROUTINE ctbt05( 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 REAL BERR( * ), FERR( * ), RESLTS( * )
200 COMPLEX AB( LDAB, * ), B( LDB, * ), X( LDX, * ),
201 $ xact( ldxact, * )
202* ..
203*
204* =====================================================================
205*
206* .. Parameters ..
207 REAL ZERO, ONE
208 parameter( zero = 0.0e+0, one = 1.0e+0 )
209* ..
210* .. Local Scalars ..
211 LOGICAL NOTRAN, UNIT, UPPER
212 INTEGER I, IFU, IMAX, J, K, NZ
213 REAL AXBI, DIFF, EPS, ERRBND, OVFL, TMP, UNFL, XNORM
214 COMPLEX ZDUM
215* ..
216* .. External Functions ..
217 LOGICAL LSAME
218 INTEGER ICAMAX
219 REAL SLAMCH
220 EXTERNAL lsame, icamax, slamch
221* ..
222* .. Intrinsic Functions ..
223 INTRINSIC abs, aimag, max, min, real
224* ..
225* .. Statement Functions ..
226 REAL CABS1
227* ..
228* .. Statement Function definitions ..
229 cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( zdum ) )
230* ..
231* .. Executable Statements ..
232*
233* Quick exit if N = 0 or NRHS = 0.
234*
235 IF( n.LE.0 .OR. nrhs.LE.0 ) THEN
236 reslts( 1 ) = zero
237 reslts( 2 ) = zero
238 RETURN
239 END IF
240*
241 eps = slamch( 'Epsilon' )
242 unfl = slamch( 'Safe minimum' )
243 ovfl = one / unfl
244 upper = lsame( uplo, 'U' )
245 notran = lsame( trans, 'N' )
246 unit = lsame( diag, 'U' )
247 nz = min( kd, n-1 ) + 1
248*
249* Test 1: Compute the maximum of
250* norm(X - XACT) / ( norm(X) * FERR )
251* over all the vectors X and XACT using the infinity-norm.
252*
253 errbnd = zero
254 DO 30 j = 1, nrhs
255 imax = icamax( n, x( 1, j ), 1 )
256 xnorm = max( cabs1( x( imax, j ) ), unfl )
257 diff = zero
258 DO 10 i = 1, n
259 diff = max( diff, cabs1( x( i, j )-xact( i, j ) ) )
260 10 CONTINUE
261*
262 IF( xnorm.GT.one ) THEN
263 GO TO 20
264 ELSE IF( diff.LE.ovfl*xnorm ) THEN
265 GO TO 20
266 ELSE
267 errbnd = one / eps
268 GO TO 30
269 END IF
270*
271 20 CONTINUE
272 IF( diff / xnorm.LE.ferr( j ) ) THEN
273 errbnd = max( errbnd, ( diff / xnorm ) / ferr( j ) )
274 ELSE
275 errbnd = one / eps
276 END IF
277 30 CONTINUE
278 reslts( 1 ) = errbnd
279*
280* Test 2: Compute the maximum of BERR / ( NZ*EPS + (*) ), where
281* (*) = NZ*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
282*
283 ifu = 0
284 IF( unit )
285 $ ifu = 1
286 DO 90 k = 1, nrhs
287 DO 80 i = 1, n
288 tmp = cabs1( b( i, k ) )
289 IF( upper ) THEN
290 IF( .NOT.notran ) THEN
291 DO 40 j = max( i-kd, 1 ), i - ifu
292 tmp = tmp + cabs1( ab( kd+1-i+j, i ) )*
293 $ cabs1( x( j, k ) )
294 40 CONTINUE
295 IF( unit )
296 $ tmp = tmp + cabs1( x( i, k ) )
297 ELSE
298 IF( unit )
299 $ tmp = tmp + cabs1( x( i, k ) )
300 DO 50 j = i + ifu, min( i+kd, n )
301 tmp = tmp + cabs1( ab( kd+1+i-j, j ) )*
302 $ cabs1( x( j, k ) )
303 50 CONTINUE
304 END IF
305 ELSE
306 IF( notran ) THEN
307 DO 60 j = max( i-kd, 1 ), i - ifu
308 tmp = tmp + cabs1( ab( 1+i-j, j ) )*
309 $ cabs1( x( j, k ) )
310 60 CONTINUE
311 IF( unit )
312 $ tmp = tmp + cabs1( x( i, k ) )
313 ELSE
314 IF( unit )
315 $ tmp = tmp + cabs1( x( i, k ) )
316 DO 70 j = i + ifu, min( i+kd, n )
317 tmp = tmp + cabs1( ab( 1+j-i, i ) )*
318 $ cabs1( x( j, k ) )
319 70 CONTINUE
320 END IF
321 END IF
322 IF( i.EQ.1 ) THEN
323 axbi = tmp
324 ELSE
325 axbi = min( axbi, tmp )
326 END IF
327 80 CONTINUE
328 tmp = berr( k ) / ( nz*eps+nz*unfl / max( axbi, nz*unfl ) )
329 IF( k.EQ.1 ) THEN
330 reslts( 2 ) = tmp
331 ELSE
332 reslts( 2 ) = max( reslts( 2 ), tmp )
333 END IF
334 90 CONTINUE
335*
336 RETURN
337*
338* End of CTBT05
339*
340 END
ctbt05
subroutine ctbt05(uplo, trans, diag, n, kd, nrhs, ab, ldab, b, ldb, x, ldx, xact, ldxact, ferr, berr, reslts)
CTBT05
Definition ctbt05.f:189