3

I have to reassign values in a big array many times and it turned out to be the most time consuming in my program. The best way I could figure out is this

nn=1024;nnn=1024/2-1;t=9.4122*10^-6;mass=RandomComplex[{1,I},{nn,nn}];

mass = ParallelTable[mass[[k, k1]] E^(I*t*(Mod[k - 1, nn, -nnn]^2 + 
   Mod[k1 - 1, nn, -nnn]^2)), {k, 1, nn}, {k1, 1, nn}];

I wonder if there is a faster functional way to change the values in an array with a function depending on their positions.

I tried to use Compile for this function

f = Compile[{{xp, _Complex, 2}}, ParallelTable[xp[[k,k1]] 
E^(I*t*(Mod[k - 1, nn, -nnn]^2 + Mod[k1 - 1, nn, -nnn]^2)), {k, 1, nn}, {k1, 1, nn}],
{{k, _Integer}, {k1, _Integer}, {nn, _Integer}, {nnn, _Integer}, {Mod[_, _, _], _Integer}}]

But I get the answer that Compiled expression should be a mashine-size Complex number. Please tell me what is wrong.

xzczd
  • 65,995
  • 9
  • 163
  • 468
user40532
  • 215
  • 1
  • 6
  • Is this calculation related to fast Fourier transform? – xzczd Jan 19 '17 at 07:46
  • 2
    Am I interpreting this correctly: you are in your code doing the mass=ParallelTable[...] statement repeatedly? If so, why on earth would you recalculate the same factors repeatedly? Simply calculate those once (i.e. remove the mass[[k, k1]] piece, give result a name), and at each iteration to change mass it becomes mass=a_name*mass... – ciao Jan 19 '17 at 09:16
  • Yes, it is related to fft. I am doing fft with CUDAFourier and the code above has to change the matrix in the momentum space, then inverse Fourier, and similar change in the coordinate space. – user40532 Jan 19 '17 at 10:35
  • @ciao, Yes, that is very good, I just was afraid of storing another big array because memory after many iterations becomes a big issue. – user40532 Jan 19 '17 at 10:41

1 Answers1

4
rotate = l \[Function] RotateLeft[Range[l] - # - 1, #] &@Floor[(l - 1)/2];    
mass3 = (lst \[Function] mass (Exp[I t ( lst^2 + #)] & /@ (lst^2)))@rotate@nn;
(* {0.519567, Null} *)

As to Compile, I don't suggest you to use it now, but if you insist, you can start from here. Here's the compiled version of the code above:

cf = Compile[{{mass, _Complex, 2}, {lst, _Real, 1}, {t, _Real}}, 
   mass (Exp[I t ( lst^2 + #)] & /@ (lst^2))];   
mass4 = cf[mass, rotate@nn, t]; // AbsoluteTiming
(* {0.103073, Null} *)

Or, let's naively compile your implementation:

cf2 = 
  Compile[{{mass, _Complex, 2}, {t, _Real}, {nn, _Integer}, {nnn, _Integer}}, 
   Table[mass[[k, k1]] E^(I t (Mod[k - 1, nn, -nnn]^2 + Mod[k1 - 1, nn, -nnn]^2)), {k, 1,
      nn}, {k1, 1, nn}]];

mass5 = cf2[mass, t, nn, nnn]; // AbsoluteTiming
(* {0.461346, Null} *)

It can be faster if we use some even more advanced technique (related post: Why is CompilationTarget -> C slower than directly writing with C?):

cf3 = 
  Unevaluated@
     Compile[{{mass, _Complex, 2}, {t, _Real}, {nn, _Integer}, {nnn, _Integer}}, 
      Table[mass[[k, k1]] E^(I t (Mod[k - 1, nn, -nnn]^2 + Mod[k1 - 1, nn, -nnn]^2)), {k,
         1, nn}, {k1, 1, nn}], CompilationTarget -> "C", RuntimeOptions -> "Speed"] /. 
    Part -> Compile`GetElement // Last;

mass6 = cf3[mass, t, nn, nnn]; // AbsoluteTiming
(* {0.115824, Null} *)
xzczd
  • 65,995
  • 9
  • 163
  • 468