Can I return lists with different dimensions from a compiled function?

Question

Sometimes I'll need several lists with different dimensions to be returned by a compiled function, but compiled function will fail if an irregular list contained:

(* The real case is of course much more complex. *)
a1 = ConstantArray[0., {10}];
a2 = ConstantArray[0., {9}];
Compile[{}, 
  (* for the real case, a lot of low level calculations will be done on a1 and a2 here *)
           {a1, a2}][]

CompiledFunction::cflist: Nontensor object generated; proceeding with uncompiled evaluation. >>

A possible (and maybe the simplest) work around is to add some extra elements to the lists and delete them(or just ignore them) later:

a1 = ConstantArray[0., {10}];
a2 = ConstantArray[0., {10}];
Delete[Compile[{}, {a1, a2}][], {2, -1}]

But it's a waste of memory (especially when lists are big and have higher dimensions).

So I wonder if other solution exists? Of course that solution should work on high dimensional lists like:

a1 = ConstantArray[0., {300, 299, 300}];
a2 = ConstantArray[0., {299, 300, 299}];

---

To avoid possible confusion I'd like to add a more specific example:

ie = 200;
ez = ConstantArray[0., {ie + 1}];
hy = ConstantArray[0., {ie}];

(* Notice the following function hasn't been fixed yet *)
fdtd1d = Compile[{{steps}}, 
   Module[{ie = ie, ez = ez, hy = hy}, 
    Do[
     ez[[2 ;; ie]] = ez[[2 ;; ie]] + (hy[[2 ;; ie]] - hy[[1 ;; ie - 1]]); 
     ez[[1]] = Sin[n/10]; 
     hy[[1 ;; ie]] = hy[[1 ;; ie]] + (ez[[2 ;; ie + 1]] - ez[[1 ;; ie]]), 
      {n, steps}]; 
    {ez, hy}]];

fdtd1d[1000]

---

Er… maybe I should add one more multidimensional sample?:

ie = 100; je = 100;
c = 1/Sqrt[2];
ez = ConstantArray[0., {ie + 1, je + 1}];
hx = ConstantArray[0., {ie + 1, je}];
hy = ConstantArray[0., {ie, je + 1}];

(* Notice the following function hasn't been fixed yet *)
fdtd2d = Compile[{{steps}}, 
   Module[{ez = ez, hx = hx, hy = hy, ie = ie, je = je, c = c},
    Do[
     ez[[2 ;; ie, 2 ;; je]] += 
      c (hx[[2 ;; ie, 1 ;; je - 1]] - hx[[2 ;; ie, 2 ;; je]] + 
         hy[[2 ;; ie, 2 ;; je]] - hy[[1 ;; ie - 1, 2 ;; je]]);
     ez[[Floor[ie/2], Floor[je/2]]] = Sin[n/10];
     hx += c (ez[[All, 1 ;; je]] - ez[[All, 2 ;; je + 1]]);
     hy += c (ez[[2 ;; ie + 1, All]] - ez[[1 ;; ie, All]]),
     {n, steps}]; 
   {ez, hx, hy}]];

AbsoluteTiming[dat = fdtd2d[1000];]
ArrayPlot /@ {ez2, hx2, hy2}

Answer 1

I suppose one could return values through an external call:

ClearAll["f`*", f];
f`compiled = Compile[{{dim1, _Integer}, {dim2, _Integer}},
   Module[{a1, a2},
    a1 = Range[dim1];
    a2 = Range[dim2];
    f`ret1 = a1;
    f`ret2 = a2;
    {Total@a1, Total@a2}
    ]];
f[dim1_, dim2_] := Block[{f`ret1, f`ret2, f`ans},
   f`ans = f`compiled[dim1, dim2];
   {f`ret1, f`ret2, f`ans}
   ];

The data is passed back through via calls to MainEvaluate. There is a cost to that, but perhaps in your actual use case, the cost is not too much.

Example

f[5, 3]

(* {{1, 2, 3, 4, 5}, {1, 2, 3}, {6, 15}} *)

Answer 2

2D update

Here is a way to make your example work with minimal overhead by joining the ragged output into a single larger list and storing the dimensions at the start and putting it back together with a post-processing function.

ie = 100; je = 100;
c = 1/Sqrt[2];
ez = ConstantArray[0., {ie + 1, je + 1}];
hx = ConstantArray[0., {ie + 1, je}];
hy = ConstantArray[0., {ie, je + 1}];

fdtd2d = Compile[{{steps, _Integer, 0}, {ie, _Integer, 
     0}, {je, _Integer, 0}, {c, _Real, 0}, {ez2, _Real, 
     2}, {hx2, _Real, 2}, {hy2, _Real, 2}}, 
   Module[{ez = ez2, hx = hx2, hy = hy2, ezdims, hxdims, hydims}, 
    Do[ez[[2 ;; ie, 2 ;; je]] += 
      c (hx[[2 ;; ie, 1 ;; je - 1]] - hx[[2 ;; ie, 2 ;; je]] + 
         hy[[2 ;; ie, 2 ;; je]] - hy[[1 ;; ie - 1, 2 ;; je]]);
     ez[[Floor[ie/2], Floor[je/2]]] = Sin[n/10];
     hx += c (ez[[All, 1 ;; je]] - ez[[All, 2 ;; je + 1]]);
     hy += c (ez[[2 ;; ie + 1, All]] - ez[[1 ;; ie, All]]), {n, 
      steps}];
    ezdims = Dimensions@ez;
    hxdims = Dimensions@hx;
    hydims = Dimensions@hy;
    Join @@ {{3, Length@ezdims, Length@hxdims, Length@hydims}, ezdims, hxdims, hydims, Flatten@ez, Flatten@hx, Flatten@hy}]];

positionindexer[data_] := 
 With[{accumulated = Accumulate@Join[{0}, data]}, 
  Table[{accumulated[[i]] + 1, accumulated[[i + 1]]}, {i, 
    Length@accumulated - 1}]]

postprocess[data_] := Module[{},
  output = data;
  numberofArrays = Round@First@output;
  output = Rest@output;
  numdims = Round@Take[output, numberofArrays];
  output = Drop[output, Length@numdims];
  dimpositions = positionindexer[numdims];
  dims = Round[Take[output, #] & /@ dimpositions];
  totallength = Times @@@ dims;
  output = Drop[output, Total@numdims];
  datapositions = positionindexer[totallength];
  dataflat = Take[output, #] & /@ datapositions;
  MapThread[ArrayReshape, {dataflat, dims}]]

postprocess@fdtd2d[1000, ie, je, c, ez, hx, hy]

Answer 3

As far as I know, the only case where this is possible is for compiled functions with RuntimeAttributes -> Listable. Example:

testComp  = 
   Compile[{{lst, _Integer, 1}}, lst^2, RuntimeAttributes -> Listable]

so that, for example:

testComp[Range /@ Range[4]]

(* {{1}, {1, 4}, {1, 4, 9}, {1, 4, 9, 16}} *)

I am not aware of other situations where CompiledFunction may return ragged lists.