diff --git a/src/solver/solver_cg.cc b/src/solver/solver_cg.cc
index c6885a1..7da51cb 100644
--- a/src/solver/solver_cg.cc
+++ b/src/solver/solver_cg.cc
@@ -1,128 +1,128 @@
/**
* @file solver_cg.cc
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 20 Dec 2017
*
* @brief Implementation of cg solver
*
* Copyright © 2017 Till Junge
*
* µSpectre is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3, or (at
* your option) any later version.
*
* µSpectre is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include "solver/solver_cg.hh"
#include "solver/solver_error.hh"
#include <iomanip>
#include <cmath>
#include <sstream>
namespace muSpectre {
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
SolverCG<DimS, DimM>::SolverCG(Cell_t& cell, Real tol, Uint maxiter,
- bool verbose)
- :Parent(cell, tol, maxiter, verbose),
+ bool verbose, Communicator comm)
+ :Parent(cell, tol, maxiter, verbose), comm{comm},
r_k{make_field<Field_t>("residual r_k", this->collection)},
p_k{make_field<Field_t>("search direction r_k", this->collection)},
Ap_k{make_field<Field_t>("Effect of tangent A*p_k", this->collection)}
{}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
void SolverCG<DimS, DimM>::solve(const Field_t & rhs,
Field_t & x_f) {
- x_f.eigenvec() = this-> solve(rhs.eigenvec(), x_f.eigenvec());
+ x_f.eigenvec() = this->solve(rhs.eigenvec(), x_f.eigenvec());
};
//----------------------------------------------------------------------------//
template <Dim_t DimS, Dim_t DimM>
typename SolverCG<DimS, DimM>::SolvVectorOut
SolverCG<DimS, DimM>::solve(const SolvVectorInC rhs, SolvVectorIn x_0) {
// Following implementation of algorithm 5.2 in Nocedal's Numerical Optimization (p. 112)
auto r = this->r_k.eigen();
auto p = this->p_k.eigen();
auto Ap = this->Ap_k.eigen();
typename Field_t::EigenMap x(x_0.data(), r.rows(), r.cols());
// initialisation of algo
r = this->cell.directional_stiffness_with_copy(x);
r -= typename Field_t::ConstEigenMap(rhs.data(), r.rows(), r.cols());
p = -r;
this->converged = false;
- Real rdr = (r*r).sum();
- Real rhs_norm2 = rhs.squaredNorm();
+ Real rdr = this->comm.sum((r*r).sum());
+ Real rhs_norm2 = this->comm.sum(rhs.squaredNorm());
Real tol2 = ipow(this->tol,2)*rhs_norm2;
size_t count_width{}; // for output formatting in verbose case
if (this->verbose) {
count_width = size_t(std::log10(this->maxiter))+1;
}
for (Uint i = 0;
i < this->maxiter && (rdr > tol2 || i == 0);
++i, ++this->counter) {
Ap = this->cell.directional_stiffness_with_copy(p);
- Real alpha = rdr/(p*Ap).sum();
+ Real alpha = rdr/this->comm.sum((p*Ap).sum());
x += alpha * p;
r += alpha * Ap;
- Real new_rdr = (r*r).sum();
+ Real new_rdr = this->comm.sum((r*r).sum());
Real beta = new_rdr/rdr;
rdr = new_rdr;
if (this->verbose) {
std::cout << " at CG step " << std::setw(count_width) << i
<< ": |r|/|b| = " << std::setw(15) << std::sqrt(rdr/rhs_norm2)
<< ", cg_tol = " << this->tol << std::endl;
}
p = -r+beta*p;
}
if (rdr < tol2) {
this->converged=true;
} else {
std::stringstream err {};
err << " After " << this->counter << " steps, the solver "
<< " FAILED with |r|/|b| = "
<< std::setw(15) << std::sqrt(rdr/rhs_norm2)
<< ", cg_tol = " << this->tol << std::endl;
throw ConvergenceError("Conjugate gradient has not converged." + err.str());
}
return x_0;
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
typename SolverCG<DimS, DimM>::Tg_req_t
SolverCG<DimS, DimM>::get_tangent_req() const {
return tangent_requirement;
}
template class SolverCG<twoD, twoD>;
//template class SolverCG<twoD, threeD>;
template class SolverCG<threeD, threeD>;
} // muSpectre
diff --git a/src/solver/solver_cg.hh b/src/solver/solver_cg.hh
index fffd4d0..6980881 100644
--- a/src/solver/solver_cg.hh
+++ b/src/solver/solver_cg.hh
@@ -1,113 +1,116 @@
/**
* @file solver_cg.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 20 Dec 2017
*
* @brief class for a simple implementation of a conjugate gradient
* solver. This follows algorithm 5.2 in Nocedal's Numerical
* Optimization (p 112)
*
* Copyright © 2017 Till Junge
*
* µSpectre is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3, or (at
* your option) any later version.
*
* µSpectre is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#ifndef SOLVER_CG_H
#define SOLVER_CG_H
#include "solver/solver_base.hh"
+#include "common/communicator.hh"
#include "common/field.hh"
#include <functional>
namespace muSpectre {
/**
* implements the `muSpectre::SolverBase` interface using a
* conjugate gradient solver. This particular class is useful for
* trouble shooting, as it can be made very verbose, but for
* production runs, it is probably better to use
* `muSpectre::SolverCGEigen`.
*/
template <Dim_t DimS, Dim_t DimM=DimS>
class SolverCG: public SolverBase<DimS, DimM>
{
public:
using Parent = SolverBase<DimS, DimM>; //!< base class
//! Input vector for solvers
using SolvVectorIn = typename Parent::SolvVectorIn;
//! Input vector for solvers
using SolvVectorInC = typename Parent::SolvVectorInC;
//! Output vector for solvers
using SolvVectorOut = typename Parent::SolvVectorOut;
using Cell_t = typename Parent::Cell_t; //!< cell type
using Ccoord = typename Parent::Ccoord; //!< cell coordinates type
//! kind of tangent that is required
using Tg_req_t = typename Parent::TangentRequirement;
//! cg only needs to handle fields that look like strain and stress
using Field_t = TensorField<
typename Parent::Collection_t, Real, secondOrder, DimM>;
//! conjugate gradient needs directional stiffness
constexpr static Tg_req_t tangent_requirement{Tg_req_t::NeedEffect};
//! Default constructor
SolverCG() = delete;
//! Constructor with domain resolutions, etc,
- SolverCG(Cell_t& cell, Real tol, Uint maxiter=0, bool verbose =false);
+ SolverCG(Cell_t& cell, Real tol, Uint maxiter=0, bool verbose=false,
+ Communicator comm=Communicator());
//! Copy constructor
SolverCG(const SolverCG &other) = delete;
//! Move constructor
SolverCG(SolverCG &&other) = default;
//! Destructor
virtual ~SolverCG() = default;
//! Copy assignment operator
SolverCG& operator=(const SolverCG &other) = delete;
//! Move assignment operator
SolverCG& operator=(SolverCG &&other) = default;
bool has_converged() const override final {return this->converged;}
//! actual solver
void solve(const Field_t & rhs,
Field_t & x);
// this simplistic implementation has no initialisation phase so the default is ok
SolvVectorOut solve(const SolvVectorInC rhs, SolvVectorIn x_0) override final;
std::string name() const override final {return "CG";}
protected:
+ Communicator comm; //!< communicator
//! returns `muSpectre::Tg_req_t::NeedEffect`
Tg_req_t get_tangent_req() const override final;
Field_t & r_k; //!< residual
Field_t & p_k; //!< search direction
Field_t & Ap_k; //!< effect of tangent on search direction
bool converged{false}; //!< whether the solver has converged
private:
};
} // muSpectre
#endif /* SOLVER_CG_H */