Speed optimization

Question

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, 2013 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, 2013 at 16:10
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, 2013 at 17:13

Daniel Lichtblau · Accepted Answer · 2013-01-04 21:08:37Z

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. — Gosere, Jan 5, 2013 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, 2013 at 20:33
I have modified the code as above. Can you can give suggestions on how to integrate it in your compile code? — Gosere, Jan 6, 2013 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, 2013 at 20:43

P. Fonseca · Accepted Answer · 2013-01-06 13:00:30Z

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]]]
  ]

To help you out, the modified parts are:

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, 2013 at 18:41
+1... looks good to me! Your compiled code is absolutely fine. — Oleksandr R., Jan 4, 2013 at 18:49
@JacobAkkerboom Compile works best with Do loops actually :) — acl, Jan 4, 2013 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, 2013 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, 2013 at 19:05

Mr.Wizard · Accepted Answer · 2013-01-04 20:34:51Z

Let me offer some advice.

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.
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.
Use vector and array operations when possible. If a function has the Listable attribute try to make use of it.
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.
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.

Stack Exchange Network

Speed optimization

3 Answers 3

Your Answer

Not the answer you're looking for? Browse other questions tagged
performance-tuning
code-review
or ask your own question.

Hot Network Questions

Speed optimization

3 Answers 3

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged performance-tuningcode-review or ask your own question.

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
performance-tuning
code-review
or ask your own question.