5#ifndef DUNE_ISTL_MULTITYPEBLOCKMATRIX_HH
6#define DUNE_ISTL_MULTITYPEBLOCKMATRIX_HH
12#include <dune/common/hybridutilities.hh>
19 template<
typename FirstRow,
typename... Args>
20 class MultiTypeBlockMatrix;
22 template<
int I,
int crow,
int remain_row>
23 class MultiTypeBlockMatrix_Solver;
43 template<
typename FirstRow,
typename... Args>
45 :
public std::tuple<FirstRow, Args...>
47 using ParentType = std::tuple<FirstRow, Args...>;
51 using ParentType::ParentType;
67 typename FieldTraits<FirstRow>::field_type,
typename FieldTraits<Args>::field_type...>;
74 using real_type = Std::detected_t<std::common_type_t,
75 typename FieldTraits<FirstRow>::real_type,
typename FieldTraits<Args>::real_type...>;
79 static_assert (
sizeof...(Args) == 0 or
80 not (std::is_same_v<field_type, Std::nonesuch> or std::is_same_v<real_type, Std::nonesuch>),
81 "No std::common_type implemented for the given field_type/real_type of the Args. Please provide an implementation and check that a FieldTraits class is present for your type.");
86 return 1+
sizeof...(Args);
92 return FirstRow::size();
111 template<
size_type index >
113 operator[] ([[maybe_unused]]
const std::integral_constant< size_type, index > indexVariable)
114 ->
decltype(std::get<index>(*
this))
116 return std::get<index>(*
this);
124 template<
size_type index >
126 operator[] ([[maybe_unused]]
const std::integral_constant< size_type, index > indexVariable)
const
127 ->
decltype(std::get<index>(*
this))
129 return std::get<index>(*
this);
137 using namespace Dune::Hybrid;
138 auto size = index_constant<1+
sizeof...(Args)>();
142 forEach(integralRange(size), [&](
auto&& i) {
152 auto size = index_constant<N()>();
153 Hybrid::forEach(Hybrid::integralRange(size), [&](
auto&& i) {
163 auto size = index_constant<N()>();
164 Hybrid::forEach(Hybrid::integralRange(size), [&](
auto&& i) {
179 auto size = index_constant<N()>();
180 Hybrid::forEach(Hybrid::integralRange(size), [&](
auto&& i) {
194 auto size = index_constant<N()>();
195 Hybrid::forEach(Hybrid::integralRange(size), [&](
auto&& i) {
204 template<
typename X,
typename Y>
205 void mv (
const X& x, Y& y)
const {
206 static_assert(X::size() ==
M(),
"length of x does not match row length");
207 static_assert(Y::size() ==
N(),
"length of y does not match row count");
214 template<
typename X,
typename Y>
215 void umv (
const X& x, Y& y)
const {
216 static_assert(X::size() ==
M(),
"length of x does not match row length");
217 static_assert(Y::size() ==
N(),
"length of y does not match row count");
218 using namespace Dune::Hybrid;
219 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
220 using namespace Dune::Hybrid;
221 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
222 (*this)[i][j].umv(x[j], y[i]);
229 template<
typename X,
typename Y>
230 void mmv (
const X& x, Y& y)
const {
231 static_assert(X::size() ==
M(),
"length of x does not match row length");
232 static_assert(Y::size() ==
N(),
"length of y does not match row count");
233 using namespace Dune::Hybrid;
234 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
235 using namespace Dune::Hybrid;
236 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
237 (*this)[i][j].mmv(x[j], y[i]);
244 template<
typename AlphaType,
typename X,
typename Y>
245 void usmv (
const AlphaType& alpha,
const X& x, Y& y)
const {
246 static_assert(X::size() ==
M(),
"length of x does not match row length");
247 static_assert(Y::size() ==
N(),
"length of y does not match row count");
248 using namespace Dune::Hybrid;
249 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
250 using namespace Dune::Hybrid;
251 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
252 (*this)[i][j].usmv(alpha, x[j], y[i]);
259 template<
typename X,
typename Y>
260 void mtv (
const X& x, Y& y)
const {
261 static_assert(X::size() ==
N(),
"length of x does not match number of rows");
262 static_assert(Y::size() ==
M(),
"length of y does not match number of columns");
269 template<
typename X,
typename Y>
270 void umtv (
const X& x, Y& y)
const {
271 static_assert(X::size() ==
N(),
"length of x does not match number of rows");
272 static_assert(Y::size() ==
M(),
"length of y does not match number of columns");
273 using namespace Dune::Hybrid;
274 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
275 using namespace Dune::Hybrid;
276 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
277 (*this)[j][i].umtv(x[j], y[i]);
284 template<
typename X,
typename Y>
285 void mmtv (
const X& x, Y& y)
const {
286 static_assert(X::size() ==
N(),
"length of x does not match number of rows");
287 static_assert(Y::size() ==
M(),
"length of y does not match number of columns");
288 using namespace Dune::Hybrid;
289 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
290 using namespace Dune::Hybrid;
291 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
292 (*this)[j][i].mmtv(x[j], y[i]);
299 template<
typename X,
typename Y>
301 static_assert(X::size() ==
N(),
"length of x does not match number of rows");
302 static_assert(Y::size() ==
M(),
"length of y does not match number of columns");
303 using namespace Dune::Hybrid;
304 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
305 using namespace Dune::Hybrid;
306 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
307 (*this)[j][i].usmtv(alpha, x[j], y[i]);
314 template<
typename X,
typename Y>
315 void umhv (
const X& x, Y& y)
const {
316 static_assert(X::size() ==
N(),
"length of x does not match number of rows");
317 static_assert(Y::size() ==
M(),
"length of y does not match number of columns");
318 using namespace Dune::Hybrid;
319 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
320 using namespace Dune::Hybrid;
321 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
322 (*this)[j][i].umhv(x[j], y[i]);
329 template<
typename X,
typename Y>
330 void mmhv (
const X& x, Y& y)
const {
331 static_assert(X::size() ==
N(),
"length of x does not match number of rows");
332 static_assert(Y::size() ==
M(),
"length of y does not match number of columns");
333 using namespace Dune::Hybrid;
334 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
335 using namespace Dune::Hybrid;
336 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
337 (*this)[j][i].mmhv(x[j], y[i]);
344 template<
typename X,
typename Y>
346 static_assert(X::size() ==
N(),
"length of x does not match number of rows");
347 static_assert(Y::size() ==
M(),
"length of y does not match number of columns");
348 using namespace Dune::Hybrid;
349 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
350 using namespace Dune::Hybrid;
351 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
352 (*this)[j][i].usmhv(alpha, x[j], y[i]);
365 auto rows = index_constant<N()>();
366 Hybrid::forEach(Hybrid::integralRange(rows), [&](
auto&& i) {
367 auto cols = index_constant<MultiTypeBlockMatrix<FirstRow, Args...>::M()>();
368 Hybrid::forEach(Hybrid::integralRange(cols), [&](
auto&& j) {
369 sum += (*this)[i][j].frobenius_norm2();
388 auto rows = index_constant<N()>();
389 Hybrid::forEach(Hybrid::integralRange(rows), [&](
auto&& i) {
391 auto cols = index_constant<MultiTypeBlockMatrix<FirstRow, Args...>::M()>();
392 Hybrid::forEach(Hybrid::integralRange(cols), [&](
auto&& j) {
393 sum += (*this)[i][j].infinity_norm();
395 norm = max(sum, norm);
407 auto rows = index_constant<N()>();
408 Hybrid::forEach(Hybrid::integralRange(rows), [&](
auto&& i) {
410 auto cols = index_constant<MultiTypeBlockMatrix<FirstRow, Args...>::M()>();
411 Hybrid::forEach(Hybrid::integralRange(cols), [&](
auto&& j) {
412 sum += (*this)[i][j].infinity_norm_real();
414 norm = max(sum, norm);
426 template <
typename... Rows>
439 template<
typename T1,
typename... Args>
441 auto N = index_constant<MultiTypeBlockMatrix<T1,Args...>::N()>();
442 auto M = index_constant<MultiTypeBlockMatrix<T1,Args...>::M()>();
443 using namespace Dune::Hybrid;
444 forEach(integralRange(N), [&](
auto&& i) {
445 using namespace Dune::Hybrid;
446 forEach(integralRange(M), [&](
auto&& j) {
447 s <<
"\t(" << i <<
", " << j <<
"): \n" << m[i][j];
455 template<
int I,
typename M>
456 struct algmeta_itsteps;
464 template<
int I,
int crow,
int ccol,
int remain_col>
470 template <
typename Trhs,
typename TVector,
typename TMatrix,
typename K>
471 static void calc_rhs(
const TMatrix& A, TVector& x, TVector& v, Trhs& b,
const K& w) {
472 std::get<ccol>( std::get<crow>(A) ).mmv( std::get<ccol>(x), b );
477 template<
int I,
int crow,
int ccol>
480 template <
typename Trhs,
typename TVector,
typename TMatrix,
typename K>
481 static void calc_rhs(
const TMatrix&, TVector&, TVector&, Trhs&,
const K&) {}
492 template<
int I,
int crow,
int remain_row>
499 template <
typename TVector,
typename TMatrix,
typename K>
500 static void dbgs(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
507 template <
typename TVector,
typename TMatrix,
typename K>
508 static void dbgs(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
509 auto rhs = std::get<crow> (b);
514 typename std::remove_cv<
515 typename std::remove_reference<
516 decltype(std::get<crow>( std::get<crow>(A)))
528 template <
typename TVector,
typename TMatrix,
typename K>
529 static void bsorf(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
535 template <
typename TVector,
typename TMatrix,
typename K>
536 static void bsorf(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
537 auto rhs = std::get<crow> (b);
542 typename std::remove_cv<
543 typename std::remove_reference<
544 decltype(std::get<crow>( std::get<crow>(A)))
548 std::get<crow>(x).axpy(w,std::get<crow>(v));
555 template <
typename TVector,
typename TMatrix,
typename K>
556 static void bsorb(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
562 template <
typename TVector,
typename TMatrix,
typename K>
563 static void bsorb(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
564 auto rhs = std::get<crow> (b);
569 typename std::remove_cv<
570 typename std::remove_reference<
571 decltype(std::get<crow>( std::get<crow>(A)))
575 std::get<crow>(x).axpy(w,std::get<crow>(v));
583 template <
typename TVector,
typename TMatrix,
typename K>
584 static void dbjac(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
590 template <
typename TVector,
typename TMatrix,
typename K>
591 static void dbjac(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
592 auto rhs = std::get<crow> (b);
597 typename std::remove_cv<
598 typename std::remove_reference<
599 decltype(std::get<crow>( std::get<crow>(A)))
610 template<
int I,
int crow>
613 template <
typename TVector,
typename TMatrix,
typename K>
614 static void dbgs(
const TMatrix&, TVector&, TVector&,
615 const TVector&,
const K&) {}
617 template <
typename TVector,
typename TMatrix,
typename K>
618 static void bsorf(
const TMatrix&, TVector&, TVector&,
619 const TVector&,
const K&) {}
621 template <
typename TVector,
typename TMatrix,
typename K>
622 static void bsorb(
const TMatrix&, TVector&, TVector&,
623 const TVector&,
const K&) {}
625 template <
typename TVector,
typename TMatrix,
typename K>
626 static void dbjac(
const TMatrix&, TVector&, TVector&,
627 const TVector&,
const K&) {}
638 template <
size_t i,
typename... Args>
639 struct tuple_element<i,
Dune::MultiTypeBlockMatrix<Args...> >
641 using type =
typename std::tuple_element<i, std::tuple<Args...> >::type;
648 template <
typename... Args>
649 struct tuple_size<
Dune::MultiTypeBlockMatrix<Args...> >
650 : std::integral_constant<std::size_t, sizeof...(Args)>
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
typename MultiTypeBlockMatrix< Rows... >::field_type field_type
Definition multitypeblockmatrix.hh:429
MultiTypeBlockMatrix< FirstRow, Args... > type
Definition multitypeblockmatrix.hh:56
static void dbjac(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition multitypeblockmatrix.hh:584
real_type infinity_norm_real() const
Bastardized version of the infinity-norm / row-sum norm.
Definition multitypeblockmatrix.hh:402
MultiTypeBlockMatrix & operator+=(const MultiTypeBlockMatrix &b)
Add the entries of another matrix to this one.
Definition multitypeblockmatrix.hh:177
void mv(const X &x, Y &y) const
y = A x
Definition multitypeblockmatrix.hh:205
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition multitypeblockmatrix.hh:285
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition multitypeblockmatrix.hh:330
Std::detected_t< std::common_type_t, typename FieldTraits< FirstRow >::real_type, typename FieldTraits< Args >::real_type... > real_type
The type used for real values.
Definition multitypeblockmatrix.hh:75
static void dbgs(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition multitypeblockmatrix.hh:614
static void dbgs(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition multitypeblockmatrix.hh:500
static void bsorb(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition multitypeblockmatrix.hh:556
void usmhv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition multitypeblockmatrix.hh:345
void usmv(const AlphaType &alpha, const X &x, Y &y) const
y += alpha A x
Definition multitypeblockmatrix.hh:245
MultiTypeBlockMatrix & operator/=(const field_type &k)
vector space division by scalar
Definition multitypeblockmatrix.hh:161
real_type infinity_norm() const
Bastardized version of the infinity-norm / row-sum norm.
Definition multitypeblockmatrix.hh:383
typename std::tuple_element< i, std::tuple< Args... > >::type type
Definition multitypeblockmatrix.hh:641
static void dbgs(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition multitypeblockmatrix.hh:508
static void bsorf(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition multitypeblockmatrix.hh:536
void umhv(const X &x, Y &y) const
y += A^H x
Definition multitypeblockmatrix.hh:315
static constexpr size_type N()
Return the number of matrix rows.
Definition multitypeblockmatrix.hh:84
void usmtv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition multitypeblockmatrix.hh:300
void operator=(const T &newval)
Definition multitypeblockmatrix.hh:136
static void calc_rhs(const TMatrix &, TVector &, TVector &, Trhs &, const K &)
Definition multitypeblockmatrix.hh:481
void umv(const X &x, Y &y) const
y += A x
Definition multitypeblockmatrix.hh:215
static void calc_rhs(const TMatrix &A, TVector &x, TVector &v, Trhs &b, const K &w)
Definition multitypeblockmatrix.hh:471
static void bsorf(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition multitypeblockmatrix.hh:529
static void dbjac(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition multitypeblockmatrix.hh:626
std::size_t size_type
Type used for sizes.
Definition multitypeblockmatrix.hh:59
real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition multitypeblockmatrix.hh:361
real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition multitypeblockmatrix.hh:377
static void bsorb(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition multitypeblockmatrix.hh:622
MultiTypeBlockMatrix & operator-=(const MultiTypeBlockMatrix &b)
Subtract the entries of another matrix from this one.
Definition multitypeblockmatrix.hh:192
typename MultiTypeBlockMatrix< Rows... >::real_type real_type
Definition multitypeblockmatrix.hh:430
auto operator[](const std::integral_constant< size_type, index > indexVariable) -> decltype(std::get< index >(*this))
Random-access operator.
Definition multitypeblockmatrix.hh:113
static void dbjac(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition multitypeblockmatrix.hh:591
static void bsorf(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition multitypeblockmatrix.hh:618
void mtv(const X &x, Y &y) const
y = A^T x
Definition multitypeblockmatrix.hh:260
static constexpr size_type M()
Return the number of matrix columns.
Definition multitypeblockmatrix.hh:90
void mmv(const X &x, Y &y) const
y -= A x
Definition multitypeblockmatrix.hh:230
Std::detected_t< std::common_type_t, typename FieldTraits< FirstRow >::field_type, typename FieldTraits< Args >::field_type... > field_type
The type used for scalars.
Definition multitypeblockmatrix.hh:67
void umtv(const X &x, Y &y) const
y += A^T x
Definition multitypeblockmatrix.hh:270
static void bsorb(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition multitypeblockmatrix.hh:563
MultiTypeBlockMatrix & operator*=(const field_type &k)
vector space multiplication with scalar
Definition multitypeblockmatrix.hh:150
Definition allocator.hh:11
std::ostream & operator<<(std::ostream &s, const BlockVector< K, A > &v)
Send BlockVector to an output stream.
Definition bvector.hh:583
A Matrix class to support different block types.
Definition multitypeblockmatrix.hh:46
static void bsorb(const M &A, X &x, const Y &b, const K &w)
Definition gsetc.hh:461
static void bsorf(const M &A, X &x, const Y &b, const K &w)
Definition gsetc.hh:418
static void dbjac(const M &A, X &x, const Y &b, const K &w)
Definition gsetc.hh:504
static void dbgs(const M &A, X &x, const Y &b, const K &w)
Definition gsetc.hh:378
solver for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition multitypeblockmatrix.hh:493
part of solvers for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition multitypeblockmatrix.hh:465