3 #ifndef DUNE_ISTL_MULTITYPEBLOCKMATRIX_HH 4 #define DUNE_ISTL_MULTITYPEBLOCKMATRIX_HH 10 #include <dune/common/hybridutilities.hh> 17 template<
typename FirstRow,
typename... Args>
18 class MultiTypeBlockMatrix;
20 template<
int I,
int crow,
int remain_row>
41 template<
typename FirstRow,
typename... Args>
43 :
public std::tuple<FirstRow, Args...>
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))
92 DUNE_UNUSED_PARAMETER(indexVariable);
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))
105 DUNE_UNUSED_PARAMETER(indexVariable);
106 return std::get<index>(*this);
114 using namespace Dune::Hybrid;
115 auto size = index_constant<1+
sizeof...(Args)>();
119 forEach(integralRange(
size), [&](
auto&& i) {
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;
141 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
142 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
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;
155 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
156 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
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;
169 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
170 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
171 (*this)[i][j].usmv(alpha, x[j], y[i]);
185 template<
typename T1,
typename... Args>
187 auto N = index_constant<MultiTypeBlockMatrix<T1,Args...>::N()>();
188 auto M = index_constant<MultiTypeBlockMatrix<T1,Args...>::M()>();
189 using namespace Dune::Hybrid;
190 forEach(integralRange(
N), [&](
auto&& i) {
191 forEach(integralRange(
M), [&](
auto&& j) {
192 s <<
"\t(" << i <<
", " << j <<
"): \n" << m[i][j];
200 template<
int I,
typename M>
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&) {}
void usmv(const AlphaType &alpha, const X &x, Y &y) const
y += alpha A x
Definition: multitypeblockmatrix.hh:165
Definition: allocator.hh:7
static constexpr std::size_t size()
Return the number of matrix rows.
Definition: multitypeblockmatrix.hh:61
void umv(const X &x, Y &y) const
y += A x
Definition: multitypeblockmatrix.hh:137
MultiTypeBlockMatrix< FirstRow, Args... > type
Definition: multitypeblockmatrix.hh:50
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
std::ostream & operator<<(std::ostream &s, const BlockVector< K, A > &v)
Send BlockVector to an output stream.
Definition: bvector.hh:620
void mv(const X &x, Y &y) const
y = A x
Definition: multitypeblockmatrix.hh:127
static constexpr std::size_t M()
Return the number of matrix columns.
Definition: multitypeblockmatrix.hh:67
auto operator[](const std::integral_constant< std::size_t, index > indexVariable) -> decltype(std::get< index >(*this))
Random-access operator.
Definition: multitypeblockmatrix.hh:90
static void bsorb(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:431
static void bsorb(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:308
static void dbjac(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition: multitypeblockmatrix.hh:371
static void dbgs(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition: multitypeblockmatrix.hh:359
static void bsorb(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition: multitypeblockmatrix.hh:367
static void dbgs(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:253
static void calc_rhs(const TMatrix &, TVector &, TVector &, Trhs &, const K &)
Definition: multitypeblockmatrix.hh:226
static void bsorb(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:301
static void dbjac(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:336
A Matrix class to support different block types.
Definition: matrixutils.hh:121
FirstRow::field_type field_type
Definition: multitypeblockmatrix.hh:52
void operator=(const T &newval)
Definition: multitypeblockmatrix.hh:113
static constexpr std::size_t N()
Return the number of matrix rows.
Definition: multitypeblockmatrix.hh:55
static void calc_rhs(const TMatrix &A, TVector &x, TVector &v, Trhs &b, const K &w)
Definition: multitypeblockmatrix.hh:216
static void dbgs(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:245
static void dbgs(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:373
part of solvers for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition: multitypeblockmatrix.hh:210
void mmv(const X &x, Y &y) const
y -= A x
Definition: multitypeblockmatrix.hh:151
static void dbjac(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:329
static void bsorf(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:401
static void dbjac(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:461
static void bsorf(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:274
solver for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition: multitypeblockmatrix.hh:21
static void bsorf(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition: multitypeblockmatrix.hh:363
static void bsorf(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:281