← Back to branch summary

~ubuntu-branches/ubuntu/wily/julia/wily

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/dintp.f

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot
  • Date: 2013-01-16 12:29:42 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130116122942-x86e42akjq31repw
Tags: 0.0.0+20130107.gitd9656f41-1
* New upstream snashot
* No longer try to rebuild helpdb.jl.
   + debian/rules: remove helpdb.jl from build-arch rule
   + debian/control: move back python-sphinx to Build-Depends-Indep
* debian/copyright: reflect upstream changes
* Add Build-Conflicts on libatlas3-base (makes linalg tests fail)
* debian/rules: replace obsolete USE_DEBIAN makeflag by a list of
  USE_SYSTEM_* flags
* debian/rules: on non-x86 systems, use libm instead of openlibm
* dpkg-buildflags.patch: remove patch, applied upstream
* Refreshed other patches

Show diffs side-by-side

added added

removed removed

Lines of Context:
deps/openlibm/slatec/dintp.f
 
1
*DECK DINTP
 
2
      SUBROUTINE DINTP (X, Y, XOUT, YOUT, YPOUT, NEQN, KOLD, PHI, IVC,
 
3
     +   IV, KGI, GI, ALPHA, OG, OW, OX, OY)
 
4
C***BEGIN PROLOGUE  DINTP
 
5
C***PURPOSE  Approximate the solution at XOUT by evaluating the
 
6
C            polynomial computed in DSTEPS at XOUT.  Must be used in
 
7
C            conjunction with DSTEPS.
 
8
C***LIBRARY   SLATEC (DEPAC)
 
9
C***CATEGORY  I1A1B
 
10
C***TYPE      DOUBLE PRECISION (SINTRP-S, DINTP-D)
 
11
C***KEYWORDS  ADAMS METHOD, DEPAC, INITIAL VALUE PROBLEMS, ODE,
 
12
C             ORDINARY DIFFERENTIAL EQUATIONS, PREDICTOR-CORRECTOR,
 
13
C             SMOOTH INTERPOLANT
 
14
C***AUTHOR  Watts, H. A., (SNLA)
 
15
C***DESCRIPTION
 
16
C
 
17
C   The methods in subroutine  DSTEPS  approximate the solution near  X
 
18
C   by a polynomial.  Subroutine  DINTP  approximates the solution at
 
19
C   XOUT  by evaluating the polynomial there.  Information defining this
 
20
C   polynomial is passed from  DSTEPS  so  DINTP  cannot be used alone.
 
21
C
 
22
C   Subroutine DSTEPS is completely explained and documented in the text
 
23
C   "Computer Solution of Ordinary Differential Equations, the Initial
 
24
C   Value Problem"  by L. F. Shampine and M. K. Gordon.
 
25
C
 
26
C   Input to DINTP --
 
27
C
 
28
C   The user provides storage in the calling program for the arrays in
 
29
C   the call list
 
30
C      DIMENSION Y(NEQN),YOUT(NEQN),YPOUT(NEQN),PHI(NEQN,16),OY(NEQN)
 
31
C                AND ALPHA(12),OG(13),OW(12),GI(11),IV(10)
 
32
C   and defines
 
33
C      XOUT -- point at which solution is desired.
 
34
C   The remaining parameters are defined in  DSTEPS  and passed to
 
35
C   DINTP  from that subroutine
 
36
C
 
37
C   Output from  DINTP --
 
38
C
 
39
C      YOUT(*) -- solution at  XOUT
 
40
C      YPOUT(*) -- derivative of solution at  XOUT
 
41
C   The remaining parameters are returned unaltered from their input
 
42
C   values.  Integration with  DSTEPS  may be continued.
 
43
C
 
44
C***REFERENCES  H. A. Watts, A smoother interpolant for DE/STEP, INTRP
 
45
C                 II, Report SAND84-0293, Sandia Laboratories, 1984.
 
46
C***ROUTINES CALLED  (NONE)
 
47
C***REVISION HISTORY  (YYMMDD)
 
48
C   840201  DATE WRITTEN
 
49
C   890831  Modified array declarations.  (WRB)
 
50
C   890831  REVISION DATE from Version 3.2
 
51
C   891214  Prologue converted to Version 4.0 format.  (BAB)
 
52
C   920501  Reformatted the REFERENCES section.  (WRB)
 
53
C***END PROLOGUE  DINTP
 
54
C
 
55
      INTEGER I, IQ, IV, IVC, IW, J, JQ, KGI, KOLD, KP1, KP2,
 
56
     1        L, M, NEQN
 
57
      DOUBLE PRECISION ALP, ALPHA, C, G, GDI, GDIF, GI, GAMMA, H, HI,
 
58
     1       HMU, OG, OW, OX, OY, PHI, RMU, SIGMA, TEMP1, TEMP2, TEMP3,
 
59
     2       W, X, XI, XIM1, XIQ, XOUT, Y, YOUT, YPOUT
 
60
C
 
61
      DIMENSION Y(*),YOUT(*),YPOUT(*),PHI(NEQN,16),OY(*)
 
62
      DIMENSION G(13),C(13),W(13),OG(13),OW(12),ALPHA(12),GI(11),IV(10)
 
63
C
 
64
C***FIRST EXECUTABLE STATEMENT  DINTP
 
65
      KP1 = KOLD + 1
 
66
      KP2 = KOLD + 2
 
67
C
 
68
      HI = XOUT - OX
 
69
      H = X - OX
 
70
      XI = HI/H
 
71
      XIM1 = XI - 1.D0
 
72
C
 
73
C   INITIALIZE W(*) FOR COMPUTING G(*)
 
74
C
 
75
      XIQ = XI
 
76
      DO 10 IQ = 1,KP1
 
77
        XIQ = XI*XIQ
 
78
        TEMP1 = IQ*(IQ+1)
 
79
 10     W(IQ) = XIQ/TEMP1
 
80
C
 
81
C   COMPUTE THE DOUBLE INTEGRAL TERM GDI
 
82
C
 
83
      IF (KOLD .LE. KGI) GO TO 50
 
84
      IF (IVC .GT. 0) GO TO 20
 
85
      GDI = 1.0D0/TEMP1
 
86
      M = 2
 
87
      GO TO 30
 
88
 20   IW = IV(IVC)
 
89
      GDI = OW(IW)
 
90
      M = KOLD - IW + 3
 
91
 30   IF (M .GT. KOLD) GO TO 60
 
92
      DO 40 I = M,KOLD
 
93
 40     GDI = OW(KP2-I) - ALPHA(I)*GDI
 
94
      GO TO 60
 
95
 50   GDI = GI(KOLD)
 
96
C
 
97
C   COMPUTE G(*) AND C(*)
 
98
C
 
99
 60   G(1) = XI
 
100
      G(2) = 0.5D0*XI*XI
 
101
      C(1) = 1.0D0
 
102
      C(2) = XI
 
103
      IF (KOLD .LT. 2) GO TO 90
 
104
      DO 80 I = 2,KOLD
 
105
        ALP = ALPHA(I)
 
106
        GAMMA = 1.0D0 + XIM1*ALP
 
107
        L = KP2 - I
 
108
        DO 70 JQ = 1,L
 
109
 70       W(JQ) = GAMMA*W(JQ) - ALP*W(JQ+1)
 
110
        G(I+1) = W(1)
 
111
 80     C(I+1) = GAMMA*C(I)
 
112
C
 
113
C   DEFINE INTERPOLATION PARAMETERS
 
114
C
 
115
 90   SIGMA = (W(2) - XIM1*W(1))/GDI
 
116
      RMU = XIM1*C(KP1)/GDI
 
117
      HMU = RMU/H
 
118
C
 
119
C   INTERPOLATE FOR THE SOLUTION -- YOUT
 
120
C   AND FOR THE DERIVATIVE OF THE SOLUTION -- YPOUT
 
121
C
 
122
      DO 100 L = 1,NEQN
 
123
        YOUT(L) = 0.0D0
 
124
 100    YPOUT(L) = 0.0D0
 
125
      DO 120 J = 1,KOLD
 
126
        I = KP2 - J
 
127
        GDIF = OG(I) - OG(I-1)
 
128
        TEMP2 = (G(I) - G(I-1)) - SIGMA*GDIF
 
129
        TEMP3 = (C(I) - C(I-1)) + RMU*GDIF
 
130
        DO 110 L = 1,NEQN
 
131
          YOUT(L) = YOUT(L) + TEMP2*PHI(L,I)
 
132
 110      YPOUT(L) = YPOUT(L) + TEMP3*PHI(L,I)
 
133
 120    CONTINUE
 
134
      DO 130 L = 1,NEQN
 
135
        YOUT(L) = ((1.0D0 - SIGMA)*OY(L) + SIGMA*Y(L)) +
 
136
     1             H*(YOUT(L) + (G(1) - SIGMA*OG(1))*PHI(L,1))
 
137
 130    YPOUT(L) = HMU*(OY(L) - Y(L)) +
 
138
     1                (YPOUT(L) + (C(1) + RMU*OG(1))*PHI(L,1))
 
139
C
 
140
      RETURN
 
141
      END