My ParallelDo does not work

Question

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?

Answer 1

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];
A=SparseArray[MapThread[({#1,#2}->RandomComplex[{-1-I,1+I}])&,
     {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.)