diff --git a/integrators/_details/interfaces/tpsPreconditioner.py b/integrators/_details/interfaces/tpsPreconditioner.py
index c8df2a7b7d21cfee2deaf1dd27fcf32bc21b0e35..9e7003f5a2cb204f0da408a6a8cbe97e129a6498 100644
--- a/integrators/_details/interfaces/tpsPreconditioner.py
+++ b/integrators/_details/interfaces/tpsPreconditioner.py
@@ -123,8 +123,8 @@ class TPSPreconditioner(_preconditionerInterface._PreconditionerInterface):
     pre = ml.aspreconditioner()
     pre1D = lambda z: pre * z
   
-    self._prePractical = lambda y: np.concatenate((
-      pre1D(y[:y.size//3]), pre1D(y[y.size//3:2*y.size//3]), pre1D(y[2*y.size//3:])))
+    self._prePractical = lambda y: \
+      pre1D(y.reshape((y.size//3,3), order='F')).reshape(y.size, order='F')
 
     practical2D = lambda x: \
       self._Q.dot(self._prePractical(self._Q.transpose().dot(x)))
@@ -154,7 +154,7 @@ class TPSPreconditioner(_preconditionerInterface._PreconditionerInterface):
     pre1D = lambda z: pre * z
 
     stationary2D = lambda x: \
-      np.vstack((pre1D(x[0::2]), pre1D(x[1::2]))).reshape((-1,),order='F')
+      pre1D(x.reshape((x.size//2,2), order='C')).reshape(x.size, order='C')
     
     self._preconditioner = scipy.sparse.linalg.LinearOperator( \
       (A.shape[0]*2, A.shape[0]*2), matvec=stationary2D)