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