diff --git a/src/Integrator/BDHI/DoublyPeriodic/DPStokesSlab.cuh b/src/Integrator/BDHI/DoublyPeriodic/DPStokesSlab.cuh
index a02a22b4..aad075a4 100644
--- a/src/Integrator/BDHI/DoublyPeriodic/DPStokesSlab.cuh
+++ b/src/Integrator/BDHI/DoublyPeriodic/DPStokesSlab.cuh
@@ -169,7 +169,7 @@ else throw std::runtime_error("[DPStokesSlab] Can only average in direction X (0
       real gw;
       real tolerance;
       WallMode mode;
-      shared_ptr<BVP::BatchedBVPHandler> bvpSolver;
+      shared_ptr<BVP::BatchedBVPHandler<real>> bvpSolver;
 
     };
 
diff --git a/src/Integrator/BDHI/DoublyPeriodic/StokesSlab/Correction.cuh b/src/Integrator/BDHI/DoublyPeriodic/StokesSlab/Correction.cuh
index f6e5f6cf..c7735e3c 100644
--- a/src/Integrator/BDHI/DoublyPeriodic/StokesSlab/Correction.cuh
+++ b/src/Integrator/BDHI/DoublyPeriodic/StokesSlab/Correction.cuh
@@ -506,7 +506,7 @@ namespace uammd{
     }
 
     auto computeZeroModeBoundaryConditions(int nz, real H, WallMode mode){
-      BVP::SchurBoundaryCondition bcs(nz, H);
+      BVP::SchurBoundaryCondition<real> bcs(nz, H);
       if(mode == WallMode::bottom){
 	correction_ns::TopBoundaryConditions top(0, H);
 	correction_ns::BottomBoundaryConditions bot(0, H);
diff --git a/src/Integrator/BDHI/DoublyPeriodic/StokesSlab/initialization.cu b/src/Integrator/BDHI/DoublyPeriodic/StokesSlab/initialization.cu
index d89ed298..d16b5bf3 100644
--- a/src/Integrator/BDHI/DoublyPeriodic/StokesSlab/initialization.cu
+++ b/src/Integrator/BDHI/DoublyPeriodic/StokesSlab/initialization.cu
@@ -144,7 +144,7 @@ namespace uammd{
       auto botBC = thrust::make_transform_iterator(thrust::make_counting_iterator<int>(0), botdispatch);
       int numberSystems = (nk.x/2+1)*nk.y;
       int nz = grid.cellDim.z;
-      this->bvpSolver = std::make_shared<BVP::BatchedBVPHandler>(klist, topBC, botBC, numberSystems, halfH, nz);
+      this->bvpSolver = std::make_shared<BVP::BatchedBVPHandler<real>>(klist, topBC, botBC, numberSystems, halfH, nz);
       CudaCheckError();
     }
 
diff --git a/src/Interactor/DoublyPeriodic/PoissonSlab/BVPPoisson.cuh b/src/Interactor/DoublyPeriodic/PoissonSlab/BVPPoisson.cuh
index 161ff7ab..b1aea296 100644
--- a/src/Interactor/DoublyPeriodic/PoissonSlab/BVPPoisson.cuh
+++ b/src/Interactor/DoublyPeriodic/PoissonSlab/BVPPoisson.cuh
@@ -90,7 +90,7 @@ namespace uammd{
     }
 
     class BVPPoissonSlab{
-      std::shared_ptr<BVP::BatchedBVPHandler> bvpSolver;
+      std::shared_ptr<BVP::BatchedBVPHandler<real>> bvpSolver;
       real2 Lxy;
       real H;
       int3 cellDim;
@@ -150,7 +150,7 @@ namespace uammd{
 						   BoundaryConditionsDispatch<BottomBoundaryConditions, decltype(klist)>(klist, H));
 
 	int numberSystems = (nk.x/2+1)*nk.y;
-	this->bvpSolver = std::make_shared<BVP::BatchedBVPHandler>(klist, topBC, bottomBC, numberSystems, H, cellDim.z);
+	this->bvpSolver = std::make_shared<BVP::BatchedBVPHandler<real>>(klist, topBC, bottomBC, numberSystems, H, cellDim.z);
 	CudaCheckError();
       }
     };
diff --git a/src/misc/BoundaryValueProblem/BVPSchurComplementMatrices.cuh b/src/misc/BoundaryValueProblem/BVPSchurComplementMatrices.cuh
index 4611db81..02ddc15c 100644
--- a/src/misc/BoundaryValueProblem/BVPSchurComplementMatrices.cuh
+++ b/src/misc/BoundaryValueProblem/BVPSchurComplementMatrices.cuh
@@ -117,6 +117,7 @@ namespace uammd{
       }
     };
 
+    template <typename T>
     class SchurBoundaryCondition{
       int nz;
       real H;
@@ -125,8 +126,8 @@ namespace uammd{
       SchurBoundaryCondition(int nz, real H):nz(nz), H(H), bcs(nz){}
 
       template<class TopBC, class BottomBC>
-      std::vector<real> computeBoundaryConditionMatrix(const TopBC &top, const BottomBC &bottom){
-	std::vector<real> CandD(2*nz+4, 0);
+      std::vector<T> computeBoundaryConditionMatrix(const TopBC &top, const BottomBC &bottom){
+	std::vector<T> CandD(2*nz+4, 0);
 	auto topRow = computeTopRow(top, bottom);
 	auto bottomRow = computeBottomRow(top, bottom);
 	std::copy(topRow.begin(), topRow.end()-2, CandD.begin());
@@ -141,8 +142,8 @@ namespace uammd{
     private:
 
       template<class TopBC, class BottomBC>
-      std::vector<real> computeTopRow(const TopBC &top, const BottomBC &bottom){
-	std::vector<real> topRow(nz+2, 0);
+      std::vector<T> computeTopRow(const TopBC &top, const BottomBC &bottom){
+	std::vector<T> topRow(nz+2, 0);
 	auto tfi = bcs.topFirstIntegral();
 	auto tsi = bcs.topSecondIntegral();
 	auto tfiFactor = top.getFirstIntegralFactor();
@@ -156,8 +157,8 @@ namespace uammd{
       }
 
       template<class TopBC, class BottomBC>
-      std::vector<real> computeBottomRow(const TopBC &top, const BottomBC &bottom){
-	std::vector<real> bottomRow(nz+2, 0);
+      std::vector<T> computeBottomRow(const TopBC &top, const BottomBC &bottom){
+	std::vector<T> bottomRow(nz+2, 0);
 	auto bfi = bcs.bottomFirstIntegral();
 	auto bsi = bcs.bottomSecondIntegral();
 	auto bfiFactor = bottom.getFirstIntegralFactor();
@@ -172,10 +173,11 @@ namespace uammd{
 
     };
 
-    std::vector<real> computeSecondIntegralMatrix(real k, real H, int nz){
-      std::vector<real> A(nz*nz, 0);
+    template<class T>
+    std::vector<T> computeSecondIntegralMatrix(T k, real H, int nz){
+      std::vector<T> A(nz*nz, 0);
       SecondIntegralMatrix sim(nz);
-      real kH2 = k*k*H*H;
+      T kH2 = k*k*H*H;
       fori(0, nz){
 	forj(0,nz){
 	  A[i+nz*j] = (i==j) - kH2*sim.getElement(i, j);
@@ -184,9 +186,10 @@ namespace uammd{
       return std::move(A);
     }
 
-    std::vector<real> computeInverseSecondIntegralMatrix(real k, real H, int nz){
+    template<class T>
+    std::vector<T> computeInverseSecondIntegralMatrix(T k, real H, int nz){
       if(k==0){
-	std::vector<real> invA(nz*nz, 0);
+	std::vector<T> invA(nz*nz, 0);
 	fori(0, nz){
 	  invA[i+nz*i] = 1;
 	}
diff --git a/src/misc/BoundaryValueProblem/BVPSolver.cuh b/src/misc/BoundaryValueProblem/BVPSolver.cuh
index 4e54b094..d2444017 100644
--- a/src/misc/BoundaryValueProblem/BVPSolver.cuh
+++ b/src/misc/BoundaryValueProblem/BVPSolver.cuh
@@ -42,39 +42,40 @@ namespace uammd{
   namespace BVP{
 
     namespace detail{
+      template <typename U>
       class SubsystemSolver{
 	int nz;
 	real H;
-	StorageHandle<real> CinvA_storage;
-	StorageHandle<real4> CinvABmD_storage;
+	StorageHandle<U> CinvA_storage;
+	StorageHandle<U> CinvABmD_storage;
 
       public:
 
 	SubsystemSolver(int nz, real H):nz(nz), H(H){}
 
 	void registerRequiredStorage(StorageRegistration &memoryManager){
-	  CinvA_storage =    memoryManager.registerStorageRequirement<real>(2*nz+2);
-	  CinvABmD_storage = memoryManager.registerStorageRequirement<real4>(1);
+	  CinvA_storage =    memoryManager.registerStorageRequirement<U>(2*nz+2);
+	  CinvABmD_storage = memoryManager.registerStorageRequirement<U>(4);
 	}
 
 	template<class TopBC, class BottomBC>
-	void precompute(real k, const TopBC &top, const BottomBC &bottom, StorageRetriever &memoryManager){
+	void precompute(U k, const TopBC &top, const BottomBC &bottom, StorageRetriever &memoryManager){
 	  auto CinvA_it = memoryManager.retrieveStorage(CinvA_storage);
 	  auto CinvABmD_it = memoryManager.retrieveStorage(CinvABmD_storage);
 	  auto invA = computeInverseSecondIntegralMatrix(k, H, nz);
-	  SchurBoundaryCondition bcs(nz, H);
+	  SchurBoundaryCondition<U> bcs(nz, H);
 	  auto CandD = bcs.computeBoundaryConditionMatrix(top, bottom);
-	  real4 D = make_real4(CandD[2*nz], CandD[2*nz+1], CandD[2*nz+2], CandD[2*nz+3]);
+	  U D[4] = {CandD[2*nz], CandD[2*nz+1], CandD[2*nz+2], CandD[2*nz+3]};
 	  auto CinvA = matmul(CandD, nz, 2, invA, nz, nz);
 	  std::copy(CinvA.begin(), CinvA.end(), CinvA_it);
-	  real4 CinvAB;
-	  real B00 = -k*k*H*H;
-	  real B11 = -k*k*H*H;
-	  CinvAB.x = CinvA[0]*B00;
-	  CinvAB.y = CinvA[1]*B11;
-	  CinvAB.z = CinvA[0+nz]*B00;
-	  CinvAB.w = CinvA[1+nz]*B11;
-	  CinvABmD_it[0] = CinvAB - D;
+	  U CinvAB[4];
+	  U B00 = -k*k*H*H;
+	  U B11 = -k*k*H*H;
+	  CinvAB[0] = CinvA[0]*B00;
+	  CinvAB[1] = CinvA[1]*B11;
+	  CinvAB[2] = CinvA[0+nz]*B00;
+	  CinvAB[3] = CinvA[1+nz]*B11;
+	  fori(0, 4) CinvABmD_it[i] = CinvAB[i] - D[i];
 	}
 
 	template<class T, class FnIterator>
@@ -83,16 +84,19 @@ namespace uammd{
 					   StorageRetriever &memoryManager){
 	  const auto CinvA = memoryManager.retrieveStorage(CinvA_storage);
 	  const auto CinvAfmab = computeRightHandSide(alpha, beta, fn, CinvA);
-	  const real4 CinvABmD = *(memoryManager.retrieveStorage(CinvABmD_storage));
+	  const auto CinvABmD_it = memoryManager.retrieveStorage(CinvABmD_storage);
+	  U CinvABmD[4] = {CinvABmD_it[0], CinvABmD_it[1],
+				       CinvABmD_it[2], CinvABmD_it[3]};
+
 	  const auto c0d0 = solveSubsystem(CinvABmD, CinvAfmab);
 	  return c0d0;
 	}
 
       private:
 
-	template<class T>
-	__device__ thrust::pair<T,T> solveSubsystem(real4 CinvABmD, thrust::pair<T,T> CinvAfmab) const{
-	  auto c0d0 = solve2x2System(CinvABmD, CinvAfmab);
+	template<class T, class T2>
+	__device__ thrust::pair<T,T> solveSubsystem(T2 CinvABmD[4], thrust::pair<T,T> CinvAfmab) const{
+	  auto c0d0 = solve2x2System<T,T2>(CinvABmD, CinvAfmab);
 	  return c0d0;
 	}
 
@@ -113,10 +117,11 @@ namespace uammd{
 
       };
 
+      template<typename U>
       class PentadiagonalSystemSolver{
 	int nz;
 	real H;
-	KBPENTA_mod pentasolve;
+	KBPENTA_mod<U> pentasolve;
       public:
 
 	PentadiagonalSystemSolver(int nz, real H):
@@ -126,12 +131,12 @@ namespace uammd{
 	  pentasolve.registerRequiredStorage(memoryManager);
 	}
 
-	void precompute(real k, StorageRetriever &memoryManager){
-	  real diagonal[nz];
-	  real diagonal_p2[nz]; diagonal_p2[nz-2] = diagonal_p2[nz-1] = 0;
-	  real diagonal_m2[nz]; diagonal_m2[0] = diagonal_m2[1] = 0;
+	void precompute(U k, StorageRetriever &memoryManager){
+	  U diagonal[nz];
+	  U diagonal_p2[nz]; diagonal_p2[nz-2] = diagonal_p2[nz-1] = 0;
+	  U diagonal_m2[nz]; diagonal_m2[0] = diagonal_m2[1] = 0;
 	  SecondIntegralMatrix sim(nz);
-	  const real kH2 = k*k*H*H;
+	  const U kH2 = k*k*H*H;
 	  for(int i = 0; i<nz; i++){
 	    diagonal[i] = 1.0 - kH2*sim.getElement(i,i);
 	    if(i<nz-2) diagonal_p2[i] = -kH2*sim.getElement(i+2, i);
@@ -148,23 +153,24 @@ namespace uammd{
       };
     }
 
+    template<typename U>
     class BoundaryValueProblemSolver{
-      detail::PentadiagonalSystemSolver pent;
-      detail::SubsystemSolver sub;
-      StorageHandle<real> waveVector;
+      detail::PentadiagonalSystemSolver<U> pent;
+      detail::SubsystemSolver<U> sub;
+      StorageHandle<U> waveVector;
       int nz;
       real H;
     public:
       BoundaryValueProblemSolver(int nz, real H): nz(nz), H(H), sub(nz, H), pent(nz, H){}
 
       void registerRequiredStorage(StorageRegistration &mem){
-	waveVector = mem.registerStorageRequirement<real>(1);
+	waveVector = mem.registerStorageRequirement<U>(1);
 	sub.registerRequiredStorage(mem);
 	pent.registerRequiredStorage(mem);
       }
 
       template<class TopBC, class BottomBC>
-      void precompute(StorageRetriever &mem, real k, const TopBC &top, const BottomBC &bot){
+      void precompute(StorageRetriever &mem, U k, const TopBC &top, const BottomBC &bot){
 	auto k_access = mem.retrieveStorage(waveVector);
 	k_access[0] = k;
 	pent.precompute(k, mem);
@@ -177,10 +183,10 @@ namespace uammd{
 			    AnIterator& an,
 			    CnIterator& cn,
 			    StorageRetriever &mem){
-        const real k = *(mem.retrieveStorage(waveVector));
+        const auto k = *(mem.retrieveStorage(waveVector));
 	T c0, d0;
 	thrust::tie(c0, d0) = sub.solve(fn, alpha, beta, mem);
-	const real kH2 = k*k*H*H;
+	const auto kH2 = k*k*H*H;
 	fn[0] += kH2*c0;
 	fn[1] += kH2*d0;
 	pent.solve(fn, an, mem);
@@ -201,15 +207,17 @@ namespace uammd{
 
     };
 
+    template<typename U>
     class BatchedBVPHandler;
 
+    template<typename U>
     struct BatchedBVPGPUSolver{
     private:
       int numberSystems;
-      BoundaryValueProblemSolver bvpSolver;
+      BoundaryValueProblemSolver<U> bvpSolver;
       char* gpuMemory;
-      friend class BatchedBVPHandler;
-      BatchedBVPGPUSolver(int numberSystems, BoundaryValueProblemSolver bvpSolver, char *raw):
+      friend class BatchedBVPHandler<U>;
+      BatchedBVPGPUSolver(int numberSystems, BoundaryValueProblemSolver<U> bvpSolver, char *raw):
 	numberSystems(numberSystems), bvpSolver(bvpSolver), gpuMemory(raw){}
     public:
 
@@ -225,9 +233,10 @@ namespace uammd{
 
     };
 
+    template <typename U>
     class BatchedBVPHandler{
       int numberSystems;
-      BoundaryValueProblemSolver bvp;
+      BoundaryValueProblemSolver<U> bvp;
       thrust::device_vector<char> gpuMemory;
     public:
 
@@ -240,9 +249,9 @@ namespace uammd{
 	precompute(klist, top, bot);
       }
 
-      BatchedBVPGPUSolver getGPUSolver(){
+      BatchedBVPGPUSolver<U> getGPUSolver(){
 	auto raw = thrust::raw_pointer_cast(gpuMemory.data());
-	BatchedBVPGPUSolver d_solver(numberSystems, bvp, raw);
+	BatchedBVPGPUSolver<U> d_solver(numberSystems, bvp, raw);
 	return d_solver;
       }
 
diff --git a/src/misc/BoundaryValueProblem/KBPENTA.cuh b/src/misc/BoundaryValueProblem/KBPENTA.cuh
index 636b217b..02873496 100644
--- a/src/misc/BoundaryValueProblem/KBPENTA.cuh
+++ b/src/misc/BoundaryValueProblem/KBPENTA.cuh
@@ -6,20 +6,21 @@ namespace uammd{
   namespace BVP{
 
 //Algorithm adapted from  http://dx.doi.org/10.1080/00207160802326507 for a special case of only three diagonals being non zero
+    template<typename U>
     class KBPENTA_mod{
-      StorageHandle<real> storageHandle;
+      StorageHandle<U> storageHandle;
       int nz;
     public:
 
       KBPENTA_mod(int nz): nz(nz){}
 
       void registerRequiredStorage(StorageRegistration &memoryManager){
-	storageHandle =  memoryManager.registerStorageRequirement<real>(3*nz+2);
+	storageHandle =  memoryManager.registerStorageRequirement<U>(3*nz+2);
       }
 
-      void store(real *diagonal, real *diagonal_p2, real *diagonal_m2, StorageRetriever &memoryManager){
+      void store(U *diagonal, U *diagonal_p2, U *diagonal_m2, StorageRetriever &memoryManager){
 	auto storage = memoryManager.retrieveStorage(storageHandle);
-	std::vector<real> beta(nz+1, 0);
+	std::vector<U> beta(nz+1, 0);
 	beta[0] = 0;
 	beta[1] = diagonal[nz-nz];
 	beta[2] = diagonal[nz-(nz-1)];
diff --git a/src/misc/BoundaryValueProblem/MatrixUtils.h b/src/misc/BoundaryValueProblem/MatrixUtils.h
index 7a1689af..282910ce 100644
--- a/src/misc/BoundaryValueProblem/MatrixUtils.h
+++ b/src/misc/BoundaryValueProblem/MatrixUtils.h
@@ -1,5 +1,5 @@
 /*Raul P. Pelaez 2019-2021. Boundary Value Problem Matrix algebra utilities
- */
+*/
 #ifndef BVP_MATRIX_UTILS_H
 #define BVP_MATRIX_UTILS_H
 #include "utils/cufftPrecisionAgnostic.h"
@@ -12,95 +12,129 @@
 #else
 #include<lapacke.h>
 #endif
-
+#include<thrust/complex.h>
 namespace uammd{
   namespace BVP{
-
+    using complex = thrust::complex<real>;
     template<class T> struct LapackeUAMMD;
     template<>
-    struct LapackeUAMMD<float>{
-      template<class ...T> static int getrf(T... args){return LAPACKE_sgetrf(args...);}
-      template<class ...T> static int getri(T... args){return LAPACKE_sgetri(args...);}
-    };
+      struct LapackeUAMMD<float>{
+        template<class ...T> static int getrf(T... args){return LAPACKE_sgetrf(args...);}
+        template<class ...T> static int getri(T... args){return LAPACKE_sgetri(args...);}
+      };
     template<>
-    struct LapackeUAMMD<double>{
-      template<class ...T> static int getrf(T... args){return LAPACKE_dgetrf(args...);}
-      template<class ...T> static int getri(T... args){return LAPACKE_dgetri(args...);}
-    };
+      struct LapackeUAMMD<double>{
+        template<class ...T> static int getrf(T... args){return LAPACKE_dgetrf(args...);}
+        template<class ...T> static int getri(T... args){return LAPACKE_dgetri(args...);}
+      };
 
     template<>
-    struct LapackeUAMMD<lapack_complex_float>{
-      template<class ...T> static int getrf(T... args){return LAPACKE_cgetrf(args...);}
-      template<class ...T> static int getri(T... args){return LAPACKE_cgetri(args...);}
-    };
+      struct LapackeUAMMD<lapack_complex_float>{
+        template<class ...T> static int getrf(T... args){return LAPACKE_cgetrf(args...);}
+        template<class ...T> static int getri(T... args){return LAPACKE_cgetri(args...);}
+      };
     template<>
-    struct LapackeUAMMD<lapack_complex_double>{
-      template<class ...T> static int getrf(T... args){return LAPACKE_zgetrf(args...);}
-      template<class ...T> static int getri(T... args){return LAPACKE_zgetri(args...);}
-    };
+      struct LapackeUAMMD<lapack_complex_double>{
+        template<class ...T> static int getrf(T... args){return LAPACKE_zgetrf(args...);}
+        template<class ...T> static int getri(T... args){return LAPACKE_zgetri(args...);}
+      };
     template<>
-    struct LapackeUAMMD<cufftComplex_t<float>>{
-      static int getrf(int matrix_layout, lapack_int n, lapack_int m,
-		       cufftComplex_t<float>* a, lapack_int lda,
-		       lapack_int* ipiv){
-	return LAPACKE_cgetrf(matrix_layout, n, m, (lapack_complex_float*)(a), lda, ipiv);
-      }
-      static int getri(int matrix_layout, lapack_int n,
-		       cufftComplex_t<float>* a, lapack_int lda,
-		       lapack_int* ipiv){
-	return LAPACKE_cgetri(matrix_layout, n, (lapack_complex_float*)(a), lda, ipiv);
-      }
-    };
+      struct LapackeUAMMD<cufftComplex_t<float>>{
+        static int getrf(int matrix_layout, lapack_int n, lapack_int m,
+            cufftComplex_t<float>* a, lapack_int lda,
+            lapack_int* ipiv){
+          return LAPACKE_cgetrf(matrix_layout, n, m, (lapack_complex_float*)(a), lda, ipiv);
+        }
+        static int getri(int matrix_layout, lapack_int n,
+            cufftComplex_t<float>* a, lapack_int lda,
+            lapack_int* ipiv){
+          return LAPACKE_cgetri(matrix_layout, n, (lapack_complex_float*)(a), lda, ipiv);
+        }
+      };
     template<>
-    struct LapackeUAMMD<cufftComplex_t<double>>{
-      static int getrf(int matrix_layout, lapack_int n, lapack_int m,
-		       cufftComplex_t<double>* a, lapack_int lda,
-		       lapack_int* ipiv){
-	return LAPACKE_zgetrf(matrix_layout, n, m, (lapack_complex_double*)(a), lda, ipiv);
-      }
-      static int getri(int matrix_layout, lapack_int n,
-		       cufftComplex_t<double>* a, lapack_int lda,
-		       lapack_int* ipiv){
-	return LAPACKE_zgetri(matrix_layout, n, (lapack_complex_double*)(a), lda, ipiv);
-      }
-    };
+      struct LapackeUAMMD<cufftComplex_t<double>>{
+        static int getrf(int matrix_layout, lapack_int n, lapack_int m,
+            cufftComplex_t<double>* a, lapack_int lda,
+            lapack_int* ipiv){
+          return LAPACKE_zgetrf(matrix_layout, n, m, (lapack_complex_double*)(a), lda, ipiv);
+        }
+        static int getri(int matrix_layout, lapack_int n,
+            cufftComplex_t<double>* a, lapack_int lda,
+            lapack_int* ipiv){
+          return LAPACKE_zgetri(matrix_layout, n, (lapack_complex_double*)(a), lda, ipiv);
+        }
+      };
+    template<>
+      struct LapackeUAMMD<thrust::complex<float>>{
+        static int getrf(int matrix_layout, lapack_int n, lapack_int m,
+            thrust::complex<float>* a, lapack_int lda,
+            lapack_int* ipiv){
+          return LAPACKE_cgetrf(matrix_layout, n, m, (lapack_complex_float*)(a), lda, ipiv);
+        }
+        static int getri(int matrix_layout, lapack_int n,
+            thrust::complex<float>* a, lapack_int lda,
+            lapack_int* ipiv){
+          return LAPACKE_cgetri(matrix_layout, n, (lapack_complex_float*)(a), lda, ipiv);
+        }
+      };
+    template<>
+      struct LapackeUAMMD<thrust::complex<double>>{
+        static int getrf(int matrix_layout, lapack_int n, lapack_int m,
+            thrust::complex<double>* a, lapack_int lda,
+            lapack_int* ipiv){
+          return LAPACKE_zgetrf(matrix_layout, n, m, (lapack_complex_double*)(a), lda, ipiv);
+        }
+        static int getri(int matrix_layout, lapack_int n,
+            thrust::complex<double>* a, lapack_int lda,
+            lapack_int* ipiv){
+          return LAPACKE_zgetri(matrix_layout, n, (lapack_complex_double*)(a), lda, ipiv);
+        }
+      };
 
 
     template<class T>
-    std::vector<T> invertSquareMatrix(const std::vector<T> &A, lapack_int N){
-      lapack_int pivotArray[N];
-      int errorHandler;
-      auto invA = A;
-      errorHandler = LapackeUAMMD<T>::getrf(LAPACK_ROW_MAJOR, N, N, invA.data(), N, pivotArray);
-      if(errorHandler){
-	throw std::runtime_error("Lapacke getrf failed with error code: " + std::to_string(errorHandler));
+      std::vector<T> invertSquareMatrix(const std::vector<T> &A, lapack_int N){
+        lapack_int pivotArray[N];
+        int errorHandler;
+        auto invA = A;
+        errorHandler = LapackeUAMMD<T>::getrf(LAPACK_ROW_MAJOR, N, N, invA.data(), N, pivotArray);
+        if(errorHandler){
+          throw std::runtime_error("Lapacke getrf failed with error code: " + std::to_string(errorHandler));
+        }
+        errorHandler = LapackeUAMMD<T>::getri(LAPACK_ROW_MAJOR, N, invA.data(), N, pivotArray);
+        if(errorHandler){
+          throw std::runtime_error("Lapacke getri failed with error code: " + std::to_string(errorHandler));
+        }
+        return invA;
       }
-      errorHandler = LapackeUAMMD<T>::getri(LAPACK_ROW_MAJOR, N, invA.data(), N, pivotArray);
-      if(errorHandler){
-	throw std::runtime_error("Lapacke getri failed with error code: " + std::to_string(errorHandler));
-      }
-      return invA;
-    }
 
     template<class T, class T2>
-    std::vector<typename T::value_type> matmul(const T &A, int ncol_a, int nrow_a,
-					       const T2 &B, int ncol_b, int nrow_b){
-      std::vector<typename T::value_type> C;
-      C.resize(ncol_b*nrow_a);
-      for(int i = 0; i<nrow_a; i++){
-	for(int j = 0; j<ncol_b; j++){
-	  real tmp = 0;
-	  for(int k = 0; k<ncol_a; k++){
-	    tmp += A[k+ncol_a*i] * B [j+ncol_b*k];
-	  }
-	  C[j+ncol_b*i] = tmp;
-	}
+      auto matmul(const T &A, int ncol_a, int nrow_a,
+          const T2 &B, int ncol_b, int nrow_b){
+        using type = typename T::value_type;
+        std::vector<type> C;
+        C.resize(ncol_b*nrow_a);
+        for(int i = 0; i<nrow_a; i++){
+          for(int j = 0; j<ncol_b; j++){
+            type tmp{};
+            for(int k = 0; k<ncol_a; k++){
+              tmp += A[k+ncol_a*i] * B [j+ncol_b*k];
+            }
+            C[j+ncol_b*i] = tmp;
+          }
+        }
+        return C;
       }
-      return C;
-    }
-
     //Solves the linear system A*x = b
-    //Given A as a real4 (2x2 matrix) and b as two numbers of an arbitrary type (could be complex, real...)
+    //Given A (2x2 matrix) and b as two numbers of an arbitrary type (could be complex, real...)
+    template<class T, class T2>
+    __device__ thrust::pair<T,T> solve2x2System(T2 A[4], thrust::pair<T,T> b){
+      const auto det  = A[0]*A[3] - A[1]*A[2];
+      const T c0 = (b.first*A[3] - A[1]*b.second)/det;
+      const T d0 = (b.second*A[0] - b.first*A[2])/det;
+      return thrust::make_pair(c0, d0);
+    }
+    
     template<class T>
     __device__ thrust::pair<T,T> solve2x2System(real4 A, thrust::pair<T,T> b){
       const real det  = A.x*A.w - A.y*A.z;
diff --git a/src/utils/vector.cuh b/src/utils/vector.cuh
index df4ed98c..43a4330f 100644
--- a/src/utils/vector.cuh
+++ b/src/utils/vector.cuh
@@ -976,4 +976,339 @@ namespace uammd{
   }
 }
 
+namespace uammd{
+////////////////////////////////////VEC2///////////////////////////////////////
+  template<class T>
+  struct vec2{
+    T x,y;
+  };
+
+  template<class T>
+  vec2<T> make_vec2(T x, T y){ return {x, y};}
+
+  template<class T>
+  VECATTR  vec2<T> operator +(const vec2<T> &a, const vec2<T> &b){
+    return {a.x + b.x, a.y + b.y};
+  }
+
+  template<class T>
+  VECATTR  vec2<T> operator +(const vec2<T> &a, const T &b){
+    return {a.x + b, a.y + b};
+  }
+
+  template<class T>
+  VECATTR  vec2<T> operator -(const vec2<T> &a, const vec2<T> &b){
+    return {a.x - b.x, a.y - b.y};
+  }
+
+
+  template<class T>
+  VECATTR  vec2<T> operator -(const vec2<T> &a, const T &b){
+    return {a.x - b, a.y - b};
+  }
+
+  template<class T>
+  VECATTR  vec2<T> operator -(const T &b, const vec2<T> &a){
+    return {b-a.x, b-a.y};
+  }
+
+
+  template<class T>
+  VECATTR  vec2<T> operator *(const vec2<T> &a, const vec2<T> &b){
+    return {a.x * b.x, a.y * b.y};
+  }
+
+
+  template<class T>
+  VECATTR  vec2<T> operator *(const vec2<T> &a, const T &b){
+    return {a.x * b, a.y * b};
+  }
+
+  template<class T>
+  VECATTR  vec2<T> operator /(const vec2<T> &a, const vec2<T> &b){
+    return {a.x / b.x, a.y / b.y};
+  }
+
+  template<class T>
+  VECATTR  vec2<T> operator /(const vec2<T> &a, const T &b){
+    return {a.x/b, a.y/b};
+  }
+
+
+  template<class T>
+  VECATTR  vec2<T> operator /(const T &b, const vec2<T> &a){
+    return {b / a.x, b / a.y};
+  }
+
+  template<class T>
+  VECATTR  void operator +=(vec2<T> &a, const vec2<T> &b){
+    a = a + b;
+  }
+
+  template<class T>
+  VECATTR  vec2<T> operator +(const T &b, const vec2<T> &a){return a+b;}
+
+  template<class T>
+  VECATTR  void operator +=(vec2<T> &a, const T &b){
+    a = a+b;
+  }
+
+  template<class T>
+  VECATTR  void operator -=(vec2<T> &a, const vec2<T> &b){
+    a = a-b;
+  }
+  template<class T>
+  VECATTR  void operator -=(vec2<T> &a, const T &b){
+    a=a-b;
+  }
+  template<class T>
+  VECATTR  void operator *=(vec2<T> &a, const vec2<T> &b){
+    a = a*b;
+  }
+
+  template<class T>
+  VECATTR  vec2<T> operator *(const T &b, const vec2<T> &a){
+    return a*b;
+  }
+
+  template<class T>
+  VECATTR  void operator *=(vec2<T> &a, const T &b){
+    a=a*b;
+  }
+
+  template<class T>
+  VECATTR  void operator /=(vec2<T> &a, const vec2<T> &b){
+    a = a/b;
+  }
+
+  template<class T>
+  VECATTR  void operator /=(vec2<T> &a, const T &b){
+    a = a/b;
+  }
+
+////////////////////////////////////VEC3///////////////////////////////////////
+  template<class T>
+  struct vec3{
+    T x,y,z;
+  };
+
+  template<class T>
+  vec3<T> make_vec3(T x, T y, T z){ return {x, y, z};}
+
+  template<class T>
+  VECATTR  vec3<T> operator +(const vec3<T> &a, const vec3<T> &b){
+    return {a.x + b.x, a.y + b.y, a.z + b.z};
+  }
+
+  template<class T>
+  VECATTR  vec3<T> operator +(const vec3<T> &a, const T &b){
+    return {a.x + b, a.y + b, a.z + b};
+  }
+
+  template<class T>
+  VECATTR  vec3<T> operator -(const vec3<T> &a, const vec3<T> &b){
+    return {a.x - b.x, a.y - b.y, a.z - b.z};
+  }
+
+
+  template<class T>
+  VECATTR  vec3<T> operator -(const vec3<T> &a, const T &b){
+    return {a.x - b, a.y - b, a.z - b};
+  }
+
+  template<class T>
+  VECATTR  vec3<T> operator -(const T &b, const vec3<T> &a){
+    return {b-a.x, b-a.y, b-a.z};
+  }
+
+
+  template<class T>
+  VECATTR  vec3<T> operator *(const vec3<T> &a, const vec3<T> &b){
+    return {a.x * b.x, a.y * b.y, a.z * b.z};
+  }
+
+
+  template<class T>
+  VECATTR  vec3<T> operator *(const vec3<T> &a, const T &b){
+    return {a.x * b, a.y * b, a.z * b};
+  }
+
+  template<class T>
+  VECATTR  vec3<T> operator /(const vec3<T> &a, const vec3<T> &b){
+    return {a.x / b.x, a.y / b.y, a.z / b.z};
+  }
+
+  template<class T>
+  VECATTR  vec3<T> operator /(const vec3<T> &a, const T &b){
+    return {a.x/b, a.y/b, a.z/b};
+  }
+
+
+  template<class T>
+  VECATTR  vec3<T> operator /(const T &b, const vec3<T> &a){
+    return {b / a.x, b / a.y, b / a.z};
+  }
+
+
+
+  template<class T>
+  VECATTR  void operator +=(vec3<T> &a, const vec3<T> &b){
+    a = a + b;
+  }
+
+  template<class T>
+  VECATTR  vec3<T> operator +(const T &b, const vec3<T> &a){return a+b;}
+
+  template<class T>
+  VECATTR  void operator +=(vec3<T> &a, const T &b){
+    a = a+b;
+  }
+
+  template<class T>
+  VECATTR  void operator -=(vec3<T> &a, const vec3<T> &b){
+    a = a-b;
+  }
+  template<class T>
+  VECATTR  void operator -=(vec3<T> &a, const T &b){
+    a=a-b;
+  }
+  template<class T>
+  VECATTR  void operator *=(vec3<T> &a, const vec3<T> &b){
+    a = a*b;
+  }
+
+  template<class T>
+  VECATTR  vec3<T> operator *(const T &b, const vec3<T> &a){
+    return a*b;
+  }
+
+  template<class T>
+  VECATTR  void operator *=(vec3<T> &a, const T &b){
+    a=a*b;
+  }
+
+  template<class T>
+  VECATTR  void operator /=(vec3<T> &a, const vec3<T> &b){
+    a = a/b;
+  }
+
+  template<class T>
+  VECATTR  void operator /=(vec3<T> &a, const T &b){
+    a = a/b;
+  }
+
+  /////////////////////////////////VEC4//////////////////////////
+
+  template<class T>
+  struct vec4{
+    T x,y,z,w;
+  };
+
+
+
+  template<class T>
+  vec4<T> make_vec4(T x, T y, T z, T w){ return {x, y, z, w};}
+
+  template<class T>
+  VECATTR  vec4<T> operator +(const vec4<T> &a, const vec4<T> &b){
+    return {a.x + b.x, a.y + b.y, a.z + b.z, a.w + b.w};
+  }
+
+  template<class T>
+  VECATTR  vec4<T> operator +(const vec4<T> &a, const T &b){
+    return {a.x + b, a.y + b, a.z + b, a.w + b};
+  }
+
+  template<class T>
+  VECATTR  vec4<T> operator -(const vec4<T> &a, const vec4<T> &b){
+    return {a.x - b.x, a.y - b.y, a.z - b.z, a.w - b.w};
+  }
+
+
+  template<class T>
+  VECATTR  vec4<T> operator -(const vec4<T> &a, const T &b){
+    return {a.x - b, a.y - b, a.z - b, a.w - b};
+  }
+
+  template<class T>
+  VECATTR  vec4<T> operator -(const T &b, const vec4<T> &a){
+    return {b-a.x, b-a.y, b-a.z, b-a.w};
+  }
+
+
+  template<class T>
+  VECATTR  vec4<T> operator *(const vec4<T> &a, const vec4<T> &b){
+    return {a.x * b.x, a.y * b.y, a.z * b.z, a.w * b.w};
+  }
+
+
+  template<class T>
+  VECATTR  vec4<T> operator *(const vec4<T> &a, const T &b){
+    return {a.x * b, a.y * b, a.z * b, a.w * b};
+  }
+
+  template<class T>
+  VECATTR  vec4<T> operator /(const vec4<T> &a, const vec4<T> &b){
+    return {a.x / b.x, a.y / b.y, a.z / b.z, a.w / b.w};
+  }
+
+  template<class T>
+  VECATTR  vec4<T> operator /(const vec4<T> &a, const T &b){
+    return {a.x/b, a.y/b, a.z/b, a.w/b};
+  }
+
+
+  template<class T>
+  VECATTR  vec4<T> operator /(const T &b, const vec4<T> &a){
+    return {b / a.x, b / a.y, b / a.z, b / a.w};
+  }
+
+
+
+  template<class T>
+  VECATTR  void operator +=(vec4<T> &a, const vec4<T> &b){
+    a = a + b;
+  }
+
+  template<class T>
+  VECATTR  vec4<T> operator +(const T &b, const vec4<T> &a){return a+b;}
+
+  template<class T>
+  VECATTR  void operator +=(vec4<T> &a, const T &b){
+    a = a+b;
+  }
+
+  template<class T>
+  VECATTR  void operator -=(vec4<T> &a, const vec4<T> &b){
+    a = a-b;
+  }
+  template<class T>
+  VECATTR  void operator -=(vec4<T> &a, const T &b){
+    a=a-b;
+  }
+  template<class T>
+  VECATTR  void operator *=(vec4<T> &a, const vec4<T> &b){
+    a = a*b;
+  }
+
+  template<class T>
+  VECATTR  vec4<T> operator *(const T &b, const vec4<T> &a){
+    return a*b;
+  }
+
+  template<class T>
+  VECATTR  void operator *=(vec4<T> &a, const T &b){
+    a=a*b;
+  }
+
+  template<class T>
+  VECATTR  void operator /=(vec4<T> &a, const vec4<T> &b){
+    a = a/b;
+  }
+
+  template<class T>
+  VECATTR  void operator /=(vec4<T> &a, const T &b){
+    a = a/b;
+  }
+}
 #endif