      SUBROUTINE SCALEG (N,MA,A,MB,B,LOW,IGH,CSCALE,CPERM,WK)
C
C     *****PARAMETERS:
      INTEGER IGH,LOW,MA,MB,N
      DOUBLE PRECISION A(MA,N),B(MB,N),CPERM(N),CSCALE(N),WK(N,6)
C
C     *****LOCAL VARIABLES:
      INTEGER I,IR,IT,J,JC,KOUNT,NR,NRP2
      DOUBLE PRECISION ALPHA,BASL,BETA,CMAX,COEF,COEF2,COEF5,COR,
     *                 EW,EWC,FI,FJ,GAMMA,PGAMMA,SUM,T,TA,TB,TC
C
C     *****FORTRAN FUNCTIONS:
      DOUBLE PRECISION DABS, DLOG10, DSIGN
C     FLOAT
C
C     *****SUBROUTINES CALLED:
C     NONE
C
C     ---------------------------------------------------------------
C
C     *****PURPOSE:
C     SCALES THE MATRICES A AND B IN THE GENERALIZED EIGENVALUE
C     PROBLEM A*X = (LAMBDA)*B*X SUCH THAT THE MAGNITUDES OF THE
C     ELEMENTS OF THE SUBMATRICES OF A AND B (AS SPECIFIED BY LOW
C     AND IGH) ARE CLOSE TO UNITY IN THE LEAST SQUARES SENSE.
C     REF.:  WARD, R. C., BALANCING THE GENERALIZED EIGENVALUE
C     PROBLEM, SIAM J. SCI. STAT. COMPUT., VOL. 2, NO. 2, JUNE 1981,
C     141-152.
C
C     *****PARAMETER DESCRIPTION:
C
C     ON INPUT:
C
C       MA,MB   INTEGER
C               ROW DIMENSIONS OF THE ARRAYS CONTAINING MATRICES
C               A AND B RESPECTIVELY, AS DECLARED IN THE MAIN CALLING
C               PROGRAM DIMENSION STATEMENT;
C
C       N       INTEGER
C               ORDER OF THE MATRICES A AND B;
C
C       A       REAL(MA,N)
C               CONTAINS THE A MATRIX OF THE GENERALIZED EIGENPROBLEM
C               DEFINED ABOVE;
C
C       B       REAL(MB,N)
C               CONTAINS THE B MATRIX OF THE GENERALIZED EIGENPROBLEM
C               DEFINED ABOVE;
C
C       LOW     INTEGER
C               SPECIFIES THE BEGINNING INDEX FOR THE ROWS AND
C               COLUMNS OF A AND B TO BE SCALED;
C
C       IGH     INTEGER
C               SPECIFIES THE ENDING INDEX FOR THE ROWS AND COLUMNS
C               OF A AND B TO BE SCALED;
C
C       CPERM   REAL(N)
C               WORK ARRAY.  ONLY LOCATIONS LOW THROUGH IGH ARE
C               REFERENCED AND ALTERED BY THIS SUBROUTINE;
C
C       WK      REAL(N,6)
C               WORK ARRAY THAT MUST CONTAIN AT LEAST 6*N LOCATIONS.
C               ONLY LOCATIONS LOW THROUGH IGH, N+LOW THROUGH N+IGH,
C               ..., 5*N+LOW THROUGH 5*N+IGH ARE REFERENCED AND
C               ALTERED BY THIS SUBROUTINE.
C
C     ON OUTPUT:
C
C       A,B     CONTAIN THE SCALED A AND B MATRICES;
C
C       CSCALE  REAL(N)
C               CONTAINS IN ITS LOW THROUGH IGH LOCATIONS THE INTEGER
C               EXPONENTS OF 2 USED FOR THE COLUMN SCALING FACTORS.
C               THE OTHER LOCATIONS ARE NOT REFERENCED;
C
C       WK      CONTAINS IN ITS LOW THROUGH IGH LOCATIONS THE INTEGER
C               EXPONENTS OF 2 USED FOR THE ROW SCALING FACTORS.
C
C     *****ALGORITHM NOTES:
C     NONE.
C
C     *****HISTORY:
C     WRITTEN BY R. C. WARD.......
C     Modified 8/86 by Bobby Bodenheimer so that if
C       SUM = 0 (corresponding to the case where the matrix
C       doesn't need to be scaled) the routine returns.
C
C     ---------------------------------------------------------------
C
      IF (LOW .EQ. IGH) GO TO 410
      DO 210 I = LOW,IGH
         WK(I,1) = 0.0D0
         WK(I,2) = 0.0D0
         WK(I,3) = 0.0D0
         WK(I,4) = 0.0D0
         WK(I,5) = 0.0D0
         WK(I,6) = 0.0D0
         CSCALE(I) = 0.0D0
         CPERM(I) = 0.0D0
  210 CONTINUE
C
C     COMPUTE RIGHT SIDE VECTOR IN RESULTING LINEAR EQUATIONS
C
      BASL = DLOG10(2.0D0)
      DO 240 I = LOW,IGH
         DO 240 J = LOW,IGH
            TB = B(I,J)
            TA = A(I,J)
            IF (TA .EQ. 0.0D0) GO TO 220
            TA = DLOG10(DABS(TA)) / BASL
  220       CONTINUE
            IF (TB .EQ. 0.0D0) GO TO 230
            TB = DLOG10(DABS(TB)) / BASL
  230       CONTINUE
            WK(I,5) = WK(I,5) - TA - TB
            WK(J,6) = WK(J,6) - TA - TB
  240 CONTINUE
      NR = IGH-LOW+1
      COEF = 1.0D0/FLOAT(2*NR)
      COEF2 = COEF*COEF
      COEF5 = 0.5D0*COEF2
      NRP2 = NR+2
      BETA = 0.0D0
      IT = 1
C
C     START GENERALIZED CONJUGATE GRADIENT ITERATION
C
  250 CONTINUE
      EW = 0.0D0
      EWC = 0.0D0
      GAMMA = 0.0D0
      DO 260 I = LOW,IGH
         GAMMA = GAMMA + WK(I,5)*WK(I,5) + WK(I,6)*WK(I,6)
         EW = EW + WK(I,5)
         EWC = EWC + WK(I,6)
  260 CONTINUE
      GAMMA = COEF*GAMMA - COEF2*(EW**2 + EWC**2)
     +        - COEF5*(EW - EWC)**2
      IF (IT .NE. 1) BETA = GAMMA / PGAMMA
      T = COEF5*(EWC - 3.0D0*EW)
      TC = COEF5*(EW - 3.0D0*EWC)
      DO 270 I = LOW,IGH
         WK(I,2) = BETA*WK(I,2) + COEF*WK(I,5) + T
         CPERM(I) = BETA*CPERM(I) + COEF*WK(I,6) + TC
  270 CONTINUE
C
C     APPLY MATRIX TO VECTOR
C
      DO 300 I = LOW,IGH
         KOUNT = 0
         SUM = 0.0D0
         DO 290 J = LOW,IGH
            IF (A(I,J) .EQ. 0.0D0) GO TO 280
            KOUNT = KOUNT+1
            SUM = SUM + CPERM(J)
  280       CONTINUE
            IF (B(I,J) .EQ. 0.0D0) GO TO 290
            KOUNT = KOUNT+1
            SUM = SUM + CPERM(J)
  290    CONTINUE
         WK(I,3) = FLOAT(KOUNT)*WK(I,2) + SUM
  300 CONTINUE
      DO 330 J = LOW,IGH
         KOUNT = 0
         SUM = 0.0D0
         DO 320 I = LOW,IGH
            IF (A(I,J) .EQ. 0.0D0) GO TO 310
            KOUNT = KOUNT+1
            SUM = SUM + WK(I,2)
  310       CONTINUE
            IF (B(I,J) .EQ. 0.0D0) GO TO 320
            KOUNT = KOUNT+1
            SUM = SUM + WK(I,2)
  320    CONTINUE
         WK(J,4) = FLOAT(KOUNT)*CPERM(J) + SUM
  330 CONTINUE
      SUM = 0.0D0
      DO 340 I = LOW,IGH
         SUM = SUM + WK(I,2)*WK(I,3) + CPERM(I)*WK(I,4)
  340 CONTINUE
      IF(SUM.EQ.0.0D0) RETURN
      ALPHA = GAMMA / SUM
C
C     DETERMINE CORRECTION TO CURRENT ITERATE
C
      CMAX = 0.0D0
      DO 350 I = LOW,IGH
         COR = ALPHA * WK(I,2)
         IF (DABS(COR) .GT. CMAX) CMAX = DABS(COR)
         WK(I,1) = WK(I,1) + COR
         COR = ALPHA * CPERM(I)
         IF (DABS(COR) .GT. CMAX) CMAX = DABS(COR)
         CSCALE(I) = CSCALE(I) + COR
  350 CONTINUE
      IF (CMAX .LT. 0.5D0) GO TO 370
      DO 360 I = LOW,IGH
         WK(I,5) = WK(I,5) - ALPHA*WK(I,3)
         WK(I,6) = WK(I,6) - ALPHA*WK(I,4)
  360 CONTINUE
      PGAMMA = GAMMA
      IT = IT+1
      IF (IT .LE. NRP2) GO TO 250
C
C     END GENERALIZED CONJUGATE GRADIENT ITERATION
C
  370 CONTINUE
      DO 380 I = LOW,IGH
         IR = WK(I,1) + DSIGN(0.5D0,WK(I,1))
         WK(I,1) = IR
         JC = CSCALE(I) + DSIGN(0.5D0,CSCALE(I))
         CSCALE(I) = JC
  380 CONTINUE
C
C     SCALE A AND B
C
      DO 400 I = 1,IGH
         IR = WK(I,1)
         FI = 2.0D0**IR
         IF (I .LT. LOW) FI = 1.0D0
         DO 400 J =LOW,N
            JC = CSCALE(J)
            FJ = 2.0D0**JC
            IF (J .LE. IGH) GO TO 390
            IF (I .LT. LOW) GO TO 400
            FJ = 1.0D0
  390       CONTINUE
            A(I,J) = A(I,J)*FI*FJ
            B(I,J) = B(I,J)*FI*FJ
  400 CONTINUE
  410 CONTINUE
      RETURN
C
C     LAST LINE OF SCALEG
C
      END
