 LAPACK  3.10.1 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.`

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: