Skip to content
Snippets Groups Projects
Commit c43f8c76 authored by Christian Engwer's avatar Christian Engwer
Browse files

[istl][multitype] fix solvers for MultiType* 

the old code didn't work (at least not with clang). The reason is that constructs like
template<typename... Arg1, typename... Arg2> don't work, as the  compiler doesn't know
where to split the list of arguments. The whole algo_itsteps got reworked to allow for an
easier specialization for MultiTypeBlockMatrix. Now the struct is also parametrized by the
matrix type. By this we can easily add a specialization.
parent d0069aa6
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/istl/gsetc.hh
+ 74
− 71
View file @ c43f8c76
......@@ -358,9 +358,6 @@ namespace Dune {
}
//============================================================
// generic steps of iteration methods
// Jacobi, Gauss-Seidel, SOR, SSOR
......@@ -368,23 +365,11 @@ namespace Dune {
// we can recurse over a fixed number of levels
//============================================================
// template meta program for iterative solver steps
template<int I>
template<int I, typename M>
struct algmeta_itsteps {
template<typename... MultiTypeMatrixArgs,
typename... MultiTypeVectorArgs,
class K>
static void dbgs (const MultiTypeBlockMatrix<MultiTypeMatrixArgs...>& A,
MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
const K& w)
{
Dune::MultiTypeBlockMatrix_Solver<I,0,A.N()>::dbgs(A, x, b, w);
}
template<class M, class X, class Y, class K>
template<class X, class Y, class K>
static void dbgs (const M& A, X& x, const Y& b, const K& w)
{
typedef typename M::ConstRowIterator rowiterator;
......@@ -405,25 +390,14 @@ namespace Dune {
coliterator diag=j++; // *diag = a_ii and increment coliterator j from a_ii to a_i+1,i to skip diagonal
for (; j != endj; ++j) // iterate over a_ij with j > i
(*j).mmv(x[j.index()],rhs); // rhs -= sum_{j>i} a_ij * xold_j
algmeta_itsteps<I-1>::dbgs(*diag,x[i.index()],rhs,w); // if I==1: xnew_i = rhs/a_ii
algmeta_itsteps<I-1,typename M::block_type>::dbgs(*diag,x[i.index()],rhs,w); // if I==1: xnew_i = rhs/a_ii
}
// next two lines: xnew_i = w / a_ii * (b_i - sum_{j<i} a_ij * xnew_j - sum_{j>=i} a_ij * xold_j) + (1-w)*xold;
x *= w;
x.axpy(K(1)-w,xold);
}
template<typename... MultiTypeMatrixArgs,
typename... MultiTypeVectorArgs,
class K>
static void bsorf (const MultiTypeBlockMatrix<MultiTypeMatrixArgs...>& A,
MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
const K& w)
{
Dune::MultiTypeBlockMatrix_Solver<I,0,A.N()>::bsorf(A, x, b, w);
}
template<class M, class X, class Y, class K>
template<class X, class Y, class K>
static void bsorf (const M& A, X& x, const Y& b, const K& w)
{
typedef typename M::ConstRowIterator rowiterator;
......@@ -448,23 +422,12 @@ namespace Dune {
coliterator diag=j; // *diag = a_ii
for (; j!=endj; ++j)
(*j).mmv(x[j.index()],rhs); // rhs -= sum_{j<i} a_ij * xnew_j
algmeta_itsteps<I-1>::bsorf(*diag,v,rhs,w); // if blocksize I==1: v = rhs/a_ii
algmeta_itsteps<I-1,typename M::block_type>::bsorf(*diag,v,rhs,w); // if blocksize I==1: v = rhs/a_ii
x[i.index()].axpy(w,v); // x_i = w / a_ii * (b_i - sum_{j<i} a_ij * xnew_j - sum_{j>=i} a_ij * xold_j)
}
}
template<typename... MultiTypeMatrixArgs,
typename... MultiTypeVectorArgs,
class K>
static void bsorb (const MultiTypeBlockMatrix<MultiTypeMatrixArgs...>& A,
MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
const K& w)
{
Dune::MultiTypeBlockMatrix_Solver<I,A.N()-1,A.N()>::bsorb(A, x, b, w);
}
template<class M, class X, class Y, class K>
template<class X, class Y, class K>
static void bsorb (const M& A, X& x, const Y& b, const K& w)
{
typedef typename M::ConstRowIterator rowiterator;
......@@ -489,23 +452,12 @@ namespace Dune {
coliterator diag=j;
for (; j!=endj; ++j)
(*j).mmv(x[j.index()],rhs);
algmeta_itsteps<I-1>::bsorb(*diag,v,rhs,w);
algmeta_itsteps<I-1,typename M::block_type>::bsorb(*diag,v,rhs,w);
x[i.index()].axpy(w,v);
}
}
template<typename... MultiTypeMatrixArgs,
typename... MultiTypeVectorArgs,
class K>
static void dbjac (const MultiTypeBlockMatrix<MultiTypeMatrixArgs...>& A,
MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
const K& w)
{
Dune::MultiTypeBlockMatrix_Solver<I,0,A.N()>::dbjac(A, x, b, w);
}
template<class M, class X, class Y, class K>
template<class X, class Y, class K>
static void dbjac (const M& A, X& x, const Y& b, const K& w)
{
typedef typename M::ConstRowIterator rowiterator;
......@@ -526,36 +478,87 @@ namespace Dune {
coliterator diag=j;
for (; j!=endj; ++j)
(*j).mmv(x[j.index()],rhs);
algmeta_itsteps<I-1>::dbjac(*diag,v[i.index()],rhs,w);
algmeta_itsteps<I-1,typename M::block_type>::dbjac(*diag,v[i.index()],rhs,w);
}
x.axpy(w,v);
}
};
// end of recursion
template<>
struct algmeta_itsteps<0> {
template<class M, class X, class Y, class K>
template<typename M>
struct algmeta_itsteps<0,M> {
template<class X, class Y, class K>
static void dbgs (const M& A, X& x, const Y& b, const K& /*w*/)
{
A.solve(x,b);
}
template<class M, class X, class Y, class K>
template<class X, class Y, class K>
static void bsorf (const M& A, X& x, const Y& b, const K& /*w*/)
{
A.solve(x,b);
}
template<class M, class X, class Y, class K>
template<class X, class Y, class K>
static void bsorb (const M& A, X& x, const Y& b, const K& /*w*/)
{
A.solve(x,b);
}
template<class M, class X, class Y, class K>
template<class X, class Y, class K>
static void dbjac (const M& A, X& x, const Y& b, const K& /*w*/)
{
A.solve(x,b);
}
};
template<int I, typename T1, typename... MultiTypeMatrixArgs>
struct algmeta_itsteps<I,MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>> {
template<
typename... MultiTypeVectorArgs,
class K>
static void dbgs (const MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>& A,
MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
const K& w)
{
static const int N = MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>::N();
Dune::MultiTypeBlockMatrix_Solver<I,0,N>::dbgs(A, x, b, w);
}
template<
typename... MultiTypeVectorArgs,
class K>
static void bsorf (const MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>& A,
MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
const K& w)
{
static const int N = MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>::N();
Dune::MultiTypeBlockMatrix_Solver<I,0,N>::bsorf(A, x, b, w);
}
template<
typename... MultiTypeVectorArgs,
class K>
static void bsorb (const MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>& A,
MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
const K& w)
{
static const int N = MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>::N();
Dune::MultiTypeBlockMatrix_Solver<I,N-1,N>::bsorb(A, x, b, w);
}
template<
typename... MultiTypeVectorArgs,
class K
>
static void dbjac (const MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>& A,
MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
const K& w)
{
static const int N = MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>::N();
Dune::MultiTypeBlockMatrix_Solver<I,0,N>::dbjac(A, x, b, w);
}
};
// user calls
......@@ -563,49 +566,49 @@ namespace Dune {
template<class M, class X, class Y, class K>
void dbgs (const M& A, X& x, const Y& b, const K& w)
{
algmeta_itsteps<1>::dbgs(A,x,b,w);
algmeta_itsteps<1,M>::dbgs(A,x,b,w);
}
//! GS step
template<class M, class X, class Y, class K, int l>
void dbgs (const M& A, X& x, const Y& b, const K& w, BL<l> /*bl*/)
{
algmeta_itsteps<l>::dbgs(A,x,b,w);
algmeta_itsteps<l,M>::dbgs(A,x,b,w);
}
//! SOR step
template<class M, class X, class Y, class K>
void bsorf (const M& A, X& x, const Y& b, const K& w)
{
algmeta_itsteps<1>::bsorf(A,x,b,w);
algmeta_itsteps<1,M>::bsorf(A,x,b,w);
}
//! SOR step
template<class M, class X, class Y, class K, int l>
void bsorf (const M& A, X& x, const Y& b, const K& w, BL<l> /*bl*/)
{
algmeta_itsteps<l>::bsorf(A,x,b,w);
algmeta_itsteps<l,M>::bsorf(A,x,b,w);
}
//! SSOR step
template<class M, class X, class Y, class K>
void bsorb (const M& A, X& x, const Y& b, const K& w)
{
algmeta_itsteps<1>::bsorb(A,x,b,w);
algmeta_itsteps<1,M>::bsorb(A,x,b,w);
}
//! Backward SOR step
template<class M, class X, class Y, class K, int l>
void bsorb (const M& A, X& x, const Y& b, const K& w, BL<l> /*bl*/)
{
algmeta_itsteps<l>::bsorb(A,x,b,w);
algmeta_itsteps<l,typename std::remove_cv<M>::type>::bsorb(A,x,b,w);
}
//! Jacobi step
template<class M, class X, class Y, class K>
void dbjac (const M& A, X& x, const Y& b, const K& w)
{
algmeta_itsteps<1>::dbjac(A,x,b,w);
algmeta_itsteps<1,M>::dbjac(A,x,b,w);
}
//! Jacobi step
template<class M, class X, class Y, class K, int l>
void dbjac (const M& A, X& x, const Y& b, const K& w, BL<l> /*bl*/)
{
algmeta_itsteps<l>::dbjac(A,x,b,w);
algmeta_itsteps<l,M>::dbjac(A,x,b,w);
}
......
dune/istl/multitypeblockmatrix.hh
+ 28
− 4
View file @ c43f8c76
......@@ -401,7 +401,13 @@ namespace Dune {
MultiTypeBlockMatrix_Solver_Col<I,crow,0, TMatrix::M()>::calc_rhs(A,x,v,rhs,w); // calculate right side of equation
//solve on blocklevel I-1
algmeta_itsteps<I-1>::dbgs(std::get<crow>( std::get<crow>(A)), std::get<crow>(x),rhs,w);
using M =
typename std::remove_cv<
typename std::remove_reference<
decltype(std::get<crow>( std::get<crow>(A)))
>::type
>::type;
algmeta_itsteps<I-1,M>::dbgs(std::get<crow>( std::get<crow>(A)), std::get<crow>(x),rhs,w);
MultiTypeBlockMatrix_Solver<I,crow+1,remain_row-1>::dbgs(A,x,v,b,w); //next row
}
......@@ -423,7 +429,13 @@ namespace Dune {
MultiTypeBlockMatrix_Solver_Col<I,crow,0,TMatrix::M()>::calc_rhs(A,x,v,rhs,w); // calculate right side of equation
//solve on blocklevel I-1
algmeta_itsteps<I-1>::bsorf(std::get<crow>( std::get<crow>(A)), std::get<crow>(v),rhs,w);
using M =
typename std::remove_cv<
typename std::remove_reference<
decltype(std::get<crow>( std::get<crow>(A)))
>::type
>::type;
algmeta_itsteps<I-1,M>::bsorf(std::get<crow>( std::get<crow>(A)), std::get<crow>(v),rhs,w);
std::get<crow>(x).axpy(w,std::get<crow>(v));
MultiTypeBlockMatrix_Solver<I,crow+1,remain_row-1>::bsorf(A,x,v,b,w); //next row
}
......@@ -444,7 +456,13 @@ namespace Dune {
MultiTypeBlockMatrix_Solver_Col<I,crow,0, TMatrix::M()>::calc_rhs(A,x,v,rhs,w); // calculate right side of equation
//solve on blocklevel I-1
algmeta_itsteps<I-1>::bsorb(std::get<crow>( std::get<crow>(A)), std::get<crow>(v),rhs,w);
using M =
typename std::remove_cv<
typename std::remove_reference<
decltype(std::get<crow>( std::get<crow>(A)))
>::type
>::type;
algmeta_itsteps<I-1,M>::bsorb(std::get<crow>( std::get<crow>(A)), std::get<crow>(v),rhs,w);
std::get<crow>(x).axpy(w,std::get<crow>(v));
MultiTypeBlockMatrix_Solver<I,crow-1,remain_row-1>::bsorb(A,x,v,b,w); //next row
}
......@@ -466,7 +484,13 @@ namespace Dune {
MultiTypeBlockMatrix_Solver_Col<I,crow,0, TMatrix::M()>::calc_rhs(A,x,v,rhs,w); // calculate right side of equation
//solve on blocklevel I-1
algmeta_itsteps<I-1>::dbjac(std::get<crow>( std::get<crow>(A)), std::get<crow>(v),rhs,w);
using M =
typename std::remove_cv<
typename std::remove_reference<
decltype(std::get<crow>( std::get<crow>(A)))
>::type
>::type;
algmeta_itsteps<I-1,M>::dbjac(std::get<crow>( std::get<crow>(A)), std::get<crow>(v),rhs,w);
MultiTypeBlockMatrix_Solver<I,crow+1,remain_row-1>::dbjac(A,x,v,b,w); //next row
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment