1
$\begingroup$

I've got some problems with my code, and I try to make it faster. Some of you suggested me to split my problem, and I'm here...

I post the same function with and without Compile.

This is a support function and does not involve my question

r=8.314472
fug1 = Compile[{v, p, t, a, b},
Module[{y, z, vbv, vb, f1, f2, f3, f4, f},
y = b/(4 v);
z = (p v)/(r t);
vbv = Log[(v + b)/v];
vb = v + b;
f1 = (4.*y - 3.*y^2.)/(1 - y)^2.;
f2 = (4.*y - 2.*y^2.)/(1 - y)^3.;
f3 = (2.*vbv)/(r t*b)*a;
f4 = (vbv/b - 1./vb)/(r t)*a;
f = f1 + f2 - f3 + f4 - Log[z];
Exp[f]]]

g2 doesen't work and I don't know why

g2 = Compile[{p, t, a0, a1, a2, b0, b1, b2},
Module[{a, b, csd, vol, sol, vliquid, vvapor, fl, fv},
a = a0*Exp[a1*t + a2*t^2];
b = b0 + b1*t + b2*t^2;
csd = a/(r*t*(b + v)) - (-(b^3/(64.*v^3)) + b^2/(16.*v^2.0) + 
    b/(4.*v) + 1.)/(1 - b/(4*v))^3 + (p*v)/(r*t);
vol = NSolve[csd == 0. && v > 0., v, Reals] // Quiet;
sol = v /. vol;
vliquid = Min[sol];
vvapor = Max[sol];
fl = fug1[vliquid, p, t, a, b];
fv = fug1[vvapor, p, t, a, b];
Print[{t, p, vol, Abs[fl - fv]}];
Abs[fl - fv]],
RuntimeAttributes -> {Listable}]

This works without Compile!

g[p_?NumericQ, t_?NumericQ, a0_?NumericQ, a1_?NumericQ, a2_?NumericQ, 
  b0_?NumericQ, b1_?NumericQ, b2_?NumericQ] := 
Module[{a, b, csd, vol, sol, vliquid, vvapor, fl, fv},
a = a0*Exp[a1*t + a2*t^2];
b = b0 + b1*t + b2*t^2;
csd = a/(r*t*(b + v)) - (-(b^3/(64.*v^3)) + b^2/(16.*v^2.0) + 
    b/(4.*v) + 1.)/(1 - b/(4*v))^3 + (p*v)/(r*t);
vol = NSolve[csd == 0. && v > 0., v, Reals];
sol = v /. vol;
vliquid = Min[sol];
vvapor = Max[sol];
fl = fug1[vliquid, p, t, a, b];
fv = fug1[vvapor, p, t, a, b];
Print[{t, p, vol, Abs[fl - fv]}];
Abs[fl - fv]];

g works very well and is the same of g2!

FindRoot[g[p, 100, 500., -4.4627562855*10^-3, -2.7625748*10^-6, 
7.30402014*10^-2, -2.2222592*10^-4, 9.42486*10^-8], {p, 
34.376}] // Timing

g2 doesen't work

FindRoot[g2[p, 100, 500., -4.4627562855*10^-3, -2.7625748*10^-6, 
7.30402014*10^-2, -2.2222592*10^-4, 9.42486*10^-8], {p, 
34.376}] // Timing
Improve this question
$\endgroup$
2
  • 5
    $\begingroup$ There are many issues with with your functions. Variables are not localized, r is not defined in fug1, NSolve is not compilable and you probably mean to use g2 but are using function g1 in the FindRoot. $\endgroup$
    – asim
    Commented Jan 27, 2013 at 16:00
  • $\begingroup$ You're right. I've corrected it, I've localized my variables and I've inserted all your corrections, but it doesn't work again $\endgroup$
    – meriens
    Commented Jan 27, 2013 at 18:16

1 Answer 1

Reset to default
3
$\begingroup$

Let's start with solving your equation with Solve:

rv = r -> 8.314472;
csd = a/(r*t*(b + v)) - (-(b^3/(64*v^3)) + b^2/(16*v^2) + b/(4*v) + 
     1)/(1 - b/(4*v))^3 + (p*v)/(r*t);
rootsols = Solve[csd == 0, v] /. rv;
v /. Last[rootsols]

Mathematica graphics

I change fug1 slightly:

fug1 = Compile[{{v, _Real, 1}, p, t, a, b}, 
   Module[{y, z, vbv, vb, f1, f2, f3, f4, f, r},
    r = 8.314472;
    y = b/(4 v);
    z = (p v)/(r t);
    vbv = Log[(v + b)/v];
    vb = v + b;
    f1 = (4.*y - 3.*y^2.)/(1 - y)^2.;
    f2 = (4.*y - 2.*y^2.)/(1 - y)^3.;
    f3 = (2.*vbv)/(r t*b)*a;
    f4 = (vbv/b - 1./vb)/(r t)*a;
    f = f1 + f2 - f3 + f4 - Log[z];
    Exp[f]]];

Now, g2, using Root instead of NSolve:

    g2 = With[{fug1 = fug1}, 
   Compile[{p, t, a0, a1, a2, b0, b1, b2}, 
    Module[{a, b, csd, vol, vliquid, vvapor, fl, fv}, 
     a = a0*Exp[a1*t + a2*t^2];
     b = b0 + b1*t + b2*t^2;
     sol = {};
     i = 0;
     While[i < 5, i++;
      soli =
       Root[-a b^3 + 
          8.314472` b^4 t + (12 a b^2 - b^4 p - 
             24.943416` b^3 t) #1 + (-48 a b + 11 b^3 p - 
             166.28944` b^2 t) #1^2 + (64 a - 36 b^2 p - 
             665.15776` b t) #1^3 + (16 b p - 532.126208` t) #1^4 + 
          64 p #1^5 &, i];
      If[Im[soli] == 0, sol = Append[sol, soli]]
      ];

     vliquid = Min[sol];
     vvapor = Max[sol];
     {fl, fv} = fug1[{vliquid, vvapor}, p, t, a, b];

     Abs[fl - fv]], {{soli, _Complex}, {sol, _Real, 1}}, 
    CompilationOptions -> {"InlineCompiledFunctions" -> True},
    RuntimeAttributes -> {Listable}
    ]];

Testing with g:

g[p_?NumericQ, t_?NumericQ, a0_?NumericQ, a1_?NumericQ, a2_?NumericQ, 
   b0_?NumericQ, b1_?NumericQ, b2_?NumericQ] := 
  Module[{a, b, csd, vol, sol, vliquid, vvapor, fl, fv}, 
   a = a0*Exp[a1*t + a2*t^2];
   b = b0 + b1*t + b2*t^2;
   csd = a/(r*t*(b + v)) - (-(b^3/(64.*v^3)) + b^2/(16.*v^2.0) + 
        b/(4.*v) + 1.)/(1 - b/(4*v))^3 + (p*v)/(r*t);
   vol = NSolve[csd == 0. && v > 0., v, Reals];
   sol = v /. vol;
   vliquid = Min[sol];
   vvapor = Max[sol];
   {fl, fv} = fug1[{vliquid, vvapor}, p, t, a, b];
   (*Print[{t,p,vol,Abs[fl-fv]}];*)
   Abs[fl - fv]];
r = 8.314472; 
FindRoot[g[p, 100, 500., -4.4627562855*10^-3, -2.7625748*10^-6, 
    7.30402014*10^-2, -2.2222592*10^-4, 9.42486*10^-8], {p, 34.376}] //
   Quiet // Timing

gives

(* 
       {0.218401, {p -> 170.441}}
*)

and with the compile function (need a simple extra function g3 here):

g3[args__?NumberQ] := g2[args];
FindRoot[g3[p, 100, 500., -4.4627562855*10^-3, -2.7625748*10^-6, 
    7.30402014*10^-2, -2.2222592*10^-4, 9.42486*10^-8], {p, 34.376}] //
   Quiet // Timing

gives

(*
    {0.093601, {p -> 170.441}}
*)

So a speed-up of a factor of 2. Why not more? Because even though g2 compiles there is still a MainEvaluate (call-back to main Mathematica) to evaluate the Root object. Maybe things will get even faster if you use a C-compiler.

Improve this answer
$\endgroup$
2
  • $\begingroup$ I think the performance improvement here is due to your choice to express the problem slightly differently rather than because of Compile. Calls out of the VM are very expensive and there are just too many of them here to make compilation of g2 worthwhile. (Although you forgot to localize i, it doesn't seem to make much difference overall.) Replacing Compile with Function actually improves the timings for me. $\endgroup$
    – Oleksandr R.
    Commented Jan 28, 2013 at 8:44
  • $\begingroup$ Thank's Oleksandr. Basically I just wanted to show how to get Compile to work at all. Did you test CompilationTarget -> "C" ? I have no time to investigate this further. Feel free to edit and extend my answer. $\endgroup$
    – Rolf Mertig
    Commented Jan 28, 2013 at 11:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.