diff --git a/bfps/cpp/full_code/symmetrize_test.cpp b/bfps/cpp/full_code/symmetrize_test.cpp
index 471e234ced35b11e85b1eea0c87d84542a15c343..13d48d6040aef829d545d364e0565230fc7cedeb 100644
--- a/bfps/cpp/full_code/symmetrize_test.cpp
+++ b/bfps/cpp/full_code/symmetrize_test.cpp
@@ -69,15 +69,12 @@ int symmetrize_test<rnumber>::do_work(void)
                 ptrdiff_t yindex,
                 ptrdiff_t zindex,
                 double k2){
-            if (k2 < kk->kM2)
-            {
-                test_field0->cval(cindex, 0, 0) = rdist(rgen);
-                test_field0->cval(cindex, 0, 1) = rdist(rgen);
-                test_field0->cval(cindex, 1, 0) = rdist(rgen);
-                test_field0->cval(cindex, 1, 1) = rdist(rgen);
-                test_field0->cval(cindex, 2, 0) = rdist(rgen);
-                test_field0->cval(cindex, 2, 1) = rdist(rgen);
-            }
+            test_field0->cval(cindex, 0, 0) = rdist(rgen);
+            test_field0->cval(cindex, 0, 1) = rdist(rgen);
+            test_field0->cval(cindex, 1, 0) = rdist(rgen);
+            test_field0->cval(cindex, 1, 1) = rdist(rgen);
+            test_field0->cval(cindex, 2, 0) = rdist(rgen);
+            test_field0->cval(cindex, 2, 1) = rdist(rgen);
             if (k2 > 0)
             {
                 test_field0->cval(cindex, 0, 0) /= sqrt(k2);
@@ -98,12 +95,11 @@ int symmetrize_test<rnumber>::do_work(void)
             }
             });
     // dealias (?!)
-    //kk->template low_pass<rnumber, THREE>(test_field0->get_cdata(), kk->kM);
     kk->template dealias<rnumber, THREE>(test_field0->get_cdata());
     // make the field divergence free
     kk->template force_divfree<rnumber>(test_field0->get_cdata());
     // apply symmetrize to test_field0
-    //test_field0->symmetrize();
+    test_field0->symmetrize();
 
 
     // make copy in test_field1
diff --git a/bfps/test/test_fftw.py b/bfps/test/test_fftw.py
index 38f82e62a96c686aead162ba8931fb1a95288a28..3de2d97df167567899fbf8b19c1123e5bf35cbe7 100644
--- a/bfps/test/test_fftw.py
+++ b/bfps/test/test_fftw.py
@@ -1,53 +1,64 @@
 #! /usr/bin/env python
 
-import os
 import numpy as np
-#import cupy as np
 import h5py
 import sys
 
 import bfps
 from bfps import TEST
 
-import matplotlib.pyplot as plt
+try:
+    import matplotlib.pyplot as plt
+except:
+    plt = None
 
 def main():
     niterations = 10
-    c = TEST()
-    c.launch(
-            ['symmetrize_test',
-             '--nx', '32',
-             '--ny', '32',
-             '--nz', '32',
-             '--np', '4',
-             '--ntpp', '1',
-             '--wd', './'] +
-             sys.argv[1:])
-
-    df = h5py.File(c.simname + '.h5', 'r')
-    df = h5py.File(c.simname + '_fields.h5', 'r')
-    field1_complex = df['field1/complex/0'].value
-    field1_real = df['field1/real/0'].value
-
-    np_field1_real = np.fft.irfftn(field1_complex, axes = (0, 1, 2)).transpose(1, 0, 2, 3)
-    L2normr = np.sqrt(np.mean(field1_real**2))
-    err = np.max(np.abs(field1_real - np_field1_real*(np_field1_real.size/3)))
-    assert(err < 1e-5)
-
-    np_field1_complex = np.fft.rfftn(field1_real.transpose(1, 0, 2, 3), axes = (0, 1, 2)) / (np_field1_real.size/3)
-
-    L2norm0 = np.sqrt(np.sum(np.abs(field1_complex[:, :, 0])**2) + 2*np.sum(np.abs(field1_complex[:, :, 1:])**2))
-    L2norm1 = np.sqrt(np.sum(np.abs(np_field1_complex[:, :, 0])**2) + 2*np.sum(np.abs(np_field1_complex[:, :, 1:])**2))
-    err = np.max(np.abs(np_field1_complex - field1_complex))
-    assert(err < 1e-5)
-    print(L2normr, L2norm0, L2norm1)
-
-    f = plt.figure()
-    a = f.add_subplot(121)
-    a.imshow(np.log(np.abs(np_field1_complex[:, :, 0, 0])), interpolation = 'nearest')
-    a = f.add_subplot(122)
-    a.imshow(np.log(np.abs(field1_complex[:, :, 0, 0])), interpolation = 'nearest')
-    f.savefig('symmetrize_test.pdf')
+    nlist = [16, 32, 48, 24, 64, 12]
+    for ii in range(len(nlist)):
+        c = TEST()
+        c.launch(
+                ['symmetrize_test',
+                 '--nx', str(nlist[ii]),
+                 '--ny', str(nlist[(ii+1)%(len(nlist))]),
+                 '--nz', str(nlist[(ii+2)%(len(nlist))]),
+                 '--Lx', str(2+np.random.random()),
+                 '--Ly', str(2+np.random.random()),
+                 '--Lz', str(2+np.random.random()),
+                 '--simname', 'fftw_vs_numpy_{0}'.format(ii),
+                 '--np', '4',
+                 '--ntpp', '1',
+                 '--wd', './'] +
+                 sys.argv[1:])
+        df = h5py.File(c.simname + '.h5', 'r')
+        df = h5py.File(c.simname + '_fields.h5', 'r')
+        field1_complex = df['field1/complex/0'].value
+        field1_real = df['field1/real/0'].value
+        npoints = field1_real.size//3
+
+        np_field1_real = np.fft.irfftn(field1_complex, axes = (0, 1, 2)).transpose(1, 0, 2, 3)
+        L2normr = np.sqrt(np.mean(np.sum(field1_real**2, axis = 3)))
+        np_L2normr = np.sqrt(np.mean(np.sum(np_field1_real**2, axis = 3)))
+        err = np.max(np.abs(field1_real - np_field1_real*npoints)) / L2normr
+        assert(err < 1e-5)
+
+        np_field1_complex = np.fft.rfftn(field1_real.transpose(1, 0, 2, 3), axes = (0, 1, 2)) / npoints
+
+        L2norm0 = np.sqrt(np.sum(np.abs(field1_complex[:, :, 0])**2) + 2*np.sum(np.abs(field1_complex[:, :, 1:])**2))
+        L2norm1 = np.sqrt(np.sum(np.abs(np_field1_complex[:, :, 0])**2) + 2*np.sum(np.abs(np_field1_complex[:, :, 1:])**2))
+        err = np.max(np.abs(np_field1_complex - field1_complex)) / L2norm0
+        assert(err < 1e-5)
+
+        err = abs(L2normr - L2norm0) / L2norm0
+        assert(err < 1e-5)
+
+        if not type(plt) == type(None):
+            f = plt.figure()
+            a = f.add_subplot(121)
+            a.imshow(np.log(np.abs(np_field1_complex[:, :, 0, 0])), interpolation = 'nearest')
+            a = f.add_subplot(122)
+            a.imshow(np.log(np.abs(field1_complex[:, :, 0, 0])), interpolation = 'nearest')
+            f.savefig(c.simname + '_complex_slice_kx0.pdf')
     return None
 
 if __name__ == '__main__':
diff --git a/setup.py b/setup.py
index 1cb8ad9387cdd8543c2b3d31d8334d595d998322..9671a3c63cd8526365c0e063ec3493b679b8511b 100644
--- a/setup.py
+++ b/setup.py
@@ -291,7 +291,8 @@ setup(
                 'bfps1 = bfps.__main__:main',
                 'bfps.test_NSVEparticles = bfps.test.test_bfps_NSVEparticles:main',
                 'bfps.test_particles = bfps.test.test_particles:main',
-                'bfps.test_Parseval = bfps.test.test_Parseval:main'],
+                'bfps.test_Parseval = bfps.test.test_Parseval:main',
+                'bfps.test_fftw = bfps.test.test_fftw:main'],
             },
         version = VERSION,
 ########################################################################
diff --git a/tests/ci-scripts/test.sh b/tests/ci-scripts/test.sh
index a2178e44c529bae0c4e4fc821aa55f4966d4bf41..e2adc661efdba6b0c918378d86d73a3c5e5411f4 100644
--- a/tests/ci-scripts/test.sh
+++ b/tests/ci-scripts/test.sh
@@ -36,10 +36,12 @@ $pythonbin setup.py install --prefix=$destdir
 ls $destdir
 ls $destdir/bin/
 
-$pythonbin $destdir/bin/bfps.test_NSVEparticles
+$pythonbin $destdir/bin/bfps.test_fftw
 
 $pythonbin $destdir/bin/bfps.test_Parseval
 
+$pythonbin $destdir/bin/bfps.test_NSVEparticles
+
 # Clean
 if [[ -d $destdir ]] ; then
     rm -rf $destdir ;