Compiling with Listable option and Parallelization

Question

I have a bunch of triangles (200k) I need to calculate the normals of. Here are five of them below, each one with its vertices in CCW order.

tris = {{{99.1175,-156.51,158.},{-63.7411,-173.073,158.},{-62.,-173.,158.}},{{-62.,-173.,3.},{-62.,-173.,158.},{-63.7411,-173.073,158.}},{{62.,-173.,158.},{99.1175,-156.51,158.},{-62.,-173.,158.}},{{-62.,-173.,3.},{62.,-173.,158.},{-62.,-173.,158.}},{{-99.1175,-156.51,158.},{-65.4575,-173.287,158.},{-63.7411,-173.073,158.}}};

Their normals can be simply calculated with

Normalize@Cross[#[[2]] - #[[1]], #[[3]] - #[[1]]] & /@ tris;

I'd like to Compile such function, therefore I reworked everything in the form

func=Compile[
{ax,ay,az,bx,by,bz,cx,cy,cz},
{
(-az*by+ay*bz+az*cy-bz*cy-ay*cz+by*cz)/Sqrt[Abs[-ay*bx+ax*by+ay*cx-by*cx-ax*cy+bx*cy]^2+Abs[az*bx-ax*bz-az*cx+bz*cx+ax*cz-bx*cz]^2+Abs[-az*by+ay*bz+az*cy-bz*cy-ay*cz+by*cz]^2],
(az*bx-ax*bz-az*cx+bz*cx+ax*cz-bx*cz)/Sqrt[Abs[-ay*bx+ax*by+ay*cx-by*cx-ax*cy+bx*cy]^2+Abs[az*bx-ax*bz-az*cx+bz*cx+ax*cz-bx*cz]^2+Abs[-az*by+ay*bz+az*cy-bz*cy-ay*cz+by*cz]^2],
(-ay*bx+ax*by+ay*cx-by*cx-ax*cy+bx*cy)/Sqrt[Abs[-ay*bx+ax*by+ay*cx-by*cx-ax*cy+bx*cy]^2+Abs[az*bx-ax*bz-az*cx+bz*cx+ax*cz-bx*cz]^2+Abs[-az*by+ay*bz+az*cy-bz*cy-ay*cz+by*cz]^2]
}
];

and run

func[Sequence @@ Flatten @ #] & /@ tris;

The speed gain is considerable (9.9 vs 1.2 seconds!), but I'd like to go further, namely using the Listable property and the Parallelization option. How can I do that? If I understand correctly parallelizability is a consequence of listability, so I added RuntimeAttributes -> {Listable}, Parallelization -> True at the end of func. But, how should I call the function then?

Answer 1

Easy: Make the input a matrix (array of rank 2); then the compiled function can thread directly over the input. No need for Map anymore. Also, you can use Mathematica to generate the symbolic code automatically from your function.

f = Normalize@Cross[#[[2]] - #[[1]], #[[3]] - #[[1]]] &;
Block[{X, XX}, 
XX = Table[Compile`GetElement[X, i, j], {i, 1, 3}, {j, 1, 3}];
cf = With[{code = f@XX},
  Compile[{{X, _Real, 2}},
    code,
    CompilationTarget -> "C",
    RuntimeAttributes -> {Listable},
    Parallelization -> True,
    RuntimeOptions -> "Speed"
    ]
   ];
 ]

Now, on my Haswell Quad Core:

tris = RandomReal[{-1, 1}, {100000, 3, 3}];
r1 = f /@ tris; // AbsoluteTiming // First
r2 = cf[tris]; // AbsoluteTiming // First
Max[Abs[r1 - r2]]

3.18289

0.006599

2.39808*10^-14

The undocumented Compile`GetElement is a read-only instruction that---in combination with RuntimeOptions -> "Speed"--- does not perform bound checks. So make sure to call this functions with matrices at least of size $3 \times 3$; otherwise, the library tries to read memory that it is not allowed to read and the operating system kills the process (together the Mathematica kernel as its parent). So when you write a program that is supposed to be used by other people, always write a Mathematica wrapper that performs all the checks before calling the CompiledFunction.