From d9c4c1495d7060b8e55f7765b2b74f259976748f Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sun, 18 Sep 2022 13:52:37 +0200 Subject: [PATCH] PreconditionRelaxation/Chebyshev fixes --- include/deal.II/lac/precondition.h | 11 ++-- tests/lac/precondition_chebyshev_07.cc | 32 +++++----- tests/lac/precondition_relaxation_01.cc | 81 ++++++++++++++++--------- 3 files changed, 77 insertions(+), 47 deletions(-) diff --git a/include/deal.II/lac/precondition.h b/include/deal.II/lac/precondition.h index d580e7d4b6..bcf4f473a7 100644 --- a/include/deal.II/lac/precondition.h +++ b/include/deal.II/lac/precondition.h @@ -1145,7 +1145,7 @@ namespace internal Assert(transposed == false, ExcNotImplemented()); - A.vmult(tmp, src); + A.vmult(tmp, dst); preconditioner.vmult( dst, @@ -1237,7 +1237,7 @@ namespace internal A.vmult( tmp, - src, + dst, [&](const unsigned int start_range, const unsigned int end_range) { // zero 'tmp' before running the vmult // operation @@ -1267,7 +1267,10 @@ namespace internal } }); - preconditioner.vmult(dst, tmp); + preconditioner.vmult(dst, tmp, [](const auto, const auto) { + // note: `dst` vector does not have to be zeroed out + // since we add the result into it + }); } } @@ -3162,7 +3165,7 @@ namespace internal DEAL_II_OPENMP_SIMD_PRAGMA for (std::size_t i = start_range; i < end_range; ++i) - tmp_ptr[i] = factor2 * (rhs_ptr[i] - tmp_ptr[i]); + tmp_ptr[i] = rhs_ptr[i] - tmp_ptr[i]; }); // 2) perform vector updates (including preconditioner application) diff --git a/tests/lac/precondition_chebyshev_07.cc b/tests/lac/precondition_chebyshev_07.cc index a7c4044ea6..2e073ae9a1 100644 --- a/tests/lac/precondition_chebyshev_07.cc +++ b/tests/lac/precondition_chebyshev_07.cc @@ -79,7 +79,10 @@ public: if (operation_before_matrix_vector_product) operation_before_matrix_vector_product(0, src.size()); - diagonal_matrix.vmult(dst, src); + if (operation_before_matrix_vector_product) + diagonal_matrix.vmult_add(dst, src); + else + diagonal_matrix.vmult(dst, src); if (operation_after_matrix_vector_product) operation_after_matrix_vector_product(0, src.size()); @@ -211,7 +214,7 @@ main() std::vector> results; { - // Test PreconditionChebyshev + DiagonalMatrix and matrix with + // 1) Test PreconditionChebyshev + DiagonalMatrix and matrix with // pre/post using PreconditionerType = DiagonalMatrix; @@ -234,7 +237,7 @@ main() } { - /// Test PreconditionChebyshev + DiagonalMatrix and matrix without + // 2) Test PreconditionChebyshev + DiagonalMatrix and matrix without // pre/post using PreconditionerType = DiagonalMatrix; @@ -253,7 +256,7 @@ main() } { - // Test PreconditionChebyshev + arbitrary preconditioner and matrix + // 3) Test PreconditionChebyshev + arbitrary preconditioner and matrix // without pre/post using PreconditionerType = MyDiagonalMatrix; @@ -272,7 +275,7 @@ main() } { - // Test PreconditionChebyshev + preconditioner with pre/post and + // 4) Test PreconditionChebyshev + preconditioner with pre/post and // matrix without pre/post using PreconditionerType = MyDiagonalMatrixWithPreAndPost; @@ -291,7 +294,7 @@ main() } { - // Test PreconditionChebyshev + preconditioner with pre/post and + // 5) Test PreconditionChebyshev + preconditioner with pre/post and // matrix with pre/post using PreconditionerType = MyDiagonalMatrixWithPreAndPost; @@ -313,18 +316,15 @@ main() results.emplace_back(test(preconditioner, src)); } - if (std::equal(results.begin(), - results.end(), - results.begin(), - [](const auto &a, const auto &b) { - if (std::abs(std::get<0>(a) - std::get<0>(b)) > 1e-6) - return false; + if (std::all_of(results.begin(), results.end(), [&](const auto &a) { + if (std::abs(std::get<0>(a) - std::get<0>(results[0])) > 1e-6) + return false; - if (std::abs(std::get<1>(a) - std::get<1>(b)) > 1e-6) - return false; + if (std::abs(std::get<1>(a) - std::get<1>(results[0])) > 1e-6) + return false; - return true; - })) + return true; + })) deallog << "OK!" << std::endl; else deallog << "ERROR!" << std::endl; diff --git a/tests/lac/precondition_relaxation_01.cc b/tests/lac/precondition_relaxation_01.cc index 8c5a83c00a..197e345b8d 100644 --- a/tests/lac/precondition_relaxation_01.cc +++ b/tests/lac/precondition_relaxation_01.cc @@ -71,15 +71,20 @@ public: vmult(VectorType & dst, const VectorType &src, const std::function - &operation_before_matrix_vector_product, + &operation_before_matrix_vector_product = {}, const std::function - &operation_after_matrix_vector_product) const + &operation_after_matrix_vector_product = {}) const { - operation_before_matrix_vector_product(0, src.size()); + if (operation_before_matrix_vector_product) + operation_before_matrix_vector_product(0, src.size()); - diagonal_matrix.vmult(dst, src); + if (operation_before_matrix_vector_product) + diagonal_matrix.vmult_add(dst, src); + else + diagonal_matrix.vmult(dst, src); - operation_after_matrix_vector_product(0, src.size()); + if (operation_after_matrix_vector_product) + operation_after_matrix_vector_product(0, src.size()); } VectorType & @@ -215,7 +220,7 @@ main() std::vector> results; { - // Test PreconditionJacobi + // 0) Test PreconditionJacobi PreconditionJacobi preconditioner; PreconditionJacobi::AdditionalData ad; @@ -228,8 +233,8 @@ main() } { - // Test PreconditionRelaxation + DiagonalMatrix: optimized path with - // lambdas + // 1) Test PreconditionRelaxation + DiagonalMatrix and matrix with + // pre/post using PreconditionerType = DiagonalMatrix; using MyMatrixType = MySparseMatrix; @@ -252,7 +257,8 @@ main() } { - // Test PreconditionRelaxation + DiagonalMatrix: optimized path + // 2) Test PreconditionRelaxation + DiagonalMatrix and matrix without + // pre/post using PreconditionerType = DiagonalMatrix; PreconditionRelaxation preconditioner; @@ -270,8 +276,8 @@ main() } { - // Test PreconditionRelaxation + wrapper around DiagonalMatrix: - // optimized path cannot be used + // 3) Test PreconditionRelaxation + arbitrary preconditioner and + // matrix without pre/post using PreconditionerType = MyDiagonalMatrix; PreconditionRelaxation preconditioner; @@ -289,8 +295,8 @@ main() } { - // Test PreconditionRelaxation + wrapper around DiagonalMatrix with - // pre and post: alternative optimized path is taken + // 4) Test PreconditionRelaxation + preconditioner with pre/post and + // matrix without pre/post using PreconditionerType = MyDiagonalMatrixWithPreAndPost; PreconditionRelaxation preconditioner; @@ -307,24 +313,45 @@ main() results.emplace_back(test(preconditioner, src, false)); } - if (std::equal(results.begin(), - results.end(), - results.begin(), - [](const auto &a, const auto &b) { - if (std::abs(std::get<0>(a) - std::get<0>(b)) > 1e-6) - return false; + { + // 5) Test PreconditionRelaxation + preconditioner with pre/post and + // matrix with pre/post + using PreconditionerType = MyDiagonalMatrixWithPreAndPost; + + using MyMatrixType = MySparseMatrix; + + MyMatrixType my_system_matrix(system_matrix); + + PreconditionRelaxation + preconditioner; + + PreconditionRelaxation::AdditionalData ad; + ad.relaxation = relaxation; + ad.n_iterations = n_iterations; + ad.preconditioner = std::make_shared(); + ad.preconditioner->get_vector() = diagonal; + + preconditioner.initialize(my_system_matrix, ad); + + results.emplace_back(test(preconditioner, src, false)); + } + + if (std::all_of(results.begin(), results.end(), [&](const auto &a) { + if (std::abs(std::get<0>(a) - std::get<0>(results[0])) > 1e-6) + return false; - if (std::abs(std::get<1>(a) - std::get<1>(b)) > 1e-6) - return false; + if (std::abs(std::get<1>(a) - std::get<1>(results[0])) > 1e-6) + return false; - if (std::abs(std::get<2>(a) - std::get<2>(b)) > 1e-6) - return false; + if (std::abs(std::get<2>(a) - std::get<2>(results[0])) > 1e-6) + return false; - if (std::abs(std::get<3>(a) - std::get<3>(b)) > 1e-6) - return false; + if (std::abs(std::get<3>(a) - std::get<3>(results[0])) > 1e-6) + return false; - return true; - })) + return true; + })) deallog << "OK!" << std::endl; else deallog << "ERROR!" << std::endl; -- 2.39.5