177 SUBROUTINE zggbal( JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE,
178 $ rscale, work, info )
187 INTEGER IHI, ILO, INFO, LDA, LDB, N
190 DOUBLE PRECISION LSCALE( * ), RSCALE( * ), WORK( * )
191 COMPLEX*16 A( lda, * ), B( ldb, * )
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 )
202 parameter ( czero = ( 0.0d+0, 0.0d+0 ) )
205 INTEGER I, ICAB, IFLOW, IP1, IR, IRAB, IT, J, JC, JP1,
206 $ k, kount, l, lcab, lm1, lrab, lsfmax, lsfmin,
208 DOUBLE PRECISION ALPHA, BASL, BETA, CAB, CMAX, COEF, COEF2,
209 $ coef5, cor, ew, ewc, gamma, pgamma, rab, sfmax,
210 $ sfmin, sum, t, ta, tb, tc
216 DOUBLE PRECISION DDOT, DLAMCH
217 EXTERNAL lsame, izamax, ddot, dlamch
223 INTRINSIC abs, dble, dimag, int, log10, max, min, sign
226 DOUBLE PRECISION CABS1
229 cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
236 IF( .NOT.lsame( job,
'N' ) .AND. .NOT.lsame( job,
'P' ) .AND.
237 $ .NOT.lsame( job,
'S' ) .AND. .NOT.lsame( job,
'B' ) )
THEN
239 ELSE IF( n.LT.0 )
THEN
241 ELSE IF( lda.LT.max( 1, n ) )
THEN
243 ELSE IF( ldb.LT.max( 1, n ) )
THEN
247 CALL xerbla(
'ZGGBAL', -info )
267 IF( lsame( job,
'N' ) )
THEN
279 IF( lsame( job,
'S' ) )
302 IF( a( i, j ).NE.czero .OR. b( i, j ).NE.czero )
310 IF( a( i, j ).NE.czero .OR. b( i, j ).NE.czero )
331 IF( a( i, j ).NE.czero .OR. b( i, j ).NE.czero )
338 IF( a( i, j ).NE.czero .OR. b( i, j ).NE.czero )
355 CALL zswap( n-k+1, a( i, k ), lda, a( m, k ), lda )
356 CALL zswap( n-k+1, b( i, k ), ldb, b( m, k ), ldb )
364 CALL zswap( l, a( 1, j ), 1, a( 1, m ), 1 )
365 CALL zswap( l, b( 1, j ), 1, b( 1, m ), 1 )
368 GO TO ( 20, 90 )iflow
374 IF( lsame( job,
'P' ) )
THEN
402 basl = log10( sclfac )
405 IF( a( i, j ).EQ.czero )
THEN
409 ta = log10( cabs1( a( i, j ) ) ) / basl
412 IF( b( i, j ).EQ.czero )
THEN
416 tb = log10( cabs1( b( i, j ) ) ) / basl
419 work( i+4*n ) = work( i+4*n ) - ta - tb
420 work( j+5*n ) = work( j+5*n ) - ta - tb
424 coef = one / dble( 2*nr )
435 gamma = ddot( nr, work( ilo+4*n ), 1, work( ilo+4*n ), 1 ) +
436 $ ddot( nr, work( ilo+5*n ), 1, work( ilo+5*n ), 1 )
441 ew = ew + work( i+4*n )
442 ewc = ewc + work( i+5*n )
445 gamma = coef*gamma - coef2*( ew**2+ewc**2 ) - coef5*( ew-ewc )**2
449 $ beta = gamma / pgamma
450 t = coef5*( ewc-three*ew )
451 tc = coef5*( ew-three*ewc )
453 CALL dscal( nr, beta, work( ilo ), 1 )
454 CALL dscal( nr, beta, work( ilo+n ), 1 )
456 CALL daxpy( nr, coef, work( ilo+4*n ), 1, work( ilo+n ), 1 )
457 CALL daxpy( nr, coef, work( ilo+5*n ), 1, work( ilo ), 1 )
460 work( i ) = work( i ) + tc
461 work( i+n ) = work( i+n ) + t
470 IF( a( i, j ).EQ.czero )
473 sum = sum + work( j )
475 IF( b( i, j ).EQ.czero )
478 sum = sum + work( j )
480 work( i+2*n ) = dble( kount )*work( i+n ) + sum
487 IF( a( i, j ).EQ.czero )
490 sum = sum + work( i+n )
492 IF( b( i, j ).EQ.czero )
495 sum = sum + work( i+n )
497 work( j+3*n ) = dble( kount )*work( j ) + sum
500 sum = ddot( nr, work( ilo+n ), 1, work( ilo+2*n ), 1 ) +
501 $ ddot( nr, work( ilo ), 1, work( ilo+3*n ), 1 )
508 cor = alpha*work( i+n )
509 IF( abs( cor ).GT.cmax )
511 lscale( i ) = lscale( i ) + cor
512 cor = alpha*work( i )
513 IF( abs( cor ).GT.cmax )
515 rscale( i ) = rscale( i ) + cor
520 CALL daxpy( nr, -alpha, work( ilo+2*n ), 1, work( ilo+4*n ), 1 )
521 CALL daxpy( nr, -alpha, work( ilo+3*n ), 1, work( ilo+5*n ), 1 )
531 sfmin = dlamch(
'S' )
533 lsfmin = int( log10( sfmin ) / basl+one )
534 lsfmax = int( log10( sfmax ) / basl )
536 irab = izamax( n-ilo+1, a( i, ilo ), lda )
537 rab = abs( a( i, irab+ilo-1 ) )
538 irab = izamax( n-ilo+1, b( i, ilo ), ldb )
539 rab = max( rab, abs( b( i, irab+ilo-1 ) ) )
540 lrab = int( log10( rab+sfmin ) / basl+one )
541 ir = int(lscale( i ) + sign( half, lscale( i ) ))
542 ir = min( max( ir, lsfmin ), lsfmax, lsfmax-lrab )
543 lscale( i ) = sclfac**ir
544 icab = izamax( ihi, a( 1, i ), 1 )
545 cab = abs( a( icab, i ) )
546 icab = izamax( ihi, b( 1, i ), 1 )
547 cab = max( cab, abs( b( icab, i ) ) )
548 lcab = int( log10( cab+sfmin ) / basl+one )
549 jc = int(rscale( i ) + sign( half, rscale( i ) ))
550 jc = min( max( jc, lsfmin ), lsfmax, lsfmax-lcab )
551 rscale( i ) = sclfac**jc
557 CALL zdscal( n-ilo+1, lscale( i ), a( i, ilo ), lda )
558 CALL zdscal( n-ilo+1, lscale( i ), b( i, ilo ), ldb )
564 CALL zdscal( ihi, rscale( j ), a( 1, j ), 1 )
565 CALL zdscal( ihi, rscale( j ), b( 1, j ), 1 )
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
subroutine zggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
ZGGBAL
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dscal(N, DA, DX, INCX)
DSCAL
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL