LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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