1

I would like to generate a set of data of matrix of size 1000 x 1000 x 1000 in a region defined by 0 outside the sphere of radius r=100 centered in the coordinate (500, 500, 500), and number 1 inside the sphere.

n = 1000;(*Matrix dimension*)
f = Compile[{i, j, k},If[(i-500)^2 + (j - 500)^2 + (k - 500)^2 <= (10)^2, 1, 0]];
ic = ParallelTable[f[N@i, N@j, N@k], {i, n}, {j, n}, {k, n}];

Next, to apply Fast Fourier Transform to data.

fftshift[dat_?ArrayQ, k : (_Integer?Positive | All) : All] := Module[{dims = Dimensions[dat]}, 
   RotateRight[dat,If[k === All, Quotient[dims, 2],Quotient[dims[[k]], 2] UnitVector[Length[dims], k]]]];
ict = Table[fftshift@ic[[i]], {i, Length[ic]}];

Finally to export the data.

Export["ic50.mat", ict]

However, the time to conclude the process is huge.

How to optimize this to reduce run time?

SAC
  • 1,335
  • 8
  • 17
  • 1
    I don't know what you mean to do with the Fast Fourier Transform. Your code does not involve any. Though you should be able to use the command Fourier for that. Applied to a multidimensional tensor, it applies FFT to each of its slots. – Henrik Schumacher Nov 16 '21 at 22:43
  • @xzczd thank you. I was looking for something like this code. – SAC Nov 20 '21 at 14:22

1 Answers1

3

The fastest solution for the array ci that I found so far is this one:

cf = Compile[{{x, _Real}, {y, _Real, 1}, {z, _Real, 1}, {r, _Real}},
   Block[{r2},
    r2 = r r - x;
    Table[
     If[Compile`GetElement[y, i] + Compile`GetElement[z, j] <= r2, 1.,0.],
     {i, 1, Length[y]},
     {j, 1, Length[z]}
     ]
    ],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

n = 1000.;
r = 100.;
x = (Range[1., n, 1.] - n/2.)^2;
ci = cf[x, x, x, r];

It takes about 3 seconds on my Quad Core machine (a dated Haswell). The idea is to perform many computations only once for the vector x and to exploint parallelization in Compile enabled by Listable. The output has size 8 GB, so performance depends critically on the available RAM. My machine has only 16 GB, so I run into heavy swapping if I generate ci and then do some further manipulations with it.

Henrik Schumacher
  • 106,770
  • 7
  • 179
  • 309