diff --git a/dune/istl/paamg/aggregates.hh b/dune/istl/paamg/aggregates.hh index 67a7886c716a1250ff46e15e001b0993e2840eca..77449dc81c5982e12ed04756407a94ee35ec5640 100644 --- a/dune/istl/paamg/aggregates.hh +++ b/dune/istl/paamg/aggregates.hh @@ -182,8 +182,8 @@ namespace Dune int row_; /** @brief The norm of the current diagonal. */ real_type diagonal_; - std::vector<typename Matrix::field_type> vals_; - typename std::vector<typename Matrix::field_type>::iterator valIter_; + std::vector<real_type> vals_; + typename std::vector<real_type>::iterator valIter_; }; @@ -201,7 +201,7 @@ namespace Dune assert(vals_.size()==row.size()); valIter_=vals_.begin(); - maxValue_ = std::min(- std::numeric_limits<typename Matrix::field_type>::max(), std::numeric_limits<typename Matrix::field_type>::min()); + maxValue_ = std::min(- std::numeric_limits<real_type>::max(), std::numeric_limits<real_type>::min()); diagonal_=norm_(row[index]); row_ = index; } @@ -381,21 +381,26 @@ namespace Dune template<class M> typename FieldTraits<typename M::field_type>::real_type operator()(const M& m) const { - return signed_abs(m[N][N]); + typedef typename M::field_type field_type; + typedef typename FieldTraits<field_type>::real_type real_type; + static_assert( std::is_convertible<field_type, real_type >::value, + "use of diagonal norm in AMG not implemented for complex field_type"); + return m[N][N]; + // possible implementation for complex types: return signed_abs(m[N][N]); } private: //! return sign * abs_value; for real numbers this is just v template<typename T> - T signed_abs(const T & v) + static T signed_abs(const T & v) { return v; } //! return sign * abs_value; for complex numbers this is csgn(v) * abs(v) template<typename T> - T signed_abs(const std::complex<T> & v) + static T signed_abs(const std::complex<T> & v) { // return sign * abs_value // in case of complex numbers this extends to using the csgn function to determine the sign @@ -404,14 +409,14 @@ namespace Dune //! sign function for complex numbers; for real numbers we assume imag(v) = 0 template<typename T> - T csgn(const T & v) + static T csgn(const T & v) { return (T(0) < v) - (v < T(0)); } //! sign function for complex numbers template<typename T> - T csgn(std::complex<T> a) + static T csgn(std::complex<T> a) { return csgn(a.real())+(a.real() == 0.0)*csgn(a.imag()); } diff --git a/dune/istl/paamg/test/amgtest.cc b/dune/istl/paamg/test/amgtest.cc index efc0817210f56f11c12e4bc0109371cc5d912306..9d09fc45bb1fc370278e529b23841844ec637c49 100644 --- a/dune/istl/paamg/test/amgtest.cc +++ b/dune/istl/paamg/test/amgtest.cc @@ -24,8 +24,10 @@ #include <dune/istl/solvers.hh> #include <cstdlib> #include <ctime> +#include <complex> typedef double XREAL; +// typedef std::complex<double> XREAL; /* typedef HPA::xreal XREAL; @@ -99,8 +101,10 @@ void testAMG(int N, int coarsenTarget, int ml) watch.reset(); Operator fop(mat); - typedef Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<BCRSMat,Dune::Amg::FirstDiagonal> > - Criterion; + typedef typename std::conditional< std::is_convertible<XREAL, typename Dune::FieldTraits<XREAL>::real_type>::value, + Dune::Amg::FirstDiagonal, Dune::Amg::RowSum >::type Norm; + typedef Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<BCRSMat,Norm> > + Criterion; typedef Dune::SeqSSOR<BCRSMat,Vector,Vector> Smoother; //typedef Dune::SeqSOR<BCRSMat,Vector,Vector> Smoother; //typedef Dune::SeqJac<BCRSMat,Vector,Vector> Smoother;