# Error from Compile with a sum over two indices

I have the following code:

compVDiag = Compile[
{{Nn, _Integer}, {M, _Integer}, {n, _Integer}, {basisvecs, _Integer,
2}, {c1m, _Integer, 2}, {c2m, _Integer, 2}},
Block[
{zpow, wpow, zlen, wlen},
zpow = Union[basisvecs[[n, 1 ;; Nn]]];
wpow = Union[basisvecs[[n, Nn + 1 ;; Nn + M]]];
zlen = Length[zpow];
wlen = Length[wpow];

Sum[
2^(-zpow[[i]] - zpow[[j]])*c1m[[n, zpow[[i]] + 1]]*
c1m[[n, zpow[[j]] + 1]]*
Binomial[zpow[[i]] + zpow[[j]], zpow[[i]]], {i, 1,
zlen - 1}, {j, i + 1, zlen}]
+
Sum[
2^(-wpow[[i]] - wpow[[j]])*c2m[[n, wpow[[i]] + 1]]*
c2m[[n, wpow[[j]] + 1]]*
Binomial[wpow[[i]] + wpow[[j]], wpow[[i]]], {i, 1,
wlen - 1}, {j, i + 1, wlen}]

]


which gives me some warnings:

Compile::cpintlt: i+1 at position 2 of zpow[[i+1]] should be either a nonzero integer or a vector of nonzero integers; evaluation will use the uncompiled function. >>

Compile::cpintlt: i+1 at position 2 of zpow[[i+1]] should be either a nonzero integer or a vector of nonzero integers; evaluation will use the uncompiled function. >>

Compile::cpintlt: i+1 at position 2 of wpow[[i+1]] should be either a nonzero integer or a vector of nonzero integers; evaluation will use the uncompiled function. >>

General::stop: Further output of Compile::cpintlt will be suppressed during this calculation. >>

CompilePrint shows the following weird snippet:

11  R2I4T(I1)3I8T(I1)4 = MainEvaluate[                                 \
\
\
-zpow[] - zpow[[i + 1]]
Function[{Nn, M, n, basisvecs, c1m, c2m, zlenCompile$1, \ wpowCompile$2, wlenCompile$3, zpowCompile$4}, Block[{zlen = \
zlenCompile$1, wpow = wpowCompile$2, wlen = wlenCompile$3, zpow = \ zpowCompile$4}, {Length[{2                           c1m[[n,zpow[] \
+ 1]] c1m[[n,zpow[[i + 1]] + 1]] Binomial[zpow[] + zpow[[i + 1]], \
zpow[]]}[]], zlen, wpow, wlen, zpow}]][ I0, I1, I2, T(I2)0, \
T(I2)1, T(I2)2, I4, T(I1)3, I8, T(I1)4]]


Could I get some help in understanding this? Seems like i is not assigned a value yet when a term in the sum is calculated, even though j starts at i+1? How can I remedy this?

• It appears that Binomial is not on the List of Compilable Functions -- You could try Gamma and perhaps Round. – Michael E2 May 3 '15 at 23:26
• @MichaelE2, thanks, after I didn't see it in the list I just implemented a compiled version myself (though not relevant to the current issue), but I overlooked Gamma in the list added for 10.0.2. I'll compare! – Marius Ladegård Meyer May 4 '15 at 4:04

It happens quite often that compilation of a correctly running Mathematica function for some or other reason fails. The reason for that is not always easy to trace. After some testing I feel that the problem you ran into is close to a bug.

First let us further simplify your problem, showing the same warnings:

fc=Compile[{}, Block[{mat},
mat=Array[1&, {10,10}];
Sum[mat[[i,j]], {i,1,10},{j,i,10}]]];

(*
During evaluation of In:= Compile::cpintlt: i at position 3 of mat[[1,i]] should be either a nonzero integer or a vector of nonzero integers; evaluation will use the uncompiled function. >>
During evaluation of In:= Compile::cpintlt: i at position 3 of mat[[1,i]] should be either a nonzero integer or a vector of nonzero integers; evaluation will use the uncompiled function. >>


*)

These two warnings are strange: there is no position 3 in mat[[1,i]], the variable i indeed is an integer, and there is no expression mat[[1, i]] in the code.

The warnings have nothing to do with the double summation. When we replace mat[[i,j]] by 1 or i+j, we get no warnings and there is no call of MainEvaluate.

The warnings only turn up when the second underlimit depends on the first index. For example,

fc=Compile[{}, Block[{mat},
mat=Array[1&, {10,10}];
Sum[mat[[i,j]], {i,1,10},{j,1,i}]]];


compiles nicely.

Somehow, your problem is a combination of using a Part-function in the first argument of your double sum and a dependency of the underlimit of the second index of the first index. In all other cases it works, and I think it should work in this case as well. So very likely a bug.

When we replace the double summation by a nested summation, there are no warnings:

fc=Compile[{}, Block[{mat},
mat=Array[1&, {10,10}];
Sum[Sum[mat[[i,j]], {j,i,10}], {i,1,10}]]];


But CompilePrint shows a call of the kernel for evaluating the nested sums, and moreover the inner sum is not recognized as returning an integer; the result is now a real:

fc[]
(* 55. *)


So this does not help.

When you want to use the double summation, the only workaound I see is to formulate it in such a way that the second underlimit does not depend on the first index:

fc=Compile[{}, Block[{mat},
mat=Array[1&, {10,10}];
Sum[mat[[i,j]], {i,1,10}, {j,10,i,-1}]]];


Ridiculous, but it works. CompilePrint does not show a MainEvaluate.

fc[]
(* 55 *)

• Thanks for the nice analysis. Do you expect the "ridiculous solution" to perform any better than a res=0;Do[res+=...] type of formulation? I would guess it's about the same. – Marius Ladegård Meyer May 3 '15 at 19:05
• No, I don't think so. My experience is that the more primitive the code, the faster the compiled function. So it might even be that your Do-construct is faster (though I expect there will be no difference), – Fred Simons May 3 '15 at 19:12