share/hompack/fortran/sintrp.f

   1       SUBROUTINE SINTRP(X,Y,XOUT,YOUT,YPOUT,NEQN,KOLD,PHI,IVC,IV,KGI,GI,
   2      1                                                ALPHA,G,W,XOLD,P)
   3 C
   4 C***BEGIN PROLOGUE  SINTRP
   5 C***DATE WRITTEN   740101   (YYMMDD)
   6 C***REVISION DATE  840201   (YYMMDD)
   7 C***CATEGORY NO.  D2A2
   8 C***KEYWORDS  INITIAL VALUE ORDINARY DIFFERENTIAL EQUATIONS,
   9 C             VARIABLE ORDER ADAMS METHODS, SMOOTH INTERPOLANT FOR
  10 C             DEABM IN THE DEPAC PACKAGE
  11 C***AUTHOR  SHAMPINE, L.F.,  SNLA
  12 C           GORDON, M.K.
  13 C             MODIFIED BY H.A. WATTS
  14 C***PURPOSE  APPROXIMATES THE SOLUTION AT XOUT BY EVALUATING THE
  15 C            POLYNOMIAL COMPUTED IN STEPS AT XOUT.  MUST BE USED IN
  16 C            CONJUNCTION WITH STEPS.
  17 C***DESCRIPTION
  18 C
  19 C   WRITTEN BY L. F. SHAMPINE AND M. K. GORDON
  20 C
  21 C   ABSTRACT
  22 C
  23 C
  24 C   THE METHODS IN SUBROUTINE  STEPS  APPROXIMATE THE SOLUTION NEAR  X
  25 C   BY A POLYNOMIAL.  SUBROUTINE  SINTRP  APPROXIMATES THE SOLUTION AT
  26 C   XOUT  BY EVALUATING THE POLYNOMIAL THERE.  INFORMATION DEFINING THIS
  27 C   POLYNOMIAL IS PASSED FROM  STEPS  SO  SINTRP  CANNOT BE USED ALONE.
  28 C
  29 C   THIS CODE IS COMPLETELY EXPLAINED AND DOCUMENTED IN THE TEXT,
  30 C   COMPUTER SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS, THE INITIAL
  31 C   VALUE PROBLEM  BY L. F. SHAMPINE AND M. K. GORDON.
  32 C   FURTHER DETAILS ON USE OF THIS CODE ARE AVAILABLE IN *SOLVING
  33 C   ORDINARY DIFFERENTIAL EQUATIONS WITH ODE, STEP, AND INTRP*,
  34 C   BY L. F. SHAMPINE AND M. K. GORDON, SLA-73-1060.
  35 C
  36 C   INPUT TO SINTRP --
  37 C
  38 C   THE USER PROVIDES STORAGE IN THE CALLING PROGRAM FOR THE ARRAYS IN
  39 C   THE CALL LIST
  40 C      DIMENSION Y(NEQN),YOUT(NEQN),YPOUT(NEQN),PHI(NEQN,16),P(NEQN),
  41 C                ALPHA(12),G(13),W(12),GI(11),IV(10)
  42 C   AND DEFINES
  43 C      XOUT -- POINT AT WHICH SOLUTION IS DESIRED.
  44 C   THE REMAINING PARAMETERS ARE DEFINED IN  STEPS  AND PASSED TO
  45 C   SINTRP  FROM THAT SUBROUTINE.
  46 C
  47 C   OUTPUT FROM  SINTRP --
  48 C
  49 C      YOUT(*) -- SOLUTION AT  XOUT
  50 C      YPOUT(*) -- DERIVATIVE OF SOLUTION AT  XOUT
  51 C
  52 C   THE REMAINING PARAMETERS ARE RETURNED UNALTERED FROM THEIR INPUT
  53 C   VALUES.  INTEGRATION WITH  STEPS  MAY BE CONTINUED.
  54 C
  55 C***REFERENCES  SHAMPINE L.F., GORDON M.K., *SOLVING ORDINARY
  56 C                 DIFFERENTIAL EQUATIONS WITH ODE, STEP, AND INTRP*,
  57 C                 SLA-73-1060, SANDIA LABORATORIES, 1973.
  58 C               WATTS H.A., SHAMPINE L.F., *A SMOOTHER INTERPOLANT FOR
  59 C                 DE/STEP,INTRP : II*, SAND84-0293, SANDIA LABORATORIES,
  60 C                 1984.
  61 C***ROUTINES CALLED  (NONE)
  62 C***END PROLOGUE  SINTRP
  63 C
  64       DOUBLE PRECISION ALP,ALPHA,C,G,GAMMA,GDI,GDIF,GI,GTEMP,
  65      1  H,HI,HMU,P,PHI,RMU,SIGMA,TEMP1,TEMP2,TEMP3,W,WTEMP,
  66      2  X,XI,XIM1,XIQ,XOLD,XOUT,Y,YOUT,YPOUT
  67       INTEGER I,IQ,IV,IVC,IW,J,JQ,KGI,KOLD,KP1,KP2,L,M,NEQN
  68 C
  69       DIMENSION Y(NEQN),YOUT(NEQN),YPOUT(NEQN),PHI(NEQN,16),P(NEQN)
  70       DIMENSION GTEMP(13),C(13),WTEMP(13),G(13),W(12),ALPHA(12),
  71      1          GI(11),IV(10)
  72 C
  73 C***FIRST EXECUTABLE STATEMENT
  74       KP1 = KOLD + 1
  75       KP2 = KOLD + 2
  76 C
  77       HI = XOUT - XOLD
  78       H = X - XOLD
  79       XI = HI/H
  80       XIM1 = XI - 1.
  81 C
  82 C   INITIALIZE WTEMP(*) FOR COMPUTING GTEMP(*)
  83 C
  84       XIQ = XI
  85       DO 10 IQ = 1,KP1
  86         XIQ = XI*XIQ
  87         TEMP1 = IQ*(IQ+1)
  88  10     WTEMP(IQ) = XIQ/TEMP1
  89 C
  90 C   COMPUTE THE DOUBLE INTEGRAL TERM GDI
  91 C
  92       IF (KOLD .LE. KGI) GO TO 50
  93       IF (IVC .GT. 0) GO TO 20
  94       GDI = 1.0/TEMP1
  95       M = 2
  96       GO TO 30
  97  20   IW = IV(IVC)
  98       GDI = W(IW)
  99       M = KOLD - IW + 3
 100  30   IF (M .GT. KOLD) GO TO 60
 101       DO 40 I = M,KOLD
 102  40     GDI = W(KP2-I) - ALPHA(I)*GDI
 103       GO TO 60
 104  50   GDI = GI(KOLD)
 105 C
 106 C   COMPUTE GTEMP(*) AND C(*)
 107 C
 108  60   GTEMP(1) = XI
 109       GTEMP(2) = 0.5*XI*XI
 110       C(1) = 1.0
 111       C(2) = XI
 112       IF (KOLD .LT. 2) GO TO 90
 113       DO 80 I = 2,KOLD
 114         ALP = ALPHA(I)
 115         GAMMA = 1.0 + XIM1*ALP
 116         L = KP2 - I
 117         DO 70 JQ = 1,L
 118  70       WTEMP(JQ) = GAMMA*WTEMP(JQ) - ALP*WTEMP(JQ+1)
 119         GTEMP(I+1) = WTEMP(1)
 120  80     C(I+1) = GAMMA*C(I)
 121 C
 122 C   DEFINE INTERPOLATION PARAMETERS
 123 C
 124  90   SIGMA = (WTEMP(2) - XIM1*WTEMP(1))/GDI
 125       RMU = XIM1*C(KP1)/GDI
 126       HMU = RMU/H
 127 C
 128 C   INTERPOLATE FOR THE SOLUTION -- YOUT
 129 C   AND FOR THE DERIVATIVE OF THE SOLUTION -- YPOUT
 130 C
 131       DO 100 L = 1,NEQN
 132         YOUT(L) = 0.0
 133  100    YPOUT(L) = 0.0
 134       DO 120 J = 1,KOLD
 135         I = KP2 - J
 136         GDIF = G(I) - G(I-1)
 137         TEMP2 = (GTEMP(I) - GTEMP(I-1)) - SIGMA*GDIF
 138         TEMP3 = (C(I) - C(I-1)) + RMU*GDIF
 139         DO 110 L = 1,NEQN
 140           YOUT(L) = YOUT(L) + TEMP2*PHI(L,I)
 141  110      YPOUT(L) = YPOUT(L) + TEMP3*PHI(L,I)
 142  120    CONTINUE
 143       DO 130 L = 1,NEQN
 144         YOUT(L) = ((1.0 - SIGMA)*P(L) + SIGMA*Y(L)) +
 145      1             H*(YOUT(L) + (GTEMP(1) - SIGMA*G(1))*PHI(L,1))
 146  130    YPOUT(L) = HMU*(P(L) - Y(L)) +
 147      1                (YPOUT(L) + (C(1) + RMU*G(1))*PHI(L,1))
 148 C
 149       RETURN
 150       END