ewaldblockjacobipreconditioner.h

Go to the documentation of this file.

#ifndef VOTCA_XTP_EWALDBLOCKJACOBIPRECONDITIONER_H
#define VOTCA_XTP_EWALDBLOCKJACOBIPRECONDITIONER_H
 
// Standard includes
#include <optional>
#include <stdexcept>
 
// Third party includes
#include <Eigen/Cholesky>
#include <Eigen/IterativeLinearSolvers>
 
// Local VOTCA includes
#include "ewaldrealspaceinteractor.h"
#include "ewaldregistry.h"
 
namespace votca {
namespace xtp {

class EwaldBlockJacobiPreconditioner {
  using Scalar = double;
  using Vector = Eigen::VectorXd;
 
 public:
  using StorageIndex = votca::Index;
  enum {
    ColsAtCompileTime = Eigen::Dynamic,
    MaxColsAtCompileTime = Eigen::Dynamic
  };
 
  EwaldBlockJacobiPreconditioner() : is_initialized_(false) {}
 
  // registry/ids/alpha_ewald/thole_a mirror EwaldPeriodicDipoleOperator's
  // own constructor exactly -- same segments, same layout, same
  // intramolecular interactor parameters -- so a caller already holding
  // those values for the operator can pass the identical ones here.
  EwaldBlockJacobiPreconditioner(EwaldRegistry& registry,
                                 std::vector<Index> ids, double alpha_ewald,
                                 double thole_a)
      : registry_(&registry), ids_(std::move(ids)) {
    intra_interactor_.emplace(alpha_ewald, thole_a);
    BuildOffsets();
    FactorizeBlocks();
    is_initialized_ = true;
  }
 
  Index rows() const { return size_; }
  Index cols() const { return size_; }
 
  // Matches DiagonalPreconditioner's own analyzePattern/factorize/compute
  // split: this class's own local blocks depend only on registry_/ids_/
  // intra_interactor_ (fixed at construction, never on the operator
  // matrix passed in here), so all of these are no-ops that just
  // return *this -- the real work already happened in the constructor.
  // Present only because Eigen::ConjugateGradient's own compute() path
  // requires them to exist on whatever preconditioner type it's given.
  template <typename MatType>
  EwaldBlockJacobiPreconditioner& analyzePattern(const MatType&) {
    return *this;
  }
  template <typename MatType>
  EwaldBlockJacobiPreconditioner& factorize(const MatType&) {
    return *this;
  }
  template <typename MatType>
  EwaldBlockJacobiPreconditioner& compute(const MatType&) {
    return *this;
  }
 
  template <typename Rhs, typename Dest>
  void _solve_impl(const Rhs& b, Dest& x) const {
    // Segments own disjoint slices of x, so they are independent.
#pragma omp parallel for schedule(dynamic, 4)
    for (Index nn = 0; nn < Index(ids_.size()); ++nn) {
      const std::size_t n = std::size_t(nn);
      const Index base = offsets_[n];
      const Index width = offsets_[n + 1] - base;
      // LDLT::solve works identically regardless of block size -- no
      // special case needed for single-site (width==3) segments; their
      // own factorizations_[n] already holds a plain 3x3 LDLT of
      // getPInv() alone (see FactorizeBlocks), the same result
      // DiagonalPreconditioner itself would give them.
      x.segment(base, width) = factorizations_[n].solve(b.segment(base, width));
    }
  }
 
  template <typename Rhs>
  inline const Eigen::Solve<EwaldBlockJacobiPreconditioner, Rhs> solve(
      const Eigen::MatrixBase<Rhs>& b) const {
    eigen_assert(is_initialized_ &&
                 "EwaldBlockJacobiPreconditioner is not initialized.");
    eigen_assert(size_ == b.rows() &&
                 "EwaldBlockJacobiPreconditioner::solve(): invalid number "
                 "of rows of the right hand side matrix b");
    return Eigen::Solve<EwaldBlockJacobiPreconditioner, Rhs>(*this,
                                                             b.derived());
  }
 
  Eigen::ComputationInfo info() const { return Eigen::Success; }
 
 private:
  void BuildOffsets() {
    offsets_.reserve(ids_.size() + 1);
    offsets_.push_back(0);
    for (Index id : ids_) {
      if (!registry_->Has(id, EwaldChargeState::Neutral)) {
        throw std::runtime_error(
            "EwaldBlockJacobiPreconditioner: segment id not registered at "
            "EwaldChargeState::Neutral");
      }
      Index n_sites = registry_->Get(id, EwaldChargeState::Neutral).size();
      offsets_.push_back(offsets_.back() + 3 * n_sites);
    }
    size_ = offsets_.back();
  }
 
  // Builds and factorizes every segment's own local block, matching
  // EwaldPeriodicDipoleOperator's own intramolecular term: the same
  // B-function/ComputeThole calls, and the same sign it enters the
  // OPERATOR with -- i.e. subtracted, since RawMultiply forms
  // A = P^-1 - C. Note this is the opposite sign to
  // AddIntraSegmentCoupling's own internal accumulation, which builds
  // the coupling field itself (a positive quantity) that RawMultiply
  // then subtracts. An earlier version of this comment claimed to
  // mirror that method's sign convention "exactly", and the code did --
  // which is precisely why this was wrong, and why it survived the
  // operator's own sign fix without being noticed. Compare against the
  // operator's assembled block, not against AddIntraSegmentCoupling in
  // isolation.
  void FactorizeBlocks() {
    // One independent block per segment: built and factorized in
    // parallel, each into its own slot.
    factorizations_.resize(ids_.size());
#pragma omp parallel for schedule(dynamic, 4)
    for (Index nn = 0; nn < Index(ids_.size()); ++nn) {
      const std::size_t n = std::size_t(nn);
      const PolarSegment& segment =
          registry_->Get(ids_[n], EwaldChargeState::Neutral);
      const Index n_sites = segment.size();
      const Index width = 3 * n_sites;
      Eigen::MatrixXd block = Eigen::MatrixXd::Zero(width, width);
 
      for (Index s = 0; s < n_sites; ++s) {
        block.block<3, 3>(3 * s, 3 * s) = segment[s].getPInv();
      }
 
      if (n_sites >= 2) {
        for (Index i = 0; i < n_sites; ++i) {
          const PolarSite& site_i = segment[i];
          for (Index j = i + 1; j < n_sites; ++j) {
            const PolarSite& site_j = segment[j];
            const Eigen::Vector3d r_vec = site_i.getPos() - site_j.getPos();
            const double r = r_vec.norm();
            const EwaldRealSpaceInteractor::BFunctions b =
                intra_interactor_->ComputeB(r);
            const EwaldRealSpaceInteractor::TholeFactors t =
                intra_interactor_->ComputeThole(r, site_j, site_i);
            const Eigen::Matrix3d coupling =
                t.l5 * b.B2 * (r_vec * r_vec.transpose()) -
                t.l3 * b.B1 * Eigen::Matrix3d::Identity();
            // BUG FIX (this session): MINUS, not plus. The operator
            // being preconditioned is A = P^-1 - C: RawMultiply builds
            // the intramolecular term with exactly the `coupling`
            // expression above and then subtracts it (result -= intra).
            // This class assembled it with a plus, so it was
            // factorizing P^-1 + C_intra -- the operator as it stood
            // BEFORE the coupling-sign fix, which this file was written
            // against and which never propagated here.
            //
            // The consequence was not a wrong answer (a preconditioner
            // cannot change the fixed point, only the path to it) but a
            // markedly worse one: on a real 5000-site solve this took
            // PCG from 16 iterations to 29, with a visibly
            // non-monotonic p.A.p curvature trace, versus a smoothly
            // decreasing one unpreconditioned.
            block.block<3, 3>(3 * i, 3 * j) -= coupling;
            block.block<3, 3>(3 * j, 3 * i) -= coupling.transpose();
          }
        }
      }
 
      factorizations_[n].compute(block);
    }
  }
 
  EwaldRegistry* registry_ = nullptr;
  std::vector<Index> ids_;
  std::optional<EwaldRealSpaceInteractor> intra_interactor_;
  std::vector<Index> offsets_;
  Index size_ = 0;
  std::vector<Eigen::LDLT<Eigen::MatrixXd>> factorizations_;
  bool is_initialized_;
};
 
}  // namespace xtp
}  // namespace votca
 
#endif  // VOTCA_XTP_EWALDBLOCKJACOBIPRECONDITIONER_H