LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.
Date
September 2012

Definition at line 211 of file sgtrfs.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: