3#ifndef DUNE_ISTL_MULTITYPEBLOCKMATRIX_HH
4#define DUNE_ISTL_MULTITYPEBLOCKMATRIX_HH
10#include <dune/common/hybridutilities.hh>
12#include "istlexception.hh"
17 template<
typename FirstRow,
typename... Args>
18 class MultiTypeBlockMatrix;
20 template<
int I,
int crow,
int remain_row>
21 class MultiTypeBlockMatrix_Solver;
41 template<
typename FirstRow,
typename... Args>
43 :
public std::tuple<FirstRow, Args...>
52 typedef typename FirstRow::field_type field_type;
55 static constexpr std::size_t
N()
57 return 1+
sizeof...(Args);
61 static constexpr std::size_t
size()
63 return 1+
sizeof...(Args);
67 static constexpr std::size_t
M()
69 return FirstRow::size();
88 template< std::
size_t index >
90 operator[] (
const std::integral_constant< std::size_t, index > indexVariable ) ->
decltype(std::get<index>(*
this))
93 return std::get<index>(*
this);
101 template< std::
size_t index >
103 operator[] (
const std::integral_constant< std::size_t, index > indexVariable )
const ->
decltype(std::get<index>(*
this))
106 return std::get<index>(*
this);
114 using namespace Dune::Hybrid;
126 template<
typename X,
typename Y>
127 void mv (
const X& x, Y& y)
const {
128 static_assert(X::size() ==
M(),
"length of x does not match row length");
129 static_assert(Y::size() ==
N(),
"length of y does not match row count");
136 template<
typename X,
typename Y>
137 void umv (
const X& x, Y& y)
const {
138 static_assert(X::size() ==
M(),
"length of x does not match row length");
139 static_assert(Y::size() ==
N(),
"length of y does not match row count");
140 using namespace Dune::Hybrid;
143 (*this)[i][j].umv(x[j], y[i]);
150 template<
typename X,
typename Y>
151 void mmv (
const X& x, Y& y)
const {
152 static_assert(X::size() ==
M(),
"length of x does not match row length");
153 static_assert(Y::size() ==
N(),
"length of y does not match row count");
154 using namespace Dune::Hybrid;
157 (*this)[i][j].mmv(x[j], y[i]);
164 template<
typename AlphaType,
typename X,
typename Y>
165 void usmv (
const AlphaType& alpha,
const X& x, Y& y)
const {
166 static_assert(X::size() ==
M(),
"length of x does not match row length");
167 static_assert(Y::size() ==
N(),
"length of y does not match row count");
168 using namespace Dune::Hybrid;
171 (*this)[i][j].usmv(alpha, x[j], y[i]);
185 template<
typename T1,
typename... Args>
189 using namespace Dune::Hybrid;
192 s <<
"\t(" << i <<
", " << j <<
"): \n" << m[i][j];
200 template<
int I,
typename M>
201 struct algmeta_itsteps;
209 template<
int I,
int crow,
int ccol,
int remain_col>
215 template <
typename Trhs,
typename TVector,
typename TMatrix,
typename K>
216 static void calc_rhs(
const TMatrix& A, TVector& x, TVector& v, Trhs& b,
const K& w) {
217 std::get<ccol>( std::get<crow>(A) ).mmv( std::get<ccol>(x), b );
222 template<
int I,
int crow,
int ccol>
225 template <
typename Trhs,
typename TVector,
typename TMatrix,
typename K>
226 static void calc_rhs(
const TMatrix&, TVector&, TVector&, Trhs&,
const K&) {}
237 template<
int I,
int crow,
int remain_row>
244 template <
typename TVector,
typename TMatrix,
typename K>
245 static void dbgs(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
252 template <
typename TVector,
typename TMatrix,
typename K>
253 static void dbgs(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
254 auto rhs = std::get<crow> (b);
259 typename std::remove_cv<
260 typename std::remove_reference<
261 decltype(std::get<crow>( std::get<crow>(A)))
273 template <
typename TVector,
typename TMatrix,
typename K>
274 static void bsorf(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
280 template <
typename TVector,
typename TMatrix,
typename K>
281 static void bsorf(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
282 auto rhs = std::get<crow> (b);
287 typename std::remove_cv<
288 typename std::remove_reference<
289 decltype(std::get<crow>( std::get<crow>(A)))
293 std::get<crow>(x).axpy(w,std::get<crow>(v));
300 template <
typename TVector,
typename TMatrix,
typename K>
301 static void bsorb(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
307 template <
typename TVector,
typename TMatrix,
typename K>
308 static void bsorb(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
309 auto rhs = std::get<crow> (b);
314 typename std::remove_cv<
315 typename std::remove_reference<
316 decltype(std::get<crow>( std::get<crow>(A)))
320 std::get<crow>(x).axpy(w,std::get<crow>(v));
328 template <
typename TVector,
typename TMatrix,
typename K>
329 static void dbjac(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
335 template <
typename TVector,
typename TMatrix,
typename K>
336 static void dbjac(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
337 auto rhs = std::get<crow> (b);
342 typename std::remove_cv<
343 typename std::remove_reference<
344 decltype(std::get<crow>( std::get<crow>(A)))
355 template<
int I,
int crow>
358 template <
typename TVector,
typename TMatrix,
typename K>
359 static void dbgs(
const TMatrix&, TVector&, TVector&,
360 const TVector&,
const K&) {}
362 template <
typename TVector,
typename TMatrix,
typename K>
363 static void bsorf(
const TMatrix&, TVector&, TVector&,
364 const TVector&,
const K&) {}
366 template <
typename TVector,
typename TMatrix,
typename K>
367 static void bsorb(
const TMatrix&, TVector&, TVector&,
368 const TVector&,
const K&) {}
370 template <
typename TVector,
typename TMatrix,
typename K>
371 static void dbjac(
const TMatrix&, TVector&, TVector&,
372 const TVector&,
const K&) {}
part of solvers for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition: multitypeblockmatrix.hh:210
solver for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition: multitypeblockmatrix.hh:238
A Matrix class to support different block types.
Definition: multitypeblockmatrix.hh:44
std::integral_constant< std::size_t, i > index_constant
An index constant with value i.
Definition: indices.hh:26
constexpr void forEach(Range &&range, F &&f)
Range based for loop.
Definition: hybridutilities.hh:314
constexpr auto integralRange(const Begin &begin, const End &end)
Create an integral range.
Definition: hybridutilities.hh:227
void bsorb(const M &A, X &x, const Y &b, const K &w)
SSOR step.
Definition: gsetc.hh:591
void dbjac(const M &A, X &x, const Y &b, const K &w)
Jacobi step.
Definition: gsetc.hh:603
void dbgs(const M &A, X &x, const Y &b, const K &w)
GS step.
Definition: gsetc.hh:567
void bsorf(const M &A, X &x, const Y &b, const K &w)
SOR step.
Definition: gsetc.hh:579
MultiTypeBlockMatrix< FirstRow, Args... > type
Definition: multitypeblockmatrix.hh:50
static void dbjac(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:329
void mv(const X &x, Y &y) const
y = A x
Definition: multitypeblockmatrix.hh:127
static constexpr std::size_t size()
Return the number of matrix rows.
Definition: multitypeblockmatrix.hh:61
static void dbgs(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:245
static void bsorb(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:301
void usmv(const AlphaType &alpha, const X &x, Y &y) const
y += alpha A x
Definition: multitypeblockmatrix.hh:165
void operator=(const T &newval)
Definition: multitypeblockmatrix.hh:113
void umv(const X &x, Y &y) const
y += A x
Definition: multitypeblockmatrix.hh:137
static void calc_rhs(const TMatrix &A, TVector &x, TVector &v, Trhs &b, const K &w)
Definition: multitypeblockmatrix.hh:216
static void bsorf(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:274
static constexpr std::size_t N()
Return the number of matrix rows.
Definition: multitypeblockmatrix.hh:55
auto operator[](const std::integral_constant< std::size_t, index > indexVariable) -> decltype(std::get< index >(*this))
Random-access operator.
Definition: multitypeblockmatrix.hh:90
void mmv(const X &x, Y &y) const
y -= A x
Definition: multitypeblockmatrix.hh:151
static constexpr std::size_t M()
Return the number of matrix columns.
Definition: multitypeblockmatrix.hh:67
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
Dune namespace.
Definition: alignment.hh:11
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentionally unused function parameters with.
Definition: unused.hh:18