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