From a67d62d9ada0fb30b1a82e2c28d807523ee27a69 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Tue, 30 Aug 2022 09:33:09 +0200 Subject: [PATCH] Enable SolverIDR for block vectors --- include/deal.II/lac/solver_idr.h | 57 ++++++++++++++++++++++++++++++-- 1 file changed, 55 insertions(+), 2 deletions(-) diff --git a/include/deal.II/lac/solver_idr.h b/include/deal.II/lac/solver_idr.h index 64bc3c6793..3c6c66daba 100644 --- a/include/deal.II/lac/solver_idr.h +++ b/include/deal.II/lac/solver_idr.h @@ -24,6 +24,7 @@ #include #include +#include #include #include #include @@ -226,6 +227,53 @@ namespace internal } return *data[i]; } + + + + template ::value, VectorType> + * = nullptr> + unsigned int + n_blocks(const VectorType &) + { + return 1; + } + + + + template ::value, VectorType> * = + nullptr> + unsigned int + n_blocks(const VectorType &vector) + { + return vector.n_blocks(); + } + + + + template ::value, VectorType> + * = nullptr> + VectorType & + block(VectorType &vector, const unsigned int b) + { + AssertDimension(b, 0); + (void)b; + return vector; + } + + + + template ::value, VectorType> * = + nullptr> + typename VectorType::BlockType & + block(VectorType &vector, const unsigned int b) + { + return vector.block(b); + } + } // namespace SolverIDRImplementation } // namespace internal @@ -328,8 +376,13 @@ SolverIDR::solve(const MatrixType & A, VectorType &tmp_q = Q(i, x); if (i != 0) { - for (auto indx : tmp_q.locally_owned_elements()) - tmp_q(indx) = normal_distribution(rng); + for (unsigned int b = 0; + b < internal::SolverIDRImplementation::n_blocks(tmp_q); + ++b) + for (auto indx : internal::SolverIDRImplementation::block(tmp_q, b) + .locally_owned_elements()) + internal::SolverIDRImplementation::block(tmp_q, b)(indx) = + normal_distribution(rng); tmp_q.compress(VectorOperation::insert); } else -- 2.39.5