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