LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ dgbrfs()

subroutine dgbrfs ( character  TRANS,
integer  N,
integer  KL,
integer  KU,
integer  NRHS,
double precision, dimension( ldab, * )  AB,
integer  LDAB,
double precision, dimension( ldafb, * )  AFB,
integer  LDAFB,
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 
)

DGBRFS

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

Purpose:
 DGBRFS improves the computed solution to a system of linear
 equations when the coefficient matrix is banded, 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]KL
          KL is INTEGER
          The number of subdiagonals within the band of A.  KL >= 0.
[in]KU
          KU is INTEGER
          The number of superdiagonals within the band of A.  KU >= 0.
[in]NRHS
          NRHS is INTEGER
          The number of right hand sides, i.e., the number of columns
          of the matrices B and X.  NRHS >= 0.
[in]AB
          AB is DOUBLE PRECISION array, dimension (LDAB,N)
          The original band matrix A, stored in rows 1 to KL+KU+1.
          The j-th column of A is stored in the j-th column of the
          array AB as follows:
          AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(n,j+kl).
[in]LDAB
          LDAB is INTEGER
          The leading dimension of the array AB.  LDAB >= KL+KU+1.
[in]AFB
          AFB is DOUBLE PRECISION array, dimension (LDAFB,N)
          Details of the LU factorization of the band matrix A, as
          computed by DGBTRF.  U is stored as an upper triangular band
          matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
          the multipliers used during the factorization are stored in
          rows KL+KU+2 to 2*KL+KU+1.
[in]LDAFB
          LDAFB is INTEGER
          The leading dimension of the array AFB.  LDAFB >= 2*KL*KU+1.
[in]IPIV
          IPIV is INTEGER array, dimension (N)
          The pivot indices from DGBTRF; for 1<=i<=N, row i of the
          matrix was interchanged with row IPIV(i).
[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 DGBTRS.
          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 202 of file dgbrfs.f.

205 *
206 * -- LAPACK computational routine --
207 * -- LAPACK is a software package provided by Univ. of Tennessee, --
208 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
209 *
210 * .. Scalar Arguments ..
211  CHARACTER TRANS
212  INTEGER INFO, KL, KU, LDAB, LDAFB, LDB, LDX, N, NRHS
213 * ..
214 * .. Array Arguments ..
215  INTEGER IPIV( * ), IWORK( * )
216  DOUBLE PRECISION AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
217  $ BERR( * ), FERR( * ), WORK( * ), X( LDX, * )
218 * ..
219 *
220 * =====================================================================
221 *
222 * .. Parameters ..
223  INTEGER ITMAX
224  parameter( itmax = 5 )
225  DOUBLE PRECISION ZERO
226  parameter( zero = 0.0d+0 )
227  DOUBLE PRECISION ONE
228  parameter( one = 1.0d+0 )
229  DOUBLE PRECISION TWO
230  parameter( two = 2.0d+0 )
231  DOUBLE PRECISION THREE
232  parameter( three = 3.0d+0 )
233 * ..
234 * .. Local Scalars ..
235  LOGICAL NOTRAN
236  CHARACTER TRANST
237  INTEGER COUNT, I, J, K, KASE, KK, NZ
238  DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
239 * ..
240 * .. Local Arrays ..
241  INTEGER ISAVE( 3 )
242 * ..
243 * .. External Subroutines ..
244  EXTERNAL daxpy, dcopy, dgbmv, dgbtrs, dlacn2, xerbla
245 * ..
246 * .. Intrinsic Functions ..
247  INTRINSIC abs, max, min
248 * ..
249 * .. External Functions ..
250  LOGICAL LSAME
251  DOUBLE PRECISION DLAMCH
252  EXTERNAL lsame, dlamch
253 * ..
254 * .. Executable Statements ..
255 *
256 * Test the input parameters.
257 *
258  info = 0
259  notran = lsame( trans, 'N' )
260  IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
261  $ lsame( trans, 'C' ) ) THEN
262  info = -1
263  ELSE IF( n.LT.0 ) THEN
264  info = -2
265  ELSE IF( kl.LT.0 ) THEN
266  info = -3
267  ELSE IF( ku.LT.0 ) THEN
268  info = -4
269  ELSE IF( nrhs.LT.0 ) THEN
270  info = -5
271  ELSE IF( ldab.LT.kl+ku+1 ) THEN
272  info = -7
273  ELSE IF( ldafb.LT.2*kl+ku+1 ) THEN
274  info = -9
275  ELSE IF( ldb.LT.max( 1, n ) ) THEN
276  info = -12
277  ELSE IF( ldx.LT.max( 1, n ) ) THEN
278  info = -14
279  END IF
280  IF( info.NE.0 ) THEN
281  CALL xerbla( 'DGBRFS', -info )
282  RETURN
283  END IF
284 *
285 * Quick return if possible
286 *
287  IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
288  DO 10 j = 1, nrhs
289  ferr( j ) = zero
290  berr( j ) = zero
291  10 CONTINUE
292  RETURN
293  END IF
294 *
295  IF( notran ) THEN
296  transt = 'T'
297  ELSE
298  transt = 'N'
299  END IF
300 *
301 * NZ = maximum number of nonzero elements in each row of A, plus 1
302 *
303  nz = min( kl+ku+2, n+1 )
304  eps = dlamch( 'Epsilon' )
305  safmin = dlamch( 'Safe minimum' )
306  safe1 = nz*safmin
307  safe2 = safe1 / eps
308 *
309 * Do for each right hand side
310 *
311  DO 140 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 dcopy( n, b( 1, j ), 1, work( n+1 ), 1 )
323  CALL dgbmv( trans, n, n, kl, ku, -one, ab, ldab, x( 1, j ), 1,
324  $ one, work( n+1 ), 1 )
325 *
326 * Compute componentwise relative backward error from formula
327 *
328 * max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
329 *
330 * where abs(Z) is the componentwise absolute value of the matrix
331 * or vector Z. If the i-th component of the denominator is less
332 * than SAFE2, then SAFE1 is added to the i-th components of the
333 * numerator and denominator before dividing.
334 *
335  DO 30 i = 1, n
336  work( i ) = abs( b( i, j ) )
337  30 CONTINUE
338 *
339 * Compute abs(op(A))*abs(X) + abs(B).
340 *
341  IF( notran ) THEN
342  DO 50 k = 1, n
343  kk = ku + 1 - k
344  xk = abs( x( k, j ) )
345  DO 40 i = max( 1, k-ku ), min( n, k+kl )
346  work( i ) = work( i ) + abs( ab( kk+i, k ) )*xk
347  40 CONTINUE
348  50 CONTINUE
349  ELSE
350  DO 70 k = 1, n
351  s = zero
352  kk = ku + 1 - k
353  DO 60 i = max( 1, k-ku ), min( n, k+kl )
354  s = s + abs( ab( kk+i, k ) )*abs( x( i, j ) )
355  60 CONTINUE
356  work( k ) = work( k ) + s
357  70 CONTINUE
358  END IF
359  s = zero
360  DO 80 i = 1, n
361  IF( work( i ).GT.safe2 ) THEN
362  s = max( s, abs( work( n+i ) ) / work( i ) )
363  ELSE
364  s = max( s, ( abs( work( n+i ) )+safe1 ) /
365  $ ( work( i )+safe1 ) )
366  END IF
367  80 CONTINUE
368  berr( j ) = s
369 *
370 * Test stopping criterion. Continue iterating if
371 * 1) The residual BERR(J) is larger than machine epsilon, and
372 * 2) BERR(J) decreased by at least a factor of 2 during the
373 * last iteration, and
374 * 3) At most ITMAX iterations tried.
375 *
376  IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
377  $ count.LE.itmax ) THEN
378 *
379 * Update solution and try again.
380 *
381  CALL dgbtrs( trans, n, kl, ku, 1, afb, ldafb, ipiv,
382  $ work( n+1 ), n, info )
383  CALL daxpy( n, one, work( n+1 ), 1, x( 1, j ), 1 )
384  lstres = berr( j )
385  count = count + 1
386  GO TO 20
387  END IF
388 *
389 * Bound error from formula
390 *
391 * norm(X - XTRUE) / norm(X) .le. FERR =
392 * norm( abs(inv(op(A)))*
393 * ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
394 *
395 * where
396 * norm(Z) is the magnitude of the largest component of Z
397 * inv(op(A)) is the inverse of op(A)
398 * abs(Z) is the componentwise absolute value of the matrix or
399 * vector Z
400 * NZ is the maximum number of nonzeros in any row of A, plus 1
401 * EPS is machine epsilon
402 *
403 * The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
404 * is incremented by SAFE1 if the i-th component of
405 * abs(op(A))*abs(X) + abs(B) is less than SAFE2.
406 *
407 * Use DLACN2 to estimate the infinity-norm of the matrix
408 * inv(op(A)) * diag(W),
409 * where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
410 *
411  DO 90 i = 1, n
412  IF( work( i ).GT.safe2 ) THEN
413  work( i ) = abs( work( n+i ) ) + nz*eps*work( i )
414  ELSE
415  work( i ) = abs( work( n+i ) ) + nz*eps*work( i ) + safe1
416  END IF
417  90 CONTINUE
418 *
419  kase = 0
420  100 CONTINUE
421  CALL dlacn2( n, work( 2*n+1 ), work( n+1 ), iwork, ferr( j ),
422  $ kase, isave )
423  IF( kase.NE.0 ) THEN
424  IF( kase.EQ.1 ) THEN
425 *
426 * Multiply by diag(W)*inv(op(A)**T).
427 *
428  CALL dgbtrs( transt, n, kl, ku, 1, afb, ldafb, ipiv,
429  $ work( n+1 ), n, info )
430  DO 110 i = 1, n
431  work( n+i ) = work( n+i )*work( i )
432  110 CONTINUE
433  ELSE
434 *
435 * Multiply by inv(op(A))*diag(W).
436 *
437  DO 120 i = 1, n
438  work( n+i ) = work( n+i )*work( i )
439  120 CONTINUE
440  CALL dgbtrs( trans, n, kl, ku, 1, afb, ldafb, ipiv,
441  $ work( n+1 ), n, info )
442  END IF
443  GO TO 100
444  END IF
445 *
446 * Normalize error.
447 *
448  lstres = zero
449  DO 130 i = 1, n
450  lstres = max( lstres, abs( x( i, j ) ) )
451  130 CONTINUE
452  IF( lstres.NE.zero )
453  $ ferr( j ) = ferr( j ) / lstres
454 *
455  140 CONTINUE
456 *
457  RETURN
458 *
459 * End of DGBRFS
460 *
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 dgbmv(TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGBMV
Definition: dgbmv.f:185
subroutine dgbtrs(TRANS, N, KL, KU, NRHS, AB, LDAB, IPIV, B, LDB, INFO)
DGBTRS
Definition: dgbtrs.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
Here is the call graph for this function:
Here is the caller graph for this function: