LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
caxpy
subroutine caxpy(n, ca, cx, incx, cy, incy)
CAXPY
Definition caxpy.f:88
ccopy
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
chbmv
subroutine chbmv(uplo, n, k, alpha, a, lda, x, incx, beta, y, incy)
CHBMV
Definition chbmv.f:187
clacn2
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
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
cpbtrs
subroutine cpbtrs(uplo, n, kd, nrhs, ab, ldab, b, ldb, info)
CPBTRS
Definition cpbtrs.f:121
Here is the call graph for this function:
Here is the caller graph for this function: