Mathematica 9 is now available
Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2006
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2006

[Date Index] [Thread Index] [Author Index]

Search the Archive

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


  • Prev by Date: Re: ndsolve & compile
  • Next by Date: Fastest method for comparing overlapping times in random time series
  • Previous by thread: Re: ndsolve & compile
  • Next by thread: Re: ndsolve & compile