From ff386db24f5f4db81a29935377b639a36a022fdc Mon Sep 17 00:00:00 2001
From: Christoph Lehmann <christoph.lehmann@ufz.de>
Date: Sun, 24 Apr 2016 00:31:23 +0200
Subject: [PATCH] [MaL] added componentwiseDivide to BLAS
---
MathLib/LinAlg/BLAS.cpp | 28 ++++++++++++++++++++++++++++
MathLib/LinAlg/BLAS.h | 5 +++++
2 files changed, 33 insertions(+)
diff --git a/MathLib/LinAlg/BLAS.cpp b/MathLib/LinAlg/BLAS.cpp
index bffbbe9a722..87550f22621 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 355af434f79..a042b045446 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);
--
GitLab