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

Definition at line 207 of file dgbrfs.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: