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