4
$\begingroup$

I need to speed up a large program that repeatedly calls the module below (I've not trimmed the code for this question to ensure there is enough calculation to show a speed improvement.) The obvious way of making a step change in speed seemed to be use of Compile, so I create the second version show below, however, there is very marginal speed improvement. Am I doing something wrong?

ZrootsAnalytic[A_Real, B_Real] :=
 Module[
  { Zroot1, Zroot2, Zroot3, term1, term2},
  term1 = (-1. + 3. A - 4. B - 10. B^2.);
  term2 = (2. - 9. A + 12. B + 36. A B - 12. B^2. - 
     56. B^3. + \[Sqrt](4. term1^3. + (2. - 9. A + 12. B + 36. A B - 
          12. B^2. - 56. B^3.)^2.))^(1./3.);

  Zroot1 = (1. - B)/3. - (2.^(1./3.) term1)/(3. term2) + 
    1./(3. 2.^(1./3.)) term2;

  Zroot2 = (1. - B)/
    3. + ((1. + I Sqrt[3.]) term1)/(3. 2.^(2./3.) term2) - 
    1./(6. 2.^(1./3.)) (1. - I Sqrt[3.]) term2;

  Zroot3 = (1. - B)/
    3. + ((1. - I Sqrt[3.]) term1)/(3. 2.^(2./3.) term2) - 
    1./(6. 2.^(1./3.)) (1. + I Sqrt[3.]) term2;
  {Zroot1, Zroot2, Zroot3 }
  ]

Compiled version of the code.

ZrootsAnalyticCC = Compile[
  {{A, _Real}, {B, _Real}}, 

  term1 = (-1. + 3. A - 4. B - 10. B^2.);
  term2 = (2. - 9. A + 12. B + 36. A B - 12. B^2. - 
     56. B^3. + \[Sqrt](4. term1^3. + (2. - 9. A + 12. B + 36. A B - 
          12. B^2. - 56. B^3.)^2.))^(1./3.);

  Zroot1 = (1. - B)/3. - (2.^(1./3.) term1)/(3. term2) + 
    1./(3. 2.^(1./3.)) term2;

  Zroot2 = (1. - B)/
    3. + ((1. + I Sqrt[3.]) term1)/(3. 2.^(2./3.) term2) - 
    1./(6. 2.^(1./3.)) (1. - I Sqrt[3.]) term2;

  Zroot3 = (1. - B)/
    3. + ((1. - I Sqrt[3.]) term1)/(3. 2.^(2./3.) term2) - 
    1./(6. 2.^(1./3.)) (1. + I Sqrt[3.]) term2;
  {Zroot1, Zroot2, Zroot3 }
  ,
  {{ Zroot1, _Complex, 0}, {Zroot2, _Complex, 0}, {Zroot3, _Complex, 
    0}, {term1, _Complex, 0}, { term2, _Complex, 0}}, 
  CompilationTarget -> "C",
  RuntimeAttributes -> {Listable},
  RuntimeOptions -> "Speed",
  CompilationOptions -> {"InlineCompiledFunctions" -> True}
  ]

Is it the use of imaginary no.'s that prevents speedup?

================== Here is a modified version of the code as suggested by Mark

ZrootsAnalyticCC =
 Compile[
  {{A, _Real}, {B, _Real}}, 
  Module[{term1, term2, Zroot1, Zroot2, Zroot3},
   term1 = (-1. + 3. A - 4. B - 10. B^2.);
   term2 = (2. - 9. A + 12. B + 36. A B - 12. B^2. - 
      56. B^3. + \[Sqrt](4. term1^3. + (2. - 9. A + 12. B + 36. A B - 
           12. B^2. - 56. B^3.)^2.))^(1./3.);

   Zroot1 = (1. - B)/3. - (2.^(1./3.) term1)/(3. term2) + 
     1./(3. 2.^(1./3.)) term2;

   Zroot2 = (1. - B)/
     3. + ((1. + I Sqrt[3.]) term1)/(3. 2.^(2./3.) term2) - 
     1./(6. 2.^(1./3.)) (1. - I Sqrt[3.]) term2;

   Zroot3 = (1. - B)/
     3. + ((1. - I Sqrt[3.]) term1)/(3. 2.^(2./3.) term2) - 
     1./(6. 2.^(1./3.)) (1. + I Sqrt[3.]) term2;
   {Zroot1, Zroot2, Zroot3 }
   ]
  ,
  {{ Zroot1, _Complex, 0}, {Zroot2, _Complex, 0}, {Zroot3, _Complex, 
    0}, { term1, _Complex, 0}, { term2, _Complex, 0}}, 
  CompilationTarget -> "C"
  ]

I get no speed up and an error message:

  ZrootsAnalyticCC[1., 1.] // AbsoluteTiming

    CompiledFunction::cfne: Numerical error encountered; proceeding with uncompiled evaluation. >>

{0.046061, {1.86081 - 2.22045*10^-16 I, -2.11491 - 4.44089*10^-16 I, 
  0.254102 + 4.44089*10^-16 I}}

However, the problem goes away when I specify {{A, _Complex}, {B, _Complex}}, the AbsoluteTiming for a single evaluation of ZrootsAnalyticCC is still about the same as the uncompiled version, but there is a 25x speed up when I evaluate it many times in table:

Table[ ZrootsAnalyticCC[1., 1.], {100000}]; // AbsoluteTiming
$\endgroup$
8
$\begingroup$

You might think that specifying the types of intermediate expressions as the third argument of Compile[] would make them local, but it does not. You need to put the expression in a Module[] with those variables in the first argument of Module[] as you did in the non-compiled code. Without that, the compiled code is effectively not compiled, since it keeps going back to Mathematica to set the global variables Zroot1, etc.

Your clue could have been that after running the compiled code, those variable names turn from blue to black. (Assuming that your front end works like mine.)

Update:

This is what I did and saw:

ZrootsAnalytic[A_Real, B_Real] := 
 Module[{Zroot1, Zroot2, Zroot3, term1, term2}, 
  term1 = (-1. + 3. A - 4. B - 10. B^2.);
  term2 = (2. - 9. A + 12. B + 36. A B - 12. B^2. - 
      56. B^3. + \[Sqrt](4. term1^3. + (2. - 9. A + 12. B + 36. A B - 
            12. B^2. - 56. B^3.)^2.))^(1./3.);
  Zroot1 = (1. - B)/3. - (2.^(1./3.) term1)/(3. term2) + 
    1./(3. 2.^(1./3.)) term2;
  Zroot2 = (1. - B)/
     3. + ((1. + I Sqrt[3.]) term1)/(3. 2.^(2./3.) term2) - 
    1./(6. 2.^(1./3.)) (1. - I Sqrt[3.]) term2;
  Zroot3 = (1. - B)/
     3. + ((1. - I Sqrt[3.]) term1)/(3. 2.^(2./3.) term2) - 
    1./(6. 2.^(1./3.)) (1. + I Sqrt[3.]) term2;
  {Zroot1, Zroot2, Zroot3}]

For this I copied the Module[] from above, and changed the arguments to _Complex to avoid a numerical error:

ZrootsAnalyticCC = 
 Compile[{{A, _Complex}, {B, _Complex}}, 
  Module[{Zroot1, Zroot2, Zroot3, term1, term2}, 
   term1 = (-1. + 3. A - 4. B - 10. B^2.);
   term2 = (2. - 9. A + 12. B + 36. A B - 12. B^2. - 
       56. B^3. + \[Sqrt](4. term1^3. + (2. - 9. A + 12. B + 
             36. A B - 12. B^2. - 56. B^3.)^2.))^(1./3.);
   Zroot1 = (1. - B)/3. - (2.^(1./3.) term1)/(3. term2) + 
     1./(3. 2.^(1./3.)) term2;
   Zroot2 = (1. - B)/
      3. + ((1. + I Sqrt[3.]) term1)/(3. 2.^(2./3.) term2) - 
     1./(6. 2.^(1./3.)) (1. - I Sqrt[3.]) term2;
   Zroot3 = (1. - B)/
      3. + ((1. - I Sqrt[3.]) term1)/(3. 2.^(2./3.) term2) - 
     1./(6. 2.^(1./3.)) (1. + I Sqrt[3.]) term2;
   {Zroot1, Zroot2, Zroot3}], {{Zroot1, _Complex, 
    0}, {Zroot2, _Complex, 0}, {Zroot3, _Complex, 
    0}, {term1, _Complex, 0}, {term2, _Complex, 0}}, 
  CompilationTarget -> "C", RuntimeAttributes -> {Listable}, 
  RuntimeOptions -> "Speed", 
  CompilationOptions -> {"InlineCompiledFunctions" -> True}]

Then running it, showing the speed-up:

Timing[Table[ZrootsAnalytic[x, y], {x, 1., 100.}, {y, 1., 1000.}];]
{7.57591, Null}
Timing[Table[ZrootsAnalyticCC[x, y], {x, 1., 100.}, {y, 1., 1000.}];]
{0.369698, Null}

Making the arguments _Complex is a bit of a hammer. It would be better to see where the computation requires complex calculations and to create a complex sub-expression at that point, in order to avoid the added work of complex calculations on real numbers.

You can see what it is compiling to by using:

Needs["CompiledFunctionTools`"]
CompilePrint[ZrootsAnalyticCC]

It will show the expressions so you see where in the calculation possible errors are occurring. It will also alert you if for some reason the code is falling back to Mathematica for some parts, indicated by MainEvaluate[]'s in the compiled code. You want to find a way to get rid of those in order to realize the benefit of compilation.

$\endgroup$
  • $\begingroup$ When I wrap the expressions in Module, and localize the variables I get an error message "Numerical error encountered; proceeding with uncompiled evaluation". Do I keep the third argument? Removing it seems to make no difference. $\endgroup$ – Robert Mar 3 '16 at 14:59
  • $\begingroup$ You need to provide in your question how you are testing it. I found no such problem, but I don't know what arguments you are giving the function. For the function itself, all I did was copy the Module[] in your first function and paste it as the second argument of Compile[] in your second function. Worked fine, and fast. $\endgroup$ – Mark Adler Mar 3 '16 at 15:28
  • $\begingroup$ Mark I've modified the post above with more explanation and some conclusions. If you are not using _Complex types for A and B then can you post what you are evaluating and the speed comparison. $\endgroup$ – Robert Mar 3 '16 at 16:22
  • $\begingroup$ Ah, yes, I forgot that I also changed the arguments to _Complex. $\endgroup$ – Mark Adler Mar 3 '16 at 17:49
  • $\begingroup$ Thanks Mark, CompilePrint was my inspiration for changing arguments to complex. There may be some overhead with treating A and B as complex but for now I'll live with that. $\endgroup$ – Robert Mar 3 '16 at 19:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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