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