@Tomas: maybe check out Numerical Recipes in C: The Art of Scientific
Computing, 2nd edition. There is also an edition for Fortran. The code that
I use in C is basically from there.
@Andrew: The xprime needs to be in the loop. I just made it ones to
simplify but normally it changes every iteration. (In the DP problem, the
loop is calculating an expecation and xprime is the possible future value
of the state variable for each state of the world). Concerning the Dierckx
package. I don't know about the general behaviour but for my particular
problem (irregular grid + cubic spline) it is very slow. Run the following
code:
using Dierckx
spacing=1.5
Nxx = 300
Naa = 350
Nalal = 200
sigma = 10
NoIter = 10000
xx=Array(Float64,Nxx)
xmin = 0.01
xmax = 400
xx[1] = xmin
for i=2:Nxx
xx[i] = xx[i-1] + (xmax-xx[i-1])/((Nxx-i+1)^spacing)
end
f_util(c) = c.^(1-sigma)/(1-sigma)
W=Array(Float64,Nxx,1)
W[:,1] = f_util(xx)
spl = Spline1D(xx,W[:,1])
function performance2(NoIter::Int64)
W_temp = Array(Float64,Nalal*Naa)
W_temp2 = Array(Float64,Nalal,Naa)
xprime=Array(Float64,Nalal,Naa)
for banana=1:NoIter
xprime=ones(Nalal,Naa)
W_temp = spl(xprime[:])
end
W_temp2 = reshape(W_temp,Nalal,Naa)
end
@time performance2(10000)
30.878093 seconds (100.01 k allocations: 15.651 GB, 2.19% gc time)
That's why I went on and asked my friend to help me out in the first place.
I think the mnspline is really (not saying it's THE fastest) fast in doing
the interpolation itself (magnitudes faster than MATLAB). But then I just
don't understand how MATLAB can catch up by just looping through the same
operation over and over. Intuitively (maybe I'm wrong) it should be
somewhat proportional. If my code in Julia is 10 times faster whitin a
loop, and then I just repeat the operation in that particular loop very
often, how can it turn out to be only equally fast as MATLAB. Again, the
mnspline uses all my threads maybe it has something to do with overhead,
whatever. I don't know, hints appreciated.
