LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
sptrfs.f
Go to the documentation of this file.
1 *> \brief \b SPTRFS
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download SPTRFS + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sptrfs.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sptrfs.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sptrfs.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SPTRFS( 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 * REAL 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 *> SPTRFS 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 REAL 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 REAL 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 REAL array, dimension (N)
76 *> The n diagonal elements of the diagonal matrix D from the
77 *> factorization computed by SPTTRF.
78 *> \endverbatim
79 *>
80 *> \param[in] EF
81 *> \verbatim
82 *> EF is REAL array, dimension (N-1)
83 *> The (n-1) subdiagonal elements of the unit bidiagonal factor
84 *> L from the factorization computed by SPTTRF.
85 *> \endverbatim
86 *>
87 *> \param[in] B
88 *> \verbatim
89 *> B is REAL 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 REAL array, dimension (LDX,NRHS)
102 *> On entry, the solution matrix X, as computed by SPTTRS.
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 REAL 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 REAL 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 REAL 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 *> \ingroup realPTcomputational
159 *
160 * =====================================================================
161  SUBROUTINE sptrfs( N, NRHS, D, E, DF, EF, B, LDB, X, LDX, FERR,
162  $ BERR, WORK, INFO )
163 *
164 * -- LAPACK computational routine --
165 * -- LAPACK is a software package provided by Univ. of Tennessee, --
166 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
167 *
168 * .. Scalar Arguments ..
169  INTEGER INFO, LDB, LDX, N, NRHS
170 * ..
171 * .. Array Arguments ..
172  REAL B( LDB, * ), BERR( * ), D( * ), DF( * ),
173  $ e( * ), ef( * ), ferr( * ), work( * ),
174  $ x( ldx, * )
175 * ..
176 *
177 * =====================================================================
178 *
179 * .. Parameters ..
180  INTEGER ITMAX
181  parameter( itmax = 5 )
182  REAL ZERO
183  parameter( zero = 0.0e+0 )
184  REAL ONE
185  parameter( one = 1.0e+0 )
186  REAL TWO
187  parameter( two = 2.0e+0 )
188  REAL THREE
189  parameter( three = 3.0e+0 )
190 * ..
191 * .. Local Scalars ..
192  INTEGER COUNT, I, IX, J, NZ
193  REAL BI, CX, DX, EPS, EX, LSTRES, S, SAFE1, SAFE2,
194  $ safmin
195 * ..
196 * .. External Subroutines ..
197  EXTERNAL saxpy, spttrs, xerbla
198 * ..
199 * .. Intrinsic Functions ..
200  INTRINSIC abs, max
201 * ..
202 * .. External Functions ..
203  INTEGER ISAMAX
204  REAL SLAMCH
205  EXTERNAL isamax, slamch
206 * ..
207 * .. Executable Statements ..
208 *
209 * Test the input parameters.
210 *
211  info = 0
212  IF( n.LT.0 ) THEN
213  info = -1
214  ELSE IF( nrhs.LT.0 ) THEN
215  info = -2
216  ELSE IF( ldb.LT.max( 1, n ) ) THEN
217  info = -8
218  ELSE IF( ldx.LT.max( 1, n ) ) THEN
219  info = -10
220  END IF
221  IF( info.NE.0 ) THEN
222  CALL xerbla( 'SPTRFS', -info )
223  RETURN
224  END IF
225 *
226 * Quick return if possible
227 *
228  IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
229  DO 10 j = 1, nrhs
230  ferr( j ) = zero
231  berr( j ) = zero
232  10 CONTINUE
233  RETURN
234  END IF
235 *
236 * NZ = maximum number of nonzero elements in each row of A, plus 1
237 *
238  nz = 4
239  eps = slamch( 'Epsilon' )
240  safmin = slamch( 'Safe minimum' )
241  safe1 = nz*safmin
242  safe2 = safe1 / eps
243 *
244 * Do for each right hand side
245 *
246  DO 90 j = 1, nrhs
247 *
248  count = 1
249  lstres = three
250  20 CONTINUE
251 *
252 * Loop until stopping criterion is satisfied.
253 *
254 * Compute residual R = B - A * X. Also compute
255 * abs(A)*abs(x) + abs(b) for use in the backward error bound.
256 *
257  IF( n.EQ.1 ) THEN
258  bi = b( 1, j )
259  dx = d( 1 )*x( 1, j )
260  work( n+1 ) = bi - dx
261  work( 1 ) = abs( bi ) + abs( dx )
262  ELSE
263  bi = b( 1, j )
264  dx = d( 1 )*x( 1, j )
265  ex = e( 1 )*x( 2, j )
266  work( n+1 ) = bi - dx - ex
267  work( 1 ) = abs( bi ) + abs( dx ) + abs( ex )
268  DO 30 i = 2, n - 1
269  bi = b( i, j )
270  cx = e( i-1 )*x( i-1, j )
271  dx = d( i )*x( i, j )
272  ex = e( i )*x( i+1, j )
273  work( n+i ) = bi - cx - dx - ex
274  work( i ) = abs( bi ) + abs( cx ) + abs( dx ) + abs( ex )
275  30 CONTINUE
276  bi = b( n, j )
277  cx = e( n-1 )*x( n-1, j )
278  dx = d( n )*x( n, j )
279  work( n+n ) = bi - cx - dx
280  work( n ) = abs( bi ) + abs( cx ) + abs( dx )
281  END IF
282 *
283 * Compute componentwise relative backward error from formula
284 *
285 * max(i) ( abs(R(i)) / ( abs(A)*abs(X) + abs(B) )(i) )
286 *
287 * where abs(Z) is the componentwise absolute value of the matrix
288 * or vector Z. If the i-th component of the denominator is less
289 * than SAFE2, then SAFE1 is added to the i-th components of the
290 * numerator and denominator before dividing.
291 *
292  s = zero
293  DO 40 i = 1, n
294  IF( work( i ).GT.safe2 ) THEN
295  s = max( s, abs( work( n+i ) ) / work( i ) )
296  ELSE
297  s = max( s, ( abs( work( n+i ) )+safe1 ) /
298  $ ( work( i )+safe1 ) )
299  END IF
300  40 CONTINUE
301  berr( j ) = s
302 *
303 * Test stopping criterion. Continue iterating if
304 * 1) The residual BERR(J) is larger than machine epsilon, and
305 * 2) BERR(J) decreased by at least a factor of 2 during the
306 * last iteration, and
307 * 3) At most ITMAX iterations tried.
308 *
309  IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
310  $ count.LE.itmax ) THEN
311 *
312 * Update solution and try again.
313 *
314  CALL spttrs( n, 1, df, ef, work( n+1 ), n, info )
315  CALL saxpy( n, one, work( n+1 ), 1, x( 1, j ), 1 )
316  lstres = berr( j )
317  count = count + 1
318  GO TO 20
319  END IF
320 *
321 * Bound error from formula
322 *
323 * norm(X - XTRUE) / norm(X) .le. FERR =
324 * norm( abs(inv(A))*
325 * ( abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) ))) / norm(X)
326 *
327 * where
328 * norm(Z) is the magnitude of the largest component of Z
329 * inv(A) is the inverse of A
330 * abs(Z) is the componentwise absolute value of the matrix or
331 * vector Z
332 * NZ is the maximum number of nonzeros in any row of A, plus 1
333 * EPS is machine epsilon
334 *
335 * The i-th component of abs(R)+NZ*EPS*(abs(A)*abs(X)+abs(B))
336 * is incremented by SAFE1 if the i-th component of
337 * abs(A)*abs(X) + abs(B) is less than SAFE2.
338 *
339  DO 50 i = 1, n
340  IF( work( i ).GT.safe2 ) THEN
341  work( i ) = abs( work( n+i ) ) + nz*eps*work( i )
342  ELSE
343  work( i ) = abs( work( n+i ) ) + nz*eps*work( i ) + safe1
344  END IF
345  50 CONTINUE
346  ix = isamax( n, work, 1 )
347  ferr( j ) = work( ix )
348 *
349 * Estimate the norm of inv(A).
350 *
351 * Solve M(A) * x = e, where M(A) = (m(i,j)) is given by
352 *
353 * m(i,j) = abs(A(i,j)), i = j,
354 * m(i,j) = -abs(A(i,j)), i .ne. j,
355 *
356 * and e = [ 1, 1, ..., 1 ]**T. Note M(A) = M(L)*D*M(L)**T.
357 *
358 * Solve M(L) * x = e.
359 *
360  work( 1 ) = one
361  DO 60 i = 2, n
362  work( i ) = one + work( i-1 )*abs( ef( i-1 ) )
363  60 CONTINUE
364 *
365 * Solve D * M(L)**T * x = b.
366 *
367  work( n ) = work( n ) / df( n )
368  DO 70 i = n - 1, 1, -1
369  work( i ) = work( i ) / df( i ) + work( i+1 )*abs( ef( i ) )
370  70 CONTINUE
371 *
372 * Compute norm(inv(A)) = max(x(i)), 1<=i<=n.
373 *
374  ix = isamax( n, work, 1 )
375  ferr( j ) = ferr( j )*abs( work( ix ) )
376 *
377 * Normalize error.
378 *
379  lstres = zero
380  DO 80 i = 1, n
381  lstres = max( lstres, abs( x( i, j ) ) )
382  80 CONTINUE
383  IF( lstres.NE.zero )
384  $ ferr( j ) = ferr( j ) / lstres
385 *
386  90 CONTINUE
387 *
388  RETURN
389 *
390 * End of SPTRFS
391 *
392  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine sptrfs(N, NRHS, D, E, DF, EF, B, LDB, X, LDX, FERR, BERR, WORK, INFO)
SPTRFS
Definition: sptrfs.f:163
subroutine spttrs(N, NRHS, D, E, B, LDB, INFO)
SPTTRS
Definition: spttrs.f:109
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:89