LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ dgtrfs()

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

DGTRFS

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

Purpose:
 DGTRFS 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 DOUBLE PRECISION array, dimension (N-1)
          The (n-1) subdiagonal elements of A.
[in]D
          D is DOUBLE PRECISION array, dimension (N)
          The diagonal elements of A.
[in]DU
          DU is DOUBLE PRECISION array, dimension (N-1)
          The (n-1) superdiagonal elements of A.
[in]DLF
          DLF is DOUBLE PRECISION array, dimension (N-1)
          The (n-1) multipliers that define the matrix L from the
          LU factorization of A as computed by DGTTRF.
[in]DF
          DF is DOUBLE PRECISION array, dimension (N)
          The n diagonal elements of the upper triangular matrix U from
          the LU factorization of A.
[in]DUF
          DUF is DOUBLE PRECISION array, dimension (N-1)
          The (n-1) elements of the first superdiagonal of U.
[in]DU2
          DU2 is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDX,NRHS)
          On entry, the solution matrix X, as computed by DGTTRS.
          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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 dgtrfs.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  DOUBLE PRECISION 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  DOUBLE PRECISION ZERO, ONE
231  parameter( zero = 0.0d+0, one = 1.0d+0 )
232  DOUBLE PRECISION TWO
233  parameter( two = 2.0d+0 )
234  DOUBLE PRECISION THREE
235  parameter( three = 3.0d+0 )
236 * ..
237 * .. Local Scalars ..
238  LOGICAL NOTRAN
239  CHARACTER TRANSN, TRANST
240  INTEGER COUNT, I, J, KASE, NZ
241  DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN
242 * ..
243 * .. Local Arrays ..
244  INTEGER ISAVE( 3 )
245 * ..
246 * .. External Subroutines ..
247  EXTERNAL daxpy, dcopy, dgttrs, dlacn2, dlagtm, xerbla
248 * ..
249 * .. Intrinsic Functions ..
250  INTRINSIC abs, max
251 * ..
252 * .. External Functions ..
253  LOGICAL LSAME
254  DOUBLE PRECISION DLAMCH
255  EXTERNAL lsame, dlamch
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( 'DGTRFS', -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 = dlamch( 'Epsilon' )
302  safmin = dlamch( '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 dcopy( n, b( 1, j ), 1, work( n+1 ), 1 )
320  CALL dlagtm( 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 dgttrs( trans, n, 1, dlf, df, duf, du2, ipiv,
392  $ work( n+1 ), n, info )
393  CALL daxpy( 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 DLACN2 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 dlacn2( 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 dgttrs( 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 dgttrs( 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 DGTRFS
470 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:89
subroutine dgttrs(TRANS, N, NRHS, DL, D, DU, DU2, IPIV, B, LDB, INFO)
DGTTRS
Definition: dgttrs.f:138
subroutine dlacn2(N, V, X, ISGN, EST, KASE, ISAVE)
DLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: dlacn2.f:136
subroutine dlagtm(TRANS, N, NRHS, ALPHA, DL, D, DU, X, LDX, BETA, B, LDB)
DLAGTM performs a matrix-matrix product of the form C = αAB+βC, where A is a tridiagonal matrix,...
Definition: dlagtm.f:145
Here is the call graph for this function:
Here is the caller graph for this function: