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...>
60 return 1+
sizeof...(Args);
64 static constexpr
size_type size() DUNE_DEPRECATED_MSG("Use method '
N' instead")
66 return 1+
sizeof...(Args);
72 return FirstRow::size();
91 template<
size_type index >
93 operator[] (
const std::integral_constant< size_type, index > indexVariable ) -> decltype(std::get<index>(*
this))
95 DUNE_UNUSED_PARAMETER(indexVariable);
96 return std::get<index>(*
this);
104 template<
size_type index >
106 operator[] (
const std::integral_constant< size_type, index > indexVariable )
const -> decltype(std::get<index>(*
this))
108 DUNE_UNUSED_PARAMETER(indexVariable);
109 return std::get<index>(*
this);
117 using namespace Dune::Hybrid;
118 auto size = index_constant<1+
sizeof...(Args)>();
122 forEach(integralRange(
size), [&](
auto&& i) {
132 auto size = index_constant<N()>();
133 Hybrid::forEach(Hybrid::integralRange(
size), [&](
auto&& i) {
143 auto size = index_constant<N()>();
144 Hybrid::forEach(Hybrid::integralRange(
size), [&](
auto&& i) {
159 auto size = index_constant<N()>();
160 Hybrid::forEach(Hybrid::integralRange(
size), [&](
auto&& i) {
174 auto size = index_constant<N()>();
175 Hybrid::forEach(Hybrid::integralRange(
size), [&](
auto&& i) {
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;
199 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
200 using namespace Dune::Hybrid;
201 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
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;
214 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
215 using namespace Dune::Hybrid;
216 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
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;
229 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
230 using namespace Dune::Hybrid;
231 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
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;
254 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
255 using namespace Dune::Hybrid;
256 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
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;
269 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
270 using namespace Dune::Hybrid;
271 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
272 (*this)[j][i].mmtv(x[j], y[i]);
279 template<
typename X,
typename Y>
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;
284 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
285 using namespace Dune::Hybrid;
286 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
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;
299 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
300 using namespace Dune::Hybrid;
301 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
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;
314 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
315 using namespace Dune::Hybrid;
316 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
317 (*this)[j][i].mmhv(x[j], y[i]);
324 template<
typename X,
typename Y>
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;
329 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
330 using namespace Dune::Hybrid;
331 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
332 (*this)[j][i].usmhv(alpha, x[j], y[i]);
344 typename FieldTraits<field_type>::real_type sum=0;
346 auto rows = index_constant<N()>();
347 Hybrid::forEach(Hybrid::integralRange(rows), [&](
auto&& i) {
348 auto cols = index_constant<MultiTypeBlockMatrix<FirstRow, Args...>::M()>();
349 Hybrid::forEach(Hybrid::integralRange(cols), [&](
auto&& j) {
350 sum += (*this)[i][j].frobenius_norm2();
368 typename FieldTraits<field_type>::real_type norm=0;
370 auto rows = index_constant<N()>();
371 Hybrid::forEach(Hybrid::integralRange(rows), [&](
auto&& i) {
373 typename FieldTraits<field_type>::real_type sum=0;
374 auto cols = index_constant<MultiTypeBlockMatrix<FirstRow, Args...>::M()>();
375 Hybrid::forEach(Hybrid::integralRange(cols), [&](
auto&& j) {
376 sum += (*this)[i][j].infinity_norm();
378 norm = max(sum, norm);
389 typename FieldTraits<field_type>::real_type norm=0;
391 auto rows = index_constant<N()>();
392 Hybrid::forEach(Hybrid::integralRange(rows), [&](
auto&& i) {
394 typename FieldTraits<field_type>::real_type sum=0;
395 auto cols = index_constant<MultiTypeBlockMatrix<FirstRow, Args...>::M()>();
396 Hybrid::forEach(Hybrid::integralRange(cols), [&](
auto&& j) {
397 sum += (*this)[i][j].infinity_norm_real();
399 norm = max(sum, norm);
412 template<
typename T1,
typename... Args>
414 auto N = index_constant<MultiTypeBlockMatrix<T1,Args...>::N()>();
415 auto M = index_constant<MultiTypeBlockMatrix<T1,Args...>::M()>();
416 using namespace Dune::Hybrid;
417 forEach(integralRange(
N), [&](
auto&& i) {
418 using namespace Dune::Hybrid;
419 forEach(integralRange(
M), [&](
auto&& j) {
420 s <<
"\t(" << i <<
", " << j <<
"): \n" << m[i][j];
428 template<
int I,
typename M>
429 struct algmeta_itsteps;
437 template<
int I,
int crow,
int ccol,
int remain_col>
443 template <
typename Trhs,
typename TVector,
typename TMatrix,
typename K>
444 static void calc_rhs(
const TMatrix& A, TVector& x, TVector& v, Trhs& b,
const K& w) {
445 std::get<ccol>( std::get<crow>(A) ).mmv( std::get<ccol>(x), b );
450 template<
int I,
int crow,
int ccol>
453 template <
typename Trhs,
typename TVector,
typename TMatrix,
typename K>
454 static void calc_rhs(
const TMatrix&, TVector&, TVector&, Trhs&,
const K&) {}
465 template<
int I,
int crow,
int remain_row>
472 template <
typename TVector,
typename TMatrix,
typename K>
473 static void dbgs(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
480 template <
typename TVector,
typename TMatrix,
typename K>
481 static void dbgs(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
482 auto rhs = std::get<crow> (b);
487 typename std::remove_cv<
488 typename std::remove_reference<
489 decltype(std::get<crow>( std::get<crow>(A)))
501 template <
typename TVector,
typename TMatrix,
typename K>
502 static void bsorf(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
508 template <
typename TVector,
typename TMatrix,
typename K>
509 static void bsorf(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
510 auto rhs = std::get<crow> (b);
515 typename std::remove_cv<
516 typename std::remove_reference<
517 decltype(std::get<crow>( std::get<crow>(A)))
521 std::get<crow>(x).axpy(w,std::get<crow>(v));
528 template <
typename TVector,
typename TMatrix,
typename K>
529 static void bsorb(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
535 template <
typename TVector,
typename TMatrix,
typename K>
536 static void bsorb(
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));
556 template <
typename TVector,
typename TMatrix,
typename K>
557 static void dbjac(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
563 template <
typename TVector,
typename TMatrix,
typename K>
564 static void dbjac(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
565 auto rhs = std::get<crow> (b);
570 typename std::remove_cv<
571 typename std::remove_reference<
572 decltype(std::get<crow>( std::get<crow>(A)))
583 template<
int I,
int crow>
586 template <
typename TVector,
typename TMatrix,
typename K>
587 static void dbgs(
const TMatrix&, TVector&, TVector&,
588 const TVector&,
const K&) {}
590 template <
typename TVector,
typename TMatrix,
typename K>
591 static void bsorf(
const TMatrix&, TVector&, TVector&,
592 const TVector&,
const K&) {}
594 template <
typename TVector,
typename TMatrix,
typename K>
595 static void bsorb(
const TMatrix&, TVector&, TVector&,
596 const TVector&,
const K&) {}
598 template <
typename TVector,
typename TMatrix,
typename K>
599 static void dbjac(
const TMatrix&, TVector&, TVector&,
600 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:50
static void dbjac(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:557
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 &, TVector &, TVector &, const TVector &, const K &)
Definition: multitypeblockmatrix.hh:587
static void dbgs(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:473
static void bsorb(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:529
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
FirstRow::field_type field_type
Definition: multitypeblockmatrix.hh:55
void usmv(const AlphaType &alpha, const X &x, Y &y) const
y += alpha A x
Definition: multitypeblockmatrix.hh:225
static void dbgs(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:481
static void bsorf(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:509
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:58
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
static void calc_rhs(const TMatrix &, TVector &, TVector &, Trhs &, const K &)
Definition: multitypeblockmatrix.hh:454
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:444
static void bsorf(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:502
static constexpr size_type size()
Return the number of matrix rows.
Definition: multitypeblockmatrix.hh:64
static void dbjac(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition: multitypeblockmatrix.hh:599
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:53
auto infinity_norm() const
Bastardized version of the infinity-norm / row-sum norm.
Definition: multitypeblockmatrix.hh:364
static void bsorb(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition: multitypeblockmatrix.hh:595
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
auto operator[](const std::integral_constant< size_type, index > indexVariable) -> decltype(std::get< index >(*this))
Random-access operator.
Definition: multitypeblockmatrix.hh:93
static void dbjac(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:564
static void bsorf(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition: multitypeblockmatrix.hh:591
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
static void bsorb(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:536
Definition: allocator.hh:7
std::ostream & operator<<(std::ostream &s, const BlockVector< K, A > &v)
Send BlockVector to an output stream.
Definition: bvector.hh:650
static void bsorb(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:457
static void bsorf(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:415
static void dbjac(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:499
static void dbgs(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:376
A Matrix class to support different block types.
Definition: multitypeblockmatrix.hh:44
solver for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition: multitypeblockmatrix.hh:466
part of solvers for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition: multitypeblockmatrix.hh:438