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