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>
41 template<
typename FirstRow,
typename... Args>
43 :
public std::tuple<FirstRow, Args...>
55 static constexpr std::size_t
N()
57 return 1+
sizeof...(Args);
61 static constexpr std::size_t
M()
63 return FirstRow::size();
82 template< std::
size_t index >
84 operator[] (
const std::integral_constant< std::size_t, index > indexVariable ) -> decltype(std::get<index>(*
this))
86 DUNE_UNUSED_PARAMETER(indexVariable);
87 return std::get<index>(*this);
95 template< std::
size_t index >
97 operator[] (
const std::integral_constant< std::size_t, index > indexVariable )
const -> decltype(std::get<index>(*
this))
99 DUNE_UNUSED_PARAMETER(indexVariable);
100 return std::get<index>(*this);
108 using namespace Dune::Hybrid;
109 auto size = index_constant<1+
sizeof...(Args)>();
113 forEach(integralRange(size), [&](
auto&& i) {
120 template<
typename X,
typename Y>
121 void mv (
const X& x, Y& y)
const {
122 static_assert(X::size() ==
M(),
"length of x does not match row length");
123 static_assert(Y::size() ==
N(),
"length of y does not match row count");
130 template<
typename X,
typename Y>
131 void umv (
const X& x, Y& y)
const {
132 static_assert(X::size() ==
M(),
"length of x does not match row length");
133 static_assert(Y::size() ==
N(),
"length of y does not match row count");
134 using namespace Dune::Hybrid;
135 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
136 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
137 (*this)[i][j].umv(x[j], y[i]);
144 template<
typename X,
typename Y>
145 void mmv (
const X& x, Y& y)
const {
146 static_assert(X::size() ==
M(),
"length of x does not match row length");
147 static_assert(Y::size() ==
N(),
"length of y does not match row count");
148 using namespace Dune::Hybrid;
149 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
150 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
151 (*this)[i][j].mmv(x[j], y[i]);
158 template<
typename AlphaType,
typename X,
typename Y>
159 void usmv (
const AlphaType& alpha,
const X& x, Y& y)
const {
160 static_assert(X::size() ==
M(),
"length of x does not match row length");
161 static_assert(Y::size() ==
N(),
"length of y does not match row count");
162 using namespace Dune::Hybrid;
163 forEach(integralRange(Hybrid::size(y)), [&](
auto&& i) {
164 forEach(integralRange(Hybrid::size(x)), [&](
auto&& j) {
165 (*this)[i][j].usmv(alpha, x[j], y[i]);
179 template<
typename T1,
typename... Args>
181 auto N = index_constant<MultiTypeBlockMatrix<T1,Args...>::N()>();
182 auto M = index_constant<MultiTypeBlockMatrix<T1,Args...>::M()>();
183 using namespace Dune::Hybrid;
184 forEach(integralRange(
N), [&](
auto&& i) {
185 forEach(integralRange(
M), [&](
auto&& j) {
186 s <<
"\t(" << i <<
", " << j <<
"): \n" << m[i][j];
194 template<
int I,
typename M>
203 template<
int I,
int crow,
int ccol,
int remain_col>
209 template <
typename Trhs,
typename TVector,
typename TMatrix,
typename K>
210 static void calc_rhs(
const TMatrix& A, TVector& x, TVector& v, Trhs& b,
const K& w) {
211 std::get<ccol>( std::get<crow>(A) ).mmv( std::get<ccol>(x), b );
216 template<
int I,
int crow,
int ccol>
219 template <
typename Trhs,
typename TVector,
typename TMatrix,
typename K>
220 static void calc_rhs(
const TMatrix&, TVector&, TVector&, Trhs&,
const K&) {}
231 template<
int I,
int crow,
int remain_row>
238 template <
typename TVector,
typename TMatrix,
typename K>
239 static void dbgs(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
246 template <
typename TVector,
typename TMatrix,
typename K>
247 static void dbgs(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
248 auto rhs = std::get<crow> (b);
253 typename std::remove_cv<
254 typename std::remove_reference<
255 decltype(std::get<crow>( std::get<crow>(A)))
267 template <
typename TVector,
typename TMatrix,
typename K>
268 static void bsorf(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
274 template <
typename TVector,
typename TMatrix,
typename K>
275 static void bsorf(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
276 auto rhs = std::get<crow> (b);
281 typename std::remove_cv<
282 typename std::remove_reference<
283 decltype(std::get<crow>( std::get<crow>(A)))
287 std::get<crow>(x).axpy(w,std::get<crow>(v));
294 template <
typename TVector,
typename TMatrix,
typename K>
295 static void bsorb(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
301 template <
typename TVector,
typename TMatrix,
typename K>
302 static void bsorb(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
303 auto rhs = std::get<crow> (b);
308 typename std::remove_cv<
309 typename std::remove_reference<
310 decltype(std::get<crow>( std::get<crow>(A)))
314 std::get<crow>(x).axpy(w,std::get<crow>(v));
322 template <
typename TVector,
typename TMatrix,
typename K>
323 static void dbjac(
const TMatrix& A, TVector& x,
const TVector& b,
const K& w) {
329 template <
typename TVector,
typename TMatrix,
typename K>
330 static void dbjac(
const TMatrix& A, TVector& x, TVector& v,
const TVector& b,
const K& w) {
331 auto rhs = std::get<crow> (b);
336 typename std::remove_cv<
337 typename std::remove_reference<
338 decltype(std::get<crow>( std::get<crow>(A)))
349 template<
int I,
int crow>
352 template <
typename TVector,
typename TMatrix,
typename K>
353 static void dbgs(
const TMatrix&, TVector&, TVector&,
354 const TVector&,
const K&) {}
356 template <
typename TVector,
typename TMatrix,
typename K>
357 static void bsorf(
const TMatrix&, TVector&, TVector&,
358 const TVector&,
const K&) {}
360 template <
typename TVector,
typename TMatrix,
typename K>
361 static void bsorb(
const TMatrix&, TVector&, TVector&,
362 const TVector&,
const K&) {}
364 template <
typename TVector,
typename TMatrix,
typename K>
365 static void dbjac(
const TMatrix&, TVector&, TVector&,
366 const TVector&,
const K&) {}
void umv(const X &x, Y &y) const
y += A x
Definition: multitypeblockmatrix.hh:131
static void calc_rhs(const TMatrix &, TVector &, TVector &, Trhs &, const K &)
Definition: multitypeblockmatrix.hh:220
A Matrix class to support different block types.
Definition: matrixutils.hh:121
static void bsorb(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:302
static void dbjac(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:330
static void bsorf(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition: multitypeblockmatrix.hh:357
static void dbgs(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:239
static void bsorb(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition: multitypeblockmatrix.hh:361
static void dbjac(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition: multitypeblockmatrix.hh:365
static void dbgs(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:247
static void calc_rhs(const TMatrix &A, TVector &x, TVector &v, Trhs &b, const K &w)
Definition: multitypeblockmatrix.hh:210
auto operator[](const std::integral_constant< std::size_t, index > indexVariable) -> decltype(std::get< index >(*this))
Random-access operator.
Definition: multitypeblockmatrix.hh:84
static void bsorb(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:295
static void bsorf(const TMatrix &A, TVector &x, TVector &v, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:275
FirstRow::field_type field_type
Definition: multitypeblockmatrix.hh:52
static constexpr std::size_t N()
Return the number of matrix rows.
Definition: multitypeblockmatrix.hh:55
part of solvers for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition: multitypeblockmatrix.hh:204
void usmv(const AlphaType &alpha, const X &x, Y &y) const
y += alpha A x
Definition: multitypeblockmatrix.hh:159
static void dbgs(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:373
static void dbjac(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:461
static void bsorf(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:401
void operator=(const T &newval)
Definition: multitypeblockmatrix.hh:107
Definition: basearray.hh:19
void mmv(const X &x, Y &y) const
y -= A x
Definition: multitypeblockmatrix.hh:145
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
static void bsorf(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:268
static void bsorb(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:431
void mv(const X &x, Y &y) const
y = A x
Definition: multitypeblockmatrix.hh:121
std::ostream & operator<<(std::ostream &s, const BlockVector< K, A > &v)
Send BlockVector to an output stream.
Definition: bvector.hh:632
static void dbjac(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:323
static constexpr std::size_t M()
Return the number of matrix columns.
Definition: multitypeblockmatrix.hh:61
MultiTypeBlockMatrix< FirstRow, Args... > type
Definition: multitypeblockmatrix.hh:50
solver for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition: multitypeblockmatrix.hh:21
static void dbgs(const TMatrix &, TVector &, TVector &, const TVector &, const K &)
Definition: multitypeblockmatrix.hh:353