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

Download SPTRFS + dependencies [TGZ] [ZIP] [TXT]

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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

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: