3#ifndef DUNE_ISTL_MULTITYPEBLOCKMATRIX_HH
4#define DUNE_ISTL_MULTITYPEBLOCKMATRIX_HH
10#include <dune/common/hybridutilities.hh>
12#include "istlexception.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;
59 typedef typename FirstRow::field_type field_type;
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;
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;
209 using namespace Dune::Hybrid;
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;
224 using namespace Dune::Hybrid;
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;
239 using namespace Dune::Hybrid;
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;
264 using namespace Dune::Hybrid;
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;
279 using namespace Dune::Hybrid;
281 (*this)[j][i].mmtv(x[j], y[i]);
288 template<
typename X,
typename Y>
289 void usmtv (
const field_type& alpha,
const X& x, Y& y)
const {
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;
294 using namespace Dune::Hybrid;
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;
309 using namespace Dune::Hybrid;
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;
324 using namespace Dune::Hybrid;
326 (*this)[j][i].mmhv(x[j], y[i]);
333 template<
typename X,
typename Y>
334 void usmhv (
const field_type& alpha,
const X& x, Y& y)
const {
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;
339 using namespace Dune::Hybrid;
341 (*this)[j][i].usmhv(alpha, x[j], y[i]);
353 typename FieldTraits<field_type_00>::real_type sum=0;
359 sum += (*this)[i][j].frobenius_norm2();
377 typename FieldTraits<field_type_00>::real_type norm=0;
382 typename FieldTraits<field_type_00>::real_type sum=0;
385 sum += (*this)[i][j].infinity_norm();
387 norm =
max(sum, norm);
398 typename FieldTraits<field_type_00>::real_type norm=0;
403 typename FieldTraits<field_type_00>::real_type sum=0;
406 sum += (*this)[i][j].infinity_norm_real();
408 norm =
max(sum, norm);
421 template<
typename T1,
typename... Args>
425 using namespace Dune::Hybrid;
427 using namespace Dune::Hybrid;
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&) {}
part of solvers for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition: multitypeblockmatrix.hh:447
solver for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition: multitypeblockmatrix.hh:475
A Matrix class to support different block types.
Definition: multitypeblockmatrix.hh:44
constexpr index_constant< 0 > _0
Compile time index with value 0.
Definition: indices.hh:51
std::integral_constant< std::size_t, i > index_constant
An index constant with value i.
Definition: indices.hh:28
constexpr void forEach(Range &&range, F &&f)
Range based for loop.
Definition: hybridutilities.hh:266
constexpr auto integralRange(const Begin &begin, const End &end)
Create an integral range.
Definition: hybridutilities.hh:182
void bsorb(const M &A, X &x, const Y &b, const K &w)
SSOR step.
Definition: gsetc.hh:644
void dbjac(const M &A, X &x, const Y &b, const K &w)
Jacobi step.
Definition: gsetc.hh:656
void dbgs(const M &A, X &x, const Y &b, const K &w)
GS step.
Definition: gsetc.hh:620
void bsorf(const M &A, X &x, const Y &b, const K &w)
SOR step.
Definition: gsetc.hh:632
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
MultiTypeBlockMatrix & operator+=(const MultiTypeBlockMatrix &b)
Add the entries of another matrix to this one.
Definition: multitypeblockmatrix.hh:166
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 &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
void usmhv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: multitypeblockmatrix.hh:334
void usmv(const AlphaType &alpha, const X &x, Y &y) const
y += alpha A x
Definition: multitypeblockmatrix.hh:234
MultiTypeBlockMatrix & operator/=(const field_type &k)
vector space division by scalar
Definition: multitypeblockmatrix.hh:150
void umhv(const X &x, Y &y) const
y += A^H x
Definition: multitypeblockmatrix.hh:304
FieldTraits< field_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: multitypeblockmatrix.hh:367
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
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
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 infinity_norm() const
Bastardized version of the infinity-norm / row-sum norm.
Definition: multitypeblockmatrix.hh:373
MultiTypeBlockMatrix & operator-=(const MultiTypeBlockMatrix &b)
Subtract the entries of another matrix from this one.
Definition: multitypeblockmatrix.hh:181
auto operator[](const std::integral_constant< size_type, index > indexVariable) -> decltype(std::get< index >(*this))
Random-access operator.
Definition: multitypeblockmatrix.hh:102
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
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
Dune namespace.
Definition: alignedallocator.hh:11