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 )
359 GO TO ( 20, 90 )iflow
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 = int(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 = int(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 daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
DGGBAL
subroutine dscal(N, DA, DX, INCX)
DSCAL