Skip to content

Refactor: clean Davidson code, remove deprecated overlap matrix scc related code #6330

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
merged 3 commits into from
Jun 24, 2025
Merged

Refactor: clean Davidson code, remove deprecated overlap matrix scc related code #6330

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
c120f6e
Remove scc calculation deprecated code of traditional Davidson algori…
Cstandardlib Jun 24, 2025
eb44e46
Remove scc completely in traditional Davidson.
Cstandardlib Jun 24, 2025
2ad90cd
Fix test after removal of scc
Cstandardlib Jun 24, 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
57 changes: 6 additions & 51 deletions source/source_hsolver/diago_david.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -76,10 +76,6 @@ DiagoDavid<T, Device>::DiagoDavid(const Real* precondition_in,
resmem_complex_op()(this->hcc, nbase_x * nbase_x, "DAV::hcc");
setmem_complex_op()(this->hcc, 0, nbase_x * nbase_x);

// scc(nbase_x, nbase_x); // Overlap on the reduced basis
// resmem_complex_op()(this->ctx, this->scc, nbase_x * nbase_x, "DAV::scc");
// setmem_complex_op()(this->ctx, this->scc, 0, nbase_x * nbase_x);

// vcc(nbase_x, nbase_x); // Eigenvectors of hcc
resmem_complex_op()(this->vcc, nbase_x * nbase_x, "DAV::vcc");
setmem_complex_op()(this->vcc, 0, nbase_x * nbase_x);
Expand All @@ -106,7 +102,6 @@ DiagoDavid<T, Device>::~DiagoDavid()
delmem_complex_op()(this->hpsi);
delmem_complex_op()(this->spsi);
delmem_complex_op()(this->hcc);
// delmem_complex_op()(this->ctx, this->scc);
delmem_complex_op()(this->vcc);
delmem_complex_op()(this->lagrange_matrix);
base_device::memory::delete_memory_op<Real, base_device::DEVICE_CPU>()(this->eigenvalue);
Expand Down Expand Up @@ -190,9 +185,9 @@ int DiagoDavid<T, Device>::diag_once(const HPsiFunc& hpsi_func,
// slice index in this piece of code is in C manner. i.e. 0:id stands for [0,id)
hpsi_func(basis, hpsi, dim, nband);

this->cal_elem(dim, nbase, nbase_x, this->notconv, this->hpsi, this->spsi, this->hcc, this->scc);
this->cal_elem(dim, nbase, nbase_x, this->notconv, this->hpsi, this->spsi, this->hcc);

this->diag_zhegvx(nbase, nband, this->hcc, this->scc, nbase_x, this->eigenvalue, this->vcc);
this->diag_zhegvx(nbase, nband, this->hcc, nbase_x, this->eigenvalue, this->vcc);

for (int m = 0; m < nband; m++)
{
Expand All @@ -218,9 +213,9 @@ int DiagoDavid<T, Device>::diag_once(const HPsiFunc& hpsi_func,
unconv.data(),
this->eigenvalue);

this->cal_elem(dim, nbase, nbase_x, this->notconv, this->hpsi, this->spsi, this->hcc, this->scc);
this->cal_elem(dim, nbase, nbase_x, this->notconv, this->hpsi, this->spsi, this->hcc);

this->diag_zhegvx(nbase, nband, this->hcc, this->scc, nbase_x, this->eigenvalue, this->vcc);
this->diag_zhegvx(nbase, nband, this->hcc, nbase_x, this->eigenvalue, this->vcc);

// check convergence and update eigenvalues
ModuleBase::timer::tick("DiagoDavid", "check_update");
Expand Down Expand Up @@ -282,7 +277,6 @@ int DiagoDavid<T, Device>::diag_once(const HPsiFunc& hpsi_func,
this->hpsi,
this->spsi,
this->hcc,
this->scc,
this->vcc);
ModuleBase::timer::tick("DiagoDavid", "last");
}
Expand Down Expand Up @@ -536,8 +530,7 @@ void DiagoDavid<T, Device>::cal_elem(const int& dim,
const int& notconv, // number of newly added basis vectors
const T* hpsi,
const T* spsi,
T* hcc,
T* scc)
T* hcc)
{
if (test_david == 1) {
ModuleBase::TITLE("DiagoDavid", "cal_elem");
Expand All @@ -562,28 +555,12 @@ void DiagoDavid<T, Device>::cal_elem(const int& dim,
this->zero,
hcc + nbase, // notconv * (nbase + notconv)
nbase_x);
// scc[nbase] = basis[nbase]' * spsi
// gemm_op<T, Device>()(this->ctx,
// 'C',
// 'N',
// notconv,
// nbase + notconv,
// dim,
// this->one,
// basis + dim*nbase, // dim * notconv
// dim,
// spsi, // dim * (nbase + notconv)
// dim,
// this->zero,
// scc + nbase, // notconv * (nbase + notconv)
// nbase_x);


#ifdef __MPI
if (diag_comm.nproc > 1)
{
ModuleBase::matrixTranspose_op<T, Device>()(nbase_x, nbase_x, hcc, hcc);
// matrixTranspose_op<T, Device>()(this->ctx, nbase_x, nbase_x, scc, scc);

auto* swap = new T[notconv * nbase_x];
syncmem_complex_op()(swap, hcc + nbase * nbase_x, notconv * nbase_x);
Expand All @@ -599,21 +576,13 @@ void DiagoDavid<T, Device>::cal_elem(const int& dim,
else {
MPI_Reduce(swap, hcc + nbase * nbase_x, notconv * nbase_x, MPI_DOUBLE_COMPLEX, MPI_SUM, 0, diag_comm.comm);
}
// syncmem_complex_op()(this->ctx, this->ctx, swap, scc + nbase * nbase_x, notconv * nbase_x);
if (base_device::get_current_precision(swap) == "single") {
// MPI_Reduce(swap, scc + nbase * nbase_x, notconv * nbase_x, MPI_COMPLEX, MPI_SUM, 0, diag_comm.comm);
}
else {
// MPI_Reduce(swap, scc + nbase * nbase_x, notconv * nbase_x, MPI_DOUBLE_COMPLEX, MPI_SUM, 0, diag_comm.comm);
}

}
delete[] swap;

// Parallel_Reduce::reduce_complex_double_pool( hcc + nbase * nbase_x, notconv * nbase_x );
// Parallel_Reduce::reduce_complex_double_pool( scc + nbase * nbase_x, notconv * nbase_x );

ModuleBase::matrixTranspose_op<T, Device>()(nbase_x, nbase_x, hcc, hcc);
// matrixTranspose_op<T, Device>()(this->ctx, nbase_x, nbase_x, scc, scc);
}
#endif

Expand All @@ -637,7 +606,6 @@ template <typename T, typename Device>
void DiagoDavid<T, Device>::diag_zhegvx(const int& nbase,
const int& nband,
const T* hcc,
const T* /*scc*/,
const int& nbase_x,
Real* eigenvalue, // in CPU
T* vcc)
Expand Down Expand Up @@ -694,7 +662,6 @@ void DiagoDavid<T, Device>::refresh(const int& dim,
T* hpsi,
T* spsi,
T* hcc,
T* scc,
T* vcc)
{
if (test_david == 1) {
Expand Down Expand Up @@ -763,42 +730,31 @@ void DiagoDavid<T, Device>::refresh(const int& dim,

setmem_complex_op()(hcc, 0, nbase_x * nbase_x);

// setmem_complex_op()(this->ctx, scc, 0, nbase_x * nbase_x);

if (this->device == base_device::GpuDevice)
{
#if defined(__CUDA) || defined(__ROCM)
T* hcc_cpu = nullptr;
// T* scc_cpu = nullptr;
T* vcc_cpu = nullptr;
base_device::memory::resize_memory_op<T, base_device::DEVICE_CPU>()(hcc_cpu,
nbase_x * nbase_x,
"DAV::hcc");
// base_device::memory::resize_memory_op<T, base_device::DEVICE_CPU>()(this->cpu_ctx,
// scc_cpu,
// nbase_x * nbase_x,
// "DAV::scc");
base_device::memory::resize_memory_op<T, base_device::DEVICE_CPU>()(vcc_cpu,
nbase_x * nbase_x,
"DAV::vcc");

syncmem_d2h_op()(hcc_cpu, hcc, nbase_x * nbase_x);
// syncmem_d2h_op()(this->cpu_ctx, this->ctx, scc_cpu, scc, nbase_x * nbase_x);
syncmem_d2h_op()(vcc_cpu, vcc, nbase_x * nbase_x);

for (int i = 0; i < nbase; i++)
{
hcc_cpu[i * nbase_x + i] = eigenvalue_in[i];
// scc_cpu[i * nbase_x + i] = this->one[0];
vcc_cpu[i * nbase_x + i] = this->one[0];
}

syncmem_h2d_op()(hcc, hcc_cpu, nbase_x * nbase_x);
// syncmem_h2d_op()(this->ctx, this->cpu_ctx, scc, scc_cpu, nbase_x * nbase_x);
syncmem_h2d_op()(vcc, vcc_cpu, nbase_x * nbase_x);

base_device::memory::delete_memory_op<T, base_device::DEVICE_CPU>()(hcc_cpu);
// base_device::memory::delete_memory_op<T, base_device::DEVICE_CPU>()(this->cpu_ctx, scc_cpu);
base_device::memory::delete_memory_op<T, base_device::DEVICE_CPU>()(vcc_cpu);
#endif
}
Expand All @@ -808,7 +764,6 @@ void DiagoDavid<T, Device>::refresh(const int& dim,
{
hcc[i * nbase_x + i] = eigenvalue_in[i];
// sc(i, i) = this->one;
// scc[i * nbase_x + i] = this->one[0];
// vc(i, i) = this->one;
vcc[i * nbase_x + i] = this->one[0];
}
Expand Down
11 changes: 2 additions & 9 deletions source/source_hsolver/diago_david.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,7 @@ class DiagoDavid
* @brief Destructor for the DiagoDavid class.
*
* This destructor releases the dynamically allocated memory used by the class members.
* It deletes the basis, hpsi, spsi, hcc, scc, vcc, lagrange_matrix, and eigenvalue arrays.
* It deletes the basis, hpsi, spsi, hcc, vcc, lagrange_matrix, and eigenvalue arrays.
*
*/
~DiagoDavid();
Expand Down Expand Up @@ -170,8 +170,6 @@ class DiagoDavid

T* hcc = nullptr; /// Hamiltonian on the reduced basis

T* scc = nullptr; /// overlap on the reduced basis

T* vcc = nullptr; /// eigenvectors of hc

T* lagrange_matrix = nullptr;
Expand Down Expand Up @@ -228,16 +226,14 @@ class DiagoDavid
* @param hpsi The output array for the Hamiltonian H times blockvector psi.
* @param spsi The output array for the overlap matrix S times blockvector psi.
* @param hcc Pointer to the array where the calculated Hamiltonian matrix elements will be stored.
* @param scc Pointer to the array where the calculated overlap matrix elements will be stored.
*/
void cal_elem(const int& dim,
int& nbase,
const int nbase_x,
const int& notconv,
const T* hpsi,
const T* spsi,
T* hcc,
T* scc);
T* hcc);

/**
* Refreshes the diagonalization solver by updating the basis and the reduced Hamiltonian.
Expand All @@ -252,7 +248,6 @@ class DiagoDavid
* @param hpsi Pointer to the output array for the updated basis set.
* @param spsi Pointer to the output array for the updated basis set (nband-th column).
* @param hcc Pointer to the output array for the updated reduced Hamiltonian.
* @param scc Pointer to the output array for the updated overlap matrix.
* @param vcc Pointer to the output array for the updated eigenvector matrix.
*
*/
Expand All @@ -266,7 +261,6 @@ class DiagoDavid
T* hpsi,
T* spsi,
T* hcc,
T* scc,
T* vcc);

/**
Expand Down Expand Up @@ -304,7 +298,6 @@ class DiagoDavid
void diag_zhegvx(const int& nbase,
const int& nband,
const T* hcc,
const T* scc,
const int& nbase_x,
Real* eigenvalue,
T* vcc);
Expand Down
1 change: 0 additions & 1 deletion source/source_hsolver/test/hsolver_pw_sup.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -142,7 +142,6 @@ DiagoDavid<T, Device>::~DiagoDavid() {
delmem_complex_op()(this->hpsi);
delmem_complex_op()(this->spsi);
delmem_complex_op()(this->hcc);
delmem_complex_op()(this->scc);
delmem_complex_op()(this->vcc);
delmem_complex_op()(this->lagrange_matrix);
base_device::memory::delete_memory_op<Real, base_device::DEVICE_CPU>()(this->eigenvalue);
Expand Down
Loading