From 2e5259107f39641af210e6c6dddbd4edd07c263f Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Tue, 22 Jan 2019 22:36:46 +0100
Subject: [PATCH] Implement AMG for scalar-valued matrices
---
dune/istl/paamg/aggregates.hh | 28 ++++++++++++++++++++++++++--
dune/istl/paamg/amg.hh | 10 +++++++++-
dune/istl/paamg/test/amgtest.cc | 7 +++++++
3 files changed, 42 insertions(+), 3 deletions(-)
diff --git a/dune/istl/paamg/aggregates.hh b/dune/istl/paamg/aggregates.hh
index cf0bdd9cd..a2a17659b 100644
--- a/dune/istl/paamg/aggregates.hh
+++ b/dune/istl/paamg/aggregates.hh
@@ -9,6 +9,7 @@
#include "properties.hh"
#include "combinedfunctor.hh"
+#include <dune/common/hybridutilities.hh>
#include <dune/common/timer.hh>
#include <dune/common/stdstreams.hh>
#include <dune/common/poolallocator.hh>
@@ -382,7 +383,8 @@ namespace Dune
* @param m The matrix to compute the norm of
*/
template<class M>
- typename FieldTraits<typename M::field_type>::real_type operator()(const M& m) const
+ typename FieldTraits<typename M::field_type>::real_type operator()(const M& m,
+ typename std::enable_if_t<!Dune::IsNumber<M>::value>* sfinae = nullptr) const
{
typedef typename M::field_type field_type;
typedef typename FieldTraits<field_type>::real_type real_type;
@@ -392,6 +394,21 @@ namespace Dune
// possible implementation for complex types: return signed_abs(m[N][N]);
}
+ /**
+ * @brief Compute the norm of a scalar
+ * @param m The scalar to compute the norm of
+ */
+ template<class M>
+ auto operator()(const M& m,
+ typename std::enable_if_t<Dune::IsNumber<M>::value>* sfinae = nullptr) const
+ {
+ typedef typename FieldTraits<M>::real_type real_type;
+ static_assert( std::is_convertible<M, real_type >::value,
+ "use of diagonal norm in AMG not implemented for complex field_type");
+ return m;
+ // possible implementation for complex types: return signed_abs(m[N][N]);
+ }
+
private:
//! return sign * abs_value; for real numbers this is just v
@@ -1813,7 +1830,14 @@ namespace Dune
for(ColIterator col = row.begin(); col != end; ++col)
if(col.index()!=*vertex) {
criterion.examine(col);
- absoffdiag = max(absoffdiag, col->frobenius_norm());
+ Hybrid::ifElse(IsNumber<typename Matrix::block_type>(),
+ [&](auto id) {
+ using std::abs;
+ absoffdiag = max(absoffdiag, abs(*id(col)));
+ },
+ [&](auto id) {
+ absoffdiag = max(absoffdiag, id(col)->frobenius_norm());
+ });
}
if(absoffdiag==0)
diff --git a/dune/istl/paamg/amg.hh b/dune/istl/paamg/amg.hh
index 44ce9c3e5..7551b8dda 100644
--- a/dune/istl/paamg/amg.hh
+++ b/dune/istl/paamg/amg.hh
@@ -584,7 +584,15 @@ namespace Dune
}
}
if(isDirichlet && hasDiagonal)
- diagonal.solve(x[row.index()], b[row.index()]);
+ {
+ Hybrid::ifElse(IsNumber<Block>(),
+ [&](auto id) {
+ x[row.index()] = id(diagonal) / b[row.index()];
+ },
+ [&](auto id) {
+ id(diagonal).solve(x[row.index()], b[row.index()]);
+ });
+ }
}
if(smoothers_->levels()>0)
diff --git a/dune/istl/paamg/test/amgtest.cc b/dune/istl/paamg/test/amgtest.cc
index 826496d29..d3be925cc 100644
--- a/dune/istl/paamg/test/amgtest.cc
+++ b/dune/istl/paamg/test/amgtest.cc
@@ -179,6 +179,13 @@ try
if(argc>3)
ml = atoi(argv[3]);
+ {
+ using Matrix = Dune::BCRSMatrix<XREAL>;
+ using Vector = Dune::BlockVector<XREAL>;
+
+ testAMG<Matrix,Vector>(N, coarsenTarget, ml);
+ }
+
{
using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<XREAL,1,1> >;
using Vector = Dune::BlockVector<Dune::FieldVector<XREAL,1> >;
--
GitLab