From b3ae9cdb1947e46d577da2478accba4e9c192e11 Mon Sep 17 00:00:00 2001 From: Alejandro Gallo Date: Fri, 17 Nov 2023 01:34:38 +0100 Subject: [PATCH] Fix the IP and EA preconditioners for mpirun and add abstract Preconditioner --- src/algorithms/CcsdPreconditioner.cxx | 21 ++++---- src/algorithms/CcsdPreconditioner.hpp | 49 ++++++++++++++++--- .../UCcsdIPEquationOfMotionDavidson.cxx | 7 ++- .../UccsdAmplitudesFromCoulombIntegrals.cxx | 3 +- 4 files changed, 56 insertions(+), 24 deletions(-) diff --git a/src/algorithms/CcsdPreconditioner.cxx b/src/algorithms/CcsdPreconditioner.cxx index 66e9f9f75..ebb03eb1a 100644 --- a/src/algorithms/CcsdPreconditioner.cxx +++ b/src/algorithms/CcsdPreconditioner.cxx @@ -392,7 +392,6 @@ void EACcsdPreconditioner::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"]; } @@ -545,21 +544,23 @@ EACcsdPreconditioner::getCorrection(const sisi4s::complex lambda, template void IPCcsdPreconditioner::calculateDiagonal() { + LOG(1, "IPCcsdPreconditioner") << "Getting No and Nv" << std::endl; int No(this->Fij->lens[0]), Nv(this->Fab->lens[0]); std::vector 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, std::vector)>( {NEW(Tensor, 1, o.data(), ns.data(), *Sisi4s::world, "Di"), NEW(Tensor, 3, voo.data(), nss.data(), *Sisi4s::world, "Daij")}), std::vector({"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"]; @@ -567,7 +568,6 @@ void IPCcsdPreconditioner::calculateDiagonal() { (*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 @@ -576,7 +576,7 @@ IPCcsdPreconditioner::getInitialBasis(const int eigenVectorsCount) { calculateDiagonal(); LOG(0, "IPCcsdPreconditioner") << "Getting initial basis " << std::endl; // find K=eigenVectorsCount lowest diagonal elements at each processor - std::vector> localElements(this->diagonalH->readLocal()); + std::vector> localElements(diagonalH->readLocal()); std::sort(localElements.begin(), localElements.end(), EomDiagonalValueComparator()); @@ -616,8 +616,9 @@ IPCcsdPreconditioner::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 basisElement(*this->diagonalH); + SDFockVector basisElement(*diagonalH); basisElement *= 0.0; std::vector> elements; if (communicator.getRank() == 0) { @@ -676,7 +677,7 @@ template SDFockVector IPCcsdPreconditioner::getCorrection(const sisi4s::complex lambda, SDFockVector &residuum) { - SDFockVector w(*this->diagonalH); + SDFockVector w(*diagonalH); // Define a singleton helping class for the diagonal correction class DiagonalCorrection { @@ -693,7 +694,7 @@ IPCcsdPreconditioner::getCorrection(const sisi4s::complex lambda, double lambda; } diagonalCorrection(std::real(lambda)); - SDFockVector correction(*this->diagonalH); + SDFockVector 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()); @@ -701,7 +702,7 @@ IPCcsdPreconditioner::getCorrection(const sisi4s::complex lambda, .contract(1.0, *residuum.get(c), indices, - *this->diagonalH->get(c), + *diagonalH->get(c), indices, 0.0, indices, diff --git a/src/algorithms/CcsdPreconditioner.hpp b/src/algorithms/CcsdPreconditioner.hpp index 2a20bc7fa..1bfa2e4c0 100644 --- a/src/algorithms/CcsdPreconditioner.hpp +++ b/src/algorithms/CcsdPreconditioner.hpp @@ -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 #include #include @@ -16,6 +23,28 @@ namespace sisi4s { +template +class Preconditioner { +public: + virtual void calculateDiagonal() = 0; + + virtual std::vector getInitialBasis(int eigenVectorsCount) = 0; + + virtual V getCorrection(const complex eigenValue, V &residuum) = 0; + + /** + * \brief Setters for the main tensors + */ + __DEFINE_SETTER(Tensor *, Tai, nullptr); + __DEFINE_SETTER(Tensor *, Tabij, nullptr); + __DEFINE_SETTER(Tensor *, Fij, nullptr); + __DEFINE_SETTER(Tensor *, Fab, nullptr); + __DEFINE_SETTER(Tensor *, Vabcd, nullptr); + __DEFINE_SETTER(Tensor *, Vijab, nullptr); + __DEFINE_SETTER(Tensor *, Viajb, nullptr); + __DEFINE_SETTER(Tensor *, 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 @@ -94,25 +123,29 @@ class CcsdPreconditioner { }; template -class IPCcsdPreconditioner : public CcsdPreconditioner { +class IPCcsdPreconditioner : public Preconditioner> { public: + using V = SDFockVector; + PTR(V) diagonalH; + void calculateDiagonal(); - std::vector> getInitialBasis(int eigenVectorsCount); + std::vector getInitialBasis(int eigenVectorsCount); - SDFockVector getCorrection(const complex eigenValue, - SDFockVector &residuum); + V getCorrection(const complex eigenValue, V &residuum); }; template -class EACcsdPreconditioner : public CcsdPreconditioner { +class EACcsdPreconditioner : public Preconditioner> { public: + using V = SDFockVector; + PTR(V) diagonalH; + void calculateDiagonal(); - std::vector> getInitialBasis(int eigenVectorsCount); + std::vector getInitialBasis(int eigenVectorsCount); - SDFockVector getCorrection(const complex eigenValue, - SDFockVector &residuum); + V getCorrection(const complex eigenValue, V &residuum); }; template diff --git a/src/algorithms/UCcsdIPEquationOfMotionDavidson.cxx b/src/algorithms/UCcsdIPEquationOfMotionDavidson.cxx index 148ab598c..13c28a580 100644 --- a/src/algorithms/UCcsdIPEquationOfMotionDavidson.cxx +++ b/src/algorithms/UCcsdIPEquationOfMotionDavidson.cxx @@ -252,7 +252,8 @@ void UCcsdIPEquationOfMotionDavidson::run() { ipH.h = &H; // INITIALIZE SIMILARITY PRECONDITIONER - IPCcsdPreconditioner P; + using _Preconditioner = IPCcsdPreconditioner; + _Preconditioner P; P.setTai(&Tai) .setTabij(&Tabij) .setFij(Fij) @@ -260,9 +261,7 @@ void UCcsdIPEquationOfMotionDavidson::run() { // Set coulomb integrals .setVijab(Vijab); - EigenSystemDavidsonMono, - SDFockVector> + EigenSystemDavidsonMono> eigenSystem(&ipH, eigenStates, &P, diff --git a/src/algorithms/UccsdAmplitudesFromCoulombIntegrals.cxx b/src/algorithms/UccsdAmplitudesFromCoulombIntegrals.cxx index 7dab7c7c7..39a7c08ea 100644 --- a/src/algorithms/UccsdAmplitudesFromCoulombIntegrals.cxx +++ b/src/algorithms/UccsdAmplitudesFromCoulombIntegrals.cxx @@ -65,9 +65,8 @@ void UccsdAmplitudesFromCoulombIntegrals::run() { "PHHPCoulombIntegrals", "UccsdDoublesAmplitudes", "UccsdSinglesAmplitudes", - "UccsdEnergy" + "UccsdEnergy", // - , "onlyPPL", "initialDoublesAmplitudes", "initialSinglesAmplitudes",