diff --git a/hysop/backend/host/fortran/operator/poisson.py b/hysop/backend/host/fortran/operator/poisson.py
index fc33f70e05659a972213a84fe5c7d71232f61e3f..6f5638ff9506973ba2f8af09ad16fa54d695e4d5 100644
--- a/hysop/backend/host/fortran/operator/poisson.py
+++ b/hysop/backend/host/fortran/operator/poisson.py
@@ -4,6 +4,7 @@ from hysop.fields.continuous_field import Field
 from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
 from hysop.backend.host.fortran.operator.fortran_fftw import fftw2py, FortranFFTWOperator
 from hysop.core.graph.graph import op_apply
+from hysop.constants import HYSOP_REAL
 
 
 class PoissonFFTW(FortranFFTWOperator):
@@ -64,11 +65,11 @@ class PoissonFFTW(FortranFFTWOperator):
         if self.discretized:
             return
         super(PoissonFFTW, self).discretize()
-        self.dFin  = self.get_input_discrete_field(self.Fin)
+        self.dFin = self.get_input_discrete_field(self.Fin)
         self.dFout = self.get_output_discrete_field(self.Fout)
         assert (self.dFin.ghosts == self.dFout.ghosts).all(), \
             "Input and output fields must have the same ghosts."
-        self.ghosts = self.dFin.ghosts.astype(np.int64) # prevent f2py copy
+        self.ghosts = self.dFin.ghosts.astype(HYSOP_REAL)  # prevent f2py copy
         self.buffers = self.dFin.buffers + self.dFout.buffers
 
     @op_apply
@@ -79,4 +80,3 @@ class PoissonFFTW(FortranFFTWOperator):
         (buf_in, buf_out) = self.buffers
         self._solve(buf_in, buf_out, self.ghosts)
         self.dFout.exchange_ghosts()
-
diff --git a/hysop/core/arrays/tests/test_array.py b/hysop/core/arrays/tests/test_array.py
index 2520b3435bcd3d58773478a7b27cb79bfcc8f62b..cfe8cfeff084fd8997b14c180ba5837342b6cb13 100644
--- a/hysop/core/arrays/tests/test_array.py
+++ b/hysop/core/arrays/tests/test_array.py
@@ -3,7 +3,7 @@ import warnings
 from contextlib import contextmanager
 from random import randint
 from hysop.testsenv import opencl_failed, iter_clenv, \
-                           __HAS_OPENCL_BACKEND__, __ENABLE_LONG_TESTS__
+    __HAS_OPENCL_BACKEND__, __ENABLE_LONG_TESTS__
 from hysop.tools.contexts import printoptions
 from hysop.tools.numerics import match_float_type, is_unsigned, is_integer, is_complex
 from hysop.tools.types import to_list
@@ -14,6 +14,7 @@ from hysop.core.arrays.all import HostArrayBackend, OpenClArrayBackend, ArrayBac
 from hysop.backend.host.host_allocator import HostAllocator
 from hysop.core.memory.mempool import MemoryPool
 
+
 class TestArray(object):
 
     @classmethod
@@ -23,7 +24,7 @@ class TestArray(object):
     @classmethod
     def teardown_class(cls):
         pass
-    
+
     def setup_method(self, method):
         pass
 
@@ -38,195 +39,204 @@ class TestArray(object):
         d = backend.full_like(c, fill_value=2)
         e = backend.empty_like(d)
         e.fill(3)
-        for buf,val in zip([a,b,c,d,e],[-1,0,1,2,3]):
+        for buf, val in zip([a, b, c, d, e], [-1, 0, 1, 2, 3]):
             assert buf.ndim == 1
-            assert buf.size  == 10
+            assert buf.size == 10
             assert buf.shape == (10,)
             assert buf.itemsize == 2
             assert buf.nbytes == 20
             assert buf.dtype == np.int16
-            assert buf.sum().get()==val*10
-            assert buf.prod()==val**10
+            assert buf.sum().get() == val*10
+            assert buf.prod() == val**10
             assert buf.is_c_contiguous()
             assert buf.is_fortran_contiguous()
             assert buf.is_hysop_contiguous()
-        
-        a = backend.empty(shape=(10,10,), dtype=np.float64,
-                order=MemoryOrdering.C_CONTIGUOUS)
+
+        a = backend.empty(shape=(10, 10,), dtype=np.float64,
+                          order=MemoryOrdering.C_CONTIGUOUS)
         a.fill(15)
         b = backend.zeros_like(a)
         c = backend.ones_like(b)
         d = backend.full_like(c, fill_value=2)
-        for buf,val in zip([a,b,c,d],[15,0,1,2]):
+        for buf, val in zip([a, b, c, d], [15, 0, 1, 2]):
             assert buf.ndim == 2
-            assert buf.size  == 100
-            assert buf.shape == (10,10,)
+            assert buf.size == 100
+            assert buf.shape == (10, 10,)
             assert buf.itemsize == 8
             assert buf.nbytes == 800
             assert buf.dtype == np.float64
-            assert np.allclose(buf.sum().get(),val*100.0)
+            assert np.allclose(buf.sum().get(), val*100.0)
             assert buf.is_c_contiguous()
-       
-        a = backend.empty(shape=(10,10,), dtype=np.float32,
-                order=MemoryOrdering.F_CONTIGUOUS)
+
+        a = backend.empty(shape=(10, 10,), dtype=np.float32,
+                          order=MemoryOrdering.F_CONTIGUOUS)
         a.fill(15)
         b = backend.zeros_like(a)
         c = backend.ones_like(b)
         d = backend.full_like(c, fill_value=2)
-        for buf,val in zip([a,b,c,d],[15,0,1,2]):
+        for buf, val in zip([a, b, c, d], [15, 0, 1, 2]):
             assert buf.ndim == 2
-            assert buf.size  == 100
-            assert buf.shape == (10,10,)
+            assert buf.size == 100
+            assert buf.shape == (10, 10,)
             assert buf.itemsize == 4
             assert buf.nbytes == 400
             assert buf.dtype == np.float32
-            assert buf.sum().get()==val*100
+            assert buf.sum().get() == val*100
             assert buf.is_fortran_contiguous()
-        
+
         a.fill(5)
         b.copy_from(a)
         c.copy_from(b.handle)
         d = c.copy()
-        for buf in [a,b,c,d]:
+        for buf in [a, b, c, d]:
             assert buf.ndim == 2
-            assert buf.size  == 100
-            assert buf.shape == (10,10,)
+            assert buf.size == 100
+            assert buf.shape == (10, 10,)
             assert buf.itemsize == 4
             assert buf.nbytes == 400
             assert buf.dtype == np.float32
-            assert buf.sum().get()==5*100
+            assert buf.sum().get() == 5*100
             assert buf.is_fortran_contiguous()
 
     def _test_transpose_routines(self, backend):
-        
+
         # ensure strides are working as intended
         # (strides are in bytes but itemsize == 1 byte here)
-        A = backend.arange(2*4*8, dtype=np.int8).reshape((2,4,8), order=MemoryOrdering.C_CONTIGUOUS)
-        B = backend.arange(2*4*8, dtype=np.int8).reshape((2,4,8), order=MemoryOrdering.F_CONTIGUOUS)
-        i1,j1,k1 = (randint(0, A.shape[i]-1) for i in xrange(3))
-        i0,j0,k0 = (randint(0, A.shape[i]-1) for i in xrange(3))
-        
-        assert A.dtype.itemsize==1
-        assert A.shape == (2,4,8)
+        A = backend.arange(2*4*8, dtype=np.int8).reshape((2, 4, 8),
+                                                         order=MemoryOrdering.C_CONTIGUOUS)
+        B = backend.arange(2*4*8, dtype=np.int8).reshape((2, 4, 8),
+                                                         order=MemoryOrdering.F_CONTIGUOUS)
+        i1, j1, k1 = (randint(0, A.shape[i]-1) for i in xrange(3))
+        i0, j0, k0 = (randint(0, A.shape[i]-1) for i in xrange(3))
+
+        assert A.dtype.itemsize == 1
+        assert A.shape == (2, 4, 8)
         assert A[0][0][0] == 0
-        assert A[0][0][1] == 1 
+        assert A[0][0][1] == 1
         assert A[0][1][0] == 8
         assert A[1][0][0] == 8*4
-        assert A.strides  == (8*4,8,1)
+        assert A.strides == (8*4, 8, 1)
         assert A[1][1][1] == np.sum(np.asarray(A.strides) / A.dtype.itemsize)
-        assert A[i1][j1][k1] ==  np.sum(np.asarray(A.strides) * (i1,j1,k1))
-        assert A[i0][j0][k0] ==  np.sum(np.asarray(A.strides) * (i0,j0,k0))
-        assert (A[i1][j1][k1]-A[i0][j0][k0]) ==  np.dot(A.strides, (i1-i0,j1-j0,k1-k0))
-        
-        assert B.dtype.itemsize==1
-        assert B.shape == (2,4,8)
+        assert A[i1][j1][k1] == np.sum(np.asarray(A.strides) * (i1, j1, k1))
+        assert A[i0][j0][k0] == np.sum(np.asarray(A.strides) * (i0, j0, k0))
+        assert (A[i1][j1][k1]-A[i0][j0][k0]) == np.dot(A.strides, (i1-i0, j1-j0, k1-k0))
+
+        assert B.dtype.itemsize == 1
+        assert B.shape == (2, 4, 8)
         assert B[0][0][0] == 0
         assert B[0][0][1] == 2*4
         assert B[0][1][0] == 2
         assert B[1][0][0] == 1
-        assert B.strides  == (1,2,2*4)
+        assert B.strides == (1, 2, 2*4)
         assert B[1][1][1] == np.sum(np.asarray(B.strides) / B.dtype.itemsize)
-        assert B[i1][j1][k1] ==  np.sum(np.asarray(B.strides) * (i1,j1,k1))
-        assert B[i0][j0][k0] ==  np.sum(np.asarray(B.strides) * (i0,j0,k0))
-        assert (B[i1][j1][k1]-B[i0][j0][k0]) ==  np.dot(B.strides, (i1-i0,j1-j0,k1-k0))
-        
+        assert B[i1][j1][k1] == np.sum(np.asarray(B.strides) * (i1, j1, k1))
+        assert B[i0][j0][k0] == np.sum(np.asarray(B.strides) * (i0, j0, k0))
+        assert (B[i1][j1][k1]-B[i0][j0][k0]) == np.dot(B.strides, (i1-i0, j1-j0, k1-k0))
+
         # ensure permutations are working as intended
-        A = backend.arange(6, dtype=np.int8).reshape((1,2,3), order=MemoryOrdering.C_CONTIGUOUS)
-        B = backend.arange(6, dtype=np.int8).reshape((1,2,3), order=MemoryOrdering.F_CONTIGUOUS)
-        for arr in (A,B):
-            # all 3d permutations 
-            assert backend.transpose(arr, axes=(0,1,2)).shape == (1,2,3)
-            assert backend.transpose(arr, axes=(0,2,1)).shape == (1,3,2)
-            assert backend.transpose(arr, axes=(1,0,2)).shape == (2,1,3)
-            assert backend.transpose(arr, axes=(2,1,0)).shape == (3,2,1)
-            assert backend.transpose(arr, axes=(2,0,1)).shape == (3,1,2)
-            assert backend.transpose(arr, axes=(1,2,0)).shape == (2,3,1)
-            
+        A = backend.arange(6, dtype=np.int8).reshape((1, 2, 3), order=MemoryOrdering.C_CONTIGUOUS)
+        B = backend.arange(6, dtype=np.int8).reshape((1, 2, 3), order=MemoryOrdering.F_CONTIGUOUS)
+        for arr in (A, B):
+            # all 3d permutations
+            assert backend.transpose(arr, axes=(0, 1, 2)).shape == (1, 2, 3)
+            assert backend.transpose(arr, axes=(0, 2, 1)).shape == (1, 3, 2)
+            assert backend.transpose(arr, axes=(1, 0, 2)).shape == (2, 1, 3)
+            assert backend.transpose(arr, axes=(2, 1, 0)).shape == (3, 2, 1)
+            assert backend.transpose(arr, axes=(2, 0, 1)).shape == (3, 1, 2)
+            assert backend.transpose(arr, axes=(1, 2, 0)).shape == (2, 3, 1)
+
             # transpositions (cycles of length 2)
-            assert backend.transpose(backend.transpose(arr, axes=(0,1,2)), axes=(0,1,2)).shape == (1,2,3)
-            assert backend.transpose(backend.transpose(arr, axes=(0,2,1)), axes=(0,2,1)).shape == (1,2,3)
-            assert backend.transpose(backend.transpose(arr, axes=(1,0,2)), axes=(1,0,2)).shape == (1,2,3)
-            assert backend.transpose(backend.transpose(arr, axes=(2,1,0)), axes=(2,1,0)).shape == (1,2,3)
-            assert backend.transpose(backend.transpose(arr, axes=(2,0,1)), axes=(1,2,0)).shape == (1,2,3)
-            assert backend.transpose(backend.transpose(arr, axes=(1,2,0)), axes=(2,0,1)).shape == (1,2,3)
+            assert backend.transpose(backend.transpose(arr, axes=(0, 1, 2)),
+                                     axes=(0, 1, 2)).shape == (1, 2, 3)
+            assert backend.transpose(backend.transpose(arr, axes=(0, 2, 1)),
+                                     axes=(0, 2, 1)).shape == (1, 2, 3)
+            assert backend.transpose(backend.transpose(arr, axes=(1, 0, 2)),
+                                     axes=(1, 0, 2)).shape == (1, 2, 3)
+            assert backend.transpose(backend.transpose(arr, axes=(2, 1, 0)),
+                                     axes=(2, 1, 0)).shape == (1, 2, 3)
+            assert backend.transpose(backend.transpose(arr, axes=(2, 0, 1)),
+                                     axes=(1, 2, 0)).shape == (1, 2, 3)
+            assert backend.transpose(backend.transpose(arr, axes=(1, 2, 0)),
+                                     axes=(2, 0, 1)).shape == (1, 2, 3)
 
             # cycles of length 3
-            assert backend.transpose(backend.transpose(backend.transpose(arr, axes=(2,0,1)), axes=(2,0,1)), axes=(2,0,1)).shape == (1,2,3)
-            assert backend.transpose(backend.transpose(backend.transpose(arr, axes=(1,2,0)), axes=(1,2,0)), axes=(1,2,0)).shape == (1,2,3)
-            
+            assert backend.transpose(backend.transpose(backend.transpose(
+                arr, axes=(2, 0, 1)), axes=(2, 0, 1)), axes=(2, 0, 1)).shape == (1, 2, 3)
+            assert backend.transpose(backend.transpose(backend.transpose(
+                arr, axes=(1, 2, 0)), axes=(1, 2, 0)), axes=(1, 2, 0)).shape == (1, 2, 3)
+
             # roll, swap and move axes
-            assert backend.rollaxis(arr, axis=0, start=0).shape == (1,2,3)
-            assert backend.rollaxis(arr, axis=1, start=0).shape == (2,1,3)
-            assert backend.rollaxis(arr, axis=2, start=0).shape == (3,1,2)
-            assert backend.rollaxis(arr, axis=0, start=1).shape == (1,2,3)
-            assert backend.rollaxis(arr, axis=1, start=1).shape == (1,2,3)
-            assert backend.rollaxis(arr, axis=2, start=1).shape == (1,3,2)
-            assert backend.rollaxis(arr, axis=0, start=2).shape == (2,1,3)
-            assert backend.rollaxis(arr, axis=1, start=2).shape == (1,2,3)
-            assert backend.rollaxis(arr, axis=2, start=2).shape == (1,2,3)
-            assert backend.rollaxis(arr, axis=0, start=3).shape == (2,3,1)
-            assert backend.rollaxis(arr, axis=1, start=3).shape == (1,3,2)
-            assert backend.rollaxis(arr, axis=2, start=3).shape == (1,2,3)
-
-            assert backend.swapaxes(arr, axis1=0, axis2=0).shape == (1,2,3)
-            assert backend.swapaxes(arr, axis1=1, axis2=0).shape == (2,1,3)
-            assert backend.swapaxes(arr, axis1=2, axis2=0).shape == (3,2,1)
-            assert backend.swapaxes(arr, axis1=0, axis2=1).shape == (2,1,3)
-            assert backend.swapaxes(arr, axis1=1, axis2=1).shape == (1,2,3)
-            assert backend.swapaxes(arr, axis1=2, axis2=1).shape == (1,3,2)
-            assert backend.swapaxes(arr, axis1=0, axis2=2).shape == (3,2,1)
-            assert backend.swapaxes(arr, axis1=1, axis2=2).shape == (1,3,2)
-            assert backend.swapaxes(arr, axis1=2, axis2=2).shape == (1,2,3)
-
-            assert backend.moveaxis(arr, source=0, destination=0).shape == (1,2,3)
-            assert backend.moveaxis(arr, source=1, destination=0).shape == (2,1,3)
-            assert backend.moveaxis(arr, source=2, destination=0).shape == (3,1,2)
-            assert backend.moveaxis(arr, source=0, destination=1).shape == (2,1,3)
-            assert backend.moveaxis(arr, source=1, destination=1).shape == (1,2,3)
-            assert backend.moveaxis(arr, source=2, destination=1).shape == (1,3,2)
-            assert backend.moveaxis(arr, source=0, destination=2).shape == (2,3,1)
-            assert backend.moveaxis(arr, source=1, destination=2).shape == (1,3,2)
-            assert backend.moveaxis(arr, source=2, destination=2).shape == (1,2,3)
-        
+            assert backend.rollaxis(arr, axis=0, start=0).shape == (1, 2, 3)
+            assert backend.rollaxis(arr, axis=1, start=0).shape == (2, 1, 3)
+            assert backend.rollaxis(arr, axis=2, start=0).shape == (3, 1, 2)
+            assert backend.rollaxis(arr, axis=0, start=1).shape == (1, 2, 3)
+            assert backend.rollaxis(arr, axis=1, start=1).shape == (1, 2, 3)
+            assert backend.rollaxis(arr, axis=2, start=1).shape == (1, 3, 2)
+            assert backend.rollaxis(arr, axis=0, start=2).shape == (2, 1, 3)
+            assert backend.rollaxis(arr, axis=1, start=2).shape == (1, 2, 3)
+            assert backend.rollaxis(arr, axis=2, start=2).shape == (1, 2, 3)
+            assert backend.rollaxis(arr, axis=0, start=3).shape == (2, 3, 1)
+            assert backend.rollaxis(arr, axis=1, start=3).shape == (1, 3, 2)
+            assert backend.rollaxis(arr, axis=2, start=3).shape == (1, 2, 3)
+
+            assert backend.swapaxes(arr, axis1=0, axis2=0).shape == (1, 2, 3)
+            assert backend.swapaxes(arr, axis1=1, axis2=0).shape == (2, 1, 3)
+            assert backend.swapaxes(arr, axis1=2, axis2=0).shape == (3, 2, 1)
+            assert backend.swapaxes(arr, axis1=0, axis2=1).shape == (2, 1, 3)
+            assert backend.swapaxes(arr, axis1=1, axis2=1).shape == (1, 2, 3)
+            assert backend.swapaxes(arr, axis1=2, axis2=1).shape == (1, 3, 2)
+            assert backend.swapaxes(arr, axis1=0, axis2=2).shape == (3, 2, 1)
+            assert backend.swapaxes(arr, axis1=1, axis2=2).shape == (1, 3, 2)
+            assert backend.swapaxes(arr, axis1=2, axis2=2).shape == (1, 2, 3)
+
+            assert backend.moveaxis(arr, source=0, destination=0).shape == (1, 2, 3)
+            assert backend.moveaxis(arr, source=1, destination=0).shape == (2, 1, 3)
+            assert backend.moveaxis(arr, source=2, destination=0).shape == (3, 1, 2)
+            assert backend.moveaxis(arr, source=0, destination=1).shape == (2, 1, 3)
+            assert backend.moveaxis(arr, source=1, destination=1).shape == (1, 2, 3)
+            assert backend.moveaxis(arr, source=2, destination=1).shape == (1, 3, 2)
+            assert backend.moveaxis(arr, source=0, destination=2).shape == (2, 3, 1)
+            assert backend.moveaxis(arr, source=1, destination=2).shape == (1, 3, 2)
+            assert backend.moveaxis(arr, source=2, destination=2).shape == (1, 2, 3)
 
     def _test_array_manipulation_routines(self, backend):
         t3d = TranspositionState[3]
-             
-        a = backend.empty(shape=(8,4,2), dtype=np.int8,
-                order=MemoryOrdering.C_CONTIGUOUS)
-        b = backend.empty(shape=(2,4,8), dtype=np.int8, 
-                order=MemoryOrdering.F_CONTIGUOUS)
-        a.copy_from(backend.arange(0,a.size,dtype=a.dtype))
-        b.copy_from(backend.arange(0,b.size,dtype=b.dtype))
-
-        assert a.transposition_state()==t3d.ZYX
-        assert b.transposition_state()==t3d.XYZ
-        assert (a.ravel()==b.ravel(order=MemoryOrdering.C_CONTIGUOUS)).all()
+
+        a = backend.empty(shape=(8, 4, 2), dtype=np.int8,
+                          order=MemoryOrdering.C_CONTIGUOUS)
+        b = backend.empty(shape=(2, 4, 8), dtype=np.int8,
+                          order=MemoryOrdering.F_CONTIGUOUS)
+        a.copy_from(backend.arange(0, a.size, dtype=a.dtype))
+        b.copy_from(backend.arange(0, b.size, dtype=b.dtype))
+
+        assert a.transposition_state() == t3d.ZYX
+        assert b.transposition_state() == t3d.XYZ
+        assert (a.ravel() == b.ravel(order=MemoryOrdering.C_CONTIGUOUS)).all()
 
         a = a.transpose_to_state(t3d.XZY)
         b = b.transpose_to_state(t3d.XZY)
         assert a.order == MemoryOrdering.OUT_OF_ORDER
         assert b.order == MemoryOrdering.OUT_OF_ORDER
-        assert a.transposition_state()==t3d.XZY
-        assert b.transposition_state()==t3d.XZY
+        assert a.transposition_state() == t3d.XZY
+        assert b.transposition_state() == t3d.XZY
 
-        a = a.transpose([2,0,1])
-        b = b.transpose([1,2,0])
-        assert a.transposition_state()==t3d.YXZ
-        assert b.transposition_state()==t3d.ZYX
+        a = a.transpose([2, 0, 1])
+        b = b.transpose([1, 2, 0])
+        assert a.transposition_state() == t3d.YXZ
+        assert b.transposition_state() == t3d.ZYX
 
         a = a.transpose_to_state(t3d.ZYX)
         b = b.transpose_to_state(t3d.XYZ)
-        assert (a.ravel()==b.ravel(order=MemoryOrdering.C_CONTIGUOUS)).all()
+        assert (a.ravel() == b.ravel(order=MemoryOrdering.C_CONTIGUOUS)).all()
 
         a = a.reshape(8*4*2)
         b = a.reshape(8*4*2)
         assert a.order == default_order
         assert b.order == default_order
-        assert (a==b).all()
+        assert (a == b).all()
 
-        a = a.reshape((8,4,2), order=MemoryOrdering.C_CONTIGUOUS)
-        b = b.reshape((8,4,2), order=MemoryOrdering.F_CONTIGUOUS)
+        a = a.reshape((8, 4, 2), order=MemoryOrdering.C_CONTIGUOUS)
+        b = b.reshape((8, 4, 2), order=MemoryOrdering.F_CONTIGUOUS)
         assert a.order == MemoryOrdering.C_CONTIGUOUS
         assert b.order == MemoryOrdering.F_CONTIGUOUS
 
@@ -251,150 +261,149 @@ class TestArray(object):
         assert a0.transposition_state() == t3d.XYZ
         assert a1.transposition_state() == t3d.ZYX
 
-        assert (a0==a).all()
-        assert (a1==a).all()
-        assert (a2==a).all()
+        assert (a0 == a).all()
+        assert (a1 == a).all()
+        assert (a2 == a).all()
 
     def _test_binary_operations(self, backend):
-        a = backend.rand((10,10))
+        a = backend.rand((10, 10))
         b = backend.rint(a)
-        a = backend.rint(backend.rand((10,10))).astype(np.uint8)
-        b = backend.rint(backend.rand((10,10))).astype(np.uint8)
-        c = backend.rint(backend.rand((10,10))).astype(np.uint8)
-        d = backend.rint(backend.rand((10,10))).astype(np.uint8)
-        
+        a = backend.rint(backend.rand((10, 10))).astype(np.uint8)
+        b = backend.rint(backend.rand((10, 10))).astype(np.uint8)
+        c = backend.rint(backend.rand((10, 10))).astype(np.uint8)
+        d = backend.rint(backend.rand((10, 10))).astype(np.uint8)
+
         assert ((~(~a)) == a).all()
-        assert (((a<<2<<3)>>5) == a).all()
-        assert ((a>>1)==0).all()
-        assert ((a|b|c|d)==(d|c|b|a)).all()
-        assert ((a&b&c&d)==(d&c&b&a)).all()
-        assert ((a^b)==(b^a)).all()
-        assert ((~(a|b))==((~a)&(~b))).all()
-        
-        a = backend.rint(10000*backend.rand((10,10))).astype(np.uint64)
-        b = backend.rint(10000*backend.rand((10,10))).astype(np.uint64)
-        c = backend.rint(10000*backend.rand((10,10))).astype(np.uint64)
-        d = backend.rint(10000*backend.rand((10,10))).astype(np.uint64)
-        
+        assert (((a << 2 << 3) >> 5) == a).all()
+        assert ((a >> 1) == 0).all()
+        assert ((a | b | c | d) == (d | c | b | a)).all()
+        assert ((a & b & c & d) == (d & c & b & a)).all()
+        assert ((a ^ b) == (b ^ a)).all()
+        assert ((~(a | b)) == ((~a) & (~b))).all()
+
+        a = backend.rint(10000*backend.rand((10, 10))).astype(np.uint64)
+        b = backend.rint(10000*backend.rand((10, 10))).astype(np.uint64)
+        c = backend.rint(10000*backend.rand((10, 10))).astype(np.uint64)
+        d = backend.rint(10000*backend.rand((10, 10))).astype(np.uint64)
+
         assert ((~(~a)) == a).all()
-        assert (((a<<2<<3)>>5) == a).all()
-        assert ((a|b|c|d)==(d|c|b|a)).all()
-        assert ((a&b&c&d)==(d&c&b&a)).all()
-        assert ((a^b)==(b^a)).all()
-        assert ((~(a|b))==((~a)&(~b))).all()
-    
+        assert (((a << 2 << 3) >> 5) == a).all()
+        assert ((a | b | c | d) == (d | c | b | a)).all()
+        assert ((a & b & c & d) == (d & c & b & a)).all()
+        assert ((a ^ b) == (b ^ a)).all()
+        assert ((~(a | b)) == ((~a) & (~b))).all()
+
     def _test_arithmetic_operations(self, backend):
-        a = backend.rand((10,10)).astype(np.float64).clip(0.1, 0.9)
+        a = backend.rand((10, 10)).astype(np.float64).clip(0.1, 0.9)
         a = (a-0.5)*10
-        
-        b = 10*backend.rand((10,10)).astype(np.float64).clip(0.1, 0.9)
-        c = 10*backend.rand((10,10)).astype(np.float64).clip(0.1, 0.9)
-        d = 10*backend.rand((10,10)).astype(np.float64).clip(0.1, 0.9)
-        
-        
-        assert backend.allclose( 4.0+a, a+4.0 )
-        assert backend.allclose( 4.0-a, -(a-4.0) )
-        assert backend.allclose( 4.0*a, a*4.0 )
-        assert backend.allclose( 4.0/a, 1.0/(a/4.0) )
-        
-        f,i = backend.modf(a)
-        assert backend.allclose( backend.trunc(a),  i )
-        f,i = backend.modf(b)
-        assert backend.allclose( backend.fmod(b,1), f )
-
-        assert backend.allclose( b//1, i )
-        assert backend.allclose( b%1,  f )
-
-        assert backend.allclose( a-b, -(b-a) )
-        assert backend.allclose( a+b-c-d, -c-d+b+a )
-        assert backend.allclose( a*b*c*d, d*c*a*b )
+
+        b = 10*backend.rand((10, 10)).astype(np.float64).clip(0.1, 0.9)
+        c = 10*backend.rand((10, 10)).astype(np.float64).clip(0.1, 0.9)
+        d = 10*backend.rand((10, 10)).astype(np.float64).clip(0.1, 0.9)
+
+        assert backend.allclose(4.0+a, a+4.0)
+        assert backend.allclose(4.0-a, -(a-4.0))
+        assert backend.allclose(4.0*a, a*4.0)
+        assert backend.allclose(4.0/a, 1.0/(a/4.0))
+
+        f, i = backend.modf(a)
+        assert backend.allclose(backend.trunc(a),  i)
+        f, i = backend.modf(b)
+        assert backend.allclose(backend.fmod(b, 1), f)
+
+        assert backend.allclose(b//1, i)
+        assert backend.allclose(b % 1,  f)
+
+        assert backend.allclose(a-b, -(b-a))
+        assert backend.allclose(a+b-c-d, -c-d+b+a)
+        assert backend.allclose(a*b*c*d, d*c*a*b)
         #assert backend.allclose( (a/b)*(c/d), (a*c)/(b*d) )
-        a = a%b 
-        a = a//b 
-        
+        a = a % b
+        a = a//b
+
         a = c.copy()
         assert backend.allclose(c, a)
-        a+=1
-        assert backend.allclose(c+1,a)
-        a-=3
-        assert backend.allclose(c+1-3,a)
-        a*=2
-        assert backend.allclose(2*(c+1-3),a)
-        a/=3
-        assert backend.allclose(2*(c+1-3)/3,a)
-        a//=4
-        assert backend.allclose((2*(c+1-3)/3)//4,a)
-        a%=2
-        assert backend.allclose(((2*(c+1-3)/3)//4)%2,a)
-        
+        a += 1
+        assert backend.allclose(c+1, a)
+        a -= 3
+        assert backend.allclose(c+1-3, a)
+        a *= 2
+        assert backend.allclose(2*(c+1-3), a)
+        a /= 3
+        assert backend.allclose(2*(c+1-3)/3, a)
+        a //= 4
+        assert backend.allclose((2*(c+1-3)/3)//4, a)
+        a %= 2
+        assert backend.allclose(((2*(c+1-3)/3)//4) % 2, a)
+
     def _test_backend_versus_numpy_operations(self, backend):
         npb = backend
-        
+
         atol = [None]
 
         # pollute array with +inf, -inf and NaNs values
         def pollute(arr):
-            mask    = lambda p: (np.random.rand(*arr.shape)<p)
+            def mask(p): return (np.random.rand(*arr.shape) < p)
             arr[mask(0.20)] = -np.inf
             arr[mask(0.20)] = +np.inf
-            arr[mask(0.20)] = np.nan 
+            arr[mask(0.20)] = np.nan
 
-        def allclose(np_array, backend_array, equal_nan=True, atol=atol, 
-                relaxed_precision=False, ignore_mask=None):
-            atol = atol[0] # 1*epsilon
+        def allclose(np_array, backend_array, equal_nan=True, atol=atol,
+                     relaxed_precision=False, ignore_mask=None):
+            atol = atol[0]  # 1*epsilon
             if relaxed_precision:
-                atol=1e-2
+                atol = 1e-2
             if (backend_array is None):
-                msg='Backend returned nothing (got None).'
+                msg = 'Backend returned nothing (got None).'
                 raise ValueError(msg)
             if not np.isscalar(np_array) and not isinstance(np_array, np.ndarray):
-                msg='first arg is not a np.ndarray (got {})'
-                msg=msg.format(np_array)
+                msg = 'first arg is not a np.ndarray (got {})'
+                msg = msg.format(np_array)
                 raise ValueError(msg)
             if isinstance(backend_array, Array):
                 backend_array = backend_array.get().handle
             if (ignore_mask is not None):
-                np_array      = np_array[~ignore_mask]
+                np_array = np_array[~ignore_mask]
                 backend_array = backend_array[~ignore_mask]
             return np.allclose(np_array, backend_array, equal_nan=equal_nan, atol=atol)
-        
+
         unary_ops = [
-                     'reciprocal', 'negative', 'absolute', 'fabs', 
-                     'sin', 'cos', 'arcsin', 'arccos', 'arctan', 'tan',
-                     'degrees', 'radians', 'deg2rad', 'rad2deg', 
-                     'sinh', 'cosh', 'arcsinh', 'arccosh', 'arctanh', 
-                      #'tanh', 'around',  #FIXME
-                     'rint', 'fix', 'floor', 'ceil', 'trunc', 
-                     'exp', 'expm1', 'exp2', 
-                     'log', 'log10', 'log2', 'log1p', 
-                     'signbit', 'reciprocal', 'negative', 
-                     'sqrt', 'cbrt', 'square', 'sign',
-                     'nan_to_num',
-                     'prod', 'sum', 'nanprod', 'nansum',
-                     'cumprod', 'cumsum', 'nancumprod', 'nancumsum',
-                     'isfinite', 'isinf', 'isnan', 'isneginf', 'isposinf',
-                     'real', 'imag', 'angle', 'conj', 'real_if_close'
-                     ]
-        
+            'reciprocal', 'negative', 'absolute', 'fabs',
+            'sin', 'cos', 'arcsin', 'arccos', 'arctan', 'tan',
+            'degrees', 'radians', 'deg2rad', 'rad2deg',
+            'sinh', 'cosh', 'arcsinh', 'arccosh', 'arctanh',
+            # 'tanh', 'around',  #FIXME
+            'rint', 'fix', 'floor', 'ceil', 'trunc',
+            'exp', 'expm1', 'exp2',
+            'log', 'log10', 'log2', 'log1p',
+            'signbit', 'reciprocal', 'negative',
+            'sqrt', 'cbrt', 'square', 'sign',
+            'nan_to_num',
+            'prod', 'sum', 'nanprod', 'nansum',
+            'cumprod', 'cumsum', 'nancumprod', 'nancumsum',
+            'isfinite', 'isinf', 'isnan', 'isneginf', 'isposinf',
+            'real', 'imag', 'angle', 'conj', 'real_if_close'
+        ]
+
         binary_ops = [
-                'minimum', 'maximum', 'fmin', 'fmax', 
-                'add', 'subtract', 'multiply', 'power', 
-                'divide', 'floor_divide', 'true_divide',
-                'equal', 'not_equal', 'less_equal', 'greater_equal', 'less', 'greater',
-                'mod', 'fmod', 'remainder', 'hypot', 
-                'arctan2',
-                'logaddexp', 'logaddexp2',
-                'copysign'
-            ]
-        
+            'minimum', 'maximum', 'fmin', 'fmax',
+            'add', 'subtract', 'multiply', 'power',
+            'divide', 'floor_divide', 'true_divide',
+            'equal', 'not_equal', 'less_equal', 'greater_equal', 'less', 'greater',
+            'mod', 'fmod', 'remainder', 'hypot',
+            'arctan2',
+            'logaddexp', 'logaddexp2',
+            'copysign'
+        ]
+
         array_unary_ops = ['__neg__', '__abs__']
 
         array_binary_ops = [
-                '__eq__', '__ne__', '__le__', '__ge__', '__lt__', '__gt__', 
-                '__add__', '__sub__', '__mul__', '__pow__', 
-                '__floordiv__', '__div__', '__mod__', 
-                '__radd__', '__rsub__', '__rmul__', '__rpow__', 
-                '__rfloordiv__', '__rdiv__', '__rmod__'
+            '__eq__', '__ne__', '__le__', '__ge__', '__lt__', '__gt__',
+            '__add__', '__sub__', '__mul__', '__pow__',
+            '__floordiv__', '__div__', '__mod__',
+            '__radd__', '__rsub__', '__rmul__', '__rpow__',
+            '__rfloordiv__', '__rdiv__', '__rmod__'
         ]
 
         # real_skip_list = ['angle', 'real', 'imag', 'conj', 'real_if_close']
@@ -405,83 +414,82 @@ class TestArray(object):
                              'logaddexp', 'logaddexp2', 'copysign', 'frexp',
                              '__mod__', '__rmod__']
 
-        splitting_ops = [ 'frexp', 'modf' ]
+        splitting_ops = ['frexp', 'modf']
 
-        def ignore_infty(ref_out,backend_out,**kargs):
-            mask  = np.isinf(ref_out)
+        def ignore_infty(ref_out, backend_out, **kargs):
+            mask = np.isinf(ref_out)
             mask |= np.isinf(backend_out)
             return mask
 
         def positive_int_rhs(variables):
-            assert 'b' in variables # b is rhs
+            assert 'b' in variables  # b is rhs
             rhs = variables['b']
             dtype = rhs[0].dtype
             if is_integer(dtype):
-                for i,v in enumerate(rhs):
+                for i, v in enumerate(rhs):
                     rhs[i] = abs(v)
 
-        def clamp(_amin,_amax):
+        def clamp(_amin, _amax):
             def _filter(variables):
-                for k,_vars in variables.iteritems():
-                    for i,var in enumerate(_vars):
+                for k, _vars in variables.iteritems():
+                    for i, var in enumerate(_vars):
                         if is_complex(var):
-                            if isinstance(var,np.ndarray):
-                                np.clip(var.real, _amin,_amax, variables[k][i].real)
-                                np.clip(var.imag, _amin,_amax, variables[k][i].imag)
+                            if isinstance(var, np.ndarray):
+                                np.clip(var.real, _amin, _amax, variables[k][i].real)
+                                np.clip(var.imag, _amin, _amax, variables[k][i].imag)
                             else:
                                 amin = _amin + 1j*_amin
                                 amax = _amax + 1j*_amax
-                                var.backend.clip_components(var,amin,amax,variables[k][i])
+                                var.backend.clip_components(var, amin, amax, variables[k][i])
                         else:
-                            if isinstance(var,np.ndarray):
-                                np.clip(var.real,_amin,_amax,variables[k][i])
+                            if isinstance(var, np.ndarray):
+                                np.clip(var.real, _amin, _amax, variables[k][i])
                             else:
-                                var.backend.clip(var,_amin,_amax,variables[k][i])
+                                var.backend.clip(var, _amin, _amax, variables[k][i])
             return _filter
-                        
-        
+
         pow_constraints = [positive_int_rhs]
-        pow_constraints.append( clamp(+0,+3) )
+        pow_constraints.append(clamp(+0, +3))
 
-        ## Extra contraints on inputs 
+        # Extra contraints on inputs
         # should be a list of functions taking variables as inputsq
         input_constraints = {
             'power':      pow_constraints,
             '__pow__':    pow_constraints,
             '__rpow__':   pow_constraints,
-            'cumprod':    clamp(0.1,1.1),
-            'nancumprod': clamp(0.1,1.1)
+            'cumprod':    clamp(0.1, 1.1),
+            'nancumprod': clamp(0.1, 1.1)
         }
-        
-        ## Extra contraints on outputs
-        # Generate a mask of values thats should not 
+
+        # Extra contraints on outputs
+        # Generate a mask of values thats should not
         # be compared to numpy solution in allclose check)
         #  all keys are operator names
-        #  all values are function of dtype and backend, 
+        #  all values are function of dtype and backend,
         output_constraints = {
             'cumprod':    [ignore_infty],
             'nancumprod': [ignore_infty]
         }
-        
+
         class TestContext(object):
             def __init__(self, opname, input_constraints, variables):
-                self.opname  = opname
-                
+                self.opname = opname
+
                 # if there is a specific constraint we copy everything
                 dtypes = {}
                 if opname in input_constraints:
-                    for vname,vargs in variables.iteritems():
-                        for i,var in enumerate(vargs):
+                    for vname, vargs in variables.iteritems():
+                        for i, var in enumerate(vargs):
                             variables[vname][i] = variables[vname][i].copy()
                     filters = to_list(input_constraints[opname])
                     for f in filters:
                         f(variables)
                 self.dtypes = dtypes
 
-                for vname,vargs in variables.iteritems():
+                for vname, vargs in variables.iteritems():
                     dtypes[vname] = variables[vname][0].dtype
-                    for i,var in enumerate(vargs):
-                        varname='{}{}'.format(vname,i)
+                    for i, var in enumerate(vargs):
+                        varname = '{}{}'.format(vname, i)
                         setattr(self, varname, var)
                     setattr(self, vname, vargs)
 
@@ -490,15 +498,14 @@ class TestArray(object):
 
             def __exit__(self, exception, e, traceback):
                 if (e is not None):
-                    msg='\nTESTING: Test failed in at {}::{}() with dtypes {}\n'
-                    msg=msg.format(backend.__class__.__name__,self.opname,self.dtypes)
+                    msg = '\nTESTING: Test failed in at {}::{}() with dtypes {}\n'
+                    msg = msg.format(backend.__class__.__name__, self.opname, self.dtypes)
                     print msg
                     raise exception, e, traceback
 
-        
         def check_inputs(name, _in):
-            isclose = np.isclose(_in[0],_in[1].get(handle=True), equal_nan=True)
-            if not isclose.all(): 
+            isclose = np.isclose(_in[0], _in[1].get(handle=True), equal_nan=True)
+            if not isclose.all():
                 print '{} inputs mismatch...'.format(name)
                 print '{} NUMPY INPUT:'.format(name.upper())
                 print _in[0][~isclose]
@@ -510,72 +517,72 @@ class TestArray(object):
             if (r0 is None) and (r1 is None):
                 return
             elif (r0 is None):
-                msg='numpy::{} returned None.'.format(opname)
+                msg = 'numpy::{} returned None.'.format(opname)
                 raise TypeError(msg)
             elif (r1 is None):
-                msg='{}::{} returned None.'.format(backend.__class__.__name__,opname)
+                msg = '{}::{} returned None.'.format(backend.__class__.__name__, opname)
                 raise TypeError(msg)
             else:
                 if isinstance(r1, Array):
                     r1 = r1.get(handle=True)
-                    if isinstance(backend,OpenClArrayBackend):
+                    if isinstance(backend, OpenClArrayBackend):
                         # FIXME OpenCl support for float16
-                        if r0.dtype==np.float16:
+                        if r0.dtype == np.float16:
                             r1 = r1.astype(np.float16)
-                    if r0.dtype==np.bool:
+                    if r0.dtype == np.bool:
                         r1 = r1.astype(np.bool)
-                
+
                 if (r0.dtype == np.bool) and (r1.dtype == np.bool):
-                    l2   = np.sqrt(np.nansum(r0^r1)) / r0.size
-                    linf = np.nanmax(r0^r1)
+                    l2 = np.sqrt(np.nansum(r0 ^ r1)) / r0.size
+                    linf = np.nanmax(r0 ^ r1)
                 else:
                     m0 = np.isfinite(r0)
                     m1 = np.isfinite(r1)
-                    if (m0!=m1).any():
-                        l2   = np.inf
+                    if (m0 != m1).any():
+                        l2 = np.inf
                         linf = np.inf
                     else:
                         try:
                             R0, R1 = r0[m0], r1[m0]
-                            l2   = np.sqrt(np.sum((R0-R1)*np.conj(R0-R1)) / R0.size)
+                            l2 = np.sqrt(np.sum((R0-R1)*np.conj(R0-R1)) / R0.size)
                             linf = np.max(np.abs(R0-R1))
                         except ValueError:
-                            l2   = 0
+                            l2 = 0
                             linf = 0
-                msg1='(l2={}, linf={}).'
-                msg1=msg1.format(l2, linf)
-                
+                msg1 = '(l2={}, linf={}).'
+                msg1 = msg1.format(l2, linf)
+
                 if (r0.dtype == r1.dtype):
-                    mask=None
+                    mask = None
                     if opname in output_constraints:
                         mask_generators = to_list(output_constraints[opname])
-                        mask = mask_generators[0](ref_out=r0,backend_out=r1)
+                        mask = mask_generators[0](ref_out=r0, backend_out=r1)
                         for mask_gen in mask_generators[1:]:
-                            mask |= mask_gen(ref_out=r0,backend_out=r1)
+                            mask |= mask_gen(ref_out=r0, backend_out=r1)
 
-                    close = allclose(r0,r1,ignore_mask=mask)
+                    close = allclose(r0, r1, ignore_mask=mask)
                     tol = atol[0]
                     if not close:
-                        close = allclose(r0,r1,relaxed_precision=True,ignore_mask=mask)
+                        close = allclose(r0, r1, relaxed_precision=True, ignore_mask=mask)
                         tol = 1e-2
                         if close:
-                            msg='WARNING: test passed with relaxed precision for {}::{}.'
-                            msg=msg.format(backend.__class__.__name__, opname)
+                            msg = 'WARNING: test passed with relaxed precision for {}::{}.'
+                            msg = msg.format(backend.__class__.__name__, opname)
                             print msg
                     if not close:
                         msg = '\n{}::{} returned dtypes did match (got {}) '
-                        msg+= 'but failed to match numpy output,'
-                        msg+='\n absolute tolerance was set to {}.'
-                        msg=msg.format(backend.__class__.__name__,opname,r1.dtype,tol)
+                        msg += 'but failed to match numpy output,'
+                        msg += '\n absolute tolerance was set to {}.'
+                        msg = msg.format(backend.__class__.__name__, opname, r1.dtype, tol)
                         print msg
                         if isinstance(r0, np.ndarray) and isinstance(r1, np.ndarray):
-                            failed=(~np.isclose(r0,r1,equal_nan=True,atol=atol[0]))
+                            failed = (~np.isclose(r0, r1, equal_nan=True, atol=atol[0]))
                             if (lhs is not None):
-                                check_inputs('lhs',lhs)
+                                check_inputs('lhs', lhs)
                                 print 'LHS_INPUT'
                                 print lhs[0][failed]
                             if (rhs is not None):
-                                check_inputs('rhs',rhs)
+                                check_inputs('rhs', rhs)
                                 print 'RHS INPUT'
                                 print rhs[0][failed]
                             print 'EXPECTED'
@@ -585,28 +592,28 @@ class TestArray(object):
                         else:
                             print 'r0 => {}'.format(r0.__class__)
                             print 'r1 => {}'.format(r1.__class__)
-                        msg0='Method {}::{} failed to match numpy output'
-                        msg0=msg0.format(backend.__class__.__name__, opname)
-                        msg=msg0+msg1
+                        msg0 = 'Method {}::{} failed to match numpy output'
+                        msg0 = msg0.format(backend.__class__.__name__, opname)
+                        msg = msg0+msg1
                         print
                         print msg
                         raise ValueError(msg)
                     else:
-                        msg0='{}::{} matched numpy output '
-                        msg0=msg0.format(backend.__class__.__name__, opname)
-                        msg=msg0+msg1
+                        msg0 = '{}::{} matched numpy output '
+                        msg0 = msg0.format(backend.__class__.__name__, opname)
+                        msg = msg0+msg1
                         print msg
                 else:
                     msg = '\n{}::{} returned dtypes didn\'t match (expected {} but got {}).'
                     msg = msg.format(backend.__class__.__name__, opname, r0.dtype, r1.dtype)
                     print msg
 
-                    msg='{}::{} returned dtypes did not match, '\
+                    msg = '{}::{} returned dtypes did not match, '\
                         'got {} but numpy returned {}.'
-                    msg=msg.format(backend.__class__.__name__, opname, r1.dtype, r0.dtype)
+                    msg = msg.format(backend.__class__.__name__, opname, r1.dtype, r0.dtype)
                     raise ValueError(msg)
 
-        def test_operators(a,b,A,B,skip=[]):
+        def test_operators(a, b, A, B, skip=[]):
             with warnings.catch_warnings():
                 warnings.simplefilter('ignore')
 
@@ -615,52 +622,52 @@ class TestArray(object):
                         continue
                     f0 = getattr(np,  opname)
                     f1 = getattr(npb, opname)
-                    with TestContext(opname, input_constraints, 
-                                        variables={'a':[a,A]}) as ctx:
+                    with TestContext(opname, input_constraints,
+                                     variables={'a': [a, A]}) as ctx:
                         r0 = f0(ctx.a0)
                         r1 = f1(ctx.a1)
-                        check_close(ctx.a,None,r0,r1,opname)
+                        check_close(ctx.a, None, r0, r1, opname)
 
                 for opname in binary_ops:
                     if opname in skip:
                         continue
                     f0 = getattr(np,  opname)
                     f1 = getattr(npb, opname)
-                    with TestContext(opname, input_constraints, 
-                            variables={'a':[a,A],'b':[b,B]}) as ctx:
-                        r0 = f0(ctx.a0,ctx.b0)
-                        r1 = f1(ctx.a1,ctx.b1)
-                        check_close(ctx.a,ctx.b,r0,r1,opname)
-                
+                    with TestContext(opname, input_constraints,
+                                     variables={'a': [a, A], 'b': [b, B]}) as ctx:
+                        r0 = f0(ctx.a0, ctx.b0)
+                        r1 = f1(ctx.a1, ctx.b1)
+                        check_close(ctx.a, ctx.b, r0, r1, opname)
+
                 for opname in splitting_ops:
                     if opname in skip:
                         continue
-                    with TestContext(opname, input_constraints, 
-                            variables={'a':[a,A],'b':[b,B]}) as ctx:
+                    with TestContext(opname, input_constraints,
+                                     variables={'a': [a, A], 'b': [b, B]}) as ctx:
                         f0 = getattr(np,  opname)
                         f1 = getattr(npb, opname)
                         r00, r01 = f0(ctx.a0)
                         r10, r11 = f1(ctx.a1)
-                        check_close(ctx.a,None,r00,r10,opname)
-                        check_close(ctx.a,None,r01,r11,opname)
-                
+                        check_close(ctx.a, None, r00, r10, opname)
+                        check_close(ctx.a, None, r01, r11, opname)
+
                 for opname in array_unary_ops:
                     if opname in skip:
                         continue
-                    with TestContext(opname, input_constraints, 
-                                        variables={'a':[a,A]}) as ctx:
+                    with TestContext(opname, input_constraints,
+                                     variables={'a': [a, A]}) as ctx:
                         f0 = getattr(ctx.a0, opname)
                         f1 = getattr(ctx.a1, opname)
                         r0 = f0()
                         r1 = f1()
-                        check_close(ctx.a,None,r0,r1,opname)
+                        check_close(ctx.a, None, r0, r1, opname)
 
                 for opname in array_binary_ops:
                     if opname in skip:
                         continue
-                    with TestContext(opname, input_constraints, 
-                            variables={'a':[a,A],'b':[b,B]}) as ctx:
-                        if opname.find('__r')==0:
+                    with TestContext(opname, input_constraints,
+                                     variables={'a': [a, A], 'b': [b, B]}) as ctx:
+                        if opname.find('__r') == 0:
                             f0 = getattr(ctx.b0, opname)
                             f1 = getattr(ctx.b1, opname)
                             r0 = f0(ctx.a0)
@@ -670,9 +677,8 @@ class TestArray(object):
                             f1 = getattr(ctx.a1, opname)
                             r0 = f0(ctx.b0)
                             r1 = f1(ctx.b1)
-                        check_close(ctx.a,ctx.b,r0,r1,opname)
+                        check_close(ctx.a, ctx.b, r0, r1, opname)
 
-        
         def make_arrays(dtype):
             ftype = match_float_type(dtype)
             atol[0] = np.finfo(ftype).eps
@@ -683,93 +689,92 @@ class TestArray(object):
             if is_unsigned(dtype):
                 a = abs(a)
                 b = abs(b)
-            a = a.astype(dtype) # <= negative number to unsigned dtype conversion wraps
+            a = a.astype(dtype)  # <= negative number to unsigned dtype conversion wraps
             b = b.astype(dtype)
             if is_complex(dtype):
-                a+= (np.random.rand(8192)-0.5)*100j
-                b+= (np.random.rand(8192)-0.5)*100j
-            
-            A, B = npb.asarray(a), npb.asarray(b) 
-            assert allclose( a, A )
-            assert allclose( b, B )
-            assert npb.allclose( npb.asarray(a), A, equal_nan=True )
-            assert npb.allclose( npb.asarray(b), B, equal_nan=True )
-
-            return a,b,A,B
-
-        #FIXME numpy quad float support (gcc __float128), not implemented yet as 
+                a += (np.random.rand(8192)-0.5)*100j
+                b += (np.random.rand(8192)-0.5)*100j
+
+            A, B = npb.asarray(a), npb.asarray(b)
+            assert allclose(a, A)
+            assert allclose(b, B)
+            assert npb.allclose(npb.asarray(a), A, equal_nan=True)
+            assert npb.allclose(npb.asarray(b), B, equal_nan=True)
+
+            return a, b, A, B
+
+        # FIXME numpy quad float support (gcc __float128), not implemented yet as
         if __ENABLE_LONG_TESTS__:
-            signed_types   = (np.int8, np.int16, np.int32, np.int64,)
+            signed_types = (np.int8, np.int16, np.int32, np.int64,)
             unsigned_types = (np.uint8, np.uint16, np.uint32, np.uint64,)
-            float_types    = (np.float16,np.float32,np.float64, np.longdouble)
-            complex_types  = (np.complex64, np.complex128, np.clongdouble)
+            float_types = (np.float16, np.float32, np.float64, np.longdouble)
+            complex_types = (np.complex64, np.complex128, np.clongdouble)
         else:
-            signed_types   = ()
+            signed_types = ()
             unsigned_types = ()
-            float_types    = (np.float32,)
-            complex_types  = (np.complex64,)
-        
+            float_types = (np.float32,)
+            complex_types = (np.complex64,)
+
         for dtype in signed_types:
             print '\n== SIGNED INTEGER OPS {} =='.format(dtype)
-            a,b,A,B = make_arrays(dtype)
-            test_operators(a,b,A,B)
-        
+            a, b, A, B = make_arrays(dtype)
+            test_operators(a, b, A, B)
+
         for dtype in unsigned_types:
             print '\n== UNSIGNED INTEGER OPS {} =='.format(dtype)
-            a,b,A,B = make_arrays(dtype)
-            test_operators(a,b,A,B)
-        
+            a, b, A, B = make_arrays(dtype)
+            test_operators(a, b, A, B)
+
         # FIXME OpenCl backend half float and long double support
         for dtype in float_types:
             print '\n== FLOAT OPS {} =='.format(dtype)
-            if isinstance(backend,OpenClArrayBackend) and (dtype in [np.float16,np.longdouble]):
+            if isinstance(backend, OpenClArrayBackend) and (dtype in [np.float16, np.longdouble]):
                 print '  -- NO SUPPORT PROVIDED BY BACKEND --'
-                continue 
-            
-            a,b,A,B = make_arrays(dtype)
-            test_operators(a,b,A,B)
-        
+                continue
+
+            a, b, A, B = make_arrays(dtype)
+            test_operators(a, b, A, B)
+
             print '\n== POLLUTED FLOAT OPS {} =='.format(dtype)
             pollute(a)
             pollute(b)
 
-            A, B = npb.asarray(a), npb.asarray(b) 
-            test_operators(a,b,A,B)
-        
-        #FIXME OpenCL complex functions: arcsin, arccos, floordix, pow, ...
+            A, B = npb.asarray(a), npb.asarray(b)
+            test_operators(a, b, A, B)
+
+        # FIXME OpenCL complex functions: arcsin, arccos, floordix, pow, ...
         for dtype in complex_types:
             print '\n== COMPLEX OPS {} =='.format(dtype)
-            if isinstance(backend,OpenClArrayBackend):
+            if isinstance(backend, OpenClArrayBackend):
                 if dtype in [np.clongdouble]:
                     print '  -- NO SUPPORT PROVIDED BY BACKEND --'
-                    continue 
+                    continue
 
                 skip_list = [x for x in complex_skip_list]
                 skip_list += ['arcsin',  'arccos', 'arctan',
                               'arcsinh', 'arccosh', 'arctanh',
-                              'exp2', 'expm1', 
+                              'exp2', 'expm1',
                               'log2', 'log10', 'log1p',
                               'floor_divide', '__floordiv__', '__rfloordiv__']
             else:
                 skip_list = complex_skip_list
-        
-            a,b,A,B = make_arrays(dtype)
-            test_operators(a,b,A,B, skip=skip_list)
-            
+
+            a, b, A, B = make_arrays(dtype)
+            test_operators(a, b, A, B, skip=skip_list)
+
             print '\n== POLLUTED COMPLEX OPS {} =='.format(dtype)
             pollute(a)
             pollute(b)
-            
-            if isinstance(backend,OpenClArrayBackend):
-                skip_list+=['power', '__rpow__', '__pow__']
-            
-            A, B = npb.asarray(a), npb.asarray(b) 
-            test_operators(a,b,A,B, skip=skip_list)
+
+            if isinstance(backend, OpenClArrayBackend):
+                skip_list += ['power', '__rpow__', '__pow__']
+
+            A, B = npb.asarray(a), npb.asarray(b)
+            test_operators(a, b, A, B, skip=skip_list)
 
     def _test_backend_versus_numpy(self, backend):
         self._test_backend_versus_numpy_operations(backend)
 
-
     def _test_backend(self, backend):
         with printoptions(linewidth=240, edgeitems=4, threshold=20):
             # self._test_array_creation_routines(backend)
@@ -777,26 +782,26 @@ class TestArray(object):
             # self._test_binary_operations(backend)
             # self._test_arithmetic_operations(backend)
             self._test_backend_versus_numpy(backend)
-            #self._test_array_manipulation_routines(backend)
-    
+            # self._test_array_manipulation_routines(backend)
+
     def test_host_array_backend_allocator(self):
         allocator = HostAllocator()
         backend = HostArrayBackend(allocator=allocator)
         self._test_backend(backend)
 
-    def test_host_array_backend_mempool(self):
-        allocator = HostAllocator()
-        pool    = allocator.memory_pool(name='host')
-        backend = HostArrayBackend(allocator=pool)
+    # def test_host_array_backend_mempool(self):
+    #     allocator = HostAllocator()
+    #     pool = allocator.memory_pool(name='host')
+    #     backend = HostArrayBackend(allocator=pool)
+
+    #     self._test_backend(backend)
+
+    #     backend.allocator.print_allocation_report()
+    #     assert backend.allocator.active_blocks == 0
+    #     backend.allocator.stop_holding()
+    #     assert backend.allocator.held_blocks == 0
+    #     backend.allocator.print_allocation_report()
 
-        self._test_backend(backend)
-        
-        backend.allocator.print_allocation_report()
-        assert backend.allocator.active_blocks == 0
-        backend.allocator.stop_holding()
-        assert backend.allocator.held_blocks == 0
-        backend.allocator.print_allocation_report()
-    
     @opencl_failed
     def test_opencl_array_backend_allocator(self):
         from hysop.backend.device.opencl.opencl_allocator import OpenClImmediateAllocator
@@ -806,27 +811,28 @@ class TestArray(object):
             allocator = OpenClImmediateAllocator(queue=cl_env.default_queue)
             backend = OpenClArrayBackend(cl_env=cl_env, allocator=allocator)
             self._test_backend(backend)
-    
-    @opencl_failed
-    def test_opencl_array_backend_pool(self):
-        from hysop.backend.device.opencl.opencl_allocator import OpenClImmediateAllocator
-        for cl_env in iter_clenv():
-            allocator = OpenClImmediateAllocator(queue=cl_env.default_queue)\
-                                        .memory_pool(name=cl_env.device.name)
-            backend = OpenClArrayBackend(cl_env=cl_env, allocator=allocator)
-            
-            self._test_backend(backend)
-            
-            backend.allocator.print_allocation_report()
-            assert backend.allocator.active_blocks == 0
-            backend.allocator.stop_holding()
-            assert backend.allocator.held_blocks == 0
-            backend.allocator.print_allocation_report()
+
+    # @opencl_failed
+    # def test_opencl_array_backend_pool(self):
+    #     from hysop.backend.device.opencl.opencl_allocator import OpenClImmediateAllocator
+    #     for cl_env in iter_clenv():
+    #         allocator = OpenClImmediateAllocator(queue=cl_env.default_queue)\
+    #                                     .memory_pool(name=cl_env.device.name)
+    #         backend = OpenClArrayBackend(cl_env=cl_env, allocator=allocator)
+
+    #         self._test_backend(backend)
+
+    #         backend.allocator.print_allocation_report()
+    #         assert backend.allocator.active_blocks == 0
+    #         backend.allocator.stop_holding()
+    #         assert backend.allocator.held_blocks == 0
+    #         backend.allocator.print_allocation_report()
+
 
 if __name__ == '__main__':
     test = TestArray()
     test.test_host_array_backend_allocator()
-    #test.test_host_array_backend_mempool()
+    # test.test_host_array_backend_mempool()
     if __HAS_OPENCL_BACKEND__:
         test.test_opencl_array_backend_allocator()
-        #test.test_opencl_array_backend_pool()
+        # test.test_opencl_array_backend_pool()
diff --git a/hysop/core/memory/tests/test_mempool.py b/hysop/core/memory/tests/test_mempool.py
index ac43f5ee301aef79d4979afb4663b40d2c212704..2646442f55e4875757aa84f569eeaf987435413c 100644
--- a/hysop/core/memory/tests/test_mempool.py
+++ b/hysop/core/memory/tests/test_mempool.py
@@ -1,31 +1,37 @@
 from hysop.deps import np
 from hysop.testsenv import opencl_failed, iter_clenv, \
-                           __HAS_OPENCL_BACKEND__, __ENABLE_LONG_TESTS__
+    __HAS_OPENCL_BACKEND__, __ENABLE_LONG_TESTS__
 from hysop.core.memory.mempool import MemoryPool
 
 import random
 max_bytes_per_alloc = 1024*1024*128  # 128MB
-free   = lambda: bool(random.random()>0.1) # 80% probability of free
-nbytes = lambda: int(2.0**(np.log2(max_bytes_per_alloc)*random.random()))
+
+
+def free(): return bool(random.random() > 0.1)  # 80% probability of free
+
+
+def nbytes(): return int(2.0**(np.log2(max_bytes_per_alloc)*random.random()))
+
 
 def test_mempool_python_allocator():
     from hysop.backend.host.host_allocator import HostAllocator
     allocator = HostAllocator()
     _test_mempool_allocator('python', allocator)
 
-@opencl_failed
-def test_mempool_opencl_immediate_allocator():
-    from hysop.backend.device.opencl.opencl_allocator import OpenClImmediateAllocator
+# @opencl_failed
+# def test_mempool_opencl_immediate_allocator():
+#     from hysop.backend.device.opencl.opencl_allocator import OpenClImmediateAllocator
+
+#     for cl_env in iter_clenv():
+#         allocator = OpenClImmediateAllocator(queue=cl_env.default_queue)
+#         _test_mempool_allocator(cl_env.platform.name, allocator)
 
-    for cl_env in iter_clenv():
-        allocator = OpenClImmediateAllocator(queue=cl_env.default_queue)
-        _test_mempool_allocator(cl_env.platform.name, allocator)
 
 def _test_mempool_allocator(name, allocator):
     pool = allocator.memory_pool(name=name)
     buffers = []
     try:
-        while True:
+        for _ in range(10000):
             size = nbytes()
             buf = pool.allocate(size)
             if not free():
@@ -33,7 +39,8 @@ def _test_mempool_allocator(name, allocator):
     except MemoryError:
         pass
 
+
 if __name__ == '__main__':
     test_mempool_python_allocator()
-    if __HAS_OPENCL_BACKEND__:
-        test_mempool_opencl_immediate_allocator()
+    # if __HAS_OPENCL_BACKEND__:
+    #     test_mempool_opencl_immediate_allocator()
diff --git a/hysop/fields/tests/test_fields.py b/hysop/fields/tests/test_fields.py
index 1f1785baed5402bebfc58e80060c89d77779a6c9..ed7a1deb8cb17c29c035f28868b651d3e5176518 100644
--- a/hysop/fields/tests/test_fields.py
+++ b/hysop/fields/tests/test_fields.py
@@ -2,16 +2,17 @@ import numpy as np
 
 from hysop import Box, CartesianTopology, CartesianDiscretization
 from hysop.constants import HYSOP_REAL
-from hysop.fields.continuous_field import Field, ScalarField,TensorField
+from hysop.fields.continuous_field import Field, ScalarField, TensorField
 from hysop.fields.discrete_field import DiscreteField,       \
-                                        DiscreteTensorField, \
-                                        DiscreteScalarFieldView
+    DiscreteTensorField, \
+    DiscreteScalarFieldView
 from hysop.fields.cartesian_discrete_field import CartesianDiscreteField, \
-                                                  CartesianDiscreteScalarFieldView
+    CartesianDiscreteScalarFieldView
 from hysop.defaults import VelocityField, VorticityField
 from hysop.tools.types import check_instance
 from hysop.testsenv import domain_boundary_iterator
 
+
 def test_field():
     domain = Box(dim=3)
     F0 = Field(domain, 'F0')
@@ -23,7 +24,7 @@ def test_field():
     F6 = F0.gradient()
     F7 = F5.field_like('F7')
 
-    D0 = CartesianDiscretization(resolution=(5,5,5), default_boundaries=True)
+    D0 = CartesianDiscretization(resolution=(5, 5, 5), default_boundaries=True)
     T0 = CartesianTopology(domain=domain, discretization=D0)
 
     DF0 = F0.discretize(T0)
@@ -35,10 +36,10 @@ def test_field():
     DF6 = F6.discretize(T0)
     DF7 = F7.discretize(T0)
 
-    requests  = DF5._dfield.memory_request(op=DF5.name,
-                    request_identifier=DF5.memory_request_id)
+    requests = DF5._dfield.memory_request(op=DF5.name,
+                                          request_identifier=DF5.memory_request_id)
     requests += DF7._dfield.memory_request(op=DF7.name,
-                    request_identifier=DF7.memory_request_id)
+                                           request_identifier=DF7.memory_request_id)
     work = requests.allocate(True)
     DF5.honor_memory_request(work)
     DF7.honor_memory_request(work)
@@ -91,7 +92,7 @@ def test_field():
     assert F4.pretty_name == 'F4'
 
     def func(data, coords, component):
-        (x,y,z) = coords
+        (x, y, z) = coords
         data[...] = x*y*z
 
     DF7.initialize(func)
@@ -117,16 +118,17 @@ def test_field():
     DF7.long_description()
     str(DF7)
 
+
 def test_tensor_field():
     domain = Box(dim=3)
-    T0 = TensorField(domain, 'T0', (3,3))
-    T0_bis = Field(domain, 'T0b', shape=(3,3))
+    T0 = TensorField(domain, 'T0', (3, 3))
+    T0_bis = Field(domain, 'T0b', shape=(3, 3))
     T1 = T0.field_like('T1', dtype=np.float16)
-    T2 = T0[1:,1:]
+    T2 = T0[1:, 1:]
     T2.rename(name='T2')
-    T3 = T0[1,1]
+    T3 = T0[1, 1]
     T4 = TensorField.from_fields(name='T4',
-            fields=(T0[0,0], T0[1,1], T1[0,1], T1[1,0]), shape=(2,2))
+                                 fields=(T0[0, 0], T0[1, 1], T1[0, 1], T1[1, 0]), shape=(2, 2))
     T5 = TensorField.from_field_array('T5', T0._fields)
     T6 = T5.field_like(name='T6')
     T7 = T6.field_like(name='T7', is_tmp=True)
@@ -145,28 +147,28 @@ def test_tensor_field():
     assert T4.name == 'T4'
     assert T5.name == 'T5'
     assert T6.name == 'T6'
-    assert np.array_equal(T0.shape, (3,3))
-    assert np.array_equal(T1.shape, (3,3))
+    assert np.array_equal(T0.shape, (3, 3))
+    assert np.array_equal(T1.shape, (3, 3))
     for (idx, t0) in T0.nd_iter():
         t1 = T1[idx]
         assert t0.domain is domain
         assert t1.domain is domain
         assert t1.dtype != t0.dtype
-        assert t1.name.replace('1','0', 1) == t0.name
-        assert t1.pretty_name.replace('1','0', 1) == t0.pretty_name
+        assert t1.name.replace('1', '0', 1) == t0.name
+        assert t1.pretty_name.replace('1', '0', 1) == t0.pretty_name
         assert t0.dim is 3
         assert t1.dim is 3
-    assert np.array_equal(T0._fields[1:,1:], T2._fields)
-    assert T0._fields[1,1] is T3
-    assert T0._fields[0,0] is T4[0,0]
-    assert T0._fields[1,1] is T4[0,1]
-    assert T1._fields[0,1] is T4[1,0]
-    assert T1._fields[1,0] is T4[1,1]
+    assert np.array_equal(T0._fields[1:, 1:], T2._fields)
+    assert T0._fields[1, 1] is T3
+    assert T0._fields[0, 0] is T4[0, 0]
+    assert T0._fields[1, 1] is T4[0, 1]
+    assert T1._fields[0, 1] is T4[1, 0]
+    assert T1._fields[1, 0] is T4[1, 1]
     assert np.array_equal(T0._fields, T5._fields)
     for f in T0:
         assert f in T0
 
-    D0 = CartesianDiscretization(resolution=(5,5,5), default_boundaries=True)
+    D0 = CartesianDiscretization(resolution=(5, 5, 5), default_boundaries=True)
     topo = CartesianTopology(domain=domain, discretization=D0)
 
     DT0 = T0.discretize(topo)
@@ -179,12 +181,12 @@ def test_tensor_field():
     DT7 = T7.discretize(topo)
 
     DT8 = DiscreteTensorField.from_dfields(name='DT8',
-            dfields=(DT0.dfields[0],
-                     DT1.dfields[0],
-                     DT2.dfields[0],
-                     DT3.dfields[0]),
-            shape=(2,2))
-    DT9 = DT0[:2,:2]
+                                           dfields=(DT0.dfields[0],
+                                                    DT1.dfields[0],
+                                                    DT2.dfields[0],
+                                                    DT3.dfields[0]),
+                                           shape=(2, 2))
+    DT9 = DT0[:2, :2]
     DT10 = DT9.clone()
     DT11, requests = DT10.tmp_dfield_like(name='DT11')
     work = requests().allocate(True)
@@ -197,7 +199,7 @@ def test_tensor_field():
         assert df in DT0
 
     def func(data, coords, component):
-        (x,y,z) = coords
+        (x, y, z) = coords
         data[...] = x*y*z
 
     DT9.rename('foo')
@@ -205,7 +207,7 @@ def test_tensor_field():
     DT9.initialize(DT10.data)
     DT9.fill(4)
     DT9.randomize()
-    DT9.copy(DT0[1:,1:])
+    DT9.copy(DT0[1:, 1:])
 
     DT9.has_ghosts()
     DT9.exchange_ghosts()
@@ -225,7 +227,8 @@ def test_tensor_field():
     DT9.integrate()
 
     assert DT0.nb_components == 9
-    check_instance(DT0.discrete_field_views(), tuple, values=CartesianDiscreteScalarFieldView, size=9)
+    check_instance(DT0.discrete_field_views(), tuple,
+                   values=CartesianDiscreteScalarFieldView, size=9)
     check_instance(DT0.dfields, tuple, values=CartesianDiscreteScalarFieldView, size=9)
     check_instance(DT0.discrete_fields(), tuple, values=CartesianDiscreteField, size=9)
     check_instance(DT0.continuous_fields(), tuple, values=Field, size=9)
@@ -260,7 +263,7 @@ def test_tensor_field():
     assert DT0.has_unique_is_at_left_boundary()
     assert DT0.has_unique_is_at_right_boundary()
 
-    dfield = DT0[1,0]
+    dfield = DT0[1, 0]
     assert DT0.backend == dfield.backend
     assert DT0.backend_kind == dfield.backend_kind
     assert DT0.domain == dfield.domain
@@ -324,18 +327,19 @@ def test_tensor_field():
     except AttributeError:
         pass
 
+
 def test_boundaries():
     """This test checks that all boundaries are compatible for velocity and vorticity."""
-    for dim in (1,2,):
-        i=0
+    for dim in (1, 2,):
+        i = 0
         for (lbd, rbd) in domain_boundary_iterator(dim):
             domain = Box(dim=dim, lboundaries=lbd,
-                                  rboundaries=rbd)
+                         rboundaries=rbd)
             V = VelocityField(domain)
             S = ScalarField(name='S0', domain=domain)
-            divV  = V.div()
+            divV = V.div()
             gradV = V.gradient()
-            lapV  = V.laplacian()
+            lapV = V.laplacian()
             print
             print 'DOMAIN BOUNDARIES:'
             print ' *boundaries=[{}]'.format(domain.format_boundaries())
@@ -352,15 +356,14 @@ def test_boundaries():
             print '{} BOUNDARIES:'.format(lapV.pretty_name)
             for lVi in lapV.fields:
                 print ' *{} boundaries=[{}]'.format(lVi.pretty_name, lVi.format_boundaries())
-            if (dim>1):
+            if (dim > 1):
                 rotV = V.curl()
                 print '{} (VORTICITY) BOUNDARIES:'.format(rotV.pretty_name)
                 for Wi in rotV.fields:
                     print ' *{} boundaries=[{}]'.format(Wi.pretty_name, Wi.format_boundaries())
 
 
-
 if __name__ == '__main__':
-    #test_field()
-    #test_tensor_field()
+    test_field()
+    test_tensor_field()
     test_boundaries()
diff --git a/hysop/numerics/fftw_f/fft3d.f90 b/hysop/numerics/fftw_f/fft3d.f90
index 56a71cda650dc1ba3514f740e84d747a26a36e9b..be54dacf7d3a19e7154d39469f0849fd6b479f0d 100755
--- a/hysop/numerics/fftw_f/fft3d.f90
+++ b/hysop/numerics/fftw_f/fft3d.f90
@@ -779,7 +779,7 @@ contains
        do k = 1, fft_resolution(c_Z)
           do i = 1, local_resolution(c_X)
              coeff = one/((kx(i)**2+ky(j)**2+kz(k)**2))
-             dataout1(i,k,j) = coeff*dataout1(i,k,j)
+             dataout1(i,k,j) = -coeff*dataout1(i,k,j)
           end do
        end do
     end do
diff --git a/hysop/operator/tests/test_poisson.py b/hysop/operator/tests/test_poisson.py
index 693ed2040518b5450f7dbf7901eb7a97aab68875..8d9406699d7a6cd2751fb55c8fba1c4bed894d23 100644
--- a/hysop/operator/tests/test_poisson.py
+++ b/hysop/operator/tests/test_poisson.py
@@ -1,4 +1,5 @@
-import random, primefac
+import random
+import primefac
 from hysop.deps import it, sm, random
 from hysop.constants import HYSOP_REAL, BoundaryCondition
 from hysop.testsenv import __ENABLE_LONG_TESTS__, __HAS_OPENCL_BACKEND__
@@ -11,17 +12,18 @@ from hysop.tools.io_utils import IO
 from hysop.tools.numpywrappers import npw
 from hysop.tools.sympy_utils import truncate_expr, round_expr
 from hysop.tools.spectral_utils import make_multivariate_trigonometric_polynomial, \
-                                       make_multivariate_polynomial
+    make_multivariate_polynomial
 from hysop.operator.poisson import Poisson, Implementation
 
 from hysop import Field, Box
 
+
 class TestPoissonOperator(object):
 
     @classmethod
     def setup_class(cls,
-            enable_extra_tests=__ENABLE_LONG_TESTS__,
-            enable_debug_mode=False):
+                    enable_extra_tests=__ENABLE_LONG_TESTS__,
+                    enable_debug_mode=False):
 
         IO.set_default_path('/tmp/hysop_tests/test_poisson')
 
@@ -29,7 +31,7 @@ class TestPoissonOperator(object):
         cls.size_max = 16
 
         cls.enable_extra_tests = enable_extra_tests
-        cls.enable_debug_mode  = enable_debug_mode
+        cls.enable_debug_mode = enable_debug_mode
 
         from hysop.tools.sympy_utils import enable_pretty_printing
         enable_pretty_printing()
@@ -40,10 +42,10 @@ class TestPoissonOperator(object):
 
     @classmethod
     def build_analytic_solutions(cls, polynomial,
-                                      dim, nb_components,
-                                      lboundaries, rboundaries,
-                                      origin, end):
-        from hysop.symbolic.base  import TensorBase
+                                 dim, nb_components,
+                                 lboundaries, rboundaries,
+                                 origin, end):
+        from hysop.symbolic.base import TensorBase
         from hysop.symbolic.frame import SymbolicFrame
         from hysop.symbolic.field import laplacian
 
@@ -55,23 +57,23 @@ class TestPoissonOperator(object):
             for i in xrange(nb_components):
                 if polynomial:
                     psi, y = make_multivariate_polynomial(origin, end,
-                                                        lboundaries, rboundaries,
-                                                        10, 4)
+                                                          lboundaries, rboundaries,
+                                                          10, 4)
                 else:
                     psi, y = make_multivariate_trigonometric_polynomial(origin, end,
-                            lboundaries, rboundaries, 2)
-                psi = psi.xreplace({yi: xi for (yi,xi) in zip(y, frame.coords)})
+                                                                        lboundaries, rboundaries, 2)
+                psi = psi.xreplace({yi: xi for (yi, xi) in zip(y, frame.coords)})
                 psis += (psi,)
             return npw.asarray(psis).view(TensorBase)
 
-        Psis  = gen_psi()
-        Ws    = npw.atleast_1d(laplacian(Psis, frame))
+        Psis = gen_psi()
+        Ws = npw.atleast_1d(laplacian(Psis, frame))
 
-        fWs   = tuple(sm.lambdify(coords, W)   for W   in Ws)
+        fWs = tuple(sm.lambdify(coords, W) for W in Ws)
         fPsis = tuple(sm.lambdify(coords, Psi) for Psi in Psis)
 
-        analytic_expressions = {'Psi':Psis,  'W':Ws}
-        analytic_functions   = {'Psi':fPsis, 'W':fWs}
+        analytic_expressions = {'Psi': Psis,  'W': Ws}
+        analytic_functions = {'Psi': fPsis, 'W': fWs}
         return (analytic_expressions, analytic_functions)
 
     @staticmethod
@@ -88,23 +90,23 @@ class TestPoissonOperator(object):
         data[...] = fn(*coords).astype(data.dtype)
 
     def _test(self, dim, dtype, max_runs=5,
-                   polynomial=False, size_min=None, size_max=None):
+              polynomial=False, size_min=None, size_max=None):
 
         if (dtype == HYSOP_REAL):
-            nb_components = 1 # enable fortran poisson test
+            nb_components = 1  # enable fortran poisson test
         else:
             nb_components = 2
 
         size_min = first_not_None(size_min, self.size_min)
         size_max = first_not_None(size_max, self.size_max)
 
-        valid_factors = {2,3,5,7,11,13}
+        valid_factors = {2, 3, 5, 7, 11, 13}
         factors = {1}
         while (factors-valid_factors):
             factors.clear()
             shape = tuple(npw.random.randint(low=size_min, high=size_max+1, size=dim).tolist())
             for Si in shape:
-                factors.update( set(primefac.primefac(int(Si))) )
+                factors.update(set(primefac.primefac(int(Si))))
 
         domain_boundaries = list(domain_boundary_iterator(dim=dim))
         periodic = domain_boundaries[0]
@@ -115,39 +117,39 @@ class TestPoissonOperator(object):
         for i, (lboundaries, rboundaries) in enumerate(domain_boundaries, 1):
             domain = Box(origin=(npw.random.rand(dim)-0.5),
                          length=(npw.random.rand(dim)+0.5)*2*npw.pi,
-                            lboundaries=lboundaries,
-                            rboundaries=rboundaries)
+                         lboundaries=lboundaries,
+                         rboundaries=rboundaries)
 
             Psi = Field(domain=domain, name='Psi', dtype=dtype,
-                    nb_components=nb_components, register_object=False)
+                        nb_components=nb_components, register_object=False)
             W = Field(domain=domain, name='W', dtype=dtype,
-                    nb_components=nb_components, register_object=False)
+                      nb_components=nb_components, register_object=False)
 
             self._test_one(shape=shape, dim=dim, dtype=dtype,
-                    domain=domain, Psi=Psi, W=W,
-                    polynomial=polynomial, nb_components=nb_components)
-            if (max_runs is not None) and (i==max_runs):
+                           domain=domain, Psi=Psi, W=W,
+                           polynomial=polynomial, nb_components=nb_components)
+            if (max_runs is not None) and (i == max_runs):
                 missing = ((4**(dim+1) - 1) / 3) - i
                 print
-                print '>> MAX RUNS ACHIEVED FOR {}D DOMAINS -- SKIPING {} OTHER BOUNDARY CONDITIONS <<'.format(dim, missing)
+                print '>> MAX RUNS ACHIEVED FOR {}D DOMAINS -- SKIPING {} OTHER BOUNDARY CONDITIONS <<'.format(
+                    dim, missing)
                 print
                 print
                 break
         else:
-            assert (i==(4**(dim+1)-1)/3), (i+1, (4**(dim+1)-1)/3)
+            assert (i == (4**(dim+1)-1)/3), (i+1, (4**(dim+1)-1)/3)
             print
             print '>> TESTED ALL {}D BOUNDARY CONDITIONS <<'.format(dim)
             print
             print
 
-
     def _test_one(self, shape, dim, dtype,
-            domain, Psi, W, polynomial, nb_components):
+                  domain, Psi, W, polynomial, nb_components):
 
         (analytic_expressions, analytic_functions) = \
             self.build_analytic_solutions(
                 dim=dim, nb_components=nb_components, polynomial=polynomial,
-                lboundaries=W.lboundaries[::-1], # => boundaries in variable order x0,...,xn
+                lboundaries=W.lboundaries[::-1],  # => boundaries in variable order x0,...,xn
                 rboundaries=W.rboundaries[::-1],
                 origin=domain.origin[::-1],
                 end=domain.end[::-1])
@@ -155,82 +157,84 @@ class TestPoissonOperator(object):
         def format_expr(e):
             return truncate_expr(round_expr(e, 3), 80)
 
-        msg='\nTesting {}D Poisson: dtype={} nb_components={} shape={} polynomial={}, bc=[{}]'.format(
-                dim, dtype.__name__, nb_components, shape, polynomial, W.domain.format_boundaries())
-        msg+='\n >Corresponding field boundary conditions are [{}].'.format(W.fields[0].format_boundaries())
-        msg+='\n >Input analytic functions are (truncated):'
+        msg = '\nTesting {}D Poisson: dtype={} nb_components={} shape={} polynomial={}, bc=[{}]'.format(
+            dim, dtype.__name__, nb_components, shape, polynomial, W.domain.format_boundaries())
+        msg += '\n >Corresponding field boundary conditions are [{}].'.format(
+            W.fields[0].format_boundaries())
+        msg += '\n >Input analytic functions are (truncated):'
         for (Wi, Wis) in zip(W.fields, analytic_expressions['W']):
-            msg+='\n   *{}(x,t) = {}'.format(Wi.pretty_name, format_expr(Wis))
-        msg+='\n >Expected output solutions:'
+            msg += '\n   *{}(x,t) = {}'.format(Wi.pretty_name, format_expr(Wis))
+        msg += '\n >Expected output solutions:'
         for (Psi_i, Psis_i) in zip(Psi.fields, analytic_expressions['Psi']):
-            msg+='\n   *{}(x,t) = {}'.format(Psi_i.pretty_name, format_expr(Psis_i))
-        msg+='\n >Testing all implementations:'
+            msg += '\n   *{}(x,t) = {}'.format(Psi_i.pretty_name, format_expr(Psis_i))
+        msg += '\n >Testing all implementations:'
         print msg
 
         implementations = Poisson.implementations()
-        variables = { Psi:shape, W:shape }
+        variables = {Psi: shape, W: shape}
+
         def iter_impl(impl):
             base_kwds = dict(Fout=Psi, Fin=W, variables=variables,
                              implementation=impl,
                              name='poisson_{}'.format(str(impl).lower()))
             if impl is Implementation.FORTRAN:
-                msg='   *Fortran FFTW: '
+                msg = '   *Fortran FFTW: '
                 print msg,
                 yield Poisson(**base_kwds)
             elif impl is Implementation.PYTHON:
-                msg='   *Python FFTW: '
+                msg = '   *Python FFTW: '
                 print msg,
                 yield Poisson(**base_kwds)
             elif impl is Implementation.OPENCL:
                 from hysop.backend.device.opencl import cl
-                msg='   *OpenCl CLFFT: '
+                msg = '   *OpenCl CLFFT: '
                 print msg
                 for cl_env in iter_clenv():
-                    msg='     |platform {}, device {}'.format(cl_env.platform.name.strip(),
-                                                          cl_env.device.name.strip())
+                    msg = '     |platform {}, device {}'.format(cl_env.platform.name.strip(),
+                                                                cl_env.device.name.strip())
                     print msg,
                     yield Poisson(cl_env=cl_env, **base_kwds)
-                msg='   *OpenCl FFTW: '
+                msg = '   *OpenCl FFTW: '
                 print msg
                 cpu_envs = tuple(iter_clenv(device_type='cpu'))
                 if cpu_envs:
                     for cl_env in cpu_envs:
-                        msg='     |platform {}, device {}'.format(cl_env.platform.name.strip(),
-                                                                  cl_env.device.name.strip())
+                        msg = '     |platform {}, device {}'.format(cl_env.platform.name.strip(),
+                                                                    cl_env.device.name.strip())
                         print msg,
                         yield Poisson(cl_env=cl_env, enforce_implementation=False, **base_kwds)
             else:
-                msg='Unknown implementation to test {}.'.format(impl)
+                msg = 'Unknown implementation to test {}.'.format(impl)
                 raise NotImplementedError(msg)
 
-        #Compare to analytic solution
+        # Compare to analytic solution
         Psiref = None
         Wref = None
         for impl in implementations:
             if (impl is Implementation.FORTRAN):
-                if ((nb_components>1) or (dim!=3) or (not dtype is HYSOP_REAL)
+                if ((nb_components > 1) or (dim != 3) or (not dtype is HYSOP_REAL)
                         or any((bd != BoundaryCondition.PERIODIC) for bd in W.lboundaries)
                         or any((bd != BoundaryCondition.PERIODIC) for bd in W.rboundaries)):
                     print '   *Fortran FFTW: NO SUPPORT'
                     continue
             for op in iter_impl(impl):
-                op   = op.build()
-                dw   = op.get_input_discrete_field(W).as_contiguous_dfield()
+                op = op.build()
+                dw = op.get_input_discrete_field(W).as_contiguous_dfield()
                 dpsi = op.get_output_discrete_field(Psi).as_contiguous_dfield()
 
                 dw.initialize(self.__analytic_init, dtype=dtype,
                               fns=analytic_functions['W'])
                 if (Psiref is None):
                     dpsi.initialize(self.__analytic_init, dtype=dtype,
-                            fns=analytic_functions['Psi'])
-                    Wref   = tuple( data.get().handle.copy() for data in dw.data   )
-                    Psiref = tuple( data.get().handle.copy() for data in dpsi.data )
+                                    fns=analytic_functions['Psi'])
+                    Wref = tuple(data.get().handle.copy() for data in dw.data)
+                    Psiref = tuple(data.get().handle.copy() for data in dpsi.data)
                 dpsi.initialize(self.__random_init, dtype=dtype)
 
                 op.apply(simulation=None)
 
-                Wout   = tuple( data.get().handle.copy() for data in dw.data   )
-                Psiout = tuple( data.get().handle.copy() for data in dpsi.data )
+                Wout = tuple(data.get().handle.copy() for data in dw.data)
+                Psiout = tuple(data.get().handle.copy() for data in dpsi.data)
                 self._check_output(impl, op, Wref, Psiref, Wout, Psiout)
                 if (impl is Implementation.FORTRAN):
                     op.finalize(clean_fftw_solver=True)
@@ -244,9 +248,9 @@ class TestPoissonOperator(object):
         check_instance(Psiout, tuple, values=npw.ndarray, size=len(Wref))
 
         msg0 = 'Reference field {} is not finite.'
-        for (fields, name) in zip((Wref, Psiref),('Wref', 'Psiref')):
-            for (i,field) in enumerate(fields):
-                iname = '{}{}'.format(name,i)
+        for (fields, name) in zip((Wref, Psiref), ('Wref', 'Psiref')):
+            for (i, field) in enumerate(fields):
+                iname = '{}{}'.format(name, i)
                 mask = npw.isfinite(field)
                 if not mask.all():
                     print
@@ -258,9 +262,9 @@ class TestPoissonOperator(object):
                     raise ValueError(msg)
 
         for (out_buffers, ref_buffers, name) in zip((Wout, Psiout),
-                                                        (Wref, Psiref), ('W', 'Psi')):
-            for i, (fout,fref) in enumerate(zip(out_buffers, ref_buffers)):
-                iname = '{}{}'.format(name,i)
+                                                    (Wref, Psiref), ('W', 'Psi')):
+            for i, (fout, fref) in enumerate(zip(out_buffers, ref_buffers)):
+                iname = '{}{}'.format(name, i)
                 assert fout.dtype == fref.dtype, iname
                 assert fout.shape == fref.shape, iname
                 assert fout.flags.c_contiguous
@@ -273,7 +277,7 @@ class TestPoissonOperator(object):
                 elif has_inf:
                     deps = 'inf'
                 else:
-                    eps  = npw.finfo(fout.dtype).eps
+                    eps = npw.finfo(fout.dtype).eps
                     dist = npw.abs(fout-fref)
                     dinf = npw.max(dist)
                     deps = int(npw.ceil(dinf/eps))
@@ -290,53 +294,66 @@ class TestPoissonOperator(object):
                 print
                 if cls.enable_debug_mode:
                     print 'REFERENCE INPUTS:'
-                    for (i,w) in enumerate(Wref):
+                    for (i, w) in enumerate(Wref):
                         print 'W{}'.format(i)
                         print w
                         print
                     if (name == 'Psi'):
                         print 'REFERENCE OUTPUT:'
-                        for (i,u) in enumerate(Psiref):
+                        for (i, u) in enumerate(Psiref):
                             print 'Psi{}'.format(i)
                             print u
                             print
                         print
                         print 'OPERATOR {} OUTPUT:'.format(op.name.upper())
                         print
-                        for (i,u) in enumerate(Psiout):
+                        for (i, u) in enumerate(Psiout):
                             print 'Psi{}'.format(i)
                             print u
                             print
                     else:
                         print 'MODIFIED INPUTS:'
-                        for (i,w) in enumerate(Wout):
+                        for (i, w) in enumerate(Wout):
                             print 'W{}'.format(i)
                             print w
                             print
                     print
 
                 msg = 'Test failed for {} on component {} for implementation {}.'
-                msg=msg.format(name, i, impl)
+                msg = msg.format(name, i, impl)
                 raise RuntimeError(msg)
 
-
     def test_1d_float32(self, **kwds):
-        self._test(dim=1, dtype=npw.float32, **kwds)
+        if HYSOP_REAL == npw.float32:
+            self._test(dim=1, dtype=npw.float32, **kwds)
+
     def test_2d_float32(self, **kwds):
-        self._test(dim=2, dtype=npw.float32, **kwds)
+        if HYSOP_REAL == npw.float32:
+            self._test(dim=2, dtype=npw.float32, **kwds)
+
     def test_3d_float32(self, **kwds):
-        self._test(dim=3, dtype=npw.float32, **kwds)
+        if HYSOP_REAL == npw.float32:
+            self._test(dim=3, dtype=npw.float32, **kwds)
+
     def test_4d_float32(self, **kwds):
-        self._test(dim=4, dtype=npw.float32, **kwds)
+        if HYSOP_REAL == npw.float32:
+            self._test(dim=4, dtype=npw.float32, **kwds)
 
     def test_1d_float64(self, **kwds):
-        self._test(dim=1, dtype=npw.float64, **kwds)
+        if HYSOP_REAL == npw.float64:
+            self._test(dim=1, dtype=npw.float64, **kwds)
+
     def test_2d_float64(self, **kwds):
-        self._test(dim=2, dtype=npw.float64, **kwds)
+        if HYSOP_REAL == npw.float64:
+            self._test(dim=2, dtype=npw.float64, **kwds)
+
     def test_3d_float64(self, **kwds):
-        self._test(dim=3, dtype=npw.float64, **kwds)
+        if HYSOP_REAL == npw.float64:
+            self._test(dim=3, dtype=npw.float64, **kwds)
+
     def test_4d_float64(self, **kwds):
-        self._test(dim=4, dtype=npw.float64, **kwds)
+        if HYSOP_REAL == npw.float64:
+            self._test(dim=4, dtype=npw.float64, **kwds)
 
     # def test_polynomial_1d_float32(self, **kwds):
     #     self._test(dim=1, dtype=npw.float32, polynomial=True, **kwds)
@@ -345,26 +362,26 @@ class TestPoissonOperator(object):
     # def test_polynomial_3d_float32(self, **kwds):
     #     self._test(dim=3, dtype=npw.float32, polynomial=True, **kwds)
 
-
     def perform_tests(self):
         max_1d_runs = None
         max_2d_runs = 2
         max_3d_runs = 2
         max_4d_runs = 2
 
-        if __ENABLE_LONG_TESTS__ or (HYSOP_REAL==npw.float32):
+        if __ENABLE_LONG_TESTS__ or (HYSOP_REAL == npw.float32):
             self.test_1d_float32(max_runs=max_1d_runs)
             self.test_2d_float32(max_runs=max_2d_runs)
             if __ENABLE_LONG_TESTS__:
                 self.test_3d_float32(max_runs=max_3d_runs)
                 self.test_4d_float32(max_runs=max_4d_runs)
-        if __ENABLE_LONG_TESTS__ or (HYSOP_REAL==npw.float64):
+        if __ENABLE_LONG_TESTS__ or (HYSOP_REAL == npw.float64):
             self.test_1d_float64(max_runs=max_1d_runs)
             self.test_2d_float64(max_runs=max_2d_runs)
             if __ENABLE_LONG_TESTS__:
                 self.test_3d_float64(max_runs=max_3d_runs)
                 self.test_4d_float32(max_runs=max_4d_runs)
 
+
 if __name__ == '__main__':
     TestPoissonOperator.setup_class(enable_extra_tests=False,
                                     enable_debug_mode=False)
diff --git a/hysop/operator/tests/test_scales_advection.py b/hysop/operator/tests/test_scales_advection.py
index fa523cc8732f0463735c7d13aa439ad99e516951..10e40fd71a51e9600033dd6cb0f3135999271752 100644
--- a/hysop/operator/tests/test_scales_advection.py
+++ b/hysop/operator/tests/test_scales_advection.py
@@ -18,12 +18,13 @@ from hysop.numerics.splitting.strang import StrangSplitting, StrangOrder
 from hysop.numerics.odesolvers.runge_kutta import Euler, RK2, RK4
 from hysop.numerics.remesh.remesh import RemeshKernel
 
+
 class TestScalesAdvectionOperator(object):
 
     @classmethod
     def setup_class(cls,
-            enable_extra_tests=__ENABLE_LONG_TESTS__,
-            enable_debug_mode=False):
+                    enable_extra_tests=__ENABLE_LONG_TESTS__,
+                    enable_debug_mode=False):
 
         IO.set_default_path('/tmp/hysop_tests/test_scales_advection')
 
@@ -35,15 +36,14 @@ class TestScalesAdvectionOperator(object):
             cls.size_max = 23
 
         cls.enable_extra_tests = enable_extra_tests
-        cls.enable_debug_mode  = enable_debug_mode
+        cls.enable_debug_mode = enable_debug_mode
 
     @classmethod
     def teardown_class(cls):
         pass
 
-
     def _test(self, dim, is_inplace,
-            size_min=None, size_max=None):
+              size_min=None, size_max=None):
         assert dim > 0
 
         # periodic boundaries removes one computational point
@@ -56,7 +56,7 @@ class TestScalesAdvectionOperator(object):
         if self.enable_extra_tests:
             flt_types = (HYSOP_REAL, )
             time_integrators = (RK2, )
-            remesh_kernels =  (Remesh.L2_1, Remesh.L4_2, Remesh.L6_4)
+            remesh_kernels = (Remesh.L2_1, Remesh.L4_2, Remesh.L6_4)
             velocity_cfls = (0.62, 1.89)
         else:
             flt_types = (HYSOP_REAL,)
@@ -66,20 +66,20 @@ class TestScalesAdvectionOperator(object):
 
         domain = Box(length=(2*npw.pi,)*dim)
         for dtype in flt_types:
-            Vin  = Field(domain=domain, name='Vin', dtype=dtype,
-                    nb_components=dim, register_object=False)
-            Sin  = Field(domain=domain, name='Sin', dtype=dtype,
-                    nb_components=5, register_object=False)
+            Vin = Field(domain=domain, name='Vin', dtype=dtype,
+                        nb_components=dim, register_object=False)
+            Sin = Field(domain=domain, name='Sin', dtype=dtype,
+                        nb_components=5, register_object=False)
             Sout = Field(domain=domain, name='Sout', dtype=dtype,
-                    nb_components=5, register_object=False)
+                         nb_components=5, register_object=False)
             for time_integrator in time_integrators:
                 for remesh_kernel in remesh_kernels:
                     for velocity_cfl in velocity_cfls:
                         print
                         self._test_one(time_integrator=time_integrator, remesh_kernel=remesh_kernel,
-                                        shape=shape, dim=dim, dtype=dtype,
-                                        is_inplace=is_inplace, domain=domain,
-                                        Vin=Vin, Sin=Sin, Sout=Sout, velocity_cfl=velocity_cfl)
+                                       shape=shape, dim=dim, dtype=dtype,
+                                       is_inplace=is_inplace, domain=domain,
+                                       Vin=Vin, Sin=Sin, Sout=Sout, velocity_cfl=velocity_cfl)
 
     @classmethod
     def __velocity_init(cls, data, coords, component, axes):
@@ -91,37 +91,34 @@ class TestScalesAdvectionOperator(object):
     @classmethod
     def __scalar_init(cls, data, coords, component, offsets=None):
         offsets = first_not_None(offsets, (0.0,)*len(coords))
-        assert len(coords)==len(offsets)
+        assert len(coords) == len(offsets)
         data[...] = 1.0/(component+1)
         for (c, o) in zip(coords, offsets):
             data[...] *= npw.cos(c+o)
 
     def _test_one(self, time_integrator, remesh_kernel,
-            shape, dim,
-            dtype, is_inplace, domain, velocity_cfl,
-            Vin, Sin, Sout):
+                  shape, dim,
+                  dtype, is_inplace, domain, velocity_cfl,
+                  Vin, Sin, Sout):
 
         print '\nTesting {}D ScalesAdvection_{}_{}: inplace={} dtype={} shape={}, cfl={}'.format(
-                dim, time_integrator.name(), remesh_kernel,
-                is_inplace, dtype.__name__, shape, velocity_cfl),
+            dim, time_integrator.name(), remesh_kernel,
+            is_inplace, dtype.__name__, shape, velocity_cfl),
         if is_inplace:
             vin = Vin
             sin, sout = Sin, Sin
-            variables = { vin : shape, sin: shape }
+            variables = {vin: shape, sin: shape}
         else:
             vin = Vin
             sin, sout = Sin, Sout
-            variables = { vin: shape, sin: shape, sout: shape }
+            variables = {vin: shape, sin: shape, sout: shape}
 
         # Use optimal timestep, ||Vi||_inf is 1 on a per-axis basis
         dt = ScalarParameter('dt', initial_value=npw.nan)
         dt.value = (0.99 * velocity_cfl) / (max(shape)-1)
 
-        ref_impl = Implementation.PYTHON
-        implementations = Advection.implementations().keys()
-        msg='Python implementation is currently treated as a directional operator.'
-        assert (ref_impl not in implementations), msg
-        implementations = [ref_impl] + implementations
+        implementations = [Implementation.FORTRAN, ]
+        ref_impl = implementations[0]
 
         method = {TimeIntegrator: time_integrator, Remesh: remesh_kernel}
 
@@ -130,8 +127,8 @@ class TestScalesAdvectionOperator(object):
                              variables=variables, implementation=impl,
                              method=method, name='advection_{}'.format(str(impl).lower()))
             force_tstate = ForceTopologyState(fields=variables.keys(),
-                                            variables=variables,
-                                            backend=Backend.HOST)
+                                              variables=variables,
+                                              backend=Backend.HOST)
             graph = ComputationalGraph(name='test_graph')
             if impl is Implementation.PYTHON:
                 da = DirectionalAdvection(velocity_cfl=velocity_cfl, **base_kwds)
@@ -141,18 +138,18 @@ class TestScalesAdvectionOperator(object):
                 graph.push_nodes(split, force_tstate)
                 yield 'default', graph
             elif impl is Implementation.FORTRAN:
-                assert dim==3, "Scales is only 3D"
+                assert dim == 3, "Scales is only 3D"
                 adv = Advection(**base_kwds)
                 graph.push_nodes(adv, force_tstate)
                 yield 'SCALES', graph
             else:
-                msg='Unknown implementation to test {}.'.format(impl)
+                msg = 'Unknown implementation to test {}.'.format(impl)
                 raise NotImplementedError(msg)
 
         # Compare to other implementations
         advec_axes = (tuple(),)
         advec_axes += tuple((x,) for x in xrange(dim))
-        if (dim>1):
+        if (dim > 1):
             advec_axes += (tuple(xrange(dim)),)
 
         reference_fields = {}
@@ -163,16 +160,15 @@ class TestScalesAdvectionOperator(object):
                 print '\n   *{}: '.format(sop),
 
                 graph.build()
-
                 for axes in advec_axes:
-                    #print 'SWITCHING TO AXES {}'.format(axes)
+                    # print 'SWITCHING TO AXES {}'.format(axes)
                     ref_outputs = reference_fields.setdefault(axes, {})
                     napplies = 10
                     Vi = npw.asarray([+1.0 if (i in axes) else +0.0
-                                        for i in xrange(dim)], dtype=dtype)
+                                      for i in xrange(dim)], dtype=dtype)
 
-                    dvin  = graph.get_input_discrete_field(vin).as_contiguous_dfield()
-                    dsin  = graph.get_input_discrete_field(sin).as_contiguous_dfield()
+                    dvin = graph.get_input_discrete_field(vin).as_contiguous_dfield()
+                    dsin = graph.get_input_discrete_field(sin).as_contiguous_dfield()
                     dsout = graph.get_output_discrete_field(sout).as_contiguous_dfield()
                     dsref = dsout.clone()
 
@@ -181,27 +177,31 @@ class TestScalesAdvectionOperator(object):
                     try:
                         dvin.initialize(self.__velocity_init, axes=axes)
                         dsin.initialize(self.__scalar_init)
+                        dsout.initialize(self.__scalar_init)
+
                         _input = tuple(dsin.data[i].get().handle.copy()
-                                for i in xrange(dsin.nb_components))
+                                       for i in xrange(dsin.nb_components))
                         S0 = dsin.integrate()
 
                         for k in xrange(napplies+1):
-                            if (k>0):
+                            if (k > 0):
                                 graph.apply()
 
                             output = tuple(dsout.data[i].get().handle.copy()[dsout.compute_slices]
-                                        for i in xrange(dsout.nb_components))
+                                           for i in xrange(dsout.nb_components))
 
                             for i in xrange(dsout.nb_components):
                                 mask = npw.isfinite(output[i][dsout.compute_slices])
                                 if not mask.all():
-                                    msg='\nFATAL ERROR: Output is not finite on axis {}.\n'.format(i)
+                                    msg = '\nFATAL ERROR: Output is not finite on axis {}.\n'.format(
+                                        i)
                                     print msg
-                                    npw.fancy_print(output[i], replace_values={(lambda a: npw.isfinite(a)): '.'})
+                                    npw.fancy_print(output[i], replace_values={
+                                                    (lambda a: npw.isfinite(a)): '.'})
                                     raise RuntimeError(msg)
 
                             if is_ref:
-                                dxk = -Vi*(k+0)*dt()
+                                dxk = -Vi*(k)*dt()
                                 dsref.initialize(self.__scalar_init, offsets=dxk.tolist())
                                 d = dsout.distance(dsref, p=2)
                                 if npw.any(d > 1e-1):
@@ -211,9 +211,11 @@ class TestScalesAdvectionOperator(object):
                                     print 'DSREF'
                                     print dsref.sdata[dsref.compute_slices]
                                     print 'DSREF - DSOUT'
-                                    print (dsout.sdata[dsout.compute_slices].get() - dsref.sdata[dsref.compute_slices].get())
-                                    msg='Test failed with V={}, k={}, dxk={}, inter-field L2 distances are {}.'
-                                    msg=msg.format(Vi, k, to_tuple(dxk, cast=float), to_tuple(d, cast=float))
+                                    print (dsout.sdata[dsout.compute_slices].get(
+                                    ) - dsref.sdata[dsref.compute_slices].get())
+                                    msg = 'Test failed with V={}, k={}, dxk={}, inter-field L2 distances are {}.'
+                                    msg = msg.format(Vi, k, to_tuple(
+                                        dxk, cast=float), to_tuple(d, cast=float))
                                     raise RuntimeError(msg)
                                 ref_outputs[k] = output
                             else:
@@ -224,10 +226,11 @@ class TestScalesAdvectionOperator(object):
                                     max_di = npw.max(di)
                                     neps = 10000
                                     max_tol = neps*npw.finfo(dsout.dtype).eps
-                                    if (max_di>max_tol):
+                                    if (max_di > max_tol):
                                         print 'FATAL ERROR: Could not match other implementation results.'
-                                        print '\nComparisson failed at step {} and component {}:'.format(k,i)
-                                        for (j,dv) in dvin.iter_fields():
+                                        print '\nComparisson failed at step {} and component {}:'.format(
+                                            k, i)
+                                        for (j, dv) in dvin.iter_fields():
                                             print 'VELOCITY INPUT {}'.format(DirectionLabels[j])
                                             print dv.sdata[dv.compute_slices]
                                         print 'SCALAR INPUT'
@@ -237,22 +240,24 @@ class TestScalesAdvectionOperator(object):
                                         print 'SCALAR OUTPUT'
                                         print output[i]
                                         print 'ABS(REF - OUT)'
-                                        npw.fancy_print(di, replace_values={(lambda a: a<max_tol): '.'})
+                                        npw.fancy_print(di, replace_values={
+                                                        (lambda a: a < max_tol): '.'})
                                         print
-                                        msg='Output did not match reference output for component {} at time step {}.'
-                                        msg+='\n > max computed distance was {}.'.format(max_di)
-                                        msg+='\n > max tolerence was set to {} ({} eps).'.format(max_tol, neps)
-                                        msg=msg.format(i, k)
+                                        msg = 'Output did not match reference output for component {} at time step {}.'
+                                        msg += '\n > max computed distance was {}.'.format(max_di)
+                                        msg += '\n > max tolerence was set to {} ({} eps).'.format(
+                                            max_tol, neps)
+                                        msg = msg.format(i, k)
                                         raise RuntimeError(msg)
                             Si = dsout.integrate()
                             if not npw.all(npw.isfinite(Si)):
-                                msg='Integral is not finite. Got {}.'.format(Si)
+                                msg = 'Integral is not finite. Got {}.'.format(Si)
                                 raise RuntimeError(msg)
-                            if (npw.abs(Si-S0)>1e-3).any():
-                                msg='Scalar was not conserved on iteration {}, expected {} but got {}.'
-                                msg=msg.format(k,
-                                        to_tuple(S0, cast=float),
-                                        to_tuple(Si, cast=float))
+                            if (npw.abs(Si-S0) > 1e-3).any():
+                                msg = 'Scalar was not conserved on iteration {}, expected {} but got {}.'
+                                msg = msg.format(k,
+                                                 to_tuple(S0, cast=float),
+                                                 to_tuple(Si, cast=float))
                                 raise RuntimeError(msg)
                     except:
                         sys.stdout.write('\bx\n\n')
@@ -271,7 +276,7 @@ class TestScalesAdvectionOperator(object):
 if __name__ == '__main__':
     import hysop
     TestScalesAdvectionOperator.setup_class(enable_extra_tests=False,
-                                      enable_debug_mode=False)
+                                            enable_debug_mode=False)
 
     test = TestScalesAdvectionOperator()
 
diff --git a/hysop/operator/tests/test_spectral_derivative.py b/hysop/operator/tests/test_spectral_derivative.py
index f69e7df2af6312056e90466fb5c7b5504acc93ae..373bf7815ff11a59caad9002580f6db99a2bb118 100644
--- a/hysop/operator/tests/test_spectral_derivative.py
+++ b/hysop/operator/tests/test_spectral_derivative.py
@@ -13,7 +13,7 @@ from hysop.tools.io_utils import IO
 from hysop.tools.numpywrappers import npw
 from hysop.tools.sympy_utils import truncate_expr, round_expr
 from hysop.tools.spectral_utils import make_multivariate_trigonometric_polynomial, \
-                                       make_multivariate_polynomial
+    make_multivariate_polynomial
 from hysop.parameters.scalar_parameter import ScalarParameter
 from hysop.operator.derivative import Implementation, SpectralSpaceDerivative
 from hysop.operator.gradient import Gradient
@@ -21,12 +21,13 @@ from hysop.operator.misc import ForceTopologyState
 
 from hysop import Field, Box
 
+
 class TestSpectralDerivative(object):
 
     @classmethod
     def setup_class(cls,
-            enable_extra_tests=__ENABLE_LONG_TESTS__,
-            enable_debug_mode=False):
+                    enable_extra_tests=__ENABLE_LONG_TESTS__,
+                    enable_debug_mode=False):
 
         IO.set_default_path('/tmp/hysop_tests/test_spectral_derivative')
 
@@ -34,13 +35,13 @@ class TestSpectralDerivative(object):
         cls.size_max = 16
 
         cls.enable_extra_tests = enable_extra_tests
-        cls.enable_debug_mode  = enable_debug_mode
+        cls.enable_debug_mode = enable_debug_mode
 
         cls.t = ScalarParameter(name='t', dtype=HYSOP_REAL)
 
     @classmethod
     def build_analytic_expressions(cls, polynomial, dim, max_derivative,
-                                        lboundaries, rboundaries, origin, end):
+                                   lboundaries, rboundaries, origin, end):
         from hysop.tools.sympy_utils import enable_pretty_printing
         from hysop.symbolic.base import TensorBase
         from hysop.symbolic.frame import SymbolicFrame
@@ -58,35 +59,34 @@ class TestSpectralDerivative(object):
                                                     10, 4)
             else:
                 f, y = make_multivariate_trigonometric_polynomial(origin, end,
-                        lboundaries, rboundaries, 2)
-            f = f.xreplace({yi: xi for (yi,xi) in zip(y, frame.coords)})
-            f *= sm.Integer(1) / (sm.Integer(1) + npw.random.randint(1,5)*cls.t.s)
+                                                                  lboundaries, rboundaries, 2)
+            f = f.xreplace({yi: xi for (yi, xi) in zip(y, frame.coords)})
+            f *= sm.Integer(1) / (sm.Integer(1) + npw.random.randint(1, 5)*cls.t.s)
             return f
 
-        F  = gen_F()
+        F = gen_F()
         fF = sm.lambdify(params, F)
 
-        dFs  = {}
+        dFs = {}
         fdFs = {}
         symbolic_dvars = {}
         for idx in it.product(range(max_derivative+1), repeat=dim):
-            if sum(idx)> max_derivative:
+            if sum(idx) > max_derivative:
                 continue
-            xvars = tuple((ci,i) for (i,ci) in zip(idx, coords))
+            xvars = tuple((ci, i) for (i, ci) in zip(idx, coords))
             symbolic_dvars[idx] = xvars
             dF = F
-            for (ci,i) in xvars:
-                if (i==0):
+            for (ci, i) in xvars:
+                if (i == 0):
                     continue
-                dF = dF.diff(ci,i)
-            dFs[idx]  = dF
+                dF = dF.diff(ci, i)
+            dFs[idx] = dF
             fdFs[idx] = sm.lambdify(params, dF)
 
-        analytic_expressions = {'F':F,  'dF':dFs}
-        analytic_functions = {'F':fF, 'dF':fdFs}
+        analytic_expressions = {'F': F,  'dF': dFs}
+        analytic_functions = {'F': fF, 'dF': fdFs}
         return (symbolic_dvars, analytic_expressions, analytic_functions)
 
-
     @classmethod
     def teardown_class(cls):
         pass
@@ -106,7 +106,7 @@ class TestSpectralDerivative(object):
         data[...] = npw.asarray(fn(*(coords+(t(),)))).astype(data.dtype)
 
     def _test(self, dim, dtype, polynomial, max_derivative=2,
-            size_min=None, size_max=None, max_runs=None):
+              size_min=None, size_max=None, max_runs=None):
         enable_extra_tests = self.enable_extra_tests
 
         size_min = first_not_None(size_min, self.size_min)
@@ -130,27 +130,27 @@ class TestSpectralDerivative(object):
             F = Field(domain=domain, name='F', dtype=dtype)
 
             self._test_one(shape=shape, dim=dim, dtype=dtype,
-                    domain=domain, F=F,
-                    polynomial=polynomial,
-                    max_derivative=max_derivative)
+                           domain=domain, F=F,
+                           polynomial=polynomial,
+                           max_derivative=max_derivative)
 
-            if (max_runs is not None) and (i==max_runs):
+            if (max_runs is not None) and (i == max_runs):
                 missing = ((4**(dim+1) - 1) / 3) - i
                 print
-                print '>> MAX RUNS ACHIEVED FOR {}D DOMAINS -- SKIPING {} OTHER BOUNDARY CONDITIONS <<'.format(dim, missing)
+                print '>> MAX RUNS ACHIEVED FOR {}D DOMAINS -- SKIPING {} OTHER BOUNDARY CONDITIONS <<'.format(
+                    dim, missing)
                 print
                 print
                 break
         else:
-            assert (i==(4**(dim+1)-1)/3), (i+1, (4**(dim+1)-1)/3)
+            assert (i == (4**(dim+1)-1)/3), (i+1, (4**(dim+1)-1)/3)
             print
             print '>> TESTED ALL {}D BOUNDARY CONDITIONS <<'.format(dim)
             print
             print
 
-
     def _test_one(self, shape, dim, dtype,
-            domain, F, polynomial, max_derivative):
+                  domain, F, polynomial, max_derivative):
 
         implementations = SpectralSpaceDerivative.implementations()
 
@@ -158,39 +158,39 @@ class TestSpectralDerivative(object):
             self.build_analytic_expressions(
                 dim=dim, polynomial=polynomial,
                 max_derivative=max_derivative,
-                lboundaries=F.lboundaries[::-1], # => boundaries in variable order x0,...,xn
+                lboundaries=F.lboundaries[::-1],  # => boundaries in variable order x0,...,xn
                 rboundaries=F.rboundaries[::-1],
                 origin=domain.origin[::-1],
                 end=domain.end[::-1])
 
-        Fs  = analytic_expressions['F']
+        Fs = analytic_expressions['F']
         fFs = analytic_functions['F']
 
         def format_expr(e):
             return truncate_expr(round_expr(e, 3), 80)
 
-        msg='\nTesting {}D SpectralDerivative: dtype={} shape={}, polynomial={}, bc=[{}]'
-        msg=msg.format(dim, dtype.__name__, shape, polynomial, F.domain.format_boundaries())
-        msg+='\n >Corresponding field boundary conditions are [{}].'.format(F.format_boundaries())
-        msg+='\n >Input analytic functions (truncated):'
-        msg+='\n   *{}(x,t) = {}'.format(F.pretty_name, format_expr(Fs))
-        msg+='\n >Testing derivatives:'
+        msg = '\nTesting {}D SpectralDerivative: dtype={} shape={}, polynomial={}, bc=[{}]'
+        msg = msg.format(dim, dtype.__name__, shape, polynomial, F.domain.format_boundaries())
+        msg += '\n >Corresponding field boundary conditions are [{}].'.format(F.format_boundaries())
+        msg += '\n >Input analytic functions (truncated):'
+        msg += '\n   *{}(x,t) = {}'.format(F.pretty_name, format_expr(Fs))
+        msg += '\n >Testing derivatives:'
         print msg
 
         for idx in sorted(symbolic_dvars.keys(), key=lambda x: sum(x)):
             xvars = symbolic_dvars[idx]
             dFe = F.s()
-            for (ci,i) in xvars:
-                if (i==0):
+            for (ci, i) in xvars:
+                if (i == 0):
                     continue
-                dFe = dFe.diff(ci,i)
+                dFe = dFe.diff(ci, i)
             dF = F.from_sympy_expression(expr=dFe,
-                    space_symbols=domain.frame.coords)
-            dFs   = analytic_expressions['dF'][idx]
-            fdFs  = analytic_functions['dF'][idx]
+                                         space_symbols=domain.frame.coords)
+            dFs = analytic_expressions['dF'][idx]
+            fdFs = analytic_functions['dF'][idx]
             print '   *{}'.format(dF.pretty_name)
 
-            variables = { F:shape, dF: shape }
+            variables = {F: shape, dF: shape}
 
             def iter_impl(impl):
                 base_kwds = dict(F=F, dF=dF, derivative=idx,
@@ -198,24 +198,24 @@ class TestSpectralDerivative(object):
                                  implementation=impl,
                                  testing=True)
                 if impl is Implementation.PYTHON:
-                    msg='     |Python: '
+                    msg = '     |Python: '
                     print msg,
                     op = SpectralSpaceDerivative(**base_kwds)
                     yield op.to_graph()
                     print
                 elif impl is Implementation.OPENCL:
-                    msg='     |Opencl: '
+                    msg = '     |Opencl: '
                     print msg
                     for cl_env in iter_clenv():
-                        msg='        >platform {}, device {}:'.format(
-                                                              cl_env.platform.name.strip(),
-                                                              cl_env.device.name.strip())
+                        msg = '        >platform {}, device {}:'.format(
+                            cl_env.platform.name.strip(),
+                            cl_env.device.name.strip())
                         print msg,
                         op = SpectralSpaceDerivative(cl_env=cl_env, **base_kwds)
                         yield op.to_graph()
                         print
                 else:
-                    msg='Unknown implementation to test {}.'.format(impl)
+                    msg = 'Unknown implementation to test {}.'.format(impl)
                     raise NotImplementedError(msg)
 
             # Compare to analytic solution
@@ -223,23 +223,23 @@ class TestSpectralDerivative(object):
             for impl in implementations:
                 for op in iter_impl(impl):
                     op.build(outputs_are_inputs=False)
-                    #op.display()
+                    # op.display()
 
-                    Fd  = op.get_input_discrete_field(F)
+                    Fd = op.get_input_discrete_field(F)
                     dFd = op.get_output_discrete_field(dF)
 
                     if (Fref is None):
                         dFd.initialize(self.__analytic_init, fns=(fdFs,), t=self.t)
-                        dFref = tuple( data.get().handle.copy() for data in dFd.data )
+                        dFref = tuple(data.get().handle.copy() for data in dFd.data)
 
                     Fd.initialize(self.__analytic_init, fns=(fFs,), t=self.t)
-                    Fref = tuple( data.get().handle.copy() for data in Fd.data )
+                    Fref = tuple(data.get().handle.copy() for data in Fd.data)
 
                     dFd.initialize(self.__random_init)
                     op.apply()
 
-                    Fout  = tuple( data.get().handle.copy() for data in Fd.data )
-                    dFout = tuple( data.get().handle.copy() for data in dFd.data )
+                    Fout = tuple(data.get().handle.copy() for data in Fd.data)
+                    dFout = tuple(data.get().handle.copy() for data in dFd.data)
 
                     self._check_output(impl, op, Fref, dFref, Fout, dFout, idx)
 
@@ -251,11 +251,11 @@ class TestSpectralDerivative(object):
         check_instance(Fout,  tuple, values=npw.ndarray, size=len(Fref))
         check_instance(dFout, tuple, values=npw.ndarray, size=len(dFref))
 
-        for j,(out_buffers, ref_buffers, name) in enumerate(zip((Fout, dFout),
-                                                                (Fref, dFref),
-                                                                ('F', 'dF'))):
-            for i, (fout,fref) in enumerate(zip(out_buffers, ref_buffers)):
-                iname = '{}{}'.format(name,i)
+        for j, (out_buffers, ref_buffers, name) in enumerate(zip((Fout, dFout),
+                                                                 (Fref, dFref),
+                                                                 ('F', 'dF'))):
+            for i, (fout, fref) in enumerate(zip(out_buffers, ref_buffers)):
+                iname = '{}{}'.format(name, i)
                 assert fout.dtype == fref.dtype, iname
                 assert fout.shape == fref.shape, iname
 
@@ -267,7 +267,7 @@ class TestSpectralDerivative(object):
                 if (has_nan or has_inf):
                     pass
                 else:
-                    eps  = npw.finfo(fout.dtype).eps
+                    eps = npw.finfo(fout.dtype).eps
                     dist = npw.abs(fout-fref)
                     dinf = npw.max(dist)
                     try:
@@ -276,7 +276,7 @@ class TestSpectralDerivative(object):
                         import numpy as np
                         deps = np.inf
                     if (deps <= 10**(nidx+2)):
-                        if (j==1):
+                        if (j == 1):
                             print '{}eps ({})'.format(deps, dinf),
                         else:
                             print '{}eps, '.format(deps),
@@ -292,57 +292,64 @@ class TestSpectralDerivative(object):
                 print
                 if cls.enable_debug_mode:
                     print 'REFERENCE INPUTS:'
-                    for (i,w) in enumerate(Fref):
+                    for (i, w) in enumerate(Fref):
                         print 'F{}'.format(i)
                         print w
                         print
                     if (name == 'dF'):
                         print 'REFERENCE OUTPUT:'
-                        for (i,u) in enumerate(dFref):
+                        for (i, u) in enumerate(dFref):
                             print 'dF{}'.format(i)
                             print u
                             print
                         print
                         print 'OPERATOR {} OUTPUT:'.format(op.name.upper())
                         print
-                        for (i,u) in enumerate(dFout):
+                        for (i, u) in enumerate(dFout):
                             print 'dF{}'.format(i)
                             print u
                             print
                     else:
                         print 'MODIFIED INPUTS:'
-                        for (i,w) in enumerate(Fout):
+                        for (i, w) in enumerate(Fout):
                             print 'F{}'.format(i)
                             print w
                             print
                     print
 
                 msg = 'Test failed for {} on component {} for implementation {}.'.format(name,
-                        i, impl)
+                                                                                         i, impl)
                 raise RuntimeError(msg)
 
-
-
     def test_1d_trigonometric_float32(self, **kwds):
         self._test(dim=1, dtype=npw.float32, polynomial=False, **kwds)
+
     def test_2d_trigonometric_float32(self, **kwds):
         self._test(dim=2, dtype=npw.float32, polynomial=False, **kwds)
+
     def test_3d_trigonometric_float32(self, **kwds):
-        self._test(dim=3, dtype=npw.float32, polynomial=False, **kwds)
+        if __ENABLE_LONG_TESTS__:
+            self._test(dim=3, dtype=npw.float32, polynomial=False, **kwds)
 
     def test_1d_trigonometric_float64(self, **kwds):
         self._test(dim=1, dtype=npw.float64, polynomial=False, **kwds)
+
     def test_2d_trigonometric_float64(self, **kwds):
         self._test(dim=2, dtype=npw.float64, polynomial=False, **kwds)
+
     def test_3d_trigonometric_float64(self, **kwds):
-        self._test(dim=3, dtype=npw.float64, polynomial=False, **kwds)
+        if __ENABLE_LONG_TESTS__:
+            self._test(dim=3, dtype=npw.float64, polynomial=False, **kwds)
 
     def test_1d_polynomial_float32(self, **kwds):
         self._test(dim=1, dtype=npw.float32, polynomial=True, **kwds)
+
     def test_2d_polynomial_float32(self, **kwds):
         self._test(dim=2, dtype=npw.float32, polynomial=True, **kwds)
+
     def test_3d_polynomial_float32(self, **kwds):
-        self._test(dim=3, dtype=npw.float32, polynomial=True, **kwds)
+        if __ENABLE_LONG_TESTS__:
+            self._test(dim=3, dtype=npw.float32, polynomial=True, **kwds)
 
     def perform_tests(self):
         self.test_1d_trigonometric_float32(max_derivative=3)
@@ -358,9 +365,10 @@ class TestSpectralDerivative(object):
             # self.test_2d_polynomial_float32(max_derivative=2)
             self.test_3d_polynomial_float32(max_derivative=1)
 
+
 if __name__ == '__main__':
     TestSpectralDerivative.setup_class(enable_extra_tests=False,
-                                     enable_debug_mode=False)
+                                       enable_debug_mode=False)
 
     test = TestSpectralDerivative()
 
diff --git a/hysop/operator/tests/test_transpose.py b/hysop/operator/tests/test_transpose.py
index ce3da5a76ea3fd88d7d025805abbebf4c66c4279..1fefbf291e425903263617573fc2922363b97353 100644
--- a/hysop/operator/tests/test_transpose.py
+++ b/hysop/operator/tests/test_transpose.py
@@ -1,4 +1,3 @@
-
 import random
 from hysop.deps import np, it
 from hysop.testsenv import __ENABLE_LONG_TESTS__, __HAS_OPENCL_BACKEND__
@@ -11,32 +10,32 @@ from hysop.operator.transpose import Transpose, Implementation
 
 from hysop import Field, Box
 
+
 class TestTransposeOperator(object):
 
     @classmethod
-    def setup_class(cls, 
-            enable_extra_tests=__ENABLE_LONG_TESTS__,
-            enable_debug_mode=False):
+    def setup_class(cls,
+                    enable_extra_tests=__ENABLE_LONG_TESTS__,
+                    enable_debug_mode=False):
 
         IO.set_default_path('/tmp/hysop_tests/test_transpose')
-        
+
         if enable_debug_mode:
             cls.size_min = 3
             cls.size_max = 4
         else:
             cls.size_min = 4
             cls.size_max = 23
-        
+
         cls.enable_extra_tests = enable_extra_tests
-        cls.enable_debug_mode  = enable_debug_mode
+        cls.enable_debug_mode = enable_debug_mode
 
     @classmethod
     def teardown_class(cls):
         pass
 
-    
     def _test(self, dim, dtype, is_inplace,
-            size_min=None, size_max=None, naxes=None):
+              size_min=None, size_max=None, naxes=None):
         enable_extra_tests = self.enable_extra_tests
         assert dim > 1
 
@@ -46,7 +45,7 @@ class TestTransposeOperator(object):
         assert (((size_max-size_min+1)**dim) >= nshapes)
 
         shapes = ((size_min,)*dim,)
-        while(len(shapes)<nshapes):
+        while(len(shapes) < nshapes):
             shape = tuple(np.random.randint(low=size_min, high=size_max+1, size=dim).tolist())
             if (shape in shapes) or all((si == shape[0]) for si in shape):
                 continue
@@ -57,129 +56,129 @@ class TestTransposeOperator(object):
         all_axes = list(all_axes)
         if (naxes is not None):
             random.shuffle(all_axes)
-            all_axes = all_axes[:min(naxes,len(all_axes))]
-        
+            all_axes = all_axes[:min(naxes, len(all_axes))]
+
         if (dtype is None):
             types = [np.float32, np.float64,
                      np.complex64, np.complex128]
             dtype = types[0]
-        
+
         domain = Box(length=(1.0,)*dim)
-        for nb_components in (2,):  
-            Fin  = Field(domain=domain, name='Fin', dtype=dtype,  
-                    nb_components=nb_components, register_object=False)
-            Fout = Field(domain=domain, name='Fout', dtype=dtype, 
-                    nb_components=nb_components, register_object=False)
+        for nb_components in (2,):
+            Fin = Field(domain=domain, name='Fin', dtype=dtype,
+                        nb_components=nb_components, register_object=False)
+            Fout = Field(domain=domain, name='Fout', dtype=dtype,
+                         nb_components=nb_components, register_object=False)
             for axes in all_axes:
                 for shape in shapes:
                     self._test_one(shape=shape, axes=axes,
-                            dim=dim, dtype=dtype, is_inplace=is_inplace,
-                            domain=domain, Fin=Fin, Fout=Fout)
+                                   dim=dim, dtype=dtype, is_inplace=is_inplace,
+                                   domain=domain, Fin=Fin, Fout=Fout)
 
     @classmethod
     def __field_init(cls, data, coords, component, dtype):
         shape = data.shape
         if is_integer(dtype):
-            data[...] = np.random.random_integers(low=0, high=255, size=shape) 
+            data[...] = np.random.random_integers(low=0, high=255, size=shape)
         elif is_fp(dtype):
-            data[...] = np.random.random(size=shape) 
+            data[...] = np.random.random(size=shape)
         elif is_complex(dtype):
-            real = np.random.random(size=shape) 
-            imag = np.random.random(size=shape) 
+            real = np.random.random(size=shape)
+            imag = np.random.random(size=shape)
             data[...] = real + 1j*imag
         else:
-            msg='Unknown dtype {}.'.format(dtype)
+            msg = 'Unknown dtype {}.'.format(dtype)
             raise NotImplementedError(msg)
 
     def _test_one(self, shape, axes,
-            dim, dtype, is_inplace,
-            domain, Fin, Fout):
+                  dim, dtype, is_inplace,
+                  domain, Fin, Fout):
 
         print '\nTesting {}D Transpose: inplace={} dtype={} shape={} axes={}'.format(
-                dim, is_inplace, dtype.__name__, shape, axes)
+            dim, is_inplace, dtype.__name__, shape, axes)
         if is_inplace:
             fin, fout = Fin, Fin
-            variables = { fin: shape }
+            variables = {fin: shape}
         else:
             fin, fout = Fin, Fout
-            variables = { fin: shape, fout: shape }
+            variables = {fin: shape, fout: shape}
 
         implementations = Transpose.implementations()
         ref_impl = Implementation.PYTHON
         assert ref_impl in implementations
-       
+
         # Compute reference solution
         print '  *reference PYTHON implementation.'
         transpose = Transpose(fields=fin, output_fields=fout,
                               variables=variables, axes=axes,
                               implementation=ref_impl).build()
-        dfin  = transpose.get_input_discrete_field(fin)
+        dfin = transpose.get_input_discrete_field(fin)
         dfout = transpose.get_output_discrete_field(fout)
         dfin.initialize(self.__field_init, dtype=dtype)
-        
+
         if is_inplace:
             refin = tuple(df.copy() for df in dfin.buffers)
         else:
             refin = tuple(df for df in dfin.buffers)
 
         transpose.apply()
-        
+
         refout = tuple(df.copy() for df in dfout.buffers)
 
-        for i,(in_,out_) in enumerate(zip(refin, refout)):
+        for i, (in_, out_) in enumerate(zip(refin, refout)):
             ref = np.transpose(in_, axes=axes)
-            if (ref!=out_).any():
+            if (ref != out_).any():
                 print
                 print np.transpose(in_, axes=axes)
                 print
                 print out_
                 msg = 'Reference did not match numpy for component {}.'.format(i)
                 raise RuntimeError(msg)
-        
+
         def iter_impl(impl):
             base_kwds = dict(fields=fin, output_fields=fout, variables=variables,
-                             axes=axes, implementation=impl, 
+                             axes=axes, implementation=impl,
                              name='test_transpose_{}'.format(str(impl).lower()))
             if impl is ref_impl:
-                return 
+                return
             elif impl is Implementation.OPENCL:
                 for cl_env in iter_clenv():
-                    msg='  *platform {}, device {}'.format(cl_env.platform.name.strip(), 
-                                                          cl_env.device.name.strip())
+                    msg = '  *platform {}, device {}'.format(cl_env.platform.name.strip(),
+                                                             cl_env.device.name.strip())
                     print msg
                     yield Transpose(cl_env=cl_env, **base_kwds)
             else:
-                msg='Unknown implementation to test {}.'.format(impl)
+                msg = 'Unknown implementation to test {}.'.format(impl)
                 raise NotImplementedError(msg)
-        
+
         # Compare to other implementations
         for impl in implementations:
             for op in iter_impl(impl):
                 op = op.build()
-                dfin  = op.get_input_discrete_field(fin)
+                dfin = op.get_input_discrete_field(fin)
                 dfout = op.get_output_discrete_field(fout)
                 dfin.copy(refin)
                 op.apply()
-                out = tuple( data.get().handle for data in dfout.data )
+                out = tuple(data.get().handle for data in dfout.data)
                 self._check_output(impl, op, refin, refout, out)
-    
+
     @classmethod
     def _check_output(cls, impl, op, refin_buffers, refout_buffers, out_buffers):
         check_instance(out_buffers, tuple, values=np.ndarray)
         check_instance(refout_buffers, tuple, values=np.ndarray)
         check_instance(refin_buffers, tuple, values=np.ndarray)
 
-        for i, (out,refin,refout) in enumerate(zip(out_buffers, refin_buffers, refout_buffers)):
+        for i, (out, refin, refout) in enumerate(zip(out_buffers, refin_buffers, refout_buffers)):
             assert refout.dtype == out.dtype
             assert refout.shape == out.shape
 
             if np.all(out == refout):
                 continue
-            
+
             if cls.enable_debug_mode:
                 has_nan = np.any(np.isnan(out))
                 has_inf = np.any(np.isinf(out))
-                
+
                 print
                 print 'Test output comparisson failed for component {}:'.format(i)
                 print ' *has_nan: {}'.format(has_nan)
@@ -195,50 +194,66 @@ class TestTransposeOperator(object):
                 print out
                 print
                 print
-            
+
             msg = 'Test failed on component {} for implementation {}.'.format(i, impl)
-            raise RuntimeError(msg) 
+            raise RuntimeError(msg)
 
+    def test_2d_out_of_place(self):
+        self._test(dim=2, is_inplace=False, dtype=np.float32)
+
+    def test_3d_out_of_place(self):
+        self._test(dim=3, is_inplace=False, dtype=np.complex64)
+
+    def test_4d_out_of_place(self):
+        if __ENABLE_LONG_TESTS__:
+            self._test(dim=4, is_inplace=False, dtype=np.float64)
 
-    def test_2d_out_of_place(self, **kwds):
-        self._test(dim=2, is_inplace=False, **kwds)
-    def test_3d_out_of_place(self, **kwds):
-        self._test(dim=3, is_inplace=False, **kwds)
-    def test_4d_out_of_place(self, **kwds):
-        self._test(dim=4, is_inplace=False, **kwds)
     def test_upper_dimensions_out_of_place(self):
-        for i in xrange(5,9):
-            self._test(dim=i, dtype=None, is_inplace=False,
-                    size_min=3, size_max=4, naxes=1)
-    
-    def test_2d_inplace(self, **kwds):
-        self._test(dim=2, is_inplace=True, **kwds)
-    def test_3d_inplace(self, **kwds):
-        self._test(dim=3, is_inplace=True, **kwds)
-    def test_4d_inplace(self, **kwds):
-        self._test(dim=4, is_inplace=True, **kwds)
+        if __ENABLE_LONG_TESTS__:
+            for i in xrange(5, 9):
+                self._test(dim=i, dtype=None, is_inplace=False,
+                           size_min=3, size_max=4, naxes=1)
+
+    def test_2d_inplace(self):
+        self._test(dim=2, is_inplace=True, dtype=np.float32)
+
+    def test_3d_inplace(self):
+        self._test(dim=3, is_inplace=True, dtype=np.float32)
+
+    def test_4d_inplace(self):
+        if __ENABLE_LONG_TESTS__:
+            self._test(dim=4, is_inplace=True, dtype=np.float32)
+
     def test_upper_dimensions_inplace(self):
-        for i in xrange(5,9):
-            self._test(dim=i, dtype=None, is_inplace=True,
-                    size_min=3, size_max=4, naxes=1)
+        if __ENABLE_LONG_TESTS__:
+            for i in xrange(5, 9):
+                self._test(dim=i, dtype=None, is_inplace=True,
+                           size_min=3, size_max=4, naxes=1)
 
     def perform_tests(self):
-        self.test_2d_out_of_place(dtype=np.float32)
-        self.test_3d_out_of_place(dtype=np.complex64)
+        self._test(dim=2, is_inplace=False, dtype=np.float32)
+        self._test(dim=3, is_inplace=False, dtype=np.float64)
+        self._test(dim=3, is_inplace=False, dtype=np.complex64)
         if __ENABLE_LONG_TESTS__:
-            self.test_4d_out_of_place()
-            self.test_upper_dimensions_out_of_place()
-        
-        self.test_2d_inplace(dtype=np.float64)
-        self.test_3d_inplace(dtype=np.complex128)
+            self._test(dim=4, is_inplace=False, dtype=np.float64)
+            for i in xrange(5, 9):
+                self._test(dim=i, dtype=None, is_inplace=False,
+                           size_min=3, size_max=4, naxes=1)
+
+        self._test(dim=2, is_inplace=True, dtype=np.float32)
+        self._test(dim=3, is_inplace=True, dtype=np.float64)
+        self._test(dim=3, is_inplace=True, dtype=np.complex128)
         if __ENABLE_LONG_TESTS__:
-            self.test_4d_inplace()
-            self.test_upper_dimensions_inplace()
-    
+            self._test(dim=4, is_inplace=True, dtype=np.float32)
+            for i in xrange(5, 9):
+                self._test(dim=i, dtype=None, is_inplace=True,
+                           size_min=3, size_max=4, naxes=1)
+
+
 if __name__ == '__main__':
-    TestTransposeOperator.setup_class(enable_extra_tests=False, 
+    TestTransposeOperator.setup_class(enable_extra_tests=False,
                                       enable_debug_mode=False)
-    
+
     test = TestTransposeOperator()
     test.perform_tests()