Why doesn't my code compile?

Question

Clear["`*"];
cf = Compile[{{A, _Real, 2}, {n, _Integer}},
   With[{a = A[[1]], b = A[[2]], c = A[[3]]},
    If[n < 1, {{a, b, c}},
      Join[cf[{a, (a + b)/2, (a + c)/2}, n - 1],
       cf[{(b + a)/2, b, (b + c)/2}, n - 1],
       cf[{(c + a)/2, (c + b)/2, c}, n - 1]
       ]
     ];
    ]
   ];

cf[{{0, 0}, {1, 0}, {.5, .8}} // N, 3]

This code I wrote can't be compiled, Mathematica returns

CompiledFunction::cfex: Could not complete external evaluation at instruction 19; proceeding with uncompiled evaluation. >>

How can I correctly compile it?

The uncompiled version:

Clear["`*"];
f[{a_, b_, c_}, n_] :=
  If[n < 1, {{a, b, c}},
   Join[f[{a, (a + c)/2, (a + b)/2}, n - 1],
    f[{(b + a)/2, b, (b + c)/2}, n - 1],
    f[{(c + a)/2, (c + b)/2, c}, n - 1]]
   ];

n = 3;
data = f[{{0, 0}, {1, 0}, {.5, .8}} // N, n];
ListLinePlot[data /. {a_, b_, c_} :> {a, b, c, a}]

Answer 1

Here is a compiled version. First note, however, that your compiled version and the top level version are not the same.

ClearAll[cf]
cf = Compile[{{A, _Real, 2}, {n, _Integer}}, 
   Block[{a = A[[1]], b = A[[2]], c = A[[3]]},
    If[n < 1, A,
     Join[cf[{a, (a + c)/2, (a + b)/2}, n - 1], 
      cf[{(b + a)/2, b, (b + c)/2}, n - 1], 
      cf[{(c + a)/2, (c + b)/2, c}, n - 1]]
     ]
    ]
   ];

This is called via (see further down):

dc = Partition[cf[input, n], 3];

A slightly different version that does not need an additional Partition

ClearAll[cf2]
cf2 = Compile[{{A, _Real, 2}, {n, _Integer}}, 
   Block[{a = A[[1]], b = A[[2]], c = A[[3]]},
    If[n < 1, {A},
     Join[cf2[{a, (a + c)/2, (a + b)/2}, n - 1], 
      cf2[{(b + a)/2, b, (b + c)/2}, n - 1], 
      cf2[{(c + a)/2, (c + b)/2, c}, n - 1]]
     ]
    ]
   , {{cf2[_, _], _Real, 3}}
   ];

For a comparison:

n = 10;
input = Developer`ToPackedArray[N[{{0, 0}, {1, 0}, {.5, .8}}]];
AbsoluteTiming[dc = Partition[cf[input, n], 3];]
AbsoluteTiming[dc2 = cf2[input, n];]
AbsoluteTiming[data = f[input, n];]
data === dc === dc2

(*
0.18
0.18
0.71
True
*)

Also, a better way to rearrange the data is:

data[[All, {1, 2, 3, 1}]]

which does not unpack.

ListLinePlot[data[[All, {1, 2, 3, 1}]]]

If you look at the compiled code with CompilePrint you will note that there are calls to CopyTensor - with some clever computed array (by direct assignment of the recursive cf calls into an array in stead of using Join) it may be possible to avoid those and save computational expense.

Hope this helps.