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

◆ sgbrfs()

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

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

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

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