> > Numerical Recipes in C: The Art of Scientific Computing, 2nd edition.
Please note that code derived from this book cannot be included in BSD, LGPL, or GPL licensed libraries (which is to say, most Julia packages). Distribution is restricted to compiled binaries only, with commercial and non-commercial license variants. Source redistribution is prohibited. Numerical Recipes should not be relied on for any implementations that are submitted to Julia core or the open source Julia package ecosystem. On Sat, Jan 30, 2016 at 8:59 AM, <pokerhontas...@googlemail.com> wrote: > @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. > >
