Re: Vector Runge-Kutta ODE solver with compilation?
- To: mathgroup at smc.vnet.net
- Subject: [mg116777] Re: Vector Runge-Kutta ODE solver with compilation?
- From: DmitryG <einschlag at gmail.com>
- Date: Sun, 27 Feb 2011 04:35:14 -0500 (EST)
- References: <ik84bo$n4e$1@smc.vnet.net> <ikamss$c8o$1@smc.vnet.net>
On 26 Feb., 06:07, DmitryG <einsch... at gmail.com> 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 w=
ith
> > >>> 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 = 500=
0;
> > >>> 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 man=
y-
> > >>> 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 bod=
y.=
> If
> > >> you look at the F you posted, there the body is RHS, which is not wh=
at
> > >> 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 insid=
e
> > >>> it. A toy example is
>
> > >>> F = Function[{t, x}, RHS = {x[[2]], 1 - x[[1]] + 4/x[[1]]^3}; R=
HS=
> ];
>
> > >>> This does not compile, although then retreating to the non-compilat=
io=
> n
> > >>> mode produces correct results. My experiments show that you can onl=
y
> > >>> 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 generaliz=
e
> > >>> 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; NSampl=
in=
> 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 Flo=
or
> > 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 machin=
e-
> > > 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
> 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
> 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
Another major problem: If in the code above one removes
CompilationTarget->"C", it is working well until NN becomes larger
then ~1000. For larger NN, say NN00, compilation breaks because of
lack of memory on my laptop. One can see how fast the memory is
getting consumed during compilation. This seems to be a deficiency of
the current version of the Mathematica compiler that probably makes it
inappropriate for solving larger problems where compilation is really
needed. And there seems to be an even bigger problem with using
external C compiler.
Dmitry