 LAPACK  3.10.1 LAPACK: Linear Algebra PACKage

## ◆ sptrfs()

 subroutine sptrfs ( integer N, integer NRHS, real, dimension( * ) D, real, dimension( * ) E, real, dimension( * ) DF, real, dimension( * ) EF, real, dimension( ldb, * ) B, integer LDB, real, dimension( ldx, * ) X, integer LDX, real, dimension( * ) FERR, real, dimension( * ) BERR, real, dimension( * ) WORK, integer INFO )

SPTRFS

Purpose:
``` SPTRFS improves the computed solution to a system of linear
equations when the coefficient matrix is symmetric positive definite
and tridiagonal, and provides error bounds and backward error
estimates for the solution.```
Parameters
 [in] N ``` N is INTEGER The order of the matrix A. N >= 0.``` [in] NRHS ``` NRHS is INTEGER The number of right hand sides, i.e., the number of columns of the matrix B. NRHS >= 0.``` [in] D ``` D is REAL array, dimension (N) The n diagonal elements of the tridiagonal matrix A.``` [in] E ``` E is REAL array, dimension (N-1) The (n-1) subdiagonal elements of the tridiagonal matrix A.``` [in] DF ``` DF is REAL array, dimension (N) The n diagonal elements of the diagonal matrix D from the factorization computed by SPTTRF.``` [in] EF ``` EF is REAL array, dimension (N-1) The (n-1) subdiagonal elements of the unit bidiagonal factor L from the factorization computed by SPTTRF.``` [in] B ``` B is REAL array, dimension (LDB,NRHS) The right hand side matrix B.``` [in] LDB ``` LDB is INTEGER The leading dimension of the array B. LDB >= max(1,N).``` [in,out] X ``` X is REAL array, dimension (LDX,NRHS) On entry, the solution matrix X, as computed by SPTTRS. On exit, the improved solution matrix X.``` [in] LDX ``` LDX is INTEGER The leading dimension of the array X. LDX >= max(1,N).``` [out] FERR ``` FERR is REAL array, dimension (NRHS) The forward error bound for each solution vector X(j) (the j-th column of the solution matrix X). If XTRUE is the true solution corresponding to X(j), FERR(j) is an estimated upper bound for the magnitude of the largest element in (X(j) - XTRUE) divided by the magnitude of the largest element in X(j).``` [out] BERR ``` BERR is REAL array, dimension (NRHS) The componentwise relative backward error of each solution vector X(j) (i.e., the smallest relative change in any element of A or B that makes X(j) an exact solution).``` [out] WORK ` WORK is REAL array, dimension (2*N)` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value```
Internal Parameters:
`  ITMAX is the maximum number of steps of iterative refinement.`

Definition at line 161 of file sptrfs.f.

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 *
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:71
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
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
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: