Re: Vector Runge-Kutta ODE solver with compilation?
- To: mathgroup at smc.vnet.net
- Subject: [mg116805] Re: Vector Runge-Kutta ODE solver with compilation?
- From: Oliver Ruebenkoenig <ruebenko at wolfram.com>
- Date: Tue, 1 Mar 2011 05:20:29 -0500 (EST)
On Mon, 28 Feb 2011, DmitryG wrote:
> On 27 Feb., 04:43, Oliver Ruebenkoenig <ruebe... at wolfram.com> wrote:
>
>> You could use CompilePrint to see that for NN=11 this generates a some 900
>> line so code - quite a lot of which are Part[]. Make sure you understand
>> what you are doing here - you in-line this sum into the code.
>>
>>
> ......................
>>
>> Weird code? - No, not really.
>> ...
>>
Hi Dmitry,
>
> Hi Oliver,
>
> In this short program I mimick my real-life program that solves a
Why do you not send the example that you classify as real-life?
> system of NN ODE's with the maximal coupling: the RHS of each of NN
> equations contains all NN variables. Any molecular dynamics (MD)
> simulation is of this kind: Acceleration of any particle depends on
> positions of all particles.
>
One method I have seen for MD programs is to use a cut-off range that
mimics for example the Lennard-Jones potential - which reduces the
complexity of the alg.
> This application of Mathematica's compilation is crucial but there are
> two problems that I have encountered:
>
> 1) Compilation in C is prohibitively slow above NN=50 or so. In this
> case Mathematica's own compiler is still doing a good job.
>
I think the in-lining of the sum is not what you want to do here. You
should either use vectorized code or another compiled function, as shown
in a previous post.
> 2) For NN=1000 or above Mathematica's compiler consumes all RAM (my
> laptop has 3 GB RAM) and quits. This is very disappointing because
> generating the list of RHS does not take much memory and the RK 4
> routine does not consume much memory as well. Thus it would be only a
> matter of computer time to solve a problem for any NN. However,
> compiler generates a blown-up code (component by component, with the
> vector structure abandoned) that does not fit into memory.
>
No it is not that compiler that generated this - you told it to in-line that
into the code and that's what happened. - A compiler is optimized to
optimize hand written code - it can not do logical simplifications.
Would you copy and past your sum function in a C file? No, you'd write a
function and call it. In Mathematica you additionally have the option to
use vectorized code. What is wrong with a vector approach?
> Thank you for suggestion to look into the compiled code with
> CompilePrint. Yes, the code is very long, and its structure is the
> same for the initial problem we considered, the problem with
>
> RHS=Quiet[{x[[2]],1-x[[1]]+4/x[[1]]^3}]
>
Just the slight difference that in the initial problem NN=2.
> Also in this case the compiled code contains Part[].
>
> Is there any way to optimize this calculation or we have an intrinsic
> limitation of Mathematica here? Below is the code again:
>
The limitation is in the code you inline. It is not that you inline, but
what you inline.
> *******************************************************************************
> (* Definition *)
> 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"*)]];
>
> (* Calculation *)
> NN=15;
> Timing[RHS=Quiet[Table[-x[[i]]/(1+300Sum[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[]
>
> tt0=AbsoluteTime[];
> Timing[RK4Comp=makeCompRK[F];]
> AbsoluteTime[]-tt0
> Print[];
>
> x0=Table[RandomReal[{0,1}],{i,1,NN}]; t0=0; tMax=50; 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}]
>
>
> *********************************
>
> Thank you for your help,
>
> Dmitry
>
>
Here is the code for NN=1000, tMax0 and n = 50000
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"]];
NN = 1000;
RHS = Quiet[-x/(1 + 300*Total[x]^2/NN^2)]
F = Function[{t, x}, Evaluate[RHS]]
tt0 = AbsoluteTime[];
Timing[RK4Comp = makeCompRK[F];]
AbsoluteTime[] - tt0
x0 = RandomReal[{0, 1}, {NN}]; t0 = 0; tMax = 200; n = 50000;
tt0 = AbsoluteTime[];
Sol = RK4Comp[x0, t0, tMax, n];
AbsoluteTime[] - tt0
runs in about 3.7 sec on my laptop.
Oliver