Optimizing operations between vectors

Question

Basically I am generating a grid in space (so a vector where v[[1]]=position 1 in x) and a random number of particles in space (so a vector where p[[1]]= position of particle 1 on space).

Afterwards, I want to see in which spaces of the grid the particles are and do somethings.

The Do is the specific part of the code I want to optimize.

rho = Compile[{{xp, _Real, 1}, {xg, _Real, 
1}, {carga, _Real}, {np, _Real}, {ng, _Real}, {dx, _Real}},


n = ng + 1;
  ρ = Table[carga*np/(dx*ng), n];
  ρ[[ng + 1]] = 0;
  
  Do[If[Abs[xp[[i]] - xg[[j]]] <= 
     dx, ρ[[j]] = ρ[[j]] - 
      carga*(1 - (Abs[xp[[i]] - xg[[j]]])/dx)/dx; r += 1, 0] , {i, 1, 
    np}, {j, 1, ng + 1}];
  ρ[[1]] = ρ[[1]] + ρ[[ng + 1]];]

Basically I am worried that I am thinking too much in terms of other languages and maybe mathematica has a better way to do this instead of using

If[Abs[xp[[i]] - xg[[j]]] <= 
     dx

That manually makes a lot of operations

Here is a part of the code I am using

eps = 8.85*10^(-12);
q = 80;
m = 50;
L = 2;
nparticles = 2000;
ngrid = 200;
vth = 2;
u = 1;

dx = L/ngrid;

x = RandomVariate[UniformDistribution[{0, L}], nparticles];
r = 0;

xgrid = Range[0, L, dx] // N;
rho[x, xgrid, q, nparticles, ngrid, dx]; // AbsoluteTiming
r

Answer 1

This is how I would write your function:

rho2 = Compile[{{p, _Real, 1}, {g, _Real, 
     1}, {c, _Real}, {np, _Integer}, {ng, _Integer}, {dx, _Real}},
   Block[{\[Rho], init, pdiv, pmod, j, a, dxinv},
    init = N[np]/N[ng];
    dxinv = 1./dx;
    a = c dxinv;
    \[Rho] = Table[0., {ng + 1}];
    
    Do[
     pdiv = Compile`GetElement[p, i] dxinv;
     j = Ceiling[pdiv];
     pmod = pdiv - j;
     \[Rho][[j]] -= Abs[pmod];
     \[Rho][[j + 1]] -= Abs[pmod + 1];
     ,
     {i, 1, np}];
    \[Rho][[1]] += \[Rho][[ng + 1]];
    Do[\[Rho][[j]] = a (Compile`GetElement[\[Rho], j] + init), {j, 1, ng}];
    \[Rho][[-1]] *= a;
    \[Rho]
    ],
   CompilationTarget -> "C",
   RuntimeOptions -> "Speed"
   ];

Note that I scoped \[Rho] within the function and handed it over as output of the function.

Here is how it compares to the orginal function rho with a substantially finer grid and more particles:

q = 80;
m = 50;
L = 2;
nparticles = 10000;
ngrid = 1000;

dx = L/ngrid;
p = RandomVariate[UniformDistribution[{0, L}], nparticles];
g = Subdivide[0., L, ngrid];

r = 0;
t1 = First@RepeatedTiming[
    rho[p, g, q, nparticles, ngrid, dx];
    ];
t2 = First@RepeatedTiming[
    \[Rho]2 = rho2[p, g, q, nparticles, ngrid, dx];
    ];
"Relative error" -> Max[Abs[(1 - \[Rho]2/\[Rho])]]
"Speedup" -> t1/t2

"Relative error" -> 4.39937*10^-11

"Speedup" -> 12270.6

Note that the relative error is quite high (I would have expected something of the order of 10^-15). And it increases with ngrid. I think this is due to a serious precision issue in your code: You initialize \[Rho] with values that are orders of magnitude greater than the ones you subtract from it. This is why I initialize \[Rho] by 0., do the subtraction loop, and add the initial values afterwards.

The key idea here is what Daniel Huber suggested: Instead of comparing each particle to all grid points, we can just compute the integers belonging to the two grid points next to the particle. This reduces the original computational complexity from $O(\mathrm{ngrid} \times \mathrm{nparticles})$ to $O(\mathrm{nparticles} + \mathrm{ngrid})$. So the speedup also grows with ngrid.

The remaining changes are just elementary refactoring:

• Instead of dividing by dx so often, I compute its reciprocal only once and multiply with it. (Floating point division is several times slower than multiplication.)

• Also, I postpone the multiplication with a until after the contibutions of all particles have been accumulated into \[Rho].

Answer 2

Instead of creating a space grid, we simply imply integer coordinates from 0 to n. E.g. {0,0,0} indicates the cell at the origin. {1,0,0} the next cell along the x axis etc. Then we can get the cell number for (nearly) "free":

To create some random points:

n = 100;
npts = 10;
SeedRandom[1];
pts = RandomReal[{0, n}, {npts, 3}]

(* {{81.7389, 11.142, 78.9526}, {18.7803, 24.1361, 6.57388}, {54.2247, 23.1155, 39.6006}, {70.0474, 21.1826, 74.8657}, {42.2851, 24.7495, 97.7172}, {82.5163, 92.5275, 57.8056}, {29.287, 20.8051, 58.0474}, {12.8821, 30.6427, 71.2012}, {39.0582, 81.9967, 32.5351}, {59.326, 51.8774, 16.9013}} *)

Now, to get the cells in wich the random points are, we only need to calculate the floor of the coordinates:

Floor[pts]

(* {{81, 11, 78}, {18, 24, 6}, {54, 23, 39}, {70, 21, 74}, {42, 24, 97}, {82, 92, 57}, {29, 20, 58}, {12, 30, 71}, {39, 81, 32}, {59, 51, 16}} *)

Answer 3

A program written using Compile will only run faster if your code does not violate it's limitations. For an explanation of what those limitations are read this.