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