relibc_openlibm/slatec/dcfod.f

*DECK DCFOD
      SUBROUTINE DCFOD (METH, ELCO, TESCO)
C***BEGIN PROLOGUE  DCFOD
C***SUBSIDIARY
C***PURPOSE  Subsidiary to DDEBDF
C***LIBRARY   SLATEC
C***TYPE      DOUBLE PRECISION (CFOD-S, DCFOD-D)
C***AUTHOR  (UNKNOWN)
C***DESCRIPTION
C
C   DCFOD defines coefficients needed in the integrator package DDEBDF
C
C***SEE ALSO  DDEBDF
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   820301  DATE WRITTEN
C   890911  Removed unnecessary intrinsics.  (WRB)
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900328  Added TYPE section.  (WRB)
C***END PROLOGUE  DCFOD
C
C
      INTEGER I, IB, METH, NQ, NQM1, NQP1
      DOUBLE PRECISION AGAMQ, ELCO, FNQ, FNQM1, PC, PINT, RAGQ,
     1      RQ1FAC, RQFAC, TESCO, TSIGN, XPIN
      DIMENSION ELCO(13,12),TESCO(3,12)
C     ------------------------------------------------------------------
C      DCFOD  IS CALLED BY THE INTEGRATOR ROUTINE TO SET COEFFICIENTS
C      NEEDED THERE.  THE COEFFICIENTS FOR THE CURRENT METHOD, AS
C      GIVEN BY THE VALUE OF METH, ARE SET FOR ALL ORDERS AND SAVED.
C      THE MAXIMUM ORDER ASSUMED HERE IS 12 IF METH = 1 AND 5 IF METH =
C      2.  (A SMALLER VALUE OF THE MAXIMUM ORDER IS ALSO ALLOWED.)
C      DCFOD  IS CALLED ONCE AT THE BEGINNING OF THE PROBLEM,
C      AND IS NOT CALLED AGAIN UNLESS AND UNTIL METH IS CHANGED.
C
C      THE ELCO ARRAY CONTAINS THE BASIC METHOD COEFFICIENTS.
C      THE COEFFICIENTS EL(I), 1 .LE. I .LE. NQ+1, FOR THE METHOD OF
C      ORDER NQ ARE STORED IN ELCO(I,NQ).  THEY ARE GIVEN BY A
C      GENERATING POLYNOMIAL, I.E.,
C          L(X) = EL(1) + EL(2)*X + ... + EL(NQ+1)*X**NQ.
C      FOR THE IMPLICIT ADAMS METHODS, L(X) IS GIVEN BY
C          DL/DX = (X+1)*(X+2)*...*(X+NQ-1)/FACTORIAL(NQ-1),    L(-1) =
C      0.  FOR THE BDF METHODS, L(X) IS GIVEN BY
C          L(X) = (X+1)*(X+2)* ... *(X+NQ)/K,
C      WHERE         K = FACTORIAL(NQ)*(1 + 1/2 + ... + 1/NQ).
C
C      THE TESCO ARRAY CONTAINS TEST CONSTANTS USED FOR THE
C      LOCAL ERROR TEST AND THE SELECTION OF STEP SIZE AND/OR ORDER.
C      AT ORDER NQ, TESCO(K,NQ) IS USED FOR THE SELECTION OF STEP
C      SIZE AT ORDER NQ - 1 IF K = 1, AT ORDER NQ IF K = 2, AND AT ORDER
C      NQ + 1 IF K = 3.
C     ------------------------------------------------------------------
      DIMENSION PC(12)
C
C***FIRST EXECUTABLE STATEMENT  DCFOD
      GO TO (10,60), METH
C
   10 CONTINUE
         ELCO(1,1) = 1.0D0
         ELCO(2,1) = 1.0D0
         TESCO(1,1) = 0.0D0
         TESCO(2,1) = 2.0D0
         TESCO(1,2) = 1.0D0
         TESCO(3,12) = 0.0D0
         PC(1) = 1.0D0
         RQFAC = 1.0D0
         DO 50 NQ = 2, 12
C           ------------------------------------------------------------
C            THE PC ARRAY WILL CONTAIN THE COEFFICIENTS OF THE
C                POLYNOMIAL P(X) = (X+1)*(X+2)*...*(X+NQ-1).
C            INITIALLY, P(X) = 1.
C           ------------------------------------------------------------
            RQ1FAC = RQFAC
            RQFAC = RQFAC/NQ
            NQM1 = NQ - 1
            FNQM1 = NQM1
            NQP1 = NQ + 1
C           FORM COEFFICIENTS OF P(X)*(X+NQ-1).
C           ----------------------------------
            PC(NQ) = 0.0D0
            DO 20 IB = 1, NQM1
               I = NQP1 - IB
               PC(I) = PC(I-1) + FNQM1*PC(I)
   20       CONTINUE
            PC(1) = FNQM1*PC(1)
C           COMPUTE INTEGRAL, -1 TO 0, OF P(X) AND X*P(X).
C           -----------------------
            PINT = PC(1)
            XPIN = PC(1)/2.0D0
            TSIGN = 1.0D0
            DO 30 I = 2, NQ
               TSIGN = -TSIGN
               PINT = PINT + TSIGN*PC(I)/I
               XPIN = XPIN + TSIGN*PC(I)/(I+1)
   30       CONTINUE
C           STORE COEFFICIENTS IN ELCO AND TESCO.
C           --------------------------------
            ELCO(1,NQ) = PINT*RQ1FAC
            ELCO(2,NQ) = 1.0D0
            DO 40 I = 2, NQ
               ELCO(I+1,NQ) = RQ1FAC*PC(I)/I
   40       CONTINUE
            AGAMQ = RQFAC*XPIN
            RAGQ = 1.0D0/AGAMQ
            TESCO(2,NQ) = RAGQ
            IF (NQ .LT. 12) TESCO(1,NQP1) = RAGQ*RQFAC/NQP1
            TESCO(3,NQM1) = RAGQ
   50    CONTINUE
      GO TO 100
C
   60 CONTINUE
         PC(1) = 1.0D0
         RQ1FAC = 1.0D0
         DO 90 NQ = 1, 5
C           ------------------------------------------------------------
C            THE PC ARRAY WILL CONTAIN THE COEFFICIENTS OF THE
C                POLYNOMIAL P(X) = (X+1)*(X+2)*...*(X+NQ).
C            INITIALLY, P(X) = 1.
C           ------------------------------------------------------------
            FNQ = NQ
            NQP1 = NQ + 1
C           FORM COEFFICIENTS OF P(X)*(X+NQ).
C           ------------------------------------
            PC(NQP1) = 0.0D0
            DO 70 IB = 1, NQ
               I = NQ + 2 - IB
               PC(I) = PC(I-1) + FNQ*PC(I)
   70       CONTINUE
            PC(1) = FNQ*PC(1)
C           STORE COEFFICIENTS IN ELCO AND TESCO.
C           --------------------------------
            DO 80 I = 1, NQP1
               ELCO(I,NQ) = PC(I)/PC(2)
   80       CONTINUE
            ELCO(2,NQ) = 1.0D0
            TESCO(1,NQ) = RQ1FAC
            TESCO(2,NQ) = NQP1/ELCO(1,NQ)
            TESCO(3,NQ) = (NQ+2)/ELCO(1,NQ)
            RQ1FAC = RQ1FAC/FNQ
   90    CONTINUE
  100 CONTINUE
      RETURN
C     ----------------------- END OF SUBROUTINE DCFOD
C     -----------------------
      END