diff --git a/MathLib/LinAlg/BLAS.cpp b/MathLib/LinAlg/BLAS.cpp
index bffbbe9a72257dd9b4eca9c707753fbd17797d8d..87550f22621bb9cf1493ab22f42e2a52eba642a6 100644
--- a/MathLib/LinAlg/BLAS.cpp
+++ b/MathLib/LinAlg/BLAS.cpp
@@ -19,6 +19,15 @@
 namespace MathLib { namespace BLAS
 {
 
+// Explicit specialization
+// Computes w = x/y componentwise.
+template<>
+void componentwiseDivide(Eigen::VectorXd& w,
+                         Eigen::VectorXd const& x, Eigen::VectorXd const& y)
+{
+    w.noalias() = x.cwiseQuotient(y);
+}
+
 // Explicit specialization
 // Computes the Manhattan norm of x
 template<>
@@ -90,6 +99,15 @@ void axpby(PETScVector& y, double const a, double const b, PETScVector const& x)
     VecAXPBY(*y.getRawVector(), a, b, *x.getRawVector());
 }
 
+// Explicit specialization
+// Computes w = x/y componentwise.
+template<>
+void componentwiseDivide(PETScVector& w,
+                         PETScVector const& x, PETScVector const& y)
+{
+    VecPointwiseDivide(*w.getRawVector(), *x.getRawVector(), *y.getRawVector());
+}
+
 // Explicit specialization
 // Computes the Manhattan norm of x
 template<>
@@ -224,6 +242,16 @@ void axpby(EigenVector& y, double const a, double const b, EigenVector const& x)
     y.getRawVector() = a*x.getRawVector() + b*y.getRawVector();
 }
 
+// Explicit specialization
+// Computes w = x/y componentwise.
+template<>
+void componentwiseDivide(EigenVector& w,
+                         EigenVector const& x, EigenVector const& y)
+{
+    w.getRawVector().noalias() =
+            x.getRawVector().cwiseQuotient(y.getRawVector());
+}
+
 // Explicit specialization
 // Computes the Manhattan norm of x
 template<>
diff --git a/MathLib/LinAlg/BLAS.h b/MathLib/LinAlg/BLAS.h
index 355af434f7958cb0f5dfad93670a8d0554e2a11a..a042b045446fb65f9bc7a48c3c5b896f8fa3cd94 100644
--- a/MathLib/LinAlg/BLAS.h
+++ b/MathLib/LinAlg/BLAS.h
@@ -66,6 +66,11 @@ void axpby(MatrixOrVector& y, double const a, double const b, MatrixOrVector con
     y = a*x + b*y;
 }
 
+//! Computes \f$w = x/y\f$ componentwise.
+template<typename MatrixOrVector>
+void componentwiseDivide(MatrixOrVector& w,
+                         MatrixOrVector const& x, MatrixOrVector const& y);
+
 //! Computes the Manhattan norm of \c x.
 template<typename MatrixOrVector>
 double norm1(MatrixOrVector const& x);