LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ sgtrfs()

subroutine sgtrfs ( character  TRANS,
integer  N,
integer  NRHS,
real, dimension( * )  DL,
real, dimension( * )  D,
real, dimension( * )  DU,
real, dimension( * )  DLF,
real, dimension( * )  DF,
real, dimension( * )  DUF,
real, dimension( * )  DU2,
integer, dimension( * )  IPIV,
real, dimension( ldb, * )  B,
integer  LDB,
real, dimension( ldx, * )  X,
integer  LDX,
real, dimension( * )  FERR,
real, dimension( * )  BERR,
real, dimension( * )  WORK,
integer, dimension( * )  IWORK,
integer  INFO 
)

SGTRFS

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

Purpose:
 SGTRFS improves the computed solution to a system of linear
 equations when the coefficient matrix is tridiagonal, and provides
 error bounds and backward error estimates for the solution.
Parameters
[in]TRANS
          TRANS is CHARACTER*1
          Specifies the form of the system of equations:
          = 'N':  A * X = B     (No transpose)
          = 'T':  A**T * X = B  (Transpose)
          = 'C':  A**H * X = B  (Conjugate transpose = Transpose)
[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]DL
          DL is REAL array, dimension (N-1)
          The (n-1) subdiagonal elements of A.
[in]D
          D is REAL array, dimension (N)
          The diagonal elements of A.
[in]DU
          DU is REAL array, dimension (N-1)
          The (n-1) superdiagonal elements of A.
[in]DLF
          DLF is REAL array, dimension (N-1)
          The (n-1) multipliers that define the matrix L from the
          LU factorization of A as computed by SGTTRF.
[in]DF
          DF is REAL array, dimension (N)
          The n diagonal elements of the upper triangular matrix U from
          the LU factorization of A.
[in]DUF
          DUF is REAL array, dimension (N-1)
          The (n-1) elements of the first superdiagonal of U.
[in]DU2
          DU2 is REAL array, dimension (N-2)
          The (n-2) elements of the second superdiagonal of U.
[in]IPIV
          IPIV is INTEGER array, dimension (N)
          The pivot indices; for 1 <= i <= n, row i of the matrix was
          interchanged with row IPIV(i).  IPIV(i) will always be either
          i or i+1; IPIV(i) = i indicates a row interchange was not
          required.
[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 SGTTRS.
          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 estimated 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).  The estimate is as reliable as
          the estimate for RCOND, and is almost always a slight
          overestimate of the true error.
[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 (3*N)
[out]IWORK
          IWORK is INTEGER array, dimension (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 206 of file sgtrfs.f.

209 *
210 * -- LAPACK computational routine --
211 * -- LAPACK is a software package provided by Univ. of Tennessee, --
212 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
213 *
214 * .. Scalar Arguments ..
215  CHARACTER TRANS
216  INTEGER INFO, LDB, LDX, N, NRHS
217 * ..
218 * .. Array Arguments ..
219  INTEGER IPIV( * ), IWORK( * )
220  REAL B( LDB, * ), BERR( * ), D( * ), DF( * ),
221  $ DL( * ), DLF( * ), DU( * ), DU2( * ), DUF( * ),
222  $ FERR( * ), WORK( * ), X( LDX, * )
223 * ..
224 *
225 * =====================================================================
226 *
227 * .. Parameters ..
228  INTEGER ITMAX
229  parameter( itmax = 5 )
230  REAL ZERO, ONE
231  parameter( zero = 0.0e+0, one = 1.0e+0 )
232  REAL TWO
233  parameter( two = 2.0e+0 )
234  REAL THREE
235  parameter( three = 3.0e+0 )
236 * ..
237 * .. Local Scalars ..
238  LOGICAL NOTRAN
239  CHARACTER TRANSN, TRANST
240  INTEGER COUNT, I, J, KASE, NZ
241  REAL EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN
242 * ..
243 * .. Local Arrays ..
244  INTEGER ISAVE( 3 )
245 * ..
246 * .. External Subroutines ..
247  EXTERNAL saxpy, scopy, sgttrs, slacn2, slagtm, xerbla
248 * ..
249 * .. Intrinsic Functions ..
250  INTRINSIC abs, max
251 * ..
252 * .. External Functions ..
253  LOGICAL LSAME
254  REAL SLAMCH
255  EXTERNAL lsame, slamch
256 * ..
257 * .. Executable Statements ..
258 *
259 * Test the input parameters.
260 *
261  info = 0
262  notran = lsame( trans, 'N' )
263  IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
264  $ lsame( trans, 'C' ) ) THEN
265  info = -1
266  ELSE IF( n.LT.0 ) THEN
267  info = -2
268  ELSE IF( nrhs.LT.0 ) THEN
269  info = -3
270  ELSE IF( ldb.LT.max( 1, n ) ) THEN
271  info = -13
272  ELSE IF( ldx.LT.max( 1, n ) ) THEN
273  info = -15
274  END IF
275  IF( info.NE.0 ) THEN
276  CALL xerbla( 'SGTRFS', -info )
277  RETURN
278  END IF
279 *
280 * Quick return if possible
281 *
282  IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
283  DO 10 j = 1, nrhs
284  ferr( j ) = zero
285  berr( j ) = zero
286  10 CONTINUE
287  RETURN
288  END IF
289 *
290  IF( notran ) THEN
291  transn = 'N'
292  transt = 'T'
293  ELSE
294  transn = 'T'
295  transt = 'N'
296  END IF
297 *
298 * NZ = maximum number of nonzero elements in each row of A, plus 1
299 *
300  nz = 4
301  eps = slamch( 'Epsilon' )
302  safmin = slamch( 'Safe minimum' )
303  safe1 = nz*safmin
304  safe2 = safe1 / eps
305 *
306 * Do for each right hand side
307 *
308  DO 110 j = 1, nrhs
309 *
310  count = 1
311  lstres = three
312  20 CONTINUE
313 *
314 * Loop until stopping criterion is satisfied.
315 *
316 * Compute residual R = B - op(A) * X,
317 * where op(A) = A, A**T, or A**H, depending on TRANS.
318 *
319  CALL scopy( n, b( 1, j ), 1, work( n+1 ), 1 )
320  CALL slagtm( trans, n, 1, -one, dl, d, du, x( 1, j ), ldx, one,
321  $ work( n+1 ), n )
322 *
323 * Compute abs(op(A))*abs(x) + abs(b) for use in the backward
324 * error bound.
325 *
326  IF( notran ) THEN
327  IF( n.EQ.1 ) THEN
328  work( 1 ) = abs( b( 1, j ) ) + abs( d( 1 )*x( 1, j ) )
329  ELSE
330  work( 1 ) = abs( b( 1, j ) ) + abs( d( 1 )*x( 1, j ) ) +
331  $ abs( du( 1 )*x( 2, j ) )
332  DO 30 i = 2, n - 1
333  work( i ) = abs( b( i, j ) ) +
334  $ abs( dl( i-1 )*x( i-1, j ) ) +
335  $ abs( d( i )*x( i, j ) ) +
336  $ abs( du( i )*x( i+1, j ) )
337  30 CONTINUE
338  work( n ) = abs( b( n, j ) ) +
339  $ abs( dl( n-1 )*x( n-1, j ) ) +
340  $ abs( d( n )*x( n, j ) )
341  END IF
342  ELSE
343  IF( n.EQ.1 ) THEN
344  work( 1 ) = abs( b( 1, j ) ) + abs( d( 1 )*x( 1, j ) )
345  ELSE
346  work( 1 ) = abs( b( 1, j ) ) + abs( d( 1 )*x( 1, j ) ) +
347  $ abs( dl( 1 )*x( 2, j ) )
348  DO 40 i = 2, n - 1
349  work( i ) = abs( b( i, j ) ) +
350  $ abs( du( i-1 )*x( i-1, j ) ) +
351  $ abs( d( i )*x( i, j ) ) +
352  $ abs( dl( i )*x( i+1, j ) )
353  40 CONTINUE
354  work( n ) = abs( b( n, j ) ) +
355  $ abs( du( n-1 )*x( n-1, j ) ) +
356  $ abs( d( n )*x( n, j ) )
357  END IF
358  END IF
359 *
360 * Compute componentwise relative backward error from formula
361 *
362 * max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
363 *
364 * where abs(Z) is the componentwise absolute value of the matrix
365 * or vector Z. If the i-th component of the denominator is less
366 * than SAFE2, then SAFE1 is added to the i-th components of the
367 * numerator and denominator before dividing.
368 *
369  s = zero
370  DO 50 i = 1, n
371  IF( work( i ).GT.safe2 ) THEN
372  s = max( s, abs( work( n+i ) ) / work( i ) )
373  ELSE
374  s = max( s, ( abs( work( n+i ) )+safe1 ) /
375  $ ( work( i )+safe1 ) )
376  END IF
377  50 CONTINUE
378  berr( j ) = s
379 *
380 * Test stopping criterion. Continue iterating if
381 * 1) The residual BERR(J) is larger than machine epsilon, and
382 * 2) BERR(J) decreased by at least a factor of 2 during the
383 * last iteration, and
384 * 3) At most ITMAX iterations tried.
385 *
386  IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
387  $ count.LE.itmax ) THEN
388 *
389 * Update solution and try again.
390 *
391  CALL sgttrs( trans, n, 1, dlf, df, duf, du2, ipiv,
392  $ work( n+1 ), n, info )
393  CALL saxpy( n, one, work( n+1 ), 1, x( 1, j ), 1 )
394  lstres = berr( j )
395  count = count + 1
396  GO TO 20
397  END IF
398 *
399 * Bound error from formula
400 *
401 * norm(X - XTRUE) / norm(X) .le. FERR =
402 * norm( abs(inv(op(A)))*
403 * ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
404 *
405 * where
406 * norm(Z) is the magnitude of the largest component of Z
407 * inv(op(A)) is the inverse of op(A)
408 * abs(Z) is the componentwise absolute value of the matrix or
409 * vector Z
410 * NZ is the maximum number of nonzeros in any row of A, plus 1
411 * EPS is machine epsilon
412 *
413 * The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
414 * is incremented by SAFE1 if the i-th component of
415 * abs(op(A))*abs(X) + abs(B) is less than SAFE2.
416 *
417 * Use SLACN2 to estimate the infinity-norm of the matrix
418 * inv(op(A)) * diag(W),
419 * where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
420 *
421  DO 60 i = 1, n
422  IF( work( i ).GT.safe2 ) THEN
423  work( i ) = abs( work( n+i ) ) + nz*eps*work( i )
424  ELSE
425  work( i ) = abs( work( n+i ) ) + nz*eps*work( i ) + safe1
426  END IF
427  60 CONTINUE
428 *
429  kase = 0
430  70 CONTINUE
431  CALL slacn2( n, work( 2*n+1 ), work( n+1 ), iwork, ferr( j ),
432  $ kase, isave )
433  IF( kase.NE.0 ) THEN
434  IF( kase.EQ.1 ) THEN
435 *
436 * Multiply by diag(W)*inv(op(A)**T).
437 *
438  CALL sgttrs( transt, n, 1, dlf, df, duf, du2, ipiv,
439  $ work( n+1 ), n, info )
440  DO 80 i = 1, n
441  work( n+i ) = work( i )*work( n+i )
442  80 CONTINUE
443  ELSE
444 *
445 * Multiply by inv(op(A))*diag(W).
446 *
447  DO 90 i = 1, n
448  work( n+i ) = work( i )*work( n+i )
449  90 CONTINUE
450  CALL sgttrs( transn, n, 1, dlf, df, duf, du2, ipiv,
451  $ work( n+1 ), n, info )
452  END IF
453  GO TO 70
454  END IF
455 *
456 * Normalize error.
457 *
458  lstres = zero
459  DO 100 i = 1, n
460  lstres = max( lstres, abs( x( i, j ) ) )
461  100 CONTINUE
462  IF( lstres.NE.zero )
463  $ ferr( j ) = ferr( j ) / lstres
464 *
465  110 CONTINUE
466 *
467  RETURN
468 *
469 * End of SGTRFS
470 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine sgttrs(TRANS, N, NRHS, DL, D, DU, DU2, IPIV, B, LDB, INFO)
SGTTRS
Definition: sgttrs.f:138
subroutine slacn2(N, V, X, ISGN, EST, KASE, ISAVE)
SLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: slacn2.f:136
subroutine slagtm(TRANS, N, NRHS, ALPHA, DL, D, DU, X, LDX, BETA, B, LDB)
SLAGTM performs a matrix-matrix product of the form C = αAB+βC, where A is a tridiagonal matrix,...
Definition: slagtm.f:145
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
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: