2
$\begingroup$

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]
  ]

enter image description here 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. enter image description here

Input file can be download from here:

r.txt

r0.txt

w.txt

xc1:24.2

yc1:31.37

$\endgroup$
  • 1
    $\begingroup$ Btw.: Questions about code without example data or a usage example are hard to answer. Please keep that in mind when you post here. $\endgroup$ – Henrik Schumacher Jul 22 at 14:26
6
$\begingroup$

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 = Compile`GetElement[r0, i, 1] - xc;
     x = Compile`GetElement[r, i, 1] - xc;
     y0 = Compile`GetElement[r0, i, 2] - yc;
     y = Compile`GetElement[r, i, 2] - yc;
     z0 = Compile`GetElement[r0, i, 3];
     z = Compile`GetElement[r, i, 3];
     (*//matrix elements//*)

     w1 = Compile`GetElement[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 Compile`GetElement, 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 Compile`GetElement. 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.

$\endgroup$
  • $\begingroup$ Thanks! It works. Now I want to proceed further. To parallelize a lot of similar simulation. Could you see the original post? $\endgroup$ – Bob Lin Jul 22 at 16:48
  • $\begingroup$ Could you please post example input data? $\endgroup$ – Henrik Schumacher Jul 22 at 16:57
  • 1
    $\begingroup$ 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. $\endgroup$ – b3m2a1 Jul 22 at 17:10
  • 1
    $\begingroup$ @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 $\endgroup$ – Henrik Schumacher Jul 22 at 17:15
  • 1
    $\begingroup$ @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. $\endgroup$ – b3m2a1 Jul 22 at 17:18

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.