Revert to storing T and R explicitly instead of computing T^-1R

Author: Roberto Di Remigio

src/solver/IEFSolver.cpp

@@ -50,28 +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)
 {
-  TinvR_  = anisotropicIEFMatrix(cav, gf_i, gf_o);
+  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(blockTinvR_, TinvR_, 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)
 {
-  TinvR_  = isotropicIEFMatrix(cav, gf_i, profiles::epsilon(gf_o.permittivity()));
+  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(blockTinvR_, TinvR_, dimBlock, nrBlocks);
+  symmetryPacking(blockTepsilon_, Tepsilon_, dimBlock, nrBlocks);
+
+  Rinfinity_ = isotropicRinfinity(cav, gf_i);
+  symmetryPacking(blockRinfinity_, Rinfinity_, dimBlock, nrBlocks);
 
   built_ = true;
 }
@@ -81,12 +87,12 @@ 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 = TinvR_.rows();
+  int fullDim = Rinfinity_.rows();
   Eigen::VectorXd charge = Eigen::VectorXd::Zero(fullDim);
-  int nrBlocks = blockTinvR_.size();
+  int nrBlocks = blockRinfinity_.size();
   int irrDim = fullDim/nrBlocks;
   charge.segment(irrep*irrDim, irrDim) =
-    - blockTinvR_[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);
 

src/solver/IEFSolver.hpp

@@ -45,7 +45,10 @@ class IGreensFunction;
  *  \author Luca Frediani, Roberto Di Remigio
  *  \date 2011, 2015, 2016
  *
- *  \note We store the unsymmetrized T^-1R matrix.
+ *  \note We store the non-Hermitian, symmetry-blocked 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.
+ *  This avoids computing and storing the inverse explicitly.
  */
 
 class IEFSolver : public PCMSolver
@@ -75,10 +78,14 @@ class IEFSolver : public PCMSolver
 private:
     /*! Whether the system matrix has to be symmetrized */
     bool hermitivitize_;
-    /*! T^-1R matrix, not symmetry blocked */
-    Eigen::MatrixXd TinvR_;
-    /*! T^-1R matrix, symmetry blocked form */
-    std::vector<Eigen::MatrixXd> blockTinvR_;
+    /*! 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