Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Re: Compilation: Avoiding inlining

• To: mathgroup at smc.vnet.net
• Subject: [mg121459] Re: Compilation: Avoiding inlining
• From: Oliver Ruebenkoenig <ruebenko at wolfram.com>
• Date: Fri, 16 Sep 2011 05:46:21 -0400 (EDT)
• Delivered-to: l-mathgroup@mail-archive0.wolfram.com
• References: <201109150841.EAA29262@smc.vnet.net>

On Thu, 15 Sep 2011, DmitryG wrote:

> This message refers to the extensive exchange here:
>
> http://groups.google.com/group/comp.soft-sys.math.mathematica/browse_thread/thread/4ea9938ad1d523c6/735f864eab188778?lnk=gst&q=dmitryg#735f864eab188778
>
> The problem is to solve a system of N ordinary differential equations
> (ODE) with Runge-Kutta-4 method with compilation. Here is an example
> of working program
>
> ***************************************************************************
> (* Runge-Kutta-4 routine *)
> makeCompRK[fIn_] :=
>  With[{f = fIn},
>   Compile[{{x0, _Real,
>      1}, {t0, _Real}, {tMax, _Real}, {n, _Integer}},
>    Module[{h, K1, K2, K3, K4, SolList, x = x0, t},
>     h = (tMax - t0)/n;
>     SolList = Table[x0, {n + 1}];
>     Do[t = t0 + k h;
>      K1 = h f[t, x];
>      K2 = h f[t + (1/2) h, x + (1/2) K1];
>      K3 = h f[t + (1/2) h, x + (1/2) K2];
>      K4 = h f[t + h, x + K3];
>      x = x + (1/6) K1 + (1/3) K2 + (1/3) K3 + (1/6) K4;
>      SolList[[k + 1]] = x
>      , {k, 1, n}];
>     SolList],(*CompilationTarget->"C"*)]];
>
> (* Definition of the equations *)
> NN = 100;
> Timing[
> RHS = Quiet[
>    Table[-x[[i]] Sin[
>        0.3 t]^2/(1 + 300 Sum[x[[j]], {j, 1, NN}]^2/NN^2), {i, 1,
> NN}]];] (* Quiet to suppress complaints *)
> F = Function[{t, x},
>   Evaluate[RHS]];  (* Evaluate to use actual form of RHS *)
> Print[]
>
> (* Compilation *)
> tt0 = AbsoluteTime[];
> Timing[RK4Comp = makeCompRK[F];]
> AbsoluteTime[] - tt0
> Print[];
>
> (* Setting parameters and Calculation *)
> x0 = Table[
>  RandomReal[{0, 1}], {i, 1, NN}];   t0 = 0;    tMax = 100;   n = 500;
> tt0 = AbsoluteTime[];
> Sol = RK4Comp[x0, t0, tMax, n];
> AbsoluteTime[] - tt0
>
> Print["Compilation: ", Developer`PackedArrayQ@Sol]
>
> (* Plotting time dependences *)
> x1List = Table[{1. t0 + (tMax - t0) k/n, Sol[[k + 1, 1]]}, {k, 0,
> n}];
> x2List = Table[{1. t0 + (tMax - t0) k/n, Sol[[k + 1, 2]]}, {k, 0, n}];
> x3List = Table[{1. t0 + (tMax - t0) k/n, Sol[[k + 1, 3]]}, {k, 0, n}];
> ListLinePlot[{x1List, x2List, x3List}, PlotMarkers -> Automatic,
> PlotStyle -> {Blue, Green, Red}, PlotRange -> {0, 1}]
>
> **********************************************************************************
>
> The drawback of this solution is what Oliver Ruebenkoenig called
> "inlining". For large NN the RHS of the equations is a lot the
> explicit code generated by Mathematica. This code then gets compiled
> and the size of the compiled code strongly increases with NN (it can
> be seen by the CompilePrint[RK4Comp] command). That is, we write a lot
> of code inline instead of letting it be generated in the course of the
> execution, with the compiled code being small.
>
> Mathematica's own compiler tolerates inlining, although the large code
> eats up the memory. If we uncomment the option "CompilationTarget-
>> "C"", we will see that for large NN compilation goes for ages. That
> is, inlining has to be avoided.
>
> Daniel Lichtblau has rewritten my similar program in the vectorized
> form that avoided inlining and could be compiled to C, and was fast.
> But vectorization is not always possible or easy. Oliver wrote that
> inlining can also be avoided by defining a function and then calling
> it, but there were no examples at that time. Now I am interested in
> the cases such as the program above where vectorization is unwanted.
> How to avoid inlining by defining a function in this case? I have
> tried to replace
>
> RHS = Quiet[Table[-x[[i]] Sin[0.3 t]^2/(1 + 300 Sum[x[[j]], {j, 1,
> NN}]^2/NN^2), {i, 1, NN}]];
>
> by
>
> FRHS[i_][x_List, t_] = -x[[i]] Sin[0.3 t]^2/(1 + 100 Sum[x[[i + j]],
> {j, 1, Min[3, NN - i]}]^2);
> RHS = Quiet[Table[FRHS[i][x, t], {i, 1, NN}]];
>
> This works correctly but does not compile, and the execution time is
> longer. Something is wrong here. Oliver, Daniel, and others, please,
> help!
>
> Thank you,
>
> Dmitry
>
>

Dmitry,

here is an idea that hopefully will get you started:

ClearAll[makeCompRK]
makeCompRK[fIn_] :=
   With[{f = fIn},
    Compile[{{x0, _Real,
       1}, {t0, _Real}, {tMax, _Real}, {n, _Integer}},
     Module[{h, K1, K2, K3, K4, SolList, x = x0, t},
      h = (tMax - t0)/n;
      SolList = Table[x0, {n + 1}];
      Do[t = t0 + k h;
       K1 = h f[t, x];
       K2 = h f[t + (1/2) h, x + (1/2) K1];
       K3 = h f[t + (1/2) h, x + (1/2) K2];
       K4 = h f[t + h, x + K3];
       x = x + (1/6) K1 + (1/3) K2 + (1/3) K3 + (1/6) K4;
       SolList[[k + 1]] = x, {k, 1, n}];
      SolList],(*CompilationTarget->"C",*)
     CompilationOptions -> {"InlineCompiledFunctions" -> True}]];

(* Have a look at InlineCompiledFunctions *)

(*Definition of the equations*)
NN = 100;
Timing[RHS =
    Quiet[Table[-x[[i]] Sin[0.3 t]^2/(1 +
          300 Sum[x[[j]], {j, 1, NN}]^2/NN^2), {i, 1,
       NN}]];] (*Quiet to suppress complaints*)

F = Function[{t, x}, Evaluate[RHS]];

>> 0.06


(*Compilation*)
tt0 = AbsoluteTime[];
Timing[RK4Comp = makeCompRK[F];]
AbsoluteTime[] - tt0

>> 0.2

cRHS = Compile[{{t, _Real, 0}, {x, _Real,
     1}}, -x*(Sin[0.3*t]^2/(1 + 300*Total[x]^2/Length[x]^2))]

(*cRHS[0.3,{1,2,3,4}//N]*)

(*Compilation*)
tt0 = AbsoluteTime[];
Timing[RK4Comp2 = makeCompRK[cRHS];]
AbsoluteTime[] - tt0

>> 0.002

(*Needs["CompiledFunctionTools`"]*)
(*CompilePrint[RK4Comp2]*)

(*
Have a look at how "InlineCompiledFunctions" -> True/False influence the 
output
*)

(*Setting parameters and Calculation*)
x0 = Table[RandomReal[{0, 1}], {i, 1, NN}]; t0 = 0; tMax = 100; n = 500;

tt0 = AbsoluteTime[];
Sol1 = RK4Comp[x0, t0, tMax, n];
AbsoluteTime[] - tt0

>> 0.43


tt0 = AbsoluteTime[];
Sol2 = RK4Comp2[x0, t0, tMax, n];
AbsoluteTime[] - tt0

>> 0.007

Norm[Sol1 - Sol2]
>> 10^-14

Sol = Sol2;
Print["Compilation: ", Developer`PackedArrayQ@Sol]

Some further notes:

- I did not try the ->"C" version, but I believe both 
compiles need to be compiled to C to be effective. Have a play. I'd be 
interested to about the results

- I understand that you wanted an non vectorized version. That in 
principal would work in the same way; you could provide an example of such 
a non-vectorizable function.

- I think this is a good illustration of why compilation to C is not 
always beneficial. It is only beneficial if the compile time vs the runtime 
is small.

Hope this gets you started.

Oliver

• References:
  • Compilation: Avoiding inlining
    • From: DmitryG <einschlag@gmail.com>