1
$\begingroup$

I would like to caculate the following codes with four variables. When the dimension is 10*10*10*10, the complile time is about 2 seconds. (just for estimate the time). When the dimension the 50*50*50*50, the complile time is about 1200 seconds, (20 mintutes), which is too long time for me, because I need to repeat this calcualtion for many times.

How can I speed up this compile time for 50*50*50*50 dimensions? Any help or suggestion will be highly appreciated!

data = Table[Exp[-((i + j - 50.)/10)^2] Exp[-((i - j)/10)^2], {i, 50}, {j, 50}];
data = Chop[data, 0.00001];
data = data/Sqrt[Sum[(data[[i, j]])^2, {i, 1, 50}, {j, 1, 50}]];


c = 2.99792458*10^8;
Δ = 0.16;(*nm*)
λs0 = 1500;
ListDensityPlot[data, InterpolationOrder -> 0, Mesh -> All, PlotRange -> All, ColorFunction -> (Blend[{Hue[2/3], Hue[0]}, #] &)]

CC1[i_, j_, k_, l_, τ_] := 1/2 (2 data[[i, l]]^2 data[[j, k]]^2 + 
 2 data[[i, k]] data[[i, l]] data[[j, k]] data[[j, l]] + 
 2 data[[i, k]]^2 data[[j, l]]^2 + 
 Cos[((2 c π)/(-4 - 0.5` Δ + 
       i Δ + λs0) + (
      2 c π)/(-4 - 0.5` Δ + 
       j Δ + λs0) - (
      2 c π)/(-4 - 0.5` Δ + 
       k Δ + λs0) - (
      2 c π)/(-4 - 0.5` Δ + 
       l Δ + λs0)) τ] (data[[i, 
      l]] data[[j, k]] + data[[i, k]] data[[j, l]])^2 + 
 2 data[[i, j]] data[[i, l]] data[[j, k]] data[[k, l]] + 
 2 data[[i, j]] data[[i, k]] data[[j, l]] data[[k, l]] + 
 2 data[[i, j]]^2 data[[k, l]]^2 + 
 2 Cos[((2 c π)/(-4 - 0.5` Δ + 
       j Δ + λs0) - (
      2 c π)/(-4 - 0.5` Δ + 
       k Δ + λs0)) τ] (data[[i, 
      l]] data[[j, k]] + 
    data[[i, k]] data[[j, l]]) (data[[i, l]] data[[j, k]] + 
    data[[i, j]] data[[k, l]]) + 
 2 Cos[((2 c π)/(-4 - 0.5` Δ + 
       i Δ + λs0) - (
      2 c π)/(-4 - 0.5` Δ + 
       l Δ + λs0)) τ] (data[[i, 
      l]] data[[j, k]] + 
    data[[i, k]] data[[j, l]]) (data[[i, l]] data[[j, k]] + 
    data[[i, j]] data[[k, l]]) + 
 Cos[((2 c π)/(-4 - 0.5` Δ + 
       i Δ + λs0) - (
      2 c π)/(-4 - 0.5` Δ + 
       j Δ + λs0) + (
      2 c π)/(-4 - 0.5` Δ + 
       k Δ + λs0) - (
      2 c π)/(-4 - 0.5` Δ + 
       l Δ + λs0)) τ] (data[[i, 
      l]] data[[j, k]] + data[[i, j]] data[[k, l]])^2 + 
 2 Cos[((2 c π)/(-4 - 0.5` Δ + 
       i Δ + λs0) - (
      2 c π)/(-4 - 0.5` Δ + 
       k Δ + λs0)) τ] (data[[i, 
      l]] data[[j, k]] + 
    data[[i, k]] data[[j, l]]) (data[[i, k]] data[[j, l]] + 
    data[[i, j]] data[[k, l]]) + 
 2 Cos[((2 c π)/(-4 - 0.5` Δ + 
       j Δ + λs0) - (
      2 c π)/(-4 - 0.5` Δ + 
       l Δ + λs0)) τ] (data[[i, 
      l]] data[[j, k]] + 
    data[[i, k]] data[[j, l]]) (data[[i, k]] data[[j, l]] + 
    data[[i, j]] data[[k, l]]) + 
 2 Cos[((2 c π)/(-4 - 0.5` Δ + 
       i Δ + λs0) - (
      2 c π)/(-4 - 0.5` Δ + 
       j Δ + λs0)) τ] (data[[i, 
      l]] data[[j, k]] + 
    data[[i, j]] data[[k, l]]) (data[[i, k]] data[[j, l]] + 
    data[[i, j]] data[[k, l]]) + 
 2 Cos[((2 c π)/(-4 - 0.5` Δ + k Δ + λs0) - (2 c π)/(-4 - 0.5` Δ + l Δ + λs0)) τ] (data[[i, l]] data[[j, k]] + 
    data[[i, j]] data[[k, l]]) (data[[i, k]] data[[j, l]] + 
    data[[i, j]] data[[k, l]]) + 
 Cos[((2 c π)/(-4 - 0.5` Δ + i Δ + λs0) - (2 c π)/(-4 - 0.5` Δ + j Δ + λs0) - (2 c π)/(-4 - 0.5` Δ + k Δ + λs0) + (2 c π)/(-4 - 0.5` Δ + l Δ + λs0)) τ] (data[[i, k]] data[[j, l]] + data[[i, j]] data[[k, l]])^2);

CC2comp = Compile[{{t, _Real}}, 
Evaluate[
 Sum[CC1[i, j, k, l, t], {i, 10}, {j, 10}, {k, 10}, {l, 10}]]]; // AbsoluteTiming

ListPlot[Table[{i, CC2comp[i*0.001]}, {i, -10, 10, 0.1}], Joined -> True, Mesh -> All, Axes -> None, PlotRange -> {All, All}, Frame -> True, AspectRatio -> .67, ImageSize -> {400, 250}] // AbsoluteTiming


CC2comp = Compile[{{t, _Real}},  Evaluate[
 Sum[CC1[i, j, k, l, t], {i, 50}, {j, 50}, {k, 50}, {l, 50}]]]; // AbsoluteTiming

ListPlot[Table[{i, CC2comp[i*0.001]}, {i, -10, 10, 0.1}], Joined -> True, Mesh -> All, Axes -> None, PlotRange -> {All, All},  Frame -> True, AspectRatio -> .67,  ImageSize -> {400, 250}] // AbsoluteTiming

enter image description here

$\endgroup$
  • $\begingroup$ @Shutao Tang, Many thanks to Shutao Tang for re-editing! How can I input Δ here? I just copy the code from a Mathematica notebook, and the result is [CapitalDelta]. $\endgroup$ – user14634 Sep 6 '16 at 15:36
  • $\begingroup$ See Mathematica Editor-buttons for StackExchange for chrome. $\endgroup$ – xyz Sep 7 '16 at 0:49
  • $\begingroup$ @Shutao Tang, Thanks a lot! $\endgroup$ – user14634 Sep 7 '16 at 9:05
3
$\begingroup$

The reason your code is slow is due to the Evaluate you have placed in the body of Compile: this forces the Sum to evaluate with a symbolic t in the place of τ, and use the resulting symbolic formula as the body of Compile. This symbolic processing is slow and unnecessary.

My proposal is to make a Compiled function of CC1 and then calling that from the compiled CC2:

CC1 = Compile[{{data, _Real, 
    2}, {i, _Integer}, {j, _Integer}, {k, _Integer}, {l, _Integer}, {\
τ, _Real}}, 
  1/2 (2 data[[i, l]]^2 data[[j, k]]^2 + 
     2 data[[i, k]] data[[i, l]] data[[j, k]] data[[j, l]] + 
     2 data[[i, k]]^2 data[[j, l]]^2 + 
     Cos[((2 c π)/(-4 - 0.5` Δ + 
             i Δ + λs0) + (2 c π)/(-4 - 
             0.5` Δ + 
             j Δ + λs0) - (2 c π)/(-4 - 
             0.5` Δ + 
             k Δ + λs0) - (2 c π)/(-4 - 
             0.5` Δ + 
             l Δ + λs0)) τ] (data[[i, 
            l]] data[[j, k]] + data[[i, k]] data[[j, l]])^2 + 
     2 data[[i, j]] data[[i, l]] data[[j, k]] data[[k, l]] + 
     2 data[[i, j]] data[[i, k]] data[[j, l]] data[[k, l]] + 
     2 data[[i, j]]^2 data[[k, l]]^2 + 
     2 Cos[((2 c π)/(-4 - 0.5` Δ + 
             j Δ + λs0) - (2 c π)/(-4 - 
             0.5` Δ + 
             k Δ + λs0)) τ] (data[[i, 
           l]] data[[j, k]] + 
        data[[i, k]] data[[j, l]]) (data[[i, l]] data[[j, k]] + 
        data[[i, j]] data[[k, l]]) + 
     2 Cos[((2 c π)/(-4 - 0.5` Δ + 
             i Δ + λs0) - (2 c π)/(-4 - 
             0.5` Δ + 
             l Δ + λs0)) τ] (data[[i, 
           l]] data[[j, k]] + 
        data[[i, k]] data[[j, l]]) (data[[i, l]] data[[j, k]] + 
        data[[i, j]] data[[k, l]]) + 
     Cos[((2 c π)/(-4 - 0.5` Δ + 
             i Δ + λs0) - (2 c π)/(-4 - 
             0.5` Δ + 
             j Δ + λs0) + (2 c π)/(-4 - 
             0.5` Δ + 
             k Δ + λs0) - (2 c π)/(-4 - 
             0.5` Δ + 
             l Δ + λs0)) τ] (data[[i, 
            l]] data[[j, k]] + data[[i, j]] data[[k, l]])^2 + 
     2 Cos[((2 c π)/(-4 - 0.5` Δ + 
             i Δ + λs0) - (2 c π)/(-4 - 
             0.5` Δ + 
             k Δ + λs0)) τ] (data[[i, 
           l]] data[[j, k]] + 
        data[[i, k]] data[[j, l]]) (data[[i, k]] data[[j, l]] + 
        data[[i, j]] data[[k, l]]) + 
     2 Cos[((2 c π)/(-4 - 0.5` Δ + 
             j Δ + λs0) - (2 c π)/(-4 - 
             0.5` Δ + 
             l Δ + λs0)) τ] (data[[i, 
           l]] data[[j, k]] + 
        data[[i, k]] data[[j, l]]) (data[[i, k]] data[[j, l]] + 
        data[[i, j]] data[[k, l]]) + 
     2 Cos[((2 c π)/(-4 - 0.5` Δ + 
             i Δ + λs0) - (2 c π)/(-4 - 
             0.5` Δ + 
             j Δ + λs0)) τ] (data[[i, 
           l]] data[[j, k]] + 
        data[[i, j]] data[[k, l]]) (data[[i, k]] data[[j, l]] + 
        data[[i, j]] data[[k, l]]) + 
     2 Cos[((2 c π)/(-4 - 0.5` Δ + 
             k Δ + λs0) - (2 c π)/(-4 - 
             0.5` Δ + 
             l Δ + λs0)) τ] (data[[i, 
           l]] data[[j, k]] + 
        data[[i, j]] data[[k, l]]) (data[[i, k]] data[[j, l]] + 
        data[[i, j]] data[[k, l]]) + 
     Cos[((2 c π)/(-4 - 0.5` Δ + 
             i Δ + λs0) - (2 c π)/(-4 - 
             0.5` Δ + 
             j Δ + λs0) - (2 c π)/(-4 - 
             0.5` Δ + 
             k Δ + λs0) + (2 c π)/(-4 - 
             0.5` Δ + 
             l Δ + λs0)) τ] (data[[i, 
            k]] data[[j, l]] + data[[i, j]] data[[k, l]])^2)
  , CompilationTarget -> "C"
  , CompilationOptions -> {"InlineExternalDefinitions" -> True}
  ]

and

CC2comp = Compile[{{t, _Real}, {dat, _Real, 2}},
    Sum[CC1[dat, i, j, k, l, t], {i, 50}, {j, 50}, {k, 50}, {l, 50}]
    ,
    CompilationTarget -> "C", 
    CompilationOptions -> {"InlineExternalDefinitions" -> 
       True}]; // AbsoluteTiming

{1.556052, Null}

You can drop the compilation to C if you don't have a working compiler but for such low-level operations as here, compiling to C should help quite a bit.

Now, the plot takes some time to create, but that's expected since each plot points needs the 50*50*50*50 sum to be evaluated, but it's still much less than waiting for the symbolic sum.

ListPlot[Table[{i, CC2comp[i*0.001, data]}, {i, -10, 10, 0.1}], 
  Joined -> True, Mesh -> All, Axes -> None, PlotRange -> {All, All}, 
  Frame -> True, AspectRatio -> .67, 
  ImageSize -> {400, 250}] // AbsoluteTiming

enter image description here

An alternative way would be to try to speed of the generation of the symbolic result of the sum. To approach this I would try to vectorize some of the sums if possible, do vector operations, and use Total instead of Sum, but I don't have time to investigate this right now :/ If I have time later, I'll give it a go.

$\endgroup$
  • $\begingroup$ Hi Marius Ladegård Meyer, Thank you so much for your good answer! $\endgroup$ – user14634 Sep 6 '16 at 15:33
  • $\begingroup$ I am also looking forward to seeing the AbsoluteTiming results for codes with Total instead of Sum. $\endgroup$ – user14634 Sep 6 '16 at 15:38
  • $\begingroup$ Can the speed be faster if I use Matlab codes instead of Mathematica codes? $\endgroup$ – user14634 Sep 7 '16 at 8:38

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.