 LAPACK  3.6.1 LAPACK: Linear Algebra PACKage
 subroutine sgbrfs ( character TRANS, integer N, integer KL, integer KU, integer NRHS, real, dimension( ldab, * ) AB, integer LDAB, real, dimension( ldafb, * ) AFB, integer LDAFB, integer, dimension( * ) IPIV, 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 )

SGBRFS

Purpose:
``` SGBRFS 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 = 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 REAL 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 REAL array, dimension (LDAFB,N) Details of the LU factorization of the band matrix A, as computed by SGBTRF. 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 SGBTRF; for 1<=i<=N, row i of the matrix was interchanged with row IPIV(i).``` [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 SGBTRS. 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.`
Date
November 2011

Definition at line 207 of file sgbrfs.f.

207 *
208 * -- LAPACK computational routine (version 3.4.0) --
209 * -- LAPACK is a software package provided by Univ. of Tennessee, --
210 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
211 * November 2011
212 *
213 * .. Scalar Arguments ..
214  CHARACTER trans
215  INTEGER info, kl, ku, ldab, ldafb, ldb, ldx, n, nrhs
216 * ..
217 * .. Array Arguments ..
218  INTEGER ipiv( * ), iwork( * )
219  REAL ab( ldab, * ), afb( ldafb, * ), b( ldb, * ),
220  \$ berr( * ), ferr( * ), work( * ), x( ldx, * )
221 * ..
222 *
223 * =====================================================================
224 *
225 * .. Parameters ..
226  INTEGER itmax
227  parameter ( itmax = 5 )
228  REAL zero
229  parameter ( zero = 0.0e+0 )
230  REAL one
231  parameter ( one = 1.0e+0 )
232  REAL two
233  parameter ( two = 2.0e+0 )
234  REAL three
235  parameter ( three = 3.0e+0 )
236 * ..
237 * .. Local Scalars ..
238  LOGICAL notran
239  CHARACTER transt
240  INTEGER count, i, j, k, kase, kk, nz
241  REAL eps, lstres, s, safe1, safe2, safmin, xk
242 * ..
243 * .. Local Arrays ..
244  INTEGER isave( 3 )
245 * ..
246 * .. External Subroutines ..
247  EXTERNAL saxpy, scopy, sgbmv, sgbtrs, slacn2, xerbla
248 * ..
249 * .. Intrinsic Functions ..
250  INTRINSIC abs, max, min
251 * ..
252 * .. External Functions ..
253  LOGICAL lsame
254  REAL slamch
255  EXTERNAL lsame, slamch
256 * ..
257 * .. Executable Statements ..
258 *
259 * Test the input parameters.
260 *
261  info = 0
262  notran = lsame( trans, 'N' )
263  IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
264  \$ lsame( trans, 'C' ) ) THEN
265  info = -1
266  ELSE IF( n.LT.0 ) THEN
267  info = -2
268  ELSE IF( kl.LT.0 ) THEN
269  info = -3
270  ELSE IF( ku.LT.0 ) THEN
271  info = -4
272  ELSE IF( nrhs.LT.0 ) THEN
273  info = -5
274  ELSE IF( ldab.LT.kl+ku+1 ) THEN
275  info = -7
276  ELSE IF( ldafb.LT.2*kl+ku+1 ) THEN
277  info = -9
278  ELSE IF( ldb.LT.max( 1, n ) ) THEN
279  info = -12
280  ELSE IF( ldx.LT.max( 1, n ) ) THEN
281  info = -14
282  END IF
283  IF( info.NE.0 ) THEN
284  CALL xerbla( 'SGBRFS', -info )
285  RETURN
286  END IF
287 *
288 * Quick return if possible
289 *
290  IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
291  DO 10 j = 1, nrhs
292  ferr( j ) = zero
293  berr( j ) = zero
294  10 CONTINUE
295  RETURN
296  END IF
297 *
298  IF( notran ) THEN
299  transt = 'T'
300  ELSE
301  transt = 'N'
302  END IF
303 *
304 * NZ = maximum number of nonzero elements in each row of A, plus 1
305 *
306  nz = min( kl+ku+2, n+1 )
307  eps = slamch( 'Epsilon' )
308  safmin = slamch( 'Safe minimum' )
309  safe1 = nz*safmin
310  safe2 = safe1 / eps
311 *
312 * Do for each right hand side
313 *
314  DO 140 j = 1, nrhs
315 *
316  count = 1
317  lstres = three
318  20 CONTINUE
319 *
320 * Loop until stopping criterion is satisfied.
321 *
322 * Compute residual R = B - op(A) * X,
323 * where op(A) = A, A**T, or A**H, depending on TRANS.
324 *
325  CALL scopy( n, b( 1, j ), 1, work( n+1 ), 1 )
326  CALL sgbmv( trans, n, n, kl, ku, -one, ab, ldab, x( 1, j ), 1,
327  \$ one, work( n+1 ), 1 )
328 *
329 * Compute componentwise relative backward error from formula
330 *
331 * max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
332 *
333 * where abs(Z) is the componentwise absolute value of the matrix
334 * or vector Z. If the i-th component of the denominator is less
335 * than SAFE2, then SAFE1 is added to the i-th components of the
336 * numerator and denominator before dividing.
337 *
338  DO 30 i = 1, n
339  work( i ) = abs( b( i, j ) )
340  30 CONTINUE
341 *
342 * Compute abs(op(A))*abs(X) + abs(B).
343 *
344  IF( notran ) THEN
345  DO 50 k = 1, n
346  kk = ku + 1 - k
347  xk = abs( x( k, j ) )
348  DO 40 i = max( 1, k-ku ), min( n, k+kl )
349  work( i ) = work( i ) + abs( ab( kk+i, k ) )*xk
350  40 CONTINUE
351  50 CONTINUE
352  ELSE
353  DO 70 k = 1, n
354  s = zero
355  kk = ku + 1 - k
356  DO 60 i = max( 1, k-ku ), min( n, k+kl )
357  s = s + abs( ab( kk+i, k ) )*abs( x( i, j ) )
358  60 CONTINUE
359  work( k ) = work( k ) + s
360  70 CONTINUE
361  END IF
362  s = zero
363  DO 80 i = 1, n
364  IF( work( i ).GT.safe2 ) THEN
365  s = max( s, abs( work( n+i ) ) / work( i ) )
366  ELSE
367  s = max( s, ( abs( work( n+i ) )+safe1 ) /
368  \$ ( work( i )+safe1 ) )
369  END IF
370  80 CONTINUE
371  berr( j ) = s
372 *
373 * Test stopping criterion. Continue iterating if
374 * 1) The residual BERR(J) is larger than machine epsilon, and
375 * 2) BERR(J) decreased by at least a factor of 2 during the
376 * last iteration, and
377 * 3) At most ITMAX iterations tried.
378 *
379  IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
380  \$ count.LE.itmax ) THEN
381 *
382 * Update solution and try again.
383 *
384  CALL sgbtrs( trans, n, kl, ku, 1, afb, ldafb, ipiv,
385  \$ work( n+1 ), n, info )
386  CALL saxpy( n, one, work( n+1 ), 1, x( 1, j ), 1 )
387  lstres = berr( j )
388  count = count + 1
389  GO TO 20
390  END IF
391 *
392 * Bound error from formula
393 *
394 * norm(X - XTRUE) / norm(X) .le. FERR =
395 * norm( abs(inv(op(A)))*
396 * ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
397 *
398 * where
399 * norm(Z) is the magnitude of the largest component of Z
400 * inv(op(A)) is the inverse of op(A)
401 * abs(Z) is the componentwise absolute value of the matrix or
402 * vector Z
403 * NZ is the maximum number of nonzeros in any row of A, plus 1
404 * EPS is machine epsilon
405 *
406 * The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
407 * is incremented by SAFE1 if the i-th component of
408 * abs(op(A))*abs(X) + abs(B) is less than SAFE2.
409 *
410 * Use SLACN2 to estimate the infinity-norm of the matrix
411 * inv(op(A)) * diag(W),
412 * where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
413 *
414  DO 90 i = 1, n
415  IF( work( i ).GT.safe2 ) THEN
416  work( i ) = abs( work( n+i ) ) + nz*eps*work( i )
417  ELSE
418  work( i ) = abs( work( n+i ) ) + nz*eps*work( i ) + safe1
419  END IF
420  90 CONTINUE
421 *
422  kase = 0
423  100 CONTINUE
424  CALL slacn2( n, work( 2*n+1 ), work( n+1 ), iwork, ferr( j ),
425  \$ kase, isave )
426  IF( kase.NE.0 ) THEN
427  IF( kase.EQ.1 ) THEN
428 *
429 * Multiply by diag(W)*inv(op(A)**T).
430 *
431  CALL sgbtrs( transt, n, kl, ku, 1, afb, ldafb, ipiv,
432  \$ work( n+1 ), n, info )
433  DO 110 i = 1, n
434  work( n+i ) = work( n+i )*work( i )
435  110 CONTINUE
436  ELSE
437 *
438 * Multiply by inv(op(A))*diag(W).
439 *
440  DO 120 i = 1, n
441  work( n+i ) = work( n+i )*work( i )
442  120 CONTINUE
443  CALL sgbtrs( trans, n, kl, ku, 1, afb, ldafb, ipiv,
444  \$ work( n+1 ), n, info )
445  END IF
446  GO TO 100
447  END IF
448 *
449 * Normalize error.
450 *
451  lstres = zero
452  DO 130 i = 1, n
453  lstres = max( lstres, abs( x( i, j ) ) )
454  130 CONTINUE
455  IF( lstres.NE.zero )
456  \$ ferr( j ) = ferr( j ) / lstres
457 *
458  140 CONTINUE
459 *
460  RETURN
461 *
462 * End of SGBRFS
463 *
subroutine sgbmv(TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGBMV
Definition: sgbmv.f:187
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
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:138
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:54
subroutine sgbtrs(TRANS, N, KL, KU, NRHS, AB, LDAB, IPIV, B, LDB, INFO)
SGBTRS
Definition: sgbtrs.f:140
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53

Here is the call graph for this function:

Here is the caller graph for this function: