From d0cd8360f5ddaaf316538f82e093434f6066eb04 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 9 May 2022 10:57:11 +0200 Subject: [PATCH] Make interface to different variants more slim by moving iteration to separate class --- doc/doxygen/references.bib | 12 + include/deal.II/lac/diagonal_matrix.h | 25 +- include/deal.II/lac/solver_cg.h | 758 ++++++++++++++++---------- tests/lac/solver.output | 10 +- tests/lapack/solver_cg.output | 2 +- 5 files changed, 514 insertions(+), 293 deletions(-) diff --git a/doc/doxygen/references.bib b/doc/doxygen/references.bib index 3f9e730389..afa41bf28a 100644 --- a/doc/doxygen/references.bib +++ b/doc/doxygen/references.bib @@ -1787,3 +1787,15 @@ url = {https://doi.org/10.1016/0045-7930(73)90027-3} doi = {10.11588/heidok.00005743}, url = {https://doi.org/10.11588/heidok.00005743} } + +@article{Chronopoulos1989, + author = {Chronopoulos, A. T. and Gear, C. W.}, + title = {S-step Iterative Methods for Symmetric Linear Systems}, + journal = {Journal of Computational and Applied Mathematics}, + volume = {25}, + number = {2}, + year = {1989}, + pages = {153--168}, + doi = {10.1016/0377-0427(89)90045-9}, + url = {https://doi.org/10.1016/0377-0427(89)90045-9} +} diff --git a/include/deal.II/lac/diagonal_matrix.h b/include/deal.II/lac/diagonal_matrix.h index 3050ce3312..46ffbf76f5 100644 --- a/include/deal.II/lac/diagonal_matrix.h +++ b/include/deal.II/lac/diagonal_matrix.h @@ -197,13 +197,16 @@ public: Tvmult_add(VectorType &dst, const VectorType &src) const; /** - * Apply the preconditioner only on a subrange of elements on the vector. + * Apply the preconditioner only to a subrange of elements of the given + * vector. To support this operation, the given `VectorType` template + * argument needs to support a method `begin()` to return the pointer to the + * start of the stored elements. */ void - apply_on_subrange(const std::size_t index_of_first_unknown, - const std::size_t length, - const value_type *src, - value_type * dst) const; + apply_to_subrange(const unsigned int index_of_first_unknown, + const unsigned int length, + const value_type * src, + value_type * dst) const; /** * Initialize vector @p dst to have the same size and partition as @@ -415,11 +418,11 @@ DiagonalMatrix::Tvmult_add(VectorType & dst, template void -DiagonalMatrix::apply_on_subrange( - const std::size_t index_of_first_unknown, - const std::size_t length, - const value_type *src, - value_type * dst) const +DiagonalMatrix::apply_to_subrange( + const unsigned int index_of_first_unknown, + const unsigned int length, + const value_type * src, + value_type * dst) const { AssertIndexRange(index_of_first_unknown, diagonal.locally_owned_elements().n_elements()); @@ -429,7 +432,7 @@ DiagonalMatrix::apply_on_subrange( const value_type *diagonal_entry = diagonal.begin() + index_of_first_unknown; DEAL_II_OPENMP_SIMD_PRAGMA - for (std::size_t i = 0; i < length; ++i) + for (unsigned int i = 0; i < length; ++i) dst[i] = diagonal_entry[i] * src[i]; } diff --git a/include/deal.II/lac/solver_cg.h b/include/deal.II/lac/solver_cg.h index 269d99acd7..f227e18301 100644 --- a/include/deal.II/lac/solver_cg.h +++ b/include/deal.II/lac/solver_cg.h @@ -34,6 +34,14 @@ DEAL_II_NAMESPACE_OPEN // forward declaration #ifndef DOXYGEN class PreconditionIdentity; +namespace LinearAlgebra +{ + namespace distributed + { + template + class Vector; + } +} // namespace LinearAlgebra #endif @@ -93,6 +101,71 @@ class PreconditionIdentity; * The solve() function of this class uses the mechanism described in the * Solver base class to determine convergence. This mechanism can also be used * to observe the progress of the iteration. + * + *

Optimized operations with specific `MatrixType` argument

+ * + * This class enables to embed the vector updates into the matrix-vector + * product in case the `MatrixType` and `PreconditionerType` support such a + * mode of operation. To this end, the `VectorType` needs to be + * LinearAlgebra::distributed::Vector, the class `MatrixType` needs to provide + * a function with the signature + * @code + * void MatrixType::vmult( + * VectorType &, + * const VectorType &, + * const std::function &, + * const std::function &) const + * @endcode + * where the two given functions run before and after the matrix-vector + * product, respectively, and the `PreconditionerType` needs to provide a + * function with the signature + * @code + * void PreconditionerType::apply_to_subrange(unsigned int start_range, + * unsigned int end_range, + * const Number* src_ptr, + * Number* dst_ptr) + * @endcode + + * where the pointers `src_ptr` and `dst_ptr` point to the location in the + * vector where the operation should be applied to. The functions passed to + * `MatrixType::vmult` take as arguments a sub-range among the locally owned + * elements of the vector, defined as half-open intervals. The intervals are + * designed to be scheduled close to the time the matrix-vector product + * touches those entries in the `src` and `dst` vectors, respectively, with + * the requirement that + *
    + *
  • the matrix-vector product may only access an entry in `src` or `dst` + * once the `operation_before_matrix_vector_product` has been run on that + * vector entry;
    • + *
  • `operation_after_matrix_vector_product` may run on a range of entries + * `[i,j)` once the matrix-vector product does not access the entries `[i,j)` + * in `src` and `dst` any more.
    • + *
+ * The motivation for this function is to increase data locality and hence + * cache usage. For the example of a class similar to the one in the step-37 + * tutorial program, the implementation is + * @code + * void + * vmult(LinearAlgebra::distributed::Vector & dst, + * const LinearAlgebra::distributed::Vector &src, + * const std::function + * &operation_before_matrix_vector_product, + * const std::function + * &operation_after_matrix_vector_product) const + * { + * data.cell_loop(&LaplaceOperator::local_apply, + * this, + * dst, + * src, + * operation_before_matrix_vector_product, + * operation_after_matrix_vector_product); + * } + * @endcode + * + * In terms of the SolverCG implementation, the operation before the loop will + * run the updates on the vectors according to a variant presented in + * Algorithm 2.2 of @cite Chronopoulos1989 (but for a preconditioner), whereas + * the operation after the loop performs a total of 7 reductions in parallel. */ template > class SolverCG : public SolverBase @@ -432,8 +505,8 @@ namespace internal template static decltype( std::declval().apply_to_subrange( - std::size_t(), - std::size_t(), + 0U, + 0U, std::declval(), std::declval())) detect_preconditioner(const PreconditionerType2 &); @@ -448,238 +521,417 @@ namespace internal bool>::value && !std::is_same())), - bool>::value; + bool>::value && + std::is_same< + VectorType, + LinearAlgebra::distributed::Vector>::value; }; + + // We need to have a separate declaration for static const members template const bool supports_vmult_with_std_functions::value; - /** - * Internal function to run one iteration of the conjugate gradient solver - * for standard matrix and preconditioner arguments. - */ - template ::value, - int>::type * = nullptr> - void - do_cg_iteration(const MatrixType & A, - const PreconditionerType &preconditioner, - const typename dealii::SolverCG::AdditionalData - & additional_data, - const unsigned int iteration, - VectorType & x, - VectorType & r, - VectorType & p, - VectorType & v, - VectorType & z, - Number & r_dot_preconditioner_dot_r, - Number & alpha, - Number & beta, - double & residual_norm, - const Number /*old_alpha*/, - const Number /*old_beta*/) + typename PreconditionerType> + struct IterationWorkerBase { - const Number old_r_dot_preconditioner_dot_r = r_dot_preconditioner_dot_r; + using Number = typename VectorType::value_type; + + const MatrixType & A; + const PreconditionerType &preconditioner; + const bool flexible; + VectorType & x; + + typename VectorMemory::Pointer r_pointer; + typename VectorMemory::Pointer p_pointer; + typename VectorMemory::Pointer v_pointer; + typename VectorMemory::Pointer z_pointer; + + // Define some aliases for simpler access, using the variables 'r' for + // the residual b - A*x, 'p' for the search direction, and 'v' for the + // auxiliary vector. This naming convention is used e.g. by the + // description on + // https://en.wikipedia.org/wiki/Conjugate_gradient_method. The variable + // 'z' gets only used for the flexible variant of the CG method. + VectorType &r; + VectorType &p; + VectorType &v; + VectorType &z; + + Number r_dot_preconditioner_dot_r; + Number alpha; + Number beta; + double residual_norm; + Number previous_alpha; + Number previous_beta; + + IterationWorkerBase(const MatrixType & A, + const PreconditionerType &preconditioner, + const bool flexible, + VectorMemory &memory, + VectorType & x) + : A(A) + , preconditioner(preconditioner) + , flexible(flexible) + , x(x) + , r_pointer(memory) + , p_pointer(memory) + , v_pointer(memory) + , z_pointer(memory) + , r(*r_pointer) + , p(*p_pointer) + , v(*v_pointer) + , z(*z_pointer) + , r_dot_preconditioner_dot_r(Number()) + , alpha(Number()) + , beta(Number()) + , residual_norm(0.0) + , previous_alpha(Number()) + , previous_beta(Number()) + {} + + void + startup(const VectorType &b) + { + // Initialize without setting the vector entries, as those would soon + // be overwritten anyway + r.reinit(x, true); + p.reinit(x, true); + v.reinit(x, true); + if (flexible) + z.reinit(x, true); + + // compute residual. if vector is zero, then short-circuit the full + // computation + if (!x.all_zero()) + { + A.vmult(r, x); + r.sadd(-1., 1., b); + } + else + r.equ(1., b); - if (std::is_same::value == - false) - { - preconditioner.vmult(v, r); - r_dot_preconditioner_dot_r = r * v; - } - else - r_dot_preconditioner_dot_r = residual_norm * residual_norm; + residual_norm = r.l2_norm(); + } - const VectorType &direction = - std::is_same::value ? r : v; + void + do_iteration(const unsigned int iteration_index) + { + const Number previous_r_dot_preconditioner_dot_r = + r_dot_preconditioner_dot_r; + previous_alpha = alpha; + previous_beta = beta; - if (iteration > 1) - { - Assert(std::abs(old_r_dot_preconditioner_dot_r) != 0., - ExcDivideByZero()); - beta = r_dot_preconditioner_dot_r / old_r_dot_preconditioner_dot_r; - if (additional_data.use_flexible_variant) - beta -= (r * z) / old_r_dot_preconditioner_dot_r; - p.sadd(beta, 1., direction); - } - else - p.equ(1., direction); + if (std::is_same::value == + false) + { + preconditioner.vmult(v, r); + r_dot_preconditioner_dot_r = r * v; + } + else + r_dot_preconditioner_dot_r = residual_norm * residual_norm; + + const VectorType &direction = + std::is_same::value ? r : v; - if (additional_data.use_flexible_variant) - z.swap(v); + if (iteration_index > 1) + { + Assert(std::abs(previous_r_dot_preconditioner_dot_r) != 0., + ExcDivideByZero()); + beta = + r_dot_preconditioner_dot_r / previous_r_dot_preconditioner_dot_r; + if (flexible) + beta -= (r * z) / previous_r_dot_preconditioner_dot_r; + p.sadd(beta, 1., direction); + } + else + p.equ(1., direction); - A.vmult(v, p); + if (flexible) + z.swap(v); - 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; + A.vmult(v, p); - x.add(alpha, p); - residual_norm = std::sqrt(std::abs(r.add_and_dot(-alpha, v, r))); - } + 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; - /** - * Internal function to run one iteration of the conjugate gradient solver - * for matrices and preconditioners that support interleaving the vector - * updates with the matrix-vector product. - */ - template ::value, - int>::type * = nullptr> - void - do_cg_iteration( - const MatrixType & A, - const PreconditionerType &preconditioner, - const typename dealii::SolverCG::AdditionalData &, - const unsigned int iteration, - VectorType & x_vector, - VectorType & r_vector, - VectorType & p_vector, - VectorType & v_vector, - VectorType &, - Number & r_dot_preconditioner_dot_r, - Number & alpha, - Number & beta, - double & residual_norm, - const Number old_alpha, - const Number old_beta) + typename = int> + struct IterationWorker + : public IterationWorkerBase { - const auto operation_before_loop = [&](const unsigned int start_range, - const unsigned int end_range) { - Number * x = x_vector.begin() + start_range; - Number * r = r_vector.begin() + start_range; - Number * p = p_vector.begin() + start_range; - Number * v = v_vector.begin() + start_range; - constexpr unsigned int grain_size = 32; - std::array prec_r; - if (iteration == 1) - { - for (unsigned int j = start_range; j < end_range; j += grain_size) - { - const unsigned int length = std::min(grain_size, end_range - j); - preconditioner.apply_on_subrange(j, length, r, prec_r.data()); - DEAL_II_OPENMP_SIMD_PRAGMA - for (unsigned int i = 0; i < length; ++i) - { - p[i] = prec_r[i]; - v[i] = Number(); - } - p += length; - r += length; - v += length; - } - } - else if (iteration % 2 == 0) - { - for (unsigned int j = start_range; j < end_range; j += grain_size) - { - const unsigned int length = std::min(grain_size, end_range - j); - DEAL_II_OPENMP_SIMD_PRAGMA - for (unsigned int i = 0; i < length; ++i) - r[i] -= alpha * v[i]; - preconditioner.apply_on_subrange(j, length, r, prec_r.data()); - DEAL_II_OPENMP_SIMD_PRAGMA - for (unsigned int i = 0; i < length; ++i) - { - p[i] = beta * p[i] + prec_r[i]; - v[i] = Number(); - } - p += length; - r += length; - v += length; - } - } + IterationWorker(const MatrixType & A, + const PreconditionerType &preconditioner, + const bool flexible, + VectorMemory &memory, + VectorType & x) + : IterationWorkerBase( + A, + preconditioner, + flexible, + memory, + x) + {} + }; + + + + // Internal function to run one iteration of the conjugate gradient solver + // for matrices and preconditioners that support interleaving the vector + // updates with the matrix-vector product. + template + struct IterationWorker< + VectorType, + MatrixType, + PreconditionerType, + typename std::enable_if< + supports_vmult_with_std_functions::value, + int>::type> + : public IterationWorkerBase + { + static constexpr unsigned int grain_size = 32; + + IterationWorker(const MatrixType & A, + const PreconditionerType &preconditioner, + const bool flexible, + VectorMemory &memory, + VectorType & x) + : IterationWorkerBase( + A, + preconditioner, + flexible, + memory, + x) + {} + + void + do_iteration(const unsigned int iteration_index) + { + using Number = typename VectorType::value_type; + + const auto operation_before_loop = [&](const unsigned int start_range, + const unsigned int end_range) { + Number * x = this->x.begin() + start_range; + Number * r = this->r.begin() + start_range; + Number * p = this->p.begin() + start_range; + Number * v = this->v.begin() + start_range; + std::array prec_r; + if (iteration_index == 1) + { + for (unsigned int j = start_range; j < end_range; j += grain_size) + { + const unsigned int length = + std::min(grain_size, end_range - j); + this->preconditioner.apply_to_subrange(j, + length, + r, + prec_r.data()); + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int i = 0; i < length; ++i) + { + p[i] = prec_r[i]; + v[i] = Number(); + } + p += length; + r += length; + v += length; + } + } + else if (iteration_index % 2 == 0) + { + for (unsigned int j = start_range; j < end_range; j += grain_size) + { + const unsigned int length = + std::min(grain_size, end_range - j); + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int i = 0; i < length; ++i) + r[i] -= this->alpha * v[i]; + this->preconditioner.apply_to_subrange(j, + length, + r, + prec_r.data()); + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int i = 0; i < length; ++i) + { + p[i] = this->beta * p[i] + prec_r[i]; + v[i] = Number(); + } + p += length; + r += length; + v += length; + } + } + else + { + const Number alpha_plus_previous_alpha_over_beta = + this->alpha + this->previous_alpha / this->previous_beta; + const Number previous_alpha_over_beta = + this->previous_alpha / this->previous_beta; + for (unsigned int j = start_range; j < end_range; j += grain_size) + { + const unsigned int length = + std::min(grain_size, end_range - j); + this->preconditioner.apply_to_subrange(j, + length, + r, + prec_r.data()); + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int i = 0; i < length; ++i) + { + x[i] += alpha_plus_previous_alpha_over_beta * p[i] - + previous_alpha_over_beta * prec_r[i]; + r[i] -= this->alpha * v[i]; + } + this->preconditioner.apply_to_subrange(j, + length, + r, + prec_r.data()); + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int i = 0; i < length; ++i) + { + p[i] = this->beta * p[i] + prec_r[i]; + v[i] = Number(); + } + p += length; + r += length; + v += length; + x += length; + } + } + }; + + std::array local_sums = {}; + const auto operation_after_loop = [&](const unsigned int start_range, + const unsigned int end_range) { + const Number * x = this->x.begin() + start_range; + const Number * r = this->r.begin() + start_range; + const Number * p = this->p.begin() + start_range; + const Number * v = this->v.begin() + start_range; + std::array prec_r; + std::array prec_v; + for (unsigned int j = start_range; j < end_range; j += grain_size) + { + const unsigned int length = std::min(grain_size, end_range - j); + this->preconditioner.apply_to_subrange(j, + length, + r, + prec_r.data()); + this->preconditioner.apply_to_subrange(j, + length, + v, + prec_v.data()); + for (unsigned int i = 0; i < length; ++i) + { + local_sums[0] += p[i] * v[i]; + local_sums[1] += v[i] * v[i]; + local_sums[2] += r[i] * v[i]; + local_sums[3] += r[i] * r[i]; + local_sums[4] += r[i] * prec_v[i]; + local_sums[5] += v[i] * prec_v[i]; + local_sums[6] += r[i] * prec_r[i]; + } + p += length; + v += length; + r += length; + } + }; + + this->A.vmult(this->v, + this->p, + operation_before_loop, + operation_after_loop); + + Utilities::MPI::sum(dealii::ArrayView(local_sums.data(), + 7), + this->r.get_mpi_communicator(), + dealii::ArrayView(local_sums.data(), 7)); + + this->previous_alpha = this->alpha; + this->previous_beta = this->beta; + + const Number p_dot_A_dot_p = local_sums[0]; + Assert(std::abs(p_dot_A_dot_p) != 0., ExcDivideByZero()); + + const Number previous_r_dot_preconditioner_dot_r = local_sums[6]; + this->alpha = previous_r_dot_preconditioner_dot_r / p_dot_A_dot_p; + this->residual_norm = std::sqrt( + local_sums[3] + + this->alpha * (-2. * local_sums[2] + this->alpha * local_sums[1])); + + this->r_dot_preconditioner_dot_r = + previous_r_dot_preconditioner_dot_r + + this->alpha * (-2. * local_sums[4] + this->alpha * local_sums[5]); + + this->beta = this->r_dot_preconditioner_dot_r / + previous_r_dot_preconditioner_dot_r; + } + + void + finalize_after_convergence(const unsigned int iteration_index) + { + if (iteration_index % 2 == 1) + this->x.add(this->alpha, this->p); else { - const Number alpha_plus_alpha_old = alpha + old_alpha / old_beta; - const Number alpha_old_beta_old = old_alpha / old_beta; - for (unsigned int j = start_range; j < end_range; j += grain_size) + 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(); + const Number alpha_plus_previous_alpha_over_beta = + this->alpha + this->previous_alpha / this->previous_beta; + const Number previous_alpha_over_beta = + this->previous_alpha / this->previous_beta; + + std::array prec_r; + for (unsigned int j = 0; j < end_range; j += grain_size) { const unsigned int length = std::min(grain_size, end_range - j); - preconditioner.apply_on_subrange(j, length, r, prec_r.data()); + this->preconditioner.apply_to_subrange(j, + length, + r, + prec_r.data()); DEAL_II_OPENMP_SIMD_PRAGMA for (unsigned int i = 0; i < length; ++i) - { - x[i] += alpha_plus_alpha_old * p[i] + - alpha_old_beta_old * prec_r[i]; - r[i] -= alpha * v[i]; - } - preconditioner.apply_on_subrange(j, length, r, prec_r.data()); - DEAL_II_OPENMP_SIMD_PRAGMA - for (unsigned int i = 0; i < length; ++i) - { - p[i] = beta * p[i] + prec_r[i]; - v[i] = 0.; - } - p += length; - r += length; - v += length; + x[i] += alpha_plus_previous_alpha_over_beta * p[i] - + previous_alpha_over_beta * prec_r[i]; + x += length; + r += length; + p += length; } } - }; - - std::array local_sums = {}; - const auto operation_after_loop = [&](const unsigned int start_range, - const unsigned int end_range) { - Number * x = x_vector.begin() + start_range; - Number * r = r_vector.begin() + start_range; - Number * p = p_vector.begin() + start_range; - Number * v = v_vector.begin() + start_range; - constexpr unsigned int grain_size = 32; - std::array prec_r; - std::array prec_v; - for (unsigned int j = start_range; j < end_range; j += grain_size) - { - const unsigned int length = std::min(grain_size, end_range - j); - preconditioner.apply_on_subrange(j, length, r, prec_r.data()); - preconditioner.apply_on_subrange(j, length, v, prec_v.data()); - for (unsigned int i = 0; i < length; ++i) - { - local_sums[0] += p[i] * v[i]; - local_sums[1] += v[i] * v[i]; - local_sums[2] += r[i] * v[i]; - local_sums[3] += r[i] * r[i]; - local_sums[4] += r[i] * prec_v[i]; - local_sums[5] += v[i] * prec_v[i]; - local_sums[6] += r[i] * prec_r[i]; - } - } - }; - - A.vmult(v_vector, p_vector, operation_before_loop, operation_after_loop); - - Utilities::MPI::sum(dealii::ArrayView(local_sums.begin(), - 7), - r_vector.get_mpi_communicator(), - dealii::ArrayView(local_sums.begin_raw(), 7)); - - const Number p_dot_A_dot_p = local_sums[0]; - Assert(std::abs(p_dot_A_dot_p) != 0., ExcDivideByZero()); - - const Number old_r_dot_preconditioner_dot_r = local_sums[6]; - alpha = old_r_dot_preconditioner_dot_r / p_dot_A_dot_p; - residual_norm = std::sqrt(local_sums[3] + 2 * local_sums[2] + - alpha * alpha * local_sums[1]); - - r_dot_preconditioner_dot_r = - old_r_dot_preconditioner_dot_r - - alpha * (2. * local_sums[4] - alpha * local_sums[5]); - beta = r_dot_preconditioner_dot_r / old_r_dot_preconditioner_dot_r; - } + } + }; } // namespace SolverCG } // namespace internal @@ -699,22 +951,6 @@ SolverCG::solve(const MatrixType & A, LogStream::Prefix prefix("cg"); - // Memory allocation - typename VectorMemory::Pointer r_pointer(this->memory); - typename VectorMemory::Pointer p_pointer(this->memory); - typename VectorMemory::Pointer v_pointer(this->memory); - typename VectorMemory::Pointer z_pointer(this->memory); - - // Define some aliases for simpler access, using the variables 'r' for the - // residual b - A*x, 'p' for the search direction, and 'v' for the auxiliary - // vector. This naming convention is used e.g. by the description on - // https://en.wikipedia.org/wiki/Conjugate_gradient_method. The variable 'z' - // gets only used for the flexible variant of the CG method. - VectorType &r = *r_pointer; - VectorType &p = *p_pointer; - VectorType &v = *v_pointer; - VectorType &z = *z_pointer; - // Should we build the matrix for eigenvalue computations? const bool do_eigenvalues = !condition_number_signal.empty() || !all_condition_numbers_signal.empty() || @@ -726,33 +962,16 @@ SolverCG::solve(const MatrixType & A, typename VectorType::value_type eigen_beta_alpha = 0; - // resize the vectors, but do not set the values since they'd be overwritten - // soon anyway. - r.reinit(x, true); - p.reinit(x, true); - v.reinit(x, true); - if (determine_beta_by_flexible_formula) - z.reinit(x, true); - - int it = 0; - number r_dot_preconditioner_dot_r = number(); - number beta = number(); - number alpha = number(); - number old_beta = number(); - number old_alpha = number(); - - // compute residual. if vector is zero, then short-circuit the full - // computation - if (!x.all_zero()) - { - A.vmult(r, x); - r.sadd(-1., 1., b); - } - else - r.equ(1., b); + int it = 0; + + internal::SolverCG:: + IterationWorker + worker( + A, preconditioner, determine_beta_by_flexible_formula, this->memory, x); + + worker.startup(b); - double residual_norm = r.l2_norm(); - solver_state = this->iteration_status(0, residual_norm, x); + solver_state = this->iteration_status(0, worker.residual_norm, x); if (solver_state != SolverControl::iterate) return; @@ -760,37 +979,22 @@ SolverCG::solve(const MatrixType & A, { it++; - internal::SolverCG::do_cg_iteration(A, - preconditioner, - additional_data, - it, - x, - r, - p, - v, - z, - r_dot_preconditioner_dot_r, - alpha, - beta, - residual_norm, - old_alpha, - old_beta); - - old_alpha = alpha; - old_beta = beta; - - print_vectors(it, x, r, p); + worker.do_iteration(it); + + print_vectors(it, x, worker.r, worker.p); if (it > 1) { - this->coefficients_signal(old_alpha, beta); - // set up the vectors containing the diagonal and the off diagonal of - // the projected matrix. + this->coefficients_signal(worker.previous_alpha, worker.beta); + // set up the vectors containing the diagonal and the off diagonal + // of the projected matrix. if (do_eigenvalues) { - diagonal.push_back(number(1.) / old_alpha + eigen_beta_alpha); - eigen_beta_alpha = beta / old_alpha; - offdiagonal.push_back(std::sqrt(beta) / old_alpha); + diagonal.push_back(number(1.) / worker.previous_alpha + + eigen_beta_alpha); + eigen_beta_alpha = worker.beta / worker.previous_alpha; + offdiagonal.push_back(std::sqrt(worker.beta) / + worker.previous_alpha); } compute_eigs_and_cond(diagonal, offdiagonal, @@ -798,16 +1002,18 @@ SolverCG::solve(const MatrixType & A, all_condition_numbers_signal); } - solver_state = this->iteration_status(it, residual_norm, x); + solver_state = this->iteration_status(it, worker.residual_norm, x); } + worker.finalize_after_convergence(it); + compute_eigs_and_cond(diagonal, offdiagonal, eigenvalues_signal, condition_number_signal); AssertThrow(solver_state == SolverControl::success, - SolverControl::NoConvergence(it, residual_norm)); + SolverControl::NoConvergence(it, worker.residual_norm)); } diff --git a/tests/lac/solver.output b/tests/lac/solver.output index c40adef8e3..cacd06ee37 100644 --- a/tests/lac/solver.output +++ b/tests/lac/solver.output @@ -63,7 +63,7 @@ DEAL:sor:Richardson::Starting value 3.000 DEAL:sor:Richardson::Convergence step 7 value 0.0004339 DEAL:sor:cg::Starting value 3.000 DEAL:sor:cg::Failure step 100 value 0.2585 -DEAL:sor::Exception: SolverControl::NoConvergence(it, residual_norm) +DEAL:sor::Exception: SolverControl::NoConvergence(it, worker.residual_norm) DEAL:sor:Bicgstab::Starting value 3.000 DEAL:sor:Bicgstab::Convergence step 5 value 4.201e-18 DEAL:sor:GMRES::Starting value 1.462 @@ -78,7 +78,7 @@ DEAL:psor:Richardson::Starting value 3.000 DEAL:psor:Richardson::Convergence step 8 value 0.0004237 DEAL:psor:cg::Starting value 3.000 DEAL:psor:cg::Failure step 100 value 0.1024 -DEAL:psor::Exception: SolverControl::NoConvergence(it, residual_norm) +DEAL:psor::Exception: SolverControl::NoConvergence(it, worker.residual_norm) DEAL:psor:Bicgstab::Starting value 3.000 DEAL:psor:Bicgstab::Convergence step 4 value 0.0007969 DEAL:psor:GMRES::Starting value 1.467 @@ -91,7 +91,7 @@ DEAL::Size 12 Unknowns 121 DEAL::SOR-diff:0.000 DEAL:no-fail:cg::Starting value 11.00 DEAL:no-fail:cg::Failure step 10 value 0.1496 -DEAL:no-fail::Exception: SolverControl::NoConvergence(it, residual_norm) +DEAL:no-fail::Exception: SolverControl::NoConvergence(it, worker.residual_norm) DEAL:no-fail:Bicgstab::Starting value 11.00 DEAL:no-fail:Bicgstab::Failure step 10 value 0.002830 DEAL:no-fail:Bicgstab::Failure step 10 value 0.001961 @@ -161,7 +161,7 @@ DEAL:sor:Richardson::Starting value 11.00 DEAL:sor:Richardson::Convergence step 88 value 0.0009636 DEAL:sor:cg::Starting value 11.00 DEAL:sor:cg::Failure step 100 value 5.643 -DEAL:sor::Exception: SolverControl::NoConvergence(it, residual_norm) +DEAL:sor::Exception: SolverControl::NoConvergence(it, worker.residual_norm) DEAL:sor:Bicgstab::Starting value 11.00 DEAL:sor:Bicgstab::Convergence step 14 value 0.0009987 DEAL:sor:GMRES::Starting value 7.322 @@ -176,7 +176,7 @@ DEAL:psor:Richardson::Starting value 11.00 DEAL:psor:Richardson::Convergence step 89 value 0.0009736 DEAL:psor:cg::Starting value 11.00 DEAL:psor:cg::Failure step 100 value 2.935 -DEAL:psor::Exception: SolverControl::NoConvergence(it, residual_norm) +DEAL:psor::Exception: SolverControl::NoConvergence(it, worker.residual_norm) DEAL:psor:Bicgstab::Starting value 11.00 DEAL:psor:Bicgstab::Convergence step 11 value 0.0005151 DEAL:psor:GMRES::Starting value 7.345 diff --git a/tests/lapack/solver_cg.output b/tests/lapack/solver_cg.output index e7c652f327..e5ae166fea 100644 --- a/tests/lapack/solver_cg.output +++ b/tests/lapack/solver_cg.output @@ -29,7 +29,7 @@ DEAL:no-fail:cg::Condition number estimate: 53.54 DEAL:no-fail:cg::Condition number estimate: 54.75 DEAL:no-fail:cg::Failure step 10 value 0.1496 DEAL:no-fail:cg::Final Eigenvalues: 0.1363 0.6565 1.499 2.357 3.131 3.994 5.392 6.576 7.464 -DEAL:no-fail::Exception: SolverControl::NoConvergence(it, residual_norm) +DEAL:no-fail::Exception: SolverControl::NoConvergence(it, worker.residual_norm) DEAL:no:cg::Starting value 11.00 DEAL:no:cg::Condition number estimate: 11.01 DEAL:no:cg::Condition number estimate: 21.10 -- 2.39.5