 LAPACK  3.10.1 LAPACK: Linear Algebra PACKage

## ◆ zgbrfs()

 subroutine zgbrfs ( character TRANS, integer N, integer KL, integer KU, integer NRHS, complex*16, dimension( ldab, * ) AB, integer LDAB, complex*16, dimension( ldafb, * ) AFB, integer LDAFB, integer, dimension( * ) IPIV, complex*16, dimension( ldb, * ) B, integer LDB, complex*16, dimension( ldx, * ) X, integer LDX, double precision, dimension( * ) FERR, double precision, dimension( * ) BERR, complex*16, dimension( * ) WORK, double precision, dimension( * ) RWORK, integer INFO )

ZGBRFS

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

Purpose:
``` ZGBRFS 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*16 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*16 array, dimension (LDAFB,N) Details of the LU factorization of the band matrix A, as computed by ZGBTRF. 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 ZGBTRF; for 1<=i<=N, row i of the matrix was interchanged with row IPIV(i).``` [in] B ``` B is COMPLEX*16 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*16 array, dimension (LDX,NRHS) On entry, the solution matrix X, as computed by ZGBTRS. 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 COMPLEX*16 array, dimension (2*N)` [out] RWORK ` RWORK is DOUBLE PRECISION 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.`

Definition at line 203 of file zgbrfs.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  DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
218  COMPLEX*16 AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
219  \$ WORK( * ), X( LDX, * )
220 * ..
221 *
222 * =====================================================================
223 *
224 * .. Parameters ..
225  INTEGER ITMAX
226  parameter( itmax = 5 )
227  DOUBLE PRECISION ZERO
228  parameter( zero = 0.0d+0 )
229  COMPLEX*16 CONE
230  parameter( cone = ( 1.0d+0, 0.0d+0 ) )
231  DOUBLE PRECISION TWO
232  parameter( two = 2.0d+0 )
233  DOUBLE PRECISION THREE
234  parameter( three = 3.0d+0 )
235 * ..
236 * .. Local Scalars ..
237  LOGICAL NOTRAN
238  CHARACTER TRANSN, TRANST
239  INTEGER COUNT, I, J, K, KASE, KK, NZ
240  DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
241  COMPLEX*16 ZDUM
242 * ..
243 * .. Local Arrays ..
244  INTEGER ISAVE( 3 )
245 * ..
246 * .. External Subroutines ..
247  EXTERNAL xerbla, zaxpy, zcopy, zgbmv, zgbtrs, zlacn2
248 * ..
249 * .. Intrinsic Functions ..
250  INTRINSIC abs, dble, dimag, max, min
251 * ..
252 * .. External Functions ..
253  LOGICAL LSAME
254  DOUBLE PRECISION DLAMCH
255  EXTERNAL lsame, dlamch
256 * ..
257 * .. Statement Functions ..
258  DOUBLE PRECISION CABS1
259 * ..
260 * .. Statement Function definitions ..
261  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( 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( 'ZGBRFS', -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 = dlamch( 'Epsilon' )
316  safmin = dlamch( '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 zcopy( n, b( 1, j ), 1, work, 1 )
334  CALL zgbmv( 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 zgbtrs( trans, n, kl, ku, 1, afb, ldafb, ipiv, work, n,
393  \$ info )
394  CALL zaxpy( 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 ZLACN2 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 zlacn2( 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 zgbtrs( 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 zgbtrs( 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 ZGBRFS
471 *
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 zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:88
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:81
subroutine zgbmv(TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGBMV
Definition: zgbmv.f:187
subroutine zgbtrs(TRANS, N, KL, KU, NRHS, AB, LDAB, IPIV, B, LDB, INFO)
ZGBTRS
Definition: zgbtrs.f:138
subroutine zlacn2(N, V, X, EST, KASE, ISAVE)
ZLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: zlacn2.f:133
Here is the call graph for this function:
Here is the caller graph for this function: