File: banded5x5.f

package info (click to toggle)
python-scipy 0.18.1-2
  • links: PTS, VCS
  • area: main
  • in suites: stretch
  • size: 75,464 kB
  • ctags: 79,406
  • sloc: python: 143,495; cpp: 89,357; fortran: 81,650; ansic: 79,778; makefile: 364; sh: 265
file content (240 lines) | stat: -rw-r--r-- 6,668 bytes parent folder | download | duplicates (6)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
c   banded5x5.f
c
c   This Fortran library contains implementations of the
c   differential equation
c       dy/dt = A*y
c   where A is a 5x5 banded matrix (see below for the actual
c   values).  These functions will be used to test
c   scipy.integrate.odeint.
c
c   The idea is to solve the system two ways: pure Fortran, and
c   using odeint.  The "pure Fortran" solver is implemented in
c   the subroutine banded5x5_solve below.  It calls LSODA to
c   solve the system.
c
c   To solve the same system using odeint, the functions in this
c   file are given a python wrapper using f2py.  Then the code
c   in test_odeint_jac.py uses the wrapper to implement the
c   equation and Jacobian functions required by odeint.  Because
c   those functions ultimately call the Fortran routines defined
c   in this file, the two method (pure Fortran and odeint) should
c   produce exactly the same results.  (That's assuming floating
c   point calculations are deterministic, which can be an
c   incorrect assumption.)  If we simply re-implemented the
c   equation and Jacobian functions using just python and numpy,
c   the floating point calculations would not be performed in
c   the same sequence as in the Fortran code, and we would obtain
c   different answers.  The answer for either method would be
c   numerically "correct", but the errors would be different,
c   and the counts of function and Jacobian evaluations would
c   likely be different.
c
      block data jacobian
      implicit none

      double precision bands
      dimension bands(4,5)
      common /jac/ bands

c     The data for a banded Jacobian stored in packed banded
c     format.  The full Jacobian is
c
c           -1,  0.25,     0,     0,     0
c         0.25,    -5,  0.25,     0,     0
c         0.10,  0.25,   -25,  0.25,     0
c            0,  0.10,  0.25,  -125,  0.25
c            0,     0,  0.10,  0.25,  -625
c
c     The columns in the following layout of numbers are
c     the upper diagonal, main diagonal and two lower diagonals
c     (i.e. each row in the layout is a column of the packed
c     banded Jacobian).  The values 0.00D0 are in the "don't
c     care" positions.

      data bands/
     +      0.00D0,   -1.0D0, 0.25D0, 0.10D0,
     +      0.25D0,   -5.0D0, 0.25D0, 0.10D0,
     +      0.25D0,  -25.0D0, 0.25D0, 0.10D0,
     +      0.25D0, -125.0D0, 0.25D0, 0.00D0,
     +      0.25D0, -625.0D0, 0.00D0, 0.00D0
     +      /

      end

      subroutine getbands(jac)
      double precision jac
      dimension jac(4, 5)
cf2py intent(out) jac

      double precision bands
      dimension bands(4,5)
      common /jac/ bands

      integer i, j
      do 5 i = 1, 4
          do 5 j = 1, 5
              jac(i, j) = bands(i, j)
 5    continue

      return
      end

c
c Differential equations, right-hand-side
c
      subroutine banded5x5(n, t, y, f)
      implicit none
      integer n
      double precision t, y, f
      dimension y(n), f(n)

      double precision bands
      dimension bands(4,5)
      common /jac/ bands

      f(1) = bands(2,1)*y(1) + bands(1,2)*y(2)
      f(2) = bands(3,1)*y(1) + bands(2,2)*y(2) + bands(1,3)*y(3)
      f(3) = bands(4,1)*y(1) + bands(3,2)*y(2) + bands(2,3)*y(3)
     +                                         + bands(1,4)*y(4)
      f(4) = bands(4,2)*y(2) + bands(3,3)*y(3) + bands(2,4)*y(4)
     +                                         + bands(1,5)*y(5)
      f(5) = bands(4,3)*y(3) + bands(3,4)*y(4) + bands(2,5)*y(5)

      return
      end

c
c Jacobian
c
c The subroutine assumes that the full Jacobian is to be computed.
c ml and mu are ignored, and nrowpd is assumed to be n.
c
      subroutine banded5x5_jac(n, t, y, ml, mu, jac, nrowpd)
      implicit none
      integer n, ml, mu, nrowpd
      double precision t, y, jac
      dimension y(n), jac(nrowpd, n)

      integer i, j

      double precision bands
      dimension bands(4,5)
      common /jac/ bands

      do 15 i = 1, 4
          do 15 j = 1, 5
              if ((i - j) .gt. 0) then
                  jac(i - j, j) = bands(i, j)
              end if
15    continue

      return
      end

c
c Banded Jacobian
c
c ml = 2, mu = 1
c
      subroutine banded5x5_bjac(n, t, y, ml, mu, bjac, nrowpd)
      implicit none
      integer n, ml, mu, nrowpd
      double precision t, y, bjac
      dimension y(5), bjac(nrowpd, n)

      integer i, j

      double precision bands
      dimension bands(4,5)
      common /jac/ bands

      do 20 i = 1, 4
          do 20 j = 1, 5
              bjac(i, j) = bands(i, j)
 20   continue

      return
      end


      subroutine banded5x5_solve(y, nsteps, dt, jt, nst, nfe, nje)

c     jt is the Jacobian type:
c         jt = 1  Use the full Jacobian.
c         jt = 4  Use the banded Jacobian.
c     nst, nfe and nje are outputs:
c         nst:  Total number of internal steps
c         nfe:  Total number of function (i.e. right-hand-side)
c               evaluations
c         nje:  Total number of Jacobian evaluations

      implicit none

      external banded5x5
      external banded5x5_jac
      external banded5x5_bjac
      external LSODA

c     Arguments...
      double precision y, dt
      integer nsteps, jt, nst, nfe, nje
cf2py intent(inout) y
cf2py intent(in) nsteps, dt, jt
cf2py intent(out) nst, nfe, nje

c     Local variables...
      double precision atol, rtol, t, tout, rwork
      integer iwork
      dimension y(5), rwork(500), iwork(500)
      integer neq, i
      integer itol, iopt, itask, istate, lrw, liw

c     Common block...
      double precision jacband
      dimension jacband(4,5)
      common /jac/ jacband

c     --- t range ---
      t = 0.0D0

c     --- Solver tolerances ---
      rtol = 1.0D-11
      atol = 1.0D-13
      itol = 1

c     --- Other LSODA parameters ---
      neq = 5
      itask = 1
      istate = 1
      iopt = 0
      iwork(1) = 2
      iwork(2) = 1
      lrw = 500
      liw = 500

c     --- Call LSODA in a loop to compute the solution ---
      do 40 i = 1, nsteps
          tout = i*dt
          if (jt .eq. 1) then
              call LSODA(banded5x5, neq, y, t, tout,
     &               itol, rtol, atol, itask, istate, iopt,
     &               rwork, lrw, iwork, liw,
     &               banded5x5_jac, jt)
          else
              call LSODA(banded5x5, neq, y, t, tout,
     &               itol, rtol, atol, itask, istate, iopt,
     &               rwork, lrw, iwork, liw,
     &               banded5x5_bjac, jt)
          end if
 40       if (istate .lt. 0) goto 80

      nst = iwork(11)
      nfe = iwork(12)
      nje = iwork(13)

      return

 80   write (6,89) istate
 89   format(1X,"Error: istate=",I3)
      return
      end