# Speed optimization

I have a part of a program that I have to repeatedly call many times (about 2000 to 10000 times)

f = Reap[Do[Sow[{0, 0, 0}], {i, 1, nparticle}];][[2, 1]];
vir = 0.;
hL = L/2.;
e = 0.;
Do[
Do[
dr = r[[i]] - r[[j]];
dr = Map[If[# > hL, # - L, If[# < -hL, # + L, #]] &, dr];
r2 = Dot[dr, dr];
If[r2 < rc2,
r6i = 1./(r2^3);
e += 4*(r6i*r6i - r6i) - ecut;
valuef = 48*(r6i*r6i - 0.5*r6i);
df = dr*valuef/r2;
f[[i]] += df;
f[[j]] -= df;
vir += valuef;
];
If[r2 <= rc^2,
gr[[IntegerPart[Sqrt[r2]/dh]]] += 2;
]
, {j, i + 1, nparticle}
];
, {i, 1, nparticle - 1}
];


r is a list of the form {{3,4,6},{5,7,8},.....}.
hL, L, rc2 are parameter I initialize at the start.

Is there any ways I can optimize the code like using compile or parallelization? I have tried these but always getting various kinds of errors or wrong calculations.
Is the two Do loops avoidable in this case?

As requested, here is the value of parameters:

rc = 2.5;
rc2 = rc^2;
nparticle = 108;
rho = 0.4;
L = (nparticle/rho)^(1./3.);
hL = 0.5*L;
rr3 = 1/(rc^3);
ecut = 4.*(rr3^4 - rr3^2);


The initialization of r:

LinCell = Round[(nparticle/4.)^(1./3.)];
lattconst = L/LinCell;
r = Reap[
Do[
Do[
Do[
Sow[{ix + 0.25, iy + 0.25, iz + 0.25}*lattconst];
Sow[{ix + 0.75, iy + 0.75, iz + 0.25}*lattconst];
Sow[{ix + 0.75, iy + 0.25, iz + 0.75}*lattconst];
Sow[{ix + 0.25, iy + 0.75, iz + 0.75}*lattconst];
, {iz, 0, LinCell - 1}
];
, {iy, 0, LinCell - 1}
];
, {ix, 0, LinCell - 1}
];
][[2, 1]];


--edit-- Now I need to add an output to the program

If[r2 <= rc^2,
gr[[IntegerPart[Sqrt[r2]/dh]]] += 2;
]


which is added to the inner the Do loop. (I have modified the above code)
gr is initialized to {0,0,0,...}.
dh=0.02

-
Since your r array is apparently unchanging, you could preprocess with Nearest and, for each i, avoid the j loop by explicitly requesting only those r elements within distance rc2 of r[[i]]. –  Daniel Lichtblau Jan 4 '13 at 16:03
Please include some sample data for all your parameters so that the code can be run as is rather than having to first create these. –  Mr.Wizard Jan 4 '13 at 16:10
ecut remains without an example/sample value –  Jacob Akkerboom Jan 4 '13 at 16:42
ecut value added –  Sheng Jan 4 '13 at 16:52
See the documentation of Do to see that/how you can put multiple iterators in one Do. Note also that the semicolons after each of your "Do's" is unnecessary. –  Jacob Akkerboom Jan 4 '13 at 17:13

I think this works. I'll slightly change the setup and also make a proper Module out of the main part.

rc = 2.5;
rc2 = rc^2;
nn = 6; (* Was 3 in original post *)
nparticle = 4*nn^3;
rho = 0.4;
L = (nparticle/rho)^(1./3.);
hL = 0.5*L;
rr3 = 1/(rc^3);
ecut = 4.*(rr3^4 - rr3^2);
LinCell = Round[(nparticle/4.)^(1./3.)];
lattconst = L/LinCell;
r = Flatten[Table[{{ix + 0.25, iy + 0.25, iz + 0.25},
{ix + 0.75, iy + 0.75, iz + 0.25},
{ix + 0.75, iy + 0.25, iz + 0.75},
{ix + 0.25, iy + 0.75, iz + 0.75}}*lattconst, {ix, 0,
LinCell - 1}, {iy, 0, LinCell - 1}, {iz, 0, LinCell - 1}], 3];
hL = L/2.;


Here is the original method.

func[r_] :=
Module[{e = 0., vir = 0., n = Length[r], f, dr, r2, r6i, valuef,
df},
f = ConstantArray[0., {n, 3}];
Do[
Do[
dr = r[[i]] - r[[j]];
dr = Map[If[# > hL, # - L, If[# < -hL, # + L, #]] &, dr];
r2 = Dot[dr, dr];
If[r2 < rc2,
r6i = 1./(r2^3);
e += 4*(r6i*r6i - r6i) - ecut;
valuef = 48*(r6i*r6i - 0.5*r6i);
df = dr*valuef/r2;
f[[i]] += df;
f[[j]] -= df;
If[df > 0, Print[df]];
vir += valuef;],
{j, i + 1, n}]
, {i, 1, n - 1}];
{e, vir, f}
]


Test:

Timing[{e, vir, f} = func[r];]

(* Out[837]= {6.750000, Null} *)


Since you use periodicity to dtermine proximity I enlarge the set by surrounding it with its 26 neighbors (allowing edge and vertex intersection as well as faces). We now look at every pair twice since I do not try to enforce an ordering (that would be more of a headache to code). In compensating for this, to get the numbers to work out I made a few other minor changes

func2[r_] := Module[
{e = 0., vir = 0., n = Length[r], nf, plusminus, enlargedr, f, dr,
r2, r6i, valuef, df, ri, near, fi},
plusminus = Map[Times[#, {L, L, L}] &, Tuples[{-1, 0, 1}, 3]];
enlargedr =
Flatten[Table[
Map[plusminus[[j]] + # &, r], {j, Length[plusminus]}], 1];
nf = Nearest[enlargedr];
f = Table[
fi = 0.;
ri = r[[i]];
near = Rest[nf[ri, {Infinity, rc}]];
Do[
dr = ri - near[[j]];
r2 = Dot[dr, dr];
r6i = 1./(r2^3);
e += 2*(r6i*r6i - r6i) - ecut/2;
valuef = 24*(r6i*r6i - 0.5*r6i);
df = dr*valuef/r2;
fi += 2 df;
vir += valuef;
, {j, Length[near]}];
fi
, {i, n}];
{e, vir, f}
]
Timing[{e2, vir2, f2} = func2[r];]

(* Out[839]= {0.400000, Null} *)


This seems to be in tolerable agreement with the basic method.

{e2 - e, vir2 - vir, Max[Abs[f2 - f]]}

(* Out[846]= {3.38332*10^-10, -1.63709*10^-10, 1.94289*10^-15} *)


If I change nn from 6 to 7, the first method takes 17.5 seconds and the second one takes 0.64 seconds.

--- edit ---

With modest work one can get this also to Compile, as with @P.Fonseca's method. There will be external calls for the NearestFunction lookups though.

ccC = Compile[{{r, _Real, 2}, {len, _Real}, {ecut, _Real}}, Module[
{e = 0., vir = 0., n = Length[r], f, dr, r2, r6i, valuef, df, ri,
near, fi},
f = Table[
fi = ConstantArray[0., 3];
ri = r[[i]];
near = Rest[nf[ri, {Infinity, rc}]];
Do[
dr = ri - near[[j]];
r2 = Dot[dr, dr];
r6i = 1./(r2^3);
e += 2*(r6i*r6i - r6i) - ecut/2;
valuef = 24*(r6i*r6i - 0.5*r6i);
df = dr*valuef/r2;
fi += 2 df;
vir += valuef;
, {j, Length[near]}];
fi
, {i, n}];
Join[{{e, vir, 0.}}, f]
], {{nf[__], _Real, 2}}];

func4[rr_, llen_, eecut_] := Module[{plusminus, enlargedr},
plusminus =
Map[Times[#, {llen, llen, llen}] &, Tuples[{-1, 0, 1}, 3]];
enlargedr =
Flatten[Table[
Map[plusminus[[j]] + # &, r], {j, Length[plusminus]}], 1];
nf = Nearest[enlargedr];
ccC[rr, llen, eecut]
]


With nn=9 this takes around .25 seconds whereas my func2 takes 1.3 or so. I guess if further speed is required one could make a LibraryFunction and compile to C. Maybe (I'm not really experienced with the LibraryFunction capabilities).

--- end edit ---

-
Why is there discrepancies between your results and the original algorithm? I notice that you have divide some parameters by 2, which should have compensated the act of looking at each pair twice. –  Sheng Jan 5 '13 at 16:22
My guess as to the discrepancies is cancellation error. You are summing values that, to close approximation, eventually give zero. In such situations it is not hard to have modest discrepancies based on order of operations. This is only a guess though and there may be something else I am missing. –  Daniel Lichtblau Jan 5 '13 at 20:33
I have modified the code as above. Can you can give suggestions on how to integrate it in your compile code? –  Sheng Jan 6 '13 at 2:45
Depends on the dimensions of gr (which I cannot really figure out). If it has the same length as f then you might hide it therein by adding a column to f. Alternatively just let it be a variable that is external to the Compile. That will cause evaluator calls to be made but that already happens with the NearestFunction so this probably will not slow things much. –  Daniel Lichtblau Jan 6 '13 at 20:43

I compiled your exact algorithm, and it seems to work OK.

I get 15 times speed up when using WVM as target, and 60 times when using C (CompilationTarget->"C").

The output is: {{e, vir, 0}, f}

test = Compile[{{nparticle, _Integer, 0}, {rho, _Real, 0}, {rc, _Real,
0}, {r, _Real, 2}, {f, _Real, 2}},
Module[{vir = 0., e = 0., dr = {0., 0., 0.}, r2 = 0., rc2 = 0.,
L = 0., hL = 0., rr3 = 0., r6i = 0., ecut = 0., valuef = 0.,
df = {0., 0., 0.}, i = 1, j = 1, faux = f},
rc2 = rc^2;
L = (nparticle/rho)^(1./3.);
hL = L*0.5;
rr3 = 1/(rc^3);
ecut = 4.*(rr3^4 - rr3^2);
Do[
dr = r[[i]] - r[[j]];
dr = Map[If[# > hL, # - L, If[# < -hL, # + L, #]] &, dr];
r2 = Dot[dr, dr];
If[r2 < rc2, r6i = 1./(r2^3);
e += 4*(r6i*r6i - r6i) - ecut;
valuef = 48*(r6i*r6i - 0.5*r6i);
df = dr*valuef/r2;
faux[[i]] += df;
faux[[j]] -= df;
vir += valuef;],
{i, 1, nparticle - 1}, {j, i + 1, nparticle}];
Join[{{e, vir, 0.}}, faux]
]
];

rc = 2.5;
nparticle = 108;
rho = 0.4;
LinCell = Round[(nparticle/4.)^(1./3.)];
lattconst = L/LinCell;
r = Reap[Do[Do[Do[Sow[{ix + 0.25, iy + 0.25, iz + 0.25}*lattconst];
Sow[{ix + 0.75, iy + 0.75, iz + 0.25}*lattconst];
Sow[{ix + 0.75, iy + 0.25, iz + 0.75}*lattconst];
Sow[{ix + 0.25, iy + 0.75, iz + 0.75}*lattconst];, {iz, 0,
LinCell - 1}];, {iy, 0, LinCell - 1}];, {ix, 0,
LinCell - 1}];][[2, 1]];
f = Reap[Do[Sow[{0, 0, 0}], {i, 1, nparticle}];][[2, 1]];

test[nparticle, rho, rc, r, f]


-------EDIT--------

I don't know if I correctly understood your edited request. Nevertheless, here is a different approach to the output that I think correctly extracts gr information, and that keeps more or less the same velocity gain:

test2 = Compile[{{nparticle, _Integer, 0}, {rho, _Real, 0}, {rc, _Real,
0}, {r, _Real, 2}, {f, _Real, 2}},
Module[{vir = 0., e = 0., dr = {0., 0., 0.}, r2 = 0., rc2 = 0.,
L = 0., hL = 0., rr3 = 0., r6i = 0., ecut = 0., valuef = 0.,
df = {0., 0., 0.}, i = 1, j = 1, faux = f, dh = 0.02, gr = {}},
rc2 = rc^2;
L = (nparticle/rho)^(1./3.);
hL = L*0.5;
rr3 = 1/(rc^3);
ecut = 4.*(rr3^4 - rr3^2);
Do[
dr = r[[i]] - r[[j]];
dr = Map[If[# > hL, # - L, If[# < -hL, # + L, #]] &, dr];
r2 = Dot[dr, dr];
If[r2 < rc2, r6i = 1./(r2^3);
e += 4*(r6i*r6i - r6i) - ecut;
valuef = 48*(r6i*r6i - 0.5*r6i);
df = dr*valuef/r2;
faux[[i]] += df;
faux[[j]] -= df;
vir += valuef;];
If[r2 <= rc^2, AppendTo[gr, IntegerPart[Sqrt[r2]/dh]]],
{i, 1, nparticle - 1}, {j, i + 1, nparticle}];
Join[{e, vir}, Flatten[faux], gr]
]
];

rc = 2.5;
nparticle = 108;
rho = 0.4;
LinCell = Round[(nparticle/4.)^(1./3.)];
lattconst = L/LinCell;
r = Reap[Do[Do[Do[Sow[{ix + 0.25, iy + 0.25, iz + 0.25}*lattconst];
Sow[{ix + 0.75, iy + 0.75, iz + 0.25}*lattconst];
Sow[{ix + 0.75, iy + 0.25, iz + 0.75}*lattconst];
Sow[{ix + 0.25, iy + 0.75, iz + 0.75}*lattconst];, {iz, 0,
LinCell - 1}];, {iy, 0, LinCell - 1}];, {ix, 0,
LinCell - 1}];][[2, 1]];
f = Reap[Do[Sow[{0, 0, 0}], {i, 1, nparticle}];][[2, 1]];

ansModified = test2[nparticle, rho, rc, r, f];

eAns = ansModified[[1]]
virAns = ansModified[[2]]
fAns = Partition[ansModified[[3 ;; nparticle*3 + 2]], 3]
If[Length[ansModified] == nparticle*3 + 2,
grAns = {},
grAns = {1, 2}*# & /@ Tally[ansModified[[nparticle*3 + 3 ;; -1]]]
]


• in Module variables: dh = 0.02, gr = {}

• your gr edit slightly modified: If[r2 <= rc^2, AppendTo[gr, IntegerPart[Sqrt[r2]/dh]]]

• a Flattened compile output: Join[{e, vir}, Flatten[faux], gr]

• an output post-processing (that is fast enough): (see eAns, virAns, fAns, grAns)

-
Ah, this seems really nice. I thought Compile would usually not be able to handle Do loops. –  Jacob Akkerboom Jan 4 '13 at 18:41
+1... looks good to me! Your compiled code is absolutely fine. –  Oleksandr R. Jan 4 '13 at 18:49
@JacobAkkerboom Compile works best with Do loops actually :) –  acl Jan 4 '13 at 18:52
Thank you all. I think my new avatar is working better :). For some reason, my minGW stopped working (and I'm currently with no time to investigate why). If someone could update my answer with the gain of using C, I would really appreciate it. –  P. Fonseca Jan 4 '13 at 18:56
Perhaps Compile worked less smoothly in 8.0.0, or equally likely, I do not understand well at all how Compile works. Is it possibly still a good idea to Reap and Sow in a compiled function? Also: Does the Module help with speed or is it just considered good practise? I should really practise this compiling stuff. –  Jacob Akkerboom Jan 4 '13 at 19:05

1. Do not start user symbols (variable names) with capital letters. These often may conflict with built-in functions. You have a full selection of Script, Gothic, Double-struck, and Greek characters to choose from to avoid these. Look at Palettes > Special Characters to see these, and hover the mouse pointer over any one of them to see the fast entry form using Esc.

2. Try to avoid explicit loops in Mathematica whenever possible. One exception is when you intend to compile to code and need maximum speed, but IMHO you should always try to get a fast and concise implementation in the native language first.

3. Use vector and array operations when possible. If a function has the Listable attribute try to make use of it.

4. Learn where to use Sow and Reap, and specifically for you where not to use them. These are powerful and useful functions, but they are usually best used when you need to gather results in an irregular fashion such as gathering only certain elements conditionally.

5. When asking for help provide complete executable code whenever possible. Also, please describe what the code does and what its ultimate application is; this is often significant in both fully understanding the code and recognizing superior approaches.

Here is an improvement upon part of your code. I realize this is the ancillary code and not the focus of your question but I have not dug into that yet. I expect Daniel is on the right track however as he usually is.

I hope that this will help you to see an alternative approach. I realize the code may be rather opaque for a new Mathematica user and I will be happy to explain any part of it you ask about but I encourage you to make use of the documentation (easily accessed by highlighting a function name or operator and pressing F1) as well.

rc = 2.5;
rc2 = rc^2;
nparticle = 108;
rho = 0.4;
\[ScriptCapitalL] = (nparticle/rho)^(1./3.);
hL = 0.5*\[ScriptCapitalL];
rr3 = 1/(rc^3);
ecut = 4.*(rr3^4 - rr3^2);

linCell = Round[(nparticle/4.)^(1./3.)];
lattconst = \[ScriptCapitalL]/linCell;

With[{nums = {
{0.25, 0.25, 0.25},
{0.75, 0.75, 0.25},
{0.75, 0.25, 0.75},
{0.25, 0.75, 0.75}}
},
r = Flatten[
Tuples /@ (lattconst * Outer[Plus, nums, Range[0, linCell - 1]]),
{2, 1}
]
];

f = ConstantArray[0, {nparticle, 3}];


Specifically look at Outer, Tuples, ConstantArray, and Flatten.

I shall try to revisit the body of the problem later but I expect it may already be well addressed by others before then.

-