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