# I can't get Compile to work

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

-
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. – asim Jan 27 '13 at 16:00
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 – meriens Jan 27 '13 at 18:16

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]


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.

-
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. – Oleksandr R. Jan 28 '13 at 8:44
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. – Rolf Mertig Jan 28 '13 at 11:55