LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dptrfs.f
Go to the documentation of this file.
1*> \brief \b DPTRFS
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download DPTRFS + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dptrfs.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dptrfs.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dptrfs.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE DPTRFS( N, NRHS, D, E, DF, EF, B, LDB, X, LDX, FERR,
20* BERR, WORK, INFO )
21*
22* .. Scalar Arguments ..
23* INTEGER INFO, LDB, LDX, N, NRHS
24* ..
25* .. Array Arguments ..
26* DOUBLE PRECISION B( LDB, * ), BERR( * ), D( * ), DF( * ),
27* $ E( * ), EF( * ), FERR( * ), WORK( * ),
28* $ X( LDX, * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> DPTRFS improves the computed solution to a system of linear
38*> equations when the coefficient matrix is symmetric positive definite
39*> and tridiagonal, and provides error bounds and backward error
40*> estimates for the solution.
41*> \endverbatim
42*
43* Arguments:
44* ==========
45*
46*> \param[in] N
47*> \verbatim
48*> N is INTEGER
49*> The order of the matrix A. N >= 0.
50*> \endverbatim
51*>
52*> \param[in] NRHS
53*> \verbatim
54*> NRHS is INTEGER
55*> The number of right hand sides, i.e., the number of columns
56*> of the matrix B. NRHS >= 0.
57*> \endverbatim
58*>
59*> \param[in] D
60*> \verbatim
61*> D is DOUBLE PRECISION array, dimension (N)
62*> The n diagonal elements of the tridiagonal matrix A.
63*> \endverbatim
64*>
65*> \param[in] E
66*> \verbatim
67*> E is DOUBLE PRECISION array, dimension (N-1)
68*> The (n-1) subdiagonal elements of the tridiagonal matrix A.
69*> \endverbatim
70*>
71*> \param[in] DF
72*> \verbatim
73*> DF is DOUBLE PRECISION array, dimension (N)
74*> The n diagonal elements of the diagonal matrix D from the
75*> factorization computed by DPTTRF.
76*> \endverbatim
77*>
78*> \param[in] EF
79*> \verbatim
80*> EF is DOUBLE PRECISION array, dimension (N-1)
81*> The (n-1) subdiagonal elements of the unit bidiagonal factor
82*> L from the factorization computed by DPTTRF.
83*> \endverbatim
84*>
85*> \param[in] B
86*> \verbatim
87*> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
88*> The right hand side matrix B.
89*> \endverbatim
90*>
91*> \param[in] LDB
92*> \verbatim
93*> LDB is INTEGER
94*> The leading dimension of the array B. LDB >= max(1,N).
95*> \endverbatim
96*>
97*> \param[in,out] X
98*> \verbatim
99*> X is DOUBLE PRECISION array, dimension (LDX,NRHS)
100*> On entry, the solution matrix X, as computed by DPTTRS.
101*> On exit, the improved solution matrix X.
102*> \endverbatim
103*>
104*> \param[in] LDX
105*> \verbatim
106*> LDX is INTEGER
107*> The leading dimension of the array X. LDX >= max(1,N).
108*> \endverbatim
109*>
110*> \param[out] FERR
111*> \verbatim
112*> FERR is DOUBLE PRECISION array, dimension (NRHS)
113*> The forward error bound for each solution vector
114*> X(j) (the j-th column of the solution matrix X).
115*> If XTRUE is the true solution corresponding to X(j), FERR(j)
116*> is an estimated upper bound for the magnitude of the largest
117*> element in (X(j) - XTRUE) divided by the magnitude of the
118*> largest element in X(j).
119*> \endverbatim
120*>
121*> \param[out] BERR
122*> \verbatim
123*> BERR is DOUBLE PRECISION array, dimension (NRHS)
124*> The componentwise relative backward error of each solution
125*> vector X(j) (i.e., the smallest relative change in
126*> any element of A or B that makes X(j) an exact solution).
127*> \endverbatim
128*>
129*> \param[out] WORK
130*> \verbatim
131*> WORK is DOUBLE PRECISION array, dimension (2*N)
132*> \endverbatim
133*>
134*> \param[out] INFO
135*> \verbatim
136*> INFO is INTEGER
137*> = 0: successful exit
138*> < 0: if INFO = -i, the i-th argument had an illegal value
139*> \endverbatim
140*
141*> \par Internal Parameters:
142* =========================
143*>
144*> \verbatim
145*> ITMAX is the maximum number of steps of iterative refinement.
146*> \endverbatim
147*
148* Authors:
149* ========
150*
151*> \author Univ. of Tennessee
152*> \author Univ. of California Berkeley
153*> \author Univ. of Colorado Denver
154*> \author NAG Ltd.
155*
156*> \ingroup ptrfs
157*
158* =====================================================================
159 SUBROUTINE dptrfs( N, NRHS, D, E, DF, EF, B, LDB, X, LDX, FERR,
160 $ BERR, WORK, INFO )
161*
162* -- LAPACK computational routine --
163* -- LAPACK is a software package provided by Univ. of Tennessee, --
164* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
165*
166* .. Scalar Arguments ..
167 INTEGER INFO, LDB, LDX, N, NRHS
168* ..
169* .. Array Arguments ..
170 DOUBLE PRECISION B( LDB, * ), BERR( * ), D( * ), DF( * ),
171 $ e( * ), ef( * ), ferr( * ), work( * ),
172 $ x( ldx, * )
173* ..
174*
175* =====================================================================
176*
177* .. Parameters ..
178 INTEGER ITMAX
179 parameter( itmax = 5 )
180 DOUBLE PRECISION ZERO
181 parameter( zero = 0.0d+0 )
182 DOUBLE PRECISION ONE
183 parameter( one = 1.0d+0 )
184 DOUBLE PRECISION TWO
185 parameter( two = 2.0d+0 )
186 DOUBLE PRECISION THREE
187 parameter( three = 3.0d+0 )
188* ..
189* .. Local Scalars ..
190 INTEGER COUNT, I, IX, J, NZ
191 DOUBLE PRECISION BI, CX, DX, EPS, EX, LSTRES, S, SAFE1, SAFE2,
192 $ safmin
193* ..
194* .. External Subroutines ..
195 EXTERNAL daxpy, dpttrs, xerbla
196* ..
197* .. Intrinsic Functions ..
198 INTRINSIC abs, max
199* ..
200* .. External Functions ..
201 INTEGER IDAMAX
202 DOUBLE PRECISION DLAMCH
203 EXTERNAL idamax, dlamch
204* ..
205* .. Executable Statements ..
206*
207* Test the input parameters.
208*
209 info = 0
210 IF( n.LT.0 ) THEN
211 info = -1
212 ELSE IF( nrhs.LT.0 ) THEN
213 info = -2
214 ELSE IF( ldb.LT.max( 1, n ) ) THEN
215 info = -8
216 ELSE IF( ldx.LT.max( 1, n ) ) THEN
217 info = -10
218 END IF
219 IF( info.NE.0 ) THEN
220 CALL xerbla( 'DPTRFS', -info )
221 RETURN
222 END IF
223*
224* Quick return if possible
225*
226 IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
227 DO 10 j = 1, nrhs
228 ferr( j ) = zero
229 berr( j ) = zero
230 10 CONTINUE
231 RETURN
232 END IF
233*
234* NZ = maximum number of nonzero elements in each row of A, plus 1
235*
236 nz = 4
237 eps = dlamch( 'Epsilon' )
238 safmin = dlamch( 'Safe minimum' )
239 safe1 = nz*safmin
240 safe2 = safe1 / eps
241*
242* Do for each right hand side
243*
244 DO 90 j = 1, nrhs
245*
246 count = 1
247 lstres = three
248 20 CONTINUE
249*
250* Loop until stopping criterion is satisfied.
251*
252* Compute residual R = B - A * X. Also compute
253* abs(A)*abs(x) + abs(b) for use in the backward error bound.
254*
255 IF( n.EQ.1 ) THEN
256 bi = b( 1, j )
257 dx = d( 1 )*x( 1, j )
258 work( n+1 ) = bi - dx
259 work( 1 ) = abs( bi ) + abs( dx )
260 ELSE
261 bi = b( 1, j )
262 dx = d( 1 )*x( 1, j )
263 ex = e( 1 )*x( 2, j )
264 work( n+1 ) = bi - dx - ex
265 work( 1 ) = abs( bi ) + abs( dx ) + abs( ex )
266 DO 30 i = 2, n - 1
267 bi = b( i, j )
268 cx = e( i-1 )*x( i-1, j )
269 dx = d( i )*x( i, j )
270 ex = e( i )*x( i+1, j )
271 work( n+i ) = bi - cx - dx - ex
272 work( i ) = abs( bi ) + abs( cx ) + abs( dx ) + abs( ex )
273 30 CONTINUE
274 bi = b( n, j )
275 cx = e( n-1 )*x( n-1, j )
276 dx = d( n )*x( n, j )
277 work( n+n ) = bi - cx - dx
278 work( n ) = abs( bi ) + abs( cx ) + abs( dx )
279 END IF
280*
281* Compute componentwise relative backward error from formula
282*
283* max(i) ( abs(R(i)) / ( abs(A)*abs(X) + abs(B) )(i) )
284*
285* where abs(Z) is the componentwise absolute value of the matrix
286* or vector Z. If the i-th component of the denominator is less
287* than SAFE2, then SAFE1 is added to the i-th components of the
288* numerator and denominator before dividing.
289*
290 s = zero
291 DO 40 i = 1, n
292 IF( work( i ).GT.safe2 ) THEN
293 s = max( s, abs( work( n+i ) ) / work( i ) )
294 ELSE
295 s = max( s, ( abs( work( n+i ) )+safe1 ) /
296 $ ( work( i )+safe1 ) )
297 END IF
298 40 CONTINUE
299 berr( j ) = s
300*
301* Test stopping criterion. Continue iterating if
302* 1) The residual BERR(J) is larger than machine epsilon, and
303* 2) BERR(J) decreased by at least a factor of 2 during the
304* last iteration, and
305* 3) At most ITMAX iterations tried.
306*
307 IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
308 $ count.LE.itmax ) THEN
309*
310* Update solution and try again.
311*
312 CALL dpttrs( n, 1, df, ef, work( n+1 ), n, info )
313 CALL daxpy( n, one, work( n+1 ), 1, x( 1, j ), 1 )
314 lstres = berr( j )
315 count = count + 1
316 GO TO 20
317 END IF
318*
319* Bound error from formula
320*
321* norm(X - XTRUE) / norm(X) .le. FERR =
322* norm( abs(inv(A))*
323* ( abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) ))) / norm(X)
324*
325* where
326* norm(Z) is the magnitude of the largest component of Z
327* inv(A) is the inverse of A
328* abs(Z) is the componentwise absolute value of the matrix or
329* vector Z
330* NZ is the maximum number of nonzeros in any row of A, plus 1
331* EPS is machine epsilon
332*
333* The i-th component of abs(R)+NZ*EPS*(abs(A)*abs(X)+abs(B))
334* is incremented by SAFE1 if the i-th component of
335* abs(A)*abs(X) + abs(B) is less than SAFE2.
336*
337 DO 50 i = 1, n
338 IF( work( i ).GT.safe2 ) THEN
339 work( i ) = abs( work( n+i ) ) + nz*eps*work( i )
340 ELSE
341 work( i ) = abs( work( n+i ) ) + nz*eps*work( i ) + safe1
342 END IF
343 50 CONTINUE
344 ix = idamax( n, work, 1 )
345 ferr( j ) = work( ix )
346*
347* Estimate the norm of inv(A).
348*
349* Solve M(A) * x = e, where M(A) = (m(i,j)) is given by
350*
351* m(i,j) = abs(A(i,j)), i = j,
352* m(i,j) = -abs(A(i,j)), i .ne. j,
353*
354* and e = [ 1, 1, ..., 1 ]**T. Note M(A) = M(L)*D*M(L)**T.
355*
356* Solve M(L) * x = e.
357*
358 work( 1 ) = one
359 DO 60 i = 2, n
360 work( i ) = one + work( i-1 )*abs( ef( i-1 ) )
361 60 CONTINUE
362*
363* Solve D * M(L)**T * x = b.
364*
365 work( n ) = work( n ) / df( n )
366 DO 70 i = n - 1, 1, -1
367 work( i ) = work( i ) / df( i ) + work( i+1 )*abs( ef( i ) )
368 70 CONTINUE
369*
370* Compute norm(inv(A)) = max(x(i)), 1<=i<=n.
371*
372 ix = idamax( n, work, 1 )
373 ferr( j ) = ferr( j )*abs( work( ix ) )
374*
375* Normalize error.
376*
377 lstres = zero
378 DO 80 i = 1, n
379 lstres = max( lstres, abs( x( i, j ) ) )
380 80 CONTINUE
381 IF( lstres.NE.zero )
382 $ ferr( j ) = ferr( j ) / lstres
383*
384 90 CONTINUE
385*
386 RETURN
387*
388* End of DPTRFS
389*
390 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
dptrfs
subroutine dptrfs(n, nrhs, d, e, df, ef, b, ldb, x, ldx, ferr, berr, work, info)
DPTRFS
Definition dptrfs.f:161
dpttrs
subroutine dpttrs(n, nrhs, d, e, b, ldb, info)
DPTTRS
Definition dpttrs.f:107