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