<jupyter_start><jupyter_text>batched 1D mode in 3 separate goes does not work! (need work-around for maximum of 2^{27} entries in plan for 3D with 8x Hockney algorithm in double) http://stackoverflow.com/questions/26918101/1d-ffts-of-columns-and-rows-of-a-3d-matrix-in-cuda  however, the following method might work for a work-around: - batch 2D x-y FFTs until plan is full, then launch next batch of x-y FFTs- launch len(y) many batches with length len(x) each of 1D z FFTs - check that transposing is not quicker for memory alignment for the last operation?  2D<jupyter_code>a2d = a3d[1].astype(np.complex128)
nx = 32#1024#*4
ny = 16#1024#*4
a2d = gpuarray.zeros((ny, nx), dtype=np.complex128) + 1
y, x = map(gpuarray.to_gpu, 
           np.meshgrid(np.arange(a2d.shape[0], dtype=np.float64), 
                       np.arange(a2d.shape[1], dtype=np.float64))
              )
a2d *= cumath.sin(2*x) * cumath.sin(y)
plt.contourf(np.real(a2d.get()))
plt.colorbar()
plan2d = fft.Plan(shape=a2d.shape, in_dtype=np.complex128, out_dtype=np.complex128)
res = gpuarray.empty(a2d.shape, dtype=np.complex128) - 42.
%%timeit
fft.fft(a2d, res, plan2d)
plt.imshow(np.abs(res.get()), interpolation='none', origin='lower')
plt.colorbar()<jupyter_output><empty_output><jupyter_text>FFT along x direction (first)<jupyter_code>plan_1 = fft.Plan(shape=a2d.shape[1], in_dtype=np.complex128, 
                   out_dtype=np.complex128, batch=a2d.shape[0],
                   inembed=np.asarray(a2d.shape), onembed=np.asarray(a2d.shape),
                   istride=1, idist=a2d.shape[1], 
                   ostride=1, odist=a2d.shape[1])
res = gpuarray.empty(a2d.shape, dtype=np.complex128) - 42.
%%timeit
fft.fft(a2d, res, plan_1)
context.synchronize()
plt.imshow(np.abs(res.get()), interpolation='none', origin='lower')
plt.colorbar()
np_res = np.fft.fft(a2d.get(), axis=1)
plt.imshow(np.abs(np_res), interpolation='none', origin='lower')
plt.colorbar()<jupyter_output><empty_output><jupyter_text>alternative first FFT along y direction<jupyter_code>plan_1a = fft.Plan(shape=a2d.shape[0], in_dtype=np.complex128, 
                   out_dtype=np.complex128, batch=a2d.shape[1],
                   inembed=np.asarray(a2d.shape), onembed=np.asarray(a2d.shape),
                   istride=a2d.shape[1], idist=1, 
                   ostride=a2d.shape[1], odist=1)
resa = gpuarray.empty(a2d.shape, dtype=np.complex128) - 42.
fft.fft(a2d, resa, plan_1a)

plt.imshow(np.abs(resa.get()), interpolation='none', origin='lower')
plt.colorbar()
np_res = np.fft.fft(a2d.get(), axis=0)
plt.imshow(np.abs(np_res), interpolation='none', origin='lower')
plt.colorbar()
np.allclose(np_res, res.get())<jupyter_output><empty_output><jupyter_text>FFT along y direction<jupyter_code>plan_2 = fft.Plan(shape=a2d.shape[0], in_dtype=np.complex128, 
                  out_dtype=np.complex128, batch=a2d.shape[1],
                  inembed=np.asarray(a2d.shape), onembed=np.asarray(a2d.shape),
                  idist=1, istride=a2d.shape[1], 
                  odist=1, ostride=a2d.shape[1])
res2 = gpuarray.empty(a2d.shape, dtype=np.complex128) - 42.
%%timeit
fft.fft(res, res2, plan_2)
plt.imshow(np.abs(res2.get()), interpolation='none', origin='lower')
plt.colorbar()
a2d_h = a2d.get()
%%timeit
np.fft.fft2(a2d_h)
np_res = np.fft.fft2(a2d.get())
plt.imshow(np.abs(np_res), interpolation='none', origin='lower')
plt.colorbar()
np.allclose(res2.get(), np_res)
fftw_1 = pyfftw.builders.fft(a2d_h, axis=0)
fftw_2 = pyfftw.builders.fft(a2d_h, axis=1)
%%timeit
fftw_res = fftw_1(a2d_h)
fftw_res = fftw_2(fftw_res)
fftw_res = fftw(a2d.get())
np.allclose(fftw_res, np_res)
pyfftw.interfaces.cache.enable()
pyfftw.interfaces.cache.set_keepalive_time(30)
%%timeit
pyfftw.interfaces.numpy_fft.fft2(a2d_h)
25.1/1.17<jupyter_output><empty_output><jupyter_text>1D<jupyter_code>a1d = a3d[1,1].astype(np.complex128)
plt.plot(np.real(a1d.get()))
plan = fft.Plan(shape=a1d.shape, in_dtype=np.complex128, out_dtype=np.complex128)
res = gpuarray.zeros(a1d.shape, dtype=np.complex128) - 42.
fft.fft(a1d, res, plan)
plt.plot(np.abs(res.get())[:(a1d.shape[0]/2 + 1)])
plt.plot(np.abs(np.fft.fft(a1d.get()))[:(a1d.shape[0]/2 + 1)])<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup