# how to Compile user-defined function like this

the following function can be compiled successfully

Clear[mat]
mat = Compile[{e, kx, {ii, _Integer}, {label, _Integer}, η},
Inverse[{{e - I η,
-E^(-I kx), 0, -E^(((I kx)/2)), 0, 0, 0, 0, 0, 0, 0, 0, 0,
0}, {-E^(I kx), e - I η, -E^(-((I kx)/2)), 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0}, {0, -E^(((I kx)/2)),
e - I η, -E^(-I kx), 0, -E^(((I kx)/2)), 0, 0, 0, 0, 0, 0,
0, 0}, {-E^(-((I kx)/2)), 0, -E^(I kx),
e - I η, -E^(-((I kx)/2)), 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0,
0, 0, -E^(((I kx)/2)), e - I η, -E^(-I kx),
0, -E^(((I kx)/2)), 0, 0, 0, 0, 0, 0}, {0, 0, -E^(-((I kx)/2)),
0, -E^(I kx), e - I η, -E^(-((I kx)/2)), 0, 0, 0, 0, 0, 0,
0}, {0, 0, 0, 0, 0, -E^(((I kx)/2)), e - I η, -E^(-I kx),
0, -E^(((I kx)/2)), 0, 0, 0, 0}, {0, 0, 0, 0, -E^(-((I kx)/2)),
0, -E^(I kx), e - I η, -E^(-((I kx)/2)), 0, 0, 0, 0,
0}, {0, 0, 0, 0, 0, 0, 0, -E^(((I kx)/2)),
e - I η, -E^(-I kx), 0, -E^(((I kx)/2)), 0, 0}, {0, 0, 0,
0, 0, 0, -E^(-((I kx)/2)), 0, -E^(I kx),
e - I η, -E^(-((I kx)/2)), 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0,
0, 0, -E^(((I kx)/2)), e - I η, -E^(-I kx),
0, -E^(((I kx)/2))}, {0, 0, 0, 0, 0, 0, 0, 0, -E^(-((I kx)/2)),
0, -E^(I kx), e - I η, -E^(-((I kx)/2)), 0}, {0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, -E^(((I kx)/2)),
e - I η, -E^(-I kx)}, {0, 0, 0, 0, 0, 0, 0, 0, 0,
0, -E^(-((I kx)/2)), 0, -E^(I kx), e - I η}}][[
2*(ii - 1) + label, 2*(ii - 1) + label]], CompilationTarget -> "C"]


as you can see the function calc the one of the matrix element of the inverse matrix.

the result of the above code is

CompiledFunction[{e,kx,ii,label,\[Eta]},Block[{Compile$498,Compile$502,Compile$503,Compile$507,Compile$508,Compile$509,Compile$510,Compile$504,Compile$505,Compile$506,Compile$514,Compile$515,Compile$516,Compile$511,Compile$512,Compile$513,Compile$518,Compile$519,Compile$520},Compile$498=I \[Eta];Compile$502=-Compile$498;Compile$503=e+Compile$502;Compile$507=1/2;Compile$508=I kx Compile$507;Compile$509=E^Compile$508;Compile$510=-Compile$509;Compile$504=-I kx;Compile$505=E^Compile$504;Compile$506=-Compile$505;Compile$514=-Compile$508;Compile$515=E^Compile$514;Compile$516=-Compile$515;Compile$511=I kx;Compile$512=E^Compile$511;Compile$513=-Compile$512;Compile$518=ii-1;Compile$519=2 Compile$518;Compile$520=Compile$519+label;Inverse[{{Compile$503,Compile$506,0,Compile$510,0,0,0,0,0,0,0,0,0,0},{Compile$513,Compile$503,Compile$516,0,0,0,0,0,0,0,0,0,0,0},{0,Compile$510,Compile$503,Compile$506,0,Compile$510,0,0,0,0,0,0,0,0},{Compile$516,0,Compile$513,Compile$503,Compile$516,0,0,0,0,0,0,0,0,0},{0,0,0,Compile$510,Compile$503,Compile$506,0,Compile$510,0,0,0,0,0,0},{0,0,Compile$516,0,Compile$513,Compile$503,Compile$516,0,0,0,0,0,0,0},{0,0,0,0,0,Compile$510,Compile$503,Compile$506,0,Compile$510,0,0,0,0},{0,0,0,0,Compile$516,0,Compile$513,Compile$503,Compile$516,0,0,0,0,0},{0,0,0,0,0,0,0,Compile$510,Compile$503,Compile$506,0,Compile$510,0,0},{0,0,0,0,0,0,Compile$516,0,Compile$513,Compile$503,Compile$516,0,0,0},{0,0,0,0,0,0,0,0,0,Compile$510,Compile$503,Compile$506,0,Compile$510},{0,0,0,0,0,0,0,0,Compile$516,0,Compile$513,Compile$503,Compile$516,0},{0,0,0,0,0,0,0,0,0,0,0,Compile$510,Compile$503,Compile$506},{0,0,0,0,0,0,0,0,0,0,Compile$516,0,Compile$513,Compile$503}}][[Compile$520,Compile$520]]],-CompiledCode-]


I want to speed it a little, So I defined a function

 Clear[invmat]
invmat[m_?MatrixQ, i_, j_] :=
Det[Drop[m, {j}, {i}]]*(-1)^(i + j)/Det[m]


which gives (i,j)th element of the inverse matrix of m directly with calc the whole matrix inverse.

and then I compile the function below

Clear[mat2]
mat2 = Compile[{e, kx, {ii, _Integer}, {label, _Integer}, η},
invmat[{{e - I η,
-E^(-I kx), 0, -E^(((I kx)/2)), 0, 0, 0, 0, 0, 0, 0, 0, 0,
0}, {-E^(I kx), e - I η, -E^(-((I kx)/2)), 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0}, {0, -E^(((I kx)/2)), e - I η, -E^(-I kx),
0, -E^(((I kx)/2)), 0, 0, 0, 0, 0, 0, 0, 0}, {-E^(-((I kx)/2)),
0, -E^(I kx), e - I η, -E^(-((I kx)/2)), 0, 0, 0, 0, 0, 0,
0, 0, 0}, {0, 0, 0, -E^(((I kx)/2)), e - I η, -E^(-I kx),
0, -E^(((I kx)/2)), 0, 0, 0, 0, 0, 0}, {0, 0, -E^(-((I kx)/2)),
0, -E^(I kx), e - I η, -E^(-((I kx)/2)), 0, 0, 0, 0, 0, 0,
0}, {0, 0, 0, 0, 0, -E^(((I kx)/2)), e - I η, -E^(-I kx),
0, -E^(((I kx)/2)), 0, 0, 0, 0}, {0, 0, 0, 0, -E^(-((I kx)/2)),
0, -E^(I kx), e - I η, -E^(-((I kx)/2)), 0, 0, 0, 0, 0}, {0,
0, 0, 0, 0, 0, 0, -E^(((I kx)/2)), e - I η, -E^(-I kx),
0, -E^(((I kx)/2)), 0, 0}, {0, 0, 0, 0, 0, 0, -E^(-((I kx)/2)),
0, -E^(I kx), e - I η, -E^(-((I kx)/2)), 0, 0, 0}, {0, 0, 0,
0, 0, 0, 0, 0, 0, -E^(((I kx)/2)), e - I η, -E^(-I kx),
0, -E^(((I kx)/2))}, {0, 0, 0, 0, 0, 0, 0, 0, -E^(-((I kx)/2)),
0, -E^(I kx), e - I η, -E^(-((I kx)/2)), 0}, {0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, -E^(((I kx)/2)),
e - I η, -E^(-I kx)}, {0, 0, 0, 0, 0, 0, 0, 0, 0,
0, -E^(-((I kx)/2)), 0, -E^(I kx), e - I η}},
2*(ii - 1) + label, 2*(ii - 1) + label], CompilationTarget -> "C"]


But mat2 cannot be compiled because the output is the same. I think at least it should compile the matrix variable in invmat

the output of mat2 is

CompiledFunction[{e,kx,ii,label,η},invmat2[{{e-I η,-E^(-I kx),0,-E^(((I kx)/2)),0,0,0,0,0,0,0,0,0,0},{-E^(I kx),e-I η,-E^(-(1/2) (I kx)),0,0,0,0,0,0,0,0,0,0,0},{0,-E^(((I kx)/2)),e-I η,-E^(-I kx),0,-E^(((I kx)/2)),0,0,0,0,0,0,0,0},{-E^(-(1/2) (I kx)),0,-E^(I kx),e-I η,-E^(-(1/2) (I kx)),0,0,0,0,0,0,0,0,0},{0,0,0,-E^(((I kx)/2)),e-I η,-E^(-I kx),0,-E^(((I kx)/2)),0,0,0,0,0,0},{0,0,-E^(-(1/2) (I kx)),0,-E^(I kx),e-I η,-E^(-(1/2) (I kx)),0,0,0,0,0,0,0},{0,0,0,0,0,-E^(((I kx)/2)),e-I η,-E^(-I kx),0,-E^(((I kx)/2)),0,0,0,0},{0,0,0,0,-E^(-(1/2) (I kx)),0,-E^(I kx),e-I η,-E^(-(1/2) (I kx)),0,0,0,0,0},{0,0,0,0,0,0,0,-E^(((I kx)/2)),e-I η,-E^(-I kx),0,-E^(((I kx)/2)),0,0},{0,0,0,0,0,0,-E^(-(1/2) (I kx)),0,-E^(I kx),e-I η,-E^(-(1/2) (I kx)),0,0,0},{0,0,0,0,0,0,0,0,0,-E^(((I kx)/2)),e-I η,-E^(-I kx),0,-E^(((I kx)/2))},{0,0,0,0,0,0,0,0,-E^(-(1/2) (I kx)),0,-E^(I kx),e-I η,-E^(-(1/2) (I kx)),0},{0,0,0,0,0,0,0,0,0,0,0,-E^(((I kx)/2)),e-I η,-E^(-I kx)},{0,0,0,0,0,0,0,0,0,0,-E^(-(1/2) (I kx)),0,-E^(I kx),e-I η}},2 (ii-1)+label,2 (ii-1)+label],-CompiledCode-]


So how to Compile mat2?

• You do know, that neither Inverse nor Det can be compiled down because they are heavily optimised? Please see this post to get a list of compilable functions. Det is a special case because although it appears in the second list of Oleksandr's answer, it cannot be compiled into a function free of a MainEvaluate call. Therefore, all your efforts probably won't lead to a speed up. – halirutan Aug 27 '13 at 8:30
• finally someone replied. Thank you!. Well actually I don't expect to compile Inverse. What I want to compile is the matrix, you see there is a lot similar term in it. And I have test it, mat is compiled, you can see it if you copy the code and paste it into Mathematica and run. But mat2 didn't compile – matheorem Aug 27 '13 at 8:39
• @halirutan I have edited my question – matheorem Aug 27 '13 at 8:44
• While mat indeed can be compiled, it WILL call MainEvaluate when reaches Inverse in the code, you can check this with CompilePrint@mat. Concerning mat2, you cannot compile without MainEvaluate because invmat still uses Det which cannot be compiled, no matter how hard one tries to mask it. – István Zachar Aug 27 '13 at 9:04
• possible duplicate of Compiling more functions that don't call MainEvaluate – István Zachar Aug 27 '13 at 9:31

The problem you face here is that Attributes[Compile] shows you that Compile holds its arguments. Therefore, nothing will be evaluated per default. So what you have to do is to inject/implant the function definition into the body of Compile.

This sounds harder than it is and there are several ways. One of them is (1) holding off the compilation, (2) replace invmat with its definition (3) release the hold and do the compilation. I will show you step by step how it works. First we take a minimal working example which you should have done in the first place, because I'm sure many users were afraid of the big expression in your example. So here is how you prevent evaluation with Hold

Hold[fc = Compile[{}, invmat[{{1, 0}, {0, 1}}]]]
(* Hold[fc = Compile[{}, invmat[{{1, 0}, {0, 1}}]]] *)


Now, you take the body of your definition and replace it

% /. invmat[m_?MatrixQ] :> Inverse[m]
(* Hold[fc = Compile[{}, Inverse[{{1, 0}, {0, 1}}]]] *)


And as last step, you release the Hold

ReleaseHold[%]


If you have already a function definition for the code you want to inject, then you can easily use the DownValues because this returns the internally used rule(s)$^1$. Using this would look like the following

invmat[m_?MatrixQ] := Inverse[m];
Hold[fc = Compile[{}, invmat[{{1, 0}, {0, 1}}]]] /. DownValues[invmat]


Now you can do the same with your original example.

$^1$ Beware, that you might run into problems with the pattern definitions and held code. One simple example is the following

Hold[Module[{m = {{1, 0}, {0, 1}}},
invmat[m]
]] /. DownValues[invmat]


here invmat is not replaced, because although m is/willBe a matrix, in the code it is currently only a symbol which doesn't fulfil MatrixQ.

• Hi, Halirutan! Thank you for your reply!. I was reading mathematica doc today. And found that there exists simpler solution. Just set CompilationOptions -> {"ExpressionOptimization" -> True} will do all the work. – matheorem Aug 28 '13 at 7:17
• @matheorem You are mistaken. This option only optimizes expressions, it does not include the body of your function in the compiled function. Therefore, you (supposed to be) fast compiled function needs to ask the kernel through a MainEvaluate call when this function is called. Please check the output of this small example to see, that f is still called and not replaced with the function definition inside the compiled function. – halirutan Aug 28 '13 at 22:47