0
$\begingroup$

I do belive that Mathematica can work as fast as other packages, but it is really sensitive to the way we write the code. Here I would like to evaluate a trace resulted from a matrix multiplication at many points as follwos

KS[x_, y_,z_] := {{0, I Sin[x] + Sin[y], 
    13/4 - Cos[x] - Cos[y] - Cos[z] - I Sin[z], -1}, {-I Sin[x] + 
     Sin[y], 0, -1, 
    13/4 - Cos[x] - Cos[y] - Cos[z] - I Sin[z]}, {13/4 - Cos[x] - 
     Cos[y] - Cos[z] + I Sin[z], -1, 0, -I Sin[x] - Sin[y]}, {-1, 
    13/4 - Cos[x] - Cos[y] - Cos[z] + I Sin[z], I Sin[x] - Sin[y], 0}};
Oz[x_, y_,z_] :=
  {{0, 0, 0, 1/2 (-Cos[z] - I Sin[z])}, {0, 0, 
    1/2 (Cos[z] + I Sin[z]), 0}, {0, 1/2 (Cos[z] - I Sin[z]), 0, 
    0}, {1/2 I (I Cos[z] + Sin[z]), 0, 0, 0}};
Ox[x_, y_,z_] := {{0, I Cos[x], Sin[x], 0}, {-I Cos[x], 0, 0, 
    Sin[x]}, {Sin[x], 0, 0, -I Cos[x]}, {0, Sin[x], I Cos[x], 0}};
F2[x_, y_, z_, 
   r_] := (Inverse[(r + I*0.01)*IdentityMatrix[4] - KS[x, y, z]]);
F1[x_, y_, z_, 
   r_] := (Inverse[(r - I*0.01)*IdentityMatrix[4] - KS[x, y, z]]);
Myfun[x_, y_, z_, r_] := 
  Re[Tr[(Ox[x, y, z].(F2[x, y, z, r] - F1[x, y, z, r]).Oz[x, y, 
       z].(F2[x, y, z, r].F2[x, y, z, r] + 
        F1[x, y, z, r].F1[x, y, z, r]))]];      

Then I used ParallelSum with Table to do the job, but it is very slow> Is there any way to speed it up?

nr=5;nc=9;
Table[ParallelSum[(nc/(2\[Pi]))^-3 1/nr Myfun[x,y,z,r],{r,-7.25,R,1/nr},{x,\[Pi]/nc,\[Pi],(2\[Pi])/nc},
{y,\[Pi]/nc,\[Pi],(2\[Pi])/nc},{z,\[Pi]/nc,\[Pi],(2\[Pi])/nc}],{R,-1,1,0.5}]//AbsoluteTiming
    {9.9159,{1.64289,-1.40237,-1.40319,-1.40314,0.10968}}     

I want to do it with nr=100 and nc=200.

$\endgroup$
9
  • 1
    $\begingroup$ (1) Use nr = 5.; nc = 9.; in other words, use inexact values for nr and nc (2) let Parallelize choose the best approach by doing Parallelize[Table[Sum[.... This brings it down to 2.4s on my machine. (3) Eliminate some repeated calculations of F1 and F2 in Myfun by doing Myfun[x_, y_, z_, r_] := Block[{f1 = F1[x, y, z, r], f2 = F2[x, y, z, r]}, Re[Tr[(Ox[x, y, z].(f2 - f1).Oz[x, y, z].(f2.f2 + f1.f1))]]]; This cuts it down to just over a second. $\endgroup$
    – flinty
    Commented Aug 14, 2020 at 21:03
  • 1
    $\begingroup$ I can get it a little lower to 0.8 seconds, if I remove F1 and F2 and roll them into a single function to reduce recalculation of KS, as F1F2[x_, y_, z_, r_] := Block[{ks = KS[x, y, z], id = IdentityMatrix[4]}, {(Inverse[(r - I*0.01)*id - ks]), (Inverse[(r + I*0.01)*id -ks])}] then Myfun[x_, y_, z_, r_] := Block[{f = F1F2[x, y, z, r]}, Re[Tr[(Ox[x, y, z].(f[[2]] - f[[1]]).Oz[x, y, z].(f[[2]].f[[2]] + f[[1]].f[[1]]))]]]; then apply (1) and (2) from my previous comment. $\endgroup$
    – flinty
    Commented Aug 14, 2020 at 21:13
  • $\begingroup$ @flinty I would really appreciate if you can post it as an answer please so I can also accept it? $\endgroup$
    – MMA13
    Commented Aug 14, 2020 at 21:16
  • 2
    $\begingroup$ Is 0.08 seconds on a couple year old laptop sufficiently quick? $\endgroup$
    – ciao
    Commented Aug 14, 2020 at 21:16
  • $\begingroup$ @ciao, how did you do it, the results I posted were performed with a workstation with 20 core and 120G as RAM-:) So, 0.08 sec would be AWESOME $\endgroup$
    – MMA13
    Commented Aug 14, 2020 at 21:18

1 Answer 1

4
$\begingroup$

You can compile your function to speed up the evaluation:

MyfunCompiled = Compile[{{x, _Real}, {y, _Real}, {z, _Real}, {r, _Real}},
    Evaluate@Re[Tr[(Ox[x, y, z].(F2[x, y, z, r] - F1[x, y, z, r]).Oz[x, y, z].
        (F2[x, y, z, r].F2[x, y, z, r] + F1[x, y, z, r].F1[x, y, z, r]))]], 
    CompilationTarget -> "C"];

nr = 5; nc = 9;
AbsoluteTiming@Sum[(nc/(2 \[Pi]))^-3 1/nr Myfun[x, y, z, r], 
    {r, -7.25, 1, 1/nr}, {x, π/nc, π, 2π/nc}, {y, π/nc, π, 2π/nc}, {z, π/nc, π, 2π/nc}]

AbsoluteTiming@Sum[(nc/(2 \[Pi]))^-3 1/nr MyfunCompiled[x, y, z, r], 
    {r, -7.25, 1, 1/nr}, {x, π/nc, π, 2π/nc}, {y, π/nc, π, 2π/nc}, {z, π/nc, π, 2π/nc}]
{4.8844, 0.10968}
{0.050282, 0.10968}

For further improvements, you can also compile the actual sum:

FunSum = Compile[{{R, _Real}, {nr, _Integer}, {nc, _Integer}},
  Sum[(nc/(2π))^-3 1/nr MyfunCompiled[x, y, z, r], 
    {r, -7.25, 1, 1/nr}, {x, π/nc, π, 2π/nc}, {y, π/nc, π, 2π/nc}, {z, π/nc, π, 2π/nc}]
,CompilationTarget -> "C"];

nr = 5; nc = 9;
AbsoluteTiming@FunSum[1, nr, nc]
{0.007274, 0.10968}

Your example would then look like this (following the suggestion of flinty in comments about parallelizing over the different values of R):

nr = 5; nc = 9;
ParallelTable[FunSum[R, nr, nc], {R, -1, 1, 0.5}] // AbsoluteTiming
{0.027832, {1.64289, -1.40237, -1.40319, -1.40314, 0.10968}}

This is now fast enough that you can try going for larger values of nr and nc. For example:

AbsoluteTiming[FunSum[1, 50, 50]]
AbsoluteTiming[FunSum[1, 100, 100]]
AbsoluteTiming[FunSum[1, 100, 200]]
{5.55616, 0.452771}
{73.7241, 0.484847}
{678.142, 0.48616}

which seems to converge (at least it looks like it does).

$\endgroup$
2
  • $\begingroup$ that is Great!, by the way, will it help if we remove F1 and F2 and roll them into a single function to reduce recalculation of KS,as suggested by "flinty" in the comments, and combine it with your approach ? $\endgroup$
    – MMA13
    Commented Aug 15, 2020 at 18:04
  • 1
    $\begingroup$ It is unfortunately not directly applicable to my method, since I am basically compiling a single, long function, and all the higher-level operations like Inverse are precomputed. I did just try a version where I only precompute F1, F2, Ox and Oz, and the compiled function then assembles the Re[Tr[Dot[]]], but there is essentially no improvement in the timings. $\endgroup$
    – Hausdorff
    Commented Aug 15, 2020 at 19:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.