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