From: Peter Munch Date: Wed, 7 Sep 2022 10:12:13 +0000 (+0200) Subject: PreconditionRelaxation: add specialization X-Git-Tag: v9.5.0-rc1~985^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=2330bb4e04a4ed80e27e96b9e54925e7b4b2d883;p=dealii.git PreconditionRelaxation: add specialization --- diff --git a/include/deal.II/lac/precondition.h b/include/deal.II/lac/precondition.h index 2972ddc408..04779e88cc 100644 --- a/include/deal.II/lac/precondition.h +++ b/include/deal.II/lac/precondition.h @@ -1086,13 +1086,16 @@ namespace internal // 1) specialized implementation with a preconditioner that accepts // ranges - template , - int> * = nullptr> + template < + typename MatrixType, + typename PreconditionerType, + typename VectorType, + std::enable_if_t< + has_vmult_with_std_functions_for_precondition && + !has_vmult_with_std_functions_for_precondition, + int> * = nullptr> void step_operations(const MatrixType & A, const PreconditionerType &preconditioner, @@ -1172,7 +1175,100 @@ namespace internal } } - // 2) specialized implementation for inverse-diagonal preconditioner + // 2) specialized implementation with a preconditioner and a matrix that + // accepts ranges + template < + typename MatrixType, + typename PreconditionerType, + typename VectorType, + std::enable_if_t< + has_vmult_with_std_functions_for_precondition && + has_vmult_with_std_functions_for_precondition, + int> * = nullptr> + void + step_operations(const MatrixType & A, + const PreconditionerType &preconditioner, + VectorType & dst, + const VectorType & src, + const double relaxation, + VectorType & tmp, + VectorType &, + const unsigned int i, + const bool transposed) + { + using Number = typename VectorType::value_type; + + if (i == 0) + { + Number * dst_ptr = dst.begin(); + const Number *src_ptr = src.begin(); + + preconditioner.vmult( + dst, + src, + [&](const unsigned int start_range, const unsigned int end_range) { + // zero 'dst' before running the vmult operation + if (end_range > start_range) + std::memset(dst.begin() + start_range, + 0, + sizeof(Number) * (end_range - start_range)); + }, + [&](const unsigned int start_range, const unsigned int end_range) { + if (relaxation == 1.0) + return; // nothing to do + + const auto src_ptr = src.begin(); + const auto dst_ptr = dst.begin(); + + DEAL_II_OPENMP_SIMD_PRAGMA + for (std::size_t i = start_range; i < end_range; ++i) + dst_ptr[i] *= relaxation; + }); + } + else + { + tmp.reinit(src, true); + + Assert(transposed == false, ExcNotImplemented()); + + A.vmult( + tmp, + src, + [&](const unsigned int start_range, const unsigned int end_range) { + // zero 'tmp' before running the vmult + // operation + if (end_range > start_range) + std::memset(tmp.begin() + start_range, + 0, + sizeof(Number) * (end_range - start_range)); + }, + [&](const unsigned int start_range, const unsigned int end_range) { + const auto src_ptr = src.begin(); + const auto tmp_ptr = tmp.begin(); + + if (relaxation == 1.0) + { + DEAL_II_OPENMP_SIMD_PRAGMA + for (std::size_t i = start_range; i < end_range; ++i) + tmp_ptr[i] = src_ptr[i] - tmp_ptr[i]; + } + else + { + // note: we scale the residual here to be able to add into + // the dst vector, which contains the solution from the last + // iteration + DEAL_II_OPENMP_SIMD_PRAGMA + for (std::size_t i = start_range; i < end_range; ++i) + tmp_ptr[i] = relaxation * (src_ptr[i] - tmp_ptr[i]); + } + }); + + preconditioner.vmult(dst, tmp); + } + } + + // 3) specialized implementation for inverse-diagonal preconditioner template