177 SUBROUTINE dggbal( JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE,
178 $ rscale, work, info )
187 INTEGER IHI, ILO, INFO, LDA, LDB, N
190 DOUBLE PRECISION A( lda, * ), B( ldb, * ), LSCALE( * ),
191 $ rscale( * ), work( * )
197 DOUBLE PRECISION ZERO, HALF, ONE
198 parameter( zero = 0.0d+0, half = 0.5d+0, one = 1.0d+0 )
199 DOUBLE PRECISION THREE, SCLFAC
200 parameter( three = 3.0d+0, sclfac = 1.0d+1 )
203 INTEGER I, ICAB, IFLOW, IP1, IR, IRAB, IT, J, JC, JP1,
204 $ k, kount, l, lcab, lm1, lrab, lsfmax, lsfmin,
206 DOUBLE PRECISION ALPHA, BASL, BETA, CAB, CMAX, COEF, COEF2,
207 $ coef5, cor, ew, ewc, gamma, pgamma, rab, sfmax,
208 $ sfmin, sum, t, ta, tb, tc
213 DOUBLE PRECISION DDOT, DLAMCH
214 EXTERNAL lsame, idamax, ddot, dlamch
220 INTRINSIC abs, dble, int, log10, max, min, sign
227 IF( .NOT.lsame( job,
'N' ) .AND. .NOT.lsame( job,
'P' ) .AND.
228 $ .NOT.lsame( job,
'S' ) .AND. .NOT.lsame( job,
'B' ) )
THEN
230 ELSE IF( n.LT.0 )
THEN
232 ELSE IF( lda.LT.max( 1, n ) )
THEN
234 ELSE IF( ldb.LT.max( 1, n ) )
THEN
238 CALL
xerbla(
'DGGBAL', -info )
258 IF( lsame( job,
'N' ) )
THEN
270 IF( lsame( job,
'S' ) )
293 IF( a( i, j ).NE.zero .OR. b( i, j ).NE.zero )
301 IF( a( i, j ).NE.zero .OR. b( i, j ).NE.zero )
322 IF( a( i, j ).NE.zero .OR. b( i, j ).NE.zero )
329 IF( a( i, j ).NE.zero .OR. b( i, j ).NE.zero )
346 CALL
dswap( n-k+1, a( i, k ), lda, a( m, k ), lda )
347 CALL
dswap( n-k+1, b( i, k ), ldb, b( m, k ), ldb )
355 CALL
dswap( l, a( 1, j ), 1, a( 1, m ), 1 )
356 CALL
dswap( l, b( 1, j ), 1, b( 1, m ), 1 )
365 IF( lsame( job,
'P' ) )
THEN
393 basl = log10( sclfac )
400 ta = log10( abs( ta ) ) / basl
404 tb = log10( abs( tb ) ) / basl
406 work( i+4*n ) = work( i+4*n ) - ta - tb
407 work( j+5*n ) = work( j+5*n ) - ta - tb
411 coef = one / dble( 2*nr )
422 gamma = ddot( nr, work( ilo+4*n ), 1, work( ilo+4*n ), 1 ) +
423 $ ddot( nr, work( ilo+5*n ), 1, work( ilo+5*n ), 1 )
428 ew = ew + work( i+4*n )
429 ewc = ewc + work( i+5*n )
432 gamma = coef*gamma - coef2*( ew**2+ewc**2 ) - coef5*( ew-ewc )**2
436 $ beta = gamma / pgamma
437 t = coef5*( ewc-three*ew )
438 tc = coef5*( ew-three*ewc )
440 CALL
dscal( nr, beta, work( ilo ), 1 )
441 CALL
dscal( nr, beta, work( ilo+n ), 1 )
443 CALL
daxpy( nr, coef, work( ilo+4*n ), 1, work( ilo+n ), 1 )
444 CALL
daxpy( nr, coef, work( ilo+5*n ), 1, work( ilo ), 1 )
447 work( i ) = work( i ) + tc
448 work( i+n ) = work( i+n ) + t
457 IF( a( i, j ).EQ.zero )
460 sum = sum + work( j )
462 IF( b( i, j ).EQ.zero )
465 sum = sum + work( j )
467 work( i+2*n ) = dble( kount )*work( i+n ) + sum
474 IF( a( i, j ).EQ.zero )
477 sum = sum + work( i+n )
479 IF( b( i, j ).EQ.zero )
482 sum = sum + work( i+n )
484 work( j+3*n ) = dble( kount )*work( j ) + sum
487 sum = ddot( nr, work( ilo+n ), 1, work( ilo+2*n ), 1 ) +
488 $ ddot( nr, work( ilo ), 1, work( ilo+3*n ), 1 )
495 cor = alpha*work( i+n )
496 IF( abs( cor ).GT.cmax )
498 lscale( i ) = lscale( i ) + cor
499 cor = alpha*work( i )
500 IF( abs( cor ).GT.cmax )
502 rscale( i ) = rscale( i ) + cor
507 CALL
daxpy( nr, -alpha, work( ilo+2*n ), 1, work( ilo+4*n ), 1 )
508 CALL
daxpy( nr, -alpha, work( ilo+3*n ), 1, work( ilo+5*n ), 1 )
518 sfmin = dlamch(
'S' )
520 lsfmin = int( log10( sfmin ) / basl+one )
521 lsfmax = int( log10( sfmax ) / basl )
523 irab = idamax( n-ilo+1, a( i, ilo ), lda )
524 rab = abs( a( i, irab+ilo-1 ) )
525 irab = idamax( n-ilo+1, b( i, ilo ), ldb )
526 rab = max( rab, abs( b( i, irab+ilo-1 ) ) )
527 lrab = int( log10( rab+sfmin ) / basl+one )
528 ir = lscale( i ) + sign( half, lscale( i ) )
529 ir = min( max( ir, lsfmin ), lsfmax, lsfmax-lrab )
530 lscale( i ) = sclfac**ir
531 icab = idamax( ihi, a( 1, i ), 1 )
532 cab = abs( a( icab, i ) )
533 icab = idamax( ihi, b( 1, i ), 1 )
534 cab = max( cab, abs( b( icab, i ) ) )
535 lcab = int( log10( cab+sfmin ) / basl+one )
536 jc = rscale( i ) + sign( half, rscale( i ) )
537 jc = min( max( jc, lsfmin ), lsfmax, lsfmax-lcab )
538 rscale( i ) = sclfac**jc
544 CALL
dscal( n-ilo+1, lscale( i ), a( i, ilo ), lda )
545 CALL
dscal( n-ilo+1, lscale( i ), b( i, ilo ), ldb )
551 CALL
dscal( ihi, rscale( j ), a( 1, j ), 1 )
552 CALL
dscal( ihi, rscale( j ), b( 1, j ), 1 )
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dscal(N, DA, DX, INCX)
DSCAL
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
subroutine dggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
DGGBAL