Re: Compilation: Avoiding inlining

*To*: mathgroup at smc.vnet.net*Subject*: [mg121459] Re: Compilation: Avoiding inlining*From*: Oliver Ruebenkoenig <ruebenko at wolfram.com>*Date*: Fri, 16 Sep 2011 05:46:21 -0400 (EDT)*Delivered-to*: l-mathgroup@mail-archive0.wolfram.com*References*: <201109150841.EAA29262@smc.vnet.net>

On Thu, 15 Sep 2011, DmitryG wrote: > This message refers to the extensive exchange here: > > http://groups.google.com/group/comp.soft-sys.math.mathematica/browse_thread/thread/4ea9938ad1d523c6/735f864eab188778?lnk=gst&q=dmitryg#735f864eab188778 > > The problem is to solve a system of N ordinary differential equations > (ODE) with Runge-Kutta-4 method with compilation. Here is an example > of working program > > *************************************************************************** > (* Runge-Kutta-4 routine *) > 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"*)]]; > > (* Definition of the equations *) > NN = 100; > Timing[ > RHS = Quiet[ > Table[-x[[i]] Sin[ > 0.3 t]^2/(1 + 300 Sum[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[] > > (* Compilation *) > tt0 = AbsoluteTime[]; > Timing[RK4Comp = makeCompRK[F];] > AbsoluteTime[] - tt0 > Print[]; > > (* Setting parameters and Calculation *) > x0 = Table[ > RandomReal[{0, 1}], {i, 1, NN}]; t0 = 0; tMax = 100; 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}] > > ********************************************************************************** > > The drawback of this solution is what Oliver Ruebenkoenig called > "inlining". For large NN the RHS of the equations is a lot the > explicit code generated by Mathematica. This code then gets compiled > and the size of the compiled code strongly increases with NN (it can > be seen by the CompilePrint[RK4Comp] command). That is, we write a lot > of code inline instead of letting it be generated in the course of the > execution, with the compiled code being small. > > Mathematica's own compiler tolerates inlining, although the large code > eats up the memory. If we uncomment the option "CompilationTarget- >> "C"", we will see that for large NN compilation goes for ages. That > is, inlining has to be avoided. > > Daniel Lichtblau has rewritten my similar program in the vectorized > form that avoided inlining and could be compiled to C, and was fast. > But vectorization is not always possible or easy. Oliver wrote that > inlining can also be avoided by defining a function and then calling > it, but there were no examples at that time. Now I am interested in > the cases such as the program above where vectorization is unwanted. > How to avoid inlining by defining a function in this case? I have > tried to replace > > RHS = Quiet[Table[-x[[i]] Sin[0.3 t]^2/(1 + 300 Sum[x[[j]], {j, 1, > NN}]^2/NN^2), {i, 1, NN}]]; > > by > > FRHS[i_][x_List, t_] = -x[[i]] Sin[0.3 t]^2/(1 + 100 Sum[x[[i + j]], > {j, 1, Min[3, NN - i]}]^2); > RHS = Quiet[Table[FRHS[i][x, t], {i, 1, NN}]]; > > This works correctly but does not compile, and the execution time is > longer. Something is wrong here. Oliver, Daniel, and others, please, > help! > > Thank you, > > Dmitry > > Dmitry, here is an idea that hopefully will get you started: ClearAll[makeCompRK] 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",*) CompilationOptions -> {"InlineCompiledFunctions" -> True}]]; (* Have a look at InlineCompiledFunctions *) (*Definition of the equations*) NN = 100; Timing[RHS = Quiet[Table[-x[[i]] Sin[0.3 t]^2/(1 + 300 Sum[x[[j]], {j, 1, NN}]^2/NN^2), {i, 1, NN}]];] (*Quiet to suppress complaints*) F = Function[{t, x}, Evaluate[RHS]]; >> 0.06 (*Compilation*) tt0 = AbsoluteTime[]; Timing[RK4Comp = makeCompRK[F];] AbsoluteTime[] - tt0 >> 0.2 cRHS = Compile[{{t, _Real, 0}, {x, _Real, 1}}, -x*(Sin[0.3*t]^2/(1 + 300*Total[x]^2/Length[x]^2))] (*cRHS[0.3,{1,2,3,4}//N]*) (*Compilation*) tt0 = AbsoluteTime[]; Timing[RK4Comp2 = makeCompRK[cRHS];] AbsoluteTime[] - tt0 >> 0.002 (*Needs["CompiledFunctionTools`"]*) (*CompilePrint[RK4Comp2]*) (* Have a look at how "InlineCompiledFunctions" -> True/False influence the output *) (*Setting parameters and Calculation*) x0 = Table[RandomReal[{0, 1}], {i, 1, NN}]; t0 = 0; tMax = 100; n = 500; tt0 = AbsoluteTime[]; Sol1 = RK4Comp[x0, t0, tMax, n]; AbsoluteTime[] - tt0 >> 0.43 tt0 = AbsoluteTime[]; Sol2 = RK4Comp2[x0, t0, tMax, n]; AbsoluteTime[] - tt0 >> 0.007 Norm[Sol1 - Sol2] >> 10^-14 Sol = Sol2; Print["Compilation: ", Developer`PackedArrayQ@Sol] Some further notes: - I did not try the ->"C" version, but I believe both compiles need to be compiled to C to be effective. Have a play. I'd be interested to about the results - I understand that you wanted an non vectorized version. That in principal would work in the same way; you could provide an example of such a non-vectorizable function. - I think this is a good illustration of why compilation to C is not always beneficial. It is only beneficial if the compile time vs the runtime is small. Hope this gets you started. Oliver

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