Re: ndsolve & compile

*To*: mathgroup at smc.vnet.net*Subject*: [mg64941] Re: ndsolve & compile*From*: Peter Pein <petsie at dordos.net>*Date*: Wed, 8 Mar 2006 01:00:13 -0500 (EST)*References*: <duh3bs$5ut$1@smc.vnet.net> <dujs8u$95d$1@smc.vnet.net>*Sender*: owner-wri-mathgroup at wolfram.com

rknapp at wolfram.com schrieb: > A big part of the speedup I obtained from Alberto Verga's code was by > using Mathematica's listable arithmetic. > > However, there is another potential pitfall with precompiling a > function that you have > fallen into, so I'll show the advantages both of working around this, > and of using listable > arithmetic. > > I couldn't find a tmax in your messge for timing comparison, so I just > took the setting of > it to be tmax = 100. With your parameter settings, your compiled > function took 5.8 Seconds on my machine versus 0.312 seconds to handle > the final case you show > where you NDSolve process the equations (and automatically compile the > function) > > The problem with the CompiledFunction cf you have generated is that it > was compiled > without knowledge of the values of the parameters. To resolve this, it > has to call the > Mathematica evaluator each time a parameter is used to get the current > value -- without > this, the CompiledFunction would not be correct if you changed values > of the parameters. > > To work around this, you want to compile the function with the specific > values of the > parameters set. This is a little tricky, since Compile is HoldAll, so > parameters do not get evaluated. A convenient way to work around this > is to use the way Mathematica's > pattern matcher substitutes values to get specific values of the > parameters into the function. > > Here is a simple example. > > In[1]:= > a=1; > c1=Compile[{{n,_Integer}},Module[{h=N[Pi]/n},Table[Sin[a h > i],{i,0,n}]]] > > Out[2]= > CompiledFunction[{n}, > > N[Pi] > Module[{h = -----}, Table[Sin[a h i], {i, 0, n}]], > n > > -CompiledCode-] > > Note that a is represented as a symbol in the output of the > CompiledFunction. > Its value is checked each time through the loop, which is very > expensive > > In[3]:= > Timing[c1[10^6];] > > Out[3]= > {5.656 Second, Null} > > Here, I've jsut moved the evaluation of the parameter outside the loop > -- and the result is > the same but much faster. > > In[4]:= > a=1; > c2=Compile[{{n,_Integer}},Module[{h=a N[Pi]/n},Table[Sin[h > i],{i,0,n}]]] > > Out[5]= > CompiledFunction[{n}, > > N[Pi] > Module[{h = a -----}, Table[Sin[h i], {i, 0, n}]], > n > > -CompiledCode-] > > In[6]:= > Timing[c2[10^6];] > > Out[6]= > {0.235 Second, Null} > > Setting up something to make sure evaluation of parameters is outside > of the loop > can be cumbersome, and for your DEs would be inconvenient. > > When you define a function like f[x_] := x^2, Mathematica's pattern > matcher looks for > a match of f[x_] (f[whatever] will be a match), then substitutes the > values of named > patterns (x) in the right hand side *before* evaluating. This is > independent of evaluation > so the substitution is done whether or not there are HoldAll functions > in the right > hand side > > In[7]:= > CompileSin[a_] := Compile[{{n,_Integer}},Module[{h= > N[Pi]/n},Table[Sin[a h i],{i,0,n}]]]; > c3 = CompileSin[a] > > Out[8]= > CompiledFunction[{n}, > > N[Pi] > Module[{h$ = -----}, > n > > Table[Sin[1 h$ i], {i, 0, n}]], -CompiledCode-] > > When I call CompileSin[a] with a set to 1, the rhs becomes > > Compile[{{n,_Integer}},Module[{h= N[Pi]/n},Table[Sin[1 h i],{i,0,n}]] > > evaluating to the result you see above. This is fast since it does not > have to check the value of a. > > In[9]:= > Timing[c3[10^6];] > > Out[9]= > {0.234 Second, Null} > > When you use substitution to set parameters, you have to keep in mind > that > whenever you change parameter values, you have to recompile: > > In[10]:= > n = 2000;hn = 1001; > > In[11]:= > {c1[n][[hn]], c2[n][[hn]], c3[n][[hn]]} > > Out[11]= > {1., 1., 1.} > > In[12]:= > a = 2; > {c1[n][[hn]], c2[n][[hn]], c3[n][[hn]]} > > Out[13]= > -16 -16 > {1.22461 10 , 1.22461 10 , 1.} > > c3 was compiled for a specific value of a, so it does not pick up the > change in a. > > > Here is how you can use this method to compile your function. > > > In[11]:= > CompileRHS[r_, k_, kp_, c_, m_] := Compile[{{n,_Real,1},{p,_Real,1}}, > > Module[{dndt,dpdt,nx=Length[n]},Transpose[Table[dndt=r*n[[x]]*(1-n[[ > x]]/k)-c*n[[x]]*p[[x]]/(n[[x]]+kp); > dpdt=c*n[[x]]*p[[x]]/(n[[x]]+kp)-m*p[[x]]; > {dndt,dpdt},{x,1,nx}]]]]; > > In[12]:= > cf1 = CompileRHS[r, k, kp, c, m]; > g1[{n_?VectorQ,p_?VectorQ}]:=cf1[n,p]; > First[Timing[sol=First[NDSolve[{np'[t]\[Equal]g1[np[ > t]],np[0]\[Equal]{Table[0.69807,{x, > 1,nx}],Table[0.32324,{x,1,nx}]}},np,{t,0,tmax}]]]] > > Out[14]= > 0.312 Second > > which is about as fast as just giving NDSolve the DEs -- this is to be > expected since > NDSolve automatically compiles the equations in a form a bit more > efficient than using Part. > > In general, using Part inside a loop is not the fastest way to evaluate > expressions > in Mathematica. Mathematica's listable arithmetic is highly optimized > and is > usually faster. > > Here I compile the function to use listable arithmetic: > > In[15]:= > CompileRHSListable[r_, k_, kp_, c_, m_] := > Compile[{{n,_Real,1},{p,_Real,1}}, > { > r*n*(1 - n/k) - c*n*p/(n + kp), > c*n*p/(n + kp) - m*p > }]; > cfL = CompileRHSListable[r, k, kp, c, m]; > gL[{n_?VectorQ,p_?VectorQ}]:=cfL[n,p]; > First[Timing[sol= > First[NDSolve[{np'[t]\[Equal]gL[np[ > t]],np[0]\[Equal]{Table[0.69807,{x,1, > nx}],Table[0.32324,{x,1,nx}]}},np,{t,0,tmax}]]]] > > Out[18]= > 0.188 Second > > Getting an extra factor of 2 out of the situtation. > > Note that the advantages of using the compiler diminish when you are > using listable > arithmetic since the arithmetic operations are the same whether you > compile or not. > You can give the equations in a vector form > > In[19]:= > First[Timing[sol = NDSolve[ > {{n'[t] \[Equal] r*n[t]*(1 - n[t]/k) - c*n[t]*p[t]/(n[t] + kp), > p'[t] \[Equal] c*n[t]*p[t]/(n[t] + kp) - m*p[t]}, > {n[0] \[Equal] Table[0.69807,{x,1,nx}], > p[0] \[Equal] Table[0.32324,{x,1,nx}]}},{n,p}, > {t,0,tmax}]]] > > Out[19]= > 0.157 Second > > This is a little faster since the separation of n and p does not have > to be handled > by the pattern matcher. > > Note that in this example, the _?VectorQ is not needed since with > symbolic n[t] and p[t] > the rhs of the equations do not evaluate. > > > > Sorry for the long reply, but this issues are not trivial and require > some patience to understand. > > Rob Knapp > Wolfram Research > > Christopher Klausmeier wrote: [...] > This is really comprehensive. Is it possible to add this knowledge to the documentation ? Peter