5
$\begingroup$

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?

$\endgroup$

1 Answer 1

5
$\begingroup$

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.

$\endgroup$
3
  • 1
    $\begingroup$ 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. :) $\endgroup$ Jun 3 at 10:57
  • 2
    $\begingroup$ 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. $\endgroup$ Jun 3 at 11:00
  • 1
    $\begingroup$ Right the ensemble averaging, sheesh brain not working today. $\endgroup$ Jun 3 at 15:11

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.