From 95fa515a5a3a3362e4869112fc62f71b366bb9a7 Mon Sep 17 00:00:00 2001
From: Markus Blatt <markus@dr-blatt.de>
Date: Fri, 22 May 2015 20:47:02 +0200
Subject: [PATCH] Adds the possibility to compute a parallel inner product.

---
 opm/core/linalg/ParallelIstlInformation.hpp | 37 +++++++++++++++++++++
 tests/test_parallelistlinformation.cpp      | 10 ++++--
 2 files changed, 44 insertions(+), 3 deletions(-)

diff --git a/opm/core/linalg/ParallelIstlInformation.hpp b/opm/core/linalg/ParallelIstlInformation.hpp
index 090fd707..d7f9d383 100644
--- a/opm/core/linalg/ParallelIstlInformation.hpp
+++ b/opm/core/linalg/ParallelIstlInformation.hpp
@@ -403,6 +403,37 @@ private:
         BinaryOperator b_;
     };
 
+    /// \brief An operator for computing a parallel inner product.
+    template<class T>
+    struct InnerProductFunctor
+    {
+        /// \brief Apply the underlying binary operator according to the mask.
+        ///
+        /// The BinaryOperator will be called with t1, and mask*t2.
+        /// \param t1 first value
+        /// \param t2 second value (might be modified).
+        /// \param mask The mask (0 or 1).
+        template<class T1>
+        T operator()(const T& t1, const T& t2, const T1& mask)
+        {
+            T masked =  maskValue(t2, mask);
+            return t1 + masked * masked;
+        }
+        template<class T1>
+        T maskValue(const T& t, const T1& mask)
+        {
+            return t*mask;
+        }
+        std::plus<T> localOperator()
+        {
+            return std::plus<T>();
+        }
+        T getInitialValue()
+        {
+            return T();
+        }
+    };
+
     /// \brief An operator that converts the values where mask is 0 to the minimum value
     ///
     /// Could be used to compute a global maximum.
@@ -549,6 +580,12 @@ private:
             (std::pointer_to_binary_function<const T&,const T&,const T&>
              ((const T&(*)(const T&, const T&))std::min<T>));
     }
+    template<class T>
+    InnerProductFunctor<T>
+    makeInnerProductFunctor()
+    {
+        return InnerProductFunctor<T>();
+    }
     } // end namespace Reduction
 } // end namespace Opm
 
diff --git a/tests/test_parallelistlinformation.cpp b/tests/test_parallelistlinformation.cpp
index ecea92cd..7c54d72b 100644
--- a/tests/test_parallelistlinformation.cpp
+++ b/tests/test_parallelistlinformation.cpp
@@ -48,16 +48,20 @@ void runSumMaxMinTest(const int offset)
     assert(comm.indexSet()->size()==x.size());
     for(auto i=comm.indexSet()->begin(), iend=comm.indexSet()->end(); i!=iend; ++i)
         x[i->local()]=i->global()+offset;
-    auto containers = std::make_tuple(x, x, x);
+    auto containers = std::make_tuple(x, x, x, x);
     auto operators  = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<int>(),
                                       Opm::Reduction::makeGlobalMaxFunctor<int>(),
-                                      Opm::Reduction::makeGlobalMinFunctor<int>());
-    auto values     = std::make_tuple(0,0,100000);
+                                      Opm::Reduction::makeGlobalMinFunctor<int>(),
+                                      Opm::Reduction::makeInnerProductFunctor<int>());
+    auto values     = std::make_tuple(0,0,100000, 0);
     auto oldvalues  = values;
+    start = offset;
+    end   = start+N;
     comm.computeReduction(containers,operators,values);
     BOOST_CHECK(std::get<0>(values)==std::get<0>(oldvalues)+((N-1+2*offset)*N)/2);
     BOOST_CHECK(std::get<1>(values)==std::max(N+offset-1, std::get<1>(oldvalues)));
     BOOST_CHECK(std::get<2>(values)==std::min(offset, std::get<2>(oldvalues)));
+    BOOST_CHECK(std::get<3>(values)==((end-1)*end*(2*end-1)-(start-1)*start*(2*start-1))/6+std::get<3>(oldvalues));
 }
 
 BOOST_AUTO_TEST_CASE(tupleReductionTest)