# My ParallelDo does not work

I have to do a lot of calculations that take lot of time, and using a Do loop is simply too long. It is the first time I am using ParallelDo and it is not working as fast as it should. I am sure it is an easy fix, likely due to how I define the variables (shared, not shared... I cannot understand how to do it right).

My code is the following:

(*Parameter to set the size of the computation. in the final version will be 15*)
Dim = 8;
(*Matrices considered*)
eHe = RandomReal[{-1, 1}, {2^Dim, 2^Dim}];
eHe = eHe + ConjugateTranspose[eHe];
(*Some parameters*)
Nst = Log2[Length[eHe]]
SysDim = N[2^Nst];

(*Matrices for the scalar product*)
PauliString = {"Id", "X", "Y", "Z"};
S0 = SparseArray[{{1, 1} -> N[1], {2, 2} -> N[1]}];
S1 = SparseArray[{{1, 2} -> N[1], {2, 1} -> N[1]}];
S2 = SparseArray[{{1, 2} -> N[-I], {2, 1} -> N[I]}];
S3 = SparseArray[{{1, 1} -> N[1], {2, 2} -> N[-1]}];
SVec = {S0, S1, S2, S3};

(*Parameters for the ParallelDo*)
Soglia = 10^-10;
Hper = eHe;
resHper = {};
SetSharedVariable[resHper];

(*ParallelDo*)
ParallelDo[

tmp = Tr[
ConjugateTranspose[
Apply[KroneckerProduct,
Table[SVec[[el[[j]]]], {j, 1, Nst}]]].Hper]/SysDim;

If[Abs[tmp] > Soglia,
AppendTo[resHper,
Flatten[
Append[{tmp},
Table[PauliString[[el[[j]]]], {j, 1, Nst}]
]
]
]
]
, {el, Tuples[{1, 2, 3, 4}, Nst]}]



Here, Hper and Apply[KroneckerProduct,Table[Vec[[el[[j]]]],{j, 1, Nst}]]] are matrices and Table[PauliString[[el[[j]]]] a table of strings. The ParallelDo loop actually works, but is seemingly much slower if compared to a normal Do loop...

I have spent really long time, and have not understood where the problem is. Any help would be really appreciated. Thanks!

UPDATE: The code above now seems to yield the correct result, but is much slower than expected. I have benchmarked the code by substituting ParallelDo with Do and the first takes up to 18 times longer if compared to Do. I have checked the CPU and it is always mostly unused when using the ParallelDo. I guess it is because Mathematica has trouble with the tmp variable inside the ParallelDo?

• tmp is a shared variable but you’re setting it in every kernel, so there’s a big problem. Dec 6, 2022 at 18:22
• Yes I fully agree with you... The point is that I do not know how to solve this issue, and wasn't able to find any tutorial. Any guess? Dec 6, 2022 at 18:26
• Your code does not run. Please make it self-contained, so that people can experiment and help. Dec 6, 2022 at 18:52
• Might be a question of coarse vs fine graining. If ParallelDo is using a bad choice by default, possibly an explicit setting will help. Dec 6, 2022 at 21:24
• See this: mathematica.stackexchange.com/a/138911/12 and also this: mathematica.stackexchange.com/q/48295/12 The general advice is to use data parallelism. If you cannot formulate your problem in terms of ParalellTable or similar (more generally: ParallelCombine), it is likely just not a good fit for parallelization. Dec 7, 2022 at 9:49

Part 1. Compare

Tr1[A_,B_]:=Tr[ConjugateTranspose[A].B];
Tr2[A_,B_]:=With[{X=Most[ArrayRules[A]]},Conjugate[X[[;;,2]]].Extract[B,X[[;;,1]]]];


They give the same results but in OP's situation and at least with Dim=12 the second is much faster:

Dim=12;
SeedRandom[1];
{Range[2^Dim],Permute[Range[2^Dim],RandomPermutation[2^Dim]]}]];
B=RandomReal[{-1,1},{2^Dim,2^Dim}];

AbsoluteTiming[Tr1[A,B]]
(* {0.110279, -4.45864+3.28819 I} *)

AbsoluteTiming[Tr2[A,B]]
(* {0.002091, -4.45864+3.28819 I} *)


For Dim=8 there is little difference, but OP mentions that they are interested in larger Dim.

Part 2. Here is a modification of OP's code.

• Instead of shared variables and ParallelDo, it uses ParallelTable.
• To collect results locally, it uses Reap-Sow and avoids AppendTo.
• Coarse graining.
• Some other changes for convenience.

Please restart the kernel before trying this:

(*Parameter to set the size of the computation.in the final version will be 15*)
Dim=8;

(*Matrices considered*)
eHe=RandomReal[{-1,1},{2^Dim,2^Dim}];
eHe=eHe+ConjugateTranspose[eHe];

(*Some parameters*)
Nst=Log2[Length[eHe]]
SysDim=N[2^Nst];

(*Matrices for the scalar product*)
PS[1]="Id"; PS[2]="X"; PS[3]="Y"; PS[4]="Z";
S[1]=SparseArray[{{1,1}->N[1],{2,2}->N[1]}];
S[2]=SparseArray[{{1,2}->N[1],{2,1}->N[1]}];
S[3]=SparseArray[{{1,2}->N[-I],{2,1}->N[I]}];
S[4]=SparseArray[{{1,1}->N[1],{2,2}->N[-1]}];

(*Parameters for the ParallelDo*)
Soglia=10^-10;
Hper=eHe;

Tr2[A_,B_]:=With[{X=Most[ArrayRules[A]]},Conjugate[X[[;;,2]]].Extract[B,X[[;;,1]]]];

NstX=3; (* should be ok if num of kernels is a divisor of 4^NstX *)

(*ParallelTable*)
resHper=Join@@ParallelTable[Reap[
Do[With[{el=Join[elX,elY]},
With[{tmp=Tr2[Apply[KroneckerProduct,Map[S,el]],Hper]/SysDim},
If[Abs[tmp]>Soglia,Sow[Prepend[Map[PS,el],tmp]]]]],
{elY,Tuples[{1,2,3,4},Nst-NstX]}]][[2,1]],
{elX,Tuples[{1,2,3,4},NstX]}];


The result in resHper should agree with OP's result up to reordering. It does give a considerable parallel speed-up when compared to Table.

(Final comment: I do not think this is optimal. Doing the calculation separately for every el is not optimal. In my code for example, one could contract Hper with just Map[S,elX] first, and only contract with Map[S,elY] in an inner loop, or something like this. Will stop here.)

• This is a massive improvement! As a small comment, it seems it is actually advantageous if matrix A is much sparser than matrix B... Works for me, so thanks a lot! Dec 6, 2022 at 21:07
• Fantastic, this works perfectly! I was ready to give up, due to a comment above, but your code works fantastically :)! Dec 8, 2022 at 14:30