Skip to content
Snippets Groups Projects
Commit 5af571d9 authored by Carsten Gräser's avatar Carsten Gräser
Browse files

Merge branch 'origin/feature/FS1622' into releases/2.4 

This brings releases/2.4 on par with "master" at 287b0ae9.
We can't do a fast forward because another merge commit was
cherry-picked in between.
parents 69c6ba89 287b0ae9
No related branches found
No related tags found
1 merge request!112Add data() accessors to BlockVector
Hide whitespace changes
Inline Side-by-side
dune/istl/ilu.hh
+ 3
− 1
View file @ 5af571d9
......@@ -188,7 +188,9 @@ namespace Dune {
coliterator kj = ILU[(*ik).first].find((*ik).first); // diagonal in k
for (++kj; kj!=endk; ++kj) // row k eliminates in row i
{
int generation = (int) firstmatrixelement(*kj);
// we misuse the storage to store an int. If the field_type is std::complex, we have to access the real part
// starting from C++11, we can use std::real to always return a real value, even if it is double/float
int generation = (int) std::real( firstmatrixelement(*kj) );
if (generation<n)
{
mapiterator ij = rowpattern.find(kj.index());
......
dune/istl/paamg/aggregates.hh
+ 65
− 20
View file @ 5af571d9
......@@ -19,6 +19,7 @@
#include <utility>
#include <set>
#include <algorithm>
#include <complex>
#include <limits>
#include <ostream>
......@@ -172,15 +173,17 @@ namespace Dune
/** @brief The matrix we work on. */
const Matrix* matrix_;
/** @brief The current max value.*/
typename Matrix::field_type maxValue_;
typedef typename Matrix::field_type field_type;
typedef typename FieldTraits<field_type>::real_type real_type;
real_type maxValue_;
/** @brief The functor for calculating the norm. */
Norm norm_;
/** @brief index of the currently evaluated row. */
int row_;
/** @brief The norm of the current diagonal. */
typename Matrix::field_type diagonal_;
std::vector<typename Matrix::field_type> vals_;
typename std::vector<typename Matrix::field_type>::iterator valIter_;
real_type diagonal_;
std::vector<real_type> vals_;
typename std::vector<real_type>::iterator valIter_;
};
......@@ -198,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;
}
......@@ -207,7 +210,7 @@ namespace Dune
inline void SymmetricMatrixDependency<M,N>::examine(const ColIter& col)
{
// skip positive offdiagonals if norm preserves sign of them.
typename Matrix::field_type eij = norm_(*col);
real_type eij = norm_(*col);
if(!N::is_sign_preserving || eij<0) // || eji<0)
{
*valIter_ = eij/diagonal_*eij/norm_(matrix_->operator[](col.index())[col.index()]);
......@@ -287,13 +290,15 @@ namespace Dune
/** @brief The matrix we work on. */
const Matrix* matrix_;
/** @brief The current max value.*/
typename Matrix::field_type maxValue_;
typedef typename Matrix::field_type field_type;
typedef typename FieldTraits<field_type>::real_type real_type;
real_type maxValue_;
/** @brief The functor for calculating the norm. */
Norm norm_;
/** @brief index of the currently evaluated row. */
int row_;
/** @brief The norm of the current diagonal. */
typename Matrix::field_type diagonal_;
real_type diagonal_;
};
/**
......@@ -346,13 +351,15 @@ namespace Dune
/** @brief The matrix we work on. */
const Matrix* matrix_;
/** @brief The current max value.*/
typename Matrix::field_type maxValue_;
typedef typename Matrix::field_type field_type;
typedef typename FieldTraits<field_type>::real_type real_type;
real_type maxValue_;
/** @brief The functor for calculating the norm. */
Norm norm_;
/** @brief index of the currently evaluated row. */
int row_;
/** @brief The norm of the current diagonal. */
typename Matrix::field_type diagonal_;
real_type diagonal_;
};
/**
......@@ -372,10 +379,48 @@ namespace Dune
* @param m The matrix ro compute the norm of.
*/
template<class M>
typename M::field_type operator()(const M& m) const
typename FieldTraits<typename M::field_type>::real_type operator()(const M& m) const
{
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>
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>
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
return csgn(v) * std::abs(v);
}
//! sign function for complex numbers; for real numbers we assume imag(v) = 0
template<typename T>
static T csgn(const T & v)
{
return (T(0) < v) - (v < T(0));
}
//! sign function for complex numbers
template<typename T>
static T csgn(std::complex<T> a)
{
return csgn(a.real())+(a.real() == 0.0)*csgn(a.imag());
}
};
/**
......@@ -401,7 +446,7 @@ namespace Dune
* @param m The matrix row to compute the norm of.
*/
template<class M>
typename M::field_type operator()(const M& m) const
typename FieldTraits<typename M::field_type>::real_type operator()(const M& m) const
{
return m.infinity_norm();
}
......@@ -418,7 +463,7 @@ namespace Dune
* @param m The matrix row to compute the norm of.
*/
template<class M>
typename M::field_type operator()(const M& m) const
typename FieldTraits<typename M::field_type>::real_type operator()(const M& m) const
{
return m.frobenius_norm();
}
......@@ -434,7 +479,7 @@ namespace Dune
* @param m The matrix row to compute the norm of.
*/
template<class M>
typename M::field_type operator()(const M& m) const
typename FieldTraits<typename M::field_type>::real_type operator()(const M& m) const
{
return 1;
}
......@@ -1362,8 +1407,8 @@ namespace Dune
template<class M, class N>
inline void SymmetricDependency<M,N>::examine(const ColIter& col)
{
typename Matrix::field_type eij = norm_(*col);
typename Matrix::field_type eji = norm_(matrix_->operator[](col.index())[row_]);
real_type eij = norm_(*col);
real_type eji = norm_(matrix_->operator[](col.index())[row_]);
// skip positive offdiagonals if norm preserves sign of them.
if(!N::is_sign_preserving || eij<0 || eji<0)
......@@ -1376,8 +1421,8 @@ namespace Dune
template<class G>
inline void SymmetricDependency<M,N>::examine(G& graph, const typename G::EdgeIterator& edge, const ColIter& col)
{
typename Matrix::field_type eij = norm_(*col);
typename Matrix::field_type eji = norm_(matrix_->operator[](col.index())[row_]);
real_type eij = norm_(*col);
real_type eji = norm_(matrix_->operator[](col.index())[row_]);
// skip positve offdiagonals if norm preserves sign of them.
if(!N::is_sign_preserving || (eij<0 || eji<0))
......@@ -1409,7 +1454,7 @@ namespace Dune
inline void Dependency<M,N>::initRow(const Row& row, int index)
{
DUNE_UNUSED_PARAMETER(row);
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());
row_ = index;
diagonal_ = norm_(matrix_->operator[](row_)[row_]);
}
......@@ -1758,7 +1803,7 @@ namespace Dune
// the calculator should know whether the vertex is isolated.
typedef typename Matrix::ConstColIterator ColIterator;
ColIterator end = row.end();
typename Matrix::field_type absoffdiag=0;
typename FieldTraits<typename Matrix::field_type>::real_type absoffdiag=0.;
if(firstlevel) {
for(ColIterator col = row.begin(); col != end; ++col)
......
dune/istl/paamg/test/amgtest.cc
+ 6
− 2
View file @ 5af571d9
......@@ -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% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment