Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: Compilation: Avoiding inlining

• To: mathgroup at smc.vnet.net
• Subject: [mg121595] Re: Compilation: Avoiding inlining
• From: DmitryG <einschlag at gmail.com>
• Date: Wed, 21 Sep 2011 05:34:56 -0400 (EDT)
• Delivered-to: l-mathgroup@mail-archive0.wolfram.com
• References: <201109160946.FAA12301@smc.vnet.net> <j4vana$dal$1@smc.vnet.net> <j59omg$l9$1@smc.vnet.net>

On Sep 20, 6:08 am, David Bailey <d... at removedbailey.co.uk> wrote:
> On 16/09/2011 12:08, Oliver Ruebenkoenig wrote:
>
> > On Fri, 16 Sep 2011, DmitryG wrote:
>
> >> Here is a program with (* Definition of the equations *) made one-
> >> step, that performs the same as in my previous post....................
>
> I tried pasting your example into Mathematica, but unfortunately there
> seems to be a variable 'x' which is undefined - presumably some input
> data. It might be worth posting a complete example, so that people can
> explore how to get decent performance.
>
> Error:
>
> Part::partd: "Part specification x[[1]] is longer than depth of object."
>
> David Baileyhttp://www.dbaileyconsultancy.co.uk

Hi David,

The RK4 procedure here works with the solution vector x whose initial
value is defined after the RK4 procedure and the equations are
defined. Mathematica does not know the lenght of x at the beginning
and this is why it complains. You can ignore these complaints. Of the
several codes posted above, that in Oliver's 19 September post is the
best because it is fully compiled. Here is this code with plotting the
solution:

***************************************
(* Runge-Kutta-4 routine *)
ClearAll[makeCompRK]
makeCompRK[f_] :=
 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](*,Parallelization->True*), CompilationTarget -> "C",
  CompilationOptions -> {"InlineCompiledFunctions" -> True}]

(* Defining equations *)
NN = 1000;
cRHS = With[{NN = NN}, Compile[{{t, _Real, 0}, {x, _Real, 1}},
    Table[-x[[i]]*
      Sin[0.1 t]^2/(1 +
         100 Sum[x[[i + j]], {j, 1, Min[3, NN - i]}]^2), {i, 1, NN}]
(*,
    CompilationTarget->"C"*)(*,
    CompilationOptions->{"InlineExternalDefinitions"->True}*)]];

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

(*Setting parameters and Calculation*)
x0 = Table[
  RandomReal[{0, 1}], {i, 1, NN}]; t0 = 0; tMax = 200; n = 500;
tt0 = AbsoluteTime[];
Sol = RK4Comp[x0, t0, tMax, n];
AbsoluteTime[] - tt0

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

(* Plotting *)
tList = Table[1. t0 + (tMax - t0) k/n, {k, 0, n}];
x1List = Transpose[{tList, Transpose[Sol][[1]]}];
x2List = Transpose[{tList, Transpose[Sol][[2]]}];
x3List = Transpose[{tList, Transpose[Sol][[3]]}];
ListLinePlot[{x1List, x2List, x3List}, PlotMarkers -> Automatic,
 PlotStyle -> {Blue, Green, Red}, PlotRange -> {0, 1}]

Best,

Dmitry

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