Skip to content

Residual compute #384

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
jeanlucf22 merged 2 commits into from
Nov 20, 2025
Merged

Residual compute #384

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
5e68097
Refactor H application and R computation
jeanlucf22 Nov 16, 2025
5a3b0bf
Rm optional args to Hamiltonian::addHlocal2matrix()
jeanlucf22 Nov 17, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 3 additions & 1 deletion src/Hamiltonian.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -141,7 +141,7 @@ void Hamiltonian<T>::applyLocal(const int ncolors, T& phi, T& hphi)
for (int i = 0; i < ncolors; i++)
{
using memory_space_type = typename T::memory_space_type;
ORBDTYPE* ihphi = hphi.getPsi(i);
auto ihphi = hphi.getPsi(i);
unsigned int const size = hphi.getNumpt();
ORBDTYPE* ihphi_host_view = MemorySpace::Memory<ORBDTYPE,
memory_space_type>::allocate_host_view(size);
Expand Down Expand Up @@ -221,6 +221,8 @@ template <>
void Hamiltonian<LocGridOrbitals>::addHlocal2matrix(LocGridOrbitals& phi1,
LocGridOrbitals& phi2, ReplicatedMatrix& hij, const bool force)
{
(void)hij;

applyLocal(phi2, force);

// phi1.addDotWithNcol2Matrix(*hlphi_, hij);
Expand Down
2 changes: 1 addition & 1 deletion src/Hamiltonian.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@ class Hamiltonian

template <class MatrixType>
void addHlocal2matrix(OrbitalsType& orbitals1, OrbitalsType& orbitals2,
MatrixType& mat, const bool force = false);
MatrixType& mat, const bool force);
void addHlocalij(OrbitalsType& orbitals1, OrbitalsType& orbitals2,
ProjectedMatricesInterface*);
void addHlocalij(OrbitalsType& orbitals1, ProjectedMatricesInterface*);
Expand Down
8 changes: 4 additions & 4 deletions src/HamiltonianMVPSolver.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -151,7 +151,7 @@ int HamiltonianMVPSolver<MatrixType, ProjMatrixType, OrbitalsType>::solve(
// compute new h11 for the current potential by adding local part to
// nonlocal components
h11 = h11nl;
hamiltonian_->addHlocal2matrix(orbitals, orbitals, h11);
hamiltonian_->addHlocal2matrix(orbitals, orbitals, h11, false);

projmatrices->assignH(h11);
projmatrices->setHB2H();
Expand Down Expand Up @@ -179,7 +179,7 @@ int HamiltonianMVPSolver<MatrixType, ProjMatrixType, OrbitalsType>::solve(

// update H and compute energy at midpoint
h11 = h11nl;
hamiltonian_->addHlocal2matrix(orbitals, orbitals, h11);
hamiltonian_->addHlocal2matrix(orbitals, orbitals, h11, false);

projmatrices->assignH(h11);
projmatrices->setHB2H();
Expand Down Expand Up @@ -214,7 +214,7 @@ int HamiltonianMVPSolver<MatrixType, ProjMatrixType, OrbitalsType>::solve(

// update H with new potential
h11 = h11nl;
hamiltonian_->addHlocal2matrix(orbitals, orbitals, h11);
hamiltonian_->addHlocal2matrix(orbitals, orbitals, h11, false);

projmatrices->assignH(h11);
projmatrices->setHB2H();
Expand Down Expand Up @@ -270,7 +270,7 @@ int HamiltonianMVPSolver<MatrixType, ProjMatrixType, OrbitalsType>::solve(

// update H
h11 = h11nl;
hamiltonian_->addHlocal2matrix(orbitals, orbitals, h11);
hamiltonian_->addHlocal2matrix(orbitals, orbitals, h11, false);

projmatrices->assignH(h11);
projmatrices->setHB2H();
Expand Down
55 changes: 16 additions & 39 deletions src/MGmol.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1171,53 +1171,43 @@ void MGmol<OrbitalsType>::projectOutKernel(OrbitalsType& phi)
}

template <class OrbitalsType>
void MGmol<OrbitalsType>::setGamma(
const pb::Lap<ORBDTYPE>& lapOper, const Potentials& pot)
void MGmol<OrbitalsType>::precond_mg(OrbitalsType& phi)
{
assert(orbitals_precond_);

Control& ct = *(Control::instance());

orbitals_precond_->setGamma(
lapOper, pot, ct.getMGlevels(), proj_matrices_.get());
}
Potentials& pot = hamiltonian_->potential();
pb::Lap<ORBDTYPE>* lapOper = hamiltonian_->lapOper();

template <class OrbitalsType>
void MGmol<OrbitalsType>::precond_mg(OrbitalsType& phi)
{
assert(orbitals_precond_);
orbitals_precond_->setGamma(
*lapOper, pot, ct.getMGlevels(), proj_matrices_.get());

orbitals_precond_->precond_mg(phi);
}

template <class OrbitalsType>
double MGmol<OrbitalsType>::computeResidual(OrbitalsType& orbitals,
OrbitalsType& work_orbitals, Ions& ions, OrbitalsType& res,
const bool print_residual, const bool norm_res)
double MGmol<OrbitalsType>::computeResidual(OrbitalsType& phi,
OrbitalsType& hphi, Ions& ions, OrbitalsType& res,
const KBPsiMatrixSparse* const kbpsi, const bool print_residual,
const bool norm_res)

{
assert(orbitals.getIterativeIndex() >= 0);

comp_res_tm_.start();
// os_<<"computeResidual()"<<endl;

Control& ct(*(Control::instance()));

proj_matrices_->computeInvB();

Potentials& pot = hamiltonian_->potential();
pb::Lap<ORBDTYPE>* lapop = hamiltonian_->lapOper();

setGamma(*lapop, pot);

// get H*psi stored in work_orbitals.psi
// get H*phi stored in hphi
// and psi^T H psi in Hij
getHpsiAndTheta(ions, orbitals, work_orbitals);
getHpsiAndTheta(ions, phi, hphi, kbpsi);

double norm2Res = computeConstraintResidual(
orbitals, work_orbitals, res, print_residual, norm_res);
double norm2Res
= computeConstraintResidual(phi, hphi, res, print_residual, norm_res);

if (ct.isSpreadFunctionalEnergy()) addResidualSpreadPenalty(orbitals, res);
if (ct.isSpreadFunctionalEnergy()) addResidualSpreadPenalty(phi, res);

comp_res_tm_.stop();

Expand Down Expand Up @@ -1344,19 +1334,8 @@ double MGmol<OrbitalsType>::computePrecondResidual(OrbitalsType& phi,
{
Control& ct = *(Control::instance());

proj_matrices_->computeInvB();

Potentials& pot = hamiltonian_->potential();
pb::Lap<ORBDTYPE>* lapop = hamiltonian_->lapOper();

setGamma(*lapop, pot);

// get H*psi stored in hphi
// and psi^T H psi in Hij
getHpsiAndTheta(ions, phi, hphi, kbpsi);

double norm2Res
= computeConstraintResidual(phi, hphi, res, print_residual, norm_res);
double norm2Res = computeResidual(
phi, hphi, ions, res, kbpsi, print_residual, norm_res);

if (ct.withPreconditioner())
{
Expand All @@ -1366,8 +1345,6 @@ double MGmol<OrbitalsType>::computePrecondResidual(OrbitalsType& phi,
orbitals_precond_->precond_mg(res);
}

// if( ct.isSpreadFunctionalActive() )addResidualSpreadPenalty(phi,res);

return norm2Res;
}

Expand Down
11 changes: 9 additions & 2 deletions src/MGmol.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -321,10 +321,17 @@ class MGmol : public MGmolInterface
void projectOutKernel(OrbitalsType& phi);

void precond_mg(OrbitalsType& orbitals);
void setGamma(const pb::Lap<ORBDTYPE>& lapOper, const Potentials& pot);
double computeResidual(OrbitalsType& orbitals, OrbitalsType& work_orbitals,
Ions& ions, OrbitalsType& res, const KBPsiMatrixSparse* const kbpsi,
const bool print_residual, const bool norm_res);
double computeResidual(OrbitalsType& orbitals, OrbitalsType& work_orbitals,
Ions& ions, OrbitalsType& res, const bool print_residual,
const bool norm_res);
const bool norm_res)
{
return computeResidual(orbitals, work_orbitals, ions, res,
g_kbpsi_.get(), print_residual, norm_res);
}

void applyAOMMprojection(OrbitalsType&);
void force(OrbitalsType& orbitals, Ions& ions)
{
Expand Down
9 changes: 4 additions & 5 deletions src/MVPSolver.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -239,7 +239,7 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
// compute h11 for the current potential by adding local part to
// nonlocal components
MatrixType h11(h11_nl);
hamiltonian_->addHlocal2matrix(orbitals, orbitals, h11);
hamiltonian_->addHlocal2matrix(orbitals, orbitals, h11, false);

current_proj_mat->assignH(h11);
current_proj_mat->setHB2H();
Expand Down Expand Up @@ -317,10 +317,9 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
energy_->saveVofRho();

// update h11
{
h11 = h11_nl;
hamiltonian_->addHlocal2matrix(orbitals, orbitals, h11);
}
h11 = h11_nl;
hamiltonian_->addHlocal2matrix(
orbitals, orbitals, h11, false);

proj_mat_work_->assignH(h11);
proj_mat_work_->setHB2H();
Expand Down
19 changes: 9 additions & 10 deletions src/computeHij.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -197,7 +197,7 @@ void MGmol<OrbitalsType>::computeHij_private(OrbitalsType& orbitals_i,
ss2dm->accumulate(submat, hij, 0.);

// add local Hamiltonian part to phi^T*H*phi
hamiltonian_->addHlocal2matrix(orbitals_i, orbitals_j, hij);
hamiltonian_->addHlocal2matrix(orbitals_i, orbitals_j, hij, false);
}

template <class OrbitalsType>
Expand Down Expand Up @@ -324,13 +324,6 @@ void MGmol<OrbitalsType>::computeHnlPhiAndAdd2HPhi(Ions& ions,
hphi.setIterativeIndex(phi.getIterativeIndex());
}

template <class OrbitalsType>
void MGmol<OrbitalsType>::getHpsiAndTheta(
Ions& ions, OrbitalsType& phi, OrbitalsType& hphi)
{
getHpsiAndTheta(ions, phi, hphi, g_kbpsi_.get());
}

template <class OrbitalsType>
void MGmol<OrbitalsType>::getHpsiAndTheta(Ions& ions, OrbitalsType& phi,
OrbitalsType& hphi, const KBPsiMatrixSparse* const kbpsi)
Expand All @@ -345,7 +338,7 @@ void MGmol<OrbitalsType>::getHpsiAndTheta(Ions& ions, OrbitalsType& phi,
os_ << " getHpsiAndTheta" << std::endl;
#endif

hphi.assign(hamiltonian_->applyLocal(phi));
hamiltonian_->applyLocal(phi.chromatic_number(), phi, hphi);

// Compute "nstates" columns of matrix
// Hij = phi**T * H_loc * phi and save in sh
Expand All @@ -370,7 +363,13 @@ void MGmol<OrbitalsType>::getHpsiAndTheta(Ions& ions, OrbitalsType& phi,
kbpsi->computeHvnlMatrix(ions, proj_matrices_.get());

// add local part of H to sh
hamiltonian_->addHlocalij(phi, proj_matrices_.get());
SquareLocalMatrices<MATDTYPE, MemorySpace::Host> slm(
phi.subdivx(), phi.chromatic_number());

phi.computeLocalProduct(hphi, slm);
proj_matrices_->setLocalMatrixElementsHl(slm);

proj_matrices_->consolidateH();

energy_->saveVofRho();

Expand Down