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