From 53544b2cfa102a3027003631976079d65be74e68 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Wed, 11 Sep 2019 11:11:30 +0200
Subject: [PATCH] Implement negation for dense vectors and matrices

In other words, with this patch you can now write

FieldVector<double, 3> a = {1,2,3};
FieldVector<double, 3> b = -a;

and similarly for matrices.
---
 CHANGELOG.md                    |  5 +++++
 dune/common/densematrix.hh      | 13 +++++++++++++
 dune/common/densevector.hh      | 12 ++++++++++++
 dune/common/test/fmatrixtest.cc |  9 +++++++++
 dune/common/test/fvectortest.cc |  3 +++
 5 files changed, 42 insertions(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index b16fba5f6..bed3db0df 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -25,6 +25,11 @@
   implementation in `dune-common` is automatically disabled, and the official
   implementation from the standard library is used instead.
 
+- By popular demand, dense vectors and matrices like `FieldVector` and `FieldMatrix`
+  now have additional operators.  In particular, there are
+  - Vector = - Vector
+  - Matrix = - Matrix
+
 - There is now (finally!) a method `power` in the file `math.hh` that computes
   powers with an integer exponent, and is usable in compile-time expressions.
   The use of the old power methods in `power.hh` is henceforth discouraged.
diff --git a/dune/common/densematrix.hh b/dune/common/densematrix.hh
index b03bbb0b5..5e511e922 100644
--- a/dune/common/densematrix.hh
+++ b/dune/common/densematrix.hh
@@ -321,6 +321,19 @@ namespace Dune
       return asImp();
     }
 
+    //! Matrix negation
+    derived_type operator- () const
+    {
+      MAT result;
+      typedef typename decltype(result)::size_type size_type;
+
+      for (size_type i = 0; i < rows(); ++i)
+        for (size_type j = 0; j < cols(); ++j)
+          result[i][j] = - asImp()[i][j];
+
+      return result;
+    }
+
     //! vector space subtraction
     template <class Other>
     derived_type &operator-= (const DenseMatrix<Other>& x)
diff --git a/dune/common/densevector.hh b/dune/common/densevector.hh
index 5c2bfacf7..e3951a312 100644
--- a/dune/common/densevector.hh
+++ b/dune/common/densevector.hh
@@ -448,6 +448,18 @@ namespace Dune {
       return (z-=b);
     }
 
+    //! Vector negation
+    derived_type operator- () const
+    {
+      V result;
+      typedef typename decltype(result)::size_type size_type;
+
+      for (size_type i = 0; i < size(); ++i)
+        result[i] = -asImp()[i];
+
+      return result;
+    }
+
     //! \brief vector space add scalar to all comps
     /**
        we use enable_if to avoid an ambiguity, if the
diff --git a/dune/common/test/fmatrixtest.cc b/dune/common/test/fmatrixtest.cc
index 9c6adcf52..b92f16690 100644
--- a/dune/common/test/fmatrixtest.cc
+++ b/dune/common/test/fmatrixtest.cc
@@ -416,6 +416,15 @@ void test_matrix()
       if (tmp.infinity_norm() > 1e-12)
         DUNE_THROW(FMatrixError, "Return value of axpy() incorrect!");
     }
+    // -Matrix
+    {
+      FM neg = -A;
+      FM ref = typename FM::field_type(-1) * A;
+
+      if ((neg-ref).infinity_norm() > 1e-12)
+        DUNE_THROW(FMatrixError, "Return value of operator-(matrix) incorrect!");
+    }
+
   }
   {
     using std::abs;
diff --git a/dune/common/test/fvectortest.cc b/dune/common/test/fvectortest.cc
index bc411c081..455816d1f 100644
--- a/dune/common/test/fvectortest.cc
+++ b/dune/common/test/fvectortest.cc
@@ -102,6 +102,9 @@ struct FieldVectorMainTestCommons
     w *= a;
     w /= a;
 
+    // Negation
+    -v;
+
     // test scalar product, axpy
     a = v * w;
     a = v.dot(w);
-- 
GitLab