# How to compile or speed up a function with eigensystem?

I have defined two functions as following:

MapQ = Compile[{{r0, _Real, 2}, {r, _Real, 2}, xc, yc, {w, _Real, 1}},
Block[{ M = ({
{0., 0., 0., 0.},
{0., 0., 0., 0.},
{0., 0., 0., 0.},
{0., 0., 0., 0.}
}), x0, y0, z0, x, y, z, w1, w2, r2, c1, c2, c3, c4, x0x, y0y,
z0z, sys, \[Theta], temp},
Do[
x0 = r0[[i, 1]] - xc; x = r[[i, 1]] - xc;
y0 = r0[[i, 2]] - yc; y = r[[i, 2]] - yc;
z0 = r0[[i, 3]]; z = r[[i, 3]];

(*// matrix elements//*)
w1 = w[[i]];
w2 = 2*w1;
r2 = x^2 + y^2 + z^2 + x0^2 + y0^2 + z0^2;
x0x = x0 x;
y0y = y0 y;
z0z = z0 z;
c1 = (r2 - 2 (x0x + y0y + z0z));
c2 = (r2 + 2 (-x0x + y0y + z0z));
c3 = (r2 + 2 (x0x - y0y + z0z));
c4 = (r2 + 2 (x0x + y0y - z0z));
M[[1, 1]] += w1 c1;
M[[1, 2]] += w2 (y z0 - z y0);
M[[1, 3]] += w2 (-x z0 + z x0);
M[[1, 4]] += w2 (x y0 - y x0);
M[[2, 2]] += w1 c2;
M[[2, 3]] -= w2 (x y0 + y x0);
M[[2, 4]] -= w2 (x z0 + z x0);
M[[3, 3]] += w1 c3;
M[[3, 4]] -= w2 (y z0 + z y0);
M[[4, 4]] += w1 c4;

, {i, Length[r0]}];
M = ({
{M[[1, 1]], M[[1, 2]], M[[1, 3]], M[[1, 4]]},
{M[[1, 2]], M[[2, 2]], M[[2, 3]], M[[2, 4]]},
{M[[1, 3]], M[[2, 3]], M[[3, 3]], M[[3, 4]]},
{M[[1, 4]], M[[2, 4]], M[[3, 4]], M[[4, 4]]}
})

], RuntimeAttributes -> {Listable}, Parallelization -> True,
CompilationTarget -> "C"];

AngleQ[M_] := Block[{sys, \[Theta]},
sys = Eigensystem[M];
If[sys[[2, 4, 1]] < 0, sys[[2, 4]] = -sys[[2, 4]], None];
sys[[2, 4]] = Normalize[sys[[2, 4]]];
\[Theta] = 2 ArcCos[sys[[2, 4, 1]]];
If[sys[[2, 4, 4]]/Sin[\[Theta]/2] < 0, \[Theta] = -\[Theta], None];
\[Theta]
]


The first function can be compiled to generate C-code and parallelizable.

My question is that how to combine these two functions together and speed up eventually?

Since Compile does not allow me to put Eigensystem inside to optimize, any suggestion?

Many thanks!

Update

In particular, I want to map this to a set of inputs. As in the following: My question is how can I parallelize this map? I try ParallelMap but it seems not work.

r.txt

r0.txt

w.txt

xc1:24.2

yc1:31.37

• Btw.: Questions about code without example data or a usage example are hard to answer. Please keep that in mind when you post here. Commented Jul 22, 2019 at 14:26

Here is my attempt on the function MapQ.

# The modified code

MapQ2 = Compile[{{r0, _Real, 2}, {r, _Real, 2}, xc, yc, {w, _Real, 1}},
Block[{M11, M12, M13, M14, M22, M23, M24, M33, M34, M44, x0, y0, z0, x, y, z, w1, w2, r2, c1, c2, c3, c4, x0x, y0y, z0z},
M11 = M12 = M13 = M14 = M22 = M23 = M24 = M33 = M34 = M44 = 0.;
Do[
x0 = CompileGetElement[r0, i, 1] - xc;
x = CompileGetElement[r, i, 1] - xc;
y0 = CompileGetElement[r0, i, 2] - yc;
y = CompileGetElement[r, i, 2] - yc;
z0 = CompileGetElement[r0, i, 3];
z = CompileGetElement[r, i, 3];
(*//matrix elements//*)

w1 = CompileGetElement[w, i];
w2 = 2. w1;
r2 = x^2 + y^2 + z^2 + x0^2 + y0^2 + z0^2;
x0x = x0 x;
y0y = y0 y;
z0z = z0 z;
c1 = (r2 - 2. (x0x + y0y + z0z));
c2 = (r2 + 2. (-x0x + y0y + z0z));
c3 = (r2 + 2. (x0x - y0y + z0z));
c4 = (r2 + 2. (x0x + y0y - z0z));
M11 += w1 c1;
M12 += w2 (y z0 - z y0);
M13 += w2 (-x z0 + z x0);
M14 += w2 (x y0 - y x0);
M22 += w1 c2;
M23 -= w2 (x y0 + y x0);
M24 -= w2 (x z0 + z x0);
M33 += w1 c3;
M34 -= w2 (y z0 + z y0);
M44 += w1 c4;
, {i, Length[r0]}];
{
{M11, M12, M13, M14},
{M12, M22, M23, M24},
{M13, M23, M33, M34},
{M14, M24, M34, M44}
}
],
CompilationTarget -> "C",
Parallelization -> True,
RuntimeAttributes -> {Listable},
RuntimeOptions -> "Speed"
];


# Timing example and test

n = 100000;
r0 = RandomReal[{-1, 1}, {n, 3}];
r = RandomReal[{-1, 1}, {n, 3}];
w = RandomReal[{-1, 1}, {n}];

aa = MapQ[r, r0, 0., 0., w]; // RepeatedTiming // First
bb = MapQ2[r, r0, 0., 0., w]; // RepeatedTiming // First
Max[Abs[aa - bb]]


0.014

0.00092

0.

So this generates exactly the same result, but 15 times faster.

# Crucial changes

• Storing the entries of the upper triangle matrix of M in scalar registers; these are faster to access than an array, because every access operation )read or write) into an array with Part in Mathematica will put a bound-check into the C code. For read operations, this can be mended by CompileGetElement, but I do not know of a way to do the same for write operations. And the frequently used += and -= are both read and write operations. This is why using real registers during the computations and building the matrix only in the end is so much faster.

• As noted above, I replaced a few occurences of Part by CompileGetElement. In order to ake this effective, one has to add the option RuntimeOptions -> "Speed"; only then the creation of C-code for the bound checking will be deactivated.

• I replaced a few occurences of integer 2 by 2.; otherwise one has a few needless runtime type-casts. Does not make a big difference in this case, though.

As it turns out, the Eigensystem is not the bottleneck here. Eigensystem effectlively calls a compiled library and cannot be sped up with Compile any further. Which does not meen that it were impossible to write a more efficient function that is dedicated to compute solely 4D eigensystems. In principle, that might be possible because polynomials of order 4 can be solved algebraically in closed form. But I would try that only if Eigensystem consitutes more than, say, 30% of your final running time. And even then I would have to think about how worthwhile this investment of programming time would be.

• Thanks! It works. Now I want to proceed further. To parallelize a lot of similar simulation. Could you see the original post? Commented Jul 22, 2019 at 16:48
• Could you please post example input data? Commented Jul 22, 2019 at 16:57
• It would be interesting to have some kind of simple Compile optimizer that does things like register unwrapping whenever it sees that a given matrix variable is only ever accessed by explicit indices. Commented Jul 22, 2019 at 17:10
• @b3m2a1 True. I am also not convinced that the unwrapping of loop constructs done by Mathematica in the creation of C code is for the better. Modern C compilers should be able to vectorize for-loops, but I do not know if they can optimize these goto` mess produced by Mathematica... oO Commented Jul 22, 2019 at 17:15
• @HenrikSchumacher if this were not WRI we're talking about, I'd imagine the new compiler could be able to make better use of modern compiler tech once it matures more. Since this is WRI...maybe there's like a 50% chance this happens. Commented Jul 22, 2019 at 17:18