Probably a hard question, but I think it's better to cry out loud.

I've hesitated for a while about whether I should post this in StackOverflow with a c tag or not, but finally decide to keep it here.

This question can be viewed as a follow up of Has this implementation of FDM touched the speed limit of Mathematica?. In the answer under that post, Daniel managed to implement a compiled Mathematica function that's almost as fast (to be more precise, 3/4 as fast) as the one directly implementing with C++, with the help of devectorization,CompilationTarget -> "C", RuntimeOptions -> "Speed" and Compile`GetElement. Since then, this combination has been tested in various samples, and turns out to be quite effective in speeding up CompiledFunction that involves a lot of array element accessing. I do benefit a lot from this technique, nevertheless in the mean time, another question never disappear in my mind, that is:

Why is the CompiledFunction created with the combination above still slower than the one directly writing with C++?

To make the question more clear and answerable, let's use a simpler example. In the answers under this post about Laplacian of a matrix, I create the following function with the technique above:

cLa = Hold@Compile[{{z, _Real, 2}}, 
     Module[{d1, d2}, {d1, d2} = Dimensions@z; 
      Table[z[[i + 1, j]] + z[[i, j + 1]] + z[[i - 1, j]] + z[[i, j - 1]] - 
        4 z[[i, j]], {i, 2, d1 - 1}, {j, 2, d2 - 1}]], CompilationTarget -> C, 
     RuntimeOptions -> "Speed"] /. Part -> Compile`GetElement // ReleaseHold;

and Shutao create one with LibraryLink (which is almost equivalent to writing code directly with C):

src = "
  #include \"WolframLibrary.h\"
  DLLEXPORT int laplacian(WolframLibraryData libData, mint Argc, MArgument *Args, \
MArgument Res) {
      MTensor tensor_A, tensor_B;
      mreal *a, *b;
      mint const *A_dims;
      mint n;
      int err;
      mint dims[2];
      mint i, j, idx;
      tensor_A = MArgument_getMTensor(Args[0]);
      a = libData->MTensor_getRealData(tensor_A);
      A_dims = libData->MTensor_getDimensions(tensor_A);
      n = A_dims[0];
      dims[0] = dims[1] = n - 2;
      err = libData->MTensor_new(MType_Real, 2, dims, &tensor_B);
      b = libData->MTensor_getRealData(tensor_B);
      for (i = 1; i <= n - 2; i++) {
          for (j = 1; j <= n - 2; j++) {
              idx = n*i + j;
              b[idx+1-2*i-n] = a[idx-n] + a[idx-1] + a[idx+n] + a[idx+1] - 4*a[idx];
      MArgument_setMTensor(Res, tensor_B);
      return LIBRARY_NO_ERROR;
lib = CreateLibrary[src, "laplacian"];

lapShutao = LibraryFunctionLoad[lib, "laplacian", {{Real, 2}}, {Real, 2}];

and the following is the benchmark by anderstood:

enter image description here

Why cLa is slower than lapShutao?

Do we really touch the speed limit of Mathematica this time?

Answer(s) addressing the reason for the inferiority of cLa or improving the speed of cLa are both welcomed.


…OK, the example above turns out to be special, as mentioned in the comment below, cLa will be as fast as lapShutao if we extract the LibraryFunction inside it:

cLaCore = cLa[[-1]];

mat = With[{n = 5000}, RandomReal[1, {n, n}]];

cLaCore@mat; // AbsoluteTiming
(* {0.269556, Null} *)

lapShutao@mat; // AbsoluteTiming
(* {0.269062, Null} *)

However, the effect of this trick is remarkable only if the output is memory consuming.

Since I've chosen such a big title for my question, I somewhat feel responsible to add a more general example. The following is the fastest 1D FDTD implementation in Mathematica so far:

fdtd1d = ReleaseHold@
   With[{ie = 200, cg = Compile`GetElement}, 
    Hold@Compile[{{steps, _Integer}}, 
        Module[{ez = Table[0., {ie + 1}], hy = Table[0., {ie}]},
          Do[ez[[j]] += hy[[j]] - hy[[j - 1]], {j, 2, ie}];
          ez[[1]] = Sin[n/10.];
          Do[hy[[j]] += ez[[j + 1]] - ez[[j]], {j, 1, ie}], {n, steps}]; ez], 
        "CompilationTarget" -> "C", "RuntimeOptions" -> "Speed"] /. Part -> cg /. 
     HoldPattern@(h : Set | AddTo)[cg@a__, b_] :> h[Part@a, b]];

fdtdcore = fdtd1d[[-1]];

and the following is an implemenation via LibraryLink (which is almost equivalent to writing code directly with C):

str = "#include \"WolframLibrary.h\"
  #include <math.h>
  DLLEXPORT int fdtd1d(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument \
    MTensor tensor_ez;
    double *ez;
    int i,t;
    const int ie=200,steps=MArgument_getInteger(Args[0]);
    const mint dimez=ie+1;
    double hy[ie];
    libData->MTensor_new(MType_Real, 1, &dimez, &tensor_ez);
    ez = libData->MTensor_getRealData(tensor_ez);
    MArgument_setMTensor(Res, tensor_ez);
    return 0;}

fdtdlib = CreateLibrary[str, "fdtd"];    
fdtdc = LibraryFunctionLoad[fdtdlib, "fdtd1d", {Integer}, {Real, 1}];

test = fdtdcore[10^6]; // AbsoluteTiming
(* {0.551254, Null} *)    
testc = fdtdc[10^6]; // AbsoluteTiming
(* {0.261192, Null} *)

As one can see, the algorithms in both pieces of code are the same, but fdtdc is twice as fast as fdtdcore. (Well, the speed difference is larger than two years ago, the reason might be I'm no longer on a 32 bit machine. )

My C compiler is TDM-GCC 4.9.2, with "SystemCompileOptions"->"-Ofast" set in Mathematica.

  • $\begingroup$ Unrelated to your question it might be easier to write With[{Part = Compile`GetElement}, Compile[. . .] ] instead of using the replacement. $\endgroup$
    – Mr.Wizard
    Commented Nov 27, 2016 at 17:52
  • 1
    $\begingroup$ @Mr.Wizard Yeah, for this specific example With is a better choice, but when things get more complicated, for example, if the code inside Compile involves something like a[[i, j]] = a[[i, j]] + 1, then ReplaceAll seems to be unavoidable. (We need to use a more complicated rule in this case of course, for example With[{cg = Compile`GetElement}, Hold@Compile[……] /. Part -> cg /. HoldPattern@(h : Set | AddTo)[cg@a__, b_] :> h[Part@a, b]]) $\endgroup$
    – xzczd
    Commented Nov 28, 2016 at 2:41
  • 5
    $\begingroup$ C-CompiledFunction expressions carry raw LibraryFunction objects as last elements. You can test speed of cLaLib = Last@cLa to factor out overhead of CompiledFunction and compare only speed difference between handwritten and auto-generated C-code. Part of slowness of CompiledFunction is related to fact that it performs one additional copy of result so this overhead grows linearly with size of result. $\endgroup$
    – jkuczm
    Commented Nov 29, 2016 at 16:57
  • 1
    $\begingroup$ @jkuczm I think you should give an answer. Though this technique doesn't speed up the FDTD example much, it makes cLa as fast as lapShutao! $\endgroup$
    – xzczd
    Commented Nov 30, 2016 at 8:20

2 Answers 2


Using Compile with CompilationTarget->"C" does generate C-Code to be compiled in a generalized way, the resulting code will contain some overhead due to that compared to hand-written code which can easily explain any difference in runtimes. Even for cases where that overhead is minimal or non-existent automatic code generation will always produce something that is very different from what manually written code would look like, so it is no surprise that the runtimes can differ quite a lot.

I did not find a reference which documents that clearly, but I think when using cLa=Compile[...,CompilationTarget->"C"] what is actually compiled is what you get with:

ExportString[cLa, "C"]

If you look at the result, you clearly find that the generated code is very different from what you would write manually, it more looks like some intermediate state on the way to compile to the WVM: loops are changed into gotos, data access is at a very low level. If you look at the code, it is easy to imagine that the compiler will have a harder time to optimize that code than it has for the nested loop in Shutaos code. The generated code also seems to switch between the one-based Mathematica indices and the zero-based C-indexing at the innermost loop level, using a block local temporary variable. That alone might explain (some of) the runtime differences. The latter is probably easiest to see for something like this:

tst = Hold@Compile[{{z, _Real, 2}}, Table[z[[k, 1]], {k, Length[z]}], 
   CompilationTarget -> C, RuntimeOptions -> "Speed"] /. 
   Part -> Compile`GetElement // ReleaseHold
ExportString[tst, "C"]
  • $\begingroup$ Thanks for exploring. I should say I'll feel a bit surprising if the code at very low level finally turns out to be the culprit. (I had a vague impression that machine-like code is easier to understand for the compiler.) $\endgroup$
    – xzczd
    Commented Nov 28, 2016 at 2:57
  • 1
    $\begingroup$ @xzczd understanding which kind of code is easy to optimize for a compiler is certainly out of my knowledge. I know it was hard 20 years ago and I think it hasn't become easier today. All I know is that it was (and still is) possible to see tremendous differences in runtimes if you manage to formulate code in a way which lets the compiler see its opportunities to optimize its output. And that means that lower level code not automatically is easier to optimize than more abstract formulations. That is by the way also true for formulating Mathematica code on a different abstraction level. $\endgroup$ Commented Nov 28, 2016 at 9:29
  • $\begingroup$ 千金市骨,hope someone will come up with a complete answer someday :) $\endgroup$
    – xzczd
    Commented Dec 15, 2016 at 11:36
  • $\begingroup$ @xzczd: I see my answer more to be a "no answer" so I feel a bit stupid to get a bounty for it, especially as the other answers did put more effort in their analysis. As for your hope - while I always am interested to learn new things I'm afraid that a more detailed analysis would be quite some effort, need high skills in understanding compiler functionality and output AND still be specific to this very code (and probably even for the target processor/platform)... $\endgroup$ Commented Dec 15, 2016 at 15:37

This is only a partial answer, but it was too long for a comment.

It seems that you can slightly increase the speed of the MMA code generated for cLa, by dealing in details with the options of Compile. Let us introduce

cLabis = Hold@Compile[{{z, _Real, 2}}, 
 Module[{d1, d2}, {d1, d2} = Dimensions@z; 
  Table[z[[i + 1, j]] + z[[i, j + 1]] + z[[i - 1, j]] + z[[i, j - 1]] - 
    4 z[[i, j]], {i, 2, d1 - 1}, {j, 2, d2 - 1}]], CompilationTarget -> "C", 
CompilationOptions -> {"ExpressionOptimization" -> True, 
"InlineCompiledFunctions" -> True, 
"InlineExternalDefinitions" -> 
True}, RuntimeOptions -> {"CatchMachineOverflow" -> False , 
"CatchMachineUnderflow" -> False, 
"CatchMachineIntegerOverflow" -> False, 
"CompareWithTolerance" -> False, "EvaluateSymbolically" -> False, 
"WarningMessages" -> False, 
"RuntimeErrorHandler" -> Function[Throw[$Failed]]}
] /. Part -> Compile`GetElement // ReleaseHold;

In particular, we note that we turned off the EvaluateSymbolically option, which is not turned off by the choice RuntimeOptions -> "Speed". I also ensured that all the parameters of CompilationOptions were specified.

One can then compare the runtimes of these various functions using (40s to run this)

compare[n_] := 
Block[{mat = RandomReal[10, {n, n}]}, 
d2 = SparseArray@
N@Sum[NDSolve`FiniteDifferenceDerivative[i, {#, #} &[Range[n]], 
    "DifferenceOrder" -> 2][
   "DifferentiationMatrix"], {i, {{2, 0}, {0, 2}}}];
{AbsoluteTiming[Array[cLa[mat] &, 10];], 
AbsoluteTiming[Array[lapShutao[mat] &, 10];], 
AbsoluteTiming[Array[cLabis[mat] &, 10];]}[[All, 1]]]

tab = Table[{Floor[1.3^i], #} & /@ compare[Floor[1.3^i]], {i, 6, 

PlotLegends -> {"cLa", "Shutao", "cLabis"}, 
AxesLabel -> {"Size", "Time"}, PlotRange -> Full]

In the end, you get a comparison of the form Time comparisons

Unfortunately, the gain only remains very minor...

  • $\begingroup$ I think the d2 inside compare is redundant? $\endgroup$
    – xzczd
    Commented Nov 30, 2016 at 2:54
  • $\begingroup$ @xzczd I naively copy-pasted the code for compare from the associated question... $\endgroup$
    – jibe
    Commented Nov 30, 2016 at 18:11

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.