204 SUBROUTINE dgbrfs( TRANS, N, KL, KU, NRHS, AB, LDAB, AFB, LDAFB,
205 $ ipiv, b, ldb, x, ldx, ferr, berr, work, iwork,
215 INTEGER info, kl, ku, ldab, ldafb, ldb, ldx, n, nrhs
218 INTEGER ipiv( * ), iwork( * )
219 DOUBLE PRECISION ab( ldab, * ), afb( ldafb, * ), b( ldb, * ),
220 $ berr( * ), ferr( * ), work( * ), x( ldx, * )
227 parameter( itmax = 5 )
228 DOUBLE PRECISION zero
229 parameter( zero = 0.0d+0 )
231 parameter( one = 1.0d+0 )
233 parameter( two = 2.0d+0 )
234 DOUBLE PRECISION three
235 parameter( three = 3.0d+0 )
240 INTEGER count, i, j, k, kase, kk, nz
241 DOUBLE PRECISION eps, lstres, s, safe1, safe2, safmin, xk
250 INTRINSIC abs, max, min
262 notran =
lsame( trans,
'N' )
263 IF( .NOT.notran .AND. .NOT.
lsame( trans,
'T' ) .AND. .NOT.
264 $
lsame( trans,
'C' ) )
THEN
266 ELSE IF( n.LT.0 )
THEN
268 ELSE IF( kl.LT.0 )
THEN
270 ELSE IF( ku.LT.0 )
THEN
272 ELSE IF( nrhs.LT.0 )
THEN
274 ELSE IF( ldab.LT.kl+ku+1 )
THEN
276 ELSE IF( ldafb.LT.2*kl+ku+1 )
THEN
278 ELSE IF( ldb.LT.max( 1, n ) )
THEN
280 ELSE IF( ldx.LT.max( 1, n ) )
THEN
284 CALL
xerbla(
'DGBRFS', -info )
290 IF( n.EQ.0 .OR. nrhs.EQ.0 )
THEN
306 nz = min( kl+ku+2, n+1 )
308 safmin =
dlamch(
'Safe minimum' )
325 CALL
dcopy( n, b( 1, j ), 1, work( n+1 ), 1 )
326 CALL
dgbmv( trans, n, n, kl, ku, -one, ab, ldab, x( 1, j ), 1,
327 $ one, work( n+1 ), 1 )
339 work( i ) = abs( b( i, j ) )
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
356 DO 60 i = max( 1, k-ku ), min( n, k+kl )
357 s = s + abs( ab( kk+i, k ) )*abs( x( i, j ) )
359 work( k ) = work( k ) + s
364 IF( work( i ).GT.safe2 )
THEN
365 s = max( s, abs( work( n+i ) ) / work( i ) )
367 s = max( s, ( abs( work( n+i ) )+safe1 ) /
368 $ ( work( i )+safe1 ) )
379 IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
380 $ count.LE.itmax )
THEN
384 CALL
dgbtrs( trans, n, kl, ku, 1, afb, ldafb, ipiv,
385 $ work( n+1 ), n, info )
386 CALL
daxpy( n, one, work( n+1 ), 1, x( 1, j ), 1 )
415 IF( work( i ).GT.safe2 )
THEN
416 work( i ) = abs( work( n+i ) ) + nz*eps*work( i )
418 work( i ) = abs( work( n+i ) ) + nz*eps*work( i ) + safe1
424 CALL
dlacn2( n, work( 2*n+1 ), work( n+1 ), iwork, ferr( j ),
431 CALL
dgbtrs( transt, n, kl, ku, 1, afb, ldafb, ipiv,
432 $ work( n+1 ), n, info )
434 work( n+i ) = work( n+i )*work( i )
441 work( n+i ) = work( n+i )*work( i )
443 CALL
dgbtrs( trans, n, kl, ku, 1, afb, ldafb, ipiv,
444 $ work( n+1 ), n, info )
453 lstres = max( lstres, abs( x( i, j ) ) )
456 $ ferr( j ) = ferr( j ) / lstres