Speed up dot and cross products for a list of vectors, in a "brick wall" pattern

Question

Yesterday I asked a question on how to speed of computation of dot and cross products for a list of vectors, for which I received a nice answer from Henrik Schumacher.

Now I have a follow up question to that. In the question, in one time step one had to perform the operation,

$f(\mathbf{\tilde{S}}_{2n-1},\mathbf{\tilde{S}}_{2n})$.

Instead if I want to perform two operations $f(\mathbf{\tilde{S}}_{2n-1},\mathbf{\tilde{S}}_{2n})$ and $f(\mathbf{\tilde{S}}_{2n},\mathbf{\tilde{S}}_{2n+1})$, i.e. the "brick wall" implementation on one time step, how should the code be implemented?

My understanding is I have to take the loop out of the compiled function in this case, and here is my implementation,

randomSpherePoint[n_] := Module[{z, ϕ, r},
   ϕ = RandomReal[{0, 2 Pi}, n];
   z = RandomReal[{-1, 1}, n];
   r = Sqrt[1. - z z];
   Transpose[{r Cos[ϕ], r Sin[ϕ], z}]
   ];

cf = Compile[{{x, _Real, 1}, {y, _Real, 1}, {a, _Real}, {b, _Real}, {c, _Real}, {iters, _Integer}},
   Block[{x1, x2, x3, y1, y2, y3, u1, u2, u3, v1, v2, v3, ufactor, 
     vfactor,dot1},
    
    x1 = Compile`GetElement[x, 1];
    x2 = Compile`GetElement[x, 2];
    x3 = Compile`GetElement[x, 3];
    
    y1 = Compile`GetElement[y, 1];
    y2 = Compile`GetElement[y, 2];
    y3 = Compile`GetElement[y, 3];
    
    u1 = u2 = u3 = v1 = v2 = v3 = 0.;
    
    Table[
    
     u1 = a x1 + b y1 + c (-x3 y2 + x2 y3);
     u2 = a x2 + b y2 + c (x3 y1 - x1 y3);
     u3 = a x3 + c (-x2 y1 + x1 y2) + b y3;
     v1 = b x1 + a y1 - c (-x3 y2 + x2 y3);
     v2 = b x2 + a y2 - c (x3 y1 - x1 y3);
     v3 = b x3 - c (-x2 y1 + x1 y2) + a y3;
     ufactor = 1./Sqrt[u1 u1 + u2 u2 + u3 u3];
     vfactor = 1./Sqrt[v1 v1 + v2 v2 + v3 v3];
     
     x1 = u1 ufactor;
     x2 = u2 ufactor;
     x3 = u3 ufactor;
     
     y1 = v1 vfactor;
     y2 = v2 vfactor;
     y3 = v3 vfactor;
     
     {{x1, x2, x3}, {y1, y2, y3}}
     
     , {i, 1, iters}]
    ],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];
a = 2.;
b = 3.;
c = 4.;
N1 = 10000;
iters = 20;

x = randomSpherePoint[N1/2];
y = randomSpherePoint[N1/2];

spins = Riffle[x, y];
result1 = Table[
     list1 = spins[[1 ;; N1 ;; 2]]; list2 = spins[[2 ;; N1 ;; 2]]; 
     list3 = MapThread[Dot, {list1, list2}]; 
     list31 = MapThread[Cross, {list1, list2}]; 
     list32 = Sqrt[a^2 + b^2 + c^2 (1 - list3^2) + 2 a b list3]; 
     list4 = (a list1 + b list2 + c list31)/(list32); 
     list5 = (a list2 + b list1 - c list31)/(list32); 
     spins = Normalize /@ Flatten[Transpose[Join[{list4, list5}]], 1];
     
     list1 = spins[[2*Range[N1/2 - 1]]]; list2 = spins[[2*Range[N1/2 - 1] + 1]]; list31 = MapThread[Cross, {list1, list2}]; list3 = MapThread[Dot, {list1, list2}]; 
     list32 = Sqrt[a^2 + b^2 + c^2 (1 - list3^2) + 2 a b list3]; 
     list4 = (a list1 + b list2 + c list31)/(list32); 
     list5 = (a list2 + b list1 - c list31)/(list32);  
     spins = Normalize /@ Join[{spins[[1]]}, Flatten[Transpose[Join[{list4, list5}]], 1], {spins[[N1]]}]
     
     
     
     , {i, 1, iters}]; // AbsoluteTiming // First

(*result2 = cf[x, y, a, b, c, iters]; // AbsoluteTiming // First*)


rs2=Table[{result2 = cf[x, y, a,b, c, 1];
x1=result2[[All,1,1]];
y1=result2[[All,1,2]];
xtemp=x1[[1]];
ytemp=y1[[N1/2]];
x2=y1[[1;;N1/2-1]];
y2=x1[[2;;N1/2]];
result2 = cf[x2, y2, a,b, c, 1];
x=Join[{xtemp},result2[[All,1,2]]],
y=Join[result2[[All,1,1]],{ytemp}]}

,{i,1,iters}];//AbsoluteTiming //First

Max[Abs[result1 - Flatten[rs2, {{1}, {3, 2}, {4}}]]]

The results in my machine are as follows,

5.16973

0.027213

4.12334*10^-11

Both the error and the speed up are orders of magnitude worse. What am I doing wrong? Is there any way to compile the entire iteration process as was the case in the last question?

Answer 1

The main problem is that you use again Table which is quite a high language construct.

The data dependicies makes it indeed harder to parallelize that. However, a single-threaded compiled version is readily constructed:

cg = Compile[{{spins0, _Real,2}, {a, _Real}, {b, _Real}, {c, _Real}, {iters, _Integer}}, 
  Block[{A, spins, n, x1, x2, x3, y1, y2, y3, z1, z2, z3, u1, u2, u3, v1, v2, v3, w1, w2, w3, ufactor, vfactor, wfactor},
   n = Length[spins0];
   spins = spins0;
   u1 = u2 = u3 = v1 = v2 = v3 = 0.;
   
   Table[
    Do[
     x1 = Compile`GetElement[spins, k, 1];
     x2 = Compile`GetElement[spins, k, 2];
     x3 = Compile`GetElement[spins, k, 3];
     y1 = Compile`GetElement[spins, k + 1, 1];
     y2 = Compile`GetElement[spins, k + 1, 2];
     y3 = Compile`GetElement[spins, k + 1, 3];
     z1 = Compile`GetElement[spins, k + 2, 1];
     z2 = Compile`GetElement[spins, k + 2, 2];
     z3 = Compile`GetElement[spins, k + 2, 3];
     
     u1 = a x1 + b y1 + c (-x3 y2 + x2 y3);
     u2 = a x2 + b y2 + c (x3 y1 - x1 y3);
     u3 = a x3 + c (-x2 y1 + x1 y2) + b y3;
     v1 = b x1 + a y1 - c (-x3 y2 + x2 y3);
     v2 = b x2 + a y2 - c (x3 y1 - x1 y3);
     v3 = b x3 - c (-x2 y1 + x1 y2) + a y3;
     ufactor = 1./Sqrt[u1 u1 + u2 u2 + u3 u3];
     vfactor = 1./Sqrt[v1 v1 + v2 v2 + v3 v3];
     x1 = u1 ufactor;
     x2 = u2 ufactor;
     x3 = u3 ufactor;
     y1 = v1 vfactor;
     y2 = v2 vfactor;
     y3 = v3 vfactor;
     spins[[k, 1]] = x1;
     spins[[k, 2]] = x2;
     spins[[k, 3]] = x3;
     spins[[k + 1, 1]] = y1;
     spins[[k + 1, 2]] = y2;
     spins[[k + 1, 3]] = y3;
     
     , {k, 1, n - 1, 2}];
    
    Do[
     x1 = Compile`GetElement[spins, k, 1];
     x2 = Compile`GetElement[spins, k, 2];
     x3 = Compile`GetElement[spins, k, 3];
     y1 = Compile`GetElement[spins, k + 1, 1];
     y2 = Compile`GetElement[spins, k + 1, 2];
     y3 = Compile`GetElement[spins, k + 1, 3];
     
     u1 = a x1 + b y1 + c (-x3 y2 + x2 y3);
     u2 = a x2 + b y2 + c (x3 y1 - x1 y3);
     u3 = a x3 + c (-x2 y1 + x1 y2) + b y3;
     v1 = b x1 + a y1 - c (-x3 y2 + x2 y3);
     v2 = b x2 + a y2 - c (x3 y1 - x1 y3);
     v3 = b x3 - c (-x2 y1 + x1 y2) + a y3;
     ufactor = 1./Sqrt[u1 u1 + u2 u2 + u3 u3];
     vfactor = 1./Sqrt[v1 v1 + v2 v2 + v3 v3];
     x1 = u1 ufactor;
     x2 = u2 ufactor;
     x3 = u3 ufactor;
     y1 = v1 vfactor;
     y2 = v2 vfactor;
     y3 = v3 vfactor;
     spins[[k, 1]] = x1;
     spins[[k, 2]] = x2;
     spins[[k, 3]] = x3;
     spins[[k + 1, 1]] = y1;
     spins[[k + 1, 2]] = y2;
     spins[[k + 1, 3]] = y3;
     
     , {k, 2, n - 1, 2}];
    
    spins
    
    , {i, 1, iters}]
   ],
  CompilationTarget -> "C",
  RuntimeOptions -> "Speed"
  ]

Now

spins0 = Riffle[x, y]
result2 = cg[spins0, a, b, c, iters]

produces (up to some small error) the same as result1. On my machine it is about 420 times faster. Multiply this with the number of cores (it is 8 on my machine), then you get a speedup in the right ball park.