feat: add PPCG eigenvalue solver implementation#7414
Open
OIerwsq wants to merge 1 commit into
Open
Conversation
There was a problem hiding this comment.
Pull request overview
Adds a new plane-wave Kohn–Sham eigensolver ppcg to ABACUS. A new hsolver::DiagoPPCG<T, Device> class is implemented (header + source), the PW dispatcher in HSolverPW is extended to construct and invoke the solver with the same hpsi/spsi/subspace callback shape as DiagoCG, the new value ppcg is added to the list of accepted ks_solver strings for basis_type=pw, the build is updated, and a GoogleTest unit test is registered. The implementation closely mirrors DiagoCG (per-band Schmidt orthogonalization, preconditioned gradient, Polak–Ribiere update, 2D line minimization, and an outer retry/subspace loop) rather than a true block PPCG algorithm.
Changes:
- New
DiagoPPCGclass withdiag()callback API, retry loop, and complex CPU+GPU template instantiations. - PW solver dispatch and input validation extended to accept
"ppcg". - New GoogleTest target
MODULE_HSOLVER_ppcgand a single parametric eigenvalue-comparison test against LAPACK.
Reviewed changes
Copilot reviewed 7 out of 7 changed files in this pull request and generated 8 comments.
Show a summary per file
| File | Description |
|---|---|
| source/source_hsolver/diago_ppcg.h | Declares DiagoPPCG class, callback typedefs, and helper methods. |
| source/source_hsolver/diago_ppcg.cpp | Implements PPCG diagonalization, including gradient/CG/line-minimization routines and template instantiations. |
| source/source_hsolver/hsolver_pw.cpp | Adds ppcg to the supported PW methods list and a dispatch branch constructing DiagoPPCG with subspace callback. |
| source/source_hsolver/CMakeLists.txt | Adds diago_ppcg.cpp to the hsolver object library. |
| source/source_hsolver/test/CMakeLists.txt | Registers the new MODULE_HSOLVER_ppcg test target under ENABLE_MPI. |
| source/source_hsolver/test/diago_ppcg_test.cpp | New unit test that compares PPCG eigenvalues against zheev_ on a random Hermitian matrix. |
| source/source_io/module_parameter/read_input_item_elec_stru.cpp | Adds "ppcg" to the validated PW ks_solver list. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| Real psi_norm = ModuleBase::dot_real_op<T, Device>()(this->n_basis_, psi_m.data<T>(), psi_m.data<T>()); | ||
| if (psi_norm <= 0.0) | ||
| { | ||
| std::cout << " DiagoPPGC: psi_norm <= 0.0, m = " << m << std::endl; |
Comment on lines
+604
to
+611
| avg = (this->n_band_ > 0) ? avg / this->n_band_ : 0; | ||
| avg_iter_ += avg; | ||
| ntry++; | ||
|
|
||
| } while (this->test_exit_cond(ntry, this->notconv_)); | ||
|
|
||
| ModuleBase::timer::end("DiagoPPCG", "diag"); | ||
| return avg; |
Comment on lines
+381
to
+386
| this->hpsi_func_ = hpsi_func; | ||
| this->spsi_func_ = spsi_func; | ||
|
|
||
| // Get convergence parameters from DiagoIterAssist | ||
| this->pw_diag_nmax_ = DiagoIterAssist<T, Device>::PW_DIAG_NMAX; | ||
| this->pw_diag_thr_ = DiagoIterAssist<T, Device>::PW_DIAG_THR; |
Comment on lines
+160
to
+164
| // g0 = P * S|grad> = prec * scg (element-wise multiply) | ||
| for (int i = 0; i < this->n_basis_; i++) | ||
| { | ||
| g0.data<T>()[i] = s_grad.data<T>()[i] * prec.data<Real>()[i]; | ||
| } |
Comment on lines
+493
to
+566
| for (int m = 0; m < this->n_band_; m++) | ||
| { | ||
| // Copy initial guess | ||
| psi_m.sync(psi[m]); | ||
|
|
||
| // Compute S|psi_m> | ||
| this->spsi_func_(psi_m.data<T>(), spsi_m.data<T>(), this->n_basis_, 1); | ||
|
|
||
| // Schmidt orthogonalize against lower bands | ||
| this->schmidt_orth(m, psi, spsi_m, psi_m); | ||
|
|
||
| // Recompute S|psi_m> after orthogonalization | ||
| this->spsi_func_(psi_m.data<T>(), spsi_m.data<T>(), this->n_basis_, 1); | ||
|
|
||
| // Compute H|psi_m> | ||
| this->hpsi_func_(psi_m.data<T>(), hpsi_m.data<T>(), this->n_basis_, 1); | ||
|
|
||
| // Initial eigenvalue estimate: λ = <ψ|H|ψ> | ||
| eigen_pack[m] = ModuleBase::dot_real_op<T, Device>()(this->n_basis_, psi_m.data<T>(), hpsi_m.data<T>()); | ||
|
|
||
| int iter = 0; | ||
| Real gg_last = 0.0; | ||
| Real cg_norm = 0.0; | ||
| Real theta = 0.0; | ||
| bool converged = false; | ||
| // g0 = P * S|grad> for CG update (persistent between iterations) | ||
| auto g0_tensor = ct::Tensor(ct::DataTypeToEnum<T>::value, | ||
| ct::DeviceTypeToEnum<ct_Device>::value, | ||
| {this->n_basis_}); | ||
|
|
||
| do | ||
| { | ||
| // Compute preconditioned gradient (ABACUS calc_grad) | ||
| this->compute_gradient(prec_tensor, hpsi_m, spsi_m, eigen_pack[m], grad, ppsi); | ||
|
|
||
| // Orthogonalize gradient against lower bands (ABACUS orth_grad) | ||
| this->orth_gradient(psi, m, grad, s_grad, lagrange); | ||
|
|
||
| // Update CG direction with persistent g0 and Polak-Ribiere (ABACUS calc_gamma_cg) | ||
| Real gamma = this->update_cg_direction(cg, scg, grad, s_grad, prec_tensor, g0_tensor, gg_last, iter); | ||
|
|
||
| // ABACUS CG correction: cg -= gamma * cg_norm * sin(theta) * psi_m | ||
| if (iter > 0 && std::abs(gamma) > std::numeric_limits<Real>::epsilon()) | ||
| { | ||
| const Real norma = gamma * cg_norm * std::sin(theta); | ||
| if (std::abs(norma) > std::numeric_limits<Real>::epsilon()) | ||
| { | ||
| T znorma = static_cast<T>(-norma); | ||
| ModuleBase::axpy_op<T, Device>()(this->n_basis_, &znorma, psi_m.data<T>(), 1, cg.data<T>(), 1); | ||
| ModuleBase::axpy_op<T, Device>()(this->n_basis_, &znorma, spsi_m.data<T>(), 1, scg.data<T>(), 1); | ||
| } | ||
| } | ||
|
|
||
| // Compute H|cg> and S|cg> | ||
| this->hpsi_func_(cg.data<T>(), ppsi.data<T>(), this->n_basis_, 1); | ||
| this->spsi_func_(cg.data<T>(), scg.data<T>(), this->n_basis_, 1); | ||
|
|
||
| // Line minimization in [ψ, cg] subspace (ABACUS update_psi) | ||
| converged = this->line_minimization(ppsi, | ||
| cg, | ||
| scg, | ||
| ethr_band[m], | ||
| cg_norm, | ||
| theta, | ||
| eigen_pack[m], | ||
| psi_m, | ||
| spsi_m, | ||
| hpsi_m); | ||
|
|
||
| iter++; | ||
| } while (!converged && iter < this->pw_diag_nmax_); | ||
|
|
||
| // Store converged eigenvector | ||
| psi[m].sync(psi_m); |
| item.check_value = [](const Input_Item& item, const Parameter& para) { | ||
| const std::string& ks_solver = para.input.ks_solver; | ||
| const std::vector<std::string> pw_solvers = {"cg", "dav", "bpcg", "dav_subspace"}; | ||
| const std::vector<std::string> pw_solvers = {"cg", "dav", "bpcg", "dav_subspace", "ppcg"}; |
| } | ||
| } | ||
|
|
||
| double *en = new double[npw]; |
Comment on lines
+202
to
+207
| INSTANTIATE_TEST_SUITE_P(VerifyPPCG, | ||
| DiagoPPCGTest, | ||
| ::testing::Values( | ||
| // nband, npw, sparsity, eps, maxiter, threshold | ||
| DiagoPPCGPrepare(4, 100, 0, 1e-5, 300, 5e-2) | ||
| )); |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Reminder
Linked Issue
Fix #...
Unit Tests and/or Case Tests for my changes
What's changed?
Any changes of core modules? (ignore if not applicable)