```      SUBROUTINE DGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
*
*  -- LAPACK routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
CHARACTER          JOB
INTEGER            IHI, ILO, INFO, LDA, N
*     ..
*     .. Array Arguments ..
DOUBLE PRECISION   A( LDA, * ), SCALE( * )
*     ..
*
*  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.
*
*  Arguments
*  =========
*
*  JOB     (input) 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.
*
*  N       (input) INTEGER
*          The order of the matrix A.  N >= 0.
*
*  A       (input/output) 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.
*
*  LDA     (input) INTEGER
*          The leading dimension of the array A.  LDA >= max(1,N).
*
*  ILO     (output) INTEGER
*  IHI     (output) 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.
*
*  SCALE   (output) 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.
*
*  INFO    (output) INTEGER
*          = 0:  successful exit.
*          < 0:  if INFO = -i, the i-th argument had an illegal value.
*
*  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
*
*  =====================================================================
*
*     .. 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
INTEGER            I, ICA, IEXC, IRA, J, K, L, M
DOUBLE PRECISION   C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
\$                   SFMIN2
*     ..
*     .. External Functions ..
LOGICAL            LSAME
INTEGER            IDAMAX
DOUBLE PRECISION   DLAMCH
EXTERNAL           LSAME, IDAMAX, DLAMCH
*     ..
*     .. External Subroutines ..
EXTERNAL           DSCAL, DSWAP, XERBLA
*     ..
*     .. Intrinsic Functions ..
INTRINSIC          ABS, MAX, MIN
*     ..
*     .. 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
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DGEBAL', -INFO )
RETURN
END IF
*
K = 1
L = N
*
IF( N.EQ.0 )
\$   GO TO 210
*
IF( LSAME( JOB, 'N' ) ) THEN
DO 10 I = 1, N
SCALE( I ) = ONE
10    CONTINUE
GO TO 210
END IF
*
IF( LSAME( JOB, 'S' ) )
\$   GO TO 120
*
*     Permutation to isolate eigenvalues if possible
*
GO TO 50
*
*     Row and column exchange.
*
20 CONTINUE
SCALE( M ) = J
IF( J.EQ.M )
\$   GO TO 30
*
CALL DSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
CALL DSWAP( N-K+1, A( J, K ), LDA, A( M, K ), LDA )
*
30 CONTINUE
GO TO ( 40, 80 )IEXC
*
*     Search for rows isolating an eigenvalue and push them down.
*
40 CONTINUE
IF( L.EQ.1 )
\$   GO TO 210
L = L - 1
*
50 CONTINUE
DO 70 J = L, 1, -1
*
DO 60 I = 1, L
IF( I.EQ.J )
\$         GO TO 60
IF( A( J, I ).NE.ZERO )
\$         GO TO 70
60    CONTINUE
*
M = L
IEXC = 1
GO TO 20
70 CONTINUE
*
GO TO 90
*
*     Search for columns isolating an eigenvalue and push them left.
*
80 CONTINUE
K = K + 1
*
90 CONTINUE
DO 110 J = K, L
*
DO 100 I = K, L
IF( I.EQ.J )
\$         GO TO 100
IF( A( I, J ).NE.ZERO )
\$         GO TO 110
100    CONTINUE
*
M = K
IEXC = 2
GO TO 20
110 CONTINUE
*
120 CONTINUE
DO 130 I = K, L
SCALE( I ) = ONE
130 CONTINUE
*
IF( LSAME( JOB, 'P' ) )
\$   GO TO 210
*
*     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
140 CONTINUE
NOCONV = .FALSE.
*
DO 200 I = K, L
C = ZERO
R = ZERO
*
DO 150 J = K, L
IF( J.EQ.I )
\$         GO TO 150
C = C + ABS( A( J, I ) )
R = R + ABS( A( I, J ) )
150    CONTINUE
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 )
\$      GO TO 200
G = R / SCLFAC
F = ONE
S = C + R
160    CONTINUE
IF( C.GE.G .OR. MAX( F, C, CA ).GE.SFMAX2 .OR.
\$       MIN( R, G, RA ).LE.SFMIN2 )GO TO 170
F = F*SCLFAC
C = C*SCLFAC
CA = CA*SCLFAC
R = R / SCLFAC
G = G / SCLFAC
RA = RA / SCLFAC
GO TO 160
*
170    CONTINUE
G = C / SCLFAC
180    CONTINUE
IF( G.LT.R .OR. MAX( R, RA ).GE.SFMAX2 .OR.
\$       MIN( F, C, G, CA ).LE.SFMIN2 )GO TO 190
F = F / SCLFAC
C = C / SCLFAC
G = G / SCLFAC
CA = CA / SCLFAC
R = R*SCLFAC
RA = RA*SCLFAC
GO TO 180
*
*        Now balance.
*
190    CONTINUE
IF( ( C+R ).GE.FACTOR*S )
\$      GO TO 200
IF( F.LT.ONE .AND. SCALE( I ).LT.ONE ) THEN
IF( F*SCALE( I ).LE.SFMIN1 )
\$         GO TO 200
END IF
IF( F.GT.ONE .AND. SCALE( I ).GT.ONE ) THEN
IF( SCALE( I ).GE.SFMAX1 / F )
\$         GO TO 200
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 )
*
200 CONTINUE
*
IF( NOCONV )
\$   GO TO 140
*
210 CONTINUE
ILO = K
IHI = L
*
RETURN
*
*     End of DGEBAL
*
END

```