Skip to content

Commit

Permalink
Fix the IP and EA preconditioners for mpirun and add abstract Precond…
Browse files Browse the repository at this point in the history 
…itioner
  • Loading branch information
@alejandrogallo
alejandrogallo committed Nov 17, 2023
1 parent f5323d2 commit b3ae9cd
Show file tree
Hide file tree
Showing 4 changed files with 56 additions and 24 deletions.
21 changes: 11 additions & 10 deletions src/algorithms/CcsdPreconditioner.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -392,7 +392,6 @@ void EACcsdPreconditioner<F>::calculateDiagonal() {
(*Da)["a"] = (+1.0) * (*Fab)["aa"];

(*Dabi)["cdi"] = (-1.0) * (*Fij)["ii"];
(*Dabi)["cdi"] += (-1.0) * (*Fij)["jj"];
(*Dabi)["cdi"] += (+1.0) * (*Fab)["cc"];
(*Dabi)["cdi"] += (+1.0) * (*Fab)["dd"];
}
Expand Down Expand Up @@ -545,29 +544,30 @@ EACcsdPreconditioner<F>::getCorrection(const sisi4s::complex lambda,

template <typename F>
void IPCcsdPreconditioner<F>::calculateDiagonal() {
LOG(1, "IPCcsdPreconditioner") << "Getting No and Nv" << std::endl;
int No(this->Fij->lens[0]), Nv(this->Fab->lens[0]);
std::vector<int> o{{No}}, voo{{Nv, No, No}}, ns{{NS}}, nss{{NS, NS, NS}};

auto &Fij = this->Fij;
auto &Fab = this->Fab;

this->diagonalH = NEW(
LOG(1, "IPCcsdPreconditioner") << "Allocate diagonal" << std::endl;
diagonalH = NEW(
SDFockVector<F>,
std::vector<PTR(Tensor<F>)>(
{NEW(Tensor<F>, 1, o.data(), ns.data(), *Sisi4s::world, "Di"),
NEW(Tensor<F>, 3, voo.data(), nss.data(), *Sisi4s::world, "Daij")}),
std::vector<std::string>({"i", "aij"}));

auto Di(this->diagonalH->get(0));
auto Daij(this->diagonalH->get(1));
auto Di(diagonalH->get(0));
auto Daij(diagonalH->get(1));

// calculate diagonal elements of H
(*Di)["i"] = (-1.0) * (*Fij)["ii"];

(*Daij)["cij"] = (-1.0) * (*Fij)["ii"];
(*Daij)["cij"] += (-1.0) * (*Fij)["jj"];
(*Daij)["cij"] += (+1.0) * (*Fab)["cc"];
(*Daij)["cij"] += (+1.0) * (*Fab)["dd"];
}

template <typename F>
Expand All @@ -576,7 +576,7 @@ IPCcsdPreconditioner<F>::getInitialBasis(const int eigenVectorsCount) {
calculateDiagonal();
LOG(0, "IPCcsdPreconditioner") << "Getting initial basis " << std::endl;
// find K=eigenVectorsCount lowest diagonal elements at each processor
std::vector<std::pair<size_t, F>> localElements(this->diagonalH->readLocal());
std::vector<std::pair<size_t, F>> localElements(diagonalH->readLocal());
std::sort(localElements.begin(),
localElements.end(),
EomDiagonalValueComparator<F>());
Expand Down Expand Up @@ -616,8 +616,9 @@ IPCcsdPreconditioner<F>::getInitialBasis(const int eigenVectorsCount) {
int currentEigenVectorCount(0);
unsigned int b(0);
int zeroVectorCount(0);
LOG(1, "IPCcsdPreconditioner") << "Starting with loop" << std::endl;
while (currentEigenVectorCount < eigenVectorsCount) {
SDFockVector<F> basisElement(*this->diagonalH);
SDFockVector<F> basisElement(*diagonalH);
basisElement *= 0.0;
std::vector<std::pair<size_t, F>> elements;
if (communicator.getRank() == 0) {
Expand Down Expand Up @@ -676,7 +677,7 @@ template <typename F>
SDFockVector<F>
IPCcsdPreconditioner<F>::getCorrection(const sisi4s::complex lambda,
SDFockVector<F> &residuum) {
SDFockVector<F> w(*this->diagonalH);
SDFockVector<F> w(*diagonalH);

// Define a singleton helping class for the diagonal correction
class DiagonalCorrection {
Expand All @@ -693,15 +694,15 @@ IPCcsdPreconditioner<F>::getCorrection(const sisi4s::complex lambda,
double lambda;
} diagonalCorrection(std::real(lambda));

SDFockVector<F> correction(*this->diagonalH);
SDFockVector<F> correction(*diagonalH);
// compute ((lambda * id - Diag(diagonal))^-1) . residuum
for (unsigned int c(0); c < w.get_components_count(); ++c) {
const char *indices(correction.component_indices[c].c_str());
(*correction.get(c))
.contract(1.0,
*residuum.get(c),
indices,
*this->diagonalH->get(c),
*diagonalH->get(c),
indices,
0.0,
indices,
Expand Down
49 changes: 41 additions & 8 deletions src/algorithms/CcsdPreconditioner.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,13 @@
} \
type name = default

#define __DEFINE_SETTER(type, name, default) \
Preconditioner &set##name(type t) { \
name = t; \
return *this; \
} \
type name = default

#include <algorithms/Algorithm.hpp>
#include <math/FockVector.hpp>
#include <vector>
Expand All @@ -16,6 +23,28 @@

namespace sisi4s {

template <typename F, typename V>
class Preconditioner {
public:
virtual void calculateDiagonal() = 0;

virtual std::vector<V> getInitialBasis(int eigenVectorsCount) = 0;

virtual V getCorrection(const complex eigenValue, V &residuum) = 0;

/**
* \brief Setters for the main tensors
*/
__DEFINE_SETTER(Tensor<F> *, Tai, nullptr);
__DEFINE_SETTER(Tensor<F> *, Tabij, nullptr);
__DEFINE_SETTER(Tensor<F> *, Fij, nullptr);
__DEFINE_SETTER(Tensor<F> *, Fab, nullptr);
__DEFINE_SETTER(Tensor<F> *, Vabcd, nullptr);
__DEFINE_SETTER(Tensor<F> *, Vijab, nullptr);
__DEFINE_SETTER(Tensor<F> *, Viajb, nullptr);
__DEFINE_SETTER(Tensor<F> *, Vijkl, nullptr);
};

/**
* \brief Implements the diagonal preconditionar for the davidson method
* \tparam F It is the field variable to be used, in general it will be
Expand Down Expand Up @@ -94,25 +123,29 @@ class CcsdPreconditioner {
};

template <typename F>
class IPCcsdPreconditioner : public CcsdPreconditioner<F> {
class IPCcsdPreconditioner : public Preconditioner<F, SDFockVector<F>> {
public:
using V = SDFockVector<F>;
PTR(V) diagonalH;

void calculateDiagonal();

std::vector<SDFockVector<F>> getInitialBasis(int eigenVectorsCount);
std::vector<V> getInitialBasis(int eigenVectorsCount);

SDFockVector<F> getCorrection(const complex eigenValue,
SDFockVector<F> &residuum);
V getCorrection(const complex eigenValue, V &residuum);
};

template <typename F>
class EACcsdPreconditioner : public CcsdPreconditioner<F> {
class EACcsdPreconditioner : public Preconditioner<F, SDFockVector<F>> {
public:
using V = SDFockVector<F>;
PTR(V) diagonalH;

void calculateDiagonal();

std::vector<SDFockVector<F>> getInitialBasis(int eigenVectorsCount);
std::vector<V> getInitialBasis(int eigenVectorsCount);

SDFockVector<F> getCorrection(const complex eigenValue,
SDFockVector<F> &residuum);
V getCorrection(const complex eigenValue, V &residuum);
};

template <typename F>
Expand Down
7 changes: 3 additions & 4 deletions src/algorithms/UCcsdIPEquationOfMotionDavidson.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -252,17 +252,16 @@ void UCcsdIPEquationOfMotionDavidson::run() {
ipH.h = &H;

// INITIALIZE SIMILARITY PRECONDITIONER
IPCcsdPreconditioner<F> P;
using _Preconditioner = IPCcsdPreconditioner<F>;
_Preconditioner P;
P.setTai(&Tai)
.setTabij(&Tabij)
.setFij(Fij)
.setFab(Fab)
// Set coulomb integrals
.setVijab(Vijab);

EigenSystemDavidsonMono<IPHamiltonian,
IPCcsdPreconditioner<F>,
SDFockVector<F>>
EigenSystemDavidsonMono<IPHamiltonian, _Preconditioner, SDFockVector<F>>
eigenSystem(&ipH,
eigenStates,
&P,
Expand Down
3 changes: 1 addition & 2 deletions src/algorithms/UccsdAmplitudesFromCoulombIntegrals.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,9 +65,8 @@ void UccsdAmplitudesFromCoulombIntegrals::run() {
"PHHPCoulombIntegrals",
"UccsdDoublesAmplitudes",
"UccsdSinglesAmplitudes",
"UccsdEnergy"
"UccsdEnergy",
//
,
"onlyPPL",
"initialDoublesAmplitudes",
"initialSinglesAmplitudes",
Expand Down

0 comments on commit b3ae9cd

Please sign in to comment.