You cannot select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

241 lines
6.5 KiB
FortranFixed

6 years ago
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