Re: Mathematica 20x slower than Java at arithmetic/special functions, is

*To*: mathgroup at smc.vnet.net*Subject*: [mg115889] Re: Mathematica 20x slower than Java at arithmetic/special functions, is*From*: Oliver Ruebenkoenig <ruebenko at wolfram.com>*Date*: Mon, 24 Jan 2011 05:23:38 -0500 (EST)

On Sun, 23 Jan 2011, Leo Alekseyev wrote: > I was playing around with JLink the other day, and noticed that Java > seems to outperform Mathematica by ~20x in an area where I'd expect > Mathematica to be rather well optimized -- arithmetic involving special > functions. In my particular example, I am simply evaluating a sum of > Bessel functions. I understand that much depends on the underlying > implementation, but I just want to run this by Mathgroup to see if > this is to be expected, or maybe if I'm doing something suboptimal in > Mathematica. Here's the code that I'm running: > > grid1dc[x_, > y_] = (With[{d = 0.1, NN = 50}, > Sum[Re[N[ > d BesselJ[1, 2 Pi d Sqrt[m^2 + n^2]]/ > Sqrt[m^2 + n^2 + 10^-7]] Exp[ > I 2.0 Pi (m x + n y)]], {m, -NN, NN, 1}, {n, -NN, NN, 1}]]) // > N > > and > > gridres1da = > With[{delta = 0.5, xlim = 2.5, ylim = 2.5}, > Table[{x, y, grid1dc[x, y]}, {x, -xlim, xlim, delta}, {y, -ylim, > ylim, delta}]] > > > Java implementation uses Colt and Apache common math libraries for the > Bessels and complex numbers, uses a double for loop, and consistently > runs 15-20 times faster. > > --Leo > > Leo, here is a speedup of about 100 and if you have multiple CPUs then there is an additional factor: tmp = (With[{d = 0.1, NN = 50}, Sum[Re[N[ d BesselJ[1, 2 Pi d Sqrt[m^2 + n^2]]/ Sqrt[m^2 + n^2 + 10^-7]] Exp[ I 2.0 Pi (m x + n y)]], {m, -NN, NN, 1}, {n, -NN, NN, 1}]]) // N // Chop; (* since tmp is a large expr, compilation to C takes a long time and does not pay *) cfGrid1D = With[{code = tmp}, Compile[{{x, _Real, 0}, {y, _Real, 0}}, code, RuntimeAttributes -> Listable, Parallelization -> True]]; gridData = Transpose[ Flatten[With[{delta = 0.5, xlim = 2.5, ylim = 2.5}, Developer`ToPackedArray@ Table[{x, y}, {x, -xlim, xlim, delta}, {y, -ylim, ylim, delta}]], 1]]; AbsoluteTiming[re = cfGrid1D @@ gridData;] Norm[ Flatten[gridres1da, 1][[All, -1]] - re] 1.88052*10^-13 Hope this helps, Oliver