diff --git a/ci/scripts/test.sh b/ci/scripts/test.sh
index be5da252b2517262dcc0c171717b0c82520f6ddb..5c3f1f9e59579f1a87984aec1f7fe2a9703e969e 100755
--- a/ci/scripts/test.sh
+++ b/ci/scripts/test.sh
@@ -52,6 +52,7 @@ python -c 'import hysop; print hysop'
python "$HYSOP_DIR/core/arrays/tests/test_array.py"
$HYSOP_DIR/fields/tests/test_cartesian.sh
+python "$HYSOP_DIR/numerics/tests/test_fft.py"
python "$HYSOP_DIR/operator/tests/test_analytic.py"
python "$HYSOP_DIR/operator/tests/test_transpose.py"
python "$HYSOP_DIR/operator/tests/test_derivative.py"
diff --git a/cmake/hysop_tests.cmake b/cmake/hysop_tests.cmake
index 8288fa8f9ffadb19afdfd797808e88f0c5de8f6e..9978521f9a1799f500eae7b2636052a09edd1609 100755
--- a/cmake/hysop_tests.cmake
+++ b/cmake/hysop_tests.cmake
@@ -52,7 +52,7 @@ endmacro()
set(testDir ${HYSOP_BUILD_PYTHONPATH})
# === Set the list of all directories which may contain tests ===
-set(py_src_dirs core/arrays fields operator)
+set(py_src_dirs core/arrays numerics fields operator)
# === Create the files list from all directories in py_src_dirs ===
diff --git a/hysop/numerics/fft/gpyfft_fft.py b/hysop/numerics/fft/gpyfft_fft.py
index e454eb68af3316fd256a1c2b738394d436485a53..b6aef6cf4b8c34b06f735e4e402634e0679063bd 100644
--- a/hysop/numerics/fft/gpyfft_fft.py
+++ b/hysop/numerics/fft/gpyfft_fft.py
@@ -4,7 +4,7 @@ FFT iterface for fast Fourier Transforms using CLFFT backend.
:class:`~hysop.numerics.GpyFFTPlan`
"""
-import warnings, struct
+import warnings, struct, primefac
import numpy as np
from abc import abstractmethod
@@ -250,6 +250,7 @@ class GpyFFTPlan(FFTPlanI):
(out_array.strides[t_axes_in[0]] == out_array.dtype.itemsize)), \
'Inplace real transforms need stride 1 for first transform axis.'
+ self.check_transform_shape(t_shape)
plan = GFFT.create_plan(self.context, t_shape[::-1])
plan.inplace = t_inplace
plan.strides_in = t_strides_in[::-1]
@@ -342,7 +343,7 @@ t_inplace, t_shape, t_precision,
layout_in, layout_out, plan.scale_forward, plan.scale_backward,
self.pre_callback_src, self.post_callback_src)
print msg
-
+
def set_callbacks(self, plan, axes, N,
in_array, out_array, fake_input, fake_output,
layout_in, layout_out, **kwds):
@@ -380,6 +381,24 @@ self.pre_callback_src, self.post_callback_src)
if (post_src is not None):
plan.set_callback('post_callback', post_src, 'post', user_data=user_data)
return (in_data, out_data)
+
+ @classmethod
+ def check_transform_shape(self, shape):
+ """Check that clFFT can handle the logical transform size."""
+ valid_factors = {2,3,5,7,11,13}
+ for Ni in shape:
+ factors = tuple( primefac.primefac(int(Ni)) )
+ invalid_factors = set(factors) - valid_factors
+ if invalid_factors:
+ factorization = ' * '.join('{}^{}'.format(factor, factors.count(factor))
+ for factor in set(factors))
+ candidates = ', '.join(str(vf) for vf in valid_factors)
+ msg ='\nInvalid transform shape {} for clFFT:'
+ msg+='\n {} = {}'
+ msg+='\nOnly {} prime factors are available.'
+ msg+='\n'
+ msg=msg.format(shape, Ni, factorization, candidates)
+ raise ValueError(msg)
@classmethod
def calculate_transform_strides(cls, taxes, array):
diff --git a/hysop/numerics/tests/test_fft.py b/hysop/numerics/tests/test_fft.py
index 0ea1a2bce9224b307350fd8b6e60bb0f85b5a20d..17e81d67a3e8709acb2074ca58b4fd869985b9df 100644
--- a/hysop/numerics/tests/test_fft.py
+++ b/hysop/numerics/tests/test_fft.py
@@ -123,7 +123,7 @@ class TestFFT(object):
print '\n FORWARD R2C: real to hermitian complex transform'
for (shape, cshape, rshape, N, Nc, Nr) in self.iter_shapes():
print ' RFFT shape={:9s} '.format(str(shape)+':'),
- Href = np.random.rand(N).astype(dtype).reshape(shape)
+ Href = np.random.rand(*shape).astype(dtype).reshape(shape)
results = {}
for (kind, implementations) in self.implementations.iteritems():
for (name, impl) in implementations.iteritems():
@@ -444,7 +444,7 @@ class TestFFT(object):
print '\n R2C-C2R transform'
for (shape, cshape, rshape, N, Nc, Nr) in self.iter_shapes():
print ' X - IRFFT(RFFT(X)) shape={:9s} '.format(str(shape)+':'),
- Href = np.random.rand(N).astype(dtype).reshape(shape)
+ Href = np.random.rand(*shape).astype(dtype).reshape(shape)
results = {}
for (kind, implementations) in self.implementations.iteritems():
for (name, impl) in implementations.iteritems():
@@ -467,8 +467,8 @@ class TestFFT(object):
results[name] = np.max(np.abs(Href-H0))
check_distances(results)
- return
print '\n R2R-R2R transforms'
+
types = ['I ','II ','III','IV ']
for (itype,stype) in enumerate(types, 1):
print '\n DCT-{}: real to real discrete cosine transform {}'.format(
@@ -476,64 +476,81 @@ class TestFFT(object):
ttype = 'COS{}'.format(itype)
for (shape, cshape, rshape, N, Nc, Nr) in self.iter_shapes():
print ' X - I{}({}(X)) shape={:9s} '.format(ttype, ttype, str(shape)+':'),
- Href = np.random.rand(N).astype(dtype).reshape(shape)
- H0 = np.empty(shape=shape, dtype=dtype)
- H1 = np.empty(shape=shape, dtype=dtype)
- H0, H1, Href = self.simd_align(H0, H1, Href)
+ if (itype==1): # real size is 2*(N-1)
+ shape = mk_shape(shape, -1, shape[-1] + 1)
+ Href = np.random.rand(*shape).astype(dtype).reshape(shape)
results = {}
- for (name, impl) in self.implementations.iteritems():
- iitype = [1,3,2,4][itype-1]
- if dtype not in impl.supported_ftypes:
- continue
- if itype not in impl.supported_cosine_transforms:
- continue
- if iitype not in impl.supported_cosine_transforms:
- continue
- forward = impl.dct(a=H0, out=H1, type=itype)
- backward = impl.idct(a=H1, out=H0, type=itype)
- H0[...] = Href
- forward.execute()
- backward.execute()
- results[name] = np.max(np.abs(Href-H0))
+ for (kind, implementations) in self.implementations.iteritems():
+ for (name, impl) in implementations.iteritems():
+ iitype = [1,3,2,4][itype-1]
+ if dtype not in impl.supported_ftypes:
+ continue
+ if itype not in impl.supported_cosine_transforms:
+ continue
+ if iitype not in impl.supported_cosine_transforms:
+ continue
+ D0 = impl.backend.empty(shape=shape, dtype=dtype)
+ D1 = impl.backend.empty(shape=shape, dtype=dtype)
+ forward = impl.dct(a=D0, out=D1, type=itype).setup()
+ backward = impl.idct(a=D1, out=D0, type=itype).setup()
+ assert forward.input_array is D0
+ assert forward.output_array is D1
+ assert backward.input_array is D1
+ assert backward.output_array is D0
+ D0[...] = Href
+ D1[...] = np.nan
+ forward.execute()
+ D0[...] = np.nan
+ backward.execute()
+ H0 = D0.get()
+ results[name] = np.max(np.abs(Href-H0))
check_distances(results)
for (itype,stype) in enumerate(types, 1):
- print '\n DST-{}: real to real discrete sine transform {}'.format(
+ print '\n DST-{}: real to real discrete sinine transform {}'.format(
stype.strip(), itype)
ttype = 'SIN{}'.format(itype)
for (shape, cshape, rshape, N, Nc, Nr) in self.iter_shapes():
print ' X - I{}({}(X)) shape={:9s} '.format(ttype, ttype, str(shape)+':'),
- Href = np.random.rand(N).astype(dtype).reshape(shape)
- H0 = np.empty(shape=shape, dtype=dtype)
- H1 = np.empty(shape=shape, dtype=dtype)
- H0, H1, Href = self.simd_align(H0, H1, Href)
+ if (itype==1): # real size is 2*(N+1)
+ shape = mk_shape(shape, -1, shape[-1] - 1)
+ Href = np.random.rand(*shape).astype(dtype).reshape(shape)
results = {}
- for (name, impl) in self.implementations.iteritems():
- iitype = [1,3,2,4][itype-1]
- if dtype not in impl.supported_ftypes:
- continue
- if itype not in impl.supported_sine_transforms:
- continue
- if iitype not in impl.supported_sine_transforms:
- continue
- forward = impl.dst(a=H0, out=H1, type=itype)
- backward = impl.idst(a=H1, out=H0, type=itype)
- H0[...] = Href
- forward.execute()
- backward.execute()
- results[name] = np.max(np.abs(Href-H0))
+ for (kind, implementations) in self.implementations.iteritems():
+ for (name, impl) in implementations.iteritems():
+ iitype = [1,3,2,4][itype-1]
+ if dtype not in impl.supported_ftypes:
+ continue
+ if itype not in impl.supported_sine_transforms:
+ continue
+ if iitype not in impl.supported_sine_transforms:
+ continue
+ D0 = impl.backend.empty(shape=shape, dtype=dtype)
+ D1 = impl.backend.empty(shape=shape, dtype=dtype)
+ forward = impl.dst(a=D0, out=D1, type=itype).setup()
+ backward = impl.idst(a=D1, out=D0, type=itype).setup()
+ assert forward.input_array is D0
+ assert forward.output_array is D1
+ assert backward.input_array is D1
+ assert backward.output_array is D0
+ D0[...] = Href
+ D1[...] = np.nan
+ forward.execute()
+ D0[...] = np.nan
+ backward.execute()
+ H0 = D0.get()
+ results[name] = np.max(np.abs(Href-H0))
check_distances(results)
def iter_shapes(self):
minj=(2,2)
- maxj=(3,3)
- # maxj=(6,6)
+ maxj=(5,5)
msg = ('EVEN', 'ODD')
for i in xrange(2):
base = 2+i
print ' '+msg[i]
for j1 in xrange(minj[i],maxj[i]):
- shape = (2,4,3,base**j1,)
+ shape = (3,2,base**j1,)
cshape = list(shape)
cshape[-1] = cshape[-1]//2 + 1
cshape = tuple(cshape)
@@ -566,15 +583,6 @@ class TestFFT(object):
ss += (s,)
failed = (Eeps>report_eps)
if failed:
- print
- print 'r0={}, r1={}'.format(r0, r1)
- print
- print results[r0]
- print
- print results[r1]
- print
- import sys
- sys.exit(1)
shape=results[r0].shape
failures.setdefault(tag, []).append((r0, r1, shape, Einf, Eeps))
@@ -585,12 +593,6 @@ class TestFFT(object):
msg='Some implementations did not agree on the result !'
raise RuntimeError(msg)
- def simd_align(self, *arrays):
- aligned_arrays = ()
- for array in arrays:
- aligned_arrays += (pyfftw.byte_align(array),)
- return aligned_arrays
-
def report_failures(self, failures):
print
print '== TEST FAILURES REPORT =='
@@ -635,9 +637,9 @@ class TestFFT(object):
# not testing np.longdouble because only one implementation supports it
# ie. we cannot compare results between different implementations
failures = {}
- for dtype in (np.float64,):
+ for dtype in (np.float32, np.float64,):
self._test_forward_backward_1d(dtype=dtype)
- # self._test_1d(dtype=dtype, failures=failures.setdefault(dtype.__name__, {}))
+ self._test_1d(dtype=dtype, failures=failures.setdefault(dtype.__name__, {}))
self.report_failures(failures)
if __name__ == '__main__':