LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ 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 *
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 cgbmv(TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGBMV
Definition: cgbmv.f:187
subroutine cgbtrs(TRANS, N, KL, KU, NRHS, AB, LDAB, IPIV, B, LDB, INFO)
CGBTRS
Definition: cgbtrs.f:138
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
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: