share/hompack/fortran/ffunp.f

   1       SUBROUTINE FFUNP(N,NUMT,MMAXT,KDEG,COEF,CL,X,
   2      $ XX,TRM,DTRM,CLX,DXNP1,F,DF)
   3 C
   4 C FFUNP  EVALUATES THE SYSTEM "F(X)=0" AND ITS PARTIAL
   5 C DERIVATIVES, USING THE "TABLEAU" INPUT: N,NUMT,KDEG,COEF.
   6 C
   7 C FFUNP  CAN BE MADE MORE EFFICIENT BY CUSTOMIZING IT TO
   8 C PARTICULAR SYSTEM TYPES.  FOR EXAMPLE,
   9 C IF X(1)**2 AND X(1)**3 ARE USED IN SEVERAL
  10 C EQUATIONS, THE CURRENT  FFUNP  RECOMPUTES BOTH OF THESE FOR
  11 C EACH EQUATION.  BUT (OF COURSE) WE CAN COMPUTE
  12 C X1SQ=X(1)**2 AND X1CU=XSQ(1)*X(1), AND
  13 C USE THESE IN EACH OF THE EQUATIONS.
  14 C
  15 C THE PART OF THE CODE BELOW LABELED "BLOCK A" CAN BE
  16 C CUSTOMIZED IN THIS WAY.   (THE CODE OUTSIDE OF
  17 C BLOCK A CONCERNS THE PROJECTIVE TRANSFORMATION AND NEED NOT
  18 C BE CHANGED.)  HOWEVER, BLOCK A REQUIRES THE HOMOGENEOUS FORM
  19 C OF THE POLYNOMIALS RATHER THAN THE STANDARD FORM.  FURTHER,
  20 C THE PARTIAL DERIVATIVES WITH RESPECT TO ALL N+1 PROJECTIVE
  21 C VARIABLES MUST BE COMPUTED.  MORE EXPLICITLY,
  22 C THE ORIGINAL SYSTEM, F(X)=0, IS GIVEN IN "NON-HOMOGENEOUS FORM" AS
  23 C DESCRIBED IN SUBROUTINE POLSYS.  F(X)  IS
  24 C REPRESENTED IN "HOMOGENEOUS FORM" AS FOLLOWS:
  25 C
  26 C              NUMT(J)
  27 C
  28 C    F(J) =     SUM   TRM(J,K)
  29 C
  30 C               K=1
  31 C
  32 C WHERE  TRM(J,K)=COEF(J,K) * XX(J,1,K)*XX(J,2,K)* ... *XX(J,N+1,K)
  33 C
  34 C WITH XX(J,L,K) = X(L)**KDEG(J,L,K) FOR J=1 TO N, L=1 TO N, AND
  35 C K=1 TO NUMT(J), AND WITH XX(J,N+1,K) = XNP1**KDEG(J,N+1,K) FOR J=1 TO
  36 C N AND K=1 TO NUMT(J), WHERE  XNP1  IS THE "HOMOGENEOUS COORDINATE,"
  37 C KDEG(J,N+1,K)=IDEG(J)-(KDEG(J,1,K)+ ... + KDEG(J,N,K)),
  38 C AND IDEG(J) THE DEGREE OF THE J-TH EQUATION.   XNP1  IS GENERATED
  39 C FROM  X  AND  CL  BEFORE BLOCK A.
  40 C
  41 C IN THIS DISCUSSION WE HAVE OMITTED, FOR SIMPLICITY OF
  42 C EXPOSITION, THE LEADING INDEX, WHICH DIFFERENTIATES THE
  43 C REAL AND IMAGINARY PARTS.  HOWEVER, THIS INDEX MUST NOT BE
  44 C OMITTED IN THE CODE.
  45 C
  46 C WE COMPLETE THE EXPOSITION OF "REPLACING BLOCK A WITH MORE EFFICIENT
  47 C CODE" WITH AN EXPLICIT EXAMPLE.  FIRST, THE SYSTEM IS DESCRIBED.
  48 C THEN THE CODE THAT SHOULD BE USED IS GIVEN (COMMENTED OUT).
  49 C IN TESTS  POLSYS  WITH THE MORE EFFICIENT  FFUNP  RAN ABOUT TWICE AS
  50 C FAST AS WITH THE GENERIC  FFUNP .
  51 C
  52 C HERE IS THE SYSTEM TO BE SOLVED:
  53 C
  54 C     F(1) = COEF(1,1) * X(1)**4
  55 C    &     + COEF(1,2) * X(1)**3 * X(2)
  56 C    &     + COEF(1,3) * X(1)**3
  57 C    &     + COEF(1,4) * X(1)
  58 C    &     + COEF(1,5)
  59 C     F(2) = COEF(2,1) * X(1)     * X(2)**2
  60 C    &     + COEF(2,2)              X(2)**2
  61 C    &     + COEF(2,3)
  62 C
  63 C THE REPLACEMENT CODE REQUIRES THE FOLLOWING DECLARATIONS:
  64 C     DOUBLE PRECISION X1SQ,X1CU,X2SQ,X3SQ,X3CU,
  65 C    &  TEMPA,TEMPB,TEMPC,TEMPD,TEMPE,TEMPF
  66 C     DIMENSION X1SQ(2),X1CU(2),X2SQ(2),X3SQ(2),X3CU(2),
  67 C    &  TEMPA(2),TEMPB(2),TEMPC(2),TEMPD(2),TEMPE(2),TEMPF(2)
  68 C
  69 C HERE IS CODE TO REPLACE BLOCK A:
  70 C
  71 C******************  BEGIN BLOCK A  *******************
  72 C
  73 C     CALL MULP(X(1,1),X(1,1),X1SQ)
  74 C     CALL MULP(X1SQ  ,X(1,1),X1CU)
  75 C     CALL MULP(X(1,2),X(1,2),X2SQ)
  76 C     CALL MULP(XNP1,  XNP1,  X3SQ)
  77 C     CALL MULP(X3SQ  ,XNP1,  X3CU)
  78 C
  79 C     DO 1 I=1,2
  80 C       TEMPA(I)=   COEF(1,1) * X(I,1)
  81 C    &            + COEF(1,2) * X(I,2)
  82 C    &            + COEF(1,3) * XNP1(I)
  83 C       TEMPB(I)=   COEF(1,4) * X(I,1)
  84 C    &            + COEF(1,5) * XNP1(I)
  85 C 1   CONTINUE
  86 C
  87 C     CALL MULP(X1SQ,  TEMPA,TEMPC)
  88 C     CALL MULP(X(1,1),TEMPC,TEMPD)
  89 C     CALL MULP(X3SQ,  TEMPB,TEMPE)
  90 C     CALL MULP(XNP1,  TEMPE,TEMPF)
  91 C
  92 C     DO 2 I=1,2
  93 C       F(I,1)=TEMPD(I) + TEMPF(I)
  94 C       DF(I,1,1)= 3. *TEMPC(I) + COEF(1,1)*X1CU(I) + COEF(1,4)*X3CU(I)
  95 C       DF(I,1,2)= COEF(1,2) * X1CU(I)
  96 C       DF(I,1,3)= COEF(1,3)*X1CU(I) + 3. *TEMPE(I) + COEF(1,5)*X3CU(I)
  97 C
  98 C       TEMPA(I) = COEF(2,1) * X(I,1) + COEF(2,2) * XNP1(I)
  99 C  2  CONTINUE
 100 C
 101 C     CALL MULP(TEMPA,X(1,2),TEMPB)
 102 C     CALL MULP(TEMPB,X(1,2),TEMPC)
 103 C
 104 C     DO 3 I=1,2
 105 C       F(I,2) = TEMPC(I) + COEF(2,3) * X3CU(I)
 106 C       DF(I,2,1) = COEF(2,1) * X2SQ(I)
 107 C       DF(I,2,2) = 2. * TEMPB(I)
 108 C       DF(I,2,3) = COEF(2,2) * X2SQ(I) + COEF(2,3) * 3. * X3SQ(I)
 109 C  3  CONTINUE
 110 C******************  END OF BLOCK A  *******************
 111 C
 112 C ON INPUT:
 113 C
 114 C N  IS THE NUMBER OF EQUATIONS AND VARIABLES.
 115 C
 116 C NUMT(J)  IS THE NUMBER OF TERMS IN THE JTH EQUATION.
 117 C
 118 C MMAXT  IS AN UPPER BOUND ON NUMT(J) FOR J=1 TO N.
 119 C
 120 C KDEG(J,L,K)  IS THE DEGREE OF THE L-TH VARIABLE IN THE K-TH TERM
 121 C   OF THE J-TH EQUATION.
 122 C
 123 C COEF(J,K)  IS THE K-TH COEFFICIENT OF THE J-TH EQUATION.
 124 C
 125 C CL  IS USED TO DEFINE THE PROJECTIVE TRANSFORMATION.  IF
 126 C   THE PROJECTIVE TRANSFORMATION IS NOT SPECIFIED, THEN  CL
 127 C   CONTAINS DUMMY VALUES.
 128 C
 129 C X(1,J), X(2,J)  ARE THE REAL AND IMAGINARY PARTS RESPECTIVELY OF
 130 C   THE J-TH INDEPENDENT VARIABLE.
 131 C
 132 C XX, TRM, DTRM, CLX, DXNP1  ARE WORKSPACE VARIABLES.
 133 C
 134 C ON OUTPUT:
 135 C
 136 C F(1,J), F(2,J)  ARE THE REAL AND IMAGINARY PARTS RESPECTIVELY OF
 137 C   THE J-TH EQUATION.
 138 C
 139 C DF(1,J,K), DF(2,J,K)  ARE THE REAL AND IMAGINARY PARTS RESPECTIVELY
 140 C   OF THE K-TH PARTIAL DERIVATIVE OF THE J-TH EQUATION.
 141 C
 142 C
 143 C VARIABLES: XNP1,TEMP1,TEMP2.
 144 C
 145 C NOTE:  XNP1(1), XNP1(2)  ARE THE REAL AND IMAGINARY PARTS,
 146 C   RESPECTIVELY, OF THE PROJECTIVE VARIABLE.  XNP1  IS UNITY
 147 C   IF THE PROJECTIVE TRANSFORMATION IS NOT SPECIFIED.
 148 C
 149 C  SUBROUTINES: MULP,POWP,DIVP.
 150 C
 151 C
 152 C DECLARATION OF INPUT AND OUTPUT:
 153       INTEGER N,NUMT,MMAXT,KDEG
 154       DOUBLE PRECISION COEF,CL,X,XX,TRM,DTRM,CLX,DXNP1,F,DF
 155       DIMENSION NUMT(N),KDEG(N,N+1,MMAXT),
 156      $  COEF(N,MMAXT),CL(2,N+1),X(2,N),
 157      $  XX(2,N,N+1,MMAXT),TRM(2,N,MMAXT),DTRM(2,N,N+1,MMAXT),
 158      $  CLX(2,N),DXNP1(2,N),F(2,N),DF(2,N,N+1)
 159 C
 160 C DECLARATION OF VARIABLES:
 161       INTEGER I,IERR,J,K,L,M,NNNN,NP1
 162       DOUBLE PRECISION TEMP1,TEMP2,XNP1
 163       DIMENSION TEMP1(2),TEMP2(2),XNP1(2)
 164 C
 165       NP1=N+1
 166 C
 167 C GENERATE XNP1, THE PROJECTIVE COORDINATE, AND ITS DERIVATIVES.
 168       DO 40 J=1,N
 169          CALL MULP(CL(1,J),X(1,J),CLX(1,J))
 170  40   CONTINUE
 171 C
 172       DO 60 I=1,2
 173           XNP1(I)=CL(I,NP1)
 174           DO 50 J=1,N
 175               XNP1(I) = XNP1(I) + CLX(I,J)
 176               DXNP1(I,J)=CL(I,J)
 177  50       CONTINUE
 178  60   CONTINUE
 179 C
 180 C******************  BEGIN BLOCK A  *******************
 181 C
 182 C "BLOCK A" TAKES  X  AND  XNP1  AS INPUT AND RETURNS  F
 183 C AND  DF  AS OUTPUT.   F  IS THE HOMOGENEOUS FORM OF THE
 184 C ORIGINAL  F, AND  DF  CONSISTS OF THE PARTIAL
 185 C DERIVATIVES OF THE HOMOGENEOUS FORM OF  F  WITH RESPECT
 186 C TO THE N+1 VARIABLES X(1), ... ,X(N), XNP1.
 187 C
 188 C BEGIN "COMPUTE F"
 189 C
 190       DO 100 J=1,N
 191         DO 100 K=1,NUMT(J)
 192           CALL POWP(KDEG(J,NP1,K),XNP1, XX(1,J,NP1,K))
 193           DO 100 L=1,N
 194             CALL POWP(KDEG(J, L,K),X(1,L),XX(1,J,  L,K))
 195  100  CONTINUE
 196       DO 200 J=1,N
 197         DO 200 K=1,NUMT(J)
 198           TRM(1,J,K)=COEF(J,K)
 199           TRM(2,J,K)=0.0
 200           DO 120 L=1,NP1
 201             CALL MULP(XX(1,J,L,K), TRM(1,J,K),TEMP1)
 202             TRM(1,J,K  )=TEMP1(1)
 203             TRM(2,J,K  )=TEMP1(2)
 204  120      CONTINUE
 205  200  CONTINUE
 206       DO 300 J=1,N
 207           F(1,J)=0.0
 208           F(2,J)=0.0
 209           DO 220 I=1,2
 210           DO 220 K=1,NUMT(J)
 211               F(I,J)= F(I,J) + TRM(I,J,K)
 212  220      CONTINUE
 213  300  CONTINUE
 214 C
 215 C END OF "COMPUTE F"
 216 C
 217 C BEGIN "COMPUTE DF"
 218 C
 219       DO 400 J=1,N
 220         DO 400 K=1,NUMT(J)
 221         DO 400 M=1,NP1
 222 C
 223 C IF TERM DOES NOT INCLUDE X(M), SET PARTIAL DERIVATIVE OF TERM
 224 C   EQUAL TO ZERO.
 225           IF(KDEG(J,M,K) .EQ. 0) THEN
 226             DTRM(1,J,M,K)=0.0
 227             DTRM(2,J,M,K)=0.0
 228           ELSE
 229 C
 230 C IF TERM DOES INCLUDE X(M), TRY COMPUTING THE PARTIAL BY DIVIDING
 231 C   THE TERM BY X(M).
 232             IF(M.LE.N) CALL DIVP(TRM(1,J,K),X(1,M),DTRM(1,J,M,K),IERR)
 233             IF(M.EQ.NP1) CALL DIVP(TRM(1,J,K),XNP1,DTRM(1,J,M,K),IERR)
 234             IF (IERR .EQ. 0) THEN
 235               DTRM(1,J,M,K)=KDEG(J,M,K)*DTRM(1,J,M,K)
 236               DTRM(2,J,M,K)=KDEG(J,M,K)*DTRM(2,J,M,K)
 237             ELSE
 238 C
 239 C IF DIVISION WOULD CAUSE OVERFLOW, GENERATE THE PARTIAL BY
 240 C   THE POLYNOMIAL FORMULA.
 241               DTRM(1,J,M,K)=COEF(J,K)
 242               DTRM(2,J,M,K)=0.0
 243               DO 320 L=1,NP1
 244                 IF (L .EQ. M) GOTO 320
 245                 CALL MULP(XX(1,J,L,K),DTRM(1,J,M,K),TEMP1)
 246                 DTRM(1,J,M,K)=TEMP1(1)
 247                 DTRM(2,J,M,K)=TEMP1(2)
 248  320          CONTINUE
 249               NNNN=KDEG(J,M,K)-1
 250               IF (M .LE. N) CALL POWP(NNNN,X(1,M),TEMP2)
 251               IF (M .EQ. NP1) CALL POWP(NNNN,XNP1 ,TEMP2)
 252               CALL MULP(TEMP2,TEMP1,DTRM(1,J,M,K))
 253               DTRM(1,J,M,K)=KDEG(J,M,K)*DTRM(1,J,M,K)
 254               DTRM(2,J,M,K)=KDEG(J,M,K)*DTRM(2,J,M,K)
 255             END IF
 256           END IF
 257  400  CONTINUE
 258       DO 600 J=1,N
 259       DO 600 M=1,NP1
 260           DF(1,J,M)=0.0
 261           DF(2,J,M)=0.0
 262           DO 420 I=1,2
 263             DO 420 K=1,NUMT(J)
 264               DF(I,J,M)= DF(I,J,M) + DTRM(I,J,M,K)
 265  420      CONTINUE
 266  600  CONTINUE
 267 C
 268 C END OF "COMPUTE DF"
 269 C*******************  END BLOCK A  ********************
 270 C
 271 C CONVERT  DF  TO BE PARTIALS WITH RESPECT TO  X(1), ... ,X(N),
 272 C BY APPLYING THE CHAIN RULE WITH  XNP1  CONSIDERED A FUNCTION OF
 273 C OF  X(1), ... ,X(N).
 274 C
 275       DO 700 J=1,N
 276         DO 700 K=1,N
 277           CALL MULP(DF(1,J,NP1),DXNP1(1,K),TEMP1)
 278           DO 700 I=1,2
 279             DF(I,J,K)=DF(I,J,K)+TEMP1(I)
 280  700  CONTINUE
 281       RETURN
 282       END