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;
61 typedef typename FirstRow::field_type field_type;
66 return 1+
sizeof...(Args);
72 return FirstRow::size();
91 template<
size_type index >
93 operator[] ([[maybe_unused]]
const std::integral_constant< size_type, index > indexVariable)
94 -> decltype(std::get<index>(*
this))
96 return std::get<index>(*
this);
104 template<
size_type index >
106 operator[] ([[maybe_unused]]
const std::integral_constant< size_type, index > indexVariable)
const
107 -> decltype(std::get<index>(*
this))
109 return std::get<index>(*
this);
117 using namespace Dune::Hybrid;
184 template<
typename X,
typename Y>
185 void mv (
const X& x, Y& y)
const {
186 static_assert(X::size() ==
M(),
"length of x does not match row length");
187 static_assert(Y::size() ==
N(),
"length of y does not match row count");
194 template<
typename X,
typename Y>
195 void umv (
const X& x, Y& y)
const {
196 static_assert(X::size() ==
M(),
"length of x does not match row length");
197 static_assert(Y::size() ==
N(),
"length of y does not match row count");
198 using namespace Dune::Hybrid;
200 using namespace Dune::Hybrid;
202 (*this)[i][j].umv(x[j], y[i]);
209 template<
typename X,
typename Y>
210 void mmv (
const X& x, Y& y)
const {
211 static_assert(X::size() ==
M(),
"length of x does not match row length");
212 static_assert(Y::size() ==
N(),
"length of y does not match row count");
213 using namespace Dune::Hybrid;
215 using namespace Dune::Hybrid;
217 (*this)[i][j].mmv(x[j], y[i]);
224 template<
typename AlphaType,
typename X,
typename Y>
225 void usmv (
const AlphaType& alpha,
const X& x, Y& y)
const {
226 static_assert(X::size() ==
M(),
"length of x does not match row length");
227 static_assert(Y::size() ==
N(),
"length of y does not match row count");
228 using namespace Dune::Hybrid;
230 using namespace Dune::Hybrid;
232 (*this)[i][j].usmv(alpha, x[j], y[i]);
239 template<
typename X,
typename Y>
240 void mtv (
const X& x, Y& y)
const {
241 static_assert(X::size() ==
N(),
"length of x does not match number of rows");
242 static_assert(Y::size() ==
M(),
"length of y does not match number of columns");
249 template<
typename X,
typename Y>
250 void umtv (
const X& x, Y& y)
const {
251 static_assert(X::size() ==
N(),
"length of x does not match number of rows");
252 static_assert(Y::size() ==
M(),
"length of y does not match number of columns");
253 using namespace Dune::Hybrid;
255 using namespace Dune::Hybrid;
257 (*this)[j][i].umtv(x[j], y[i]);
264 template<
typename X,
typename Y>
265 void mmtv (
const X& x, Y& y)
const {
266 static_assert(X::size() ==
N(),
"length of x does not match number of rows");
267 static_assert(Y::size() ==
M(),
"length of y does not match number of columns");
268 using namespace Dune::Hybrid;
270 using namespace Dune::Hybrid;
272 (*this)[j][i].mmtv(x[j], y[i]);
279 template<
typename X,
typename Y>
280 void usmtv (
const field_type& alpha,
const X& x, Y& y)
const {
281 static_assert(X::size() ==
N(),
"length of x does not match number of rows");
282 static_assert(Y::size() ==
M(),
"length of y does not match number of columns");
283 using namespace Dune::Hybrid;
285 using namespace Dune::Hybrid;
287 (*this)[j][i].usmtv(alpha, x[j], y[i]);
294 template<
typename X,
typename Y>
295 void umhv (
const X& x, Y& y)
const {
296 static_assert(X::size() ==
N(),
"length of x does not match number of rows");
297 static_assert(Y::size() ==
M(),
"length of y does not match number of columns");
298 using namespace Dune::Hybrid;
300 using namespace Dune::Hybrid;
302 (*this)[j][i].umhv(x[j], y[i]);
309 template<
typename X,
typename Y>
310 void mmhv (
const X& x, Y& y)
const {
311 static_assert(X::size() ==
N(),
"length of x does not match number of rows");
312 static_assert(Y::size() ==
M(),
"length of y does not match number of columns");
313 using namespace Dune::Hybrid;
315 using namespace Dune::Hybrid;
317 (*this)[j][i].mmhv(x[j], y[i]);
324 template<
typename X,
typename Y>
325 void usmhv (
const field_type& alpha,
const X& x, Y& y)
const {
326 static_assert(X::size() ==
N(),
"length of x does not match number of rows");
327 static_assert(Y::size() ==
M(),
"length of y does not match number of columns");
328 using namespace Dune::Hybrid;
330 using namespace Dune::Hybrid;
332 (*this)[j][i].usmhv(alpha, x[j], y[i]);
344 typename FieldTraits<field_type_00>::real_type sum=0;
350 sum += (*this)[i][j].frobenius_norm2();
368 typename FieldTraits<field_type_00>::real_type norm=0;
373 typename FieldTraits<field_type_00>::real_type sum=0;
376 sum += (*this)[i][j].infinity_norm();
378 norm =
max(sum, norm);
389 typename FieldTraits<field_type_00>::real_type norm=0;
394 typename FieldTraits<field_type_00>::real_type sum=0;
397 sum += (*this)[i][j].infinity_norm_real();
399 norm =
max(sum, norm);
411 template <
typename... Rows>
412 struct FieldTraits< MultiTypeBlockMatrix<Rows...> >
414 using field_type =
typename MultiTypeBlockMatrix<Rows...>::field_type;
415 using real_type =
typename FieldTraits<field_type>::real_type;
424 template<
typename T1,
typename... Args>
428 using namespace Dune::Hybrid;
430 using namespace Dune::Hybrid;
432 s <<
"\t(" << i <<
", " << j <<
"): \n" << m[i][j];
440 template<
int I,
typename M>
441 struct algmeta_itsteps;
449 template<
int I,
int crow,
int ccol,
int remain_col>
455 template <
typename Trhs,
typename TVector,
typename TMatrix,
typename K>
456 static void calc_rhs(
const TMatrix& A, TVector& x, TVector& v, Trhs& b,
const K& w) {
457 std::get<ccol>( std::get<crow>(A) ).mmv( std::get<ccol>(x), b );
462 template<
int I,
int crow,
int ccol>
463 class MultiTypeBlockMatrix_Solver_Col<I,crow,ccol,0> {
465 template <
typename Trhs,
typename TVector,
typename TMatrix,
typename K>
466 static void calc_rhs(
const TMatrix&, TVector&, TVector&, Trhs&,
const K&) {}
477 template<
int I,
int crow,
int remain_row>
484 template <
typename TVector,
typename TMatrix,
typename K>
485 static void dbgs(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
492 template <
typename TVector,
typename TMatrix,
typename K>
493 static void dbgs(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
494 auto rhs = std::get<crow> (b);
499 typename std::remove_cv<
500 typename std::remove_reference<
501 decltype(std::get<crow>( std::get<crow>(A)))
513 template <
typename TVector,
typename TMatrix,
typename K>
514 static void bsorf(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
520 template <
typename TVector,
typename TMatrix,
typename K>
521 static void bsorf(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
522 auto rhs = std::get<crow> (b);
527 typename std::remove_cv<
528 typename std::remove_reference<
529 decltype(std::get<crow>( std::get<crow>(A)))
533 std::get<crow>(x).axpy(w,std::get<crow>(v));
540 template <
typename TVector,
typename TMatrix,
typename K>
541 static void bsorb(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
547 template <
typename TVector,
typename TMatrix,
typename K>
548 static void bsorb(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
549 auto rhs = std::get<crow> (b);
554 typename std::remove_cv<
555 typename std::remove_reference<
556 decltype(std::get<crow>( std::get<crow>(A)))
560 std::get<crow>(x).axpy(w,std::get<crow>(v));
568 template <
typename TVector,
typename TMatrix,
typename K>
569 static void dbjac(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
575 template <
typename TVector,
typename TMatrix,
typename K>
576 static void dbjac(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
577 auto rhs = std::get<crow> (b);
582 typename std::remove_cv<
583 typename std::remove_reference<
584 decltype(std::get<crow>( std::get<crow>(A)))
595 template<
int I,
int crow>
596 class MultiTypeBlockMatrix_Solver<I,crow,0> {
598 template <
typename TVector,
typename TMatrix,
typename K>
599 static void dbgs(
const TMatrix&, TVector&, TVector&,
600 const TVector&,
const K&) {}
602 template <
typename TVector,
typename TMatrix,
typename K>
603 static void bsorf(
const TMatrix&, TVector&, TVector&,
604 const TVector&,
const K&) {}
606 template <
typename TVector,
typename TMatrix,
typename K>
607 static void bsorb(
const TMatrix&, TVector&, TVector&,
608 const TVector&,
const K&) {}
610 template <
typename TVector,
typename TMatrix,
typename K>
611 static void dbjac(
const TMatrix&, TVector&, TVector&,
612 const TVector&,
const K&) {}
623 template <
size_t i,
typename... Args>
624 struct tuple_element<i,
Dune::MultiTypeBlockMatrix<Args...> >
626 using type =
typename std::tuple_element<i, std::tuple<Args...> >
::type;
633 template <
typename... Args>
634 struct tuple_size<
Dune::MultiTypeBlockMatrix<Args...> >
635 : std::integral_constant<std::size_t, sizeof...(Args)>
part of solvers for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition: multitypeblockmatrix.hh:450
solver for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition: multitypeblockmatrix.hh:478
A Matrix class to support different block types.
Definition: multitypeblockmatrix.hh:46
constexpr index_constant< 0 > _0
Compile time index with value 0.
Definition: indices.hh:52
std::integral_constant< std::size_t, i > index_constant
An index constant with value i.
Definition: indices.hh:29
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:569
void mv(const X &x, Y &y) const
y = A x
Definition: multitypeblockmatrix.hh:185
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: multitypeblockmatrix.hh:265
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition: multitypeblockmatrix.hh:310
static void dbgs(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:485
static void bsorb(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:541
MultiTypeBlockMatrix & operator/=(const field_type &k)
vector space division by scalar
Definition: multitypeblockmatrix.hh:141
void usmhv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: multitypeblockmatrix.hh:325
void usmv(const AlphaType &alpha, const X &x, Y &y) const
y += alpha A x
Definition: multitypeblockmatrix.hh:225
void umhv(const X &x, Y &y) const
y += A^H x
Definition: multitypeblockmatrix.hh:295
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:280
void operator=(const T &newval)
Definition: multitypeblockmatrix.hh:116
void umv(const X &x, Y &y) const
y += A x
Definition: multitypeblockmatrix.hh:195
static void calc_rhs(const TMatrix &A, TVector &x, TVector &v, Trhs &b, const K &w)
Definition: multitypeblockmatrix.hh:456
static void bsorf(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:514
auto infinity_norm_real() const
Bastardized version of the infinity-norm / row-sum norm.
Definition: multitypeblockmatrix.hh:385
auto frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: multitypeblockmatrix.hh:341
std::size_t size_type
Type used for sizes.
Definition: multitypeblockmatrix.hh:59
auto operator[]([[maybe_unused]] const std::integral_constant< size_type, index > indexVariable) -> decltype(std::get< index >(*this))
Random-access operator.
Definition: multitypeblockmatrix.hh:93
auto infinity_norm() const
Bastardized version of the infinity-norm / row-sum norm.
Definition: multitypeblockmatrix.hh:364
FieldTraits< field_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: multitypeblockmatrix.hh:358
MultiTypeBlockMatrix & operator+=(const MultiTypeBlockMatrix &b)
Add the entries of another matrix to this one.
Definition: multitypeblockmatrix.hh:157
void mtv(const X &x, Y &y) const
y = A^T x
Definition: multitypeblockmatrix.hh:240
static constexpr size_type M()
Return the number of matrix columns.
Definition: multitypeblockmatrix.hh:70
void mmv(const X &x, Y &y) const
y -= A x
Definition: multitypeblockmatrix.hh:210
void umtv(const X &x, Y &y) const
y += A^T x
Definition: multitypeblockmatrix.hh:250
MultiTypeBlockMatrix & operator*=(const field_type &k)
vector space multiplication with scalar
Definition: multitypeblockmatrix.hh:130
MultiTypeBlockMatrix & operator-=(const MultiTypeBlockMatrix &b)
Subtract the entries of another matrix from this one.
Definition: multitypeblockmatrix.hh:172
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