LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ dpbrfs()

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.

Definition at line 187 of file dpbrfs.f.

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