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...>
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>
466 class MultiTypeBlockMatrix_Solver {
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&) {}