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>
21 class MultiTypeBlockMatrix_Solver;
41 template<
typename FirstRow,
typename... Args>
43 :
public std::tuple<FirstRow, Args...>
45 using ParentType = std::tuple<FirstRow, Args...>;
49 using ParentType::ParentType;
64 return 1+
sizeof...(Args);
72 [[deprecated(
"Use method 'N' instead")]]
75 return 1+
sizeof...(Args);
81 return FirstRow::size();
100 template<
size_type index >
102 operator[] ([[maybe_unused]]
const std::integral_constant< size_type, index > indexVariable)
103 -> decltype(std::get<index>(*
this))
105 return std::get<index>(*
this);
113 template<
size_type index >
115 operator[] ([[maybe_unused]]
const std::integral_constant< size_type, index > indexVariable)
const
116 -> decltype(std::get<index>(*
this))
118 return std::get<index>(*
this);
126 using namespace Dune::Hybrid;
127 auto size = index_constant<1+
sizeof...(Args)>();
131 forEach(integralRange(
size), [&](
auto&& i) {
141 auto size = index_constant<N()>();
142 Hybrid::forEach(Hybrid::integralRange(
size), [&](
auto&& i) {
152 auto size = index_constant<N()>();
153 Hybrid::forEach(Hybrid::integralRange(
size), [&](
auto&& i) {
168 auto size = index_constant<N()>();
169 Hybrid::forEach(Hybrid::integralRange(
size), [&](
auto&& i) {
183 auto size = index_constant<N()>();
184 Hybrid::forEach(Hybrid::integralRange(
size), [&](
auto&& i) {
193 template<
typename X,
typename Y>
194 void mv (
const X& x, Y& y)
const {
195 static_assert(X::size() ==
M(),
"length of x does not match row length");
196 static_assert(Y::size() ==
N(),
"length of y does not match row count");
203 template<
typename X,
typename Y>
204 void umv (
const X& x, Y& y)
const {
205 static_assert(X::size() ==
M(),
"length of x does not match row length");
206 static_assert(Y::size() ==
N(),
"length of y does not match row count");
207 using namespace Dune::Hybrid;
208 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
209 using namespace Dune::Hybrid;
210 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
211 (*this)[i][j].umv(x[j], y[i]);
218 template<
typename X,
typename Y>
219 void mmv (
const X& x, Y& y)
const {
220 static_assert(X::size() ==
M(),
"length of x does not match row length");
221 static_assert(Y::size() ==
N(),
"length of y does not match row count");
222 using namespace Dune::Hybrid;
223 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
224 using namespace Dune::Hybrid;
225 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
226 (*this)[i][j].mmv(x[j], y[i]);
233 template<
typename AlphaType,
typename X,
typename Y>
234 void usmv (
const AlphaType& alpha,
const X& x, Y& y)
const {
235 static_assert(X::size() ==
M(),
"length of x does not match row length");
236 static_assert(Y::size() ==
N(),
"length of y does not match row count");
237 using namespace Dune::Hybrid;
238 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
239 using namespace Dune::Hybrid;
240 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
241 (*this)[i][j].usmv(alpha, x[j], y[i]);
248 template<
typename X,
typename Y>
249 void mtv (
const X& x, Y& y)
const {
250 static_assert(X::size() ==
N(),
"length of x does not match number of rows");
251 static_assert(Y::size() ==
M(),
"length of y does not match number of columns");
258 template<
typename X,
typename Y>
259 void umtv (
const X& x, Y& y)
const {
260 static_assert(X::size() ==
N(),
"length of x does not match number of rows");
261 static_assert(Y::size() ==
M(),
"length of y does not match number of columns");
262 using namespace Dune::Hybrid;
263 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
264 using namespace Dune::Hybrid;
265 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
266 (*this)[j][i].umtv(x[j], y[i]);
273 template<
typename X,
typename Y>
274 void mmtv (
const X& x, Y& y)
const {
275 static_assert(X::size() ==
N(),
"length of x does not match number of rows");
276 static_assert(Y::size() ==
M(),
"length of y does not match number of columns");
277 using namespace Dune::Hybrid;
278 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
279 using namespace Dune::Hybrid;
280 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
281 (*this)[j][i].mmtv(x[j], y[i]);
288 template<
typename X,
typename Y>
290 static_assert(X::size() ==
N(),
"length of x does not match number of rows");
291 static_assert(Y::size() ==
M(),
"length of y does not match number of columns");
292 using namespace Dune::Hybrid;
293 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
294 using namespace Dune::Hybrid;
295 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
296 (*this)[j][i].usmtv(alpha, x[j], y[i]);
303 template<
typename X,
typename Y>
304 void umhv (
const X& x, Y& y)
const {
305 static_assert(X::size() ==
N(),
"length of x does not match number of rows");
306 static_assert(Y::size() ==
M(),
"length of y does not match number of columns");
307 using namespace Dune::Hybrid;
308 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
309 using namespace Dune::Hybrid;
310 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
311 (*this)[j][i].umhv(x[j], y[i]);
318 template<
typename X,
typename Y>
319 void mmhv (
const X& x, Y& y)
const {
320 static_assert(X::size() ==
N(),
"length of x does not match number of rows");
321 static_assert(Y::size() ==
M(),
"length of y does not match number of columns");
322 using namespace Dune::Hybrid;
323 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
324 using namespace Dune::Hybrid;
325 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
326 (*this)[j][i].mmhv(x[j], y[i]);
333 template<
typename X,
typename Y>
335 static_assert(X::size() ==
N(),
"length of x does not match number of rows");
336 static_assert(Y::size() ==
M(),
"length of y does not match number of columns");
337 using namespace Dune::Hybrid;
338 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
339 using namespace Dune::Hybrid;
340 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
341 (*this)[j][i].usmhv(alpha, x[j], y[i]);
352 using field_type_00 =
typename std::decay_t<decltype((*
this)[Indices::_0][Indices::_0])>
::field_type;
353 typename FieldTraits<field_type_00>::real_type sum=0;
355 auto rows = index_constant<N()>();
356 Hybrid::forEach(Hybrid::integralRange(rows), [&](
auto&& i) {
357 auto cols = index_constant<MultiTypeBlockMatrix<FirstRow, Args...>::M()>();
358 Hybrid::forEach(Hybrid::integralRange(cols), [&](
auto&& j) {
359 sum += (*this)[i][j].frobenius_norm2();
375 using field_type_00 =
typename std::decay_t<decltype((*
this)[Indices::_0][Indices::_0])>
::field_type;
377 typename FieldTraits<field_type_00>::real_type norm=0;
379 auto rows = index_constant<N()>();
380 Hybrid::forEach(Hybrid::integralRange(rows), [&](
auto&& i) {
382 typename FieldTraits<field_type_00>::real_type sum=0;
383 auto cols = index_constant<MultiTypeBlockMatrix<FirstRow, Args...>::M()>();
384 Hybrid::forEach(Hybrid::integralRange(cols), [&](
auto&& j) {
385 sum += (*this)[i][j].infinity_norm();
387 norm = max(sum, norm);
396 using field_type_00 =
typename std::decay_t<decltype((*
this)[Indices::_0][Indices::_0])>
::field_type;
398 typename FieldTraits<field_type_00>::real_type norm=0;
400 auto rows = index_constant<N()>();
401 Hybrid::forEach(Hybrid::integralRange(rows), [&](
auto&& i) {
403 typename FieldTraits<field_type_00>::real_type sum=0;
404 auto cols = index_constant<MultiTypeBlockMatrix<FirstRow, Args...>::M()>();
405 Hybrid::forEach(Hybrid::integralRange(cols), [&](
auto&& j) {
406 sum += (*this)[i][j].infinity_norm_real();
408 norm = max(sum, norm);
421 template<
typename T1,
typename... Args>
423 auto N = index_constant<MultiTypeBlockMatrix<T1,Args...>::N()>();
424 auto M = index_constant<MultiTypeBlockMatrix<T1,Args...>::M()>();
425 using namespace Dune::Hybrid;
426 forEach(integralRange(
N), [&](
auto&& i) {
427 using namespace Dune::Hybrid;
428 forEach(integralRange(
M), [&](
auto&& j) {
429 s <<
"\t(" << i <<
", " << j <<
"): \n" << m[i][j];
437 template<
int I,
typename M>
438 struct algmeta_itsteps;
446 template<
int I,
int crow,
int ccol,
int remain_col>
452 template <
typename Trhs,
typename TVector,
typename TMatrix,
typename K>
453 static void calc_rhs(
const TMatrix& A, TVector& x, TVector& v, Trhs& b,
const K& w) {
454 std::get<ccol>( std::get<crow>(A) ).mmv( std::get<ccol>(x), b );
459 template<
int I,
int crow,
int ccol>
462 template <
typename Trhs,
typename TVector,
typename TMatrix,
typename K>
463 static void calc_rhs(
const TMatrix&, TVector&, TVector&, Trhs&,
const K&) {}
474 template<
int I,
int crow,
int remain_row>
481 template <
typename TVector,
typename TMatrix,
typename K>
482 static void dbgs(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
489 template <
typename TVector,
typename TMatrix,
typename K>
490 static void dbgs(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
491 auto rhs = std::get<crow> (b);
496 typename std::remove_cv<
497 typename std::remove_reference<
498 decltype(std::get<crow>( std::get<crow>(A)))
510 template <
typename TVector,
typename TMatrix,
typename K>
511 static void bsorf(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
517 template <
typename TVector,
typename TMatrix,
typename K>
518 static void bsorf(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
519 auto rhs = std::get<crow> (b);
524 typename std::remove_cv<
525 typename std::remove_reference<
526 decltype(std::get<crow>( std::get<crow>(A)))
530 std::get<crow>(x).axpy(w,std::get<crow>(v));
537 template <
typename TVector,
typename TMatrix,
typename K>
538 static void bsorb(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
544 template <
typename TVector,
typename TMatrix,
typename K>
545 static void bsorb(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
546 auto rhs = std::get<crow> (b);
551 typename std::remove_cv<
552 typename std::remove_reference<
553 decltype(std::get<crow>( std::get<crow>(A)))
557 std::get<crow>(x).axpy(w,std::get<crow>(v));
565 template <
typename TVector,
typename TMatrix,
typename K>
566 static void dbjac(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
572 template <
typename TVector,
typename TMatrix,
typename K>
573 static void dbjac(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
574 auto rhs = std::get<crow> (b);
579 typename std::remove_cv<
580 typename std::remove_reference<
581 decltype(std::get<crow>( std::get<crow>(A)))
592 template<
int I,
int crow>
595 template <
typename TVector,
typename TMatrix,
typename K>
596 static void dbgs(
const TMatrix&, TVector&, TVector&,
597 const TVector&,
const K&) {}
599 template <
typename TVector,
typename TMatrix,
typename K>
600 static void bsorf(
const TMatrix&, TVector&, TVector&,
601 const TVector&,
const K&) {}
603 template <
typename TVector,
typename TMatrix,
typename K>
604 static void bsorb(
const TMatrix&, TVector&, TVector&,
605 const TVector&,
const K&) {}
607 template <
typename TVector,
typename TMatrix,
typename K>
608 static void dbjac(
const TMatrix&, TVector&, TVector&,
609 const TVector&,
const K&) {}
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
MultiTypeBlockMatrix< FirstRow, Args... > type
Definition: multitypeblockmatrix.hh:54
static void dbjac(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:566
void mv(const X &x, Y &y) const
y = A x
Definition: multitypeblockmatrix.hh:194
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: multitypeblockmatrix.hh:274
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition: multitypeblockmatrix.hh:319
static void dbgs(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition: multitypeblockmatrix.hh:596
static void dbgs(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:482
static void bsorb(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:538
MultiTypeBlockMatrix & operator/=(const field_type &k)
vector space division by scalar
Definition: multitypeblockmatrix.hh:150
void usmhv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: multitypeblockmatrix.hh:334
FirstRow::field_type field_type
Definition: multitypeblockmatrix.hh:59
void usmv(const AlphaType &alpha, const X &x, Y &y) const
y += alpha A x
Definition: multitypeblockmatrix.hh:234
static void dbgs(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:490
static void bsorf(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:518
void umhv(const X &x, Y &y) const
y += A^H x
Definition: multitypeblockmatrix.hh:304
static constexpr size_type N()
Return the number of matrix rows.
Definition: multitypeblockmatrix.hh:62
void usmtv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition: multitypeblockmatrix.hh:289
void operator=(const T &newval)
Definition: multitypeblockmatrix.hh:125
static void calc_rhs(const TMatrix &, TVector &, TVector &, Trhs &, const K &)
Definition: multitypeblockmatrix.hh:463
void umv(const X &x, Y &y) const
y += A x
Definition: multitypeblockmatrix.hh:204
static void calc_rhs(const TMatrix &A, TVector &x, TVector &v, Trhs &b, const K &w)
Definition: multitypeblockmatrix.hh:453
static void bsorf(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:511
static constexpr size_type size()
Return the number of matrix rows.
Definition: multitypeblockmatrix.hh:73
static void dbjac(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition: multitypeblockmatrix.hh:608
auto infinity_norm_real() const
Bastardized version of the infinity-norm / row-sum norm.
Definition: multitypeblockmatrix.hh:394
auto frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: multitypeblockmatrix.hh:350
std::size_t size_type
Type used for sizes.
Definition: multitypeblockmatrix.hh:57
auto operator[]([[maybe_unused]] const std::integral_constant< size_type, index > indexVariable) -> decltype(std::get< index >(*this))
Random-access operator.
Definition: multitypeblockmatrix.hh:102
auto infinity_norm() const
Bastardized version of the infinity-norm / row-sum norm.
Definition: multitypeblockmatrix.hh:373
static void bsorb(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition: multitypeblockmatrix.hh:604
FieldTraits< field_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: multitypeblockmatrix.hh:367
MultiTypeBlockMatrix & operator+=(const MultiTypeBlockMatrix &b)
Add the entries of another matrix to this one.
Definition: multitypeblockmatrix.hh:166
static void dbjac(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:573
static void bsorf(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition: multitypeblockmatrix.hh:600
void mtv(const X &x, Y &y) const
y = A^T x
Definition: multitypeblockmatrix.hh:249
static constexpr size_type M()
Return the number of matrix columns.
Definition: multitypeblockmatrix.hh:79
void mmv(const X &x, Y &y) const
y -= A x
Definition: multitypeblockmatrix.hh:219
void umtv(const X &x, Y &y) const
y += A^T x
Definition: multitypeblockmatrix.hh:259
MultiTypeBlockMatrix & operator*=(const field_type &k)
vector space multiplication with scalar
Definition: multitypeblockmatrix.hh:139
MultiTypeBlockMatrix & operator-=(const MultiTypeBlockMatrix &b)
Subtract the entries of another matrix from this one.
Definition: multitypeblockmatrix.hh:181
static void bsorb(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:545
Definition: allocator.hh:9
std::ostream & operator<<(std::ostream &s, const BlockVector< K, A > &v)
Send BlockVector to an output stream.
Definition: bvector.hh:588
A Matrix class to support different block types.
Definition: multitypeblockmatrix.hh:44
static void bsorb(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:459
static void bsorf(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:416
static void dbjac(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:502
static void dbgs(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:376
solver for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition: multitypeblockmatrix.hh:475
part of solvers for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition: multitypeblockmatrix.hh:447