From 80b516d1840488446bf0ceeae97302b8486edaca Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Thu, 12 May 2022 15:40:54 +0200 Subject: [PATCH] Fix two bugs in interleaved CG method --- include/deal.II/lac/solver_cg.h | 58 +++++++++----- include/deal.II/matrix_free/matrix_free.h | 35 +++++++- .../parallel_multigrid_interleave.cc | 2 - ...4est=true.with_lapack=true.mpirun=2.output | 74 ++++++++--------- .../parallel_multigrid_interleave_renumber.cc | 2 - ...est=true.with_lapack=true.mpirun=2.output} | 80 +++++++++---------- 6 files changed, 147 insertions(+), 104 deletions(-) rename tests/matrix_free/{parallel_multigrid_interleave_renumber.wih_p4est=true.with_lapack=true.mpirun=2.output => parallel_multigrid_interleave_renumber.with_p4est=true.with_lapack=true.mpirun=2.output} (83%) diff --git a/include/deal.II/lac/solver_cg.h b/include/deal.II/lac/solver_cg.h index 81f50070aa..667be6679e 100644 --- a/include/deal.II/lac/solver_cg.h +++ b/include/deal.II/lac/solver_cg.h @@ -511,7 +511,6 @@ namespace internal Number beta; double residual_norm; Number previous_alpha; - Number previous_beta; IterationWorkerBase(const MatrixType & A, const PreconditionerType &preconditioner, @@ -535,7 +534,6 @@ namespace internal , beta(Number()) , residual_norm(0.0) , previous_alpha(Number()) - , previous_beta(Number()) {} void @@ -604,8 +602,6 @@ namespace internal const Number previous_r_dot_preconditioner_dot_r = r_dot_preconditioner_dot_r; - this->previous_alpha = alpha; - this->previous_beta = beta; if (std::is_same::value == false) @@ -639,7 +635,9 @@ namespace internal const Number p_dot_A_dot_p = p * v; Assert(std::abs(p_dot_A_dot_p) != 0., ExcDivideByZero()); - alpha = r_dot_preconditioner_dot_r / p_dot_A_dot_p; + + this->previous_alpha = alpha; + alpha = r_dot_preconditioner_dot_r / p_dot_A_dot_p; x.add(alpha, p); residual_norm = std::sqrt(std::abs(r.add_and_dot(-alpha, v, r))); @@ -717,6 +715,9 @@ namespace internal { using Number = typename VectorType::value_type; + Number next_r_dot_preconditioner_dot_r; + Number previous_beta; + IterationWorker(const MatrixType & A, const PreconditionerType &preconditioner, const bool flexible, @@ -728,6 +729,8 @@ namespace internal flexible, memory, x) + , next_r_dot_preconditioner_dot_r(0.) + , previous_beta(0.) {} // This is the main iteration function, that will use some of the @@ -735,6 +738,13 @@ namespace internal void do_iteration(const unsigned int iteration_index) { + if (iteration_index > 1) + { + previous_beta = this->beta; + this->beta = next_r_dot_preconditioner_dot_r / + this->r_dot_preconditioner_dot_r; + } + std::array, 7> vectorized_sums = {}; this->A.vmult( @@ -759,24 +769,23 @@ namespace internal this->r.get_mpi_communicator(), dealii::ArrayView(scalar_sums.data(), 7)); - this->previous_alpha = this->alpha; - this->previous_beta = this->beta; + this->r_dot_preconditioner_dot_r = scalar_sums[6]; const Number p_dot_A_dot_p = scalar_sums[0]; Assert(std::abs(p_dot_A_dot_p) != 0., ExcDivideByZero()); - const Number previous_r_dot_preconditioner_dot_r = scalar_sums[6]; - this->alpha = previous_r_dot_preconditioner_dot_r / p_dot_A_dot_p; - this->residual_norm = std::sqrt( - scalar_sums[3] + - this->alpha * (-2. * scalar_sums[2] + this->alpha * scalar_sums[1])); + this->previous_alpha = this->alpha; + this->alpha = this->r_dot_preconditioner_dot_r / p_dot_A_dot_p; - this->r_dot_preconditioner_dot_r = - previous_r_dot_preconditioner_dot_r + - this->alpha * (-2. * scalar_sums[4] + this->alpha * scalar_sums[5]); + // Round-off errors near zero might yield negative values, so take + // the absolute value in the next two formulas + this->residual_norm = std::sqrt(std::abs( + scalar_sums[3] + + this->alpha * (-2. * scalar_sums[2] + this->alpha * scalar_sums[1]))); - this->beta = this->r_dot_preconditioner_dot_r / - previous_r_dot_preconditioner_dot_r; + next_r_dot_preconditioner_dot_r = std::abs( + this->r_dot_preconditioner_dot_r + + this->alpha * (-2. * scalar_sums[4] + this->alpha * scalar_sums[5])); } // Function that we use if the PreconditionerType implements an apply() @@ -943,16 +952,21 @@ namespace internal typename std::enable_if, U>::type finalize_after_convergence(const unsigned int iteration_index) { - if (iteration_index % 2 == 1) + if (iteration_index % 2 == 1 || iteration_index == 2) this->x.add(this->alpha, this->p); else { using Number = typename VectorType::value_type; const unsigned int end_range = this->x.locally_owned_size(); - Number * x = this->x.begin(); - Number * r = this->r.begin(); - Number * p = this->p.begin(); + Number *x = this->x.begin(); + Number *r = this->r.begin(); + Number *p = this->p.begin(); + + // Note that we use 'beta' here rather than 'previous_beta' in the + // formula above, which is because the shift in beta -> + // previous_beta has not been applied at this stage, allowing us + // to recover the previous search direction const Number alpha_plus_previous_alpha_over_beta = this->alpha + this->previous_alpha / this->previous_beta; const Number previous_alpha_over_beta = @@ -1147,7 +1161,7 @@ namespace internal typename std::enable_if, U>::type finalize_after_convergence(const unsigned int iteration_index) { - if (iteration_index % 2 == 1) + if (iteration_index % 2 == 1 || iteration_index == 2) this->x.add(this->alpha, this->p); else { diff --git a/include/deal.II/matrix_free/matrix_free.h b/include/deal.II/matrix_free/matrix_free.h index 402b52dc81..f19e2c9432 100644 --- a/include/deal.II/matrix_free/matrix_free.h +++ b/include/deal.II/matrix_free/matrix_free.h @@ -1031,7 +1031,7 @@ public: * * @param dof_handler_index_pre_post Since MatrixFree can be initialized * with a vector of DoFHandler objects, each of them will in general have - * vector sizes and thus different ranges returned to + * different vector sizes and thus different ranges returned to * `operation_before_loop` and `operation_after_loop`. Use this variable to * specify which one of the DoFHandler objects the index range should be * associated to. Defaults to the `dof_handler_index` 0. @@ -4621,6 +4621,28 @@ namespace internal + // Apply a unit matrix operation to constrained DoFs: Default cases where we + // cannot detect a LinearAlgebra::distributed::Vector, we do not do + // anything, else we apply the constraints as a unit operation + template + inline void + apply_operation_to_constrained_dofs(const std::vector &, + const VectorStruct1 &, + VectorStruct2 &) + {} + + template + inline void + apply_operation_to_constrained_dofs( + const std::vector & constrained_dofs, + const LinearAlgebra::distributed::Vector &src, + LinearAlgebra::distributed::Vector & dst) + { + for (const unsigned int i : constrained_dofs) + dst.local_element(i) = src.local_element(i); + } + + namespace MatrixFreeFunctions { // struct to select between a const interface and a non-const interface @@ -4861,6 +4883,17 @@ namespace internal { if (operation_after_loop) { + // Run unit matrix operation on constrained dofs if we are at the + // last range + const std::vector &partition_row_index = + matrix_free.get_task_info().partition_row_index; + if (range_index == + partition_row_index[partition_row_index.size() - 2] - 1) + apply_operation_to_constrained_dofs( + matrix_free.get_constrained_dofs(dof_handler_index_pre_post), + src, + dst); + const internal::MatrixFreeFunctions::DoFInfo &dof_info = matrix_free.get_dof_info(dof_handler_index_pre_post); if (range_index == numbers::invalid_unsigned_int) diff --git a/tests/matrix_free/parallel_multigrid_interleave.cc b/tests/matrix_free/parallel_multigrid_interleave.cc index 37fd4f7d75..afebfe7e32 100644 --- a/tests/matrix_free/parallel_multigrid_interleave.cc +++ b/tests/matrix_free/parallel_multigrid_interleave.cc @@ -153,8 +153,6 @@ public: src, operation_before_loop, operation_after_loop); - for (unsigned int i : data.get_constrained_dofs()) - dst.local_element(i) = src.local_element(i); } void diff --git a/tests/matrix_free/parallel_multigrid_interleave.with_p4est=true.with_lapack=true.mpirun=2.output b/tests/matrix_free/parallel_multigrid_interleave.with_p4est=true.with_lapack=true.mpirun=2.output index f33407f3f0..ba30099891 100644 --- a/tests/matrix_free/parallel_multigrid_interleave.with_p4est=true.with_lapack=true.mpirun=2.output +++ b/tests/matrix_free/parallel_multigrid_interleave.with_p4est=true.with_lapack=true.mpirun=2.output @@ -3,63 +3,63 @@ DEAL::Testing FE_Q<2>(1) DEAL::Number of degrees of freedom: 81 DEAL:cg::Starting value 7.00000 DEAL:cg::Convergence step 3 value 1.44113e-07 -DEAL::Number of calls to special vmult for Operator of size 4: 12 -DEAL::Number of calls to special vmult for Operator of size 9: 27 -DEAL::Number of calls to special vmult for Operator of size 25: 27 -DEAL::Number of calls to special vmult for Operator of size 81: 27 +DEAL::Number of calls to special vmult for Operator of size 4: 13 +DEAL::Number of calls to special vmult for Operator of size 9: 28 +DEAL::Number of calls to special vmult for Operator of size 25: 33 +DEAL::Number of calls to special vmult for Operator of size 81: 42 DEAL::Testing FE_Q<2>(1) DEAL::Number of degrees of freedom: 289 DEAL:cg::Starting value 15.0000 DEAL:cg::Convergence step 3 value 1.49408e-06 -DEAL::Number of calls to special vmult for Operator of size 4: 12 -DEAL::Number of calls to special vmult for Operator of size 9: 27 -DEAL::Number of calls to special vmult for Operator of size 25: 27 -DEAL::Number of calls to special vmult for Operator of size 81: 27 -DEAL::Number of calls to special vmult for Operator of size 289: 27 +DEAL::Number of calls to special vmult for Operator of size 4: 13 +DEAL::Number of calls to special vmult for Operator of size 9: 28 +DEAL::Number of calls to special vmult for Operator of size 25: 33 +DEAL::Number of calls to special vmult for Operator of size 81: 42 +DEAL::Number of calls to special vmult for Operator of size 289: 42 DEAL::Testing FE_Q<2>(3) DEAL::Number of degrees of freedom: 625 DEAL:cg::Starting value 23.0000 -DEAL:cg::Convergence step 4 value 6.30996e-08 -DEAL::Number of calls to special vmult for Operator of size 16: 16 -DEAL::Number of calls to special vmult for Operator of size 49: 36 -DEAL::Number of calls to special vmult for Operator of size 169: 36 -DEAL::Number of calls to special vmult for Operator of size 625: 36 +DEAL:cg::Convergence step 4 value 6.60805e-08 +DEAL::Number of calls to special vmult for Operator of size 16: 19 +DEAL::Number of calls to special vmult for Operator of size 49: 48 +DEAL::Number of calls to special vmult for Operator of size 169: 51 +DEAL::Number of calls to special vmult for Operator of size 625: 51 DEAL::Testing FE_Q<2>(3) DEAL::Number of degrees of freedom: 2401 DEAL:cg::Starting value 47.0000 -DEAL:cg::Convergence step 4 value 1.99754e-07 -DEAL::Number of calls to special vmult for Operator of size 16: 16 -DEAL::Number of calls to special vmult for Operator of size 49: 36 -DEAL::Number of calls to special vmult for Operator of size 169: 36 -DEAL::Number of calls to special vmult for Operator of size 625: 36 -DEAL::Number of calls to special vmult for Operator of size 2401: 36 +DEAL:cg::Convergence step 4 value 2.29801e-07 +DEAL::Number of calls to special vmult for Operator of size 16: 19 +DEAL::Number of calls to special vmult for Operator of size 49: 48 +DEAL::Number of calls to special vmult for Operator of size 169: 51 +DEAL::Number of calls to special vmult for Operator of size 625: 51 +DEAL::Number of calls to special vmult for Operator of size 2401: 51 DEAL::Testing FE_Q<3>(1) DEAL::Number of degrees of freedom: 125 DEAL:cg::Starting value 5.19615 DEAL:cg::Convergence step 3 value 5.16655e-09 -DEAL::Number of calls to special vmult for Operator of size 8: 12 -DEAL::Number of calls to special vmult for Operator of size 27: 27 -DEAL::Number of calls to special vmult for Operator of size 125: 27 +DEAL::Number of calls to special vmult for Operator of size 8: 13 +DEAL::Number of calls to special vmult for Operator of size 27: 28 +DEAL::Number of calls to special vmult for Operator of size 125: 34 DEAL::Testing FE_Q<3>(1) DEAL::Number of degrees of freedom: 729 DEAL:cg::Starting value 18.5203 -DEAL:cg::Convergence step 3 value 1.01773e-06 -DEAL::Number of calls to special vmult for Operator of size 8: 12 -DEAL::Number of calls to special vmult for Operator of size 27: 27 -DEAL::Number of calls to special vmult for Operator of size 125: 27 -DEAL::Number of calls to special vmult for Operator of size 729: 27 +DEAL:cg::Convergence step 3 value 1.00122e-06 +DEAL::Number of calls to special vmult for Operator of size 8: 13 +DEAL::Number of calls to special vmult for Operator of size 27: 28 +DEAL::Number of calls to special vmult for Operator of size 125: 34 +DEAL::Number of calls to special vmult for Operator of size 729: 42 DEAL::Testing FE_Q<3>(2) DEAL::Number of degrees of freedom: 729 DEAL:cg::Starting value 18.5203 -DEAL:cg::Convergence step 4 value 5.65213e-09 -DEAL::Number of calls to special vmult for Operator of size 27: 16 -DEAL::Number of calls to special vmult for Operator of size 125: 36 -DEAL::Number of calls to special vmult for Operator of size 729: 36 +DEAL:cg::Convergence step 4 value 5.65178e-09 +DEAL::Number of calls to special vmult for Operator of size 27: 17 +DEAL::Number of calls to special vmult for Operator of size 125: 43 +DEAL::Number of calls to special vmult for Operator of size 729: 51 DEAL::Testing FE_Q<3>(2) DEAL::Number of degrees of freedom: 4913 DEAL:cg::Starting value 58.0948 -DEAL:cg::Convergence step 4 value 3.37746e-08 -DEAL::Number of calls to special vmult for Operator of size 27: 16 -DEAL::Number of calls to special vmult for Operator of size 125: 36 -DEAL::Number of calls to special vmult for Operator of size 729: 36 -DEAL::Number of calls to special vmult for Operator of size 4913: 36 +DEAL:cg::Convergence step 4 value 6.26463e-08 +DEAL::Number of calls to special vmult for Operator of size 27: 17 +DEAL::Number of calls to special vmult for Operator of size 125: 43 +DEAL::Number of calls to special vmult for Operator of size 729: 51 +DEAL::Number of calls to special vmult for Operator of size 4913: 51 diff --git a/tests/matrix_free/parallel_multigrid_interleave_renumber.cc b/tests/matrix_free/parallel_multigrid_interleave_renumber.cc index 430265140d..0a26bbdfe6 100644 --- a/tests/matrix_free/parallel_multigrid_interleave_renumber.cc +++ b/tests/matrix_free/parallel_multigrid_interleave_renumber.cc @@ -194,8 +194,6 @@ public: src, operation_before_loop, operation_after_loop); - for (unsigned int i : data.get_constrained_dofs()) - dst.local_element(i) = src.local_element(i); } void diff --git a/tests/matrix_free/parallel_multigrid_interleave_renumber.wih_p4est=true.with_lapack=true.mpirun=2.output b/tests/matrix_free/parallel_multigrid_interleave_renumber.with_p4est=true.with_lapack=true.mpirun=2.output similarity index 83% rename from tests/matrix_free/parallel_multigrid_interleave_renumber.wih_p4est=true.with_lapack=true.mpirun=2.output rename to tests/matrix_free/parallel_multigrid_interleave_renumber.with_p4est=true.with_lapack=true.mpirun=2.output index f33407f3f0..d51f7fbac3 100644 --- a/tests/matrix_free/parallel_multigrid_interleave_renumber.wih_p4est=true.with_lapack=true.mpirun=2.output +++ b/tests/matrix_free/parallel_multigrid_interleave_renumber.with_p4est=true.with_lapack=true.mpirun=2.output @@ -2,64 +2,64 @@ DEAL::Testing FE_Q<2>(1) DEAL::Number of degrees of freedom: 81 DEAL:cg::Starting value 7.00000 -DEAL:cg::Convergence step 3 value 1.44113e-07 -DEAL::Number of calls to special vmult for Operator of size 4: 12 -DEAL::Number of calls to special vmult for Operator of size 9: 27 -DEAL::Number of calls to special vmult for Operator of size 25: 27 -DEAL::Number of calls to special vmult for Operator of size 81: 27 +DEAL:cg::Convergence step 3 value 2.50859e-07 +DEAL::Number of calls to special vmult for Operator of size 4: 13 +DEAL::Number of calls to special vmult for Operator of size 9: 29 +DEAL::Number of calls to special vmult for Operator of size 25: 33 +DEAL::Number of calls to special vmult for Operator of size 81: 42 DEAL::Testing FE_Q<2>(1) DEAL::Number of degrees of freedom: 289 DEAL:cg::Starting value 15.0000 -DEAL:cg::Convergence step 3 value 1.49408e-06 -DEAL::Number of calls to special vmult for Operator of size 4: 12 -DEAL::Number of calls to special vmult for Operator of size 9: 27 -DEAL::Number of calls to special vmult for Operator of size 25: 27 -DEAL::Number of calls to special vmult for Operator of size 81: 27 -DEAL::Number of calls to special vmult for Operator of size 289: 27 +DEAL:cg::Convergence step 4 value 2.05438e-08 +DEAL::Number of calls to special vmult for Operator of size 4: 17 +DEAL::Number of calls to special vmult for Operator of size 9: 38 +DEAL::Number of calls to special vmult for Operator of size 25: 42 +DEAL::Number of calls to special vmult for Operator of size 81: 51 +DEAL::Number of calls to special vmult for Operator of size 289: 51 DEAL::Testing FE_Q<2>(3) DEAL::Number of degrees of freedom: 625 DEAL:cg::Starting value 23.0000 -DEAL:cg::Convergence step 4 value 6.30996e-08 -DEAL::Number of calls to special vmult for Operator of size 16: 16 -DEAL::Number of calls to special vmult for Operator of size 49: 36 -DEAL::Number of calls to special vmult for Operator of size 169: 36 -DEAL::Number of calls to special vmult for Operator of size 625: 36 +DEAL:cg::Convergence step 4 value 2.11999e-07 +DEAL::Number of calls to special vmult for Operator of size 16: 19 +DEAL::Number of calls to special vmult for Operator of size 49: 48 +DEAL::Number of calls to special vmult for Operator of size 169: 51 +DEAL::Number of calls to special vmult for Operator of size 625: 51 DEAL::Testing FE_Q<2>(3) DEAL::Number of degrees of freedom: 2401 DEAL:cg::Starting value 47.0000 -DEAL:cg::Convergence step 4 value 1.99754e-07 -DEAL::Number of calls to special vmult for Operator of size 16: 16 -DEAL::Number of calls to special vmult for Operator of size 49: 36 -DEAL::Number of calls to special vmult for Operator of size 169: 36 -DEAL::Number of calls to special vmult for Operator of size 625: 36 -DEAL::Number of calls to special vmult for Operator of size 2401: 36 +DEAL:cg::Convergence step 4 value 4.14038e-07 +DEAL::Number of calls to special vmult for Operator of size 16: 19 +DEAL::Number of calls to special vmult for Operator of size 49: 48 +DEAL::Number of calls to special vmult for Operator of size 169: 51 +DEAL::Number of calls to special vmult for Operator of size 625: 51 +DEAL::Number of calls to special vmult for Operator of size 2401: 51 DEAL::Testing FE_Q<3>(1) DEAL::Number of degrees of freedom: 125 DEAL:cg::Starting value 5.19615 -DEAL:cg::Convergence step 3 value 5.16655e-09 -DEAL::Number of calls to special vmult for Operator of size 8: 12 -DEAL::Number of calls to special vmult for Operator of size 27: 27 -DEAL::Number of calls to special vmult for Operator of size 125: 27 +DEAL:cg::Convergence step 3 value 4.30693e-08 +DEAL::Number of calls to special vmult for Operator of size 8: 13 +DEAL::Number of calls to special vmult for Operator of size 27: 28 +DEAL::Number of calls to special vmult for Operator of size 125: 34 DEAL::Testing FE_Q<3>(1) DEAL::Number of degrees of freedom: 729 DEAL:cg::Starting value 18.5203 -DEAL:cg::Convergence step 3 value 1.01773e-06 -DEAL::Number of calls to special vmult for Operator of size 8: 12 -DEAL::Number of calls to special vmult for Operator of size 27: 27 -DEAL::Number of calls to special vmult for Operator of size 125: 27 -DEAL::Number of calls to special vmult for Operator of size 729: 27 +DEAL:cg::Convergence step 3 value 3.62205e-07 +DEAL::Number of calls to special vmult for Operator of size 8: 13 +DEAL::Number of calls to special vmult for Operator of size 27: 28 +DEAL::Number of calls to special vmult for Operator of size 125: 34 +DEAL::Number of calls to special vmult for Operator of size 729: 42 DEAL::Testing FE_Q<3>(2) DEAL::Number of degrees of freedom: 729 DEAL:cg::Starting value 18.5203 -DEAL:cg::Convergence step 4 value 5.65213e-09 -DEAL::Number of calls to special vmult for Operator of size 27: 16 -DEAL::Number of calls to special vmult for Operator of size 125: 36 -DEAL::Number of calls to special vmult for Operator of size 729: 36 +DEAL:cg::Convergence step 4 value 1.13788e-08 +DEAL::Number of calls to special vmult for Operator of size 27: 17 +DEAL::Number of calls to special vmult for Operator of size 125: 43 +DEAL::Number of calls to special vmult for Operator of size 729: 51 DEAL::Testing FE_Q<3>(2) DEAL::Number of degrees of freedom: 4913 DEAL:cg::Starting value 58.0948 -DEAL:cg::Convergence step 4 value 3.37746e-08 -DEAL::Number of calls to special vmult for Operator of size 27: 16 -DEAL::Number of calls to special vmult for Operator of size 125: 36 -DEAL::Number of calls to special vmult for Operator of size 729: 36 -DEAL::Number of calls to special vmult for Operator of size 4913: 36 +DEAL:cg::Convergence step 4 value 9.17189e-08 +DEAL::Number of calls to special vmult for Operator of size 27: 17 +DEAL::Number of calls to special vmult for Operator of size 125: 43 +DEAL::Number of calls to special vmult for Operator of size 729: 51 +DEAL::Number of calls to special vmult for Operator of size 4913: 51 -- 2.39.5