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

◆ zpbrfs()

subroutine zpbrfs ( character  uplo,
integer  n,
integer  kd,
integer  nrhs,
complex*16, dimension( ldab, * )  ab,
integer  ldab,
complex*16, dimension( ldafb, * )  afb,
integer  ldafb,
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 
)

ZPBRFS

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

Purpose:
 ZPBRFS 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*16 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*16 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
          ZPBTRF, 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*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 ZPBTRS.
          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 187 of file zpbrfs.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 DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
200 COMPLEX*16 AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
201 $ WORK( * ), X( LDX, * )
202* ..
203*
204* =====================================================================
205*
206* .. Parameters ..
207 INTEGER ITMAX
208 parameter( itmax = 5 )
209 DOUBLE PRECISION ZERO
210 parameter( zero = 0.0d+0 )
211 COMPLEX*16 ONE
212 parameter( one = ( 1.0d+0, 0.0d+0 ) )
213 DOUBLE PRECISION TWO
214 parameter( two = 2.0d+0 )
215 DOUBLE PRECISION THREE
216 parameter( three = 3.0d+0 )
217* ..
218* .. Local Scalars ..
219 LOGICAL UPPER
220 INTEGER COUNT, I, J, K, KASE, L, NZ
221 DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
222 COMPLEX*16 ZDUM
223* ..
224* .. Local Arrays ..
225 INTEGER ISAVE( 3 )
226* ..
227* .. External Subroutines ..
228 EXTERNAL xerbla, zaxpy, zcopy, zhbmv, zlacn2, zpbtrs
229* ..
230* .. Intrinsic Functions ..
231 INTRINSIC abs, dble, dimag, max, min
232* ..
233* .. External Functions ..
234 LOGICAL LSAME
235 DOUBLE PRECISION DLAMCH
236 EXTERNAL lsame, dlamch
237* ..
238* .. Statement Functions ..
239 DOUBLE PRECISION CABS1
240* ..
241* .. Statement Function definitions ..
242 cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( 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( 'ZPBRFS', -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 = dlamch( 'Epsilon' )
286 safmin = dlamch( '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 zcopy( n, b( 1, j ), 1, work, 1 )
303 CALL zhbmv( 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( dble( 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( dble( 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 zpbtrs( uplo, n, kd, 1, afb, ldafb, work, n, info )
369 CALL zaxpy( 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 ZLACN2 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 zlacn2( 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 zpbtrs( 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 zpbtrs( 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 ZPBRFS
444*
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
zhbmv
subroutine zhbmv(uplo, n, k, alpha, a, lda, x, incx, beta, y, incy)
ZHBMV
Definition zhbmv.f:187
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:133
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
zpbtrs
subroutine zpbtrs(uplo, n, kd, nrhs, ab, ldab, b, ldb, info)
ZPBTRS
Definition zpbtrs.f:121
Here is the call graph for this function:
Here is the caller graph for this function: