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