dgebal()

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

DGEBAL

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

Purpose:

 DGEBAL balances a general real 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 DOUBLE PRECISION 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 BALANC.

  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 dgebal.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   A( LDA, * ), SCALE( * )
*     ..
*
*  =====================================================================
*
*     .. 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            IDAMAX
      DOUBLE PRECISION   DLAMCH, DNRM2
      EXTERNAL           disnan, lsame, idamax, dlamch, dnrm2
*     ..
*     .. External Subroutines ..
      EXTERNAL           dscal, dswap, xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, 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( 'DGEBAL', -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. 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 dswap( l, a( 1, i ), 1, a( 1, l ), 1 )
                     CALL dswap( 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. 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 dswap( l, a( 1, j ), 1, a( 1, k ), 1 )
                     CALL dswap( 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 = dnrm2( l-k+1, a( k, i ), 1 )
            r = dnrm2( l-k+1, a( i, k ), lda )
            ica = idamax( l, a( 1, i ), 1 )
            ca = abs( a( ica, i ) )
            ira = idamax( 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( 'DGEBAL', -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 dscal( n-k+1, g, a( i, k ), lda )
            CALL dscal( l, f, a( 1, i ), 1 )
*
         END DO
*
      END DO
*
      ilo = k
      ihi = l
*
      RETURN
*
*     End of DGEBAL
*

Here is the call graph for this function:

Here is the caller graph for this function: