LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dgttrs
subroutine dgttrs(trans, n, nrhs, dl, d, du, du2, ipiv, b, ldb, info)
DGTTRS
Definition dgttrs.f:138
dlacn2
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
dlagtm
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
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
Here is the call graph for this function:
Here is the caller graph for this function: