zggbal.f

Go to the documentation of this file.

*> \brief \b ZGGBAL
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at 
*            http://www.netlib.org/lapack/explore-html/ 
*
*> \htmlonly
*> Download ZGGBAL + dependencies 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zggbal.f"> 
*> [TGZ]</a> 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zggbal.f"> 
*> [ZIP]</a> 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zggbal.f"> 
*> [TXT]</a>
*> \endhtmlonly 
*
*  Definition:
*  ===========
*
*       SUBROUTINE ZGGBAL( JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE,
*                          RSCALE, WORK, INFO )
* 
*       .. Scalar Arguments ..
*       CHARACTER          JOB
*       INTEGER            IHI, ILO, INFO, LDA, LDB, N
*       ..
*       .. Array Arguments ..
*       DOUBLE PRECISION   LSCALE( * ), RSCALE( * ), WORK( * )
*       COMPLEX*16         A( LDA, * ), B( LDB, * )
*       ..
*  
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> ZGGBAL balances a pair of general complex matrices (A,B).  This
*> involves, first, permuting A and B by similarity transformations to
*> isolate eigenvalues in the first 1 to ILO$-$1 and last IHI+1 to N
*> elements on the diagonal; and second, applying a diagonal similarity
*> transformation to rows and columns ILO to IHI to make the rows
*> and columns as close in norm as possible. Both steps are optional.
*>
*> Balancing may reduce the 1-norm of the matrices, and improve the
*> accuracy of the computed eigenvalues and/or eigenvectors in the
*> generalized eigenvalue problem A*x = lambda*B*x.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] JOB
*> \verbatim
*>          JOB is CHARACTER*1
*>          Specifies the operations to be performed on A and B:
*>          = 'N':  none:  simply set ILO = 1, IHI = N, LSCALE(I) = 1.0
*>                  and RSCALE(I) = 1.0 for i=1,...,N;
*>          = 'P':  permute only;
*>          = 'S':  scale only;
*>          = 'B':  both permute and scale.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrices A and B.  N >= 0.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*>          A is COMPLEX*16 array, dimension (LDA,N)
*>          On entry, the input matrix A.
*>          On exit, A is overwritten by the balanced matrix.
*>          If JOB = 'N', A is not referenced.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of the array A. LDA >= max(1,N).
*> \endverbatim
*>
*> \param[in,out] B
*> \verbatim
*>          B is COMPLEX*16 array, dimension (LDB,N)
*>          On entry, the input matrix B.
*>          On exit, B is overwritten by the balanced matrix.
*>          If JOB = 'N', B is not referenced.
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*>          LDB is INTEGER
*>          The leading dimension of the array B. LDB >= max(1,N).
*> \endverbatim
*>
*> \param[out] ILO
*> \verbatim
*>          ILO is INTEGER
*> \endverbatim
*>
*> \param[out] IHI
*> \verbatim
*>          IHI is INTEGER
*>          ILO and IHI are set to integers such that on exit
*>          A(i,j) = 0 and B(i,j) = 0 if i > j and
*>          j = 1,...,ILO-1 or i = IHI+1,...,N.
*>          If JOB = 'N' or 'S', ILO = 1 and IHI = N.
*> \endverbatim
*>
*> \param[out] LSCALE
*> \verbatim
*>          LSCALE is DOUBLE PRECISION array, dimension (N)
*>          Details of the permutations and scaling factors applied
*>          to the left side of A and B.  If P(j) is the index of the
*>          row interchanged with row j, and D(j) is the scaling factor
*>          applied to row j, then
*>            LSCALE(j) = P(j)    for J = 1,...,ILO-1
*>                      = D(j)    for J = ILO,...,IHI
*>                      = P(j)    for J = IHI+1,...,N.
*>          The order in which the interchanges are made is N to IHI+1,
*>          then 1 to ILO-1.
*> \endverbatim
*>
*> \param[out] RSCALE
*> \verbatim
*>          RSCALE is DOUBLE PRECISION array, dimension (N)
*>          Details of the permutations and scaling factors applied
*>          to the right side of A and B.  If P(j) is the index of the
*>          column interchanged with column j, and D(j) is the scaling
*>          factor applied to column j, then
*>            RSCALE(j) = P(j)    for J = 1,...,ILO-1
*>                      = D(j)    for J = ILO,...,IHI
*>                      = P(j)    for J = IHI+1,...,N.
*>          The order in which the interchanges are made is N to IHI+1,
*>          then 1 to ILO-1.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is REAL array, dimension (lwork)
*>          lwork must be at least max(1,6*N) when JOB = 'S' or 'B', and
*>          at least 1 when JOB = 'N' or 'P'.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          = 0:  successful exit
*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee 
*> \author Univ. of California Berkeley 
*> \author Univ. of Colorado Denver 
*> \author NAG Ltd. 
*
*> \date November 2011
*
*> \ingroup complex16GBcomputational
*
*> \par Further Details:
*  =====================
*>
*> \verbatim
*>
*>  See R.C. WARD, Balancing the generalized eigenvalue problem,
*>                 SIAM J. Sci. Stat. Comp. 2 (1981), 141-152.
*> \endverbatim
*>
*  =====================================================================
      SUBROUTINE zggbal( JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE,
     $                   rscale, work, info )
*
*  -- LAPACK computational routine (version 3.4.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2011
*
*     .. Scalar Arguments ..
      CHARACTER          job
      INTEGER            ihi, ilo, info, lda, ldb, n
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   lscale( * ), rscale( * ), work( * )
      COMPLEX*16         a( lda, * ), b( ldb, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   zero, half, one
      parameter( zero = 0.0d+0, half = 0.5d+0, one = 1.0d+0 )
      DOUBLE PRECISION   three, sclfac
      parameter( three = 3.0d+0, sclfac = 1.0d+1 )
      COMPLEX*16         czero
      parameter( czero = ( 0.0d+0, 0.0d+0 ) )
*     ..
*     .. Local Scalars ..
      INTEGER            i, icab, iflow, ip1, ir, irab, it, j, jc, jp1,
     $                   k, kount, l, lcab, lm1, lrab, lsfmax, lsfmin,
     $                   m, nr, nrp2
      DOUBLE PRECISION   alpha, basl, beta, cab, cmax, coef, coef2,
     $                   coef5, cor, ew, ewc, gamma, pgamma, rab, sfmax,
     $                   sfmin, sum, t, ta, tb, tc
      COMPLEX*16         cdum
*     ..
*     .. External Functions ..
      LOGICAL            lsame
      INTEGER            izamax
      DOUBLE PRECISION   ddot, dlamch
      EXTERNAL           lsame, izamax, ddot, dlamch
*     ..
*     .. External Subroutines ..
      EXTERNAL           daxpy, dscal, xerbla, zdscal, zswap
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, dble, dimag, int, log10, max, min, sign
*     ..
*     .. Statement Functions ..
      DOUBLE PRECISION   cabs1
*     ..
*     .. Statement Function definitions ..
      cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters
*
      info = 0
      IF( .NOT.lsame( job, 'N' ) .AND. .NOT.lsame( job, 'P' ) .AND.
     $    .NOT.lsame( job, 'S' ) .AND. .NOT.lsame( job, 'B' ) ) THEN
         info = -1
      ELSE IF( n.LT.0 ) THEN
         info = -2
      ELSE IF( lda.LT.max( 1, n ) ) THEN
         info = -4
      ELSE IF( ldb.LT.max( 1, n ) ) THEN
         info = -6
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'ZGGBAL', -info )
         return
      END IF
*
*     Quick return if possible
*
      IF( n.EQ.0 ) THEN
         ilo = 1
         ihi = n
         return
      END IF
*
      IF( n.EQ.1 ) THEN
         ilo = 1
         ihi = n
         lscale( 1 ) = one
         rscale( 1 ) = one
         return
      END IF
*
      IF( lsame( job, 'N' ) ) THEN
         ilo = 1
         ihi = n
         DO 10 i = 1, n
            lscale( i ) = one
            rscale( i ) = one
   10    continue
         return
      END IF
*
      k = 1
      l = n
      IF( lsame( job, 'S' ) )
     $   go to 190
*
      go to 30
*
*     Permute the matrices A and B to isolate the eigenvalues.
*
*     Find row with one nonzero in columns 1 through L
*
   20 continue
      l = lm1
      IF( l.NE.1 )
     $   go to 30
*
      rscale( 1 ) = 1
      lscale( 1 ) = 1
      go to 190
*
   30 continue
      lm1 = l - 1
      DO 80 i = l, 1, -1
         DO 40 j = 1, lm1
            jp1 = j + 1
            IF( a( i, j ).NE.czero .OR. b( i, j ).NE.czero )
     $         go to 50
   40    continue
         j = l
         go to 70
*
   50    continue
         DO 60 j = jp1, l
            IF( a( i, j ).NE.czero .OR. b( i, j ).NE.czero )
     $         go to 80
   60    continue
         j = jp1 - 1
*
   70    continue
         m = l
         iflow = 1
         go to 160
   80 continue
      go to 100
*
*     Find column with one nonzero in rows K through N
*
   90 continue
      k = k + 1
*
  100 continue
      DO 150 j = k, l
         DO 110 i = k, lm1
            ip1 = i + 1
            IF( a( i, j ).NE.czero .OR. b( i, j ).NE.czero )
     $         go to 120
  110    continue
         i = l
         go to 140
  120    continue
         DO 130 i = ip1, l
            IF( a( i, j ).NE.czero .OR. b( i, j ).NE.czero )
     $         go to 150
  130    continue
         i = ip1 - 1
  140    continue
         m = k
         iflow = 2
         go to 160
  150 continue
      go to 190
*
*     Permute rows M and I
*
  160 continue
      lscale( m ) = i
      IF( i.EQ.m )
     $   go to 170
      CALL zswap( n-k+1, a( i, k ), lda, a( m, k ), lda )
      CALL zswap( n-k+1, b( i, k ), ldb, b( m, k ), ldb )
*
*     Permute columns M and J
*
  170 continue
      rscale( m ) = j
      IF( j.EQ.m )
     $   go to 180
      CALL zswap( l, a( 1, j ), 1, a( 1, m ), 1 )
      CALL zswap( l, b( 1, j ), 1, b( 1, m ), 1 )
*
  180 continue
      go to( 20, 90 )iflow
*
  190 continue
      ilo = k
      ihi = l
*
      IF( lsame( job, 'P' ) ) THEN
         DO 195 i = ilo, ihi
            lscale( i ) = one
            rscale( i ) = one
  195    continue
         return
      END IF
*
      IF( ilo.EQ.ihi )
     $   return
*
*     Balance the submatrix in rows ILO to IHI.
*
      nr = ihi - ilo + 1
      DO 200 i = ilo, ihi
         rscale( i ) = zero
         lscale( i ) = zero
*
         work( i ) = zero
         work( i+n ) = zero
         work( i+2*n ) = zero
         work( i+3*n ) = zero
         work( i+4*n ) = zero
         work( i+5*n ) = zero
  200 continue
*
*     Compute right side vector in resulting linear equations
*
      basl = log10( sclfac )
      DO 240 i = ilo, ihi
         DO 230 j = ilo, ihi
            IF( a( i, j ).EQ.czero ) THEN
               ta = zero
               go to 210
            END IF
            ta = log10( cabs1( a( i, j ) ) ) / basl
*
  210       continue
            IF( b( i, j ).EQ.czero ) THEN
               tb = zero
               go to 220
            END IF
            tb = log10( cabs1( b( i, j ) ) ) / basl
*
  220       continue
            work( i+4*n ) = work( i+4*n ) - ta - tb
            work( j+5*n ) = work( j+5*n ) - ta - tb
  230    continue
  240 continue
*
      coef = one / dble( 2*nr )
      coef2 = coef*coef
      coef5 = half*coef2
      nrp2 = nr + 2
      beta = zero
      it = 1
*
*     Start generalized conjugate gradient iteration
*
  250 continue
*
      gamma = ddot( nr, work( ilo+4*n ), 1, work( ilo+4*n ), 1 ) +
     $        ddot( nr, work( ilo+5*n ), 1, work( ilo+5*n ), 1 )
*
      ew = zero
      ewc = zero
      DO 260 i = ilo, ihi
         ew = ew + work( i+4*n )
         ewc = ewc + work( i+5*n )
  260 continue
*
      gamma = coef*gamma - coef2*( ew**2+ewc**2 ) - coef5*( ew-ewc )**2
      IF( gamma.EQ.zero )
     $   go to 350
      IF( it.NE.1 )
     $   beta = gamma / pgamma
      t = coef5*( ewc-three*ew )
      tc = coef5*( ew-three*ewc )
*
      CALL dscal( nr, beta, work( ilo ), 1 )
      CALL dscal( nr, beta, work( ilo+n ), 1 )
*
      CALL daxpy( nr, coef, work( ilo+4*n ), 1, work( ilo+n ), 1 )
      CALL daxpy( nr, coef, work( ilo+5*n ), 1, work( ilo ), 1 )
*
      DO 270 i = ilo, ihi
         work( i ) = work( i ) + tc
         work( i+n ) = work( i+n ) + t
  270 continue
*
*     Apply matrix to vector
*
      DO 300 i = ilo, ihi
         kount = 0
         sum = zero
         DO 290 j = ilo, ihi
            IF( a( i, j ).EQ.czero )
     $         go to 280
            kount = kount + 1
            sum = sum + work( j )
  280       continue
            IF( b( i, j ).EQ.czero )
     $         go to 290
            kount = kount + 1
            sum = sum + work( j )
  290    continue
         work( i+2*n ) = dble( kount )*work( i+n ) + sum
  300 continue
*
      DO 330 j = ilo, ihi
         kount = 0
         sum = zero
         DO 320 i = ilo, ihi
            IF( a( i, j ).EQ.czero )
     $         go to 310
            kount = kount + 1
            sum = sum + work( i+n )
  310       continue
            IF( b( i, j ).EQ.czero )
     $         go to 320
            kount = kount + 1
            sum = sum + work( i+n )
  320    continue
         work( j+3*n ) = dble( kount )*work( j ) + sum
  330 continue
*
      sum = ddot( nr, work( ilo+n ), 1, work( ilo+2*n ), 1 ) +
     $      ddot( nr, work( ilo ), 1, work( ilo+3*n ), 1 )
      alpha = gamma / sum
*
*     Determine correction to current iteration
*
      cmax = zero
      DO 340 i = ilo, ihi
         cor = alpha*work( i+n )
         IF( abs( cor ).GT.cmax )
     $      cmax = abs( cor )
         lscale( i ) = lscale( i ) + cor
         cor = alpha*work( i )
         IF( abs( cor ).GT.cmax )
     $      cmax = abs( cor )
         rscale( i ) = rscale( i ) + cor
  340 continue
      IF( cmax.LT.half )
     $   go to 350
*
      CALL daxpy( nr, -alpha, work( ilo+2*n ), 1, work( ilo+4*n ), 1 )
      CALL daxpy( nr, -alpha, work( ilo+3*n ), 1, work( ilo+5*n ), 1 )
*
      pgamma = gamma
      it = it + 1
      IF( it.LE.nrp2 )
     $   go to 250
*
*     End generalized conjugate gradient iteration
*
  350 continue
      sfmin = dlamch( 'S' )
      sfmax = one / sfmin
      lsfmin = int( log10( sfmin ) / basl+one )
      lsfmax = int( log10( sfmax ) / basl )
      DO 360 i = ilo, ihi
         irab = izamax( n-ilo+1, a( i, ilo ), lda )
         rab = abs( a( i, irab+ilo-1 ) )
         irab = izamax( n-ilo+1, b( i, ilo ), ldb )
         rab = max( rab, abs( b( i, irab+ilo-1 ) ) )
         lrab = int( log10( rab+sfmin ) / basl+one )
         ir = lscale( i ) + sign( half, lscale( i ) )
         ir = min( max( ir, lsfmin ), lsfmax, lsfmax-lrab )
         lscale( i ) = sclfac**ir
         icab = izamax( ihi, a( 1, i ), 1 )
         cab = abs( a( icab, i ) )
         icab = izamax( ihi, b( 1, i ), 1 )
         cab = max( cab, abs( b( icab, i ) ) )
         lcab = int( log10( cab+sfmin ) / basl+one )
         jc = rscale( i ) + sign( half, rscale( i ) )
         jc = min( max( jc, lsfmin ), lsfmax, lsfmax-lcab )
         rscale( i ) = sclfac**jc
  360 continue
*
*     Row scaling of matrices A and B
*
      DO 370 i = ilo, ihi
         CALL zdscal( n-ilo+1, lscale( i ), a( i, ilo ), lda )
         CALL zdscal( n-ilo+1, lscale( i ), b( i, ilo ), ldb )
  370 continue
*
*     Column scaling of matrices A and B
*
      DO 380 j = ilo, ihi
         CALL zdscal( ihi, rscale( j ), a( 1, j ), 1 )
         CALL zdscal( ihi, rscale( j ), b( 1, j ), 1 )
  380 continue
*
      return
*
*     End of ZGGBAL
*
      END