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