# (Too) Slow Compiled Function To Sum Waves

I have some code that sums planar waves and outputs a graphical result, as per the mock up below :

waveZoomFactor = 1;
waveHorizontalScaling = 1;
waveVerticalScaling = 1;
imageSize = 100;

createWaves =
Compile[{{x, _Real}, {y, _Real}, {sym, _Real}, {phaseShift, _Real}},
Table[Cos[
waveZoomFactor*(waveHorizontalScaling*x *Cos[2*Pi*(k/sym)] +
waveVerticalScaling*y *Sin[2*Pi (k/sym)]) + phaseShift], {k,
0, sym - 1}]
, CompilationTarget -> "C", RuntimeAttributes -> {Listable},
Parallelization -> True, RuntimeOptions -> "Speed"];

waves[phaseShift_] :=
ParallelTable[Total[createWaves[x, y, 5  , phaseShift]],
{x, -imageSize , imageSize , 0.2 }, {y, -imageSize , imageSize , 0.2 }];

CreatePattern := Table[Image[waves[t]], {t, 0, 2*Pi, 1}]
CreatePattern


I can't figure out why the code is so incredibly slow, in particular for larger 'imageSize'. I am a bit rusty with Mathematica so I am wondering if I overlooked something obvious or something dumb is going on. Thanks in advance.

• Wouldn't it be better to use Sum inside the createWaves function instead of returning list with Table, since you're going to sum it anyway with Total? Dec 29, 2022 at 16:39
• Definitively made a difference, thanks for the help. Dec 29, 2022 at 19:48
• Have you read this post?: mathematica.stackexchange.com/a/104031/1871 Dec 30, 2022 at 2:27

## 1 Answer

If you've got OpenCL support, here's an OpenCL implementation which is almost instantaneous in returning fairly high resolution results:

Remove["Global*"];
Needs["OpenCLLink"];

Check[OpenCLQ[], "OpenCL not available!"];

cwSource =
"__kernel void waves(__global float * out, float sym, float phaseShift, mint xres, mint yres, float x0, float y0, float x1, float y1, float zoom, float waveXscale, float waveYscale) {
int xi = get_global_id(0);
int yi = get_global_id(1);
float xfrac = ((float)xi)/xres;
float yfrac = ((float)yi)/yres;
float x = mix(x0,x1,xfrac);
float y = mix(y0,y1,yfrac);
float result = 0.0f;
for(int k = 0; k < sym; ++k)
{
float f = 2.0f*M_PI_F*k/sym;
result += cos(zoom*(waveXscale*x*cos(f) + waveYscale*y*sin(f)) + phaseShift);
}
out[yi * xres + xi] = result;
}
";

createWavesOpenCL = OpenCLFunctionLoad[cwSource, "waves",
{{"Float", _, "Output"}, "Float", "Float", _Integer, _Integer,
"Float", "Float", "Float", "Float", "Float", "Float",
"Float"}, {8, 8}, "ShellOutputFunction" -> Print];

waveZoomFactor = 1;
waveHorizontalScaling = 1;
waveVerticalScaling = 1;
resolution = 500;
minx = -100;
maxx = 100;
miny = -100;
maxy = 100;
output = OpenCLMemoryAllocate["Float", {resolution, resolution}];
Table[
createWavesOpenCL[output, 5, phase, resolution, resolution, minx,
miny, maxx, maxy, waveZoomFactor, waveHorizontalScaling,
waveVerticalScaling];
Image[OpenCLMemoryGet[output], "Float"],
{phase, 0, 2 Pi, 1}]


This is fast enough, that I can generate about 256 300x300 image slices and stack them in a 3D image using Image3D. It takes just a few seconds on an nVidia 1080ti.

• Thanks a lot. this is a great approach, I'll try to up my skills and use similar code moving forward. Now if I could only get CUDA set up... Dec 29, 2022 at 19:48