diff --git a/lib/simd/Grid_qpx.h b/lib/simd/Grid_qpx.h
index 07933f52..bc86291d 100644
--- a/lib/simd/Grid_qpx.h
+++ b/lib/simd/Grid_qpx.h
@@ -244,7 +244,22 @@ namespace Optimization {
       return a*b;
     }
   };
-  
+
+  struct Div{
+    // Real double
+    inline vector4double operator()(vector4double a, vector4double b){
+      return vec_swdiv(a, b);
+    }
+
+    // Real float
+    FLOAT_WRAP_2(operator(), inline)
+
+    // Integer
+    inline int operator()(int a, int b){
+      return a/b;
+    }
+  };
+
   struct Conj{
     // Complex double
     inline vector4double operator()(vector4double v){
@@ -413,6 +428,7 @@ template <typename S, typename T> using ReduceSIMD = Optimization::Reduce<S,T>;
 typedef Optimization::Sum         SumSIMD;
 typedef Optimization::Sub         SubSIMD;
 typedef Optimization::Mult        MultSIMD;
+typedef Optimization::Div         DivSIMD;
 typedef Optimization::MultComplex MultComplexSIMD;
 typedef Optimization::Conj        ConjSIMD;
 typedef Optimization::TimesMinusI TimesMinusISIMD;
diff --git a/tests/Test_simd.cc b/tests/Test_simd.cc
index 189f0559..92f9bcd8 100644
--- a/tests/Test_simd.cc
+++ b/tests/Test_simd.cc
@@ -50,6 +50,12 @@ public:
   template<class vec> void operator()(vec &rr,vec &i1,vec &i2) const { rr = i1*i2;}
   std::string name(void) const { return std::string("Times"); }
 };
+class funcDivide {
+public:
+  funcDivide() {};
+  template<class vec> void operator()(vec &rr,vec &i1,vec &i2) const { rr = i1/i2;}
+  std::string name(void) const { return std::string("Divide"); }
+};
 class funcConj {
 public:
   funcConj() {};
@@ -341,6 +347,7 @@ int main (int argc, char ** argv)
   Tester<RealF,vRealF>(funcPlus());
   Tester<RealF,vRealF>(funcMinus());
   Tester<RealF,vRealF>(funcTimes());
+  Tester<RealF,vRealF>(funcDivide());
   Tester<RealF,vRealF>(funcAdj());
   Tester<RealF,vRealF>(funcConj());
   Tester<RealF,vRealF>(funcInnerProduct());
@@ -371,6 +378,7 @@ int main (int argc, char ** argv)
   Tester<RealD,vRealD>(funcPlus());
   Tester<RealD,vRealD>(funcMinus());
   Tester<RealD,vRealD>(funcTimes());
+  Tester<RealD,vRealD>(funcDivide());
   Tester<RealD,vRealD>(funcAdj());
   Tester<RealD,vRealD>(funcConj());
   Tester<RealD,vRealD>(funcInnerProduct());