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

◆ spbrfs()

subroutine spbrfs ( character  uplo,
integer  n,
integer  kd,
integer  nrhs,
real, dimension( ldab, * )  ab,
integer  ldab,
real, dimension( ldafb, * )  afb,
integer  ldafb,
real, dimension( ldb, * )  b,
integer  ldb,
real, dimension( ldx, * )  x,
integer  ldx,
real, dimension( * )  ferr,
real, dimension( * )  berr,
real, dimension( * )  work,
integer, dimension( * )  iwork,
integer  info 
)

SPBRFS

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

Purpose:
 SPBRFS improves the computed solution to a system of linear
 equations when the coefficient matrix is symmetric 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 REAL array, dimension (LDAB,N)
          The upper or lower triangle of the symmetric 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 REAL array, dimension (LDAFB,N)
          The triangular factor U or L from the Cholesky factorization
          A = U**T*U or A = L*L**T of the band matrix A as computed by
          SPBTRF, 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 REAL 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 REAL array, dimension (LDX,NRHS)
          On entry, the solution matrix X, as computed by SPBTRS.
          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 REAL array, dimension (3*N)
[out]IWORK
          IWORK is INTEGER 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 spbrfs.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 INTEGER IWORK( * )
200 REAL AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
201 $ BERR( * ), FERR( * ), WORK( * ), X( LDX, * )
202* ..
203*
204* =====================================================================
205*
206* .. Parameters ..
207 INTEGER ITMAX
208 parameter( itmax = 5 )
209 REAL ZERO
210 parameter( zero = 0.0e+0 )
211 REAL ONE
212 parameter( one = 1.0e+0 )
213 REAL TWO
214 parameter( two = 2.0e+0 )
215 REAL THREE
216 parameter( three = 3.0e+0 )
217* ..
218* .. Local Scalars ..
219 LOGICAL UPPER
220 INTEGER COUNT, I, J, K, KASE, L, NZ
221 REAL EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
222* ..
223* .. Local Arrays ..
224 INTEGER ISAVE( 3 )
225* ..
226* .. External Subroutines ..
227 EXTERNAL saxpy, scopy, slacn2, spbtrs, ssbmv, xerbla
228* ..
229* .. Intrinsic Functions ..
230 INTRINSIC abs, max, min
231* ..
232* .. External Functions ..
233 LOGICAL LSAME
234 REAL SLAMCH
235 EXTERNAL lsame, slamch
236* ..
237* .. Executable Statements ..
238*
239* Test the input parameters.
240*
241 info = 0
242 upper = lsame( uplo, 'U' )
243 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
244 info = -1
245 ELSE IF( n.LT.0 ) THEN
246 info = -2
247 ELSE IF( kd.LT.0 ) THEN
248 info = -3
249 ELSE IF( nrhs.LT.0 ) THEN
250 info = -4
251 ELSE IF( ldab.LT.kd+1 ) THEN
252 info = -6
253 ELSE IF( ldafb.LT.kd+1 ) THEN
254 info = -8
255 ELSE IF( ldb.LT.max( 1, n ) ) THEN
256 info = -10
257 ELSE IF( ldx.LT.max( 1, n ) ) THEN
258 info = -12
259 END IF
260 IF( info.NE.0 ) THEN
261 CALL xerbla( 'SPBRFS', -info )
262 RETURN
263 END IF
264*
265* Quick return if possible
266*
267 IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
268 DO 10 j = 1, nrhs
269 ferr( j ) = zero
270 berr( j ) = zero
271 10 CONTINUE
272 RETURN
273 END IF
274*
275* NZ = maximum number of nonzero elements in each row of A, plus 1
276*
277 nz = min( n+1, 2*kd+2 )
278 eps = slamch( 'Epsilon' )
279 safmin = slamch( 'Safe minimum' )
280 safe1 = nz*safmin
281 safe2 = safe1 / eps
282*
283* Do for each right hand side
284*
285 DO 140 j = 1, nrhs
286*
287 count = 1
288 lstres = three
289 20 CONTINUE
290*
291* Loop until stopping criterion is satisfied.
292*
293* Compute residual R = B - A * X
294*
295 CALL scopy( n, b( 1, j ), 1, work( n+1 ), 1 )
296 CALL ssbmv( uplo, n, kd, -one, ab, ldab, x( 1, j ), 1, one,
297 $ work( n+1 ), 1 )
298*
299* Compute componentwise relative backward error from formula
300*
301* max(i) ( abs(R(i)) / ( abs(A)*abs(X) + abs(B) )(i) )
302*
303* where abs(Z) is the componentwise absolute value of the matrix
304* or vector Z. If the i-th component of the denominator is less
305* than SAFE2, then SAFE1 is added to the i-th components of the
306* numerator and denominator before dividing.
307*
308 DO 30 i = 1, n
309 work( i ) = abs( b( i, j ) )
310 30 CONTINUE
311*
312* Compute abs(A)*abs(X) + abs(B).
313*
314 IF( upper ) THEN
315 DO 50 k = 1, n
316 s = zero
317 xk = abs( x( k, j ) )
318 l = kd + 1 - k
319 DO 40 i = max( 1, k-kd ), k - 1
320 work( i ) = work( i ) + abs( ab( l+i, k ) )*xk
321 s = s + abs( ab( l+i, k ) )*abs( x( i, j ) )
322 40 CONTINUE
323 work( k ) = work( k ) + abs( ab( kd+1, k ) )*xk + s
324 50 CONTINUE
325 ELSE
326 DO 70 k = 1, n
327 s = zero
328 xk = abs( x( k, j ) )
329 work( k ) = work( k ) + abs( ab( 1, k ) )*xk
330 l = 1 - k
331 DO 60 i = k + 1, min( n, k+kd )
332 work( i ) = work( i ) + abs( ab( l+i, k ) )*xk
333 s = s + abs( ab( l+i, k ) )*abs( x( i, j ) )
334 60 CONTINUE
335 work( k ) = work( k ) + s
336 70 CONTINUE
337 END IF
338 s = zero
339 DO 80 i = 1, n
340 IF( work( i ).GT.safe2 ) THEN
341 s = max( s, abs( work( n+i ) ) / work( i ) )
342 ELSE
343 s = max( s, ( abs( work( n+i ) )+safe1 ) /
344 $ ( work( i )+safe1 ) )
345 END IF
346 80 CONTINUE
347 berr( j ) = s
348*
349* Test stopping criterion. Continue iterating if
350* 1) The residual BERR(J) is larger than machine epsilon, and
351* 2) BERR(J) decreased by at least a factor of 2 during the
352* last iteration, and
353* 3) At most ITMAX iterations tried.
354*
355 IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
356 $ count.LE.itmax ) THEN
357*
358* Update solution and try again.
359*
360 CALL spbtrs( uplo, n, kd, 1, afb, ldafb, work( n+1 ), n,
361 $ info )
362 CALL saxpy( n, one, work( n+1 ), 1, x( 1, j ), 1 )
363 lstres = berr( j )
364 count = count + 1
365 GO TO 20
366 END IF
367*
368* Bound error from formula
369*
370* norm(X - XTRUE) / norm(X) .le. FERR =
371* norm( abs(inv(A))*
372* ( abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) ))) / norm(X)
373*
374* where
375* norm(Z) is the magnitude of the largest component of Z
376* inv(A) is the inverse of A
377* abs(Z) is the componentwise absolute value of the matrix or
378* vector Z
379* NZ is the maximum number of nonzeros in any row of A, plus 1
380* EPS is machine epsilon
381*
382* The i-th component of abs(R)+NZ*EPS*(abs(A)*abs(X)+abs(B))
383* is incremented by SAFE1 if the i-th component of
384* abs(A)*abs(X) + abs(B) is less than SAFE2.
385*
386* Use SLACN2 to estimate the infinity-norm of the matrix
387* inv(A) * diag(W),
388* where W = abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) )))
389*
390 DO 90 i = 1, n
391 IF( work( i ).GT.safe2 ) THEN
392 work( i ) = abs( work( n+i ) ) + nz*eps*work( i )
393 ELSE
394 work( i ) = abs( work( n+i ) ) + nz*eps*work( i ) + safe1
395 END IF
396 90 CONTINUE
397*
398 kase = 0
399 100 CONTINUE
400 CALL slacn2( n, work( 2*n+1 ), work( n+1 ), iwork, ferr( j ),
401 $ kase, isave )
402 IF( kase.NE.0 ) THEN
403 IF( kase.EQ.1 ) THEN
404*
405* Multiply by diag(W)*inv(A**T).
406*
407 CALL spbtrs( uplo, n, kd, 1, afb, ldafb, work( n+1 ), n,
408 $ info )
409 DO 110 i = 1, n
410 work( n+i ) = work( n+i )*work( i )
411 110 CONTINUE
412 ELSE IF( kase.EQ.2 ) THEN
413*
414* Multiply by inv(A)*diag(W).
415*
416 DO 120 i = 1, n
417 work( n+i ) = work( n+i )*work( i )
418 120 CONTINUE
419 CALL spbtrs( uplo, n, kd, 1, afb, ldafb, work( n+1 ), n,
420 $ info )
421 END IF
422 GO TO 100
423 END IF
424*
425* Normalize error.
426*
427 lstres = zero
428 DO 130 i = 1, n
429 lstres = max( lstres, abs( x( i, j ) ) )
430 130 CONTINUE
431 IF( lstres.NE.zero )
432 $ ferr( j ) = ferr( j ) / lstres
433*
434 140 CONTINUE
435*
436 RETURN
437*
438* End of SPBRFS
439*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
saxpy
subroutine saxpy(n, sa, sx, incx, sy, incy)
SAXPY
Definition saxpy.f:89
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
ssbmv
subroutine ssbmv(uplo, n, k, alpha, a, lda, x, incx, beta, y, incy)
SSBMV
Definition ssbmv.f:184
slacn2
subroutine slacn2(n, v, x, isgn, est, kase, isave)
SLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition slacn2.f:136
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
spbtrs
subroutine spbtrs(uplo, n, kd, nrhs, ab, ldab, b, ldb, info)
SPBTRS
Definition spbtrs.f:121
Here is the call graph for this function:
Here is the caller graph for this function: