Why is CompilationTarget -> C slower than directly writing with C?

Question

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;
  }
  ";
Needs["CCompilerDriver`"]
lib = CreateLibrary[src, "laplacian"];

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

and the following is the benchmark by anderstood:

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.

---

Update

…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[
          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 \
Res){
    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);
    
    for(i=0;i<ie+1;i++){ez[i]=0;}   
    for(i=0;i<ie;i++){hy[i]=0;}
  
    for(t=1;t<=steps;t++){
        for(i=1;i<ie;i++){ez[i]+=(hy[i]-hy[i-1]);}
        ez[0]=sin(t/10.);
        for(i=0;i<ie;i++){hy[i]+=(ez[i+1]-ez[i]);}
    }
    
    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.

Answer 1

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"]

Answer 2

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, 
 31}];

ListLinePlot[Transpose@tab, 
PlotLegends -> {"cLa", "Shutao", "cLabis"}, 
AxesLabel -> {"Size", "Time"}, PlotRange -> Full]

In the end, you get a comparison of the form

Unfortunately, the gain only remains very minor...