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...>
55 typedef typename FirstRow::field_type field_type;
60 return 1+
sizeof...(Args);
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))
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))
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>::real_type sum=0;
350 sum += (*this)[i][j].frobenius_norm2();
368 typename FieldTraits<field_type>::real_type norm=0;
373 typename FieldTraits<field_type>::real_type sum=0;
376 sum += (*this)[i][j].infinity_norm();
378 norm =
max(sum, norm);
389 typename FieldTraits<field_type>::real_type norm=0;
394 typename FieldTraits<field_type>::real_type sum=0;
397 sum += (*this)[i][j].infinity_norm_real();
399 norm =
max(sum, norm);
412 template<
typename T1,
typename... Args>
416 using namespace Dune::Hybrid;
418 using namespace Dune::Hybrid;
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&) {}
part of solvers for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition: multitypeblockmatrix.hh:438
solver for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition: multitypeblockmatrix.hh:466
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
#define DUNE_DEPRECATED_MSG(text)
Mark some entity as deprecated.
Definition: deprecated.hh:169
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentionally unused function parameters with.
Definition: unused.hh:25
constexpr void forEach(Range &&range, F &&f)
Range based for loop.
Definition: hybridutilities.hh:267
constexpr auto integralRange(const Begin &begin, const End &end)
Create an integral range.
Definition: hybridutilities.hh:183
void bsorb(const M &A, X &x, const Y &b, const K &w)
SSOR step.
Definition: gsetc.hh:640
void dbjac(const M &A, X &x, const Y &b, const K &w)
Jacobi step.
Definition: gsetc.hh:652
void dbgs(const M &A, X &x, const Y &b, const K &w)
GS step.
Definition: gsetc.hh:616
void bsorf(const M &A, X &x, const Y &b, const K &w)
SOR step.
Definition: gsetc.hh:628
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
MultiTypeBlockMatrix & operator+=(const MultiTypeBlockMatrix &b)
Add the entries of another matrix to this one.
Definition: multitypeblockmatrix.hh:157
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:473
static void bsorb(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:529
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
MultiTypeBlockMatrix & operator/=(const field_type &k)
vector space division by scalar
Definition: multitypeblockmatrix.hh:141
void umhv(const X &x, Y &y) const
y += A^H x
Definition: multitypeblockmatrix.hh:295
FieldTraits< field_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: multitypeblockmatrix.hh:358
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
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
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
MultiTypeBlockMatrix & operator-=(const MultiTypeBlockMatrix &b)
Subtract the entries of another matrix from this one.
Definition: multitypeblockmatrix.hh:172
auto operator[](const std::integral_constant< size_type, index > indexVariable) -> decltype(std::get< index >(*this))
Random-access operator.
Definition: multitypeblockmatrix.hh:93
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
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:14