LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
cgbt05.f
Go to the documentation of this file.
1 *> \brief \b CGBT05
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 CGBT05( TRANS, N, KL, KU, NRHS, AB, LDAB, B, LDB, X,
12 * LDX, XACT, LDXACT, FERR, BERR, RESLTS )
13 *
14 * .. Scalar Arguments ..
15 * CHARACTER TRANS
16 * INTEGER KL, KU, LDAB, LDB, LDX, LDXACT, N, NRHS
17 * ..
18 * .. Array Arguments ..
19 * REAL BERR( * ), FERR( * ), RESLTS( * )
20 * COMPLEX AB( LDAB, * ), B( LDB, * ), X( LDX, * ),
21 * $ XACT( LDXACT, * )
22 * ..
23 *
24 *
25 *> \par Purpose:
26 * =============
27 *>
28 *> \verbatim
29 *>
30 *> CGBT05 tests the error bounds from iterative refinement for the
31 *> computed solution to a system of equations op(A)*X = B, where A is a
32 *> general band matrix of order n with kl subdiagonals and ku
33 *> superdiagonals and op(A) = A, A**T, or A**H, depending on TRANS.
34 *>
35 *> RESLTS(1) = test of the error bound
36 *> = norm(X - XACT) / ( norm(X) * FERR )
37 *>
38 *> A large value is returned if this ratio is not less than one.
39 *>
40 *> RESLTS(2) = residual from the iterative refinement routine
41 *> = the maximum of BERR / ( NZ*EPS + (*) ), where
42 *> (*) = NZ*UNFL / (min_i (abs(op(A))*abs(X) +abs(b))_i )
43 *> and NZ = max. number of nonzeros in any row of A, plus 1
44 *> \endverbatim
45 *
46 * Arguments:
47 * ==========
48 *
49 *> \param[in] TRANS
50 *> \verbatim
51 *> TRANS is CHARACTER*1
52 *> Specifies the form of the system of equations.
53 *> = 'N': A * X = B (No transpose)
54 *> = 'T': A**T * X = B (Transpose)
55 *> = 'C': A**H * X = B (Conjugate transpose = Transpose)
56 *> \endverbatim
57 *>
58 *> \param[in] N
59 *> \verbatim
60 *> N is INTEGER
61 *> The number of rows of the matrices X, B, and XACT, and the
62 *> order of the matrix A. N >= 0.
63 *> \endverbatim
64 *>
65 *> \param[in] KL
66 *> \verbatim
67 *> KL is INTEGER
68 *> The number of subdiagonals within the band of A. KL >= 0.
69 *> \endverbatim
70 *>
71 *> \param[in] KU
72 *> \verbatim
73 *> KU is INTEGER
74 *> The number of superdiagonals within the band of A. KU >= 0.
75 *> \endverbatim
76 *>
77 *> \param[in] NRHS
78 *> \verbatim
79 *> NRHS is INTEGER
80 *> The number of columns of the matrices X, B, and XACT.
81 *> NRHS >= 0.
82 *> \endverbatim
83 *>
84 *> \param[in] AB
85 *> \verbatim
86 *> AB is COMPLEX array, dimension (LDAB,N)
87 *> The original band matrix A, stored in rows 1 to KL+KU+1.
88 *> The j-th column of A is stored in the j-th column of the
89 *> array AB as follows:
90 *> AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(n,j+kl).
91 *> \endverbatim
92 *>
93 *> \param[in] LDAB
94 *> \verbatim
95 *> LDAB is INTEGER
96 *> The leading dimension of the array AB. LDAB >= KL+KU+1.
97 *> \endverbatim
98 *>
99 *> \param[in] B
100 *> \verbatim
101 *> B is COMPLEX array, dimension (LDB,NRHS)
102 *> The right hand side vectors for the system of linear
103 *> equations.
104 *> \endverbatim
105 *>
106 *> \param[in] LDB
107 *> \verbatim
108 *> LDB is INTEGER
109 *> The leading dimension of the array B. LDB >= max(1,N).
110 *> \endverbatim
111 *>
112 *> \param[in] X
113 *> \verbatim
114 *> X is COMPLEX array, dimension (LDX,NRHS)
115 *> The computed solution vectors. Each vector is stored as a
116 *> column of the matrix X.
117 *> \endverbatim
118 *>
119 *> \param[in] LDX
120 *> \verbatim
121 *> LDX is INTEGER
122 *> The leading dimension of the array X. LDX >= max(1,N).
123 *> \endverbatim
124 *>
125 *> \param[in] XACT
126 *> \verbatim
127 *> XACT is COMPLEX array, dimension (LDX,NRHS)
128 *> The exact solution vectors. Each vector is stored as a
129 *> column of the matrix XACT.
130 *> \endverbatim
131 *>
132 *> \param[in] LDXACT
133 *> \verbatim
134 *> LDXACT is INTEGER
135 *> The leading dimension of the array XACT. LDXACT >= max(1,N).
136 *> \endverbatim
137 *>
138 *> \param[in] FERR
139 *> \verbatim
140 *> FERR is REAL array, dimension (NRHS)
141 *> The estimated forward error bounds for each solution vector
142 *> X. If XTRUE is the true solution, FERR bounds the magnitude
143 *> of the largest entry in (X - XTRUE) divided by the magnitude
144 *> of the largest entry in X.
145 *> \endverbatim
146 *>
147 *> \param[in] BERR
148 *> \verbatim
149 *> BERR is REAL array, dimension (NRHS)
150 *> The componentwise relative backward error of each solution
151 *> vector (i.e., the smallest relative change in any entry of A
152 *> or B that makes X an exact solution).
153 *> \endverbatim
154 *>
155 *> \param[out] RESLTS
156 *> \verbatim
157 *> RESLTS is REAL array, dimension (2)
158 *> The maximum over the NRHS solution vectors of the ratios:
159 *> RESLTS(1) = norm(X - XACT) / ( norm(X) * FERR )
160 *> RESLTS(2) = BERR / ( NZ*EPS + (*) )
161 *> \endverbatim
162 *
163 * Authors:
164 * ========
165 *
166 *> \author Univ. of Tennessee
167 *> \author Univ. of California Berkeley
168 *> \author Univ. of Colorado Denver
169 *> \author NAG Ltd.
170 *
171 *> \ingroup complex_lin
172 *
173 * =====================================================================
174  SUBROUTINE cgbt05( TRANS, N, KL, KU, NRHS, AB, LDAB, B, LDB, X,
175  $ LDX, XACT, LDXACT, FERR, BERR, RESLTS )
176 *
177 * -- LAPACK test routine --
178 * -- LAPACK is a software package provided by Univ. of Tennessee, --
179 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
180 *
181 * .. Scalar Arguments ..
182  CHARACTER TRANS
183  INTEGER KL, KU, LDAB, LDB, LDX, LDXACT, N, NRHS
184 * ..
185 * .. Array Arguments ..
186  REAL BERR( * ), FERR( * ), RESLTS( * )
187  COMPLEX AB( LDAB, * ), B( LDB, * ), X( LDX, * ),
188  $ xact( ldxact, * )
189 * ..
190 *
191 * =====================================================================
192 *
193 * .. Parameters ..
194  REAL ZERO, ONE
195  parameter( zero = 0.0e+0, one = 1.0e+0 )
196 * ..
197 * .. Local Scalars ..
198  LOGICAL NOTRAN
199  INTEGER I, IMAX, J, K, NZ
200  REAL AXBI, DIFF, EPS, ERRBND, OVFL, TMP, UNFL, XNORM
201  COMPLEX ZDUM
202 * ..
203 * .. External Functions ..
204  LOGICAL LSAME
205  INTEGER ICAMAX
206  REAL SLAMCH
207  EXTERNAL lsame, icamax, slamch
208 * ..
209 * .. Intrinsic Functions ..
210  INTRINSIC abs, aimag, max, min, real
211 * ..
212 * .. Statement Functions ..
213  REAL CABS1
214 * ..
215 * .. Statement Function definitions ..
216  cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( zdum ) )
217 * ..
218 * .. Executable Statements ..
219 *
220 * Quick exit if N = 0 or NRHS = 0.
221 *
222  IF( n.LE.0 .OR. nrhs.LE.0 ) THEN
223  reslts( 1 ) = zero
224  reslts( 2 ) = zero
225  RETURN
226  END IF
227 *
228  eps = slamch( 'Epsilon' )
229  unfl = slamch( 'Safe minimum' )
230  ovfl = one / unfl
231  notran = lsame( trans, 'N' )
232  nz = min( kl+ku+2, n+1 )
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 = icamax( 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 / ( NZ*EPS + (*) ), where
266 * (*) = NZ*UNFL / (min_i (abs(op(A))*abs(X) +abs(b))_i )
267 *
268  DO 70 k = 1, nrhs
269  DO 60 i = 1, n
270  tmp = cabs1( b( i, k ) )
271  IF( notran ) THEN
272  DO 40 j = max( i-kl, 1 ), min( i+ku, n )
273  tmp = tmp + cabs1( ab( ku+1+i-j, j ) )*
274  $ cabs1( x( j, k ) )
275  40 CONTINUE
276  ELSE
277  DO 50 j = max( i-ku, 1 ), min( i+kl, n )
278  tmp = tmp + cabs1( ab( ku+1+j-i, i ) )*
279  $ cabs1( x( j, k ) )
280  50 CONTINUE
281  END IF
282  IF( i.EQ.1 ) THEN
283  axbi = tmp
284  ELSE
285  axbi = min( axbi, tmp )
286  END IF
287  60 CONTINUE
288  tmp = berr( k ) / ( nz*eps+nz*unfl / max( axbi, nz*unfl ) )
289  IF( k.EQ.1 ) THEN
290  reslts( 2 ) = tmp
291  ELSE
292  reslts( 2 ) = max( reslts( 2 ), tmp )
293  END IF
294  70 CONTINUE
295 *
296  RETURN
297 *
298 * End of CGBT05
299 *
300  END
subroutine cgbt05(TRANS, N, KL, KU, NRHS, AB, LDAB, B, LDB, X, LDX, XACT, LDXACT, FERR, BERR, RESLTS)
CGBT05
Definition: cgbt05.f:176