SUBROUTINE SCLGNP(N,NN,MMAXT,NUMT,DEG,MODE,EPS0,COEF,
     $                    NNUMT,DDEG,CCOEF,ALPHA,BETA,RWORK,XWORK,
     $                                      FACV,FACE,COESCL,IERR)
C
C SCLGNP  SCALES THE COEFFICIENTS OF A POLYNOMIAL SYSTEM OF N
C EQUATIONS IN N UNKNOWNS, F(X)=0, WHERE THE JTH TERM OF
C THE ITH EQUATION LOOKS LIKE:
C
C    COEF(I,J) * X(1)**DEG(I,1,J) ... X(N)**DEG(I,N,J)
C
C THE ITH EQUATION IS SCALED BY 10**FACE(I).  THE KTH
C VARIABLE IS SCALED BY 10**FACV(K).  IN OTHER WORDS, X(K) =
C 10**FACV(K) * Y(K), WHERE Y SOLVES THE SCALED EQUATION.
C THE SCALED EQUATION HAS THE SAME FORM AS THE ORIGINAL
C EQUATION, EXCEPT THAT COESCL(I,J) REPLACES COEF(I,J), WHERE
C
C COESCL(I,J)=COEF(I,J)* 10**( FACE(I) + FACV(1)*DEG(I,1,J)+ ...
C                                       +FACV(N)*DEG(I,N,J) )
C
C THE CRITERION FOR GENERATING FACE AND FACV IS THAT OF
C MINIMIZING THE SUM OF SQUARES OF THE EXPONENTS OF THE SCALED
C COEFFICIENTS.  IT TURNS OUT THAT THIS CRITERION REDUCES TO
C SOLVING A SINGLE LINEAR SYSTEM, ALPHA*X = BETA, AS DEFINED
C IN THE CODE BELOW.  FURTHER, THE FORM OF THE POLYNOMIAL
C SYSTEM ALONE DETERMINES THE MATRIX ALPHA.  THUS, IN CASES
C IN WHICH MANY SYSTEMS OF THE SAME FORM, BUT WITH DIFFERENT
C COEFFICIENTS, ARE TO BE SCALED, THE MATRIX ALPHA IS
C UNCHANGED AND MAY BE FACTORED ONLY ONCE (BY  QRFAQF).  WHEN
C SCLGNP  IS CALLED WITH MODE=1,  SCLGNP  DOES NOT RECOMPUTE OR
C REFACTOR THE MATRIX ALPHA.  SEE MEINTJES AND MORGAN "A
C METHODOLOGY FOR SOLVING CHEMICAL EQUILIBRIUM SYSTEMS"
C (GENERAL MOTORS RESEARCH LABORATORIES TECHNICAL REPORT
C GMR-4971).
C
C SUBROUTINES CALLED DIRECTLY:  QRFAQF, QRSLQF.
C SUBROUTINES CALLED INDIRECTLY:  DAXPY, DCOPY, DDOT, DNRM2, DSCAL.
C
C ON INPUT:
C
C N  IS THE NUMBER OF EQUATIONS AND THE NUMBER OF VARIABLES.
C
C NN  IS THE DECLARED DIMENSION OF SEVERAL ARRAY INDICES.
C
C MMAXT  IS AN UPPER BOUND ON THE SET NUMT(I), I=1 TO N.
C
C NUMT(I)  IS THE NUMBER OF TERMS IN THE I-TH EQUATION FOR I=1 TO N.
C
C DEG(I,K,J)  IS THE DEGREE OF THE K-TH VARIABLE IN THE
C   J-TH TERM OF THE I-TH EQUATION FOR I=1 TO N, J=1 TO NUMT(I), AND
C   K=1 TO N.
C
C MODE
C  =1  THIS IS NOT THE FIRST CALL TO  SCLGNP, AND THE FORM OF THE
C      SYSTEM HAS NOT CHANGED.
C  =0  THIS IS THE FIRST CALL TO  SCLGNP.
C
C EPS0  ZERO-EPSILON FOR TERMS (TERMS LESS THAN  EPS0  IN MAGNITUDE
C   ARE TREATED AS ZERO BY THE SCALING ALGORITHM).
C
C COEF(I,J)  IS THE COEFFICIENT OF THE JTH TERM OF THE ITH EQUATION
C   FOR I=1 TO N AND J=1 TO NUMT(N).  (COEF(I,J) MAY BE ZERO.)
C
C NNUMT, DDEG, CCOEF, ALPHA, BETA, RWORK, AND  XWORK  ARE WORKSPACES.
C
C ON OUTPUT:
C
C N, NUMT, DEG, MODE, EPS0, AND  COEF  ARE UNCHANGED.
C
C FACV(I)  IS THE VARIABLE SCALE FACTOR FOR THE I-TH VARIABLE, FOR
C   I=1 TO N.
C
C FACE(I)  IS THE EQUATION SCALE FACTOR FOR THE I-TH EQUATION, FOR
C   I=1 TO N.
C
C COESCL(I,J)  IS THE SCALED VERSION OF COEFFICIENT  COEF(I,J), FOR
C   I=1 TO N, J=1 TO NUMT(I), UNLESS IERR=1.
C
C IERR
C   =0  IF SCALING MATRIX, ALPHA, IS WELL CONDITIONED.
C   =1  OTHERWISE.  IN THIS CASE, ALPHA IS "REPAIRED" AND A
C         SCALING IS COMPUTED.
C
C
C DECLARATION OF INPUT
      INTEGER N,NN,MMAXT,NUMT(NN),DEG(NN,NN+1,MMAXT)
      INTEGER MODE
      DOUBLE PRECISION EPS0,COEF
      DIMENSION COEF(NN,MMAXT)
C
C DECLARATION OF WORKSPACE
      INTEGER NNUMT,DDEG
      DOUBLE PRECISION CCOEF,ALPHA,BETA,RWORK,XWORK
      DIMENSION NNUMT(N),DDEG(N,N+1,MMAXT)
      DIMENSION CCOEF(N,MMAXT),ALPHA(2*N,2*N),BETA(2*N),
     $ RWORK(N*(2*N+1)),XWORK(2*N)
C
C DECLARATION OF OUTPUT
      INTEGER IERR
      DOUBLE PRECISION FACV,FACE,COESCL
      DIMENSION FACV(N),FACE(N),COESCL(N,MMAXT)
C
C DECLARATION OF VARIABLES
      INTEGER I,IDAMAX,IFLAG,INDEX,IRMAX,J,JJ,K,LENR,N2,S
      DOUBLE PRECISION D1MACH,DUM,LMFPN,NTUR,RTOL,SUM
C
      SAVE
C
      IERR=0
      N2=2*N
      LMFPN=D1MACH(2)
      NTUR=D1MACH(4)*N
      LENR=N*(N+1)/2
C
C  DELETE NEAR ZERO TERMS
      DO  60 I=1,N
         JJ=0
         NNUMT(I)=0
         DO 40 J=1,NUMT(I)
             IF(ABS(COEF(I,J)) .GT. EPS0) THEN
                JJ=JJ+1
                NNUMT(I)=NNUMT(I)+1
                  CCOEF(I,JJ)=COEF(I,J)
                  DO 20 K=1,N
                    DDEG(I,K,JJ)=DEG(I,K,J)
  20              CONTINUE
             END IF
  40     CONTINUE
  60  CONTINUE
      DO 90 I=1,N
          DO 80 J=1,NNUMT(I)
              COESCL(I,J)=LOG10(ABS(CCOEF(I,J)))
  80      CONTINUE
  90  CONTINUE
C
C SKIP OVER THE GENERATION AND DECOMPOSITON OF MATRIX ALPHA IF MODE=1
      IF (MODE .EQ. 0) THEN
C
C GENERATE THE MATRIX ALPHA
      DO 110 S=1,N
      DO 110 K=1,N
          ALPHA(S,K)=0
 110  CONTINUE
      DO 200 S=1,N
        ALPHA(S,S)=NNUMT(S)
 200  CONTINUE
      DO 300 S=1,N
        DO 300 I=1,N
          SUM=0
          DO 220 J=1,NNUMT(I)
            SUM=SUM+DDEG(I,S,J)
 220      CONTINUE
          ALPHA(N+S,I)=SUM
 300  CONTINUE
      DO 400 S=1,N
        DO 330 K=1,N
          SUM=0
          DO 320 I=1,N
            DO 310 J=1,NNUMT(I)
              SUM=SUM+DDEG(I,S,J)*DDEG(I,K,J)
 310        CONTINUE
 320      CONTINUE
          ALPHA(N+S,N+K)=SUM
 330    CONTINUE
 400  CONTINUE
      DO 500 S=1,N
        DO 500 K=1,N
          SUM=0
          DO 420 J=1,NNUMT(S)
            SUM=SUM+DDEG(S,K,J)
 420      CONTINUE
          ALPHA(S,N+K)=SUM
 500  CONTINUE
C
C COMPUTE QR FACTORIZATION OF MATRIX ALPHA
      CALL QRFAQF(ALPHA,RWORK,2*N,IFLAG)
C
C REPAIR ILL CONDITIONED SCALING MATRIX
      IRMAX=IDAMAX(LENR,RWORK,1)
      RTOL=RWORK(IRMAX)*NTUR
      INDEX=1
      DO 510 I=N,2,-1
        IF (ABS(RWORK(INDEX)) .LT. RTOL) THEN
          RWORK(INDEX)=LMFPN
          IERR=1
        ENDIF
      INDEX=INDEX+I
 510  CONTINUE
      IF (ABS(RWORK(INDEX)) .LT. RTOL) THEN
        RWORK(INDEX)=LMFPN
        IERR=1
      ENDIF
C
      ENDIF
C
C CONTROL PASSES HERE IF MODE=1
C
C
C GENERATE THE COLUMN BETA
      DO 600 S=1,N
        SUM=0
        DO 550 J=1,NNUMT(S)
          SUM=SUM+COESCL(S,J)
 550    CONTINUE
        BETA(S)=-SUM
 600  CONTINUE
      DO 700 S=1,N
        SUM=0
        DO 620 I=1,N
          DO 610 J=1,NNUMT(I)
            SUM=SUM+COESCL(I,J)*DDEG(I,S,J)
 610      CONTINUE
 620    CONTINUE
        BETA(N+S)=-SUM
 700  CONTINUE
C
C SOLVE THE LINEAR SYSTEM ALPHA * X = BETA
      CALL QRSLQF(ALPHA,RWORK,BETA,XWORK,2*N)
C
C GENERATE FACE, FACV, AND THE MATRIX COESCL
      DO 800 I=1,N
        FACE(I)=BETA(I)
        FACV(I)=BETA(N+I)
 800  CONTINUE
      DO 900 I=1,N
        DO 820 J=1,NUMT(I)
          DUM = ABS(COEF(I,J))
          IF (DUM .EQ. 0.0) THEN
            COESCL(I,J) = 0.0
          ELSE
            SUM = FACE(I) + LOG10( DUM )
            DO 810 K=1,N
              SUM = SUM + FACV(K)*DEG(I,K,J)
 810        CONTINUE
            COESCL(I,J) = SIGN(10.0**(SUM), COEF(I,J))
          ENDIF
 820    CONTINUE
 900  CONTINUE
      RETURN
      END

.