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;
66 return 1+
sizeof...(Args);
74 [[deprecated(
"Use method 'N' instead")]]
77 return 1+
sizeof...(Args);
83 return FirstRow::size();
102 template<
size_type index >
104 operator[] ([[maybe_unused]]
const std::integral_constant< size_type, index > indexVariable)
105 ->
decltype(std::get<index>(*
this))
107 return std::get<index>(*
this);
115 template<
size_type index >
117 operator[] ([[maybe_unused]]
const std::integral_constant< size_type, index > indexVariable)
const
118 ->
decltype(std::get<index>(*
this))
120 return std::get<index>(*
this);
128 using namespace Dune::Hybrid;
129 auto size = index_constant<1+
sizeof...(Args)>();
133 forEach(integralRange(
size), [&](
auto&& i) {
143 auto size = index_constant<N()>();
144 Hybrid::forEach(Hybrid::integralRange(
size), [&](
auto&& i) {
154 auto size = index_constant<N()>();
155 Hybrid::forEach(Hybrid::integralRange(
size), [&](
auto&& i) {
170 auto size = index_constant<N()>();
171 Hybrid::forEach(Hybrid::integralRange(
size), [&](
auto&& i) {
185 auto size = index_constant<N()>();
186 Hybrid::forEach(Hybrid::integralRange(
size), [&](
auto&& i) {
195 template<
typename X,
typename Y>
196 void mv (
const X& x, Y& y)
const {
197 static_assert(X::size() ==
M(),
"length of x does not match row length");
198 static_assert(Y::size() ==
N(),
"length of y does not match row count");
205 template<
typename X,
typename Y>
206 void umv (
const X& x, Y& y)
const {
207 static_assert(X::size() ==
M(),
"length of x does not match row length");
208 static_assert(Y::size() ==
N(),
"length of y does not match row count");
209 using namespace Dune::Hybrid;
210 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
211 using namespace Dune::Hybrid;
212 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
213 (*this)[i][j].umv(x[j], y[i]);
220 template<
typename X,
typename Y>
221 void mmv (
const X& x, Y& y)
const {
222 static_assert(X::size() ==
M(),
"length of x does not match row length");
223 static_assert(Y::size() ==
N(),
"length of y does not match row count");
224 using namespace Dune::Hybrid;
225 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
226 using namespace Dune::Hybrid;
227 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
228 (*this)[i][j].mmv(x[j], y[i]);
235 template<
typename AlphaType,
typename X,
typename Y>
236 void usmv (
const AlphaType& alpha,
const X& x, Y& y)
const {
237 static_assert(X::size() ==
M(),
"length of x does not match row length");
238 static_assert(Y::size() ==
N(),
"length of y does not match row count");
239 using namespace Dune::Hybrid;
240 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
241 using namespace Dune::Hybrid;
242 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
243 (*this)[i][j].usmv(alpha, x[j], y[i]);
250 template<
typename X,
typename Y>
251 void mtv (
const X& x, Y& y)
const {
252 static_assert(X::size() ==
N(),
"length of x does not match number of rows");
253 static_assert(Y::size() ==
M(),
"length of y does not match number of columns");
260 template<
typename X,
typename Y>
261 void umtv (
const X& x, Y& y)
const {
262 static_assert(X::size() ==
N(),
"length of x does not match number of rows");
263 static_assert(Y::size() ==
M(),
"length of y does not match number of columns");
264 using namespace Dune::Hybrid;
265 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
266 using namespace Dune::Hybrid;
267 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
268 (*this)[j][i].umtv(x[j], y[i]);
275 template<
typename X,
typename Y>
276 void mmtv (
const X& x, Y& y)
const {
277 static_assert(X::size() ==
N(),
"length of x does not match number of rows");
278 static_assert(Y::size() ==
M(),
"length of y does not match number of columns");
279 using namespace Dune::Hybrid;
280 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
281 using namespace Dune::Hybrid;
282 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
283 (*this)[j][i].mmtv(x[j], y[i]);
290 template<
typename X,
typename Y>
292 static_assert(X::size() ==
N(),
"length of x does not match number of rows");
293 static_assert(Y::size() ==
M(),
"length of y does not match number of columns");
294 using namespace Dune::Hybrid;
295 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
296 using namespace Dune::Hybrid;
297 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
298 (*this)[j][i].usmtv(alpha, x[j], y[i]);
305 template<
typename X,
typename Y>
306 void umhv (
const X& x, Y& y)
const {
307 static_assert(X::size() ==
N(),
"length of x does not match number of rows");
308 static_assert(Y::size() ==
M(),
"length of y does not match number of columns");
309 using namespace Dune::Hybrid;
310 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
311 using namespace Dune::Hybrid;
312 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
313 (*this)[j][i].umhv(x[j], y[i]);
320 template<
typename X,
typename Y>
321 void mmhv (
const X& x, Y& y)
const {
322 static_assert(X::size() ==
N(),
"length of x does not match number of rows");
323 static_assert(Y::size() ==
M(),
"length of y does not match number of columns");
324 using namespace Dune::Hybrid;
325 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
326 using namespace Dune::Hybrid;
327 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
328 (*this)[j][i].mmhv(x[j], y[i]);
335 template<
typename X,
typename Y>
337 static_assert(X::size() ==
N(),
"length of x does not match number of rows");
338 static_assert(Y::size() ==
M(),
"length of y does not match number of columns");
339 using namespace Dune::Hybrid;
340 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
341 using namespace Dune::Hybrid;
342 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
343 (*this)[j][i].usmhv(alpha, x[j], y[i]);
354 using field_type_00 =
typename std::decay_t<
decltype((*this)[Indices::_0][Indices::_0])>
::field_type;
355 typename FieldTraits<field_type_00>::real_type sum=0;
357 auto rows = index_constant<N()>();
358 Hybrid::forEach(Hybrid::integralRange(rows), [&](
auto&& i) {
359 auto cols = index_constant<MultiTypeBlockMatrix<FirstRow, Args...>::M()>();
360 Hybrid::forEach(Hybrid::integralRange(cols), [&](
auto&& j) {
361 sum += (*this)[i][j].frobenius_norm2();
377 using field_type_00 =
typename std::decay_t<
decltype((*this)[Indices::_0][Indices::_0])>
::field_type;
379 typename FieldTraits<field_type_00>::real_type norm=0;
381 auto rows = index_constant<N()>();
382 Hybrid::forEach(Hybrid::integralRange(rows), [&](
auto&& i) {
384 typename FieldTraits<field_type_00>::real_type sum=0;
385 auto cols = index_constant<MultiTypeBlockMatrix<FirstRow, Args...>::M()>();
386 Hybrid::forEach(Hybrid::integralRange(cols), [&](
auto&& j) {
387 sum += (*this)[i][j].infinity_norm();
389 norm = max(sum, norm);
398 using field_type_00 =
typename std::decay_t<
decltype((*this)[Indices::_0][Indices::_0])>
::field_type;
400 typename FieldTraits<field_type_00>::real_type norm=0;
402 auto rows = index_constant<N()>();
403 Hybrid::forEach(Hybrid::integralRange(rows), [&](
auto&& i) {
405 typename FieldTraits<field_type_00>::real_type sum=0;
406 auto cols = index_constant<MultiTypeBlockMatrix<FirstRow, Args...>::M()>();
407 Hybrid::forEach(Hybrid::integralRange(cols), [&](
auto&& j) {
408 sum += (*this)[i][j].infinity_norm_real();
410 norm = max(sum, norm);
423 template<
typename T1,
typename... Args>
425 auto N = index_constant<MultiTypeBlockMatrix<T1,Args...>::N()>();
426 auto M = index_constant<MultiTypeBlockMatrix<T1,Args...>::M()>();
427 using namespace Dune::Hybrid;
428 forEach(integralRange(N), [&](
auto&& i) {
429 using namespace Dune::Hybrid;
430 forEach(integralRange(M), [&](
auto&& j) {
431 s <<
"\t(" << i <<
", " << j <<
"): \n" << m[i][j];
439 template<
int I,
typename M>
440 struct algmeta_itsteps;
448 template<
int I,
int crow,
int ccol,
int remain_col>
454 template <
typename Trhs,
typename TVector,
typename TMatrix,
typename K>
455 static void calc_rhs(
const TMatrix& A, TVector& x, TVector& v, Trhs& b,
const K& w) {
456 std::get<ccol>( std::get<crow>(A) ).mmv( std::get<ccol>(x), b );
461 template<
int I,
int crow,
int ccol>
464 template <
typename Trhs,
typename TVector,
typename TMatrix,
typename K>
465 static void calc_rhs(
const TMatrix&, TVector&, TVector&, Trhs&,
const K&) {}
476 template<
int I,
int crow,
int remain_row>
483 template <
typename TVector,
typename TMatrix,
typename K>
484 static void dbgs(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
491 template <
typename TVector,
typename TMatrix,
typename K>
492 static void dbgs(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
493 auto rhs = std::get<crow> (b);
498 typename std::remove_cv<
499 typename std::remove_reference<
500 decltype(std::get<crow>( std::get<crow>(A)))
512 template <
typename TVector,
typename TMatrix,
typename K>
513 static void bsorf(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
519 template <
typename TVector,
typename TMatrix,
typename K>
520 static void bsorf(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
521 auto rhs = std::get<crow> (b);
526 typename std::remove_cv<
527 typename std::remove_reference<
528 decltype(std::get<crow>( std::get<crow>(A)))
532 std::get<crow>(x).axpy(w,std::get<crow>(v));
539 template <
typename TVector,
typename TMatrix,
typename K>
540 static void bsorb(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
546 template <
typename TVector,
typename TMatrix,
typename K>
547 static void bsorb(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
548 auto rhs = std::get<crow> (b);
553 typename std::remove_cv<
554 typename std::remove_reference<
555 decltype(std::get<crow>( std::get<crow>(A)))
559 std::get<crow>(x).axpy(w,std::get<crow>(v));
567 template <
typename TVector,
typename TMatrix,
typename K>
568 static void dbjac(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
574 template <
typename TVector,
typename TMatrix,
typename K>
575 static void dbjac(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
576 auto rhs = std::get<crow> (b);
581 typename std::remove_cv<
582 typename std::remove_reference<
583 decltype(std::get<crow>( std::get<crow>(A)))
594 template<
int I,
int crow>
597 template <
typename TVector,
typename TMatrix,
typename K>
598 static void dbgs(
const TMatrix&, TVector&, TVector&,
599 const TVector&,
const K&) {}
601 template <
typename TVector,
typename TMatrix,
typename K>
602 static void bsorf(
const TMatrix&, TVector&, TVector&,
603 const TVector&,
const K&) {}
605 template <
typename TVector,
typename TMatrix,
typename K>
606 static void bsorb(
const TMatrix&, TVector&, TVector&,
607 const TVector&,
const K&) {}
609 template <
typename TVector,
typename TMatrix,
typename K>
610 static void dbjac(
const TMatrix&, TVector&, TVector&,
611 const TVector&,
const K&) {}
622 template <
size_t i,
typename... Args>
623 struct tuple_element<i,
Dune::MultiTypeBlockMatrix<Args...> >
625 using type =
typename std::tuple_element<i, std::tuple<Args...> >::type;
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
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:568
MultiTypeBlockMatrix & operator+=(const MultiTypeBlockMatrix &b)
Add the entries of another matrix to this one.
Definition: multitypeblockmatrix.hh:168
void mv(const X &x, Y &y) const
y = A x
Definition: multitypeblockmatrix.hh:196
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: multitypeblockmatrix.hh:276
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition: multitypeblockmatrix.hh:321
static void dbgs(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition: multitypeblockmatrix.hh:598
static void dbgs(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:484
static void bsorb(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:540
void usmhv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: multitypeblockmatrix.hh:336
FirstRow::field_type field_type
Definition: multitypeblockmatrix.hh:61
void usmv(const AlphaType &alpha, const X &x, Y &y) const
y += alpha A x
Definition: multitypeblockmatrix.hh:236
MultiTypeBlockMatrix & operator/=(const field_type &k)
vector space division by scalar
Definition: multitypeblockmatrix.hh:152
typename std::tuple_element< i, std::tuple< Args... > >::type type
Definition: multitypeblockmatrix.hh:625
static void dbgs(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:492
static void bsorf(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:520
void umhv(const X &x, Y &y) const
y += A^H x
Definition: multitypeblockmatrix.hh:306
FieldTraits< field_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: multitypeblockmatrix.hh:369
static constexpr size_type N()
Return the number of matrix rows.
Definition: multitypeblockmatrix.hh:64
void usmtv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition: multitypeblockmatrix.hh:291
void operator=(const T &newval)
Definition: multitypeblockmatrix.hh:127
static void calc_rhs(const TMatrix &, TVector &, TVector &, Trhs &, const K &)
Definition: multitypeblockmatrix.hh:465
void umv(const X &x, Y &y) const
y += A x
Definition: multitypeblockmatrix.hh:206
static void calc_rhs(const TMatrix &A, TVector &x, TVector &v, Trhs &b, const K &w)
Definition: multitypeblockmatrix.hh:455
static void bsorf(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:513
static constexpr size_type size()
Return the number of matrix rows.
Definition: multitypeblockmatrix.hh:75
static void dbjac(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition: multitypeblockmatrix.hh:610
auto infinity_norm_real() const
Bastardized version of the infinity-norm / row-sum norm.
Definition: multitypeblockmatrix.hh:396
auto frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: multitypeblockmatrix.hh:352
std::size_t size_type
Type used for sizes.
Definition: multitypeblockmatrix.hh:59
auto infinity_norm() const
Bastardized version of the infinity-norm / row-sum norm.
Definition: multitypeblockmatrix.hh:375
static void bsorb(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition: multitypeblockmatrix.hh:606
MultiTypeBlockMatrix & operator-=(const MultiTypeBlockMatrix &b)
Subtract the entries of another matrix from this one.
Definition: multitypeblockmatrix.hh:183
auto operator[](const std::integral_constant< size_type, index > indexVariable) -> decltype(std::get< index >(*this))
Random-access operator.
Definition: multitypeblockmatrix.hh:104
static void dbjac(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:575
static void bsorf(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition: multitypeblockmatrix.hh:602
void mtv(const X &x, Y &y) const
y = A^T x
Definition: multitypeblockmatrix.hh:251
static constexpr size_type M()
Return the number of matrix columns.
Definition: multitypeblockmatrix.hh:81
void mmv(const X &x, Y &y) const
y -= A x
Definition: multitypeblockmatrix.hh:221
void umtv(const X &x, Y &y) const
y += A^T x
Definition: multitypeblockmatrix.hh:261
static void bsorb(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:547
MultiTypeBlockMatrix & operator*=(const field_type &k)
vector space multiplication with scalar
Definition: multitypeblockmatrix.hh:141
Definition: allocator.hh:11
std::ostream & operator<<(std::ostream &s, const BlockVector< K, A > &v)
Send BlockVector to an output stream.
Definition: bvector.hh:590
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:477
part of solvers for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition: multitypeblockmatrix.hh:449