LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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 *> \date November 2011
178 *
179 *> \ingroup complex16_lin
180 *
181 * =====================================================================
182  SUBROUTINE ztrt05( UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, X,
183  $ ldx, xact, ldxact, ferr, berr, reslts )
184 *
185 * -- LAPACK test routine (version 3.4.0) --
186 * -- LAPACK is a software package provided by Univ. of Tennessee, --
187 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
188 * November 2011
189 *
190 * .. Scalar Arguments ..
191  CHARACTER DIAG, TRANS, UPLO
192  INTEGER LDA, LDB, LDX, LDXACT, N, NRHS
193 * ..
194 * .. Array Arguments ..
195  DOUBLE PRECISION BERR( * ), FERR( * ), RESLTS( * )
196  COMPLEX*16 A( lda, * ), B( ldb, * ), X( ldx, * ),
197  $ xact( ldxact, * )
198 * ..
199 *
200 * =====================================================================
201 *
202 * .. Parameters ..
203  DOUBLE PRECISION ZERO, ONE
204  parameter ( zero = 0.0d+0, one = 1.0d+0 )
205 * ..
206 * .. Local Scalars ..
207  LOGICAL NOTRAN, UNIT, UPPER
208  INTEGER I, IFU, IMAX, J, K
209  DOUBLE PRECISION AXBI, DIFF, EPS, ERRBND, OVFL, TMP, UNFL, XNORM
210  COMPLEX*16 ZDUM
211 * ..
212 * .. External Functions ..
213  LOGICAL LSAME
214  INTEGER IZAMAX
215  DOUBLE PRECISION DLAMCH
216  EXTERNAL lsame, izamax, dlamch
217 * ..
218 * .. Intrinsic Functions ..
219  INTRINSIC abs, dble, dimag, max, min
220 * ..
221 * .. Statement Functions ..
222  DOUBLE PRECISION CABS1
223 * ..
224 * .. Statement Function definitions ..
225  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
226 * ..
227 * .. Executable Statements ..
228 *
229 * Quick exit if N = 0 or NRHS = 0.
230 *
231  IF( n.LE.0 .OR. nrhs.LE.0 ) THEN
232  reslts( 1 ) = zero
233  reslts( 2 ) = zero
234  RETURN
235  END IF
236 *
237  eps = dlamch( 'Epsilon' )
238  unfl = dlamch( 'Safe minimum' )
239  ovfl = one / unfl
240  upper = lsame( uplo, 'U' )
241  notran = lsame( trans, 'N' )
242  unit = lsame( diag, 'U' )
243 *
244 * Test 1: Compute the maximum of
245 * norm(X - XACT) / ( norm(X) * FERR )
246 * over all the vectors X and XACT using the infinity-norm.
247 *
248  errbnd = zero
249  DO 30 j = 1, nrhs
250  imax = izamax( n, x( 1, j ), 1 )
251  xnorm = max( cabs1( x( imax, j ) ), unfl )
252  diff = zero
253  DO 10 i = 1, n
254  diff = max( diff, cabs1( x( i, j )-xact( i, j ) ) )
255  10 CONTINUE
256 *
257  IF( xnorm.GT.one ) THEN
258  GO TO 20
259  ELSE IF( diff.LE.ovfl*xnorm ) THEN
260  GO TO 20
261  ELSE
262  errbnd = one / eps
263  GO TO 30
264  END IF
265 *
266  20 CONTINUE
267  IF( diff / xnorm.LE.ferr( j ) ) THEN
268  errbnd = max( errbnd, ( diff / xnorm ) / ferr( j ) )
269  ELSE
270  errbnd = one / eps
271  END IF
272  30 CONTINUE
273  reslts( 1 ) = errbnd
274 *
275 * Test 2: Compute the maximum of BERR / ( (n+1)*EPS + (*) ), where
276 * (*) = (n+1)*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
277 *
278  ifu = 0
279  IF( unit )
280  $ ifu = 1
281  DO 90 k = 1, nrhs
282  DO 80 i = 1, n
283  tmp = cabs1( b( i, k ) )
284  IF( upper ) THEN
285  IF( .NOT.notran ) THEN
286  DO 40 j = 1, i - ifu
287  tmp = tmp + cabs1( a( j, i ) )*cabs1( x( j, k ) )
288  40 CONTINUE
289  IF( unit )
290  $ tmp = tmp + cabs1( x( i, k ) )
291  ELSE
292  IF( unit )
293  $ tmp = tmp + cabs1( x( i, k ) )
294  DO 50 j = i + ifu, n
295  tmp = tmp + cabs1( a( i, j ) )*cabs1( x( j, k ) )
296  50 CONTINUE
297  END IF
298  ELSE
299  IF( notran ) THEN
300  DO 60 j = 1, i - ifu
301  tmp = tmp + cabs1( a( i, j ) )*cabs1( x( j, k ) )
302  60 CONTINUE
303  IF( unit )
304  $ tmp = tmp + cabs1( x( i, k ) )
305  ELSE
306  IF( unit )
307  $ tmp = tmp + cabs1( x( i, k ) )
308  DO 70 j = i + ifu, n
309  tmp = tmp + cabs1( a( j, i ) )*cabs1( x( j, k ) )
310  70 CONTINUE
311  END IF
312  END IF
313  IF( i.EQ.1 ) THEN
314  axbi = tmp
315  ELSE
316  axbi = min( axbi, tmp )
317  END IF
318  80 CONTINUE
319  tmp = berr( k ) / ( ( n+1 )*eps+( n+1 )*unfl /
320  $ max( axbi, ( n+1 )*unfl ) )
321  IF( k.EQ.1 ) THEN
322  reslts( 2 ) = tmp
323  ELSE
324  reslts( 2 ) = max( reslts( 2 ), tmp )
325  END IF
326  90 CONTINUE
327 *
328  RETURN
329 *
330 * End of ZTRT05
331 *
332  END
ztrt05
subroutine ztrt05(UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, X, LDX, XACT, LDXACT, FERR, BERR, RESLTS)
ZTRT05
Definition: ztrt05.f:184