* (for the requirements on matrices and vectors in order to work with this
* class, see the documentation of the Solver base class). The type of the
* solution vector must be passed as template argument, and defaults to
- * dealii::Vector<double>.
+ * dealii::Vector<double>. The AdditionalData structure allows to control the
+ * type of residual for the stopping condition.
*
* @note The CG method requires a symmetric preconditioner (i.e., for example,
* SOR is not a possible choice). There is a variant of the solver,
* 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.
+ *
+ * <h3>Preconditioned residual</h3>
+ *
+ * @p AdditionalData allows you to choose between using the explicit
+ * or implicit residual as a stopping condition for the iterative
+ * solver. This behavior can be overridden by using the flag
+ * AdditionalData::use_default_residual. A <tt>true</tt> value refers to the
+ * implicit residual, while <tt>false</tt> reverts
+ * it. The former uses the result of the matrix-vector
+ * product already computed in other algorithm steps to derive the residual by
+ * a mere vector update, whereas the latter explicitly calculates the system
+ * residual with an additional matrix-vector product.
+ * More information on explicit and implicit residual stopping
+ * criteria can be found
+ * <a
+ * href="https://en.wikipedia.org/wiki/Conjugate_gradient_method#Explicit_residual_calculation">link
+ * here</a>.
*/
template <typename VectorType = Vector<double>>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
*/
using size_type = types::global_dof_index;
+
/**
* Standardized data struct to pipe additional data to the solver.
- * Here, it does not store anything but just exists for consistency
- * with the other solver classes.
*/
struct AdditionalData
- {};
+ {
+ /**
+ * Constructor. By default, set the residual of the stopping criterion
+ * to the implicit residual. A <tt>true</tt> value of
+ * AdditionalData::use_default_residual refers to the
+ * implicit residual, while <tt>false</tt> reverts
+ * it. The former uses the result of the matrix-vector
+ * product already computed in other algorithm steps to derive the residual
+ * by a mere vector update, whereas the latter explicitly calculates the
+ * system residual with an additional matrix-vector product. More
+ * information on explicit and implicit residual stopping criteria can be
+ * found <a
+ * href="https://en.wikipedia.org/wiki/Conjugate_gradient_method#Explicit_residual_calculation">link
+ * here</a>.
+ */
+ explicit AdditionalData(const bool use_default_residual = true);
+
+ /**
+ * Flag for the default residual that is used to measure convergence.
+ */
+ bool use_default_residual;
+ };
+
/**
* Constructor.
*/
using size_type = types::global_dof_index;
+
/**
* Standardized data struct to pipe additional data to the solver.
- * Here, it does not store anything but just exists for consistency
- * with the other solver classes.
*/
struct AdditionalData
- {};
+ {
+ /**
+ * Constructor. By default, set the residual of the stopping criterion
+ * to the implicit residual. A <tt>true</tt> value of
+ * AdditionalData::use_default_residual refers to the
+ * implicit residual, while <tt>false</tt> reverts
+ * it. The former uses the result of the matrix-vector
+ * product already computed in other algorithm steps to derive the residual
+ * by a mere vector update, whereas the latter explicitly calculates the
+ * system residual with an additional matrix-vector product. More
+ * information on explicit and implicit residual stopping criteria can be
+ * found <a
+ * href="https://en.wikipedia.org/wiki/Conjugate_gradient_method#Explicit_residual_calculation">link
+ * here</a>.
+ */
+ explicit AdditionalData(const bool use_default_residual = true);
+
+ /**
+ * Flag for the default residual that is used to measure convergence.
+ */
+ bool use_default_residual;
+ };
+
/**
* Constructor.
*/
SolverFlexibleCG(SolverControl &cn,
const AdditionalData &data = AdditionalData());
+
+protected:
+ /**
+ * Additional parameters.
+ */
+ AdditionalData additional_data;
};
#ifndef DOXYGEN
+template <typename VectorType>
+DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
+inline SolverCG<VectorType>::AdditionalData::AdditionalData(
+ const bool use_default_residual)
+ : use_default_residual(use_default_residual)
+{}
template <typename VectorType>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
namespace internal
{
+
namespace SolverCG
{
// This base class is used to select different variants of the conjugate
const PreconditionerType &preconditioner;
const bool flexible;
VectorType &x;
+ const VectorType &b;
typename VectorMemory<VectorType>::Pointer r_pointer;
typename VectorMemory<VectorType>::Pointer p_pointer;
typename VectorMemory<VectorType>::Pointer v_pointer;
typename VectorMemory<VectorType>::Pointer z_pointer;
+ typename VectorMemory<VectorType>::Pointer explicit_r_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
VectorType &p;
VectorType &v;
VectorType &z;
+ VectorType &explicit_r;
+
+ Number r_dot_preconditioner_dot_r;
+ Number alpha;
+ Number beta;
+ double residual_norm;
+ Number previous_alpha;
+ const bool use_default_residual;
- Number r_dot_preconditioner_dot_r;
- Number alpha;
- Number beta;
- double residual_norm;
- Number previous_alpha;
IterationWorkerBase(const MatrixType &A,
const PreconditionerType &preconditioner,
const bool flexible,
VectorMemory<VectorType> &memory,
- VectorType &x)
+ VectorType &x,
+ const VectorType &b,
+ const bool &use_default_residual)
: A(A)
, preconditioner(preconditioner)
, flexible(flexible)
, x(x)
+ , b(b)
, r_pointer(memory)
, p_pointer(memory)
, v_pointer(memory)
, z_pointer(memory)
+ , explicit_r_pointer(memory)
, r(*r_pointer)
, p(*p_pointer)
, v(*v_pointer)
, z(*z_pointer)
+ , explicit_r(*explicit_r_pointer)
, r_dot_preconditioner_dot_r(Number())
, alpha(Number())
, beta(Number())
, residual_norm(0.0)
, previous_alpha(Number())
+ , use_default_residual(use_default_residual)
{}
void
- startup(const VectorType &b)
+ startup()
{
// Initialize without setting the vector entries, as those would soon
// be overwritten anyway
v.reinit(x, true);
if (flexible)
z.reinit(x, true);
+ if (!use_default_residual)
+ explicit_r.reinit(x, true);
// compute residual. if vector is zero, then short-circuit the full
// computation
using BaseClass =
IterationWorkerBase<VectorType, MatrixType, PreconditionerType>;
+
IterationWorker(const MatrixType &A,
const PreconditionerType &preconditioner,
const bool flexible,
VectorMemory<VectorType> &memory,
- VectorType &x)
- : BaseClass(A, preconditioner, flexible, memory, x)
+ VectorType &x,
+ const VectorType &b,
+ const bool &use_default_residual)
+ : BaseClass(A,
+ preconditioner,
+ flexible,
+ memory,
+ x,
+ b,
+ use_default_residual)
{}
using BaseClass::A;
using BaseClass::alpha;
+ using BaseClass::b;
using BaseClass::beta;
+ using BaseClass::explicit_r;
using BaseClass::p;
using BaseClass::preconditioner;
using BaseClass::r;
using BaseClass::r_dot_preconditioner_dot_r;
using BaseClass::residual_norm;
+ using BaseClass::use_default_residual;
using BaseClass::v;
using BaseClass::x;
using BaseClass::z;
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)));
+
+ // compute the residual norm with implicit residual
+ if (use_default_residual)
+ {
+ residual_norm = std::sqrt(std::abs(r.add_and_dot(-alpha, v, r)));
+ }
+ // compute the residual norm with the explicit residual, i.e.
+ // compute l2 norm of Ax - b.
+ else
+ {
+ // compute the residual conjugate gradient update
+ r.add(-alpha, v);
+ // compute explicit residual
+ A.vmult(explicit_r, x);
+ explicit_r.add(-1, b);
+ residual_norm = explicit_r.l2_norm();
+ }
}
void
const PreconditionerType &preconditioner,
const bool flexible,
VectorMemory<VectorType> &memory,
- VectorType &x)
+ VectorType &x,
+ const VectorType &b,
+ const bool &use_default_residual)
: IterationWorkerBase<VectorType, MatrixType, PreconditionerType>(
A,
preconditioner,
flexible,
memory,
- x)
+ x,
+ b,
+ use_default_residual)
, next_r_dot_preconditioner_dot_r(0.)
, previous_beta(0.)
{}
internal::SolverCG::
IterationWorker<VectorType, MatrixType, PreconditionerType>
- worker(
- A, preconditioner, determine_beta_by_flexible_formula, this->memory, x);
+ worker(A,
+ preconditioner,
+ determine_beta_by_flexible_formula,
+ this->memory,
+ x,
+ b,
+ additional_data.use_default_residual);
- worker.startup(b);
+ worker.startup();
solver_state = this->iteration_status(0, worker.residual_norm, x);
if (solver_state != SolverControl::iterate)
+template <typename VectorType>
+DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
+SolverFlexibleCG<VectorType>::AdditionalData::AdditionalData(
+ const bool use_default_residual)
+ : use_default_residual(use_default_residual)
+{}
+
+
template <typename VectorType>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
SolverFlexibleCG<VectorType>::SolverFlexibleCG(SolverControl &cn,
VectorMemory<VectorType> &mem,
- const AdditionalData &)
+ const AdditionalData &data)
: SolverCG<VectorType>(cn, mem)
{
this->determine_beta_by_flexible_formula = true;
+ this->additional_data = data;
}
template <typename VectorType>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
-SolverFlexibleCG<VectorType>::SolverFlexibleCG(SolverControl &cn,
- const AdditionalData &)
+SolverFlexibleCG<VectorType>::SolverFlexibleCG(SolverControl &cn,
+ const AdditionalData &data)
: SolverCG<VectorType>(cn)
{
this->determine_beta_by_flexible_formula = true;
+ this->additional_data = data;
}
-
#endif // DOXYGEN
DEAL_II_NAMESPACE_CLOSE
SolverControl control(100, 1.e-3);
SolverControl verbose_control(100, 1.e-3, true);
SolverCG<> cg(control, mem);
+ SolverCG<>::AdditionalData data0(false);
+ SolverCG<> cg_add_data(control, mem, data0);
SolverGMRES<>::AdditionalData data1(6);
SolverGMRES<> gmres(control, mem, data1);
SolverGMRES<>::AdditionalData data2(6, true);
control.set_max_steps(10);
check_solve(cg, A, u, f, prec_no);
+ check_solve(cg_add_data, A, u, f, prec_no);
check_solve(bicgstab, A, u, f, prec_no);
check_solve(gmres, A, u, f, prec_no);
check_solve(gmresright, A, u, f, prec_no);
rich.set_omega(1. / A.diag_element(0));
check_solve(rich, A, u, f, prec_no);
check_solve(cg, A, u, f, prec_no);
+ check_solve(cg_add_data, A, u, f, prec_no);
check_solve(bicgstab, A, u, f, prec_no);
check_solve(gmres, A, u, f, prec_no);
check_solve(gmresright, A, u, f, prec_no);
rich.set_omega(1. / A.diag_element(0));
check_solve(rich, A, u, f, prec_richardson);
check_solve(cg, A, u, f, prec_richardson);
+ check_solve(cg_add_data, A, u, f, prec_richardson);
check_solve(bicgstab, A, u, f, prec_richardson);
check_solve(gmres, A, u, f, prec_richardson);
check_solve(gmresright, A, u, f, prec_richardson);
check_Tsolve(rich, A, u, f, prec_ssor);
check_solve(rich, A, u, f, prec_ssor);
check_solve(cg, A, u, f, prec_ssor);
+ check_solve(cg_add_data, A, u, f, prec_ssor);
check_solve(bicgstab, A, u, f, prec_ssor);
check_solve(gmres, A, u, f, prec_ssor);
check_solve(gmresright, A, u, f, prec_ssor);
check_Tsolve(rich, A, u, f, prec_sor);
check_solve(rich, A, u, f, prec_sor);
check_solve(cg, A, u, f, prec_sor);
+ check_solve(cg_add_data, A, u, f, prec_sor);
check_solve(bicgstab, A, u, f, prec_sor);
check_solve(gmres, A, u, f, prec_sor);
check_solve(gmresright, A, u, f, prec_sor);
check_Tsolve(rich, A, u, f, prec_psor);
check_solve(rich, A, u, f, prec_psor);
check_solve(cg, A, u, f, prec_psor);
+ check_solve(cg_add_data, A, u, f, prec_psor);
check_solve(bicgstab, A, u, f, prec_psor);
check_solve(gmres, A, u, f, prec_psor);
check_solve(gmresright, A, u, f, prec_psor);
DEAL::SOR-diff:0.000
DEAL:no-fail:cg::Starting value 3.000
DEAL:no-fail:cg::Convergence step 3 value 4.807e-17
+DEAL:no-fail:cg::Starting value 3.000
+DEAL:no-fail:cg::Convergence step 3 value 9.679e-16
DEAL:no-fail:Bicgstab::Starting value 3.000
DEAL:no-fail:Bicgstab::Convergence step 3 value 1.104e-17
DEAL:no-fail:GMRES::Starting value 3.000
-DEAL:no-fail:GMRES::Convergence step 3 value 4.551e-16
+DEAL:no-fail:GMRES::Convergence step 3 value 3.512e-09
DEAL:no-fail:GMRES::Starting value 3.000
-DEAL:no-fail:GMRES::Convergence step 3 value 4.551e-16
+DEAL:no-fail:GMRES::Convergence step 3 value 3.512e-09
DEAL:no-fail:GMRES::Starting value 3.000
DEAL:no-fail:GMRES::Convergence step 3 value 4.414e-16
DEAL:no-fail:SQMR::Starting value 3.000
DEAL:no:Richardson::Convergence step 24 value 0.0007118
DEAL:no:cg::Starting value 3.000
DEAL:no:cg::Convergence step 3 value 4.807e-17
+DEAL:no:cg::Starting value 3.000
+DEAL:no:cg::Convergence step 3 value 9.679e-16
DEAL:no:Bicgstab::Starting value 3.000
DEAL:no:Bicgstab::Convergence step 3 value 1.104e-17
DEAL:no:GMRES::Starting value 3.000
-DEAL:no:GMRES::Convergence step 3 value 4.551e-16
+DEAL:no:GMRES::Convergence step 3 value 3.512e-09
DEAL:no:GMRES::Starting value 3.000
-DEAL:no:GMRES::Convergence step 3 value 4.551e-16
+DEAL:no:GMRES::Convergence step 3 value 3.512e-09
DEAL:no:GMRES::Starting value 3.000
DEAL:no:GMRES::Convergence step 3 value 4.414e-16
DEAL:no:SQMR::Starting value 3.000
DEAL:rich:Richardson::Convergence step 42 value 0.0008696
DEAL:rich:cg::Starting value 3.000
DEAL:rich:cg::Convergence step 3 value 2.871e-16
+DEAL:rich:cg::Starting value 3.000
+DEAL:rich:cg::Convergence step 3 value 1.691e-15
DEAL:rich:Bicgstab::Starting value 3.000
DEAL:rich:Bicgstab::Convergence step 3 value 3.238e-18
DEAL:rich:GMRES::Starting value 1.800
-DEAL:rich:GMRES::Convergence step 3 value 3.294e-16
+DEAL:rich:GMRES::Convergence step 3 value 3.512e-09
DEAL:rich:GMRES::Starting value 3.000
-DEAL:rich:GMRES::Convergence step 3 value 7.148e-16
+DEAL:rich:GMRES::Convergence step 3 value 4.139e-09
DEAL:rich:GMRES::Starting value 1.800
-DEAL:rich:GMRES::Convergence step 3 value 1.934e-16
+DEAL:rich:GMRES::Convergence step 3 value 1.976e-16
DEAL:rich:SQMR::Starting value 3.000
DEAL:rich:SQMR::Convergence step 3 value 1.422e-15
DEAL:rich:FIRE::Starting value 9.000
DEAL:ssor:Richardson::Convergence step 7 value 0.0006128
DEAL:ssor:cg::Starting value 3.000
DEAL:ssor:cg::Convergence step 4 value 6.018e-05
+DEAL:ssor:cg::Starting value 3.000
+DEAL:ssor:cg::Convergence step 4 value 6.018e-05
DEAL:ssor:Bicgstab::Starting value 3.000
DEAL:ssor:Bicgstab::Convergence step 2 value 0.0003670
DEAL:ssor:GMRES::Starting value 1.920
DEAL:sor:cg::Starting value 3.000
DEAL:sor:cg::Failure step 100 value 0.2585
DEAL:sor::Exception: SolverControl::NoConvergence(it, worker.residual_norm)
+DEAL:sor:cg::Starting value 3.000
+DEAL:sor:cg::Failure step 100 value 0.2585
+DEAL:sor::Exception: SolverControl::NoConvergence(it, worker.residual_norm)
DEAL:sor:Bicgstab::Starting value 3.000
DEAL:sor:Bicgstab::Convergence step 5 value 4.331e-18
DEAL:sor:GMRES::Starting value 1.462
-DEAL:sor:GMRES::Convergence step 5 value 9.441e-19
+DEAL:sor:GMRES::Convergence step 5 value 3.488e-10
DEAL:sor:GMRES::Starting value 3.000
-DEAL:sor:GMRES::Convergence step 5 value 1.393e-17
+DEAL:sor:GMRES::Convergence step 5 value 1.520e-09
DEAL:sor:GMRES::Starting value 1.462
DEAL:sor:GMRES::Convergence step 21 value 0.0007019
DEAL:sor:FIRE::Starting value 9.000
DEAL:psor:cg::Starting value 3.000
DEAL:psor:cg::Failure step 100 value 0.1024
DEAL:psor::Exception: SolverControl::NoConvergence(it, worker.residual_norm)
+DEAL:psor:cg::Starting value 3.000
+DEAL:psor:cg::Failure step 100 value 0.1024
+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
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, worker.residual_norm)
+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, 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
DEAL:no::Exception: SolverControl::NoConvergence(iter, last_criterion)
DEAL:no:cg::Starting value 11.00
DEAL:no:cg::Convergence step 15 value 0.0005794
+DEAL:no:cg::Starting value 11.00
+DEAL:no:cg::Convergence step 15 value 0.0005794
DEAL:no:Bicgstab::Starting value 11.00
DEAL:no:Bicgstab::Convergence step 11 value 0.0006357
DEAL:no:GMRES::Starting value 11.00
DEAL:rich::Exception: SolverControl::NoConvergence(iter, last_criterion)
DEAL:rich:cg::Starting value 11.00
DEAL:rich:cg::Convergence step 15 value 0.0005794
+DEAL:rich:cg::Starting value 11.00
+DEAL:rich:cg::Convergence step 15 value 0.0005794
DEAL:rich:Bicgstab::Starting value 11.00
DEAL:rich:Bicgstab::Convergence step 11 value 0.0006357
DEAL:rich:GMRES::Starting value 6.600
DEAL:ssor:Richardson::Convergence step 48 value 0.0009502
DEAL:ssor:cg::Starting value 11.00
DEAL:ssor:cg::Convergence step 8 value 0.0005816
+DEAL:ssor:cg::Starting value 11.00
+DEAL:ssor:cg::Convergence step 8 value 0.0005816
DEAL:ssor:Bicgstab::Starting value 11.00
DEAL:ssor:Bicgstab::Convergence step 5 value 0.0002498
DEAL:ssor:GMRES::Starting value 13.27
DEAL:sor:cg::Starting value 11.00
DEAL:sor:cg::Failure step 100 value 5.643
DEAL:sor::Exception: SolverControl::NoConvergence(it, worker.residual_norm)
+DEAL:sor:cg::Starting value 11.00
+DEAL:sor:cg::Failure step 100 value 5.643
+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
DEAL:psor:cg::Starting value 11.00
DEAL:psor:cg::Failure step 100 value 2.935
DEAL:psor::Exception: SolverControl::NoConvergence(it, worker.residual_norm)
+DEAL:psor:cg::Starting value 11.00
+DEAL:psor:cg::Failure step 100 value 2.935
+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
projects / dealii.git / commitdiff