184
$\begingroup$

What are the best practices of compiling functions? I understand that this is a vague question, but let me list some aspects that might trigger useful answers. Some of these have already been answered (these don't need answers) but I included them for sake of clarity and completeness:

  1. Which functions are allowed in Compile (see these lists: functions, distributions)?
  2. There are structures which are not allowed by the Mathematica virtual machine, e.g. boolean tensors are not supported (see this question, particularly Oleksandr's comment). Are there other disallowed objects?
  3. Which operations should be left for the main evaluation to do? (e.g. x = x + c1[y] or x = c2[x, y], where the compiled function c2 contains the addition inside).
  4. Which practice is better: partition a large computation into smaller compiled functions, or wrap the whole into one big function? What factors should be considered when making a decision? One such factor is for example that a large compiled function cannot easily return intermediate values (at least it is not trivial how to do this, see this post, and perhaps it also decreases performance).
  5. Type-conversions: Compile silently converts Integer/Real input to the appropriate type, and it also packs non packed arrays, but I'm not sure what the cost of it is. Is it efficient/safe to leave the compiler to deal with type-conversions?
  6. How to call for custom functions in the compiled code? This involves calls for external user functions, recursive calls (perhaps recursive calls passing arguments by reference), or directly injecting code into a held Compile.
  7. Anything else...?

As I'm in the process of understanding Compile, I have a very basic understanding at the moment. Please feel free to add any knowledge to this post.

$\endgroup$
2
  • 1
    $\begingroup$ One thing which I feel compelled to mention, but doesn't seem worth a whole answer, is that Goto won't work right in Compile, at least as of version 8. It's kind of annoying since otherwise Compile's limited subset of Mathematica would be a great fit for crufty old Fortran code. $\endgroup$
    – Pillsy
    Commented Feb 27, 2012 at 20:26
  • 2
    $\begingroup$ Here is a talk from wolfram research conference I think very useful: Effective Use of the Mathematica Compiler and Code Generation wolfram.com/broadcast/video.php?channel=101&video=751 $\endgroup$ Commented Mar 1, 2013 at 23:00

6 Answers 6

154
$\begingroup$

I'll just throw in a few random thoughts in no particular order, but this will be a rather high-level view on things. This is necessarily a subjective exposition, so treat it as such.

Typical use cases

In my opinion, Compile as an efficiency-boosting device is effective in two kinds of situations (and their mixes):

  • The problem is solved most efficiently with a procedural style, because for example an efficient algorithm for it is formulated procedurally and does not have a simple / efficient functional counterpart (note also that functional programming in Mathematica is peculiar in many respects, reflecting the fact that functional layer is a thin one on top of the rule-based engine. So, some algorithms which are efficient in other functional languages may be inefficient in Mathematica). A very clear sign of it is when you have to do array indexing in a loop.

  • The problem can be solved by joining several Compile-able built-in functions together, but there are (perhaps several) "joints" where you face the performance-hit if using the top-level code, because it stays general and can not use specialized versions of these functions, and for a few other reasons. In such cases, Compile merely makes the code more efficient by effectively type-specializing to numerical arguments and not using the main evaluator. One example that comes to mind is when we compile Select with a custom (compilable) predicate and can get a substantial performance boost (here is one example).

I use this rule of thumb when determining whether or not I will benefit from Compile: the more my code inside Compile looks like C code I'd write otherwise, the more I benefit from it (strictly speaking, this is only true for the compilation to C, not MVM).

It may happen that some portions of top-level code will be the major bottleneck and can not be recast into a more efficient form, for a given approach to the problem. In such a case, Compile may not really help, and it is best to rethink the whole approach and try to find another formulation for the problem. In other words, it often saves time and effort to do some profiling and get a good idea about the real places where the bottlenecks are, before turning to Compile.

Limitations of Compile

Here is an (incomplete) list of limitations, most of which you mentioned yourself

  • Can only accept regular arrays (tensors) of numerical or boolean types. This excludes ragged arrays and more general Mathematica expressions.
  • In most cases, can only return a single tensor of some type
  • Only machine-precision arithmetic
  • From the user-defined functions, only pure functions are compilable, plus one can inline other compiled functions. Rules and "functions" defined with rules are inherently not compilable.
  • No way to create functions with memory (a-la static variables in C)
  • Only a small subset of built-in functions can be compiled to byte-code (or C)
  • Possibilities for writing recursive compiled functions seem to be very limited, and most interesting cases seem to be ruled out
  • No decent pass-by-reference semantics, which is a big deal (to me anyways)
  • You can not really use indexed variables in Compile, although it may appear that you can.
  • ...

Whether or not to compile to C?

I think this depends on the circumstances. Compilation to C is expensive, so this makes sense only for performance-critical code to be used many times. There are also many cases when compilation to MVM will give similar performance, while being much faster. One such example can be found in this answer, where the just-in-time compilation to MVM target led to a major speed-up, while compilation to C would have likely destroyed the purpose of it - in that particular case.

Another class of situations when compiling to C is may not be the best option is when you want to "serialize" the CompiledFunction object, and distribute it to others, for example in a package, and you don't want to count on a C compiler being installed on the user's machine. As far as I know, there is no automatic mechanism yet to grab the generated shared library and package it together with the CompiledFunction, and also one would have to cross-compile for all platforms and automatically dispatch to the right library to load. All this is possible but complicated, so, unless the speed gain can justify such complications for a given problem, it may be not worth it, while compilation to MVM target creates the top-level CompiledFunction object, which is automatically cross-platform, and does not require anything (except Mathematica) to be installed.

So, it really depends, although more often than not compilation to C will lead to faster execution and, if you at all decide to use Compile, will be justified.

What to include in Compile

I share an opinion that, unless you have some specific requirements, it is best to only use Compile on minimal code fragments which would benefit from it the most, rather than have one big Compile. This is good because:

  • It allows you to better understand where the real bottlenecks are
  • It makes your compiled code more testable and composable
  • If you really need it, you can then combine these pieces and use "InlineCompiledFunctions" -> True option setting, to get all the benefits that one large Compile would give you
  • Since Compile is limited in what it can take, you will have less headaches on how to include some uncompilable pieces, plus less chances to overlook a callback to the main evaluator

That said, you may benefit from one large Compile in some situations, including:

  • Cases when you want to grab the resulting C code and use it stand-alone (linked against Wolfram RTL)
  • Cases when you want to run your compiled code in parallel on several kernels and don't want to think about possible definitions distribution issues etc (this was noted by @halirutan)

Listable functions

When you can, it may be a good idea to use the RuntimeAttributes -> Listable option, so that your code can be executed on (all or some) available cores in parallel. I will give one example which I think is rather interesting, because it represents a problem which may not initially look like one directly amenable to this (although it is surely not at all hard to realize that parallelization may work here) - computation of Pi as a partial sum, of a well-known infinite sum representation. Here is a single-core function:

Clear[numpi1];
numpi1 = 
   Compile[{{nterms, _Integer}}, 
      4*Sum[(-1)^k/(2 k + 1), {k, 0, nterms}], 
        CompilationTarget -> "C", RuntimeOptions -> "Speed"];

Here is a parallel version:

numpiParallelC = 
  Compile[{{start, _Integer}, {end, _Integer}}, 
    4*Sum[(-1)^k/(2 k + 1), {k, start, end}], CompilationTarget -> "C",
       RuntimeAttributes -> Listable, RuntimeOptions -> "Speed"];

Clear[numpiParallel];
numpiParallel[nterms_, nkernels_] := 
  Total@Apply[numpiParallelC, 
     MapAt[# + 1 &, {Most[#], Rest[#] - 1}, {2, -1}] &@
        IntegerPart[Range[0, nterms, nterms/nkernels]]];

Now, some benchmarks (on a 6-core machine):

(res0=numpiParallel[10000000,1])//AbsoluteTiming
(res1=numpiParallel[10000000,6])//AbsoluteTiming
(res2=numpi1[10000000])//AbsoluteTiming
Chop[{res0-res2,res0-res1,res2-res1}]

(*
 ==>
 {0.0722656,3.14159}
 {0.0175781,3.14159}
 {0.0566406,3.14159}
 {0,0,0}
*)

A few points to note here:

  • It may happen that the time it takes to prepare the data to be fed into a Listable compiled function, will be much more than the time the function runs (e.g. when we use Transpose or Partition etc on huge lists), which then sort of destroys the purpose. So, it is good to make an estimate whether or not that will be the case.
  • A more "coarse-grained" alternative to this is to run a single-threaded compiled function in parallel on several Mathematica kernels, using the built-in parallel functionality (ParallelEvaluate, ParallelMap, etc). These two possibilities are useful in different situations.

Auto-compilation

While this is not directly related to the explicit use of Compile, this topic logically belongs here. There are a number of built-in (higher-order) functions, such as Map, which can auto-compile. What this means is that when we execute

Map[f, list]

the function f is analyzed by Map, which attempts to automatically call Compile on it (this is not done at the top-level, so using Trace won't show an explicit call to Compile). To benefit from this, the function f must be compilable. As a rule of thumb, it has to be a pure function for that (which is not by itself a sufficient condition) - and generally the question of whether or not a function is compilable is answered here in the same way as for explicit Compile. In particular, functions defined by patterns will not benefit from auto-compilation, which is something to keep in mind.

Here is a little contrived but simple example to illustrate the point:

sumThousandNumbers[n_] := 
   Module[{sum = 0}, Do[sum += i, {i, n, n + 1000}]; sum]

sumThousandNumbersPF = 
   Module[{sum = 0}, Do[sum += i, {i, #, # + 1000}]; sum] &

Now, we try:

Map[sumThousandNumbers, Range[3000]]//Short//Timing
Map[sumThousandNumbersPF, Range[3000]]//Short//Timing

(*
  ==> {3.797,{501501,502502,503503,504504,505505,<<2990>>,3499496,
               3500497,3501498,3502499,3503500}}

      {0.094,{501501,502502,503503,504504,505505,<<2990>>,3499496,
               3500497,3501498,3502499,3503500}}
*)

which shows a 40-times speedup in this particular case, due to auto-compilation.

There are in fact many cases when this is important, and not all of them are as obvious as the above example. One such case was considered in a recent answer to the question of extracting numbers from a sorted list belonging to some window. The solution is short and I will reproduce it here:

window[list_, {xmin_, xmax_}] := 
    Pick[list, Boole[xmin <= # <= xmax] & /@ list, 1]

What may look like a not particularly efficient solution, is actually quite fast due to the auto-compilation of the predicate Boole[...] inside Map, plus Pick being optimized on packed arrays. See the aforementioned question for more context and discussion.

This shows us another benefit of auto-compilation: not only does it often make the code run much faster, but it also does not unpack, allowing surrounding functions to also benefit from packed arrays when they can.

Which functions can auto-compile? One way to find out is to inspect SystemOptions["CompileOptions"]:

Cases["CompileOptions"/.SystemOptions["CompileOptions"],
      opt:(s_String->_)/;StringMatchQ[s,__~~"Length"]]

{"ApplyCompileLength" -> \[Infinity], "ArrayCompileLength" -> 250, 
 "FoldCompileLength" -> 100, "ListableFunctionCompileLength" -> 250, 
 "MapCompileLength" -> 100, "NestCompileLength" -> 100, 
 "ProductCompileLength" -> 250, "SumCompileLength" -> 250, 
 "TableCompileLength" -> 250}

This also tells you the threshold lengths of the list beyond which the auto-compilation is turned on. You can also change these values. Setting the value of ...CompileLength to Infinity is effectively disabling the auto-compilation. You can see that "ApplyCompileLength" has this value. This is because it can only compile 3 heads: Times, Plus, and List. If you have one of those in your code, however, you can reset this value, to benefit from auto-compilation. Generally, the default values are pretty meaningful, so it is rarely necessary to change these defaults.

A few more techniques

There are a number of techniques involving Compile, which are perhaps somewhat more advanced, but which sometimes allow one to solve problems for which plain Compile is not flexible enough. Some which I am aware of:

  • Sometimes you can trade memory for speed, and, having a nested ragged list, pad it with zeros to form a tensor, and pass that to Compile.

  • Sometimes your list is general and you can not directly process it in Compile to do what you want, however, you can reformulate a problem such that you can instead process a list of element positions, which are integers. I call it "element-position duality". One example of this technique in action is here, for a larger application of this idea see my last post in this thread (I hesitated to include this reference because my first several posts there are incorrect solutions. Note that for that particular problem, a far more elegant and short, but somewhat less efficient solution was given in the end of that thread).

  • Sometimes you may need some structural operations to prepare the input data for Compile, and the data contains lists (or, generally, tensors), of different types (say, integer positions and real values). To keep the list packed, it may make sense to convert integers to reals (in this example), converting them back to integers with IntegerPart inside Compile. One such example is here

  • Run-time generation of compiled functions, where certain run-time parameters get embedded. This may be combined with memoization. One example is here, another very good example is here

  • One can emulate pass-by-reference and have a way of composing larger compiled functions out of smaller ones with parameters (well, sort of), without a loss of efficiency. This technique is showcased for example here

  • A common wisdom is that since neither linked-lists, nor Sow-Reap are compilable, one has to pre-allocate large arrays most of the time, to store the intermediate results. There are at least two other options:

    • Use Internal`Bag, which is compilable (the problem however is that it can not be returned as a result of Compile as of now, AFAIK).
    • It is quite easy to implement an analog of a dynamic array inside your compiled code, by setting up a variable which gives the current size limit, and copy your array to a new larger array once more space is needed. In this way, you only allocate (at the end) as much space as is really needed, for a price of some overhead, which is often negligible.
  • One may often be able to use vectorized operations like UnitStep, Clip, Unitize etc, to replace the if-else control flow in inner loops, also inside Compile. This may give a huge speed-up, particularly when compiling to MVM target. Some examples are in my comments in this and this blog posts, and one other pretty illustrative example of a vectorized binary search in my answer in this thread

  • Using additional list of integers as "pointers" to some lists you may have. Here, I will make an exception for this post, and give an explicit example, illustrating the point. The following is a fairly efficient function to find a longest increasing subsequence of a list of numbers. It was developed jointly by DrMajorBob, Fred Simons and myself, in an on and off-line MathGroup discussion (so this final form is not available publicly AFAIK, thus including it here)

Here is the code

Clear[btComp];
btComp = 
Compile[{{lst, _Integer, 1}}, 
   Module[{refs, result, endrefs = {1}, ends = {First@lst}, 
      len = Length@lst, n0 = 1, n1 = 1, i = 1, n, e}, 
     refs = result = 0 lst;
     For[i = 2, i <= len, i++, 
        Which[
          lst[[i]] < First@ends, 
             (ends[[1]] = lst[[i]]; endrefs[[1]] = i; refs[[i]] = 0),
          lst[[i]] > Last@ends, 
             (refs[[i]] = Last@endrefs;AppendTo[ends, lst[[i]]]; AppendTo[endrefs, i]), 
          First@ends < lst[[i]] < Last@ends, 
             (n0 = 1; n1 = Length@ends;  
              While[n1 - n0 > 1, 
                n = Floor[(n0 + n1)/2];
                If[ends[[n]] < lst[[i]], n0 = n, n1 = n]];
                ends[[n1]] = lst[[i]];
                endrefs[[n1]] = i;
                refs[[i]] = endrefs[[n1 - 1]])
        ]];
        For[i = 1; e = Last@endrefs, e != 0, (i++; e = refs[[e]]), 
            result[[i]] = lst[[e]]];
        Reverse@Take[result, i - 1]], CompilationTarget -> "C"];

Here is an example of use (list should not contain duplicates):

test = RandomSample[#, Length[#]] &@ Union@RandomInteger[{1, 1000000}, 1000000];

btComp[test] // Length // Timing

The fastest solution based on built-ins, which is indeed very fast, is still about 6 times slower for this size of the list:

LongestCommonSequence[test, Sort@test] // Short // Timing

Anyways, the point here is that this was possible because of extra variables refs and endrefs, the use of which allowed to only manipulate single integers (representing positions of sub-lists in a larger list) instead of large integer lists.

A few assorted remarks

  • Things to watch out for: see this discussion for some tips on that. Basically, you should avoid

    • Callbacks to the main evaluator
    • Excessive copying of tensors (CopyTensor instruction)
    • Accidental unpacking happening in top-level functions preparing input for Compile or processing its output. This is not related to Compile proper, but it happens that Compile does not help at all, because the bottleneck is in the top-level code.
  • Type conversion I would not worry about performance hit, but sometimes wrong types may lead to run-time errors, or unanticipated callbacks to MainEvaluate in the compiled code.

  • Certain functions (e.g. Sort with the default comparison function, but not only), don't benefit from compilation much or at all.

  • It is not clear how Compile handles Hold- attributes in compiled code, but there are indications that it does not fully preserve the standard semantics we are used to in the top-level.

  • How to see whether or not you can effectively use Compile for a given problem. My experience is that with Compile in Mathematica you have to be "proactive" (with all my dislike for the word, I know of nothing better here). What I mean is that to use it effectively, you have to search the structure of your problem / program for places where you could transform the (parts of) data into a form which can be used in Compile. In most cases (at least in my experience), except obvious ones where you already have a procedural algorithm in pseudo-code, you have to reformulate the problem, so you have to actively ask: what should I do to use Compile here.

$\endgroup$
17
  • $\begingroup$ I did not at first realize how essential the use of Boole was in this solution. The same thing could be more easily written using a Select, but that is much slower. It seems Boole allowed the auto-compilation to happen, or perhaps the speedup is only because Boole enabled packed arrays as you mentioned in your comment. Maybe it's worth investigating this further and including it in this answer of your, or perhaps here: mathematica.stackexchange.com/q/2/12 $\endgroup$
    – Szabolcs
    Commented Feb 27, 2012 at 12:57
  • $\begingroup$ @Szabolcs Added a section to reflect this. See if you like it. $\endgroup$ Commented Feb 27, 2012 at 13:45
  • $\begingroup$ My point (just didn't have time to test it fully when I wrote the first comment) was that if we take cf = Compile[{{r, _Real, 1}}, Boole[1. < # < 10.] & /@ r], it compiles fine. But if we remove Boole, it does not (it calls back to MainEvaluate). It seems a compiled function can't return a vector of booleans, but Boole is an effective way to work around this. I have never thought of it this way before. (Of course this specific function is better vectorized, as in ruebenko's recent answer on my question, but it was just for illustration) $\endgroup$
    – Szabolcs
    Commented Feb 27, 2012 at 15:07
  • $\begingroup$ @Szabolcs I knew about this, and generally that True and False are not so well integrated into Compile as numeric types. But this is a good point. This may be worth including into my answer above - please feel free to do it, if you are inclined, or I will try to find some time later to do this. $\endgroup$ Commented Feb 27, 2012 at 15:12
  • $\begingroup$ @LeonidShifrin I thought I'd ask you this, given your choice of example here and your academic past — do you have fast implementations of the TW distribution? I know it's not a simple problem to integrate the Painlevé equation and my current implementation takes about 0.4 seconds per point on a 2.6Ghz core i7 processor. This was mostly based on an older notebook, I believe written by Tracy. Since I don't require high precision, but frequent evaluations at arbitrary points memoization isn't mostly helpful. I use interpolation now for speed, but thought I'd check if you had smarter ways... $\endgroup$
    – rm -rf
    Commented Aug 31, 2012 at 19:39
36
$\begingroup$

Setting SetSystemOptions[ "CompileOptions" -> "CompileReportExternal"->True] will emit a message when parts of your function do not get compiled. After compilation, Needs["CompiledFunctionTools`"] followed by CompilePrint[cF] (with cF the function you have compiled will display some bytecode; looking for CopyTensor or MainEvaluate in that helps locate inefficiencies (MainEvaluate calls the main mma kernel, so it effectively means not compiling).

Other useful options may be found by evaluating SystemOptions["CompileOptions"]

EDIT: As suggested by Oleksandr in a comment, On[Compile::noinfo] is also useful.

$\endgroup$
0
36
$\begingroup$

Leonid Shifrin has already given an excellent answer for the question but it's so… long and may be frustrating for someone just beginning to learn the usage of Compile so I decided to post this as an answer.

Recently (OK… actually it's more than a year ago) I found that Ted Ersek's Mathematica tricks(.nb version can be found here) contains a brief but masterly summarization for the limitation of Compile. Though it's written based on Version 4, the 6 rules listed by him are still valid and if you fully understand them, you can handle most compilation issues. For convenience, I'll paste the relevant part below (with modifications for part of the commentaries and examples).

BTW, based on the 6 rules, I wrote a Chinese tutorial for compilation here.


(1) Only a subset of kernel functions can be used in compiled evaluation.

We already have this post for the issue.


(2) Compiled evaluation can't use global variables.

In the next piece of code a compiled function uses a global variable width. Later we see that the function works, but timing tests will show that it's slower than if the function was defined without using Compile.

width=2.5;
g1 = Function[{x}, x + width];
g2 = Compile[{x}, x + width];

Do[g1[0.3], {10^6}]; // AbsoluteTiming
Do[g2[0.3], {10^6}]; // AbsoluteTiming
{0.786027, Null}
{1.303057, Null}

The next piece of code defines the same function as the previous example and takes the global variable as an argument. This is the way global variables should be used in a compiled function.

g3 = Compile[{x, y}, x + y];

Do[g3[0.3, width], {10^6}]; // AbsoluteTiming
{0.415907, Null}

The next piece of code shows what happens if we try to set the value of a global variable inside a compiled function. Here again the function works, but timing tests show that the function is slower than the same function defined without using Compile.

h1 = Function[{x}, (temp = 2.5; x + temp)];
h2 = Compile[{x}, (temp = 2.5; x + temp)];

Do[h1[0.3], {10^6}]; // AbsoluteTiming
Do[h2[0.3], {10^6}]; // AbsoluteTiming
{1.471703, Null}
{2.391085, Null}

In the last example the value of temp was set inside Compile. When something like this is needed, Block, Module, or With should be used to make temp a local variable. Compiled functions are defined below that do the same thing as the last example using Block, Module, With and they all use compiled evaluation which runs much faster. Besides the advantage in speed the approach below ensures that evaluating say ha[0.3], hb[0.3], hc[0.3] will not change the value temp might have outside the compiled function.

ha = Compile[{x}, Block[{temp = 2.5}, x + temp]];
hb = Compile[{x}, Module[{temp = 2.5}, x + temp]];
hc = Compile[{x}, With[{temp = 2.5}, x + temp]];

Do[ha[0.3], {10^6}]; // AbsoluteTiming
Do[hb[0.3], {10^6}]; // AbsoluteTiming
Do[hc[0.3], {10^6}]; // AbsoluteTiming
{0.364617, Null}
{0.352145, Null}
{0.347725, Null}

(3) Compiled evaluation can't work with patterns.

This rule can be viewed as a deduction of rule (1).

The function in next piece of code takes a list of real numbers, and determines if 0.5 is in the list. This function is fully compiled.

f4=Compile[{{lst,_Real,1}},FreeQ[lst,0.5]];
<<CompiledFunctionTools`
CompilePrint@f4
…………
1 B0 = FreeQ[ T(R1)0, T(R0)0, R1]]
2 Return

The function below also takes a list of real numbers, and determines if the list is free of negative numbers. This function works, but timing tests will show that it's no faster than the same function defined without using Compile. The problem in this case is that the compiled function uses the pattern _?Negative. Compiled evaluation isn't used if the second argument in Compile includes patterns of any kind. When faced with a problem like this we should use an algorithm that doesn't require patterns, or define a function without using Compile.

f5=Compile[{{lst,_Real,1}},FreeQ[lst,_?Negative]];
CompilePrint@f5
…………
1 R0 = MainEvaluate[ Function[{lst}, _?Negative][ T(R1)0]]
2 B0 = FreeQ[ T(R1)0, T(R0)0, R1]]
3 Return

The next piece of code shows one way the function from the last example can be written to avoid the use of patterns. This version will use compiled evaluation.

f5fixed=Compile[{{lst,_Real,1}}, FreeQ[Sign[lst],-1]  ];
CompilePrint@f5fixed
…………
1 T(I1)1 = Sign[ T(R1)0]
2 B0 = FreeQ[ T(I1)1, T(I0)0]]
3 Return

Notice that function definition like f[x_] = … or f[x_] := … also uses pattern so they can't be compiled!


(4) Local variables in compiled evaluation must always have the same type.

The next piece of code defines a function which uses a local variable temp where temp starts with an integer value. Later in the function temp is changed to a real number. In this case the function works, but timing tests will show it's no faster than the same function defined without using Compile. Local variables used in a compiled function should always have the same type (Real, Integer, Complex, True|False) to ensure compiled evaluation is used.

f6 = Compile[{x}, Module[{temp = 5}, (temp = temp + x; Round[temp])]];
CompilePrint@f6
…………
1 I1 = I0
2 R1 = I1
3 R1 = R1 + R0
4 V17I1 = MainEvaluate[ Function[{x, tempCompile$2}, 
    Block[{temp = tempCompile$2}, {temp 
= temp + x, temp}]][ R0, I1]]
5 Return

So one can fix f6 in the following ways:

f6fixed1 = Compile[{x}, Module[{temp = 5}, temp = Round[temp + x]]];
CompilePrint@f6fixed1
…………
1 I1 = I0
2 R1 = I1
3 R1 = R1 + R0
4 I2 = Round[ R1]
5 I1 = I2
6 Return
f6fixed2 = Compile[{x}, Module[{temp = 5.}, temp = temp + x; Round[temp]]];
CompilePrint@f6fixed2
…………
1 R2 = R1
2 R3 = R2 + R0
3 R2 = R3
4 I1 = Round[ R2]
5 Return

(5) A compiled function can't change the value of it's argument.

The next piece of code shows an attempt to define a function. If you try to use this function you will see it doesn't work. The problem here is that the function tries to change the value of it's argument, but a compiled function can't change the value of it's argument for the same reason we can't evaluate 3=5.

f7=Compile[{x},
Do[x=Cos[x],{4}];
x];

Some times we want to use an algorithm that starts with some value and replaces it with another value. If we want to do this in a compiled function we should initialize a local variable to the value of the functions argument, and then change the value of the local variable. The function in the next cell does this and uses compiled evaluation.

f7fixed=Compile[{x},
Module[{temp=x},
Do[temp=Cos[temp],{4}];
temp
]];

(6) Compiled evaluation can't work with all possible list structures.

The only list structures compiled evaluation can work with are vectors, matrices, and other tensors. The function defined in the next piece of code returns a list structure that doesn't meet these restrictions so compiled evaluation can't be used. In a case like this there is no advantage to using Compile to define the function.

f8=Compile[{x},{x,{2x,3x}}];
CompilePrint@f8
…………
11    T(R1)2 = MainEvaluate[ Hold[List][ R0, T(R1)1]]
12    Return
$\endgroup$
1
  • 1
    $\begingroup$ sha fa!!! : ) +1 $\endgroup$
    – matheorem
    Commented Jan 14, 2016 at 8:56
13
$\begingroup$

I do a simple test on how Compile translate the list generating function to pseudo-code C. For example, in Mathematica, Range is more efficient than Table, and Table from a known list is better than Table when i runs from 1 to imax. But in Compile, thing changes. Table with i runs from 1 to imax has the same code with Range, so has the same efficiency.

<< "CompiledFunctionTools`"
f1 = Compile[{{n, _Integer}}, Table[i, {i, 1, n}]]
CompilePrint[f1]


        1 argument
        6 Integer registers
        1 Tensor register
        Underflow checking off
        Overflow checking off
        Integer overflow checking on
        RuntimeAttributes -> {}

        I0 = A1
        I1 = 1
        I4 = 0
        Result = T(I1)0

1   I2 = I0
2   I3 = I4
3   T(I1)0 = Table[ I2]
4   I5 = I4
5   goto 7
6   Element[ T(I1)0, I3] = I5
7   if[ ++ I5 < I2] goto 6
8   Return


f2 = Compile[{{n, _Integer}}, Range[n]]
CompilePrint[f2]

        1 argument
        6 Integer registers
        1 Tensor register
        Underflow checking off
        Overflow checking off
        Integer overflow checking on
        RuntimeAttributes -> {}

        I0 = A1
        I2 = 1
        I4 = 0
        Result = T(I1)0

1   I1 = I0
2   I3 = I4
3   T(I1)0 = Table[ I1]
4   I5 = I4
5   goto 7
6   Element[ T(I1)0, I3] = I5
7   if[ ++ I5 < I1] goto 6
8   Return

Compile 'translates' the ideas of Range and Table to the same pseudo-C code. So not always the faster function in Mathematica could run faster in the C-translated code.

f3 = Compile[{{n, _Integer}}, Table[i, {i, Range[n]}]]
CompilePrint[f3]

        1 argument
        7 Integer registers
        4 Tensor registers
        Underflow checking off
        Overflow checking off
        Integer overflow checking on
        RuntimeAttributes -> {}

        I0 = A1
        I2 = 1
        I4 = 0
        Result = T(I1)3

1   I1 = I0
2   I3 = I4
3   T(I1)2 = Table[ I1]
4   I5 = I4
5   goto 7
6   Element[ T(I1)2, I3] = I5
7   if[ ++ I5 < I1] goto 6
8   I3 = Length[ T(I1)2]
9   I1 = I4
10  T(I1)3 = Table[ I3]
11  I5 = I4
12  goto 15
13  I6 = GetElement[ T(I1)2, I5]
14  Element[ T(I1)3, I1] = I6
15  if[ ++ I5 < I3] goto 13
16  Return

In this case, Table from a known list runs faster in Mathematica, but run slower in Compiled code. f3 takes 16 steps while f1, f2 takes only 8 steps.

The more native way to generate a list is shown in f4. It uses only 11 steps, but f4 uses CopyTensor in a loop so should be the slowest.

f4 = Compile[{{n, _Integer}}, 
  Block[{x = {}}, Do[AppendTo[x, i], {i, n}]; x]]
CompilePrint[f4]

    1 argument
    5 Integer registers
    1 Real register
    5 Tensor registers
    Underflow checking off
    Overflow checking off
    Integer overflow checking on
    RuntimeAttributes -> {}

    I0 = A1
    I4 = -1
    I2 = 0
    Result = T(R1)3

1   T(R1)3 ={ }
2   I1 = I0
3   I3 = I2
4   goto 10
5   T(I1)1 ={ I4 }
6   T(I2)2 ={ T(I1)1 }
7   R0 = I3
8   T(R1)4 = Insert[ T(R1)3, T(R0)0, T(I2)2]]
9   T(R1)3 = CopyTensor[ T(R1)4]]
10  if[ ++ I3 < I1] goto 5
11  Return

One can use a higher function in Mathematica to get the same result. f5 uses CopyTensor and takes 30 steps.

f5 = Compile[{{n, _Integer}}, Array[# &, n]]
CompilePrint[f5]

        1 argument
        1 Boolean register
        9 Integer registers
        3 Tensor registers
        Underflow checking off
        Overflow checking off
        Integer overflow checking on
        RuntimeAttributes -> {}

        I0 = A1
        I3 = 1
        I1 = 0
        Result = T(I1)2

1   B0 = I0 <= I1
2   if[ !B0] goto 5
3   Return Error
4   goto 5
5   I2 = I0
6   I4 = I1
7   T(I1)0 = Table[ I2]
8   I5 = I1
9   goto 11
10  Element[ T(I1)0, I4] = I5
11  if[ ++ I5 < I2] goto 10
12  I4 = I0
13  I2 = I1
14  T(I1)1 = Table[ I4]
15  I5 = I1
16  goto 18
17  Element[ T(I1)1, I2] = I5
18  if[ ++ I5 < I4] goto 17
19  T(I1)0 = CopyTensor[ T(I1)1]]
20  T(I1)1 = CopyTensor[ T(I1)0]]
21  I2 = Length[ T(I1)1]
22  I4 = I1
23  T(I1)2 = Table[ I2]
24  I6 = I1
25  goto 29
26  I7 = GetElement[ T(I1)1, I6]
27  I8 = I7
28  Element[ T(I1)2, I4] = I8
29  if[ ++ I6 < I2] goto 26
30  Return

A simple function to counts the number of line in C-code:

FindCCompileLine[f_CompiledFunction] :=
 Length@StringPosition[CompilePrint[f], "\n"];

Check, generate a list from {1,...,10}

f1[10] == f2[10] == f3[10] == f4[10] == f5[10]
True

FindCCompileLine[#] & /@ {f1, f2, f3, f4, f5}

{22, 22, 30, 26, 45}

MyTiming[#] & /@ {f1, f2, f3, f4, f5}

0.780005  0.780005  1.076407   5.054432  1.435209

Conclusions

  • Simple function like Table, Range is more efficient, translated to shorter c-code. f1 f2 uses 8 steps to finish the job. f1, f2 are the fastest.
  • Don't use high order function for a simple task: f5 vs f1 vs f2 vs f3. f5 relies on Array to generate the list {1,2,3,...10} , and needs 30 steps. f3 uses Table with a known list, and needs 16 steps.
  • CopyTensor can be used, but avoid to put CopyTensor in a loop. f4 is slowest, f5 uses CopyTensor too, but not as slow as f4. f4 put the AppendTo (which use CopyTensor to copy to a new list each time) in a Do loop.
  • In general, with a good pseudo code (without CopyTensor in a loop, without MainEvaluate) the longer C code will run slower. f3 vs f1 vs f2.
  • Avoid Append, AppendTo to insert a value in a list, try to use Table to generate a new list with the inserted value.
  • The function translated by Compile, and the function evaluated in Mathematica in the different way. Compile is not fully optimized and could be improved better. Maybe, one can try the add on package Mathcode: http://www.wolfram.com/products/applications/mathcode/ to get a better compiled C code?
  • For list building, use Internal`Bag as mentioned by István Zachar in the comments.
$\endgroup$
3
  • 3
    $\begingroup$ This is a very long post considering how little novel insight it offers. Leonid (and the documentation) has already warned us about using CopyTensor excessively. The number of lines in C does not really correlate with evaluation speed: sometimes it is more fast if steps of a loop are explicitly written (more lines) then put in a loop. For successive list building, Internal`Bag is one really fast solution. A better, more concise summary would be more useful than all the CompilePrint output. $\endgroup$ Commented Feb 15, 2014 at 15:03
  • 2
    $\begingroup$ You are right. I'm trying to do it shorter and more concise, but as you are an expert so maybe your point of view is different with the beginner's view. $\endgroup$
    – Nam Nguyen
    Commented Feb 16, 2014 at 10:28
  • 3
    $\begingroup$ Sure DaoTRINH, you're right about the viewpoints (though I don't consider myself an expert), and it is most important that you did a systematic investigation of the problem! Sorry for the blunt comment, keep it up, such approaches will do good use for the community :) $\endgroup$ Commented Feb 16, 2014 at 14:56
10
$\begingroup$

Not a complete answer by any means but if you are new to Compile I highly recommend that you read the tutorial on the Mathematica compiler. It is very useful.

$\endgroup$
10
$\begingroup$

This answer I gave illustrates how to define the type of local variables used in a Module inside Compile which helps solve some types of issues as compiled functions sometimes are more strongly typed than non compiled ones.

Internal`Bag inside Compile

Inlining external functions inside Compile can sometimes be problematic (because of the need of definining the type of some local variables like above for example).
The post Compiling more functions that don't call MainEvaluate answers this problem for many situations (for all that I have seen) by expanding Mathematica code recursively until it is only composed of compilable functions. The code will then be automatically optimized by Compile.
This is quite interesting because it allows to divide your numerical code more easily into simpler functions in order to do complex things.

$\endgroup$

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.