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