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
|