diff --git a/MathLib/LinAlg/BLAS.cpp b/MathLib/LinAlg/BLAS.cpp
index 93813d361386587e5e37b1169e2ea004285b4fa3..bf2d8e3aba5c0f2b3b03255b54dfff822c2c1d90 100644
--- a/MathLib/LinAlg/BLAS.cpp
+++ b/MathLib/LinAlg/BLAS.cpp
@@ -10,6 +10,27 @@
 #include "BLAS.h"
 
 
+// Dense Eigen matrices and vectors ////////////////////////////////////////
+#ifdef OGS_USE_EIGEN
+
+#include <Eigen/Core>
+
+namespace MathLib { namespace BLAS
+{
+
+// Explicit specialization
+// Computes the Euclidean norm of x
+template<>
+double norm2(Eigen::VectorXd const& x)
+{
+    return x.norm();
+}
+
+} } // namespaces
+
+#endif
+
+
 // Global PETScMatrix/PETScVector //////////////////////////////////////////
 #ifdef USE_PETSC
 
@@ -52,6 +73,14 @@ void axpby(PETScVector& y, double const a, double const b, PETScVector const& x)
     VecAXPBY(y.getRawVector(), a, b, x.getRawVector());
 }
 
+// Explicit specialization
+// Computes the Euclidean norm of x
+template<>
+double norm2(PETScVector const& x)
+{
+    return x.getNorm(MathLib::VecNormType::NORM2);
+}
+
 
 // Matrix
 
@@ -150,6 +179,14 @@ void axpby(EigenVector& y, double const a, double const b, EigenVector const& x)
     y.getRawVector() = a*x.getRawVector() + b*y.getRawVector();
 }
 
+// Explicit specialization
+// Euclidean norm
+template<>
+double norm2(EigenVector const& x)
+{
+    return x.getRawVector().norm();
+}
+
 
 // Matrix
 
diff --git a/MathLib/LinAlg/BLAS.h b/MathLib/LinAlg/BLAS.h
index b8b96dd8d7053a0b38461bb39361630803c99f4b..3c3623b2a8ae7ef793297763b765d49faeeffe86 100644
--- a/MathLib/LinAlg/BLAS.h
+++ b/MathLib/LinAlg/BLAS.h
@@ -66,6 +66,10 @@ void axpby(MatrixOrVector& y, double const a, double const b, MatrixOrVector con
     y = a*x + b*y;
 }
 
+//! Computes the Euclidean norm of \c x.
+template<typename MatrixOrVector>
+double norm2(MatrixOrVector const& x);
+
 
 // Matrix and Vector