Do not calculate inverse of T explicitly in IEFSolver

Roberto Di Remigio · Roberto Di Remigio · commit 4b142ea0f5e1 · 2016-07-03T15:59:48.000+02:00
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -7,11 +7,23 @@
 - The `PEDRA.OUT` cavity generator log now reports the initial _and_ final
   lists of spheres. The final list only contains those spheres that were
   actually tesselated and, possibly, the added spheres.
+- For all solvers, the symmetrization needed to obtain the polarization weights
+  happens directly on the computed charges, instead of symmetrizing the system
+  matrices.
+- The `IEFSolver` object now stores the unsymmetrized T and R matrices.
+  The explicit calculation of the inverse of T is thus avoided, the drawback is
+  that two square matrices of size equal to the cavity size are stored instead
+  of just one. A [partially pivoted LU decomposition](http://eigen.tuxfamily.org/dox/classEigen_1_1PartialPivLU.html)
+  is used to solve the linear system.
+  In the isotropic case no inverses are ever explicitly calculated. In the
+  anisotropic case, the inverse of the internal S matrix is still needed to
+  form R.
 - The `CPCMSolver` object now stores the scaled, unsymmetrized S matrix. The
-  explicit calculation and storage of its inverse is thus avoided. A robust
-  Cholesky decomposition is used to solve the linear equation system. The
-  symmetrization needed to obtain the polarization weights happens directly on
-  the computed charges.
+  explicit calculation and storage of its inverse is thus avoided.
+  A [robust Cholesky decomposition](http://eigen.tuxfamily.org/dox/classEigen_1_1LDLT.html) is
+  used to solve the linear equation system.
+  The symmetrization needed to obtain the polarization weights happens directly
+  on the computed charges.
 - The `hermitivitize` function can now correctly symmetrize vectors.
 
 ### Fixed
@@ -25,6 +37,8 @@
 ### Removed
 
 - The function `CPCMMatrix` in the `SolverImpl.hpp` header file is no longer available.
+- The functions `anisotropicIEFMatrix` and `isotropicIEFMatrix` in the
+  `SolverImpl.hpp` header file are no longer available.
 
 ## [v1.1.2] (2016-05-31)
 
diff --git a/src/solver/IEFSolver.cpp b/src/solver/IEFSolver.cpp
@@ -50,36 +50,34 @@ void IEFSolver::buildSystemMatrix_impl(const Cavity & cavity, const IGreensFunct
 
 void IEFSolver::buildAnisotropicMatrix(const Cavity & cav, const IGreensFunction & gf_i, const IGreensFunction & gf_o)
 {
-  fullPCMMatrix_ = anisotropicIEFMatrix(cav, gf_i, gf_o);
-  // Symmetrize K := (K + K+)/2
-  if (hermitivitize_) {
-    hermitivitize(fullPCMMatrix_);
-  }
+  Tepsilon_  = anisotropicTEpsilon(cav, gf_i, gf_o);
   // Pack into a block diagonal matrix
   // The number of irreps in the group
   int nrBlocks = cav.pointGroup().nrIrrep();
   // The size of the irreducible portion of the cavity
   int dimBlock = cav.irreducible_size();
   // For the moment just packs into a std::vector<Eigen::MatrixXd>
-  symmetryPacking(blockPCMMatrix_, fullPCMMatrix_, dimBlock, nrBlocks);
+  symmetryPacking(blockTepsilon_, Tepsilon_, dimBlock, nrBlocks);
+
+  Rinfinity_ = anisotropicRinfinity(cav, gf_i, gf_o);
+  symmetryPacking(blockRinfinity_, Rinfinity_, dimBlock, nrBlocks);
 
   built_ = true;
 }
 
 void IEFSolver::buildIsotropicMatrix(const Cavity & cav, const IGreensFunction & gf_i, const IGreensFunction & gf_o)
 {
-  fullPCMMatrix_ = isotropicIEFMatrix(cav, gf_i, profiles::epsilon(gf_o.permittivity()));
-  // Symmetrize K := (K + K+)/2
-  if (hermitivitize_) {
-    hermitivitize(fullPCMMatrix_);
-  }
+  Tepsilon_  = isotropicTEpsilon(cav, gf_i, profiles::epsilon(gf_o.permittivity()));
   // Pack into a block diagonal matrix
   // The number of irreps in the group
   int nrBlocks = cav.pointGroup().nrIrrep();
   // The size of the irreducible portion of the cavity
   int dimBlock = cav.irreducible_size();
   // For the moment just packs into a std::vector<Eigen::MatrixXd>
-  symmetryPacking(blockPCMMatrix_, fullPCMMatrix_, dimBlock, nrBlocks);
+  symmetryPacking(blockTepsilon_, Tepsilon_, dimBlock, nrBlocks);
+
+  Rinfinity_ = isotropicRinfinity(cav, gf_i);
+  symmetryPacking(blockRinfinity_, Rinfinity_, dimBlock, nrBlocks);
 
   built_ = true;
 }
@@ -89,12 +87,15 @@ Eigen::VectorXd IEFSolver::computeCharge_impl(const Eigen::VectorXd & potential,
   // The potential and charge vector are of dimension equal to the
   // full dimension of the cavity. We have to select just the part
   // relative to the irrep needed.
-  int fullDim = fullPCMMatrix_.rows();
+  int fullDim = Rinfinity_.rows();
   Eigen::VectorXd charge = Eigen::VectorXd::Zero(fullDim);
-  int nrBlocks = blockPCMMatrix_.size();
+  int nrBlocks = blockRinfinity_.size();
   int irrDim = fullDim/nrBlocks;
   charge.segment(irrep*irrDim, irrDim) =
-    - blockPCMMatrix_[irrep] * potential.segment(irrep*irrDim, irrDim);
+    - blockTepsilon_[irrep].partialPivLu().solve(blockRinfinity_[irrep] * potential.segment(irrep*irrDim, irrDim));
+  // Symmetrize ASC charge := (charge + charge*)/2
+  if (hermitivitize_) hermitivitize(charge);
+
   return charge;
 }
 
diff --git a/src/solver/IEFSolver.hpp b/src/solver/IEFSolver.hpp
@@ -2,22 +2,22 @@
 /*
  *     PCMSolver, an API for the Polarizable Continuum Model
  *     Copyright (C) 2013-2015 Roberto Di Remigio, Luca Frediani and contributors
- *     
+ *
  *     This file is part of PCMSolver.
- *     
+ *
  *     PCMSolver is free software: you can redistribute it and/or modify
  *     it under the terms of the GNU Lesser General Public License as published by
  *     the Free Software Foundation, either version 3 of the License, or
  *     (at your option) any later version.
- *     
+ *
  *     PCMSolver 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 Lesser General Public License for more details.
- *     
+ *
  *     You should have received a copy of the GNU Lesser General Public License
  *     along with PCMSolver.  If not, see <http://www.gnu.org/licenses/>.
- *     
+ *
  *     For information on the complete list of contributors to the
  *     PCMSolver API, see: <http://pcmsolver.readthedocs.org/>
  */
@@ -43,7 +43,13 @@ class IGreensFunction;
  *  \class IEFSolver
  *  \brief IEFPCM, collocation-based solver
  *  \author Luca Frediani, Roberto Di Remigio
- *  \date 2011, 2015
+ *  \date 2011, 2015, 2016
+ *
+ *  \note We store the unsymmetrized T(epsilon) and Rinfinity matrices.
+ *  The ASC is obtained by multiplying the MEP by Rinfinity and then using a partially
+ *  pivoted LU decomposition of T(epsilon) on the resulting vector.
+ *  The ASC is then symmetrized. This avoids computing and storing the inverse
+ *  explicitly.
  */
 
 class IEFSolver : public PCMSolver
@@ -73,10 +79,14 @@ class IEFSolver : public PCMSolver
 private:
     /*! Whether the system matrix has to be symmetrized */
     bool hermitivitize_;
-    /*! PCM matrix, not symmetry blocked */
-    Eigen::MatrixXd fullPCMMatrix_;
-    /*! PCM matrix, symmetry blocked form */
-    std::vector<Eigen::MatrixXd> blockPCMMatrix_;
+    /*! T(epsilon) matrix, not symmetry blocked */
+    Eigen::MatrixXd Tepsilon_;
+    /*! T(epsilon) matrix, symmetry blocked form */
+    std::vector<Eigen::MatrixXd> blockTepsilon_;
+    /*! R_infinity matrix, not symmetry blocked */
+    Eigen::MatrixXd Rinfinity_;
+    /*! R_infinity matrix, symmetry blocked form */
+    std::vector<Eigen::MatrixXd> blockRinfinity_;
 
     /*! \brief Calculation of the PCM matrix
      *  \param[in] cavity the cavity to be used
diff --git a/src/solver/SolverImpl.hpp b/src/solver/SolverImpl.hpp
@@ -41,157 +41,6 @@
  *  \date 2015
  */
 
-/*! \brief Builds the **anisotropic** IEFPCM matrix
- *  \param[in] cav the discretized cavity
- *  \param[in] gf_i Green's function inside the cavity
- *  \param[in] gf_o Green's function outside the cavity
- *  \return the \f$ \mathbf{K} = \mathbf{T}^{-1}\mathbf{R}\mathbf{A} \f$ matrix
- *
- *  This function calculates the PCM matrix. We use the following definitions:
- *  \f[
- *     \begin{align}
- *       \mathbf{T} &=
- *      \left(2\pi\mathbf{I} - \mathbf{D}_\mathrm{e}\mathbf{A}\right)\mathbf{S}_\mathrm{i}
- *      +\mathbf{S}_\mathrm{e}\left(2\pi\mathbf{I} +
- *      \mathbf{A}\mathbf{D}_\mathrm{i}^\dagger\right) \\
- *      \mathbf{R} &=
- *      \left(2\pi\mathbf{A}^{-1} - \mathbf{D}_\mathrm{e}\right) -
- *      \mathbf{S}_\mathrm{e}\mathbf{S}^{-1}_\mathrm{i}\left(2\pi\mathbf{A}^{-1}-\mathbf{D}_\mathrm{i}\right)
- *     \end{align}
- *  \f]
- *  The matrix is not symmetrized and is not symmetry packed.
- */
-inline Eigen::MatrixXd anisotropicIEFMatrix(const Cavity & cav, const IGreensFunction & gf_i, const IGreensFunction & gf_o)
-{
-  // The total size of the cavity
-  PCMSolverIndex cavitySize = cav.size();
-  // The number of irreps in the group
-  int nrBlocks = cav.pointGroup().nrIrrep();
-  // The size of the irreducible portion of the cavity
-  int dimBlock = cav.irreducible_size();
-
-  // Compute SI, DI and SE, DE on the whole cavity, regardless of symmetry
-  TIMER_ON("Computing SI");
-  Eigen::MatrixXd SI = gf_i.singleLayer(cav.elements());
-  TIMER_OFF("Computing SI");
-  TIMER_ON("Computing DI");
-  Eigen::MatrixXd DI = gf_i.doubleLayer(cav.elements());
-  TIMER_OFF("Computing DI");
-  TIMER_ON("Computing SE");
-  Eigen::MatrixXd SE = gf_o.singleLayer(cav.elements());
-  TIMER_OFF("Computing SE");
-  TIMER_ON("Computing DE");
-  Eigen::MatrixXd DE = gf_o.doubleLayer(cav.elements());
-  TIMER_OFF("Computing DE");
-
-  // Perform symmetry blocking
-  // If the group is C1 avoid symmetry blocking, we will just pack the fullPCMMatrix
-  // into "block diagonal" when all other manipulations are done.
-  if (cav.pointGroup().nrGenerators() != 0) {
-    TIMER_ON("Symmetry blocking");
-    symmetryBlocking(DI, cavitySize, dimBlock, nrBlocks);
-    symmetryBlocking(SI, cavitySize, dimBlock, nrBlocks);
-    symmetryBlocking(DE, cavitySize, dimBlock, nrBlocks);
-    symmetryBlocking(SE, cavitySize, dimBlock, nrBlocks);
-    TIMER_OFF("Symmetry blocking");
-  }
-
-  Eigen::MatrixXd a = cav.elementArea().asDiagonal();
-  Eigen::MatrixXd Id = Eigen::MatrixXd::Identity(cavitySize, cavitySize);
-
-  // 1. Form T
-  TIMER_ON("Assemble T matrix");
-  Eigen::MatrixXd fullPCMMatrix = ((2 * M_PI * Id - DE * a) * SI + SE * (2 * M_PI * Id + a * DI.adjoint().eval()));
-  TIMER_OFF("Assemble T matrix");
-  // 2. Invert T using LU decomposition with full pivoting
-  //    This is a rank-revealing LU decomposition, this allows us
-  //    to test if T is invertible before attempting to invert it.
-  TIMER_ON("Invert T matrix");
-  Eigen::FullPivLU<Eigen::MatrixXd> T_LU(fullPCMMatrix);
-  if (!(T_LU.isInvertible())) PCMSOLVER_ERROR("T matrix is not invertible!", BOOST_CURRENT_FUNCTION);
-  fullPCMMatrix = T_LU.inverse();
-  TIMER_OFF("Invert T matrix");
-  Eigen::FullPivLU<Eigen::MatrixXd> SI_LU(SI);
-  if (!(SI_LU.isInvertible())) PCMSOLVER_ERROR("SI matrix is not invertible!", BOOST_CURRENT_FUNCTION);
-  TIMER_ON("Assemble T^-1R matrix");
-  fullPCMMatrix *= ((2 * M_PI * Id - DE * a) - SE * SI_LU.inverse() * (2 * M_PI * Id - DI * a));
-  TIMER_OFF("Assemble T^-1R matrix");
-
-  return fullPCMMatrix;
-}
-
-/*! \brief Builds the **isotropic** IEFPCM matrix
- *  \param[in] cav the discretized cavity
- *  \param[in] gf_i Green's function inside the cavity
- *  \param[in] epsilon permittivity outside the cavity
- *  \return the \f$ \mathbf{K} = \mathbf{T}^{-1}\mathbf{R}\mathbf{A} \f$ matrix
- *
- *  This function calculates the PCM matrix. We use the following definitions:
- *  \f[
- *     \begin{align}
- *       \mathbf{T} &=
- *      \left(2\pi\frac{\varepsilon+1}{\varepsilon-1}\mathbf{I} - \mathbf{D}_\mathrm{i}\mathbf{A}\right)\mathbf{S}_\mathrm{i} \\
- *      \mathbf{R} &=
- *      \left(2\pi\mathbf{A}^{-1} - \mathbf{D}_\mathrm{i}\right)
- *     \end{align}
- *  \f]
- *  The matrix is not symmetrized and is not symmetry packed.
- */
-inline Eigen::MatrixXd isotropicIEFMatrix(const Cavity & cav, const IGreensFunction & gf_i, double epsilon)
-{
-  // The total size of the cavity
-  PCMSolverIndex cavitySize = cav.size();
-  // The number of irreps in the group
-  int nrBlocks = cav.pointGroup().nrIrrep();
-  // The size of the irreducible portion of the cavity
-  int dimBlock = cav.irreducible_size();
-
-  // Compute SI and DI on the whole cavity, regardless of symmetry
-  TIMER_ON("Computing SI");
-  Eigen::MatrixXd SI = gf_i.singleLayer(cav.elements());
-  TIMER_OFF("Computing SI");
-  TIMER_ON("Computing DI");
-  Eigen::MatrixXd DI = gf_i.doubleLayer(cav.elements());
-  TIMER_OFF("Computing DI");
-
-  // Perform symmetry blocking
-  // If the group is C1 avoid symmetry blocking, we will just pack the fullPCMMatrix
-  // into "block diagonal" when all other manipulations are done.
-  if (cav.pointGroup().nrGenerators() != 0) {
-    TIMER_ON("Symmetry blocking");
-    symmetryBlocking(DI, cavitySize, dimBlock, nrBlocks);
-    symmetryBlocking(SI, cavitySize, dimBlock, nrBlocks);
-    TIMER_OFF("Symmetry blocking");
-  }
-
-  Eigen::MatrixXd a = cav.elementArea().asDiagonal();
-  Eigen::MatrixXd Id = Eigen::MatrixXd::Identity(cavitySize, cavitySize);
-
-  // Tq = -Rv -> q = -(T^-1 * R)v = -Kv
-  // T = (2 * M_PI * fact * aInv - DI) * a * SI; R = (2 * M_PI * aInv - DI)
-  // fullPCMMatrix_ = K = T^-1 * R * a
-  // 1. Form T
-  double fact = (epsilon + 1.0)/(epsilon - 1.0);
-  TIMER_ON("Assemble T matrix");
-  Eigen::MatrixXd fullPCMMatrix = (2 * M_PI * fact * Id - DI * a) * SI;
-  TIMER_OFF("Assemble T matrix");
-  // 2. Invert T using LU decomposition with full pivoting
-  //    This is a rank-revealing LU decomposition, this allows us
-  //    to test if T is invertible before attempting to invert it.
-  TIMER_ON("Invert T matrix");
-  Eigen::FullPivLU<Eigen::MatrixXd> T_LU(fullPCMMatrix);
-  if (!(T_LU.isInvertible()))
-    PCMSOLVER_ERROR("T matrix is not invertible!", BOOST_CURRENT_FUNCTION);
-  fullPCMMatrix = T_LU.inverse();
-  TIMER_OFF("Invert T matrix");
-  // 3. Multiply T^-1 and R
-  TIMER_ON("Assemble T^-1R matrix");
-  fullPCMMatrix *= (2 * M_PI * Id - DI * a);
-  TIMER_OFF("Assemble T^-1R matrix");
-
-  return fullPCMMatrix;
-}
-
 /*! \brief Builds the **anisotropic** \f$ \mathbf{T}_\varepsilon \f$ matrix
  *  \param[in] cav the discretized cavity
  *  \param[in] gf_i Green's function inside the cavity
@@ -322,9 +171,8 @@ inline Eigen::MatrixXd anisotropicRinfinity(const Cavity & cav, const IGreensFun
   Eigen::MatrixXd a = cav.elementArea().asDiagonal();
   Eigen::MatrixXd Id = Eigen::MatrixXd::Identity(cavitySize, cavitySize);
 
-  // Form T
-  Eigen::FullPivLU<Eigen::MatrixXd> SI_LU(SI);
-  if (!(SI_LU.isInvertible())) PCMSOLVER_ERROR("SI matrix is not invertible!", BOOST_CURRENT_FUNCTION);
+  // Form R
+  Eigen::PartialPivLU<Eigen::MatrixXd> SI_LU(SI);
   return ((2 * M_PI * Id - DE * a) - SE * SI_LU.inverse() * (2 * M_PI * Id - DI * a));
 }
 
