Merged branch master into master

rveltz · Nov 27, 2016 · bcfdcdf · bcfdcdf
2 parents 82143be + 913574f
commit bcfdcdf
Show file tree

Hide file tree

Showing 5 changed files with 58 additions and 32 deletions.
diff --git a/README.md b/README.md
@@ -1 +1,8 @@
-# LSODA.jl## Introduction**LSODA.jl** is a Julia package that interfaces to the [liblsoda](https://github.com/sdwfrost/liblsoda) library, developped by [Simon Frost](http://www.vet.cam.ac.uk/directory/sdf22@cam.ac.uk) ([@sdwfrost](http://github.com/sdwfrost)), thereby providing a way to use the LSODA algorithm from Linda Petzold and Alan Hindmarsh from [Julia](http://julialang.org/). **Clang.jl** has been used to write the library.## InstallationTo install this package, run the command `Pkg.clone("https://github.com/rveltz/LSODA.jl.git")`## Simplified Functions```juliafunction rhs!(t, x, ydot, data)	ydot[1]=1.0E4 * x[2] * x[3] - .04E0 * x[1]	ydot[3]=3.0E7 * x[2] * x[2]	ydot[2]=-ydot[1] - ydot[3]  nothingendy0 = [1.,0.,0.]tspan = [0., 0.4]@time lsoda_0(rhs!, y0, tspan, reltol= 1e-4, abstol = Vector([1.e-6,1.e-10,1.e-6]))```This should give the following.```at t =   4.0000e-01 y=   9.851712e-01   3.386380e-05   1.479493e-02at t =   4.0000e+00 y=   9.055333e-01   2.240655e-05   9.444430e-02at t =   4.0000e+01 y=   7.158403e-01   9.186334e-06   2.841505e-01at t =   4.0000e+02 y=   4.505250e-01   3.222964e-06   5.494717e-01at t =   4.0000e+03 y=   1.831976e-01   8.941774e-07   8.168016e-01at t =   4.0000e+04 y=   3.898729e-02   1.621940e-07   9.610125e-01at t =   4.0000e+05 y=   4.936362e-03   1.984221e-08   9.950636e-01at t =   4.0000e+06 y=   5.161832e-04   2.065786e-09   9.994838e-01at t =   4.0000e+07 y=   5.179811e-05   2.072030e-10   9.999482e-01at t =   4.0000e+08 y=   5.283524e-06   2.113420e-11   9.999947e-01at t =   4.0000e+09 y=   4.658945e-07   1.863579e-12   9.999995e-01at t =   4.0000e+10 y=   1.423392e-08   5.693574e-14   1.000000e+00```
+[![Build Status](https://travis-ci.org/rveltz/LSODA.jl.svg?branch=master)](https://travis-ci.org/rveltz/LSODA)
+[![Coverage Status](https://coveralls.io/repos/github/rveltz/LSODA.jl/badge.svg?branch=master)](https://coveralls.io/github/rveltz/LSODA.jl?branch=master)
+[![Build status](https://ci.appveyor.com/api/projects/status/p879qachs4c52y4u/branch/master?svg=true)](https://ci.appveyor.com/project/rveltz/lsoda-jl/branch/master)
+
+
+
+
+# LSODA.jl## Introduction**LSODA.jl** is a Julia package that interfaces to the [liblsoda](https://github.com/sdwfrost/liblsoda) library, developped by [Simon Frost](http://www.vet.cam.ac.uk/directory/sdf22@cam.ac.uk) ([@sdwfrost](http://github.com/sdwfrost)), thereby providing a way to use the LSODA algorithm from Linda Petzold and Alan Hindmarsh from [Julia](http://julialang.org/). **Clang.jl** has been used to write the library and **Sundial.jl** was a inspiring source.## InstallationTo install this package, run the command `Pkg.clone("https://github.com/rveltz/LSODA.jl.git")`## Simplified FunctionsTo solve an ODE, one can call the simplified solver:```juliafunction rhs!(t, x, ydot, data)	ydot[1]=1.0E4 * x[2] * x[3] - .04E0 * x[1]	ydot[3]=3.0E7 * x[2] * x[2]	ydot[2]=-ydot[1] - ydot[3]  nothingendy0 = [1.,0.,0.]tspan = [0., 0.4]res =  lsoda(rhs!, y0, tspan, reltol= 1e-4, abstol = Vector([1.e-6,1.e-10,1.e-6]))```To reproduce the test example from liblsoda, on can use:```julialsoda_0(rhs!, y0, tspan, reltol= 1e-4, abstol = Vector([1.e-6,1.e-10,1.e-6]))```This should give the following.```at t =   4.0000e-01 y=   9.851712e-01   3.386380e-05   1.479493e-02at t =   4.0000e+00 y=   9.055333e-01   2.240655e-05   9.444430e-02at t =   4.0000e+01 y=   7.158403e-01   9.186334e-06   2.841505e-01at t =   4.0000e+02 y=   4.505250e-01   3.222964e-06   5.494717e-01at t =   4.0000e+03 y=   1.831976e-01   8.941774e-07   8.168016e-01at t =   4.0000e+04 y=   3.898729e-02   1.621940e-07   9.610125e-01at t =   4.0000e+05 y=   4.936362e-03   1.984221e-08   9.950636e-01at t =   4.0000e+06 y=   5.161832e-04   2.065786e-09   9.994838e-01at t =   4.0000e+07 y=   5.179811e-05   2.072030e-10   9.999482e-01at t =   4.0000e+08 y=   5.283524e-06   2.113420e-11   9.999947e-01at t =   4.0000e+09 y=   4.658945e-07   1.863579e-12   9.999995e-01at t =   4.0000e+10 y=   1.423392e-08   5.693574e-14   1.000000e+00```

diff --git a/src/LSODA.jl b/src/LSODA.jl
@@ -10,7 +10,7 @@ else
     error("LSODA is not properly installed. Please run Pkg.build(\"LSODA\")")
 end
 
-export lsoda, lsoda_0, lsoda_opt_t, lsoda_context_t, lsoda_prepare, lsoda_opt_t, lsoda_free
+export lsoda, lsoda_0, lsoda_opt_t, lsoda_context_t, lsoda_prepare, lsoda_opt_t, lsoda_free, lsoda_evolve!
 
 export LSODAAlgorithm, LSODAAlg, solve
 

diff --git a/src/solver.jl b/src/solver.jl
@@ -1,16 +1,11 @@
 function lsodafun{T1,T2,T3}(t::T1,y::T2,yp::T3,userfun::UserFunctionAndData)
-  y_ = pointer_to_array(y,userfun.neq)
-  ydot_ = pointer_to_array(yp,userfun.neq)
+  y_ = unsafe_wrap(Array,y,userfun.neq)
+  ydot_ = unsafe_wrap(Array,yp,userfun.neq)
   userfun.func(t, y_, ydot_,userfun.data)
   return Int32(0)
 end
 
-function lsodafun{T1,T2,T3}(t::T1,y::T2,yp::T3, userfun::Function)
-    y_ = pointer_to_array(y,3)
-    ydot_ = pointer_to_array(yp,3)
-    userfun(t, y_, ydot_)
-    return Int32(0)
-end
+const fex_c = cfunction(lsodafun,Cint,(Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ref{UserFunctionAndData}))
 
 function lsoda_0(f::Function, y0::Vector{Float64}, tspan::Vector{Float64}; userdata::Any=nothing, reltol::Union{Float64,Vector}=1e-4, abstol::Union{Float64,Vector}=1e-10)
   neq = Int32(length(y0))
@@ -34,9 +29,7 @@ function lsoda_0(f::Function, y0::Vector{Float64}, tspan::Vector{Float64}; userd
   t = Array{Float64}(1)
   tout = Array{Float64}(1)
   y = Array{Float64}(neq)
-  for i=1:neq
-    y[i] = y0[i]
-  end
+  y = copy(y0)
 
   t[1] = tspan[1]
   tout[1] = tspan[2]
@@ -48,8 +41,6 @@ function lsoda_0(f::Function, y0::Vector{Float64}, tspan::Vector{Float64}; userd
     opt.itask = 1
   #
 
-  const fex_c = cfunction(lsodafun,Cint,(Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ref{typeof(userfun)}))
-
   ctx = lsoda_context_t()
     ctx.function_ = fex_c
     ctx.neq = neq
@@ -68,7 +59,6 @@ end
 function lsoda(f::Function, y0::Vector{Float64}, tspan::Vector{Float64}; userdata::Any=nothing, reltol::Union{Float64,Vector}=1e-4, abstol::Union{Float64,Vector}=1e-10)
   neq = Int32(length(y0))
   userfun = UserFunctionAndData(f, userdata,neq)
-  # println("-->userfun = ",userfun,typeof(userfun))
 
   atol = ones(Float64,neq)
   rtol = ones(Float64,neq)
@@ -87,24 +77,19 @@ function lsoda(f::Function, y0::Vector{Float64}, tspan::Vector{Float64}; userdat
     rtol = copy(reltol)
   end
 
-  t = Array{Float64}(1)
+  t    = Array{Float64}(1)
   tout = Array{Float64}(1)
-  y = Array{Float64}(neq)
-  for i=1:neq
-    y[i] = y0[i]
-  end
+  y    = Array{Float64}(neq)
+  y = copy(y0)
 
-  t[1] = tspan[1]
+  t[1]    = tspan[1]
   tout[1] = tspan[2]
-  #
+
   opt = lsoda_opt_t()
     opt.ixpr = 0
     opt.rtol = pointer(rtol)
     opt.atol = pointer(atol)
     opt.itask = 1
-  #
-
-  const fex_c = cfunction(lsodafun,Cint,(Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ref{typeof(userfun)}))
 
   ctx = lsoda_context_t()
     ctx.function_ = fex_c
@@ -113,10 +98,24 @@ function lsoda(f::Function, y0::Vector{Float64}, tspan::Vector{Float64}; userdat
     ctx.data = pointer_from_objref(userfun)
 
   lsoda_prepare(ctx,opt)
-  for i=1:12
+  yres[1,:] = y0
+
+  for k in 2:length(tspan)
+	tout[1] = tspan[k]  
     lsoda(ctx,y,t,tout[1])
-    @printf("at t = %12.4e y= %14.6e %14.6e %14.6e\n",t[1],y[1], y[2], y[3])
-    tout[1] *= 10.0E0
+	yres[k,:] = copy(y)
   end
-  return ctx
-end
+  return ctx, yres
+end
+
+function lsoda_evolve!(ctx::lsoda_context_t,y::Vector{Float64},tspan::Vector{Float64};data=nothing)
+	if data != nothing
+		ctx.data.userdata = data
+	end
+	t    = Array{Float64}(1)
+	tout = Array{Float64}(1)
+	t[1] = tspan[1]
+	tout[1] = tspan[2]
+	lsoda(ctx,y,t,tout[1])
+end
+
diff --git a/test/runtests.jl b/test/runtests.jl
@@ -2,5 +2,7 @@ using LSODA
 using Base.Test
 
 # write your own tests here
+println("--> test1 =============")
 include("test1.jl")
+println("\n--> test2 =============")
 include("test2.jl")
diff --git a/test/test2.jl b/test/test2.jl
@@ -5,9 +5,27 @@ function rhs!(t, x, ydot, data)
   nothing
 end
 
-# rhs!(t, x, ydot) =  rhs!(t, x, ydot, nothing)
+rhs!(t, x, ydot) =  rhs!(t, x, ydot, nothing)
 
 y0 = [1.,0.,0.]
+println("\n####################################\n--> Use of a old wrapper for speed comparison")
 tspan = [0., 0.4]
 @time lsoda_0(rhs!, y0, tspan, reltol= 1e-4,abstol = Vector([1.e-6,1.e-10,1.e-6]))
 @time lsoda_0(rhs!, y0, tspan, reltol= 1e-4,abstol = Vector([1.e-6,1.e-10,1.e-6]))
+
+# case with a vector
+println("\n####################################\n--> Use of a vector of times where output is required")
+tspan = [4.*10.0^k for k=-1:10]
+_, res = lsoda(rhs!, y0, tspan, reltol= 1e-4,abstol = Vector([1.e-6,1.e-10,1.e-6]))
+_, res = @time lsoda(rhs!, y0, tspan, reltol= 1e-4,abstol = Vector([1.e-6,1.e-10,1.e-6]))
+println(res)
+
+#case where we don't have to declare a new context
+println("\n####################################\n--> Use of a lsoda_evolve!")
+y0 = [1.,0.,0.]
+tspan = [4.*10.0^k for k=-1:10]
+ctx, _ = lsoda(rhs!, y0, tspan[1:2], reltol= 1e-4,abstol = Vector([1.e-6,1.e-10,1.e-6]))
+@time for k=2:length(tspan)
+	lsoda_evolve!(ctx,y0,tspan[k-1:k])
+	@printf("at t = %12.4e y= %14.6e %14.6e %14.6e\n",tspan[k],y0[1], y0[2], y0[3])
+end