zgebal()

subroutine zgebal	(	character	job,
		integer	n,
		complex*16, dimension( lda, * )	a,
		integer	lda,
		integer	ilo,
		integer	ihi,
		double precision, dimension( * )	scale,
		integer	info )

ZGEBAL

Download ZGEBAL + dependencies [TGZ] [ZIP] [TXT]

Purpose:

!>
!> ZGEBAL balances a general complex matrix A.  This involves, first,
!> permuting A by a similarity transformation 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 matrix, and improve the
!> accuracy of the computed eigenvalues and/or eigenvectors.
!> 

Parameters

[in]	JOB	!> JOB is CHARACTER*1 !> Specifies the operations to be performed on A: !> = 'N': none: simply set ILO = 1, IHI = N, SCALE(I) = 1.0 !> for i = 1,...,N; !> = 'P': permute only; !> = 'S': scale only; !> = 'B': both permute and scale. !>
[in]	N	!> N is INTEGER !> The order of the matrix A. N >= 0. !>
[in,out]	A	!> 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. !> See Further Details. !>
[in]	LDA	!> LDA is INTEGER !> The leading dimension of the array A. LDA >= max(1,N). !>
[out]	ILO	!> ILO is INTEGER !>
[out]	IHI	!> IHI is INTEGER !> ILO and IHI are set to integers such that on exit !> A(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. !>
[out]	SCALE	!> SCALE is DOUBLE PRECISION array, dimension (N) !> Details of the permutations and scaling factors applied to !> A. If P(j) is the index of the row and column interchanged !> with row and column j and D(j) is the scaling factor !> applied to row and column j, then !> SCALE(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. !>
[out]	INFO	!> INFO is INTEGER !> = 0: successful exit. !> < 0: if INFO = -i, the i-th argument had an illegal value. !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details:

!>
!>  The permutations consist of row and column interchanges which put
!>  the matrix in the form
!>
!>             ( T1   X   Y  )
!>     P A P = (  0   B   Z  )
!>             (  0   0   T2 )
!>
!>  where T1 and T2 are upper triangular matrices whose eigenvalues lie
!>  along the diagonal.  The column indices ILO and IHI mark the starting
!>  and ending columns of the submatrix B. Balancing consists of applying
!>  a diagonal similarity transformation inv(D) * B * D to make the
!>  1-norms of each row of B and its corresponding column nearly equal.
!>  The output matrix is
!>
!>     ( T1     X*D          Y    )
!>     (  0  inv(D)*B*D  inv(D)*Z ).
!>     (  0      0           T2   )
!>
!>  Information about the permutations P and the diagonal matrix D is
!>  returned in the vector SCALE.
!>
!>  This subroutine is based on the EISPACK routine CBAL.
!>
!>  Modified by Tzu-Yi Chen, Computer Science Division, University of
!>    California at Berkeley, USA
!>
!>  Refactored by Evert Provoost, Department of Computer Science,
!>    KU Leuven, Belgium
!> 

Definition at line 162 of file zgebal.f.

*
*  -- LAPACK computational routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER          JOB
      INTEGER            IHI, ILO, INFO, LDA, N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   SCALE( * )
      COMPLEX*16         A( LDA, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      parameter( zero = 0.0d+0, one = 1.0d+0 )
      DOUBLE PRECISION   SCLFAC
      parameter( sclfac = 2.0d+0 )
      DOUBLE PRECISION   FACTOR
      parameter( factor = 0.95d+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            NOCONV, CANSWAP
      INTEGER            I, ICA, IRA, J, K, L
      DOUBLE PRECISION   C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
     $                   SFMIN2
*     ..
*     .. External Functions ..
      LOGICAL            DISNAN, LSAME
      INTEGER            IZAMAX
      DOUBLE PRECISION   DLAMCH, DZNRM2
      EXTERNAL           disnan, lsame, izamax, dlamch,
     $                   dznrm2
*     ..
*     .. External Subroutines ..
      EXTERNAL           xerbla, zdscal, zswap
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, dble, dimag, max, min
*
*     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
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'ZGEBAL', -info )
         RETURN
      END IF
*
*     Quick returns.
*
      IF( n.EQ.0 ) THEN
         ilo = 1
         ihi = 0
         RETURN
      END IF
*
      IF( lsame( job, 'N' ) ) THEN
         DO i = 1, n
            scale( i ) = one
         END DO
         ilo = 1
         ihi = n
         RETURN
      END IF
*
*     Permutation to isolate eigenvalues if possible.
*
      k = 1
      l = n
*
      IF( .NOT.lsame( job, 'S' ) ) THEN
*
*        Row and column exchange.
*
         noconv = .true.
         DO WHILE( noconv )
*
*           Search for rows isolating an eigenvalue and push them down.
*
            noconv = .false.
            DO i = l, 1, -1
               canswap = .true.
               DO j = 1, l
                  IF( i.NE.j .AND. ( dble( a( i, j ) ).NE.zero .OR.
     $                dimag( a( i, j ) ).NE.zero ) ) THEN
                     canswap = .false.
                     EXIT
                  END IF
               END DO
*
               IF( canswap ) THEN
                  scale( l ) = i
                  IF( i.NE.l ) THEN
                     CALL zswap( l, a( 1, i ), 1, a( 1, l ), 1 )
                     CALL zswap( n-k+1, a( i, k ), lda, a( l, k ),
     $                           lda )
                  END IF
                  noconv = .true.
*
                  IF( l.EQ.1 ) THEN
                     ilo = 1
                     ihi = 1
                     RETURN
                  END IF
*
                  l = l - 1
               END IF
            END DO
*
         END DO
 
         noconv = .true.
         DO WHILE( noconv )
*
*           Search for columns isolating an eigenvalue and push them left.
*
            noconv = .false.
            DO j = k, l
               canswap = .true.
               DO i = k, l
                  IF( i.NE.j .AND. ( dble( a( i, j ) ).NE.zero .OR.
     $                dimag( a( i, j ) ).NE.zero ) ) THEN
                     canswap = .false.
                     EXIT
                  END IF
               END DO
*
               IF( canswap ) THEN
                  scale( k ) = j
                  IF( j.NE.k ) THEN
                     CALL zswap( l, a( 1, j ), 1, a( 1, k ), 1 )
                     CALL zswap( n-k+1, a( j, k ), lda, a( k, k ),
     $                           lda )
                  END IF
                  noconv = .true.
*
                  k = k + 1
               END IF
            END DO
*
         END DO
*
      END IF
*
*     Initialize SCALE for non-permuted submatrix.
*
      DO i = k, l
         scale( i ) = one
      END DO
*
*     If we only had to permute, we are done.
*
      IF( lsame( job, 'P' ) ) THEN
         ilo = k
         ihi = l
         RETURN
      END IF
*
*     Balance the submatrix in rows K to L.
*
*     Iterative loop for norm reduction.
*
      sfmin1 = dlamch( 'S' ) / dlamch( 'P' )
      sfmax1 = one / sfmin1
      sfmin2 = sfmin1*sclfac
      sfmax2 = one / sfmin2
*
      noconv = .true.
      DO WHILE( noconv )
         noconv = .false.
*
         DO i = k, l
*
            c = dznrm2( l-k+1, a( k, i ), 1 )
            r = dznrm2( l-k+1, a( i, k ), lda )
            ica = izamax( l, a( 1, i ), 1 )
            ca = abs( a( ica, i ) )
            ira = izamax( n-k+1, a( i, k ), lda )
            ra = abs( a( i, ira+k-1 ) )
*
*           Guard against zero C or R due to underflow.
*
            IF( c.EQ.zero .OR. r.EQ.zero ) cycle
*
*           Exit if NaN to avoid infinite loop
*
            IF( disnan( c+ca+r+ra ) ) THEN
               info = -3
               CALL xerbla( 'ZGEBAL', -info )
               RETURN
            END IF
*
            g = r / sclfac
            f = one
            s = c + r
*
            DO WHILE( c.LT.g .AND. max( f, c, ca ).LT.sfmax2 .AND.
     $                min( r, g, ra ).GT.sfmin2 )
               f = f*sclfac
               c = c*sclfac
               ca = ca*sclfac
               r = r / sclfac
               g = g / sclfac
               ra = ra / sclfac
            END DO
*
            g = c / sclfac
*
            DO WHILE( g.GE.r .AND. max( r, ra ).LT.sfmax2 .AND.
     $                min( f, c, g, ca ).GT.sfmin2 )
               f = f / sclfac
               c = c / sclfac
               g = g / sclfac
               ca = ca / sclfac
               r = r*sclfac
               ra = ra*sclfac
            END DO
*
*           Now balance.
*
            IF( ( c+r ).GE.factor*s ) cycle
            IF( f.LT.one .AND. scale( i ).LT.one ) THEN
               IF( f*scale( i ).LE.sfmin1 ) cycle
            END IF
            IF( f.GT.one .AND. scale( i ).GT.one ) THEN
               IF( scale( i ).GE.sfmax1 / f ) cycle
            END IF
            g = one / f
            scale( i ) = scale( i )*f
            noconv = .true.
*
            CALL zdscal( n-k+1, g, a( i, k ), lda )
            CALL zdscal( l, f, a( 1, i ), 1 )
*
         END DO
*
      END DO
*
      ilo = k
      ihi = l
*
      RETURN
*
*     End of ZGEBAL
*

Here is the call graph for this function:

Here is the caller graph for this function: