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

◆ 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.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 201 of file zgbrfs.f.

205*
206* -- LAPACK computational routine --
207* -- LAPACK is a software package provided by Univ. of Tennessee, --
208* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
209*
210* .. Scalar Arguments ..
211 CHARACTER TRANS
212 INTEGER INFO, KL, KU, LDAB, LDAFB, LDB, LDX, N, NRHS
213* ..
214* .. Array Arguments ..
215 INTEGER IPIV( * )
216 DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
217 COMPLEX*16 AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
218 $ WORK( * ), X( LDX, * )
219* ..
220*
221* =====================================================================
222*
223* .. Parameters ..
224 INTEGER ITMAX
225 parameter( itmax = 5 )
226 DOUBLE PRECISION ZERO
227 parameter( zero = 0.0d+0 )
228 COMPLEX*16 CONE
229 parameter( cone = ( 1.0d+0, 0.0d+0 ) )
230 DOUBLE PRECISION TWO
231 parameter( two = 2.0d+0 )
232 DOUBLE PRECISION THREE
233 parameter( three = 3.0d+0 )
234* ..
235* .. Local Scalars ..
236 LOGICAL NOTRAN
237 CHARACTER TRANSN, TRANST
238 INTEGER COUNT, I, J, K, KASE, KK, NZ
239 DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
240 COMPLEX*16 ZDUM
241* ..
242* .. Local Arrays ..
243 INTEGER ISAVE( 3 )
244* ..
245* .. External Subroutines ..
246 EXTERNAL xerbla, zaxpy, zcopy, zgbmv, zgbtrs,
247 $ 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 ),
335 $ 1,
336 $ cone, work, 1 )
337*
338* Compute componentwise relative backward error from formula
339*
340* max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
341*
342* where abs(Z) is the componentwise absolute value of the matrix
343* or vector Z. If the i-th component of the denominator is less
344* than SAFE2, then SAFE1 is added to the i-th components of the
345* numerator and denominator before dividing.
346*
347 DO 30 i = 1, n
348 rwork( i ) = cabs1( b( i, j ) )
349 30 CONTINUE
350*
351* Compute abs(op(A))*abs(X) + abs(B).
352*
353 IF( notran ) THEN
354 DO 50 k = 1, n
355 kk = ku + 1 - k
356 xk = cabs1( x( k, j ) )
357 DO 40 i = max( 1, k-ku ), min( n, k+kl )
358 rwork( i ) = rwork( i ) + cabs1( ab( kk+i, k ) )*xk
359 40 CONTINUE
360 50 CONTINUE
361 ELSE
362 DO 70 k = 1, n
363 s = zero
364 kk = ku + 1 - k
365 DO 60 i = max( 1, k-ku ), min( n, k+kl )
366 s = s + cabs1( ab( kk+i, k ) )*cabs1( x( i, j ) )
367 60 CONTINUE
368 rwork( k ) = rwork( k ) + s
369 70 CONTINUE
370 END IF
371 s = zero
372 DO 80 i = 1, n
373 IF( rwork( i ).GT.safe2 ) THEN
374 s = max( s, cabs1( work( i ) ) / rwork( i ) )
375 ELSE
376 s = max( s, ( cabs1( work( i ) )+safe1 ) /
377 $ ( rwork( i )+safe1 ) )
378 END IF
379 80 CONTINUE
380 berr( j ) = s
381*
382* Test stopping criterion. Continue iterating if
383* 1) The residual BERR(J) is larger than machine epsilon, and
384* 2) BERR(J) decreased by at least a factor of 2 during the
385* last iteration, and
386* 3) At most ITMAX iterations tried.
387*
388 IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
389 $ count.LE.itmax ) THEN
390*
391* Update solution and try again.
392*
393 CALL zgbtrs( trans, n, kl, ku, 1, afb, ldafb, ipiv, work,
394 $ n,
395 $ info )
396 CALL zaxpy( n, cone, work, 1, x( 1, j ), 1 )
397 lstres = berr( j )
398 count = count + 1
399 GO TO 20
400 END IF
401*
402* Bound error from formula
403*
404* norm(X - XTRUE) / norm(X) .le. FERR =
405* norm( abs(inv(op(A)))*
406* ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
407*
408* where
409* norm(Z) is the magnitude of the largest component of Z
410* inv(op(A)) is the inverse of op(A)
411* abs(Z) is the componentwise absolute value of the matrix or
412* vector Z
413* NZ is the maximum number of nonzeros in any row of A, plus 1
414* EPS is machine epsilon
415*
416* The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
417* is incremented by SAFE1 if the i-th component of
418* abs(op(A))*abs(X) + abs(B) is less than SAFE2.
419*
420* Use ZLACN2 to estimate the infinity-norm of the matrix
421* inv(op(A)) * diag(W),
422* where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
423*
424 DO 90 i = 1, n
425 IF( rwork( i ).GT.safe2 ) THEN
426 rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i )
427 ELSE
428 rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i ) +
429 $ safe1
430 END IF
431 90 CONTINUE
432*
433 kase = 0
434 100 CONTINUE
435 CALL zlacn2( n, work( n+1 ), work, ferr( j ), kase, isave )
436 IF( kase.NE.0 ) THEN
437 IF( kase.EQ.1 ) THEN
438*
439* Multiply by diag(W)*inv(op(A)**H).
440*
441 CALL zgbtrs( transt, n, kl, ku, 1, afb, ldafb, ipiv,
442 $ work, n, info )
443 DO 110 i = 1, n
444 work( i ) = rwork( i )*work( i )
445 110 CONTINUE
446 ELSE
447*
448* Multiply by inv(op(A))*diag(W).
449*
450 DO 120 i = 1, n
451 work( i ) = rwork( i )*work( i )
452 120 CONTINUE
453 CALL zgbtrs( transn, n, kl, ku, 1, afb, ldafb, ipiv,
454 $ work, n, info )
455 END IF
456 GO TO 100
457 END IF
458*
459* Normalize error.
460*
461 lstres = zero
462 DO 130 i = 1, n
463 lstres = max( lstres, cabs1( x( i, j ) ) )
464 130 CONTINUE
465 IF( lstres.NE.zero )
466 $ ferr( j ) = ferr( j ) / lstres
467*
468 140 CONTINUE
469*
470 RETURN
471*
472* End of ZGBRFS
473*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zaxpy
subroutine zaxpy(n, za, zx, incx, zy, incy)
ZAXPY
Definition zaxpy.f:88
zcopy
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
zgbmv
subroutine zgbmv(trans, m, n, kl, ku, alpha, a, lda, x, incx, beta, y, incy)
ZGBMV
Definition zgbmv.f:190
zgbtrs
subroutine zgbtrs(trans, n, kl, ku, nrhs, ab, ldab, ipiv, b, ldb, info)
ZGBTRS
Definition zgbtrs.f:137
zlacn2
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:131
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
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: