LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ cpbrfs()

subroutine cpbrfs ( character  UPLO,
integer  N,
integer  KD,
integer  NRHS,
complex, dimension( ldab, * )  AB,
integer  LDAB,
complex, dimension( ldafb, * )  AFB,
integer  LDAFB,
complex, dimension( ldb, * )  B,
integer  LDB,
complex, dimension( ldx, * )  X,
integer  LDX,
real, dimension( * )  FERR,
real, dimension( * )  BERR,
complex, dimension( * )  WORK,
real, dimension( * )  RWORK,
integer  INFO 
)

CPBRFS

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

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