Skip to content
Snippets Groups Projects
Commit 6fdf16cb authored by Andreas Dedner's avatar Andreas Dedner
Browse files

fix some minor issues with amg and complex

extended the amgtest to also use field_type = complex by switching the
definition of XREAL
parent dcf41959
No related branches found
No related tags found
No related merge requests found
......@@ -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());
}
......
......@@ -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;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment