Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Matrix multiplication in JIT compile #28

Open
mbkoltai opened this issue Nov 28, 2020 · 1 comment
Open

Matrix multiplication in JIT compile #28

mbkoltai opened this issue Nov 28, 2020 · 1 comment

Comments

@mbkoltai
Copy link

mbkoltai commented Nov 28, 2020

Hi,
I'm trying to do vector equations with diffeqr+jit_compile:
de_jul <- diffeqr::diffeq_setup()
p <- c(10,28,8/3); u0 <- c(1,0.2,0.3); tspan <- c(0,50); p <- c(10,28,8/3); K_lor=matrix(c(-p[1],p[1],0,p[2],-1,0,0,0,-p[3]),3,3)
f_ode_julia_vect <- function(u,x,t) {du=x%*%u + c(0,-u[1]*u[3],u[1]*u[2]); return(du)}
prob<-de_jul$ODEProblem(f_ode_julia_vect, u0,tspan,K_lor); fastprob=diffeqr::jitoptimize_ode(de_jul,prob)

but I get the following error:
Error in x %*% u : requires numeric/complex matrix/vector arguments
Trying without jit:
prob <- de_jul$ODEProblem(f_ode_julia_vect, u0, tspan, K_lor); sol <- de_jul$solve(prob)
produces the error:
MethodError: no method matching Array{Float64,1}(::Array{Float64,2})
Closest candidates are:
Array{Float64,1}(::AbstractArray{S,N}) where {T, N, S} at array.jl:562
Array{Float64,1}() where T at boot.jl:425
Array{Float64,1}(!Matched::UndefInitializer, !Matched::Int64) where T at boot.jl:406

However if I input the matrix but don't do a matrix multiplication:
f_ode_julia_vect<-function(u,x,t){du=c(x[1,1]*u[1]+x[1,2]*u[2],x[2,1]*u[1]+x[2,2]*u[2],x[3,3]*u[3])+c(0,-u[1]*u[3],u[1]*u[2]); return(du)}
then it works, but since I work with relatively large ODE systems that are much more convenient to write as vector equations this solution doesn't really help.
Thanks a lot for any help.

@mbkoltai
Copy link
Author

update:
without JIT
f_ode_julia_vect <- function(u,x,t) {du=c(x %*% u) + c(0,-u[1]*u[3],u[1]*u[2]); return(du)}
prob<-de_jul$ODEProblem(f_ode_julia_vect, u0,tspan,K_lor)

works (so writing the matrix product as c(x %*% u)), but not with JIT compiler.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant