SUBROUTINE FFUNP(N,NUMT,MAXT,KDEG,COEF,CL,X,
& XX,TRM,DTRM,CLX,DXNP1,F,DF)
C
C FFUNP EVALUATES THE SYSTEM "F(X)=0" AND ITS PARTIAL
C DERIVATIVES, USING THE "TABLEAU" INPUT: N,NUMT,KDEG,COEF.
C
C FFUNP CAN BE MADE MORE EFFICIENT BY CUSTOMIZING IT TO
C PARTICULAR SYSTEM TYPES. FOR EXAMPLE,
C IF X(1)**2 AND X(1)**3 ARE USED IN SEVERAL
C EQUATIONS, THE CURRENT FFUNP RECOMPUTES BOTH OF THESE FOR
C EACH EQUATION. BUT (OF COURSE) WE CAN COMPUTE
C X1SQ=X(1)**2 AND X1CU=XSQ(1)*X(1), AND
C USE THESE IN EACH OF THE EQUATIONS.
C
C THE PART OF THE CODE BELOW LABELED "BLOCK A" CAN BE
C CUSTOMIZED IN THIS WAY. (THE CODE OUTSIDE OF
C BLOCK A CONCERNS THE PROJECTIVE TRANSFORMATION AND NEED NOT
C BE CHANGED.) HOWEVER, BLOCK A REQUIRES THE HOMOGENEOUS FORM
C OF THE POLYNOMIALS RATHER THAN THE STANDARD FORM. FURTHER,
C THE PARTIAL DERIVATIVES WITH RESPECT TO ALL N+1 PROJECTIVE
C VARIABLES MUST BE COMPUTED. MORE EXPLICITLY,
C THE ORIGINAL SYSTEM, F(X)=0, IS GIVEN IN "NON-HOMOGENEOUS FORM" AS
C DESCRIBED IN SUBROUTINE POLSYS. F(X) IS
C REPRESENTED IN "HOMOGENEOUS FORM" AS FOLLOWS:
C
C NUMT(J)
C
C F(J) = SUM TRM(J,K)
C
C K=1
C
C WHERE TRM(J,K)=COEF(J,K) * XX(J,1,K)*XX(J,2,K)* ... *XX(J,N+1,K)
C
C WITH XX(J,L,K) = X(L)**KDEG(J,L,K) FOR J=1 TO N, L=1 TO N, AND
C K=1 TO NUMT(J), AND WITH XX(J,N+1,K) = XNP1**KDEG(J,N+1,K) FOR J=1 TO
C N AND K=1 TO NUMT(J), WHERE XNP1 IS THE "HOMOGENEOUS COORDINATE,"
C KDEG(J,N+1,K)=IDEG(J)-(KDEG(J,1,K)+ ... + KDEG(J,N,K)),
C AND IDEG(J) THE DEGREE OF THE J-TH EQUATION. XNP1 IS GENERATED
C FROM X AND CL BEFORE BLOCK A.
C
C IN THIS DISCUSSION WE HAVE OMITTED, FOR SIMPLICITY OF
C EXPOSITION, THE LEADING INDEX, WHICH DIFFERENTIATES THE
C REAL AND IMAGINARY PARTS. HOWEVER, THIS INDEX MUST NOT BE
C OMITTED IN THE CODE.
C
C WE COMPLETE THE EXPOSITION OF "REPLACING BLOCK A WITH MORE EFFICIENT
C CODE" WITH AN EXPLICIT EXAMPLE. FIRST, THE SYSTEM IS DESCRIBED.
C THEN THE CODE THAT SHOULD BE USED IS GIVEN (COMMENTED OUT).
C IN TESTS POLSYS WITH THE MORE EFFICIENT FFUNP RAN ABOUT TWICE AS
C FAST AS WITH THE GENERIC FFUNP .
C
C HERE IS THE SYSTEM TO BE SOLVED:
C
C F(1) = COEF(1,1) * X(1)**4
C & + COEF(1,2) * X(1)**3 * X(2)
C & + COEF(1,3) * X(1)**3
C & + COEF(1,4) * X(1)
C & + COEF(1,5)
C F(2) = COEF(2,1) * X(1) * X(2)**2
C & + COEF(2,2) X(2)**2
C & + COEF(2,3)
C
C THE REPLACEMENT CODE REQUIRES THE FOLLOWING DECLARATIONS:
C DOUBLE PRECISION X1SQ,X1CU,X2SQ,X3SQ,X3CU,
C & TEMPA,TEMPB,TEMPC,TEMPD,TEMPE,TEMPF
C DIMENSION X1SQ(2),X1CU(2),X2SQ(2),X3SQ(2),X3CU(2),
C & TEMPA(2),TEMPB(2),TEMPC(2),TEMPD(2),TEMPE(2),TEMPF(2)
C
C HERE IS CODE TO REPLACE BLOCK A:
C
C****************** BEGIN BLOCK A *******************
C
C CALL MULP(X(1,1),X(1,1),X1SQ)
C CALL MULP(X1SQ ,X(1,1),X1CU)
C CALL MULP(X(1,2),X(1,2),X2SQ)
C CALL MULP(XNP1, XNP1, X3SQ)
C CALL MULP(X3SQ ,XNP1, X3CU)
C
C DO 1 I=1,2
C TEMPA(I)= COEF(1,1) * X(I,1)
C & + COEF(1,2) * X(I,2)
C & + COEF(1,3) * XNP1(I)
C TEMPB(I)= COEF(1,4) * X(I,1)
C & + COEF(1,5) * XNP1(I)
C 1 CONTINUE
C
C CALL MULP(X1SQ, TEMPA,TEMPC)
C CALL MULP(X(1,1),TEMPC,TEMPD)
C CALL MULP(X3SQ, TEMPB,TEMPE)
C CALL MULP(XNP1, TEMPE,TEMPF)
C
C DO 2 I=1,2
C F(I,1)=TEMPD(I) + TEMPF(I)
C DF(I,1,1)= 3. *TEMPC(I) + COEF(1,1)*X1CU(I) + COEF(1,4)*X3CU(I)
C DF(I,1,2)= COEF(1,2) * X1CU(I)
C DF(I,1,3)= COEF(1,3)*X1CU(I) + 3. *TEMPE(I) + COEF(1,5)*X3CU(I)
C
C TEMPA(I) = COEF(2,1) * X(I,1) + COEF(2,2) * XNP1(I)
C 2 CONTINUE
C
C CALL MULP(TEMPA,X(1,2),TEMPB)
C CALL MULP(TEMPB,X(1,2),TEMPC)
C
C DO 3 I=1,2
C F(I,2) = TEMPC(I) + COEF(2,3) * X3CU(I)
C DF(I,2,1) = COEF(2,1) * X2SQ(I)
C DF(I,2,2) = 2. * TEMPB(I)
C DF(I,2,3) = COEF(2,2) * X2SQ(I) + COEF(2,3) * 3. * X3SQ(I)
C 3 CONTINUE
C****************** END OF BLOCK A *******************
C
C ON INPUT:
C
C N IS THE NUMBER OF EQUATIONS AND VARIABLES.
C
C NUMT(J) IS THE NUMBER OF TERMS IN THE JTH EQUATION.
C
C MAXT IS AN UPPER BOUND ON NUMT(J) FOR J=1 TO N.
C
C KDEG(J,L,K) IS THE DEGREE OF THE L-TH VARIABLE IN THE K-TH TERM
C OF THE J-TH EQUATION.
C
C COEF(J,K) IS THE K-TH COEFFICIENT OF THE J-TH EQUATION.
C
C CL IS USED TO DEFINE THE PROJECTIVE TRANSFORMATION. IF
C THE PROJECTIVE TRANSFORMATION IS NOT SPECIFIED, THEN CL
C CONTAINS DUMMY VALUES.
C
C X(1,J), X(2,J) ARE THE REAL AND IMAGINARY PARTS RESPECTIVELY OF
C THE J-TH INDEPENDENT VARIABLE.
C
C XX, TRM, DTRM, CLX, DXNP1 ARE WORKSPACE VARIABLES.
C
C ON OUTPUT:
C
C F(1,J), F(2,J) ARE THE REAL AND IMAGINARY PARTS RESPECTIVELY OF
C THE J-TH EQUATION.
C
C DF(1,J,K), DF(2,J,K) ARE THE REAL AND IMAGINARY PARTS RESPECTIVELY
C OF THE K-TH PARTIAL DERIVATIVE OF THE J-TH EQUATION.
C
C
C VARIABLES: XNP1,TEMP1,TEMP2.
C
C NOTE: XNP1(1), XNP1(2) ARE THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF THE PROJECTIVE VARIABLE. XNP1 IS UNITY
C IF THE PROJECTIVE TRANSFORMATION IS NOT SPECIFIED.
C
C SUBROUTINES: MULP,POWP,DIVP.
C
USE REAL_PRECISION
C
INTERFACE
SUBROUTINE DIVP(XXXX,YYYY,ZZZZ,IERR)
USE REAL_PRECISION
REAL (KIND=R8), DIMENSION(2), INTENT(IN):: XXXX,YYYY
INTEGER, INTENT(OUT):: IERR
REAL (KIND=R8), DIMENSION(2), INTENT(OUT):: ZZZZ
END SUBROUTINE DIVP
SUBROUTINE MULP(XXXX,YYYY,ZZZZ)
USE REAL_PRECISION
REAL (KIND=R8), DIMENSION(2), INTENT(IN):: XXXX,YYYY
REAL (KIND=R8), DIMENSION(2), INTENT(OUT):: ZZZZ
END SUBROUTINE MULP
SUBROUTINE POWP(NNNN,XXXX,YYYY)
USE REAL_PRECISION
INTEGER, INTENT(IN):: NNNN
REAL (KIND=R8), DIMENSION(2), INTENT(IN):: XXXX
REAL (KIND=R8), DIMENSION(2), INTENT(IN OUT):: YYYY
END SUBROUTINE POWP
END INTERFACE
C
C DECLARATION OF INPUT AND OUTPUT:
INTEGER, INTENT(IN):: N,NUMT(N),MAXT,KDEG(N,N+1,MAXT)
REAL (KIND=R8), INTENT(IN):: COEF(N,MAXT),CL(2,N+1),X(2,N)
REAL (KIND=R8), INTENT(IN OUT)::
& XX(2,N,N+1,MAXT),TRM(2,N,MAXT),DTRM(2,N,N+1,MAXT)
REAL (KIND=R8), INTENT(OUT)::
& CLX(2,N),DXNP1(2,N),F(2,N),DF(2,N,N+1)
C
C DECLARATION OF LOCAL VARIABLES:
INTEGER:: IERR,J,K,L,M,NNNN,NP1
REAL (KIND=R8), DIMENSION(2):: TEMP1,TEMP2,XNP1
C
NP1=N+1
C
C GENERATE XNP1, THE PROJECTIVE COORDINATE, AND ITS DERIVATIVES.
DO J=1,N
CALL MULP(CL(1,J),X(1,J),CLX(1,J))
END DO
C
XNP1(1:2)=CL(1:2,NP1) + SUM(CLX(1:2,1:N),DIM=2)
DXNP1(1:2,1:N)=CL(1:2,1:N)
C
C****************** BEGIN BLOCK A *******************
C
C "BLOCK A" TAKES X AND XNP1 AS INPUT AND RETURNS F
C AND DF AS OUTPUT. F IS THE HOMOGENEOUS FORM OF THE
C ORIGINAL F, AND DF CONSISTS OF THE PARTIAL
C DERIVATIVES OF THE HOMOGENEOUS FORM OF F WITH RESPECT
C TO THE N+1 VARIABLES X(1), ... ,X(N), XNP1.
C
C BEGIN "COMPUTE F"
C
DO J=1,N
DO K=1,NUMT(J)
CALL POWP(KDEG(J,NP1,K),XNP1, XX(1,J,NP1,K))
DO L=1,N
CALL POWP(KDEG(J, L,K),X(1,L),XX(1,J, L,K))
END DO
END DO
END DO
TRM = 0.0
DO J=1,N
DO K=1,NUMT(J)
TRM(1,J,K)=COEF(J,K)
DO L=1,NP1
CALL MULP(XX(1,J,L,K), TRM(1,J,K),TEMP1)
TRM(1:2,J,K ) = TEMP1(1:2)
END DO
END DO
END DO
F(1:2,1:N) = SUM(TRM(1:2,1:N,:),DIM=3)
C
C END OF "COMPUTE F"
C
C BEGIN "COMPUTE DF"
C
J_LOOP: DO J=1,N
K_LOOP: DO K=1,NUMT(J)
M_LOOP: DO M=1,NP1
C
C IF TERM DOES NOT INCLUDE X(M), SET PARTIAL DERIVATIVE OF TERM
C EQUAL TO ZERO.
IF(KDEG(J,M,K) .EQ. 0) THEN
DTRM(1:2,J,M,K) = 0.0
ELSE
C
C IF TERM DOES INCLUDE X(M), TRY COMPUTING THE PARTIAL BY DIVIDING
C THE TERM BY X(M).
IF(M.LE.N) CALL DIVP(TRM(1,J,K),X(1,M),DTRM(1,J,M,K),IERR)
IF(M.EQ.NP1) CALL DIVP(TRM(1,J,K),XNP1,DTRM(1,J,M,K),IERR)
IF (IERR .EQ. 0) THEN
DTRM(1:2,J,M,K)=KDEG(J,M,K)*DTRM(1:2,J,M,K)
ELSE
C
C IF DIVISION WOULD CAUSE OVERFLOW, GENERATE THE PARTIAL BY
C THE POLYNOMIAL FORMULA.
DTRM(1,J,M,K)=COEF(J,K)
DTRM(2,J,M,K)=0.0
DO L=1,NP1
IF (L .EQ. M) CYCLE
CALL MULP(XX(1,J,L,K),DTRM(1,J,M,K),TEMP1)
DTRM(1:2,J,M,K)=TEMP1(1:2)
END DO
NNNN=KDEG(J,M,K)-1
IF (M .LE. N) CALL POWP(NNNN,X(1,M),TEMP2)
IF (M .EQ. NP1) CALL POWP(NNNN,XNP1 ,TEMP2)
CALL MULP(TEMP2,TEMP1,DTRM(1,J,M,K))
DTRM(1:2,J,M,K)=KDEG(J,M,K)*DTRM(1:2,J,M,K)
END IF
END IF
END DO M_LOOP
END DO K_LOOP
END DO J_LOOP
DO J=1,N
DF(1:2,J,1:NP1) = SUM(DTRM(1:2,J,1:NP1,1:NUMT(J)), DIM=3)
END DO
C
C END OF "COMPUTE DF"
C******************* END BLOCK A ********************
C
C CONVERT DF TO BE PARTIALS WITH RESPECT TO X(1), ... ,X(N),
C BY APPLYING THE CHAIN RULE WITH XNP1 CONSIDERED A FUNCTION OF
C OF X(1), ... ,X(N).
C
DO J=1,N
DO K=1,N
CALL MULP(DF(1,J,NP1),DXNP1(1,K),TEMP1)
DF(1:2,J,K)=DF(1:2,J,K)+TEMP1(1:2)
END DO
END DO
RETURN
END SUBROUTINE FFUNP