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

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 = CompileGetElement[x, 1];
x2 = CompileGetElement[x, 2];
x3 = CompileGetElement[x, 3];

y1 = CompileGetElement[y, 1];
y2 = CompileGetElement[y, 2];
y3 = CompileGetElement[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]];
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?

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 = CompileGetElement[spins, k, 1];
x2 = CompileGetElement[spins, k, 2];
x3 = CompileGetElement[spins, k, 3];
y1 = CompileGetElement[spins, k + 1, 1];
y2 = CompileGetElement[spins, k + 1, 2];
y3 = CompileGetElement[spins, k + 1, 3];
z1 = CompileGetElement[spins, k + 2, 1];
z2 = CompileGetElement[spins, k + 2, 2];
z3 = CompileGetElement[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 = CompileGetElement[spins, k, 1];
x2 = CompileGetElement[spins, k, 2];
x3 = CompileGetElement[spins, k, 3];
y1 = CompileGetElement[spins, k + 1, 1];
y2 = CompileGetElement[spins, k + 1, 2];
y3 = CompileGetElement[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.

• Thanks. This style seems pretty similar to how I would code this in FORTRAN. I did get that Table was the bottleneck again but couldn't figure out how to tackle it. Anyway I shall compare speeds in our server for this and the other way around, though probably that is not parallelizing as well. Thanks for the free lesson. :) Commented Jun 3, 2022 at 10:57
• You're welcome. Regarding parallelization: It's typically the best strategy to parallelize the most outer loop. So if you have to run many similar experiments, just run the single threaded code on each core of the server. Commented Jun 3, 2022 at 11:00
• Right the ensemble averaging, sheesh brain not working today. Commented Jun 3, 2022 at 15:11