5#ifndef DUNE_ISTL_MULTITYPEBLOCKMATRIX_HH
6#define DUNE_ISTL_MULTITYPEBLOCKMATRIX_HH
12#include <dune/common/hybridutilities.hh>
14#include "istlexception.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...>;
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;
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;
220 using namespace Dune::Hybrid;
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;
235 using namespace Dune::Hybrid;
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;
250 using namespace Dune::Hybrid;
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;
275 using namespace Dune::Hybrid;
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;
290 using namespace Dune::Hybrid;
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;
305 using namespace Dune::Hybrid;
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;
320 using namespace Dune::Hybrid;
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;
335 using namespace Dune::Hybrid;
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;
350 using namespace Dune::Hybrid;
352 (*this)[j][i].usmhv(alpha, x[j], y[i]);
369 sum += (*this)[i][j].frobenius_norm2();
393 sum += (*this)[i][j].infinity_norm();
395 norm =
max(sum, norm);
412 sum += (*this)[i][j].infinity_norm_real();
414 norm =
max(sum, norm);
426 template <
typename... Rows>
439 template<
typename T1,
typename... Args>
443 using namespace Dune::Hybrid;
445 using namespace Dune::Hybrid;
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>
641 using type =
typename std::tuple_element<i, std::tuple<Args...> >::type;
648 template <
typename... Args>
650 : std::integral_constant<std::size_t, sizeof...(Args)>
part of solvers for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition: multitypeblockmatrix.hh:465
solver for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition: multitypeblockmatrix.hh:493
A Matrix class to support different block types.
Definition: multitypeblockmatrix.hh:46
std::integral_constant< std::size_t, i > index_constant
An index constant with value i.
Definition: indices.hh:29
typename detected_or< nonesuch, Op, Args... >::type detected_t
Returns Op<Args...> if that is valid; otherwise returns nonesuch.
Definition: type_traits.hh:174
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 &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
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
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
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
MultiTypeBlockMatrix & operator-=(const MultiTypeBlockMatrix &b)
Subtract the entries of another matrix from this one.
Definition: multitypeblockmatrix.hh:192
auto operator[](const std::integral_constant< size_type, index > indexVariable) -> decltype(std::get< index >(*this))
Random-access operator.
Definition: multitypeblockmatrix.hh:113
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
MultiTypeBlockMatrix & operator*=(const field_type &k)
vector space multiplication with scalar
Definition: multitypeblockmatrix.hh:150
constexpr void forEach(Range &&range, F &&f)
Range based for loop.
Definition: hybridutilities.hh:256
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
constexpr auto integralRange(const Begin &begin, const End &end)
Create an integral range.
Definition: hybridutilities.hh:172
void bsorb(const M &A, X &x, const Y &b, const K &w)
SSOR step.
Definition: gsetc.hh:646
void dbjac(const M &A, X &x, const Y &b, const K &w)
Jacobi step.
Definition: gsetc.hh:658
void dbgs(const M &A, X &x, const Y &b, const K &w)
GS step.
Definition: gsetc.hh:622
void bsorf(const M &A, X &x, const Y &b, const K &w)
SOR step.
Definition: gsetc.hh:634
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75