173 SUBROUTINE dggbal( JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE,
174 $ RSCALE, WORK, INFO )
182 INTEGER IHI, ILO, INFO, LDA, LDB, N
185 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), LSCALE( * ),
186 $ rscale( * ), work( * )
192 DOUBLE PRECISION ZERO, HALF, ONE
193 parameter( zero = 0.0d+0, half = 0.5d+0, one = 1.0d+0 )
194 DOUBLE PRECISION THREE, SCLFAC
195 parameter( three = 3.0d+0, sclfac = 1.0d+1 )
198 INTEGER I, ICAB, IFLOW, IP1, IR, IRAB, IT, J, JC, JP1,
199 $ k, kount, l, lcab, lm1, lrab, lsfmax, lsfmin,
201 DOUBLE PRECISION ALPHA, BASL, BETA, CAB, CMAX, COEF, COEF2,
202 $ coef5, cor, ew, ewc, gamma, pgamma, rab, sfmax,
203 $ sfmin, sum, t, ta, tb, tc
208 DOUBLE PRECISION DDOT, DLAMCH
209 EXTERNAL lsame, idamax, ddot, dlamch
215 INTRINSIC abs, dble, int, log10, max, min, sign
222 IF( .NOT.lsame( job,
'N' ) .AND.
223 $ .NOT.lsame( job,
'P' ) .AND.
224 $ .NOT.lsame( job,
'S' ) .AND.
225 $ .NOT.lsame( job,
'B' ) )
THEN
227 ELSE IF( n.LT.0 )
THEN
229 ELSE IF( lda.LT.max( 1, n ) )
THEN
231 ELSE IF( ldb.LT.max( 1, n ) )
THEN
235 CALL xerbla(
'DGGBAL', -info )
255 IF( lsame( job,
'N' ) )
THEN
267 IF( lsame( job,
'S' ) )
290 IF( a( i, j ).NE.zero .OR. b( i, j ).NE.zero )
298 IF( a( i, j ).NE.zero .OR. b( i, j ).NE.zero )
319 IF( a( i, j ).NE.zero .OR. b( i, j ).NE.zero )
326 IF( a( i, j ).NE.zero .OR. b( i, j ).NE.zero )
343 CALL dswap( n-k+1, a( i, k ), lda, a( m, k ), lda )
344 CALL dswap( n-k+1, b( i, k ), ldb, b( m, k ), ldb )
352 CALL dswap( l, a( 1, j ), 1, a( 1, m ), 1 )
353 CALL dswap( l, b( 1, j ), 1, b( 1, m ), 1 )
356 GO TO ( 20, 90 )iflow
362 IF( lsame( job,
'P' ) )
THEN
390 basl = log10( sclfac )
397 ta = log10( abs( ta ) ) / basl
401 tb = log10( abs( tb ) ) / basl
403 work( i+4*n ) = work( i+4*n ) - ta - tb
404 work( j+5*n ) = work( j+5*n ) - ta - tb
408 coef = one / dble( 2*nr )
419 gamma = ddot( nr, work( ilo+4*n ), 1, work( ilo+4*n ), 1 ) +
420 $ ddot( nr, work( ilo+5*n ), 1, work( ilo+5*n ), 1 )
425 ew = ew + work( i+4*n )
426 ewc = ewc + work( i+5*n )
429 gamma = coef*gamma - coef2*( ew**2+ewc**2 ) - coef5*( ew-ewc )**2
433 $ beta = gamma / pgamma
434 t = coef5*( ewc-three*ew )
435 tc = coef5*( ew-three*ewc )
437 CALL dscal( nr, beta, work( ilo ), 1 )
438 CALL dscal( nr, beta, work( ilo+n ), 1 )
440 CALL daxpy( nr, coef, work( ilo+4*n ), 1, work( ilo+n ), 1 )
441 CALL daxpy( nr, coef, work( ilo+5*n ), 1, work( ilo ), 1 )
444 work( i ) = work( i ) + tc
445 work( i+n ) = work( i+n ) + t
454 IF( a( i, j ).EQ.zero )
457 sum = sum + work( j )
459 IF( b( i, j ).EQ.zero )
462 sum = sum + work( j )
464 work( i+2*n ) = dble( kount )*work( i+n ) + sum
471 IF( a( i, j ).EQ.zero )
474 sum = sum + work( i+n )
476 IF( b( i, j ).EQ.zero )
479 sum = sum + work( i+n )
481 work( j+3*n ) = dble( kount )*work( j ) + sum
484 sum = ddot( nr, work( ilo+n ), 1, work( ilo+2*n ), 1 ) +
485 $ ddot( nr, work( ilo ), 1, work( ilo+3*n ), 1 )
492 cor = alpha*work( i+n )
493 IF( abs( cor ).GT.cmax )
495 lscale( i ) = lscale( i ) + cor
496 cor = alpha*work( i )
497 IF( abs( cor ).GT.cmax )
499 rscale( i ) = rscale( i ) + cor
504 CALL daxpy( nr, -alpha, work( ilo+2*n ), 1, work( ilo+4*n ),
506 CALL daxpy( nr, -alpha, work( ilo+3*n ), 1, work( ilo+5*n ),
517 sfmin = dlamch(
'S' )
519 lsfmin = int( log10( sfmin ) / basl+one )
520 lsfmax = int( log10( sfmax ) / basl )
522 irab = idamax( n-ilo+1, a( i, ilo ), lda )
523 rab = abs( a( i, irab+ilo-1 ) )
524 irab = idamax( n-ilo+1, b( i, ilo ), ldb )
525 rab = max( rab, abs( b( i, irab+ilo-1 ) ) )
526 lrab = int( log10( rab+sfmin ) / basl+one )
527 ir = int(lscale( i ) + sign( half, lscale( i ) ))
528 ir = min( max( ir, lsfmin ), lsfmax, lsfmax-lrab )
529 lscale( i ) = sclfac**ir
530 icab = idamax( ihi, a( 1, i ), 1 )
531 cab = abs( a( icab, i ) )
532 icab = idamax( ihi, b( 1, i ), 1 )
533 cab = max( cab, abs( b( icab, i ) ) )
534 lcab = int( log10( cab+sfmin ) / basl+one )
535 jc = int(rscale( i ) + sign( half, rscale( i ) ))
536 jc = min( max( jc, lsfmin ), lsfmax, lsfmax-lcab )
537 rscale( i ) = sclfac**jc
543 CALL dscal( n-ilo+1, lscale( i ), a( i, ilo ), lda )
544 CALL dscal( n-ilo+1, lscale( i ), b( i, ilo ), ldb )
550 CALL dscal( ihi, rscale( j ), a( 1, j ), 1 )
551 CALL dscal( ihi, rscale( j ), b( 1, j ), 1 )