Re: Vector Runge-Kutta ODE solver with compilation?
- To: mathgroup at smc.vnet.net
- Subject: [mg116782] Re: Vector Runge-Kutta ODE solver with compilation?
- From: Oliver Ruebenkoenig <ruebenko at wolfram.com>
- Date: Sun, 27 Feb 2011 04:36:08 -0500 (EST)
On Sat, 26 Feb 2011, DmitryG wrote:
> On 25 Feb., 06:38, Oliver Ruebenkoenig <ruebe... at wolfram.com> wrote:
>> On Thu, 24 Feb 2011, DmitryG wrote:
>>> On 23 Feb., 06:25, Oliver Ruebenkoenig <ruebe... at wolfram.com> wrote:
>>>> On Wed, 23 Feb 2011, DmitryG wrote:
>>
>>>> [....]
>>
>>>>> Dear Oliver and All interested,
>>
>>>>> Today I tried to apply the above very efficient vector RK4 method with
>>>>> compilation to my real problems that I was solving with NDSolve
>>>>> before. Unfortunately, it has proven to be more complicated than it
>>>>> seemed. To illustrate the difficulties, I will use the same working
>>>>> model as above and will try to generalize it. So, the working code is
>>
>>>>> *************************************************
>>
>>>>> (* Definition of the RK4 solver *)
>>>>> 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 *)
>>>>> F = Function[{t, x}, {x[[2]], 1 - x[[1]] + 4/x[[1]]^3}];
>>>>> RK4Comp = makeCompRK[F];
>>>>> x0 = {1.0, 1.0}; t0 = 0; tMax = 200; n = 5000;
>>>>> tt0 = AbsoluteTime[];
>>>>> Sol = RK4Comp[x0, t0, tMax, n];
>>>>> AbsoluteTime[] - tt0
>>>>> Developer`PackedArrayQ@Sol
>>
>>>>> (* Plotting *)
>>>>> nTake = 200;
>>>>> ListLinePlot[Take[Sol, nTake], PlotMarkers -> Automatic, PlotStyle ->
>>>>> {Red}] (* Trajectory *)
>>
>>>>> (* 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=
> }]=
>>> ;
>>>>> ListLinePlot[{Take[x1List, nTake], Take[x2List, nTake]},
>>>>> PlotMarkers -> Automatic, PlotStyle -> {Blue, Green}]
>>>>> *********************************************************************=
> ****** > > **************************
>>
>>>>> In this example, the vector of ODE's RHS' is defined by
>>
>>>>> F = Function[{t, x}, {x[[2]], 1 - x[[1]] + 4/x[[1]]^3}];
>>
>>>>> In real life, this is something more complicated, a result of a many-
>>>>> step calculation. One can try to generalize the code as
>>
>>>>> RHS={x[[2]], 1 - x[[1]] + 4/x[[1]]^3};
>>>>> F = Function[{t, x}, RHS];
>>
>>>> Dmitry,
>>
>>>> your friend is - in princial -
>>
>>>> F = Function[{t, x}, Evaluate[RHS]]
>>
>>>> That has a small problem that we used Part in the RHS and Part gives a
>>>> message that you try to access elements that the _symbol_ x does not h=
> ave.
>>
>>>> One way to work around this is
>>
>>>> RHS = Quiet[{x[[2]], 1 - x[[1]] + 4/x[[1]]^3}];
>>>> F = Function[{t, x}, Evaluate[RHS]];
>>
>>>> If you evaluate F you'll see that the RHS in now in the function body.=
> If
>>>> you look at the F you posted, there the body is RHS, which is not what
>>>> you want. All this hast to do with
>>
>>>> Attributes[Function]
>>
>>>> HoldAll.
>>
>>>> Oliver
>>
>>>>> (that should be equivalent to the above) and then substitute RHS by
>>>>> whatever you have for your problem. However, Mathematica complains an=
> d
>>>>> produces crap results. My suspicion is that it is only F =
>>>>> Function[{t, x}, RHS]; that is injected into the compiled RK4 code
>>>>> whereas RHS={x[[2]], 1 - x[[1]] + 4/x[[1]]^3}; is not injected.
>>>>> However, I cannot be sure because it was never explained.
>>
>>>>> Further idea is to define F as a procedure and put everything inside
>>>>> it. A toy example is
>>
>>>>> F = Function[{t, x}, RHS = {x[[2]], 1 - x[[1]] + 4/x[[1]]^3}; RHS=
> ];
>>
>>>>> This does not compile, although then retreating to the non-compilatio=
> n
>>>>> mode produces correct results. My experiments show that you can only
>>>>> use a single expression in the body of F that does not use any
>>>>> definitions. For instance,
>>
>>>>> F = Function[{t, x}, Table[-x[[i]]^i, {i, 1, 2}]];
>>
>>>>> compiles and works fine but this is, again, no more than a toy
>>>>> example.
>>
>>>>> To summarize, I am very excited by the prospect to speed up my
>>>>> computations by a factor about 100 but I don't see how to generalize
>>>>> the approach for real-life problems.
>>
>>> Thank you, Oliver, this was a great help.
>>
>>> Now, I am trying to introduce sampling of the results to reduce the
>>> number of output points. This is also very important. However, there
>>> is a new problem. Here is the code:
>>
>>> ***********************************************************************=
> **** ***
>>> (* Vector RK4 routine *)
>>> makeCompRK[fIn_]:=With[{f=fIn},Compile[{{x0,_Real,1},{t0,_Real},
>>> {ttMax,_Real},{n,_Integer},{nsamp,_Integer}},
>>> Module[{h,K1,K2,K3,K4,SolList,x=x0,t,ksamp},
>>> h=(ttMax-t0)/n;
>>> SolList=Table[x0,{nsamp+1}];
>>> Do[t=t0+k h; ksamp=k nsamp/n;
>>> 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;
>>
>>> If[FractionalPart[ksamp]==0,SolList[[ksamp+1]]=x]
>>
>>> ,{k,1,n}];
>>> SolList],CompilationTarget->"C"]];
>>
>>> (* Calculation *)
>>> RHS=Quiet[{x[[2]],1-x[[1]]+4/x[[1]]^3}]; (* Quiet to suppress
>>> complaints *)
>>> F=Function[{t,x},Evaluate[RHS]]; (* Evaluate to use actual form of
>>> RHS *)
>>> RK4Comp=makeCompRK[F];
>>> x0={1.0,1.0}; t0=0; tMax; NSteps=5000; NSamplin=
> g0;
>>> tt0=AbsoluteTime[];
>>> Sol=RK4Comp[x0,t0,tMax,NSteps,NSampling];
>>> AbsoluteTime[]-tt0
>>> Developer`PackedArrayQ@Sol
>>> ***********************************************************************=
> **** ***
>>
>>> The problem is in the line
>>
>>> If[FractionalPart[ksamp]==0,SolList[[ksamp+1]]=x]
>>
>>> In the first execution, Mathematica complains
>>
>>> Compile::cpintlt: "ksamp$+1 at position 2 of SolList$[[ksamp$+1]]
>>> should be either a nonzero integer or a vector of nonzero integers;
>>> evaluation will use the uncompiled function."
>>
>> Of what type is ksamp?
>>
>> ksamp = k nsamp/n
>>
>> For example
>>
>> Head[5 3/2]
>>
>> In worst case this is a rational and you check that with
>> FractionalPart. But Part needs an integer at position 2. With using Floor
>> you can tell the compiler that ksamp is always going to be an integer.
>>
>> SolList[[Floor[ksamp] + 1]]
>>
>> Oliver
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>> I do not understand this. Still Mathematica compiles and the results
>>> are correct.
>>
>>> In the second execution, Mathematica complains more and does not
>>> compile, although the results are still correct.
>>
>>> If I use the command
>>
>>> If[IntegerQ[kSampling],SolList[[kSampling+1]]=x
>>
>>> Mathematica complains again:
>>
>>> CompiledFunction::cfse: "Compiled expression False should be a machine-
>>> size real number)."
>>
>>> does not compile, but the results are correct.
>>
>>> What is wrong here? I believe I am doing everything in a natural way.
>>> But compiling is too tricky..
>>
>>> Dmitry
>
> Thank you for the response, Oliver! I have already fugured it out.
>
> There is now a problem with compilation in C in the case of a system
> of N equations. For N above 10, it takes a long time to compile, and
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.
> the compilation time is dramatically increasing with N, so that one
> cannot solve any real-life problems with large N using compilation in
> C. On the other hand, using Mathematca's own compiler is non-
> problematic. In both cases the solution is correct. Here is the code:
>
> ***************************************************************************=
> ***************
> (* 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=30;
> RHS=Quiet[Table[-x[[i]]/(1+Sum[x[[j]],{j,1,NN}]^2),{i,1,NN}]]; (*
> Quiet to suppress complaints *)
> F=Function[{t,x},Evaluate[RHS]]; (* Evaluate to use actual form of
> RHS *)
>
> tt0=AbsoluteTime[];
> Timing[RK4Comp=makeCompRK[F];]
> AbsoluteTime[]-tt0
> Print[];
>
> x0=Table[RandomReal[{0,1}],{i,1,NN}]; t0=0; tMax0; 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}]
>
>
> Compile::nogen: A library could not be generated from the compiled
> function. >>
> Out[4]= {11.172,Null}
> Out[5]= 142.0625000
>
>
> Out[10]= 0.0468750
> Compilation: True
>
> ***************************************************************************=
> **********************************
>
> Whereas execution is fast, 0.0468750, compilation takes a long time. I
> measure both Timing and AbsoluteTime. The former, I believe, is the
> time used by Mathematica, whereas the latter is the total time used by
> Mathematica and my C compiler, Microsoft Visual C++ 2010 Express. Both
> times are long. Presumably, Mathematica generates some weird code that
Weird code? - No, not really.
> gives a hard time to C compiler. Here I also have an error message
> concerning a library but I haven' had this message as I ran this
> program before.
>
> Any help?
>
> Dmitry
>
>
You have several options
1)
F = With[{NN = NN}, Compile[{{t, _Real, 0}, {x, _Real, 1}},
Table[-x[[i]]/(1 + Sum[x[[j]], {j, 1, NN}]^2), {i, 1, NN}],
CompilationTarget -> C]];
Use a compiled function and insert that into the RK4 code. This is like
another function call.
- If you look at your code - that is does a lot of computations
repetitively, the sum is computed for every i - I do not think that any
compiler can optimize that.
2) Use vector code
F = Compile[{{t, _Real, 0}, {x, _Real, 1}}, -x/(1 + Total[x]^2),
CompilationTarget -> C];
3) To save the additional function call you can the inline this piece
code.
RHS = Quiet[(-x/(1 + Total[x]^2))]
F = Function[{t, x}, Evaluate[RHS]]
NN = 30;
x0 = Table[
RandomReal[{0, 1}], {i, 1, NN}]; t0 = 0; tMax = 2000; n = 50000;
(* n = 50000 *)
tt0 = AbsoluteTime[];
Timing[RK4Comp = makeCompRK[F];]
AbsoluteTime[] - tt0
{0.104007, Null}
0.257110
tt0 = AbsoluteTime[];
Sol = RK4Comp[x0, t0, tMax, n];
AbsoluteTime[] - tt0
1.481468
Oliver