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

◆ cgbrfs()

subroutine cgbrfs ( character  trans,
integer  n,
integer  kl,
integer  ku,
integer  nrhs,
complex, dimension( ldab, * )  ab,
integer  ldab,
complex, dimension( ldafb, * )  afb,
integer  ldafb,
integer, dimension( * )  ipiv,
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 
)

CGBRFS

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

Purpose:
 CGBRFS 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)
[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 COMPLEX 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 COMPLEX array, dimension (LDAFB,N)
          Details of the LU factorization of the band matrix A, as
          computed by CGBTRF.  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 CGBTRF; for 1<=i<=N, row i of the
          matrix was interchanged with row IPIV(i).
[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 CGBTRS.
          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 203 of file cgbrfs.f.

206*
207* -- LAPACK computational routine --
208* -- LAPACK is a software package provided by Univ. of Tennessee, --
209* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
210*
211* .. Scalar Arguments ..
212 CHARACTER TRANS
213 INTEGER INFO, KL, KU, LDAB, LDAFB, LDB, LDX, N, NRHS
214* ..
215* .. Array Arguments ..
216 INTEGER IPIV( * )
217 REAL BERR( * ), FERR( * ), RWORK( * )
218 COMPLEX AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
219 $ WORK( * ), X( LDX, * )
220* ..
221*
222* =====================================================================
223*
224* .. Parameters ..
225 INTEGER ITMAX
226 parameter( itmax = 5 )
227 REAL ZERO
228 parameter( zero = 0.0e+0 )
229 COMPLEX CONE
230 parameter( cone = ( 1.0e+0, 0.0e+0 ) )
231 REAL TWO
232 parameter( two = 2.0e+0 )
233 REAL THREE
234 parameter( three = 3.0e+0 )
235* ..
236* .. Local Scalars ..
237 LOGICAL NOTRAN
238 CHARACTER TRANSN, TRANST
239 INTEGER COUNT, I, J, K, KASE, KK, NZ
240 REAL EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
241 COMPLEX ZDUM
242* ..
243* .. Local Arrays ..
244 INTEGER ISAVE( 3 )
245* ..
246* .. External Subroutines ..
247 EXTERNAL caxpy, ccopy, cgbmv, cgbtrs, clacn2, xerbla
248* ..
249* .. Intrinsic Functions ..
250 INTRINSIC abs, aimag, max, min, real
251* ..
252* .. External Functions ..
253 LOGICAL LSAME
254 REAL SLAMCH
255 EXTERNAL lsame, slamch
256* ..
257* .. Statement Functions ..
258 REAL CABS1
259* ..
260* .. Statement Function definitions ..
261 cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( zdum ) )
262* ..
263* .. Executable Statements ..
264*
265* Test the input parameters.
266*
267 info = 0
268 notran = lsame( trans, 'N' )
269 IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
270 $ lsame( trans, 'C' ) ) THEN
271 info = -1
272 ELSE IF( n.LT.0 ) THEN
273 info = -2
274 ELSE IF( kl.LT.0 ) THEN
275 info = -3
276 ELSE IF( ku.LT.0 ) THEN
277 info = -4
278 ELSE IF( nrhs.LT.0 ) THEN
279 info = -5
280 ELSE IF( ldab.LT.kl+ku+1 ) THEN
281 info = -7
282 ELSE IF( ldafb.LT.2*kl+ku+1 ) THEN
283 info = -9
284 ELSE IF( ldb.LT.max( 1, n ) ) THEN
285 info = -12
286 ELSE IF( ldx.LT.max( 1, n ) ) THEN
287 info = -14
288 END IF
289 IF( info.NE.0 ) THEN
290 CALL xerbla( 'CGBRFS', -info )
291 RETURN
292 END IF
293*
294* Quick return if possible
295*
296 IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
297 DO 10 j = 1, nrhs
298 ferr( j ) = zero
299 berr( j ) = zero
300 10 CONTINUE
301 RETURN
302 END IF
303*
304 IF( notran ) THEN
305 transn = 'N'
306 transt = 'C'
307 ELSE
308 transn = 'C'
309 transt = 'N'
310 END IF
311*
312* NZ = maximum number of nonzero elements in each row of A, plus 1
313*
314 nz = min( kl+ku+2, n+1 )
315 eps = slamch( 'Epsilon' )
316 safmin = slamch( 'Safe minimum' )
317 safe1 = nz*safmin
318 safe2 = safe1 / eps
319*
320* Do for each right hand side
321*
322 DO 140 j = 1, nrhs
323*
324 count = 1
325 lstres = three
326 20 CONTINUE
327*
328* Loop until stopping criterion is satisfied.
329*
330* Compute residual R = B - op(A) * X,
331* where op(A) = A, A**T, or A**H, depending on TRANS.
332*
333 CALL ccopy( n, b( 1, j ), 1, work, 1 )
334 CALL cgbmv( trans, n, n, kl, ku, -cone, ab, ldab, x( 1, j ), 1,
335 $ cone, work, 1 )
336*
337* Compute componentwise relative backward error from formula
338*
339* max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
340*
341* where abs(Z) is the componentwise absolute value of the matrix
342* or vector Z. If the i-th component of the denominator is less
343* than SAFE2, then SAFE1 is added to the i-th components of the
344* numerator and denominator before dividing.
345*
346 DO 30 i = 1, n
347 rwork( i ) = cabs1( b( i, j ) )
348 30 CONTINUE
349*
350* Compute abs(op(A))*abs(X) + abs(B).
351*
352 IF( notran ) THEN
353 DO 50 k = 1, n
354 kk = ku + 1 - k
355 xk = cabs1( x( k, j ) )
356 DO 40 i = max( 1, k-ku ), min( n, k+kl )
357 rwork( i ) = rwork( i ) + cabs1( ab( kk+i, k ) )*xk
358 40 CONTINUE
359 50 CONTINUE
360 ELSE
361 DO 70 k = 1, n
362 s = zero
363 kk = ku + 1 - k
364 DO 60 i = max( 1, k-ku ), min( n, k+kl )
365 s = s + cabs1( ab( kk+i, k ) )*cabs1( x( i, j ) )
366 60 CONTINUE
367 rwork( k ) = rwork( k ) + s
368 70 CONTINUE
369 END IF
370 s = zero
371 DO 80 i = 1, n
372 IF( rwork( i ).GT.safe2 ) THEN
373 s = max( s, cabs1( work( i ) ) / rwork( i ) )
374 ELSE
375 s = max( s, ( cabs1( work( i ) )+safe1 ) /
376 $ ( rwork( i )+safe1 ) )
377 END IF
378 80 CONTINUE
379 berr( j ) = s
380*
381* Test stopping criterion. Continue iterating if
382* 1) The residual BERR(J) is larger than machine epsilon, and
383* 2) BERR(J) decreased by at least a factor of 2 during the
384* last iteration, and
385* 3) At most ITMAX iterations tried.
386*
387 IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
388 $ count.LE.itmax ) THEN
389*
390* Update solution and try again.
391*
392 CALL cgbtrs( trans, n, kl, ku, 1, afb, ldafb, ipiv, work, n,
393 $ info )
394 CALL caxpy( n, cone, work, 1, x( 1, j ), 1 )
395 lstres = berr( j )
396 count = count + 1
397 GO TO 20
398 END IF
399*
400* Bound error from formula
401*
402* norm(X - XTRUE) / norm(X) .le. FERR =
403* norm( abs(inv(op(A)))*
404* ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
405*
406* where
407* norm(Z) is the magnitude of the largest component of Z
408* inv(op(A)) is the inverse of op(A)
409* abs(Z) is the componentwise absolute value of the matrix or
410* vector Z
411* NZ is the maximum number of nonzeros in any row of A, plus 1
412* EPS is machine epsilon
413*
414* The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
415* is incremented by SAFE1 if the i-th component of
416* abs(op(A))*abs(X) + abs(B) is less than SAFE2.
417*
418* Use CLACN2 to estimate the infinity-norm of the matrix
419* inv(op(A)) * diag(W),
420* where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
421*
422 DO 90 i = 1, n
423 IF( rwork( i ).GT.safe2 ) THEN
424 rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i )
425 ELSE
426 rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i ) +
427 $ safe1
428 END IF
429 90 CONTINUE
430*
431 kase = 0
432 100 CONTINUE
433 CALL clacn2( n, work( n+1 ), work, ferr( j ), kase, isave )
434 IF( kase.NE.0 ) THEN
435 IF( kase.EQ.1 ) THEN
436*
437* Multiply by diag(W)*inv(op(A)**H).
438*
439 CALL cgbtrs( transt, n, kl, ku, 1, afb, ldafb, ipiv,
440 $ work, n, info )
441 DO 110 i = 1, n
442 work( i ) = rwork( i )*work( i )
443 110 CONTINUE
444 ELSE
445*
446* Multiply by inv(op(A))*diag(W).
447*
448 DO 120 i = 1, n
449 work( i ) = rwork( i )*work( i )
450 120 CONTINUE
451 CALL cgbtrs( transn, n, kl, ku, 1, afb, ldafb, ipiv,
452 $ work, n, info )
453 END IF
454 GO TO 100
455 END IF
456*
457* Normalize error.
458*
459 lstres = zero
460 DO 130 i = 1, n
461 lstres = max( lstres, cabs1( x( i, j ) ) )
462 130 CONTINUE
463 IF( lstres.NE.zero )
464 $ ferr( j ) = ferr( j ) / lstres
465*
466 140 CONTINUE
467*
468 RETURN
469*
470* End of CGBRFS
471*
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
cgbmv
subroutine cgbmv(trans, m, n, kl, ku, alpha, a, lda, x, incx, beta, y, incy)
CGBMV
Definition cgbmv.f:190
cgbtrs
subroutine cgbtrs(trans, n, kl, ku, nrhs, ab, ldab, ipiv, b, ldb, info)
CGBTRS
Definition cgbtrs.f:138
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
Here is the call graph for this function:
Here is the caller graph for this function: