LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dpbrfs ( character  UPLO,
integer  N,
integer  KD,
integer  NRHS,
double precision, dimension( ldab, * )  AB,
integer  LDAB,
double precision, dimension( ldafb, * )  AFB,
integer  LDAFB,
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 
)

DPBRFS

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

Purpose: 
 DPBRFS improves the computed solution to a system of linear
 equations when the coefficient matrix is symmetric positive definite
 and banded, and provides error bounds and backward error estimates
 for the solution.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in]KD
          KD is INTEGER
          The number of superdiagonals of the matrix A if UPLO = 'U',
          or the number of subdiagonals if UPLO = 'L'.  KD >= 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 upper or lower triangle of the symmetric band matrix A,
          stored in the first KD+1 rows of the array.  The j-th column
          of A is stored in the j-th column of the array AB as follows:
          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
[in]LDAB
          LDAB is INTEGER
          The leading dimension of the array AB.  LDAB >= KD+1.
[in]AFB
          AFB is DOUBLE PRECISION array, dimension (LDAFB,N)
          The triangular factor U or L from the Cholesky factorization
          A = U**T*U or A = L*L**T of the band matrix A as computed by
          DPBTRF, in the same storage format as A (see AB).
[in]LDAFB
          LDAFB is INTEGER
          The leading dimension of the array AFB.  LDAFB >= KD+1.
[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 DPBTRS.
          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 191 of file dpbrfs.f.

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