5
$\begingroup$

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?

$\endgroup$
7
$\begingroup$

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you for showcasing Compile``GetElement as well. Can you please explain why I need to use two localisation constructs to create code instead of writing directly f@Table[Compile``GetElement[X, i, j], {i, 1, 3}, {j, 1, 3}]? $\endgroup$ – Marco Apr 1 at 14:12
  • $\begingroup$ You're welcome. The Block shields X and XX from potential global definitions and the With is a trick to enforce that the symbolic expression is evaluated and supplied in verbatim to Compile. One could also try Compile[{...}, Evaluate[f@XX],...], but that is not very robust. $\endgroup$ – Henrik Schumacher Apr 1 at 14:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.