Hello, I have decided to take a leap and use Julia instead of Matlab, hence I have converted a previous Matlab script in Julia.
The only thing the script does is to build a Jacobian matrix of a map. Everything is vectorized (no loop, only matrix operation: repeat / .(operator) / sparse.) and yet it takes at least twice the time in Julia vs Matlab (using @elpased function(…)).
This is the code:
module jac
using SparseArrays
using LinearAlgebra
export DF
function DF(x,m0,m,lambda,n,l)
#############################
# masses on the ith circle #
#############################
theta = repeat(360/l*collect(0:1:l1),1,n)
xi_ext = repeat(x,1,l)'
m_ext = repeat(m,1,l)'
A =  m_ext
B = similar(Array{Float64},axes(xi_ext))
B .= sqrt(2)*(xi_ext.^3).*(ones(l,n)  cosd.(theta)).^(1/2)
@views A[1,:] = zeros(n)
@views B[1,:] = ones(n)
# combine
DF = sparse(collect(1:1:n),collect(1:1:n),vec(lambda*ones(n,1)2*m0./x.^3 + sum(A./B,dims=1)'))
#############################
# masses on the jth circle #
#############################
# angles associated with each mass
theta = 360/l*collect(0:1:l1)
theta = repeat(theta,1,(n1)*n) # extend theta
# extend xi
# l=0 > [ x_1 ... x_1  ...  x_n ... x_n ]
# l=1 > [ x_1 ... x_1  ...  x_n ... x_n ]
# . [ .  ...  . ]
# l=l1 > [ x_1 ... x_1  ...  x_n ... x_n ]
xi_ext = repeat(reshape(repeat(x,1,n1)',n*(n1),1),1,l)'
# extend xj
# l=0 > [ x_2 x_3 ... x_n  ...  x_1 x_2 x_3 ... x_{n1} ]
# l=1 > [ x_2 x_3 ... x_n  ...  x_1 x_2 x_3 ... x_{n1} ]
# . [ .  ...  . ]
# l=l1 > [ x_2 x_3 ... x_n  ...  x_1 x_2 x_3 ... x_{n1} ]
xj_ext = repeat(x,n,l)'
@views xj_ext = xj_ext[:,setdiff(1:end, (0:1:n1)+(1:n:n*n))]
# extend mj the same way as xj
mj_ext = repeat(m,n,l)'
@views mj_ext = mj_ext[:,setdiff(1:end, (0:1:n1)+(1:n:n*n))]
##### Computate the derivatives
xi_carre = xi_ext.^2
xj_carre = xj_ext.^2
prod_xixjcos = xi_ext.*xj_ext.*cosd.(theta)
# diagonal elements of DF
A, B = similar(Array{Float64},axes(xi_ext)), similar(Array{Float64},axes(xi_ext))
A .= mj_ext.*(4*xi_carre + xj_carre  8*prod_xixjcos + 3*xj_carre.*cosd.(2*theta))
B .= 2*(xi_carre + xj_carre  2*prod_xixjcos).^(5/2)
sum_j=sum(reshape(sum(A./B,dims=1),n1,n),dims=1)
# nondiagonal elements of DF
A .= mj_ext.*(4*(xi_carre + xj_carre).*cosd.(theta) + xi_ext.*xj_ext.*(7*ones(l,(n1)*n) + cosd.(2*theta)))
sum_theta = sum(A./B,dims=1)
# combine
return DF + sparse(collect(1:1:n),collect(1:1:n),vec(sum_j)) + [tril(reshape(sum_theta,n1,n)',1) zeros(n,1)]+[zeros(n,1) triu(reshape(sum_theta,n1,n)')]
end
end
Three possible issues:

the way I run the code is not efficient. What I did: I created a module, say X, in a .jl. This module contains the function creating the Jacobin matrix. Then, I created a main.jl file which executes the function using include("/path/to/MODULE/X.jl") and X.FUNCTION_NAME(…). I launch from the Julia interface the script main.jl with include("/path/to/SCRIPT/X.jl").

the functions repeat / sparse are not as efficient in Julia.

it has to do with “Multithreading / parallelism computing” but I still don’t understand how to use/implement it in my code.
Could anyone help me understand what I have done wrong in the code?
Thank you for your help!