adding helper class to collect L,R,S components of self-energy

FlexibleSUSY · Mar 28, 2017 · 2363afa · 2363afa
1 parent e2de61a
commit 2363afa
Show file tree

Hide file tree

Showing 2 changed files with 95 additions and 0 deletions.
diff --git a/src/eigen_utils.hpp b/src/eigen_utils.hpp
@@ -28,6 +28,46 @@
 
 namespace flexiblesusy {
 
+template <int N, typename Scalar = double>
+struct LRS_tensor {
+public:
+   using Matrix_t = Eigen::Matrix<Scalar,N,N>;
+   using Tuple_t = std::tuple<Matrix_t,Matrix_t,Matrix_t>;
+
+   LRS_tensor()
+      : lrs(std::make_tuple(Matrix_t::Zero(),Matrix_t::Zero(),Matrix_t::Zero())) {}
+   explicit LRS_tensor(const Tuple_t& t)
+      : lrs(t) {}
+
+   Matrix_t L() const { return std::get<0>(lrs); }
+   Matrix_t R() const { return std::get<1>(lrs); }
+   Matrix_t S() const { return std::get<2>(lrs); }
+   Tuple_t tuple() const { return lrs; }
+   static LRS_tensor PL() { return LRS_tensor(std::make_tuple(Matrix_t::Ones(),Matrix_t::Zero(),Matrix_t::Zero())); }
+   static LRS_tensor PR() { return LRS_tensor(std::make_tuple(Matrix_t::Zero(),Matrix_t::Ones(),Matrix_t::Zero())); }
+   static LRS_tensor PS() { return LRS_tensor(std::make_tuple(Matrix_t::Zero(),Matrix_t::Zero(),Matrix_t::Ones())); }
+
+   LRS_tensor operator+(const LRS_tensor& rhs)
+   {
+      return LRS_tensor(std::make_tuple((L()+rhs.L()).eval(), (R()+rhs.R()).eval(), (S()+rhs.S()).eval()));
+   }
+
+   template <typename T>
+   LRS_tensor operator*(T rhs)
+   {
+      return LRS_tensor(std::make_tuple((rhs*L()).eval(), (rhs*R()).eval(), (rhs*S()).eval()));
+   }
+
+   template <typename T>
+   friend LRS_tensor operator*(T lhs, const LRS_tensor& rhs)
+   {
+      return LRS_tensor(std::make_tuple((lhs*rhs.L()).eval(), (lhs*rhs.R()).eval(), (lhs*rhs.S()).eval()));
+   }
+
+private:
+   Tuple_t lrs; ///< L, R, S components
+};
+
 template <typename Derived>
 int closest_index(double mass, const Eigen::ArrayBase<Derived>& v)
 {

diff --git a/test/test_eigen_utils.cpp b/test/test_eigen_utils.cpp
@@ -22,6 +22,7 @@
 #include <boost/test/unit_test.hpp>
 #include <typeinfo>
 #include "eigen_utils.hpp"
+#include "wrappers.hpp"
 #include <complex>
 
 using namespace flexiblesusy;
@@ -146,3 +147,57 @@ BOOST_AUTO_TEST_CASE(test_Eval)
    BOOST_CHECK_EQUAL(typeid(Eval(i)).hash_code(), typeid(i).hash_code());
 
 }
+
+BOOST_AUTO_TEST_CASE(test_LRS_tensor_sum)
+{
+   using Tensor_t = LRS_tensor<2>;
+   auto t = Tensor_t::PL() + Tensor_t::PR() + Tensor_t::PS();
+   BOOST_CHECK_EQUAL(t.L(), Tensor_t::Matrix_t::Ones());
+   BOOST_CHECK_EQUAL(t.R(), Tensor_t::Matrix_t::Ones());
+   BOOST_CHECK_EQUAL(t.S(), Tensor_t::Matrix_t::Ones());
+}
+
+BOOST_AUTO_TEST_CASE(test_LRS_tensor_product_int)
+{
+   using Tensor_t = LRS_tensor<2>;
+   const int i = 2;
+   auto t = i*Tensor_t::PL() + Tensor_t::PR()*2*i + Tensor_t::PS();
+   BOOST_CHECK_EQUAL(t.L(), 2*Tensor_t::Matrix_t::Ones());
+   BOOST_CHECK_EQUAL(t.R(), 4*Tensor_t::Matrix_t::Ones());
+   BOOST_CHECK_EQUAL(t.S(), 1*Tensor_t::Matrix_t::Ones());
+}
+
+BOOST_AUTO_TEST_CASE(test_LRS_tensor_product_double)
+{
+   using Tensor_t = LRS_tensor<2>;
+   const double d = 2.;
+   auto t = d*Tensor_t::PL() + Tensor_t::PR()*2*d + Tensor_t::PS();
+   BOOST_CHECK_EQUAL(t.L(), 2*Tensor_t::Matrix_t::Ones());
+   BOOST_CHECK_EQUAL(t.R(), 4*Tensor_t::Matrix_t::Ones());
+   BOOST_CHECK_EQUAL(t.S(), 1*Tensor_t::Matrix_t::Ones());
+}
+
+BOOST_AUTO_TEST_CASE(test_LRS_tensor_product_complex)
+{
+   using Tensor_t = LRS_tensor<2,std::complex<double> >;
+   const auto c = std::complex<double>(2.,2.);
+   auto t = c*Tensor_t::PL() + Tensor_t::PR()*2*c + Tensor_t::PS();
+   BOOST_CHECK_EQUAL(t.L(), std::complex<double>(2.,2.)*  Tensor_t::Matrix_t::Ones());
+   BOOST_CHECK_EQUAL(t.R(), std::complex<double>(2.,2.)*2*Tensor_t::Matrix_t::Ones());
+   BOOST_CHECK_EQUAL(t.S(), Tensor_t::Matrix_t::Ones());
+}
+
+BOOST_AUTO_TEST_CASE(test_LRS_tensor_product_Matrix)
+{
+   using Tensor_t = LRS_tensor<2,std::complex<double> >;
+   using Matrix_t = Tensor_t::Matrix_t;
+   const auto c = std::complex<double>(2.,2.);
+   Matrix_t m;
+   m << std::complex<double>(2.,0.), std::complex<double>(0.,0.),
+        std::complex<double>(0.,0.), std::complex<double>(2.,0.);
+
+   auto t = c*m*Tensor_t::PL() + Tensor_t::PR()*2*m + Tensor_t::PS()*m;
+   BOOST_CHECK_EQUAL(t.L(), std::complex<double>(4.,4.)*  Tensor_t::Matrix_t::Ones());
+   BOOST_CHECK_EQUAL(t.R(), std::complex<double>(2.,0.)*2*Tensor_t::Matrix_t::Ones());
+   BOOST_CHECK_EQUAL(t.S(), 2.*Tensor_t::Matrix_t::Ones());
+}