diff --git a/docs/advanced/scf/hsolver.md b/docs/advanced/scf/hsolver.md index cc85b4d87be..ba384eef217 100644 --- a/docs/advanced/scf/hsolver.md +++ b/docs/advanced/scf/hsolver.md @@ -4,7 +4,15 @@ Method of explicit solving KS-equation can be chosen by variable "ks_solver" in INPUT file. -When "basis_type = pw", `ks_solver` can be `cg`, `bpcg` or `dav`. The default setting `cg` is recommended, which is band-by-band conjugate gradient diagonalization method. There is a large probability that the use of setting of `dav` , which is block Davidson diagonalization method, can be tried to improve performance. +When "basis_type = pw", `ks_solver` can be `cg`, `bpcg`, `dav`, `dav_subspace`, or `ppcg`. The default setting `cg` is recommended, which is a band-by-band conjugate-gradient diagonalization method. The `dav` and `dav_subspace` settings use Davidson-style subspace diagonalization and can be tried to improve performance. + +The `ppcg` setting uses the projection preconditioned conjugate-gradient method, a restarted block method with bounded workspace. It is optimized and validated for CPU plane-wave calculations; non-CPU devices use a transitional host/device bridge. PPCG is unavailable for `esolver_type = sdft`, including mixed stochastic-deterministic calculations. + +PPCG uses `pw_diag_ndim` as its block size and performs the conservative full Rayleigh--Ritz correction on every iteration. Cached Hamiltonian and overlap images are refreshed independently, with additional refreshes when numerical cancellation threatens the residual. + +If PPCG reaches `pw_diag_nmax` before satisfying its residual criterion, the running log reports the active-band count and maximum residual. An SCF calculation can continue to its next electronic iteration with a warning; an NSCF calculation stops with an error because it has no subsequent SCF iteration to improve these eigenpairs. A non-finite residual is an error in either calculation. + +For diagnostics, `ABACUS_PPCG_RESIDUAL_TRACE` selects an optional CSV trace. Only the MPI world root writes this file, and its resolved full path is recorded in the running log each time the file is created or overwritten. When "basis_type = lcao", `ks_solver` can be `genelpa` or `scalapack_gvx`. The default setting `genelpa` is recommended, which is based on ELPA (EIGENVALUE SOLVERS FOR PETAFLOP APPLICATIONS) (https://elpa.mpcdf.mpg.de/) and the kernel is auto choosed by GENELPA(https://github.com/pplab/GenELPA), usually faster than the setting of "scalapack_gvx", which is based on ScaLAPACK(Scalable Linear Algebra PACKage) diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 243b1ee57f4..9e4dcaa0511 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -503,6 +503,7 @@ OBJS_HSOLVER=diago_cg.o\ diago_david.o\ diago_dav_subspace.o\ diago_bpcg.o\ + diago_ppcg.o\ para_lin_tf.o\ hsolver.o\ hsolver_pw.o\ diff --git a/source/source_base/parallel_reduce.cpp b/source/source_base/parallel_reduce.cpp index 741b38f36d5..2d67255b19d 100644 --- a/source/source_base/parallel_reduce.cpp +++ b/source/source_base/parallel_reduce.cpp @@ -13,6 +13,29 @@ const MPI_Datatype Parallel_Reduce::MPI_Type>::value = MPI_C const MPI_Datatype Parallel_Reduce::MPI_Type::value = MPI_LONG_LONG; #endif +bool Parallel_Reduce::is_world_root() +{ +#ifdef __MPI + int initialized = 0; + MPI_Initialized(&initialized); + if (!initialized) + { + return true; + } + int finalized = 0; + MPI_Finalized(&finalized); + if (finalized) + { + return false; + } + int rank = 0; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + return rank == 0; +#else + return true; +#endif +} + template void Parallel_Reduce::reduce_all(T& object) { diff --git a/source/source_base/parallel_reduce.h b/source/source_base/parallel_reduce.h index 2e2f842f3a2..9130a002104 100644 --- a/source/source_base/parallel_reduce.h +++ b/source/source_base/parallel_reduce.h @@ -11,6 +11,9 @@ namespace Parallel_Reduce { +/// Return whether this process is the MPI world root (or is running serially). +bool is_world_root(); + #ifdef __MPI template struct MPI_Type; diff --git a/source/source_estate/elecstate_print.cpp b/source/source_estate/elecstate_print.cpp index 4c4a2d9edf5..6acf3c95fa7 100644 --- a/source/source_estate/elecstate_print.cpp +++ b/source/source_estate/elecstate_print.cpp @@ -57,6 +57,7 @@ void print_scf_iterinfo(const std::string& ks_solver, {"scalapack_gvx", "GV"}, {"cusolver", "CU"}, {"bpcg", "BP"}, + {"ppcg", "PP"}, {"pexsi", "PE"}, {"cusolvermp", "CM"}, {"sdft", "CT"}}; // CT = Chebyshev Trace, for pure SDFT (nbands=0) where no H diagonalization is performed diff --git a/source/source_estate/unittests/test_elecstate_print.cpp b/source/source_estate/unittests/test_elecstate_print.cpp index d2a8e0d058f..6a956881818 100644 --- a/source/source_estate/unittests/test_elecstate_print.cpp +++ b/source/source_estate/unittests/test_elecstate_print.cpp @@ -150,7 +150,8 @@ TEST_F(ElecStatePrintTest, PrintEtot) } // iteration of different ks_solver - std::vector ks_solvers = {"cg", "lapack", "genelpa", "dav", "scalapack_gvx", "cusolver"}; + std::vector ks_solvers + = {"cg", "lapack", "genelpa", "dav", "scalapack_gvx", "cusolver", "ppcg"}; for (int i = 0; i < ks_solvers.size(); i++) { inp.ks_solver = ks_solvers[i]; @@ -184,6 +185,10 @@ TEST_F(ElecStatePrintTest, PrintEtot) { EXPECT_THAT(output, testing::HasSubstr("CU")); } + else if (inp.ks_solver == "ppcg") + { + EXPECT_THAT(output, testing::HasSubstr("PP")); + } } GlobalV::ofs_running.close(); ifs.open("test.dat", std::ios::in); diff --git a/source/source_hsolver/CMakeLists.txt b/source/source_hsolver/CMakeLists.txt index 8a4b36cc149..29496029778 100644 --- a/source/source_hsolver/CMakeLists.txt +++ b/source/source_hsolver/CMakeLists.txt @@ -6,6 +6,7 @@ list( diago_david.cpp diago_dav_subspace.cpp diago_bpcg.cpp + diago_ppcg.cpp para_lin_tf.cpp hsolver_pw.cpp hsolver_lcaopw.cpp diff --git a/source/source_hsolver/diago_params.cpp b/source/source_hsolver/diago_params.cpp index 28e1040a974..6a49fa52a2f 100644 --- a/source/source_hsolver/diago_params.cpp +++ b/source/source_hsolver/diago_params.cpp @@ -16,7 +16,8 @@ void setup_diago_params_pw(const int istep, DiagoIterAssist::SCF_ITER = iter; DiagoIterAssist::PW_DIAG_THR = ethr; - if (inp.calculation != "nscf") + // PPCG also uses the configured iteration budget for NSCF calculations. + if (inp.calculation != "nscf" || inp.ks_solver == "ppcg") { DiagoIterAssist::PW_DIAG_NMAX = inp.pw_diag_nmax; } diff --git a/source/source_hsolver/diago_ppcg.cpp b/source/source_hsolver/diago_ppcg.cpp new file mode 100644 index 00000000000..92849abb7fe --- /dev/null +++ b/source/source_hsolver/diago_ppcg.cpp @@ -0,0 +1,2415 @@ +#include "diago_ppcg.h" +#include "source_base/parallel_reduce.h" +#include "source_base/timer.h" +#include +#include +#include +#include +#include +#include "source_base/kernels/math_kernel_op.h" +#include +#include +#include +#include +#include +#include +#include + +namespace hsolver { +namespace { + +const int ppcg_openmp_work_threshold = 4096; +const int ppcg_openmp_column_threshold = 16; +const double ppcg_minimum_diagonalization_threshold = 1.0e-14; +const double ppcg_preconditioner_threshold = 1.0e-12; +const double ppcg_numerical_threshold = 1.0e-30; +const double ppcg_scaling_threshold = 1.0e-15; +const int ppcg_stagnation_restart_interval = 15; +// Operator images are transformed with their vectors. Refresh their accumulated +// roundoff independently of the user-selected Ritz diagonalization interval. +const int ppcg_operator_refresh_interval = 16; +// Relative cutoff for overlap eigenvalues in a local [X, W, P] subspace. +// Directions below this level are numerically dependent; retaining them can +// make the reduced generalized solve stagnate instead of improving a Ritz pair. +const double ppcg_subspace_rank_threshold = 1.0e-8; + +// Diagonal shifts used by the small projected eigensolve fallback. +const double ppcg_subspace_shifts[] = {0.0, 1.0e-10, 1.0e-8, 1.0e-6}; + +// A whole-function timer also closes correctly on early return or exception. +class PPCGFunctionTimer +{ +public: + explicit PPCGFunctionTimer(const char* function) : function_(function) + { + ModuleBase::timer::start("DiagoPPCG", function_); + } + + ~PPCGFunctionTimer() + { + ModuleBase::timer::end("DiagoPPCG", function_); + } + +private: + const char* function_; +}; + +} // namespace +} // namespace hsolver + + + +namespace hsolver { +namespace { + +template +void reduce_pool_values(Value& value) +{ + Parallel_Reduce::reduce_pool(value); +} + +template +void reduce_pool_values(Value* value, const int n) +{ + Parallel_Reduce::reduce_pool(value, n); +} + +bool all_pool_operations_succeeded(const bool local_success) +{ + int failure_count = local_success ? 0 : 1; + reduce_pool_values(&failure_count, 1); + return failure_count == 0; +} + +template +Real max_generalized_residual( + const T* hpsi, + const T* spsi, + const Real* eigenvalue, + int ld, + int n_dim, + int ncol) +{ + const PPCGFunctionTimer timer("max_generalized_residual"); + Real max_res = 0; + std::vector nrm2_all(ncol, 0.0); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) if (n_dim * ncol > ppcg_openmp_work_threshold) +#endif + for (int j = 0; j < ncol; ++j) + { + double nrm2 = 0.0; + for (int ig = 0; ig < n_dim; ++ig) + { + const T r = hpsi[ig + j * ld] - T(eigenvalue[j]) * spsi[ig + j * ld]; + nrm2 += double(std::norm(r)); + } + nrm2_all[j] = nrm2; + } + reduce_pool_values(nrm2_all.data(), ncol); + for (int j = 0; j < ncol; ++j) + { + const Real residual_norm = std::sqrt(Real(nrm2_all[j])); + if (!std::isfinite(residual_norm)) + { + return std::numeric_limits::infinity(); + } + max_res = std::max(max_res, residual_norm); + } + return max_res; +} + +template +inline void set_zero(std::vector& x) +{ + const int n = int(x.size()); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) if (n > ppcg_openmp_work_threshold) +#endif + for (int i = 0; i < n; ++i) + { + x[i] = T(0); + } +} + +} // anonymous namespace +} // namespace hsolver + + + +namespace hsolver { +namespace { + +template +struct HermitianLapack +{ + using Real = typename container::GetTypeReal::type; + using Device = container::DEVICE_CPU; + + static void sygvd(int n, Scalar* a, Scalar* b, Real* w) + { + std::vector eigenvectors(n * n); + container::kernels::lapack_hegvd()( + n, n, a, b, w, eigenvectors.data()); + std::copy(eigenvectors.begin(), eigenvectors.end(), a); + } +}; + +} // anonymous namespace +} // namespace hsolver + + + +namespace hsolver { +namespace { + +inline bool ppcg_contiguous_cols(const std::vector& cols, int& first) +{ + if (cols.empty()) + { + return false; + } + + first = cols.front(); + for (int j = 0; j < int(cols.size()); ++j) + { + if (cols[j] != first + j) + { + return false; + } + } + return true; +} + +} // anonymous namespace + +bool ppcg_workspace_dimensions_valid(const int ld_psi, + const int nband, + const int dim, + const int sbsize) +{ + if (ld_psi <= 0 || nband <= 0 || dim <= 0 || dim > ld_psi) + { + return false; + } + const long long limit = std::numeric_limits::max(); + const int bounded_sbsize = std::min(nband, sbsize); + const int block_size = std::max(1, bounded_sbsize); + if (static_cast(nband) > limit / ld_psi || block_size > limit / 3 || nband > limit / 3) + { + return false; + } + // The largest local basis holds [X, W, P]. Startup correction can hold + // [X, W], while discarded directions can add up to 2*nband columns during + // root recovery. The expanded recovery solve therefore reaches 3*nband. + const long long local_dim = 3LL * block_size; + const long long recovery_dim = 3LL * nband; + if (local_dim > limit / ld_psi || recovery_dim > limit / ld_psi) + { + return false; + } + long long projected_dim = std::max(local_dim, recovery_dim); + if (block_size < nband) + { + const long long startup_dim = 2LL * nband; + projected_dim = std::max(projected_dim, startup_dim); + } + const long long matrix_size = projected_dim * projected_dim; + const long long linear_workspace = 1 + 6 * projected_dim; + if (matrix_size > limit || linear_workspace > limit) + { + return false; + } + // lapack_heevd/hegvd additionally form 1+6*d+2*d*d in int arithmetic. + return matrix_size <= (limit - linear_workspace) / 2; +} + +// ============================================================================= +// Constructor +// ============================================================================= +template +DiagoPPCG::DiagoPPCG(const Real& diag_thr, + const int& diag_iter_max, + const int& sbsize, + const int& rr_step) + : maxiter_(diag_iter_max), + sbsize_(sbsize), + rr_step_(rr_step), + diag_thr_(diag_thr) +{ + sbsize_ = std::max(1, sbsize_); + rr_step_ = std::max(1, rr_step_); + diag_thr_ = std::max(diag_thr_, Real(ppcg_minimum_diagonalization_threshold)); +} + +// ============================================================================= +// Input validation +// ============================================================================= +template +void DiagoPPCG::validate_input( + const HPsiFunc& hpsi_func, + const T* psi_in, + const Real* eigenvalue_in, + const std::vector& ethr_band, + const Real* prec) const +{ + std::string local_error; + try + { + if (!hpsi_func) + { + throw std::invalid_argument("PPCG: H operator is empty."); + } + if (psi_in == nullptr || eigenvalue_in == nullptr) + { + throw std::invalid_argument("PPCG: psi/eigenvalue pointer is null."); + } + if (prec == nullptr) + { + throw std::invalid_argument("PPCG: preconditioner pointer is null."); + } + if (ld_psi_ <= 0 || n_band_ <= 0 || n_dim_ <= 0) + { + throw std::invalid_argument("PPCG: invalid dimensions."); + } + if (n_dim_ > ld_psi_) + { + throw std::invalid_argument("PPCG: dim must not exceed ld_psi."); + } + if (!ppcg_workspace_dimensions_valid(ld_psi_, n_band_, n_dim_, sbsize_)) + { + throw std::invalid_argument("PPCG: workspace dimensions exceed integer indexing or LAPACK limits."); + } + if (ethr_band.size() < size_t(n_band_)) + { + throw std::invalid_argument("PPCG: ethr_band size is smaller than nband."); + } + for (int i = 0; i < n_band_; ++i) + { + if (!std::isfinite(ethr_band[i])) + { + throw std::invalid_argument("PPCG: ethr_band contains non-finite value."); + } + } + for (int i = 0; i < n_dim_; ++i) + { + if (!std::isfinite(prec[i])) + { + throw std::invalid_argument("PPCG: preconditioner contains non-finite value."); + } + } + } + catch (const std::invalid_argument& error) + { + local_error = error.what(); + } + // All pool ranks must reject the call before any rank enters an operator + // callback or Gram-matrix collective, even if only one input is invalid. + const bool local_success = local_error.empty(); + if (!all_pool_operations_succeeded(local_success)) + { + if (local_error.empty()) + { + local_error = "PPCG: invalid input on another rank in the MPI pool."; + } + throw std::invalid_argument(local_error); + } +} + +// ============================================================================= +// Operator application +// ============================================================================= +template +void DiagoPPCG::apply_h(const HPsiFunc& hpsi_func, + T* psi_in, T* hpsi_out, + int ncol) const +{ + const PPCGFunctionTimer timer("apply_h"); + hpsi_func(psi_in, hpsi_out, ld_psi_, ncol); +} + +template +void DiagoPPCG::apply_s(const SPsiFunc& spsi_func, + T* psi_in, T* spsi_out, + int ncol) const +{ + const PPCGFunctionTimer timer("apply_s"); + if (spsi_func) + { + spsi_func(psi_in, spsi_out, ld_psi_, ncol); + } + else + { +#ifdef _OPENMP +#pragma omp parallel for schedule(static) if (ld_psi_ * ncol > ppcg_openmp_work_threshold) +#endif + for (int j = 0; j < ncol; ++j) + { + std::copy(psi_in + j * ld_psi_, psi_in + (j + 1) * ld_psi_, + spsi_out + j * ld_psi_); + } + } +} + +template +void DiagoPPCG::apply_s_current(T* psi_in, T* spsi_out, + int ncol) const +{ + apply_s(spsi_func_, psi_in, spsi_out, ncol); +} + +// ============================================================================= +// Inner product (real part only, for Hermitian operators) +// ============================================================================= +template +typename DiagoPPCG::Real +DiagoPPCG::gamma_dot(const T* x, const T* y) const +{ + Real result = ModuleBase::dot_real_op()(n_dim_, x, y, false); + reduce_pool_values(result); + return result; +} + +// ============================================================================= +// Gram matrix: out[i, j] = +// ============================================================================= +template +void DiagoPPCG::gram(const T* mat_a, const T* mat_b, + int ncol_a, int ncol_b, + std::vector& out, + int ld_out) const +{ + const int reduction_count = ld_out * ncol_b; + out.resize(reduction_count); + const T one = T(1); + const T zero = T(0); + ModuleBase::gemm_op()('C', + 'N', + ncol_a, + ncol_b, + n_dim_, + &one, + mat_a, + ld_psi_, + mat_b, + ld_psi_, + &zero, + out.data(), + ld_out); + reduce_pool_values(out.data(), reduction_count); +} + +// ============================================================================= +// Column gather: extract selected columns into contiguous storage +// ============================================================================= +template +void DiagoPPCG::copy_cols(const T* src, + const std::vector& cols, + std::vector& dst) const +{ + const int ncols = int(cols.size()); + dst.resize(ld_psi_ * ncols); + if (ncols == 0) + { + return; + } + + int first = 0; + if (ppcg_contiguous_cols(cols, first)) + { + std::copy(src + first * ld_psi_, + src + (first + ncols) * ld_psi_, + dst.begin()); + return; + } + +#ifdef _OPENMP +#pragma omp parallel for schedule(static) if (ld_psi_ * ncols > ppcg_openmp_work_threshold) +#endif + for (int j = 0; j < ncols; ++j) + { + const int c = cols[j]; + std::copy(src + c * ld_psi_, src + c * ld_psi_ + ld_psi_, + dst.begin() + j * ld_psi_); + } +} + +// ============================================================================= +// Column scatter: write contiguous storage back into selected columns +// ============================================================================= +template +void DiagoPPCG::scatter_cols( + T* dst, + const std::vector& cols, + const std::vector& src) const +{ + const int ncols = int(cols.size()); + if (ncols == 0) + { + return; + } + + int first = 0; + if (ppcg_contiguous_cols(cols, first)) + { + std::copy(src.begin(), + src.begin() + ld_psi_ * ncols, + dst + first * ld_psi_); + return; + } + +#ifdef _OPENMP +#pragma omp parallel for schedule(static) if (ld_psi_ * ncols > ppcg_openmp_work_threshold) +#endif + for (int j = 0; j < ncols; ++j) + { + const int c = cols[j]; + std::copy(src.begin() + j * ld_psi_, + src.begin() + (j + 1) * ld_psi_, + dst + c * ld_psi_); + } +} + +// ============================================================================= +// Project x onto vectors orthogonal to S-orthonormal basis +// ============================================================================= +template +void DiagoPPCG::project_against( + const T* basis, const T* sbasis, const T* hbasis, + const std::vector& basis_cols, + std::vector& x, std::vector& sx, std::vector* hx, + const std::vector& x_cols) const +{ + const PPCGFunctionTimer timer("project_against"); + if (basis_cols.empty() || x_cols.empty()) + { + return; + } + + const int nbasis = int(basis_cols.size()); + const int nx = int(x_cols.size()); + + int x_first = 0; + const bool contiguous_x = ppcg_contiguous_cols(x_cols, x_first); + + std::vector x_l; + std::vector sx_l; + std::vector hx_l; + T* x_data = x.data() + x_first * ld_psi_; + T* sx_data = sx.data() + x_first * ld_psi_; + T* hx_data = hx == nullptr ? nullptr : hx->data() + x_first * ld_psi_; + if (!contiguous_x) + { + x_l.reserve(ld_psi_ * nx); + sx_l.reserve(ld_psi_ * nx); + copy_cols(x.data(), x_cols, x_l); + copy_cols(sx.data(), x_cols, sx_l); + x_data = x_l.data(); + sx_data = sx_l.data(); + if (hx != nullptr) + { + copy_cols(hx->data(), x_cols, hx_l); + hx_data = hx_l.data(); + } + } + + int basis_first = 0; + const bool contiguous_basis = + ppcg_contiguous_cols(basis_cols, basis_first); + + std::vector basis_l; + std::vector sbasis_l; + std::vector hbasis_l; + const T* basis_data = basis + basis_first * ld_psi_; + const T* sbasis_data = sbasis + basis_first * ld_psi_; + const T* hbasis_data = hbasis == nullptr ? nullptr : hbasis + basis_first * ld_psi_; + if (!contiguous_basis) + { + basis_l.reserve(ld_psi_ * nbasis); + sbasis_l.reserve(ld_psi_ * nbasis); + copy_cols(basis, basis_cols, basis_l); + copy_cols(sbasis, basis_cols, sbasis_l); + basis_data = basis_l.data(); + sbasis_data = sbasis_l.data(); + if (hx != nullptr) + { + copy_cols(hbasis, basis_cols, hbasis_l); + hbasis_data = hbasis_l.data(); + } + } + + std::vector coeff(nbasis * nx, T(0)); + gram(basis_data, sx_data, nbasis, nx, coeff, nbasis); + + const T minus_one = T(-1); + const T one = T(1); + ModuleBase::gemm_op()('N', + 'N', + n_dim_, + nx, + nbasis, + &minus_one, + basis_data, + ld_psi_, + coeff.data(), + nbasis, + &one, + x_data, + ld_psi_); + ModuleBase::gemm_op()('N', + 'N', + n_dim_, + nx, + nbasis, + &minus_one, + sbasis_data, + ld_psi_, + coeff.data(), + nbasis, + &one, + sx_data, + ld_psi_); + if (hx != nullptr) + { + // H(x-X*C) = Hx-HX*C; the same coefficients update all three images. + ModuleBase::gemm_op()('N', + 'N', + n_dim_, + nx, + nbasis, + &minus_one, + hbasis_data, + ld_psi_, + coeff.data(), + nbasis, + &one, + hx_data, + ld_psi_); + } + + if (!contiguous_x) + { + scatter_cols(x.data(), x_cols, x_l); + scatter_cols(sx.data(), x_cols, sx_l); + if (hx != nullptr) + { + scatter_cols(hx->data(), x_cols, hx_l); + } + } +} + +// ============================================================================= +// Project P and refresh its operator images when subtraction loses precision. +// ============================================================================= +template +void DiagoPPCG::project_search_directions( + const HPsiFunc& hpsi_func, + const T* psi, + const std::vector& all_cols, + const std::vector& active_cols) +{ + const PPCGFunctionTimer timer("project_search_directions"); + const int nactive = int(active_cols.size()); + std::vector norms(2 * nactive, 0.0); + for (int j = 0; j < nactive; ++j) + { + const int column = active_cols[j]; + for (int i = 0; i < n_dim_; ++i) + { + norms[j] += std::norm(p_[idx(i, column, ld_psi_)]); + } + } + project_against(psi, spsi_.data(), hpsi_.data(), all_cols, p_, sp_, &hp_, active_cols); + for (int j = 0; j < nactive; ++j) + { + const int column = active_cols[j]; + for (int i = 0; i < n_dim_; ++i) + { + norms[nactive + j] += std::norm(p_[idx(i, column, ld_psi_)]); + } + } + const int norm_count = int(norms.size()); + reduce_pool_values(norms.data(), norm_count); + std::vector refresh_cols; + for (int j = 0; j < nactive; ++j) + { + if (norms[j] == 0.0 || norms[nactive + j] == 0.0) + { + const int offset = active_cols[j] * ld_psi_; + std::fill(p_.begin() + offset, p_.begin() + offset + ld_psi_, T(0)); + std::fill(sp_.begin() + offset, sp_.begin() + offset + ld_psi_, T(0)); + std::fill(hp_.begin() + offset, hp_.begin() + offset + ld_psi_, T(0)); + continue; + } + // A strongly reduced P suffers cancellation in HP-HX*C and SP-SX*C. + // Normalizing that tiny P would amplify the absolute image error. + // Refresh only these columns; reduction gives every rank one selection. + if (norms[nactive + j] < 0.5 * norms[j]) + { + refresh_cols.push_back(active_cols[j]); + } + } + if (!refresh_cols.empty()) + { + std::vector projected; + copy_cols(p_.data(), refresh_cols, projected); + std::vector hprojected(projected.size(), T(0)); + std::vector sprojected(projected.size(), T(0)); + const int refresh_count = int(refresh_cols.size()); + apply_h(hpsi_func, projected.data(), hprojected.data(), refresh_count); + apply_s_current(projected.data(), sprojected.data(), refresh_count); + scatter_cols(hp_.data(), refresh_cols, hprojected); + scatter_cols(sp_.data(), refresh_cols, sprojected); + } +} + +// ============================================================================= +// Preconditioner: x[c] /= max(prec, eps) for each active column c +// ============================================================================= +template +void DiagoPPCG::divide_by_preconditioner( + const std::vector& active_cols, + const Real* prec, + std::vector& x) const +{ + const int ncols = int(active_cols.size()); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) if (n_dim_ * ncols > ppcg_openmp_work_threshold) +#endif + for (int j = 0; j < ncols; ++j) + { + const int c = active_cols[j]; + for (int ig = 0; ig < n_dim_; ++ig) + { + x[idx(ig, c, ld_psi_)] /= + std::max(prec[ig], Real(ppcg_preconditioner_threshold)); + } + } +} + +} // namespace hsolver + + +namespace hsolver { + +//============================================================================== +// BLOCK_SUBSPACE STRATEGY +//============================================================================== + +// --------------------------------------------------------------------------- +// Lock only a contiguous prefix of converged eigenpairs. A small residual +// proves that a vector is an eigenpair, but it does not prove that it belongs +// to the requested lowest roots. Keeping every band after the first +// unconverged root active lets an evolving lower root replace a later exact +// high-energy state. +// --------------------------------------------------------------------------- +template +void DiagoPPCG::lock_epairs( + const std::vector& residual, + const std::vector& ethr_band, + std::vector& active_cols) +{ + active_cols.clear(); + active_cols.reserve(n_band_); + std::vector nrm2_all(n_band_, 0.0); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) if (n_dim_ * n_band_ > ppcg_openmp_work_threshold) +#endif + for (int j = 0; j < n_band_; ++j) + { + double nrm2 = 0.0; + for (int ig = 0; ig < n_dim_; ++ig) + { + nrm2 += double(std::norm(residual[idx(ig, j, ld_psi_)])); + } + nrm2_all[j] = nrm2; + } + reduce_pool_values(nrm2_all.data(), n_band_); + bool found_unconverged = false; + max_residual_ = Real(0); + for (int j = 0; j < n_band_; ++j) + { + const Real rnrm = std::sqrt(std::max(Real(nrm2_all[j]), Real(0))); + const Real thr = std::max(Real(ethr_band[j]), diag_thr_); + if (std::isfinite(rnrm)) + { + max_residual_ = std::max(max_residual_, rnrm); + } + else + { + max_residual_ = std::numeric_limits::infinity(); + } + if (!std::isfinite(rnrm) || rnrm > thr) + { + found_unconverged = true; + } + if (found_unconverged) + { + active_cols.push_back(j); + } + } +} + +template +void DiagoPPCG::scale_to_unit_snorm(std::vector& x, + std::vector& sx, + std::vector& hx, + const int ncols) const +{ + std::vector inverse_snorm(ncols, 0.0); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) if (n_dim_ * ncols > ppcg_openmp_work_threshold) +#endif + for (int j = 0; j < ncols; ++j) + { + double snorm_squared = 0.0; + for (int ig = 0; ig < n_dim_; ++ig) + { + snorm_squared += double(std::real( + std::conj(x[idx(ig, j, ld_psi_)]) + * sx[idx(ig, j, ld_psi_)])); + } + inverse_snorm[j] = snorm_squared; + } + reduce_pool_values(inverse_snorm.data(), ncols); + for (int j = 0; j < ncols; ++j) + { + const Real snorm = std::sqrt(std::max( + Real(inverse_snorm[j]), Real(ppcg_numerical_threshold))); + inverse_snorm[j] = snorm > Real(ppcg_scaling_threshold) + ? double(Real(1) / snorm) + : 1.0; + } +#ifdef _OPENMP +#pragma omp parallel for collapse(2) schedule(static) if (n_dim_ * ncols > ppcg_openmp_work_threshold) +#endif + for (int j = 0; j < ncols; ++j) + { + for (int ig = 0; ig < n_dim_; ++ig) + { + const Real scale = Real(inverse_snorm[j]); + x[idx(ig, j, ld_psi_)] *= scale; + sx[idx(ig, j, ld_psi_)] *= scale; + hx[idx(ig, j, ld_psi_)] *= scale; + } + } +} + +template +void DiagoPPCG::hermitize_projected(std::vector& matrix, + const int dim) const +{ + for (int j = 0; j < dim; ++j) + { + matrix[j + j * dim] = T(std::real(matrix[j + j * dim]), 0); + for (int i = j + 1; i < dim; ++i) + { + const T value = (matrix[i + j * dim] + + std::conj(matrix[j + i * dim])) + * Real(0.5); + matrix[i + j * dim] = value; + matrix[j + i * dim] = std::conj(value); + } + } +} + +template +void DiagoPPCG::insert_startup_gram_block( + const T* left, + const T* right, + const int row_count, + const int column_count, + const int row_offset, + const int column_offset, + const int matrix_dim, + std::vector& matrix, + std::vector& workspace) const +{ + gram(left, right, row_count, column_count, workspace, row_count); + for (int j = 0; j < column_count; ++j) + { + for (int i = 0; i < row_count; ++i) + { + const int matrix_index = row_offset + i + + (column_offset + j) * matrix_dim; + matrix[matrix_index] = workspace[i + j * row_count]; + } + } +} + +template +void DiagoPPCG::build_startup_global_subspace( + const T* psi, + std::vector& w_active, + std::vector& sw_active, + std::vector& hw_active, + const int active_count, + SmallSubspace& subspace) +{ + scale_to_unit_snorm(w_active, sw_active, hw_active, active_count); + + const int state_count = n_band_; + const int correction_offset = state_count; + const int local_dim = state_count + active_count; + const int matrix_size = local_dim * local_dim; + subspace.k.assign(matrix_size, T(0)); + subspace.m.assign(matrix_size, T(0)); + subspace.eval.resize(local_dim); + std::vector gram_workspace; + + insert_startup_gram_block(psi, + hpsi_.data(), + state_count, + state_count, + 0, + 0, + local_dim, + subspace.k, + gram_workspace); + insert_startup_gram_block(psi, + hw_active.data(), + state_count, + active_count, + 0, + correction_offset, + local_dim, + subspace.k, + gram_workspace); + insert_startup_gram_block(w_active.data(), + hpsi_.data(), + active_count, + state_count, + correction_offset, + 0, + local_dim, + subspace.k, + gram_workspace); + insert_startup_gram_block(w_active.data(), + hw_active.data(), + active_count, + active_count, + correction_offset, + correction_offset, + local_dim, + subspace.k, + gram_workspace); + insert_startup_gram_block(psi, + spsi_.data(), + state_count, + state_count, + 0, + 0, + local_dim, + subspace.m, + gram_workspace); + insert_startup_gram_block(psi, + sw_active.data(), + state_count, + active_count, + 0, + correction_offset, + local_dim, + subspace.m, + gram_workspace); + insert_startup_gram_block(w_active.data(), + spsi_.data(), + active_count, + state_count, + correction_offset, + 0, + local_dim, + subspace.m, + gram_workspace); + insert_startup_gram_block(w_active.data(), + sw_active.data(), + active_count, + active_count, + correction_offset, + correction_offset, + local_dim, + subspace.m, + gram_workspace); + + hermitize_projected(subspace.k, local_dim); + hermitize_projected(subspace.m, local_dim); +} + +template +void DiagoPPCG::update_startup_global_subspace( + T* psi, + const std::vector& w_active, + const std::vector& sw_active, + const std::vector& hw_active, + const int active_count, + const SmallSubspace& subspace) +{ + const int local_dim = n_band_ + active_count; + const T* state_coefficients = subspace.k.data(); + const T* correction_coefficients = state_coefficients + n_band_; + const T one = T(1); + const T zero = T(0); + + set_zero(w_); + set_zero(sw_); + set_zero(hw_); + set_zero(p_); + set_zero(sp_); + set_zero(hp_); + + auto combine_state = [&](const T* state, + const T* correction, + T* output) + { + ModuleBase::gemm_op()('N', + 'N', + n_dim_, + n_band_, + n_band_, + &one, + state, + ld_psi_, + state_coefficients, + local_dim, + &zero, + output, + ld_psi_); + ModuleBase::gemm_op()('N', + 'N', + n_dim_, + n_band_, + active_count, + &one, + correction, + ld_psi_, + correction_coefficients, + local_dim, + &one, + output, + ld_psi_); + }; + auto combine_direction = [&](const T* correction, T* output) + { + ModuleBase::gemm_op()('N', + 'N', + n_dim_, + n_band_, + active_count, + &one, + correction, + ld_psi_, + correction_coefficients, + local_dim, + &zero, + output, + ld_psi_); + }; + + combine_state(psi, w_active.data(), w_.data()); + combine_state(spsi_.data(), sw_active.data(), sw_.data()); + combine_state(hpsi_.data(), hw_active.data(), hw_.data()); + combine_direction(w_active.data(), p_.data()); + combine_direction(sw_active.data(), sp_.data()); + combine_direction(hw_active.data(), hp_.data()); + + std::copy(w_.begin(), w_.end(), psi); + std::copy(sw_.begin(), sw_.end(), spsi_.begin()); + std::copy(hw_.begin(), hw_.end(), hpsi_.begin()); +} + +// --------------------------------------------------------------------------- +// Build K = V^H H V and M = V^H S V where V = [psi, w] +// --------------------------------------------------------------------------- +template +void DiagoPPCG::build_small_subspace( + const T* psi, + const std::vector& cols, + int nblk, + SmallSubspace& subspace) const +{ + const PPCGFunctionTimer timer("build_small_subspace"); + const int l = int(cols.size()); + const int dim = nblk * l; + subspace.k.resize(dim * dim); + subspace.m.resize(dim * dim); + subspace.eval.resize(dim); + + copy_cols(psi, cols, subspace.psi_l); + copy_cols(spsi_.data(), cols, subspace.spsi_l); + copy_cols(hpsi_.data(), cols, subspace.hpsi_l); + copy_cols(w_.data(), cols, subspace.w_l); + copy_cols(sw_.data(), cols, subspace.sw_l); + copy_cols(hw_.data(), cols, subspace.hw_l); + if (nblk >= 3) + { + copy_cols(p_.data(), cols, subspace.p_l); + copy_cols(sp_.data(), cols, subspace.sp_l); + copy_cols(hp_.data(), cols, subspace.hp_l); + } + + // --------------------------------------------------------------------------- + // Normalize w columns to unit S-norm for numerical stability. + // + // The w block of the Gram matrix M has entries O(||w||^2) which become + // tiny when residuals are small, making M nearly singular and causing + // sygvd to produce garbage eigenvectors. + // + // Scaling to unit S-norm keeps M well-conditioned (diagonal ~1) without + // changing the subspace. The same scaled basis is reused in update_one_block. + // --------------------------------------------------------------------------- + scale_to_unit_snorm(subspace.w_l, + subspace.sw_l, + subspace.hw_l, + l); + if (nblk >= 3) + { + scale_to_unit_snorm(subspace.p_l, + subspace.sp_l, + subspace.hp_l, + l); + } + + auto copy_block = [&](const std::vector& src, + const int col0, + std::vector& dst) + { +#ifdef _OPENMP +#pragma omp parallel for schedule(static) if (ld_psi_ * l > ppcg_openmp_work_threshold) +#endif + for (int j = 0; j < l; ++j) + { + const auto source_begin = src.begin() + j * ld_psi_; + const auto source_end = source_begin + ld_psi_; + const auto destination_begin = dst.begin() + (col0 + j) * ld_psi_; + std::copy(source_begin, source_end, destination_begin); + } + }; + + const int basis_size = ld_psi_ * dim; + subspace.basis.resize(basis_size); + subspace.hbasis.resize(basis_size); + subspace.sbasis.resize(basis_size); + copy_block(subspace.psi_l, 0, subspace.basis); + copy_block(subspace.hpsi_l, 0, subspace.hbasis); + copy_block(subspace.spsi_l, 0, subspace.sbasis); + copy_block(subspace.w_l, l, subspace.basis); + copy_block(subspace.hw_l, l, subspace.hbasis); + copy_block(subspace.sw_l, l, subspace.sbasis); + if (nblk >= 3) + { + const int p_block_offset = 2 * l; + copy_block(subspace.p_l, p_block_offset, subspace.basis); + copy_block(subspace.hp_l, p_block_offset, subspace.hbasis); + copy_block(subspace.sp_l, p_block_offset, subspace.sbasis); + } + + gram(subspace.basis.data(), subspace.hbasis.data(), dim, dim, subspace.k, dim); + gram(subspace.basis.data(), subspace.sbasis.data(), dim, dim, subspace.m, dim); + hermitize_projected(subspace.k, dim); + hermitize_projected(subspace.m, dim); +} + +// --------------------------------------------------------------------------- +// Solve K v = λ M v (small generalized eigenvalue problem) +// --------------------------------------------------------------------------- +template +int DiagoPPCG::solve_small_generalized( + const int dim, + const int nstates, + SmallSubspace& subspace) const +{ + const PPCGFunctionTimer timer("solve_small_generalized"); + const std::vector k0 = subspace.k; + const std::vector m0 = subspace.m; + + // A nearly dependent local basis makes the generalized solve unstable + // even when LAPACK accepts M as positive definite. Work in the retained + // eigenspace of M and map the Ritz vectors back to the original basis. + std::vector overlap_eigenvectors = m0; + std::vector overlap_eigenvalues(dim); + bool local_rank_truncation_succeeded = false; + int retained_count = 0; + try + { + container::kernels::lapack_heevd()( + dim, + overlap_eigenvectors.data(), + dim, + overlap_eigenvalues.data()); + + const Real largest_overlap_eigenvalue = std::max(overlap_eigenvalues.back(), Real(0)); + const Real precision_floor = Real(10) * std::numeric_limits::epsilon(); + const Real relative_rank_threshold = std::max( + Real(ppcg_subspace_rank_threshold), + precision_floor); + const Real overlap_cutoff = largest_overlap_eigenvalue + * relative_rank_threshold; + std::vector retained_indices; + retained_indices.reserve(dim); + for (int i = 0; i < dim; ++i) + { + if (overlap_eigenvalues[i] > overlap_cutoff) + { + retained_indices.push_back(i); + } + } + + const int retained_dim = int(retained_indices.size()); + if (retained_dim >= nstates && retained_dim < dim) + { + std::vector retained_basis(dim * retained_dim, T(0)); + for (int j = 0; j < retained_dim; ++j) + { + const int source_index = retained_indices[j]; + const Real scale = Real(1) + / std::sqrt(overlap_eigenvalues[source_index]); + for (int row = 0; row < dim; ++row) + { + retained_basis[row + j * dim] = overlap_eigenvectors[row + source_index * dim] + * scale; + } + } + + // Q = U_r Lambda_r^(-1/2) orthonormalizes the retained overlap + // subspace. Use BLAS for Q^H K Q rather than a four-index + // contraction: the latter becomes prohibitively expensive when + // the first global correction contains many nearly dependent + // residual directions. + const T one = T(1); + const T zero = T(0); + std::vector projected_hamiltonian(dim * retained_dim, T(0)); + ModuleBase::gemm_op()('N', + 'N', + dim, + retained_dim, + dim, + &one, + k0.data(), + dim, + retained_basis.data(), + dim, + &zero, + projected_hamiltonian.data(), + dim); + std::vector reduced_hamiltonian(retained_dim * retained_dim, T(0)); + ModuleBase::gemm_op()('C', + 'N', + retained_dim, + retained_dim, + dim, + &one, + retained_basis.data(), + dim, + projected_hamiltonian.data(), + dim, + &zero, + reduced_hamiltonian.data(), + retained_dim); + for (int j = 0; j < retained_dim; ++j) + { + const Real diagonal_value = std::real( + reduced_hamiltonian[j + j * retained_dim]); + reduced_hamiltonian[j + j * retained_dim] = T(diagonal_value, Real(0)); + for (int i = 0; i < j; ++i) + { + const T value = Real(0.5) * ( + reduced_hamiltonian[i + j * retained_dim] + + std::conj(reduced_hamiltonian[j + i * retained_dim])); + reduced_hamiltonian[i + j * retained_dim] = value; + reduced_hamiltonian[j + i * retained_dim] = std::conj(value); + } + } + + std::vector reduced_eigenvalues(retained_dim); + container::kernels::lapack_heevd()( + retained_dim, + reduced_hamiltonian.data(), + retained_dim, + reduced_eigenvalues.data()); + + std::fill(subspace.k.begin(), subspace.k.end(), T(0)); + for (int state = 0; state < retained_dim; ++state) + { + subspace.eval[state] = reduced_eigenvalues[state]; + } + ModuleBase::gemm_op()('N', + 'N', + dim, + retained_dim, + retained_dim, + &one, + retained_basis.data(), + dim, + reduced_hamiltonian.data(), + retained_dim, + &zero, + subspace.k.data(), + dim); + local_rank_truncation_succeeded = true; + retained_count = retained_dim; + } + } + catch (const std::runtime_error&) + { + // The shifted generalized solve below remains the fallback. + } + if (all_pool_operations_succeeded(local_rank_truncation_succeeded)) + { + return retained_count; + } + + // Try with increasing diagonal shifts; fall back to identity (no update) + // if the subspace is too ill-conditioned. sygvd modifies both matrices + // in-place before it may fail. + const Real shifts[] = {Real(ppcg_subspace_shifts[0]), + Real(ppcg_subspace_shifts[1]), + Real(ppcg_subspace_shifts[2]), + Real(ppcg_subspace_shifts[3])}; + for (const Real shift : shifts) + { + subspace.k = k0; + subspace.m = m0; + for (int i = 0; i < dim; ++i) + { + subspace.m[i + i * dim] += T(shift); + } + + bool local_sygvd_succeeded = false; + try + { + HermitianLapack::sygvd(dim, subspace.k.data(), + subspace.m.data(), + subspace.eval.data()); + local_sygvd_succeeded = true; + } + catch (const std::runtime_error&) + { + // Try the next diagonal shift. + } + if (all_pool_operations_succeeded(local_sygvd_succeeded)) + { + return dim; + } + } + // All attempts failed — set eigenvectors to identity (no update). + std::fill(subspace.k.begin(), subspace.k.end(), T(0)); + for (int i = 0; i < dim; ++i) + { + subspace.k[i + i * dim] = T(1); + subspace.eval[i] = Real(std::real(k0[i + i * dim])) + / std::max(Real(std::real(m0[i + i * dim])), + Real(ppcg_numerical_threshold)); + } + return 0; +} + +// Keep valid local Ritz directions that a block-size restriction would discard. +template +void DiagoPPCG::collect_discarded_ritz( + const SmallSubspace& subspace, + const int kept_count, + const int valid_count, + OperatorBlock& candidates) const +{ + const PPCGFunctionTimer timer("collect_discarded_ritz"); + const int count = valid_count - kept_count; + if (count <= 0) + { + return; + } + const int dim = int(subspace.eval.size()); + const int offset = int(candidates.vectors.size()); + const T* coefficients = subspace.k.data() + kept_count * dim; + const T one(1); + const T zero(0); + const auto append = [&](const std::vector& basis, std::vector& output) { + output.resize(offset + ld_psi_ * count, T(0)); + T* destination = output.data() + offset; + ModuleBase::gemm_op()('N', 'N', n_dim_, count, dim, + &one, basis.data(), ld_psi_, coefficients, dim, + &zero, destination, ld_psi_); + }; + append(subspace.basis, candidates.vectors); + append(subspace.hbasis, candidates.hvectors); + append(subspace.sbasis, candidates.svectors); +} + +// Preserve X exactly as a subspace: rank decisions apply only to its complement. +template +void DiagoPPCG::project_candidate_complement( + const HPsiFunc& hpsi_func, const T* psi, OperatorBlock& candidates) const +{ + const PPCGFunctionTimer timer("project_candidate_complement"); + const int count = int(candidates.vectors.size()) / ld_psi_; + std::vector basis_cols(n_band_); + std::iota(basis_cols.begin(), basis_cols.end(), 0); + std::vector candidate_cols(count); + std::iota(candidate_cols.begin(), candidate_cols.end(), 0); + std::vector norms(2 * count, 0.0); + for (int column = 0; column < count; ++column) + { + for (int row = 0; row < n_dim_; ++row) + { + norms[column] += std::norm(candidates.vectors[idx(row, column, ld_psi_)]); + } + } + // Two passes suppress overlap left by cancellation against the retained X. + for (int pass = 0; pass < 2; ++pass) + { + project_against(psi, spsi_.data(), hpsi_.data(), basis_cols, + candidates.vectors, candidates.svectors, &candidates.hvectors, candidate_cols); + } + for (int column = 0; column < count; ++column) + { + for (int row = 0; row < n_dim_; ++row) + { + norms[count + column] += std::norm(candidates.vectors[idx(row, column, ld_psi_)]); + } + } + const int norm_count = int(norms.size()); + reduce_pool_values(norms.data(), norm_count); + const Real precision_floor = Real(10) * std::numeric_limits::epsilon(); + const Real rank_threshold = std::max(Real(ppcg_subspace_rank_threshold), precision_floor); + std::vector retained; + std::vector refresh; + for (int column = 0; column < count; ++column) + { + const double before = norms[column]; + const double after = norms[count + column]; + if (!std::isfinite(before) || !std::isfinite(after) || after <= rank_threshold * before) + { + continue; + } + retained.push_back(column); + if (after < 0.5 * before) + { + refresh.push_back(column); + } + } + // The GPU bridge owns only nband columns, even if the complement is larger. + for (int first = 0; first < int(refresh.size()); first += n_band_) + { + const int last = std::min(first + n_band_, int(refresh.size())); + const std::vector columns(refresh.begin() + first, refresh.begin() + last); + std::vector vectors; + copy_cols(candidates.vectors.data(), columns, vectors); + std::vector hvectors(vectors.size(), T(0)); + std::vector svectors(vectors.size(), T(0)); + const int batch_count = int(columns.size()); + apply_h(hpsi_func, vectors.data(), hvectors.data(), batch_count); + apply_s_current(vectors.data(), svectors.data(), batch_count); + scatter_cols(candidates.hvectors.data(), columns, hvectors); + scatter_cols(candidates.svectors.data(), columns, svectors); + } + OperatorBlock projected; + copy_cols(candidates.vectors.data(), retained, projected.vectors); + copy_cols(candidates.hvectors.data(), retained, projected.hvectors); + copy_cols(candidates.svectors.data(), retained, projected.svectors); + candidates = std::move(projected); +} + +template +bool DiagoPPCG::recover_competing_ritz( + const HPsiFunc& hpsi_func, T* psi, Real* eigenvalue, + const bool ordered_ritz, std::vector& active_cols, + const std::vector& ethr_band, OperatorBlock& candidates) +{ + const PPCGFunctionTimer timer("recover_competing_ritz"); + if (candidates.vectors.empty()) + { + return false; + } + Real upper_bound = eigenvalue[n_band_ - 1]; + if (!ordered_ritz) + { + // Between Ritz steps a column Rayleigh quotient is not an upper bound + // on the largest target Ritz value. Use the Hermitian row-sum bound. + upper_bound = std::numeric_limits::lowest(); + for (int row = 0; row < n_band_; ++row) + { + Real row_bound = std::real(rr_hsub_[row + row * n_band_]); + for (int column = 0; column < n_band_; ++column) + { + if (column != row) + { + row_bound += std::abs(rr_hsub_[row + column * n_band_]); + } + } + upper_bound = std::max(upper_bound, row_bound); + } + } + project_candidate_complement(hpsi_func, psi, candidates); + const int count = int(candidates.vectors.size()) / ld_psi_; + if (count == 0) + { + return false; + } + std::vector moments(2 * count, 0.0); + for (int column = 0; column < count; ++column) + { + for (int row = 0; row < n_dim_; ++row) + { + const int offset = idx(row, column, ld_psi_); + const T conjugate = std::conj(candidates.vectors[offset]); + moments[column] += std::real(conjugate * candidates.svectors[offset]); + moments[count + column] += std::real(conjugate * candidates.hvectors[offset]); + } + } + const int moment_count = int(moments.size()); + reduce_pool_values(moments.data(), moment_count); + const Real scale = std::max(Real(1), std::abs(upper_bound)); + const Real margin = Real(32) * std::numeric_limits::epsilon() * scale; + std::vector competing; + for (int column = 0; column < count; ++column) + { + const double norm = moments[column]; + const double energy = moments[count + column]; + if (!std::isfinite(norm) || !std::isfinite(energy) || norm <= 0.0) + { + continue; + } + // Use the physical quotient, not an eigenvalue of the shifted-M + // fallback; they need not represent the same projected problem. + if (energy / norm < upper_bound - margin) + { + competing.push_back(column); + } + } + if (competing.empty()) + { + return false; + } + OperatorBlock selected; + copy_cols(candidates.vectors.data(), competing, selected.vectors); + copy_cols(candidates.hvectors.data(), competing, selected.hvectors); + copy_cols(candidates.svectors.data(), competing, selected.svectors); + // The selected copies own the surviving directions; release the larger + // candidate buffers before constructing their orthogonal basis. + candidates = OperatorBlock{}; + const int selected_count = int(competing.size()); + scale_to_unit_snorm(selected.vectors, selected.svectors, selected.hvectors, selected_count); + std::vector overlap; + gram(selected.vectors.data(), selected.svectors.data(), selected_count, selected_count, + overlap, selected_count); + hermitize_projected(overlap, selected_count); + std::vector norms(selected_count); + bool local_success = false; + try + { + container::kernels::lapack_heevd()( + selected_count, overlap.data(), selected_count, norms.data()); + local_success = true; + } + catch (const std::runtime_error&) + { + } + if (!all_pool_operations_succeeded(local_success)) + { + throw std::runtime_error("PPCG: competing-direction overlap diagonalization failed."); + } + const Real precision_floor = Real(10) * std::numeric_limits::epsilon(); + const Real rank_threshold = std::max(Real(ppcg_subspace_rank_threshold), precision_floor); + const Real cutoff = norms.back() * rank_threshold; + int retained_count = 0; + for (int column = 0; column < selected_count; ++column) + { + if (norms[column] <= cutoff) + { + continue; + } + const Real inverse_norm = Real(1) / std::sqrt(norms[column]); + for (int row = 0; row < selected_count; ++row) + { + overlap[row + retained_count * selected_count] = overlap[row + column * selected_count] * inverse_norm; + } + ++retained_count; + } + if (retained_count == 0) + { + throw std::runtime_error("PPCG: competing-direction complement has zero numerical rank."); + } + OperatorBlock orthogonal; + const T one(1); + const T zero(0); + const auto transform = [&](const std::vector& input, std::vector& output) { + output.assign(ld_psi_ * retained_count, T(0)); + ModuleBase::gemm_op()('N', 'N', n_dim_, retained_count, selected_count, + &one, input.data(), ld_psi_, overlap.data(), selected_count, + &zero, output.data(), ld_psi_); + }; + transform(selected.vectors, orthogonal.vectors); + transform(selected.hvectors, orthogonal.hvectors); + transform(selected.svectors, orthogonal.svectors); + selected = OperatorBlock{}; + SmallSubspace expanded; + build_startup_global_subspace(psi, orthogonal.vectors, orthogonal.svectors, + orthogonal.hvectors, retained_count, expanded); + const int expanded_dim = n_band_ + retained_count; + local_success = false; + try + { + // No rank truncation of [X,Q]: X is retained intact and Q has already + // been rank-revealed independently in its S-orthogonal complement. + HermitianLapack::sygvd(expanded_dim, expanded.k.data(), expanded.m.data(), expanded.eval.data()); + local_success = true; + } + catch (const std::runtime_error&) + { + } + if (!all_pool_operations_succeeded(local_success)) + { + throw std::runtime_error("PPCG: competing-direction Rayleigh-Ritz solve failed."); + } + update_startup_global_subspace(psi, orthogonal.vectors, orthogonal.svectors, + orthogonal.hvectors, retained_count, expanded); + std::copy(expanded.eval.begin(), expanded.eval.begin() + n_band_, eigenvalue); + set_zero(p_); + set_zero(sp_); + set_zero(hp_); + compute_eigenpair_residual(eigenvalue); + lock_epairs(w_, ethr_band, active_cols); + ++rayleigh_ritz_count_; + ++root_recovery_count_; + return true; +} + +// --------------------------------------------------------------------------- +// Update wavefunctions from small subspace eigenvectors +// --------------------------------------------------------------------------- +template +void DiagoPPCG::update_one_block( + T* psi, + const std::vector& cols, + int l, + int nblk, + SmallSubspace& subspace) +{ + const PPCGFunctionTimer timer("update_one_block"); + const int dim = nblk * l; + const T* eigvec = subspace.k.data(); + + subspace.psi_new.assign(ld_psi_ * l, T(0)); + subspace.spsi_new.assign(ld_psi_ * l, T(0)); + subspace.hpsi_new.assign(ld_psi_ * l, T(0)); + subspace.p_new.assign(ld_psi_ * l, T(0)); + subspace.sp_new.assign(ld_psi_ * l, T(0)); + subspace.hp_new.assign(ld_psi_ * l, T(0)); + + // coeff_state: full Ritz-vector rows [psi, w, p] -> new iterate. + // coeff_p: the [w, p] rows (psi rows zero) -> new search direction. + subspace.coeff_state.resize(dim * l); + subspace.coeff_p.resize(dim * l); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) if (l * l > ppcg_openmp_work_threshold) +#endif + for (int j = 0; j < l; ++j) + { + for (int i = 0; i < l; ++i) + { + const T c_psi = eigvec[i + j * dim]; + const T c_w = eigvec[(l + i) + j * dim]; + subspace.coeff_state[i + j * dim] = c_psi; + subspace.coeff_state[(l + i) + j * dim] = c_w; + subspace.coeff_p[i + j * dim] = T(0); + subspace.coeff_p[(l + i) + j * dim] = c_w; + if (nblk >= 3) + { + const T c_p = eigvec[(2 * l + i) + j * dim]; + subspace.coeff_state[(2 * l + i) + j * dim] = c_p; + subspace.coeff_p[(2 * l + i) + j * dim] = c_p; + } + } + } + + auto fill_basis = [&](const std::vector& a, + const std::vector& b, + const std::vector& c, + std::vector& basis) + { + basis.resize(ld_psi_ * dim); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) if (ld_psi_ * l > ppcg_openmp_work_threshold) +#endif + for (int j = 0; j < l; ++j) + { + std::copy(a.begin() + j * ld_psi_, + a.begin() + (j + 1) * ld_psi_, + basis.begin() + j * ld_psi_); + std::copy(b.begin() + j * ld_psi_, + b.begin() + (j + 1) * ld_psi_, + basis.begin() + (l + j) * ld_psi_); + if (nblk >= 3) + { + std::copy(c.begin() + j * ld_psi_, + c.begin() + (j + 1) * ld_psi_, + basis.begin() + (2 * l + j) * ld_psi_); + } + } + }; + + auto combine = [&](const std::vector& basis, + const std::vector& coeff, + std::vector& out) + { + const T one = T(1); + const T zero = T(0); + ModuleBase::gemm_op()('N', + 'N', + n_dim_, + l, + dim, + &one, + basis.data(), + ld_psi_, + coeff.data(), + dim, + &zero, + out.data(), + ld_psi_); + }; + + fill_basis(subspace.psi_l, subspace.w_l, subspace.p_l, subspace.basis); + fill_basis(subspace.spsi_l, subspace.sw_l, subspace.sp_l, subspace.sbasis); + fill_basis(subspace.hpsi_l, subspace.hw_l, subspace.hp_l, subspace.hbasis); + + combine(subspace.basis, subspace.coeff_state, subspace.psi_new); + combine(subspace.sbasis, subspace.coeff_state, subspace.spsi_new); + combine(subspace.hbasis, subspace.coeff_state, subspace.hpsi_new); + combine(subspace.basis, subspace.coeff_p, subspace.p_new); + combine(subspace.sbasis, subspace.coeff_p, subspace.sp_new); + combine(subspace.hbasis, subspace.coeff_p, subspace.hp_new); + + scatter_cols(psi, cols, subspace.psi_new); + scatter_cols(spsi_.data(), cols, subspace.spsi_new); + scatter_cols(hpsi_.data(), cols, subspace.hpsi_new); + scatter_cols(p_.data(), cols, subspace.p_new); + scatter_cols(sp_.data(), cols, subspace.sp_new); + scatter_cols(hp_.data(), cols, subspace.hp_new); +} + +// --------------------------------------------------------------------------- +// Rayleigh-Ritz: full subspace diagonalization + residual computation +// --------------------------------------------------------------------------- +template +bool DiagoPPCG::rayleigh_ritz( + T* psi, Real* eigenvalue, + std::vector& active_cols, + const std::vector& ethr_band) +{ + const PPCGFunctionTimer timer("rayleigh_ritz"); + ++rayleigh_ritz_count_; + gram(psi, hpsi_.data(), n_band_, n_band_, rr_hsub_, n_band_); + gram(psi, spsi_.data(), n_band_, n_band_, rr_ssub_, n_band_); + + bool local_sygvd_ok = false; + try + { + HermitianLapack::sygvd(n_band_, rr_hsub_.data(), rr_ssub_.data(), + rr_eval_.data()); + local_sygvd_ok = true; + } + catch (const std::runtime_error&) + { + // The collective result below selects one common fallback on all ranks. + } + const bool sygvd_ok = all_pool_operations_succeeded(local_sygvd_ok); + if (!sygvd_ok) + { + // LAPACK may overwrite the projected matrices before failing. Re-form + // them on every rank so all ranks enter the same fallback path. + gram(psi, hpsi_.data(), n_band_, n_band_, rr_hsub_, n_band_); + gram(psi, spsi_.data(), n_band_, n_band_, rr_ssub_, n_band_); + for (int ii = 0; ii < n_band_; ++ii) + { + rr_eval_[ii] = Real(std::real(rr_hsub_[ii + ii * n_band_])) + / std::max(Real( + std::real(rr_ssub_[ii + ii * n_band_])), + Real(ppcg_numerical_threshold)); + } + } + + if (sygvd_ok) + { + const int sz = ld_psi_ * n_band_; + std::copy(psi, psi + sz, w_.begin()); + std::copy(spsi_.begin(), spsi_.end(), sw_.begin()); + std::copy(hpsi_.begin(), hpsi_.end(), hw_.begin()); + + std::fill(psi, psi + ld_psi_ * n_band_, T(0)); + set_zero(spsi_); + set_zero(hpsi_); + + const T one = T(1); + const T zero = T(0); + ModuleBase::gemm_op()('N', + 'N', + n_dim_, + n_band_, + n_band_, + &one, + w_.data(), + ld_psi_, + rr_hsub_.data(), + n_band_, + &zero, + psi, + ld_psi_); + ModuleBase::gemm_op()('N', + 'N', + n_dim_, + n_band_, + n_band_, + &one, + sw_.data(), + ld_psi_, + rr_hsub_.data(), + n_band_, + &zero, + spsi_.data(), + ld_psi_); + ModuleBase::gemm_op()('N', + 'N', + n_dim_, + n_band_, + n_band_, + &one, + hw_.data(), + ld_psi_, + rr_hsub_.data(), + n_band_, + &zero, + hpsi_.data(), + ld_psi_); + + for (int j = 0; j < n_band_; ++j) + { + eigenvalue[j] = rr_eval_[j]; + } + // Keep each conjugate direction attached to the same Ritz column + // when the global solve mixes or permutes the retained states. + const auto rotate_direction = [&](std::vector& values) { + std::copy(values.begin(), values.end(), w_.begin()); + set_zero(values); + ModuleBase::gemm_op()('N', 'N', n_dim_, n_band_, n_band_, + &one, w_.data(), ld_psi_, rr_hsub_.data(), n_band_, + &zero, values.data(), ld_psi_); + }; + rotate_direction(p_); + rotate_direction(sp_); + rotate_direction(hp_); + } + else + { + // No rotation: just update eigenvalues with Rayleigh quotients. + for (int j = 0; j < n_band_; ++j) + { + eigenvalue[j] = rr_eval_[j]; + } + } + + compute_eigenpair_residual(eigenvalue); + + if (sygvd_ok) + { + lock_epairs(w_, ethr_band, active_cols); + } + else + { + // Without a successful global Rayleigh-Ritz solve, the columns are + // not known to be ordered by eigenvalue. Do not hard-lock an + // arbitrary column based on its residual alone. + active_cols.resize(n_band_); + std::iota(active_cols.begin(), active_cols.end(), 0); + max_residual_ = max_generalized_residual(hpsi_.data(), + spsi_.data(), + eigenvalue, + ld_psi_, + n_dim_, + n_band_); + } + return sygvd_ok; +} + +template +void DiagoPPCG::compute_eigenpair_residual(const Real* eigenvalue) +{ + const PPCGFunctionTimer timer("compute_eigenpair_residual"); + set_zero(w_); +#ifdef _OPENMP +#pragma omp parallel for collapse(2) schedule(static) if (n_dim_ * n_band_ > ppcg_openmp_work_threshold) +#endif + for (int j = 0; j < n_band_; ++j) + { + for (int ig = 0; ig < n_dim_; ++ig) + { + w_[idx(ig, j, ld_psi_)] = hpsi_[idx(ig, j, ld_psi_)] + - spsi_[idx(ig, j, ld_psi_)] * eigenvalue[j]; + } + } +} + +template +bool DiagoPPCG::orthonormalize(T* psi) +{ + const PPCGFunctionTimer timer("orthonormalize"); + gram(psi, spsi_.data(), n_band_, n_band_, rr_ssub_, n_band_); + bool local_success = false; + try + { + container::kernels::lapack_potrf()( + 'U', n_band_, rr_ssub_.data(), n_band_); + container::kernels::lapack_trtri()( + 'U', 'N', n_band_, rr_ssub_.data(), n_band_); + local_success = true; + } + catch (const std::runtime_error&) + { + // The caller performs a common Rayleigh-Ritz fallback on all ranks. + } + if (!all_pool_operations_succeeded(local_success)) + { + return false; + } + for (int j = 0; j < n_band_; ++j) + { + for (int i = j + 1; i < n_band_; ++i) + { + rr_ssub_[i + j * n_band_] = T(0); + } + } + // X^H SX = R^H R. Right multiplication by the upper triangular R^{-1} + // keeps a previously locked prefix independent of the active columns. + const int size = ld_psi_ * n_band_; + const T one = T(1); + const T zero = T(0); + const auto rotate = [&](T* values) { + std::copy(values, values + size, w_.begin()); + std::fill(values, values + size, T(0)); + ModuleBase::gemm_op()('N', + 'N', + n_dim_, + n_band_, + n_band_, + &one, + w_.data(), + ld_psi_, + rr_ssub_.data(), + n_band_, + &zero, + values, + ld_psi_); + }; + rotate(psi); + rotate(hpsi_.data()); + rotate(spsi_.data()); + // The non-Ritz sweep changes the column basis too. Transform the + // conjugate directions and both cached images with the same coefficients. + rotate(p_.data()); + rotate(sp_.data()); + rotate(hp_.data()); + return true; +} + +template +bool DiagoPPCG::update_block_gradient( + const T* psi, + Real* eigenvalue, + const std::vector& ethr_band) +{ + const PPCGFunctionTimer timer("update_block_gradient"); + gram(psi, hpsi_.data(), n_band_, n_band_, rr_hsub_, n_band_); + hermitize_projected(rr_hsub_, n_band_); + for (int j = 0; j < n_band_; ++j) + { + // These are Rayleigh quotients for tracing, not ordered Ritz values. + eigenvalue[j] = std::real(rr_hsub_[j + j * n_band_]); + } + std::copy(hpsi_.begin(), hpsi_.end(), w_.begin()); + const T minus_one = T(-1); + const T one = T(1); + ModuleBase::gemm_op()('N', + 'N', + n_dim_, + n_band_, + n_band_, + &minus_one, + spsi_.data(), + ld_psi_, + rr_hsub_.data(), + n_band_, + &one, + w_.data(), + ld_psi_); + std::vector candidates; + lock_epairs(w_, ethr_band, candidates); + // Do not change the actual locked prefix without a successful Ritz solve. + return candidates.empty(); +} + +} // namespace hsolver + + +namespace hsolver { + +//============================================================================== +// MAIN DIAGONALIZATION ROUTINE +//============================================================================== +template +double DiagoPPCG::diag(const HPsiFunc& hpsi_func, + const SPsiFunc& spsi_func, + int ld_psi, + int nband, + int dim, + T* psi_in, + Real* eigenvalue_in, + const std::vector& ethr_band, + const Real* prec) +{ + return diag(hpsi_func, spsi_func, ld_psi, nband, dim, + psi_in, eigenvalue_in, ethr_band, prec, std::clog); +} + +template +double DiagoPPCG::diag(const HPsiFunc& hpsi_func, + const SPsiFunc& spsi_func, + int ld_psi, + int nband, + int dim, + T* psi_in, + Real* eigenvalue_in, + const std::vector& ethr_band, + const Real* prec, + std::ostream& log) +{ + const PPCGFunctionTimer timer("diag"); + ld_psi_ = ld_psi; + n_band_ = nband; + n_dim_ = dim; + converged_ = false; + active_band_count_ = 0; + iteration_limit_reached_ = false; + max_residual_ = Real(0); + first_active_band_ = -1; + first_active_threshold_ = Real(0); + rayleigh_ritz_count_ = 0; + root_recovery_count_ = 0; + + validate_input(hpsi_func, psi_in, eigenvalue_in, ethr_band, prec); + spsi_func_ = spsi_func; + + // Allocate working storage. + const int ncol = n_band_; + const int sz = ld_psi_ * ncol; + + hpsi_.assign(sz, T(0)); + spsi_.assign(sz, T(0)); + w_.assign(sz, T(0)); + sw_.assign(sz, T(0)); + hw_.assign(sz, T(0)); + p_.assign(sz, T(0)); + sp_.assign(sz, T(0)); + hp_.assign(sz, T(0)); + rr_hsub_.resize(ncol * ncol); + rr_ssub_.resize(ncol * ncol); + rr_eval_.resize(ncol); + + std::vector all_cols(ncol); + std::iota(all_cols.begin(), all_cols.end(), 0); + + apply_h(hpsi_func, psi_in, hpsi_.data(), ncol); + apply_s_current(psi_in, spsi_.data(), ncol); + + double avg_iter = 1.0; + int iter = 1; + std::vector active_cols; + active_cols.reserve(ncol); + + std::ofstream residual_trace; + const char* residual_trace_path = std::getenv("ABACUS_PPCG_RESIDUAL_TRACE"); + if (residual_trace_path != nullptr && Parallel_Reduce::is_world_root()) + { + // A single writer avoids file corruption when several MPI pools solve + // different k-points concurrently. + residual_trace.open(residual_trace_path); + if (residual_trace) + { +#ifdef _WIN32 + char* full_path = ::_fullpath(nullptr, residual_trace_path, 0); +#else + char* full_path = ::realpath(residual_trace_path, nullptr); +#endif + if (full_path != nullptr) + { + const std::string trace_filename(full_path); + std::free(full_path); + log << "PPCG residual trace written: " << trace_filename << std::endl; + } + else + { + log << "WARNING: PPCG residual trace opened, but the full path could not be resolved: " + << residual_trace_path << std::endl; + } + residual_trace << std::setprecision(std::numeric_limits::max_digits10); + residual_trace << "iteration,stage,max_residual"; + // Test-side tools compare these Ritz values with analytical or + // LAPACK reference eigenvalues without coupling them to the solver. + for (int ib = 0; ib < ncol; ++ib) + { + residual_trace << ",eigenvalue_" << ib; + } + residual_trace << ",local_overlap_min,local_overlap_max,local_overlap_condition\n"; + } + } + const bool write_residual_trace = residual_trace.is_open() && residual_trace.good(); + int residual_trace_writer_count = write_residual_trace ? 1 : 0; + reduce_pool_values(&residual_trace_writer_count, 1); + const bool residual_trace_enabled = residual_trace_writer_count > 0; + Real local_overlap_min = std::numeric_limits::quiet_NaN(); + Real local_overlap_max = std::numeric_limits::quiet_NaN(); + Real local_overlap_condition = std::numeric_limits::quiet_NaN(); + auto record_residual = [&](int iteration, const char* stage) { + if (!residual_trace_enabled) + { + return; + } + const Real max_residual = max_generalized_residual(hpsi_.data(), + spsi_.data(), + eigenvalue_in, + ld_psi_, + n_dim_, + ncol); + if (!write_residual_trace) + { + return; + } + residual_trace << iteration << ',' << stage << ',' << max_residual; + for (int ib = 0; ib < ncol; ++ib) + { + residual_trace << ',' << eigenvalue_in[ib]; + } + residual_trace << ',' << local_overlap_min + << ',' << local_overlap_max + << ',' << local_overlap_condition + << '\n'; + }; + + auto measure_local_overlap = [&](const SmallSubspace& local_subspace, + const int local_dim) { + if (!write_residual_trace) + { + return; + } + std::vector overlap = local_subspace.m; + std::vector overlap_eigenvalues(local_dim); + try + { + container::kernels::lapack_heevd()( + local_dim, + overlap.data(), + local_dim, + overlap_eigenvalues.data()); + } + catch (const std::exception&) + { + // Trace diagnostics must not change the solver's MPI control flow. + local_overlap_min = std::numeric_limits::quiet_NaN(); + local_overlap_max = std::numeric_limits::quiet_NaN(); + local_overlap_condition = std::numeric_limits::infinity(); + return; + } + const Real min_eigenvalue = overlap_eigenvalues.front(); + const Real max_eigenvalue = overlap_eigenvalues.back(); + const Real denominator = std::max(std::abs(min_eigenvalue), + Real(ppcg_numerical_threshold)); + const Real condition = std::abs(max_eigenvalue) / denominator; + if (std::isnan(local_overlap_condition) || condition > local_overlap_condition) + { + local_overlap_min = min_eigenvalue; + local_overlap_max = max_eigenvalue; + local_overlap_condition = condition; + } + }; + + // Initialize with Rayleigh-Ritz. + bool last_rr_succeeded = rayleigh_ritz( + psi_in, eigenvalue_in, active_cols, ethr_band); + const bool needs_startup_global_correction = !active_cols.empty() + && sbsize_ < ncol; + // The Ritz transformation updates X, HX and SX with identical coefficients. + record_residual(0, "initial_rr"); + + std::vector active_vector; + std::vector s_active_vector; + std::vector h_active_vector; + std::vector cols; + cols.reserve(std::min(sbsize_, ncol)); + SmallSubspace subspace; + bool use_p = false; // previous search direction becomes available after + // the first block update. + Real best_res = std::numeric_limits::max(); // best residual seen + int no_improve = 0; // iterations since the residual last decreased + + while (!active_cols.empty() && iter <= maxiter_) + { + OperatorBlock discarded; + // Near the floating-point cancellation scale, form a reliable + // residual before using it as a new search direction. This safeguard + // also covers cancellation in the W/P recurrence, not only in the + // explicit projection of P. The final fresh residual remains the + // convergence certificate; this scale is not an error bound. + Real spectral_scale = Real(1); + for (int ib = 0; ib < ncol; ++ib) + { + spectral_scale = std::max(spectral_scale, std::abs(eigenvalue_in[ib])); + } + const Real reliable_update_threshold = std::sqrt(std::numeric_limits::epsilon()) + * spectral_scale; + const bool near_precision = max_residual_ <= reliable_update_threshold; + if (near_precision) + { + apply_h(hpsi_func, psi_in, hpsi_.data(), ncol); + apply_s_current(psi_in, spsi_.data(), ncol); + apply_h(hpsi_func, p_.data(), hp_.data(), ncol); + apply_s_current(p_.data(), sp_.data(), ncol); + if (last_rr_succeeded) + { + compute_eigenpair_residual(eigenvalue_in); + } + else + { + update_block_gradient(psi_in, eigenvalue_in, ethr_band); + } + } + local_overlap_min = std::numeric_limits::quiet_NaN(); + local_overlap_max = std::numeric_limits::quiet_NaN(); + local_overlap_condition = std::numeric_limits::quiet_NaN(); + const int nact = int(active_cols.size()); + const int nsb = 1 + (nact - 1) / sbsize_; + // Keep the same number of local blocks and the same size limit, but + // distribute the remainder evenly. When the limit permits it, this + // avoids a trailing one-column block that cannot exchange search + // directions with a neighboring degenerate Ritz vector. + const int base_block_size = nact / nsb; + const int larger_block_count = nact % nsb; + + // Precondition the residual. + divide_by_preconditioner(active_cols, prec, w_); + copy_cols(w_.data(), active_cols, active_vector); + s_active_vector.assign(ld_psi_ * nact, T(0)); + apply_s_current(active_vector.data(), s_active_vector.data(), nact); + scatter_cols(sw_.data(), active_cols, s_active_vector); + project_against(psi_in, spsi_.data(), nullptr, all_cols, w_, sw_, nullptr, active_cols); + + // Apply H to the search direction. + copy_cols(w_.data(), active_cols, active_vector); + h_active_vector.assign(ld_psi_ * nact, T(0)); + copy_cols(sw_.data(), active_cols, s_active_vector); + apply_h(hpsi_func, active_vector.data(), h_active_vector.data(), nact); + scatter_cols(hw_.data(), active_cols, h_active_vector); + + // Project the cached P, SP and HP together. Linearity preserves their + // consistency without another Hamiltonian/overlap application except + // when strong cancellation requires a fresh image of selected columns. + if (use_p) + { + project_search_directions(hpsi_func, psi_in, all_cols, active_cols); + } + + avg_iter += double(nact) / double(ncol); + + // LOBPCG-style block subspace. On the first sweep only [psi, w] is + // available; afterwards the previous search direction p is added as a + // third block, which restores the conjugate-gradient acceleration. + // The w/p blocks are normalized to unit S-norm before building the + // Gram matrix (see build_small_subspace), keeping M well-conditioned. + + // The startup correction needs a global [X, W_active] comparison to + // retain the requested lowest roots. All current Ritz vectors remain + // available, while a converged prefix contributes no correction + // vectors and stays outside the active iteration count. + const bool startup_correction = iter == 1 && needs_startup_global_correction; + if (startup_correction) + { + const int local_dim = ncol + nact; + build_startup_global_subspace(psi_in, + active_vector, + s_active_vector, + h_active_vector, + nact, + subspace); + measure_local_overlap(subspace, local_dim); + solve_small_generalized(local_dim, ncol, subspace); + update_startup_global_subspace(psi_in, + active_vector, + s_active_vector, + h_active_vector, + nact, + subspace); + } + else + { + const int nblk = use_p ? 3 : 2; + int i0 = 0; + for (int isb = 0; isb < nsb; ++isb) + { + int l = base_block_size; + if (isb < larger_block_count) + { + ++l; + } + cols.assign(active_cols.begin() + i0, + active_cols.begin() + i0 + l); + + build_small_subspace(psi_in, cols, nblk, subspace); + const int local_dim = nblk * l; + measure_local_overlap(subspace, local_dim); + const int valid_ritz_count = solve_small_generalized(local_dim, l, subspace); + collect_discarded_ritz(subspace, l, valid_ritz_count, discarded); + update_one_block(psi_in, cols, l, nblk, subspace); + i0 += l; + } + } + use_p = true; + + // Refresh cached images independently of the full Ritz cadence, so + // rr_step=1 does not add two full H/S applications on every sweep. + // Intermediate sweeps maintain S-orthonormality without diagonalizing + // the full projected Hamiltonian. + const bool scheduled_rr = (iter % rr_step_) == 0; + const bool refresh_operators = (iter % ppcg_operator_refresh_interval) == 0 + && !near_precision; + if (refresh_operators) + { + apply_h(hpsi_func, psi_in, hpsi_.data(), ncol); + apply_s_current(psi_in, spsi_.data(), ncol); + apply_h(hpsi_func, p_.data(), hp_.data(), ncol); + apply_s_current(p_.data(), sp_.data(), ncol); + } + // The startup subspace can insert a newly discovered lower root into + // an old locked prefix. Reclassify it before freezing any column. + bool perform_rr = scheduled_rr || startup_correction || iter == maxiter_; + if (!perform_rr) + { + perform_rr = !orthonormalize(psi_in); + } + if (!perform_rr) + { + perform_rr = update_block_gradient(psi_in, eigenvalue_in, ethr_band); + } + last_rr_succeeded = false; + if (perform_rr) + { + last_rr_succeeded = rayleigh_ritz( + psi_in, eigenvalue_in, active_cols, ethr_band); + } + // A residual only certifies the retained eigenpairs, not that no lower + // state was discarded by another local block. Inspect those available + // directions before accepting a newly converged set of columns. + if (!perform_rr || last_rr_succeeded) + { + const bool recovered = recover_competing_ritz(hpsi_func, psi_in, eigenvalue_in, + last_rr_succeeded, active_cols, + ethr_band, discarded); + if (recovered) + { + last_rr_succeeded = true; + perform_rr = true; + use_p = false; + } + } + if (last_rr_succeeded && active_cols.empty()) + { + // A small cached residual is only tentative. Refresh the + // operator images before stopping, and keep iterating if + // roundoff had made the cached criterion overly optimistic. + apply_h(hpsi_func, psi_in, hpsi_.data(), ncol); + apply_s_current(psi_in, spsi_.data(), ncol); + compute_eigenpair_residual(eigenvalue_in); + lock_epairs(w_, ethr_band, active_cols); + if (!active_cols.empty()) + { + set_zero(p_); + set_zero(sp_); + set_zero(hp_); + use_p = false; + } + } + // Restart the search direction if the residual keeps rising. With a + // poor preconditioner the LOBPCG recurrence can stagnate (or slowly + // diverge) instead of reducing the residual. Requiring several + // consecutive rises avoids resetting on a transient bump, and a reset + // falls back to a steepest-descent step to recover the low eigenpairs. + { + const Real cur_res = max_residual_; + if (cur_res < best_res) + { + best_res = cur_res; + no_improve = 0; + } + else + { + ++no_improve; + } + if (no_improve >= ppcg_stagnation_restart_interval) + { + std::fill(p_.begin(), p_.end(), T(0)); + std::fill(sp_.begin(), sp_.end(), T(0)); + std::fill(hp_.begin(), hp_.end(), T(0)); + use_p = false; + no_improve = 0; + best_res = cur_res; + } + } + const char* stage = perform_rr ? "rayleigh_ritz" : "block_sweep"; + record_residual(iter, stage); + + ++iter; + } + + // Final consistency: ensure H|psi> and S|psi> match the returned vectors, + // then classify convergence from those final residuals rather than from a + // cached Rayleigh-Ritz state. + apply_h(hpsi_func, psi_in, hpsi_.data(), ncol); + apply_s_current(psi_in, spsi_.data(), ncol); + compute_eigenpair_residual(eigenvalue_in); + if (last_rr_succeeded) + { + lock_epairs(w_, ethr_band, active_cols); + } + else + { + active_cols = all_cols; + max_residual_ = max_generalized_residual(hpsi_.data(), + spsi_.data(), + eigenvalue_in, + ld_psi_, + n_dim_, + ncol); + } + active_band_count_ = int(active_cols.size()); + if (!active_cols.empty()) + { + first_active_band_ = active_cols.front(); + first_active_threshold_ = std::max(Real(ethr_band[first_active_band_]), + diag_thr_); + } + converged_ = last_rr_succeeded && active_cols.empty(); + iteration_limit_reached_ = !converged_ && iter > maxiter_; + if (!active_cols.empty()) + { + const char* failure_stage = iteration_limit_reached_ + ? "max_iterations" + : "final_residual_failure"; + record_residual(iter - 1, failure_stage); + } + record_residual(iter - 1, "final"); + return avg_iter; +} + +} // namespace hsolver + +namespace hsolver { + +template class DiagoPPCG, base_device::DEVICE_CPU>; +template class DiagoPPCG, base_device::DEVICE_CPU>; + +} // namespace hsolver diff --git a/source/source_hsolver/diago_ppcg.h b/source/source_hsolver/diago_ppcg.h new file mode 100644 index 00000000000..399c7d5057b --- /dev/null +++ b/source/source_hsolver/diago_ppcg.h @@ -0,0 +1,333 @@ +#ifndef DIAGO_PPCG_H +#define DIAGO_PPCG_H + +#include "source_base/module_device/types.h" + +#include +#include +#include +#include +#include + +namespace hsolver { + +// ----------------------------------------------------------------------------- +// DiagoPPCG: Projection Preconditioned solver +// ----------------------------------------------------------------------------- +// +// Implements the block-subspace diagonalization strategy used by +// ks_solver=ppcg. The solver updates groups of active Ritz pairs in a bounded +// local subspace with configurable global Rayleigh-Ritz steps. Operator images +// are cached, with independent periodic and cancellation-triggered refreshes. +// ----------------------------------------------------------------------------- + +namespace base_device = ::base_device; + +// Check all possible CPU workspaces before allocating or copying through the +// GPU bridge. Dense kernels and pool reductions use int element counts. +bool ppcg_workspace_dimensions_valid(int ld_psi, int nband, int dim, int sbsize); + +template +class DiagoPPCG +{ +public: + // ------------------------------------------------------------------------- + // Type aliases + // ------------------------------------------------------------------------- + using Real = typename std::conditional< + std::is_same>::value, double, + float>::type; + using HPsiFunc = std::function; + using SPsiFunc = std::function; + + // ------------------------------------------------------------------------- + // Constructor + // rr_step is the explicit Ritz interval; PW INPUT defaults are supplied by + // the caller. The internal H/S cache-refresh interval is independent of it. + // ------------------------------------------------------------------------- + DiagoPPCG(const Real& diag_thr, + const int& diag_iter_max, + const int& sbsize, + const int& rr_step); + + // ------------------------------------------------------------------------- + // Main entry point + // + // Returns average number of subspace iterations per band. + // Optional ABACUS_PPCG_RESIDUAL_TRACE CSV: initial_rr, rayleigh_ritz and + // final rows contain ordered Ritz values after a successful Ritz solve. + // block_sweep rows instead contain the current column Rayleigh quotients; + // exclude those rows when plotting ordered eigenvalue errors. + // ------------------------------------------------------------------------- + double diag(const HPsiFunc& hpsi_func, + const SPsiFunc& spsi_func, + int ld_psi, + int nband, + int dim, + T* psi_in, + Real* eigenvalue_in, + const std::vector& ethr_band, + const Real* prec); + + // The PW caller supplies its running log; standalone callers of the + // overload above use std::clog for optional trace-file diagnostics. + double diag(const HPsiFunc& hpsi_func, + const SPsiFunc& spsi_func, + int ld_psi, + int nband, + int dim, + T* psi_in, + Real* eigenvalue_in, + const std::vector& ethr_band, + const Real* prec, + std::ostream& log); + + bool converged() const + { + return converged_; + } + + int active_band_count() const + { + return active_band_count_; + } + + bool iteration_limit_reached() const + { + return iteration_limit_reached_; + } + + Real max_residual() const + { + return max_residual_; + } + + int first_active_band() const + { + return first_active_band_; + } + + Real first_active_threshold() const + { + return first_active_threshold_; + } + + int rayleigh_ritz_count() const + { + return rayleigh_ritz_count_; + } + + int root_recovery_count() const + { + return root_recovery_count_; + } + +private: + friend class DiagoPPCGProjectionTest; + // ------------------------------------------------------------------------- + // Data members + // ------------------------------------------------------------------------- + int maxiter_; + int sbsize_; + int rr_step_; + Real diag_thr_; + + // Problem dimensions (set in diag()) + int ld_psi_ = 0; + int n_band_ = 0; + int n_dim_ = 0; + bool converged_ = false; + int active_band_count_ = 0; + bool iteration_limit_reached_ = false; + Real max_residual_ = Real(0); + int first_active_band_ = -1; + Real first_active_threshold_ = Real(0); + int rayleigh_ritz_count_ = 0; + int root_recovery_count_ = 0; + + // Cached S-operator (null if identity). + SPsiFunc spsi_func_; + + // Working storage (column-major: ld_psi_ rows, n_band_ columns). + std::vector hpsi_; + std::vector spsi_; + std::vector w_; // residual / preconditioned residual + std::vector sw_; // S * w + std::vector hw_; // H * w + std::vector p_; // previous search direction (LOBPCG) + std::vector sp_; // S * p + std::vector hp_; // H * p + std::vector rr_hsub_; + std::vector rr_ssub_; + std::vector rr_eval_; + + // ------------------------------------------------------------------------- + // Internal helpers + // ------------------------------------------------------------------------- + static inline int idx(int row, int col, int ld) + { + return row + col * ld; + } + + void validate_input(const HPsiFunc& hpsi_func, + const T* psi_in, const Real* eigenvalue_in, + const std::vector& ethr_band, + const Real* prec) const; + + // S-application (identity fallback if spsi_func is null). + void apply_h(const HPsiFunc& hpsi_func, T* psi_in, T* hpsi_out, + int ncol) const; + void apply_s(const SPsiFunc& spsi_func, T* psi_in, T* spsi_out, + int ncol) const; + void apply_s_current(T* psi_in, T* spsi_out, int ncol) const; + + // Inner product (real part only). + Real gamma_dot(const T* x, const T* y) const; + + // Gram matrix: out[i, j] = . + void gram(const T* mat_a, const T* mat_b, + int ncol_a, int ncol_b, + std::vector& out, int ld_out) const; + + // Gather / scatter columns. + void copy_cols(const T* src, const std::vector& cols, + std::vector& dst) const; + void scatter_cols(T* dst, const std::vector& cols, + const std::vector& src) const; + + // Project x onto vectors orthogonal to the S-orthonormal basis. + void project_against(const T* basis, const T* sbasis, const T* hbasis, + const std::vector& basis_cols, + std::vector& x, std::vector& sx, std::vector* hx, + const std::vector& x_cols) const; + + void project_search_directions(const HPsiFunc& hpsi_func, const T* psi, + const std::vector& all_cols, + const std::vector& active_cols); + + // x[c] /= max(prec, eps) for each active column c. + void divide_by_preconditioner(const std::vector& active_cols, + const Real* prec, + std::vector& x) const; + + // ------------------------------------------------------------------------- + // Block-subspace strategy helpers (File 1 style) + // ------------------------------------------------------------------------- + struct SmallSubspace + { + std::vector k; // K matrix (projected H) + std::vector m; // M matrix (projected S) + std::vector eval; // eigenvalues + std::vector psi_l; + std::vector spsi_l; + std::vector hpsi_l; + std::vector w_l; + std::vector sw_l; + std::vector hw_l; + std::vector p_l; + std::vector sp_l; + std::vector hp_l; + std::vector basis; + std::vector hbasis; + std::vector sbasis; + std::vector coeff_state; + std::vector coeff_p; + std::vector psi_new; + std::vector spsi_new; + std::vector hpsi_new; + std::vector p_new; + std::vector sp_new; + std::vector hp_new; + }; + + struct OperatorBlock + { + std::vector vectors; + std::vector hvectors; + std::vector svectors; + }; + + void scale_to_unit_snorm(std::vector& x, + std::vector& sx, + std::vector& hx, + int ncols) const; + + void hermitize_projected(std::vector& matrix, int dim) const; + + void insert_startup_gram_block(const T* left, + const T* right, + int row_count, + int column_count, + int row_offset, + int column_offset, + int matrix_dim, + std::vector& matrix, + std::vector& workspace) const; + + void build_startup_global_subspace(const T* psi, + std::vector& w_active, + std::vector& sw_active, + std::vector& hw_active, + int active_count, + SmallSubspace& subspace); + + void update_startup_global_subspace(T* psi, + const std::vector& w_active, + const std::vector& sw_active, + const std::vector& hw_active, + int active_count, + const SmallSubspace& subspace); + + void lock_epairs(const std::vector& residual, + const std::vector& ethr_band, + std::vector& active_cols); + + void build_small_subspace(const T* psi, + const std::vector& cols, + int nblk, + SmallSubspace& subspace) const; + + // Returns the number of valid Ritz vectors, including discarded states. + // An identity fallback reports zero rather than pretending it is a solve. + int solve_small_generalized(const int dim, + const int nstates, + SmallSubspace& subspace) const; + + void collect_discarded_ritz(const SmallSubspace& subspace, + int kept_count, int valid_count, + OperatorBlock& candidates) const; + + void project_candidate_complement(const HPsiFunc& hpsi_func, + const T* psi, OperatorBlock& candidates) const; + + bool recover_competing_ritz(const HPsiFunc& hpsi_func, + T* psi, Real* eigenvalue, + bool ordered_ritz, + std::vector& active_cols, + const std::vector& ethr_band, + OperatorBlock& candidates); + + void update_one_block(T* psi, + const std::vector& cols, + int l, + int nblk, + SmallSubspace& subspace); + + bool rayleigh_ritz(T* psi, Real* eigenvalue, + std::vector& active_cols, + const std::vector& ethr_band); + + void compute_eigenpair_residual(const Real* eigenvalue); + + bool orthonormalize(T* psi); + + // Between Ritz solves the columns need not be individual eigenvectors. + // Use the block gradient, and request a Ritz solve before declaring them + // converged or changing the locked prefix. + bool update_block_gradient(const T* psi, Real* eigenvalue, + const std::vector& ethr_band); +}; + +} // namespace hsolver + +#endif // DIAGO_PPCG_H diff --git a/source/source_hsolver/hsolver_pw.cpp b/source/source_hsolver/hsolver_pw.cpp index 1b316b567e0..9fd620aba0b 100644 --- a/source/source_hsolver/hsolver_pw.cpp +++ b/source/source_hsolver/hsolver_pw.cpp @@ -1,6 +1,8 @@ #include "hsolver_pw.h" #include "source_base/parallel_comm.h" +#include "source_base/parallel_reduce.h" +#include "source_base/module_device/memory_op.h" #include "source_base/timer.h" #include "source_base/tool_quit.h" #include "source_estate/elecstate_pw.h" @@ -10,16 +12,175 @@ #include "source_hsolver/diago_cg.h" #include "source_hsolver/diago_dav_subspace.h" #include "source_hsolver/diago_david.h" +#include "source_hsolver/diago_ppcg.h" #include "source_hsolver/diago_iter_assist.h" #include "source_psi/psi.h" #include +#include +#include +#include #include +#include +#include #include namespace hsolver { +namespace +{ +struct PPCGRunResult +{ + double average_iterations; + bool converged; + int active_band_count; + bool iteration_limit_reached; + double max_residual; +}; + +template +PPCGRunResult run_ppcg_pw(const HPsiFunc& hpsi_func, + const SPsiFunc& spsi_func, + const int ld_psi, + const int nband, + const int dim, + T* psi, + Real* eigenvalue, + const std::vector& ethr_band, + const Real* pre_condition, + const double diag_thr, + const int diag_iter_max, + const int pw_diag_ndim, + const int rr_step, + std::ostream& log, + std::true_type) +{ + const int sbsize = std::max(1, std::min(nband, pw_diag_ndim)); + const int rr_step_safe = std::max(1, rr_step); + const Real solver_diag_thr = static_cast(diag_thr); + DiagoPPCG ppcg(solver_diag_thr, diag_iter_max, sbsize, rr_step_safe); + PPCGRunResult result; + result.average_iterations = ppcg.diag(hpsi_func, + spsi_func, + ld_psi, + nband, + dim, + psi, + eigenvalue, + ethr_band, + pre_condition, + log); + result.converged = ppcg.converged(); + result.active_band_count = ppcg.active_band_count(); + result.iteration_limit_reached = ppcg.iteration_limit_reached(); + result.max_residual = double(ppcg.max_residual()); + return result; +} + +template +PPCGRunResult run_ppcg_pw(const HPsiFunc& hpsi_func, + const SPsiFunc& spsi_func, + const int ld_psi, + const int nband, + const int dim, + T* psi, + Real* eigenvalue, + const std::vector& ethr_band, + const Real* pre_condition, + const double diag_thr, + const int diag_iter_max, + const int pw_diag_ndim, + const int rr_step, + std::ostream& log, + std::false_type) +{ + const int sbsize = std::max(1, std::min(nband, pw_diag_ndim)); + const int rr_step_safe = std::max(1, rr_step); + const bool valid_dimensions = ppcg_workspace_dimensions_valid(ld_psi, nband, dim, sbsize); + int invalid_dimensions = valid_dimensions ? 0 : 1; + Parallel_Reduce::reduce_pool(&invalid_dimensions, 1); + if (invalid_dimensions != 0) + { + ModuleBase::WARNING_QUIT("PPCG GPU bridge", "Invalid or oversized workspace dimensions on a pool rank."); + } + const int nelem = ld_psi * nband; + const Real solver_diag_thr = static_cast(diag_thr); + struct DeviceBuffer + { + T* ptr = nullptr; + explicit DeviceBuffer(const int size) + { + base_device::memory::resize_memory_op()(ptr, size, "PPCG device bridge"); + } + ~DeviceBuffer() + { + if (ptr != nullptr) + { + base_device::memory::delete_memory_op()(ptr); + } + } + DeviceBuffer(const DeviceBuffer&) = delete; + DeviceBuffer& operator=(const DeviceBuffer&) = delete; + }; + std::vector psi_host(nelem, T(0)); + T* psi_host_ptr = psi_host.data(); + base_device::memory::synchronize_memory_op()(psi_host_ptr, psi, nelem); + DeviceBuffer psi_dev(nelem); + DeviceBuffer out_dev(nelem); + const auto copy_to_device = [&](T* host_ptr, const int ld, const int nvec) { + int invalid_block = ld <= 0 || nvec <= 0 || nvec > nelem / ld ? 1 : 0; + Parallel_Reduce::reduce_pool(&invalid_block, 1); + if (invalid_block != 0) + { + ModuleBase::WARNING_QUIT("PPCG GPU bridge", "Input block exceeds allocated workspace on a pool rank."); + } + const int count = ld * nvec; + base_device::memory::synchronize_memory_op()( + psi_dev.ptr, host_ptr, count); + }; + const auto copy_from_device = [&](T* host_ptr, const int ld, const int nvec) { + int invalid_block = ld <= 0 || nvec <= 0 || nvec > nelem / ld ? 1 : 0; + Parallel_Reduce::reduce_pool(&invalid_block, 1); + if (invalid_block != 0) + { + ModuleBase::WARNING_QUIT("PPCG GPU bridge", "Output block exceeds allocated workspace on a pool rank."); + } + const int count = ld * nvec; + base_device::memory::synchronize_memory_op()( + host_ptr, out_dev.ptr, count); + }; + const auto bridge_hpsi = [&](T* psi_in, T* hpsi_out, const int ld, const int nvec) { + copy_to_device(psi_in, ld, nvec); + hpsi_func(psi_dev.ptr, out_dev.ptr, ld, nvec); + copy_from_device(hpsi_out, ld, nvec); + }; + const auto bridge_spsi = [&](T* psi_in, T* spsi_out, const int ld, const int nvec) { + copy_to_device(psi_in, ld, nvec); + spsi_func(psi_dev.ptr, out_dev.ptr, ld, nvec); + copy_from_device(spsi_out, ld, nvec); + }; + DiagoPPCG ppcg(solver_diag_thr, diag_iter_max, sbsize, rr_step_safe); + PPCGRunResult result; + result.average_iterations = ppcg.diag(bridge_hpsi, + bridge_spsi, + ld_psi, + nband, + dim, + psi_host_ptr, + eigenvalue, + ethr_band, + pre_condition, + log); + result.converged = ppcg.converged(); + result.active_band_count = ppcg.active_band_count(); + result.iteration_limit_reached = ppcg.iteration_limit_reached(); + result.max_residual = double(ppcg.max_residual()); + base_device::memory::synchronize_memory_op()(psi, psi_host_ptr, nelem); + return result; +} +} // namespace + template void HSolverPW::cal_smooth_ethr(const double& wk, const double* wg, @@ -79,7 +240,7 @@ void HSolverPW::solve(HSOperator& op, this->nproc_in_pool = nproc_in_pool_in; // report if the specified diagonalization method is not supported - const std::initializer_list _methods = {"cg", "dav", "dav_subspace", "bpcg"}; + const std::initializer_list _methods = {"cg", "dav", "dav_subspace", "bpcg", "ppcg"}; if (std::find(std::begin(_methods), std::end(_methods), this->method) == std::end(_methods)) { ModuleBase::WARNING_QUIT("HSolverPW::solve", "This type of eigensolver is not supported!"); @@ -131,7 +292,8 @@ void HSolverPW::solve(HSOperator& op, // solve eigenvector and eigenvalue for H(k) - this->hamiltSolvePsiK(op, psi, precondition, eigenvalues.data() + ik * psi.get_nbands(), this->wfc_basis->nks); + Real* k_eigenvalues = eigenvalues.data() + ik * psi.get_nbands(); + this->hamiltSolvePsiK(op, psi, precondition, k_eigenvalues, this->wfc_basis->nks, log); if (skip_charge) { @@ -169,7 +331,8 @@ void HSolverPW::solve(HSOperator& op, // solve eigenvector and eigenvalue for H(k) - this->hamiltSolvePsiK(op, psi, precondition, eigenvalues.data() + ik * psi.get_nbands(), this->wfc_basis->nks); + Real* k_eigenvalues = eigenvalues.data() + ik * psi.get_nbands(); + this->hamiltSolvePsiK(op, psi, precondition, k_eigenvalues, this->wfc_basis->nks, log); // output iteration information and reset avg_iter if (skip_charge) @@ -228,7 +391,8 @@ void HSolverPW::hamiltSolvePsiK(const HSOperator& op, psi::Psi& psi, std::vector& pre_condition, Real* eigenvalue, - const int& nk_nums) + const int& nk_nums, + std::ostream& log) { ModuleBase::timer::start("HSolverPW", "solve_psik"); #ifdef __MPI @@ -347,6 +511,67 @@ void HSolverPW::hamiltSolvePsiK(const HSOperator& op, ntry_max, notconv_max)); } + else if (this->method == "ppcg") + { + const auto hpsi_func = [&op](T* psi_in, T* hpsi_out, const int ld, const int nvec) { + op.hpsi(psi_in, hpsi_out, ld, nvec); + }; + const auto spsi_func = [&op](T* psi_in, T* spsi_out, const int ld, const int nvec) { + op.spsi(psi_in, spsi_out, ld, nvec); + }; + const std::is_same cpu_path; + const int ld_psi = psi.get_nbasis(); + const int nband = psi.get_nbands(); + const int dim = psi.get_current_ngk(); + T* psi_pointer = psi.get_pointer(); + const Real* precondition_pointer = pre_condition.data(); + try + { + const PPCGRunResult ppcg_result = run_ppcg_pw( + hpsi_func, + spsi_func, + ld_psi, + nband, + dim, + psi_pointer, + eigenvalue, + this->ethr_band, + precondition_pointer, + this->diag_thr, + this->diag_iter_max, + this->pw_diag_ndim, + 1, + log, + cpu_path); + if (!ppcg_result.converged) + { + std::ostringstream message; + message << "PPCG did not converge: active bands=" << ppcg_result.active_band_count + << ", max residual=" << ppcg_result.max_residual; + if (ppcg_result.iteration_limit_reached) + { + message << ", iteration limit reached (pw_diag_nmax=" << this->diag_iter_max << ")"; + } + const std::string diagnostic = message.str(); + if (this->rank_in_pool == 0) + { + log << "WARNING: " << diagnostic << std::endl; + // Other k-point pools may not own an open running log. + std::cerr << "WARNING: " << diagnostic << std::endl; + } + if (this->calculation_type == "nscf" || !std::isfinite(ppcg_result.max_residual)) + { + ModuleBase::WARNING_QUIT("HSolverPW::hamiltSolvePsiK", diagnostic); + } + } + DiagoIterAssist::avg_iter += ppcg_result.average_iterations; + } + catch (const std::exception& error) + { + const std::string message = std::string("PPCG diagonalization failed: ") + error.what(); + ModuleBase::WARNING_QUIT("HSolverPW::hamiltSolvePsiK", message); + } + } ModuleBase::timer::end("HSolverPW", "solve_psik"); return; } diff --git a/source/source_hsolver/hsolver_pw.h b/source/source_hsolver/hsolver_pw.h index db6764154c8..e2725794a70 100644 --- a/source/source_hsolver/hsolver_pw.h +++ b/source/source_hsolver/hsolver_pw.h @@ -74,7 +74,8 @@ class HSolverPW psi::Psi& psi, std::vector& pre_condition, Real* eigenvalue, - const int& nk_nums); + const int& nk_nums, + std::ostream& log); // calculate the precondition array for diagonalization in PW base void update_precondition(std::vector& h_diag, const int ik, const int npw, const Real vl_of_0); diff --git a/source/source_hsolver/unittests/CMakeLists.txt b/source/source_hsolver/unittests/CMakeLists.txt index b46d9bf1ad7..fefabbde345 100644 --- a/source/source_hsolver/unittests/CMakeLists.txt +++ b/source/source_hsolver/unittests/CMakeLists.txt @@ -66,9 +66,12 @@ if(ENABLE_MPI) ../hsolver_lcaopw.cpp ../diago_bpcg.cpp ../diago_dav_subspace.cpp + ../diago_ppcg.cpp ../diag_const_nums.cpp + ../diago_params.cpp ../diago_iter_assist.cpp ../para_lin_tf.cpp + ../../source_basis/module_pw/test/test_tool.cpp ../../source_estate/elecstate_tools.cpp ../../source_estate/occupy.cpp ../../source_base/module_fft/fft_bundle.cpp @@ -110,6 +113,19 @@ if(ENABLE_MPI) target_compile_definitions(MODULE_HSOLVER_LCAO_cusolver PRIVATE __CUDA) endif() endif() + +addtest( + TARGET MODULE_HSOLVER_ppcg + LIBS parameter base psi device container + SOURCES test_diago_ppcg.cpp ../diago_ppcg.cpp + ../../source_basis/module_pw/test/test_tool.cpp) + +addtest( + TARGET MODULE_HSOLVER_ppcg_float + LIBS parameter base psi device container + SOURCES test_diago_ppcg_float.cpp ../diago_ppcg.cpp + ../../source_basis/module_pw/test/test_tool.cpp) + abacus_copy_test_files( H-KPoints-Si2.dat H-GammaOnly-Si2.dat @@ -126,6 +142,7 @@ abacus_copy_test_files( diago_cg_parallel_test.sh diago_david_parallel_test.sh diago_lcao_parallel_test.sh + diago_ppcg_parallel_test.sh PEXSI-H-GammaOnly-Si2.dat PEXSI-S-GammaOnly-Si2.dat PEXSI-DM-GammaOnly-Si2.dat @@ -158,6 +175,12 @@ if(ENABLE_MPI) LIBS parameter base device MPI::MPI_CXX SOURCES test_para_lin_tf.cpp ../para_lin_tf.cpp) + addtest( + TARGET MODULE_HSOLVER_ppcg_parallel + LIBS parameter base psi device container MPI::MPI_CXX + SOURCES test_diago_ppcg_parallel.cpp ../diago_ppcg.cpp + ../../source_basis/module_pw/test/test_tool.cpp) + add_test( NAME MODULE_HSOLVER_para_linear_trans COMMAND mpirun -np 4 "$" @@ -182,6 +205,11 @@ if(ENABLE_MPI) ${BASH} diago_david_parallel_test.sh "$" "$" WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + add_test( + NAME MODULE_HSOLVER_ppcg_parallel_test + COMMAND ${BASH} diago_ppcg_parallel_test.sh + "$" + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) if(ENABLE_LCAO) add_test( NAME MODULE_HSOLVER_LCAO_parallel diff --git a/source/source_hsolver/unittests/diago_ppcg_parallel_test.sh b/source/source_hsolver/unittests/diago_ppcg_parallel_test.sh new file mode 100644 index 00000000000..6c2cd2625a0 --- /dev/null +++ b/source/source_hsolver/unittests/diago_ppcg_parallel_test.sh @@ -0,0 +1,44 @@ +#!/bin/bash + +ppcg_binary=${1:-./MODULE_HSOLVER_ppcg_parallel} +np=$(nproc --all) +echo "nprocs in this machine is $np" + +runs=0 +for i in 2 3 6; do + if [[ $i -gt $np ]]; then + continue + fi + echo "TEST DIAGO PPCG in parallel, nprocs=$i" + OMP_NUM_THREADS=1 mpirun -np $i "$ppcg_binary" + e=$? + if [[ $e -ne 0 ]]; then + echo -e "\e[1;33m [ FAILED ] \e[0m"\ + "execute UT with $i cores error." + exit 1 + fi + runs=$((runs + 1)) +done + +if [[ $np -ge 4 ]]; then + echo "TEST DIAGO PPCG with two pools, nprocs=4" + OMP_NUM_THREADS=1 mpirun -np 4 "$ppcg_binary" 2 + e=$? + if [[ $e -ne 0 ]]; then + echo -e "\e[1;33m [ FAILED ] \e[0m"\ + "execute two-pool UT with 4 cores error." + exit 1 + fi + runs=$((runs + 1)) +fi + +if [[ $runs -eq 0 ]]; then + echo "No 2/3/6-process case is supported; running the 1-process fallback." + OMP_NUM_THREADS=1 mpirun -np 1 "$ppcg_binary" + e=$? + if [[ $e -ne 0 ]]; then + echo -e "\e[1;33m [ FAILED ] \e[0m"\ + "execute UT with 1 core error." + exit 1 + fi +fi diff --git a/source/source_hsolver/unittests/test_diago_ppcg.cpp b/source/source_hsolver/unittests/test_diago_ppcg.cpp new file mode 100644 index 00000000000..073ce4a1a31 --- /dev/null +++ b/source/source_hsolver/unittests/test_diago_ppcg.cpp @@ -0,0 +1,2862 @@ +/** + * test_diago_ppcg.cpp — correctness and regression tests for DiagoPPCG + * + * Test matrices include both S = I and non-trivial overlap operators: + * 1. Tridiagonal Laplacian (1D particle-in-a-box): H[i,i]=2, H[i,i±1]=-1 + * Exact λ_k = 2 - 2·cos(k·π/(n+1)). Realistic but sparse. + * 2. Diagonal matrix: H = diag(1, 2, 3, 4, 5) + * Exact eigenvalues are the diagonal entries. Simplest possible + * smoke test — should converge in very few iterations. + * + * Analytical matrix cases share setup below; numerical failure regressions keep + * their explicit initial vectors and assertions. + */ + +#include "../diago_ppcg.h" +#include "source_base/timer.h" +#ifdef __MPI +#include "source_basis/module_pw/test/test_tool.h" +#endif + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +static constexpr double PPCG_PI = 3.14159265358979323846; + +using T = std::complex; +using Real = double; + +// ----------------------------------------------------------------------------- +// Helper: dense H-matrix times a set of column vectors +// H is stored column-major: H(row, col) = H_data[row + col * n_dim] +// ----------------------------------------------------------------------------- +static void dense_h_multiply(const T* H_data, int n_dim, const T* in, T* out, int ld, int ncol) +{ + for (int j = 0; j < ncol; ++j) + { + for (int i = 0; i < n_dim; ++i) + { + T sum = 0; + for (int k = 0; k < n_dim; ++k) + { + sum += H_data[i + k * n_dim] * in[k + j * ld]; + } + out[i + j * ld] = sum; + } + } +} + +// Use the same seeded real initial vectors for the analytical matrix fixtures. +static void initialize_vectors(std::vector& psi, int n_dim, int nband, int ld, unsigned seed) +{ + std::mt19937 rng(seed); + std::uniform_real_distribution dist(-1.0, 1.0); + const int size = ld * nband; + psi.assign(size, T(0)); + for (int j = 0; j < nband; ++j) + { + for (int i = 0; i < n_dim; ++i) + { + psi[i + j * ld] = T(dist(rng), 0.0); + } + } + for (int j = 0; j < nband; ++j) + { + for (int k = 0; k < j; ++k) + { + T dot = 0; + for (int i = 0; i < n_dim; ++i) + { + dot += std::conj(psi[i + k * ld]) * psi[i + j * ld]; + } + for (int i = 0; i < n_dim; ++i) + { + psi[i + j * ld] -= dot * psi[i + k * ld]; + } + } + Real norm_squared = 0; + for (int i = 0; i < n_dim; ++i) + { + norm_squared += std::norm(psi[i + j * ld]); + } + const Real norm = std::sqrt(norm_squared); + for (int i = 0; i < n_dim; ++i) + { + psi[i + j * ld] /= norm; + } + } +} + +// Shared storage and checks for analytical matrix cases. Regression tests below +// keep their own deliberately constructed vectors and operator callbacks. +class DiagoPPCGMatrixTest : public ::testing::Test +{ + protected: + void initialize(int dimension, int bands, unsigned seed, double threshold) + { + n_dim = dimension; + nband = bands; + ld = dimension; + ethr.assign(nband, threshold); + initialize_vectors(psi, n_dim, nband, ld, seed); + } + + void set_laplacian(double scale) + { + const int size = n_dim * n_dim; + H_mat.assign(size, T(0)); + for (int i = 0; i < n_dim; ++i) + { + H_mat[i + i * n_dim] = 2.0 * scale; + if (i > 0) + { + H_mat[i + (i - 1) * n_dim] = -scale; + H_mat[i - 1 + i * n_dim] = -scale; + } + } + const double diagonal = 2.0 * scale; + prec.assign(n_dim, diagonal); + exact.resize(nband); + for (int k = 0; k < nband; ++k) + { + const double angle = Real(k + 1) * PPCG_PI / Real(n_dim + 1); + exact[k] = scale * (2.0 - 2.0 * std::cos(angle)); + } + } + + void check_eigenvalues(double diag_threshold, int max_iterations, int block_size, + int ritz_interval, double tolerance) + { + std::vector eval(nband, 0.0); + hsolver::DiagoPPCG solver( + diag_threshold, max_iterations, block_size, ritz_interval); + const auto h_op = [this](T* in, T* out, int leading_dimension, int columns) { + dense_h_multiply(H_mat.data(), n_dim, in, out, leading_dimension, columns); + }; + const double iterations = solver.diag( + h_op, nullptr, ld, nband, n_dim, psi.data(), eval.data(), ethr, prec.data()); + for (int i = 0; i < nband; ++i) + { + EXPECT_NEAR(eval[i], exact[i], tolerance) << "eigenvalue " << i; + } + const double iteration_budget = max_iterations; + EXPECT_LE(iterations, iteration_budget); + } + + int n_dim; + int nband; + int ld; + std::vector H_mat; + std::vector prec; + std::vector exact; + std::vector ethr; + std::vector psi; +}; + +namespace hsolver +{ + +// Narrow friend access exercises cancellation in the projection helper without +// relying on a particular sequence of LAPACK rotations in an end-to-end solve. +class DiagoPPCGProjectionTest : public ::testing::Test +{ +protected: + void check_direction_rotation(const bool ritz) + { + const int dimension = 3; + const int leading_dimension = 4; + const int bands = 2; + const int size = leading_dimension * bands; + const std::vector h = {T(1), T(0, -0.4), T(0), + T(0, 0.4), T(4), T(0), T(0), T(0), T(10)}; + const std::vector s = {T(1), T(0), T(0), T(0), T(2), T(0), T(0), T(0), T(3)}; + std::vector x(size, T(0)); + x[0] = T(1); + x[leading_dimension] = T(0.3); + x[leading_dimension + 1] = T(1); + DiagoPPCG solver(1e-12, 20, 1, 4); + solver.ld_psi_ = leading_dimension; + solver.n_dim_ = dimension; + solver.n_band_ = bands; + solver.hpsi_.assign(size, T(0)); + solver.spsi_.assign(size, T(0)); + solver.w_.assign(size, T(0)); + solver.sw_.assign(size, T(0)); + solver.hw_.assign(size, T(0)); + solver.p_.assign(size, T(0)); + solver.p_[2] = T(0.4, 0.3); + solver.p_[leading_dimension + 2] = T(-0.2, 0.6); + solver.hp_.assign(size, T(0)); + solver.sp_.assign(size, T(0)); + solver.rr_eval_.resize(bands); + dense_h_multiply(h.data(), dimension, x.data(), solver.hpsi_.data(), leading_dimension, bands); + dense_h_multiply(s.data(), dimension, x.data(), solver.spsi_.data(), leading_dimension, bands); + dense_h_multiply(h.data(), dimension, solver.p_.data(), solver.hp_.data(), leading_dimension, bands); + dense_h_multiply(s.data(), dimension, solver.p_.data(), solver.sp_.data(), leading_dimension, bands); + const std::vector original_p = solver.p_; + if (ritz) + { + std::vector eigenvalues(bands); + const std::vector thresholds(bands, 1e-12); + std::vector active; + ASSERT_TRUE(solver.rayleigh_ritz(x.data(), eigenvalues.data(), active, thresholds)); + } + else + { + ASSERT_TRUE(solver.orthonormalize(x.data())); + } + double change = 0.0; + for (int column = 0; column < bands; ++column) + { + // Recover C from X_new = X_old C using the known triangular + // leading 2-by-2 block, independently of the solver's workspace. + const int first = column * leading_dimension; + const T c0 = x[first] - 0.3 * x[first + 1]; + const T c1 = x[first + 1]; + for (int row = 0; row < dimension; ++row) + { + const int index = first + row; + const T expected = original_p[row] * c0 + original_p[row + leading_dimension] * c1; + const T direction_difference = solver.p_[index] - expected; + const double error = std::abs(direction_difference); + EXPECT_LT(error, 1e-13); + const T direction_change = solver.p_[index] - original_p[index]; + change += std::norm(direction_change); + } + EXPECT_EQ(solver.p_[first + dimension], T(0)); + } + EXPECT_GT(change, 1e-3); + std::vector hp(size, T(0)); + std::vector sp(size, T(0)); + dense_h_multiply(h.data(), dimension, solver.p_.data(), hp.data(), leading_dimension, bands); + dense_h_multiply(s.data(), dimension, solver.p_.data(), sp.data(), leading_dimension, bands); + for (int index = 0; index < size; ++index) + { + const T h_difference = hp[index] - solver.hp_[index]; + const T s_difference = sp[index] - solver.sp_[index]; + const double h_error = std::abs(h_difference); + const double s_error = std::abs(s_difference); + EXPECT_LT(h_error, 1e-13); + EXPECT_LT(s_error, 1e-13); + } + } + + void check_rank_truncated_outputs(const int kept_count) + { + const int dimension = 8; + const int rank = 7; + const int matrix_size = dimension * dimension; + DiagoPPCG solver(1.0e-12, 20, kept_count, 1); + DiagoPPCG::SmallSubspace subspace; + subspace.k.assign(matrix_size, T(0)); + subspace.m.assign(matrix_size, T(0)); + subspace.eval.assign(dimension, 0.0); + for (int column = 0; column < rank; ++column) + { + const double weight = column + 1.0; + subspace.m[column + column * dimension] = T(weight); + subspace.k[column + column * dimension] = T(weight * weight); + if (column + 1 < rank) + { + const double scale = std::sqrt(weight * (weight + 1.0)); + const T coupling = T(-0.3, 0.2) * scale; + subspace.k[column + (column + 1) * dimension] = coupling; + subspace.k[column + 1 + column * dimension] = std::conj(coupling); + } + } + const std::vector hamiltonian = subspace.k; + const std::vector overlap = subspace.m; + const int valid_count = solver.solve_small_generalized(dimension, kept_count, subspace); + ASSERT_EQ(valid_count, rank); + // Check every exported vector, including those not used for the local + // update: the recovery path must receive valid, normalized Ritz pairs. + for (int column = 0; column < valid_count; ++column) + { + double residual_squared = 0.0; + for (int row = 0; row < dimension; ++row) + { + T residual(0); + for (int inner = 0; inner < dimension; ++inner) + { + const int offset = row + inner * dimension; + residual += (hamiltonian[offset] - subspace.eval[column] * overlap[offset]) + * subspace.k[inner + column * dimension]; + } + residual_squared += std::norm(residual); + } + const double residual_norm = std::sqrt(residual_squared); + EXPECT_LT(residual_norm, 1.0e-12); + for (int other = 0; other < valid_count; ++other) + { + T inner_product(0); + for (int row = 0; row < rank; ++row) + { + inner_product += std::conj(subspace.k[row + column * dimension]) + * (row + 1.0) * subspace.k[row + other * dimension]; + } + const T expected = column == other ? T(1) : T(0); + const double error = std::abs(inner_product - expected); + EXPECT_LT(error, 1.0e-12); + } + } + } + + void check_projection(const bool nearly_parallel, const bool zero_direction) + { + const int dimension = 2; + const int bands = 1; + const std::vector h_diagonal = {1.0, 1000.0}; + const std::vector s_diagonal = {1.0, 3.0}; + const std::vector x = {T(0.5), T(0.5)}; + const double scale = 1.0 / std::sqrt(12.0); + const std::vector z = {T(3.0 * scale), T(-scale)}; + DiagoPPCG solver(1.0e-12, 20, bands, 1); + solver.ld_psi_ = dimension; + solver.n_dim_ = dimension; + solver.n_band_ = bands; + solver.hpsi_.resize(dimension); + solver.spsi_.resize(dimension); + solver.p_.resize(dimension); + solver.hp_.resize(dimension); + solver.sp_.resize(dimension); + for (int i = 0; i < dimension; ++i) + { + solver.hpsi_[i] = h_diagonal[i] * x[i]; + solver.spsi_[i] = s_diagonal[i] * x[i]; + solver.p_[i] = z[i]; + if (nearly_parallel) + { + solver.p_[i] = x[i] + 1.0e-10 * z[i]; + } + solver.hp_[i] = h_diagonal[i] * solver.p_[i]; + solver.sp_[i] = s_diagonal[i] * solver.p_[i]; + if (zero_direction) + { + solver.p_[i] = T(0); + solver.hp_[i] = T(1.0e-15); + solver.sp_[i] = T(1.0e-15); + } + } + int h_calls = 0; + int s_calls = 0; + const auto h_op = [&h_calls, &h_diagonal, dimension](T* in, T* out, int ld, int ncol) { + ++h_calls; + for (int j = 0; j < ncol; ++j) + { + for (int i = 0; i < dimension; ++i) + { + out[i + j * ld] = h_diagonal[i] * in[i + j * ld]; + } + } + }; + solver.spsi_func_ = [&s_calls, &s_diagonal, dimension](T* in, T* out, int ld, int ncol) { + ++s_calls; + for (int j = 0; j < ncol; ++j) + { + for (int i = 0; i < dimension; ++i) + { + out[i + j * ld] = s_diagonal[i] * in[i + j * ld]; + } + } + }; + const std::vector columns = {0}; + solver.project_search_directions(h_op, x.data(), columns, columns); + double direction_norm_squared = 0.0; + double h_error_squared = 0.0; + double s_error_squared = 0.0; + for (int i = 0; i < dimension; ++i) + { + direction_norm_squared += std::norm(solver.p_[i]); + h_error_squared += std::norm(solver.hp_[i] - h_diagonal[i] * solver.p_[i]); + s_error_squared += std::norm(solver.sp_[i] - s_diagonal[i] * solver.p_[i]); + } + if (zero_direction) + { + EXPECT_DOUBLE_EQ(direction_norm_squared, 0.0); + EXPECT_DOUBLE_EQ(h_error_squared, 0.0); + EXPECT_DOUBLE_EQ(s_error_squared, 0.0); + EXPECT_EQ(h_calls, 0); + EXPECT_EQ(s_calls, 0); + return; + } + ASSERT_GT(direction_norm_squared, 0.0); + const double h_error = std::sqrt(h_error_squared / direction_norm_squared); + const double s_error = std::sqrt(s_error_squared / direction_norm_squared); + EXPECT_LT(h_error, 1.0e-12); + EXPECT_LT(s_error, 1.0e-12); + const int expected_calls = nearly_parallel ? 1 : 0; + EXPECT_EQ(h_calls, expected_calls); + EXPECT_EQ(s_calls, expected_calls); + } +}; + +TEST_F(DiagoPPCGProjectionTest, RefreshesImagesAfterSevereProjectionCancellation) +{ + check_projection(true, false); +} + +TEST_F(DiagoPPCGProjectionTest, RitzRotationPreservesDirectionColumnAssociation) +{ + check_direction_rotation(true); +} + +TEST_F(DiagoPPCGProjectionTest, OrthonormalizationRotatesDirectionImagesTogether) +{ + check_direction_rotation(false); +} + +TEST_F(DiagoPPCGProjectionTest, ReusesImagesForWellSeparatedDirection) +{ + check_projection(false, false); +} + +TEST_F(DiagoPPCGProjectionTest, ClearsStaleImagesForZeroDirection) +{ + check_projection(false, true); +} + +TEST_F(DiagoPPCGProjectionTest, RankTruncationExportsKeptAndDiscardedRitzVectors) +{ + const int kept_counts[] = {2, 3, 7}; + for (const int kept_count : kept_counts) + { + SCOPED_TRACE(kept_count); + check_rank_truncated_outputs(kept_count); + } +} + +} // namespace hsolver + +// ============================================================================= +// Test fixture: 1D particle-in-a-box (tridiagonal Laplacian) +// ============================================================================= +class DiagoPPCGTridiagTest : public DiagoPPCGMatrixTest +{ + protected: + void SetUp() override + { + initialize(10, 3, 42, 1e-10); + set_laplacian(1.0); + } +}; + +TEST_F(DiagoPPCGTridiagTest, BlockSubspace) +{ + std::vector psi_run = psi; + std::vector eval(nband, 0.0); + + hsolver::DiagoPPCG solver( + /* diag_thr = */ 1e-12, + /* max_iter = */ 100, + /* sbsize = */ 4, + /* rr_step = */ 4); + + auto h_op = [this](T* in, T* out, int ld_in, int ncol) { + dense_h_multiply(H_mat.data(), n_dim, in, out, ld_in, ncol); + }; + + double avg_iter = solver.diag(h_op, nullptr, ld, nband, n_dim, psi_run.data(), eval.data(), ethr, prec.data()); + + for (int i = 0; i < nband; ++i) + { + EXPECT_NEAR(eval[i], exact[i], 1e-8) << "Tridiag BLOCK: eigenvalue[" << i << "] mismatch"; + } + EXPECT_LE(avg_iter, double(100)) << "Tridiag BLOCK: too many iterations"; +} + +TEST_F(DiagoPPCGTridiagTest, ReportsIterationLimit) +{ + std::vector psi_run = psi; + std::vector eval(nband, 0.0); + std::vector strict_ethr(nband, 1.0e-14); + hsolver::DiagoPPCG solver( + /* diag_thr = */ 1.0e-14, + /* max_iter = */ 1, + /* sbsize = */ 3, + /* rr_step = */ 4); + + auto h_op = [this](T* in, T* out, int ld_in, int ncol) { + dense_h_multiply(H_mat.data(), n_dim, in, out, ld_in, ncol); + }; + + solver.diag(h_op, nullptr, ld, nband, n_dim, psi_run.data(), eval.data(), strict_ethr, prec.data()); + + EXPECT_FALSE(solver.converged()); + EXPECT_GT(solver.active_band_count(), 0); + EXPECT_TRUE(solver.iteration_limit_reached()); + EXPECT_GE(solver.first_active_band(), 0); + EXPECT_GT(solver.max_residual(), Real(1.0e-14)); +} + +TEST_F(DiagoPPCGTridiagTest, ReusesProjectedOperatorImages) +{ + std::vector psi_run = psi; + std::vector eval(nband, 0.0); + const std::vector strict_ethr(nband, 1.0e-14); + int h_calls = 0; + int s_calls = 0; + const auto h_op = [this, &h_calls](T* in, T* out, int ld_in, int ncol) { + ++h_calls; + dense_h_multiply(H_mat.data(), n_dim, in, out, ld_in, ncol); + }; + const auto s_op = [this, &s_calls](T* in, T* out, int ld_in, int ncol) { + ++s_calls; + for (int j = 0; j < ncol; ++j) + { + for (int i = 0; i < n_dim; ++i) + { + out[i + j * ld_in] = (1.0 + 0.1 * i) * in[i + j * ld_in]; + } + } + }; + hsolver::DiagoPPCG solver(1.0e-14, 2, 3, 16); + solver.diag(h_op, s_op, ld, nband, n_dim, + psi_run.data(), eval.data(), strict_ethr, prec.data()); + + ASSERT_TRUE(solver.iteration_limit_reached()); + // Initial/final X applications, one W application per sweep, and one + // cancellation-triggered P refresh. No unconditional HP or projected SW + // recomputation is needed, including when S is not I. + EXPECT_EQ(h_calls, 5); + EXPECT_EQ(s_calls, 5); + double max_residual = 0.0; + for (int j = 0; j < nband; ++j) + { + double residual_squared = 0.0; + for (int i = 0; i < n_dim; ++i) + { + T hpsi = T(0); + for (int k = 0; k < n_dim; ++k) + { + hpsi += H_mat[i + k * n_dim] * psi_run[k + j * ld]; + } + const T spsi = (1.0 + 0.1 * i) * psi_run[i + j * ld]; + residual_squared += std::norm(hpsi - eval[j] * spsi); + } + const Real residual_norm = std::sqrt(residual_squared); + max_residual = std::max(max_residual, residual_norm); + } + EXPECT_NEAR(solver.max_residual(), max_residual, 1.0e-13); +} + +TEST(DiagoPPCGReliableUpdateTest, RefreshesBeforeUsingNearPrecisionResidual) +{ + const int dimension = 2; + const std::vector prec(dimension, 1.0); + const std::vector thresholds(1, 1.0e-14); + const double perturbations[] = {1.0e-5, 1.0e-10}; + for (const double perturbation : perturbations) + { + SCOPED_TRACE(perturbation); + const double normalization = std::sqrt(1.0 + perturbation * perturbation); + std::vector psi = {T(1.0 / normalization), T(perturbation / normalization)}; + std::vector eigenvalues(1, 0.0); + int h_calls = 0; + const auto h_op = [&h_calls](T* in, T* out, int ld, int ncol) { + ++h_calls; + for (int j = 0; j < ncol; ++j) + { + out[j * ld] = in[j * ld]; + out[1 + j * ld] = 3.0 * in[1 + j * ld]; + } + }; + hsolver::DiagoPPCG solver( + 1.0e-14, 1, 1, 1); + solver.diag(h_op, nullptr, dimension, 1, dimension, + psi.data(), eigenvalues.data(), thresholds, prec.data()); + EXPECT_TRUE(solver.converged()); + EXPECT_LT(solver.max_residual(), 1.0e-14); + EXPECT_NEAR(eigenvalues[0], 1.0, 1.0e-14); + const T residual0 = (1.0 - eigenvalues[0]) * psi[0]; + const T residual1 = (3.0 - eigenvalues[0]) * psi[1]; + const double residual = std::sqrt(std::norm(residual0) + std::norm(residual1)); + EXPECT_LT(residual, 1.0e-14); + // Initial X, W, tentative-convergence validation, and final validation. + // Only the near-precision solve adds a reliable X/P refresh before W. + const int expected_calls = perturbation < 1.0e-8 ? 6 : 4; + EXPECT_EQ(h_calls, expected_calls); + } +} + +TEST(DiagoPPCGPeriodicRitzTest, HonorsIntervalAndChecksFinalRitzResidual) +{ + const int dimension = 48; + const int bands = 8; + std::mt19937 random(184); + std::uniform_real_distribution uniform(-1.0, 1.0); + std::vector initial(dimension * bands); + for (T& value : initial) + { + const double real = uniform(random); + const double imaginary = uniform(random); + value = T(real, imaginary); + } + const std::vector prec(dimension, 2.0); + const std::vector thresholds(bands, 1.0e-14); + struct Configuration + { + int sweeps; + int interval; + }; + // Ten sweeps need no periodic refresh; eighteen exercise the one at 16. + // Each configuration additionally encounters one projection cancellation. + // Each run encounters one competing discarded direction and + // therefore add one conditional expanded Ritz solve, without another H call. + const Configuration configurations[] = {{10, 1}, {10, 4}, {10, 16}, + {18, 1}, {18, 4}, {18, 16}}; + for (const Configuration& configuration : configurations) + { + const int sweeps = configuration.sweeps; + const int interval = configuration.interval; + SCOPED_TRACE(sweeps); + SCOPED_TRACE(interval); + std::vector psi = initial; + std::vector eigenvalues(bands, 0.0); + int applications = 0; + const auto h_op = [&applications, dimension](T* in, T* out, int ld, int ncol) { + ++applications; + for (int j = 0; j < ncol; ++j) + { + for (int i = 0; i < dimension; ++i) + { + T value = 2.0 * in[i + j * ld]; + if (i > 0) + { + value -= in[i - 1 + j * ld]; + } + if (i + 1 < dimension) + { + value -= in[i + 1 + j * ld]; + } + out[i + j * ld] = value; + } + } + }; + hsolver::DiagoPPCG solver( + 1.0e-14, sweeps, 2, interval); + solver.diag(h_op, nullptr, dimension, bands, dimension, + psi.data(), eigenvalues.data(), thresholds, prec.data()); + ASSERT_TRUE(solver.iteration_limit_reached()); + const int scheduled = sweeps / interval; + const int final_extra = sweeps % interval == 0 ? 0 : 1; + const int startup_extra = interval == 1 ? 0 : 1; + const int competing_extra = 1; + const int expected_ritz = 1 + scheduled + final_extra + startup_extra + competing_extra; + const int operator_refreshes = sweeps / 16; + const int cancellation_refreshes = 1; + const int expected_applications = 2 + sweeps + 2 * operator_refreshes + + cancellation_refreshes; + EXPECT_EQ(solver.root_recovery_count(), competing_extra); + EXPECT_EQ(solver.rayleigh_ritz_count(), expected_ritz); + EXPECT_EQ(applications, expected_applications); + EXPECT_TRUE(std::is_sorted(eigenvalues.begin(), eigenvalues.end())); + std::vector hpsi(psi.size()); + h_op(psi.data(), hpsi.data(), dimension, bands); + double residual = 0.0; + for (int j = 0; j < bands; ++j) + { + double norm_squared = 0.0; + for (int i = 0; i < dimension; ++i) + { + const T value = hpsi[i + j * dimension] - eigenvalues[j] * psi[i + j * dimension]; + norm_squared += std::norm(value); + } + const Real residual_norm = std::sqrt(norm_squared); + residual = std::max(residual, residual_norm); + for (int k = 0; k < bands; ++k) + { + T overlap = T(0); + for (int i = 0; i < dimension; ++i) + { + overlap += std::conj(psi[i + k * dimension]) * psi[i + j * dimension]; + } + const T expected = j == k ? T(1) : T(0); + const double error = std::abs(overlap - expected); + EXPECT_LT(error, 1.0e-11); + } + } + EXPECT_NEAR(solver.max_residual(), residual, 1.0e-13); + } +} + +TEST(DiagoPPCGPeriodicRitzTest, StartupReclassifiesDisplacedLockedPrefix) +{ + const int dimension = 5; + const int bands = 2; + std::vector psi(dimension * bands, T(0)); + // The initial exact eigenvalue 2 appears converged, but it is not the + // lowest root. The startup correction discovers a lower approximate + // root, which must replace and unlock this column before the next sweep. + psi[1] = T(1); + const double scale = 1.0 / std::sqrt(3.0); + psi[dimension] = T(scale); + psi[2 + dimension] = T(scale); + psi[4 + dimension] = T(scale); + std::vector eigenvalues(bands, 0.0); + const std::vector prec(dimension, 1.0); + const std::vector thresholds(bands, 1.0e-10); + const auto h_op = [dimension](T* in, T* out, int ld, int ncol) { + for (int j = 0; j < ncol; ++j) + { + for (int i = 0; i < dimension; ++i) + { + out[i + j * ld] = double(i + 1) * in[i + j * ld]; + } + } + }; + hsolver::DiagoPPCG solver(1.0e-12, 5, 1, 16); + solver.diag(h_op, nullptr, dimension, bands, dimension, + psi.data(), eigenvalues.data(), thresholds, prec.data()); + EXPECT_TRUE(solver.converged()); + EXPECT_NEAR(eigenvalues[0], 1.0, 1.0e-10); + EXPECT_NEAR(eigenvalues[1], 2.0, 1.0e-10); + EXPECT_LE(solver.max_residual(), 1.0e-10); +} + +TEST(DiagoPPCGPeriodicRitzTest, ClusteredGeneralizedSpectrumWithLockedPrefix) +{ + const int dimension = 16; + const int bands = 4; + std::vector spectrum(dimension); + spectrum[0] = -1.0; + spectrum[1] = 0.0; + spectrum[2] = 1.0e-6; + spectrum[3] = 1.0e-6; + spectrum[4] = 2.0e-6; + std::vector overlap(dimension); + std::vector prec(dimension, 1.0); + for (int i = 0; i < dimension; ++i) + { + overlap[i] = std::pow(10.0, 3.0 * i / (dimension - 1)); + if (i >= 5) + { + spectrum[i] = 1.0 + i; + } + prec[i] = overlap[i]; + } + std::mt19937 random(392); + std::uniform_real_distribution uniform(-1.0, 1.0); + std::vector initial(dimension * bands, T(0)); + initial[0] = T(1.0 / std::sqrt(overlap[0])); + for (int j = 1; j < bands; ++j) + { + for (int i = 1; i < dimension; ++i) + { + const double real = uniform(random); + const double imaginary = uniform(random); + initial[i + j * dimension] = T(real, imaginary) / std::sqrt(overlap[i]); + } + } + const auto h_op = [&spectrum, &overlap, dimension](T* in, T* out, int ld, int ncol) { + for (int j = 0; j < ncol; ++j) + { + for (int i = 0; i < dimension; ++i) + { + out[i + j * ld] = spectrum[i] * overlap[i] * in[i + j * ld]; + } + } + }; + const auto s_op = [&overlap, dimension](T* in, T* out, int ld, int ncol) { + for (int j = 0; j < ncol; ++j) + { + for (int i = 0; i < dimension; ++i) + { + out[i + j * ld] = overlap[i] * in[i + j * ld]; + } + } + }; + const std::vector thresholds(bands, 1.0e-10); + const int intervals[] = {1, 16}; + for (const int interval : intervals) + { + SCOPED_TRACE(interval); + std::vector psi = initial; + std::vector eigenvalues(bands, 0.0); + hsolver::DiagoPPCG solver( + 1.0e-12, 600, 2, interval); + solver.diag(h_op, s_op, dimension, bands, dimension, + psi.data(), eigenvalues.data(), thresholds, prec.data()); + EXPECT_TRUE(solver.converged()); + EXPECT_LE(solver.max_residual(), 1.0e-10); + for (int j = 0; j < bands; ++j) + { + EXPECT_NEAR(eigenvalues[j], spectrum[j], 1.0e-10); + for (int k = 0; k < bands; ++k) + { + T product = T(0); + for (int i = 0; i < dimension; ++i) + { + product += std::conj(psi[i + k * dimension]) * overlap[i] * psi[i + j * dimension]; + } + const T expected = j == k ? T(1) : T(0); + const double error = std::abs(product - expected); + EXPECT_LT(error, 1.0e-10); + } + } + const double prefix_weight = std::norm(psi[0]) * overlap[0]; + EXPECT_NEAR(prefix_weight, 1.0, 1.0e-12); + } +} + +TEST(DiagoPPCGStatusTest, ReportsActiveTailStart) +{ + const int n_dim = 4; + const int nband = 3; + const int ld = n_dim; + std::vector h_mat(n_dim * n_dim, T(0)); + for (int i = 0; i < n_dim; ++i) + { + h_mat[i + i * ld] = T(Real(i + 1), 0); + } + const Real inverse_sqrt_two = Real(1) / std::sqrt(Real(2)); + std::vector psi(ld * nband, T(0)); + psi[0] = T(1, 0); + psi[1 + ld] = T(1, 0); + psi[2 + 2 * ld] = T(inverse_sqrt_two, 0); + psi[3 + 2 * ld] = T(inverse_sqrt_two, 0); + std::vector eval(nband, Real(0)); + const std::vector ethr(nband, 1.0e-12); + const std::vector prec = {1, 2, 3, 4}; + const auto h_op = [&h_mat, n_dim](T* in, T* out, int ld_in, int ncol) { + dense_h_multiply(h_mat.data(), n_dim, in, out, ld_in, ncol); + }; + hsolver::DiagoPPCG solver( + /* diag_thr = */ 1.0e-12, + /* max_iter = */ 0, + /* sbsize = */ 3, + /* rr_step = */ 1); + + solver.diag(h_op, nullptr, ld, nband, n_dim, + psi.data(), eval.data(), ethr, prec.data()); + + EXPECT_FALSE(solver.converged()); + EXPECT_EQ(solver.first_active_band(), 2); + EXPECT_EQ(solver.active_band_count(), 1); + EXPECT_DOUBLE_EQ(solver.first_active_threshold(), Real(1.0e-12)); + EXPECT_GT(solver.max_residual(), solver.first_active_threshold()); +} + +TEST(DiagoPPCGStatusTest, NonFiniteFinalResidualIsNotConverged) +{ + const int n_dim = 2; + const int nband = 2; + const int ld = n_dim; + std::vector psi(ld * nband, T(0)); + psi[0] = T(1, 0); + psi[1 + ld] = T(1, 0); + std::vector eval(nband, Real(0)); + const std::vector ethr(nband, 1e-12); + const std::vector prec = {1, 2}; + int h_apply_count = 0; + const auto h_op = [&h_apply_count](T* in, T* out, int ld_in, int ncol) { + ++h_apply_count; + for (int j = 0; j < ncol; ++j) + { + for (int i = 0; i < ld_in; ++i) + { + const Real eigenvalue = Real(i + 1); + out[i + j * ld_in] = eigenvalue * in[i + j * ld_in]; + } + } + if (h_apply_count >= 2) + { + out[0] = T(std::numeric_limits::quiet_NaN(), 0); + } + }; + hsolver::DiagoPPCG solver( + /* diag_thr = */ 1e-12, + /* max_iter = */ 10, + /* sbsize = */ 2, + /* rr_step = */ 2); + + solver.diag(h_op, nullptr, ld, nband, n_dim, + psi.data(), eval.data(), ethr, prec.data()); + + EXPECT_FALSE(solver.converged()); + EXPECT_EQ(solver.active_band_count(), nband); + EXPECT_FALSE(solver.iteration_limit_reached()); + EXPECT_EQ(solver.first_active_band(), 0); + EXPECT_TRUE(std::isinf(solver.max_residual())); +} + +TEST(DiagoPPCGStatusTest, FailedRayleighRitzCannotCertifyLowestRoots) +{ + const int n_dim = 3; + const int nband = 2; + const int ld = n_dim; + std::vector h_mat(n_dim * n_dim, T(0)); + h_mat[0] = T(1, 0); + h_mat[1 + ld] = T(2, 0); + h_mat[2 + 2 * ld] = T(3, 0); + std::vector psi(ld * nband, T(0)); + psi[2] = T(1, 0); + psi[2 + ld] = T(1, 0); + std::vector eval(nband, Real(0)); + const std::vector ethr(nband, 1e-12); + const std::vector prec = {1, 2, 3}; + const auto h_op = [&h_mat, n_dim](T* in, T* out, int ld_in, int ncol) { + dense_h_multiply(h_mat.data(), n_dim, in, out, ld_in, ncol); + }; + hsolver::DiagoPPCG solver( + /* diag_thr = */ 1e-12, + /* max_iter = */ 0, + /* sbsize = */ 2, + /* rr_step = */ 2); + + solver.diag(h_op, nullptr, ld, nband, n_dim, + psi.data(), eval.data(), ethr, prec.data()); + + EXPECT_FALSE(solver.converged()); + EXPECT_EQ(solver.active_band_count(), nband); + EXPECT_TRUE(solver.iteration_limit_reached()); +} + +TEST_F(DiagoPPCGTridiagTest, ResidualTraceWritesCsv) +{ + const char* env_name = "ABACUS_PPCG_RESIDUAL_TRACE"; + const char* old_env = std::getenv(env_name); + const bool had_old_env = old_env != nullptr; + const std::string old_env_value = had_old_env ? old_env : ""; + const std::string trace_path = testing::TempDir() + "./ppcg_residual_trace_test.csv"; + std::remove(trace_path.c_str()); + + std::vector psi_run = psi; + std::vector eval(nband, 0.0); + hsolver::DiagoPPCG solver( + /* diag_thr = */ 1e-12, + /* max_iter = */ 100, + /* sbsize = */ 4, + /* rr_step = */ 4); + + auto h_op = [this](T* in, T* out, int ld_in, int ncol) { + dense_h_multiply(H_mat.data(), n_dim, in, out, ld_in, ncol); + }; + + std::ostringstream running_log; + auto& residual_timer = ModuleBase::timer::timer_pool["DiagoPPCG"]["max_generalized_residual"]; + const size_t before_disabled = residual_timer.calls; + ASSERT_EQ(::unsetenv(env_name), 0); + solver.diag(h_op, nullptr, ld, nband, n_dim, psi_run.data(), eval.data(), ethr, prec.data(), running_log); + const size_t disabled_calls = residual_timer.calls - before_disabled; + EXPECT_TRUE(solver.converged()); + + // An unset variable and a failed open must both skip trace-only work. + // Opening the existing temporary directory fails without creating a file. + const std::string directory_path = testing::TempDir(); + ASSERT_EQ(::setenv(env_name, directory_path.c_str(), 1), 0); + psi_run = psi; + const size_t before_failed_open = residual_timer.calls; + solver.diag(h_op, nullptr, ld, nband, n_dim, psi_run.data(), eval.data(), ethr, prec.data(), running_log); + const size_t failed_open_calls = residual_timer.calls - before_failed_open; + EXPECT_EQ(disabled_calls, failed_open_calls); + EXPECT_TRUE(solver.converged()); + EXPECT_TRUE(running_log.str().empty()); + + ASSERT_EQ(::setenv(env_name, trace_path.c_str(), 1), 0); + psi_run = psi; + const size_t before_enabled = residual_timer.calls; + solver.diag(h_op, nullptr, ld, nband, n_dim, psi_run.data(), eval.data(), ethr, prec.data(), running_log); + const size_t enabled_calls = residual_timer.calls - before_enabled; + EXPECT_GT(enabled_calls, disabled_calls); + + if (had_old_env) + { + ASSERT_EQ(::setenv(env_name, old_env_value.c_str(), 1), 0); + } + else + { + ASSERT_EQ(::unsetenv(env_name), 0); + } + + std::ifstream trace(trace_path); + ASSERT_TRUE(trace.good()); + char* resolved_path = ::realpath(trace_path.c_str(), nullptr); + ASSERT_NE(resolved_path, nullptr); + const std::string full_path(resolved_path); + std::free(resolved_path); + const std::string expected_log = "PPCG residual trace written: " + full_path; + EXPECT_NE(running_log.str().find(expected_log), std::string::npos); + std::string header; + std::string first_record; + std::getline(trace, header); + std::getline(trace, first_record); + EXPECT_EQ(header, + "iteration,stage,max_residual,eigenvalue_0,eigenvalue_1,eigenvalue_2," + "local_overlap_min,local_overlap_max,local_overlap_condition"); + EXPECT_NE(first_record.find("initial_rr"), std::string::npos); + EXPECT_EQ(std::count(first_record.begin(), first_record.end(), ','), nband + 5); + + double previous_max_error = std::numeric_limits::max(); + std::string record = first_record; + do + { + std::stringstream record_stream(record); + std::string field; + for (int field_index = 0; field_index < 3; ++field_index) + { + std::getline(record_stream, field, ','); + } + double max_error = 0.0; + for (int ib = 0; ib < nband; ++ib) + { + std::getline(record_stream, field, ','); + const double eigenvalue = std::stod(field); + const double difference = eigenvalue - exact[ib]; + const double error = std::abs(difference); + max_error = std::max(max_error, error); + } + const double allowed_error = previous_max_error + 1.0e-10; + EXPECT_LE(max_error, allowed_error); + previous_max_error = max_error; + } while (std::getline(trace, record)); + + trace.close(); + std::remove(trace_path.c_str()); +} + +// ============================================================================= +// Test fixture: diagonal matrix (simplest possible Hamiltonian) +// ============================================================================= +class DiagoPPCGDiagonalTest : public ::testing::Test +{ + protected: + void SetUp() override + { + n_dim = 5; + nband = 3; + ld = n_dim; + + // Build diagonal H: H[i,i] = i+1 + H_mat.assign(n_dim * n_dim, T(0)); + for (int i = 0; i < n_dim; ++i) + { + H_mat[i + i * n_dim] = T(Real(i + 1), 0); + } + + // Preconditioner — diagonal of H + prec.resize(n_dim); + for (int i = 0; i < n_dim; ++i) + { + prec[i] = Real(i + 1); + } + + // Lowest 3 eigenvalues: 1, 2, 3 + exact = {1.0, 2.0, 3.0}; + + // Convergence thresholds + ethr.assign(nband, 1e-10); + + // Random initial guess (fixed seed) + initialize_vectors(psi, n_dim, nband, ld, 42); + } + + int n_dim; + int nband; + int ld; + std::vector H_mat; + std::vector prec; + std::vector exact; + std::vector ethr; + std::vector psi; +}; + +TEST_F(DiagoPPCGDiagonalTest, BlockSubspace) +{ + std::vector psi_run = psi; + std::vector eval(nband, 0.0); + + hsolver::DiagoPPCG solver( + /* diag_thr = */ 1e-12, + /* max_iter = */ 50, + /* sbsize = */ 3, + /* rr_step = */ 3); + + auto h_op = [this](T* in, T* out, int ld_in, int ncol) { + dense_h_multiply(H_mat.data(), n_dim, in, out, ld_in, ncol); + }; + + double avg_iter = solver.diag(h_op, nullptr, ld, nband, n_dim, psi_run.data(), eval.data(), ethr, prec.data()); + + for (int i = 0; i < nband; ++i) + { + EXPECT_NEAR(eval[i], exact[i], 1e-8) << "Diagonal BLOCK: eigenvalue[" << i << "] mismatch"; + } + EXPECT_LE(avg_iter, double(50)) << "Diagonal BLOCK: too many iterations"; +} + +TEST(DiagoPPCGRootSelectionTest, StartupCorrectionKeepsLowestRoots) +{ + const int n_dim = 8; + const int nband = 4; + const int ld = n_dim; + const Real inv_sqrt_two = Real(1) / std::sqrt(Real(2)); + + std::vector h_mat(n_dim * n_dim, T(0)); + std::vector prec(n_dim); + for (int i = 0; i < n_dim; ++i) + { + const Real diagonal = Real(i + 1); + h_mat[i + i * n_dim] = T(diagonal, 0); + prec[i] = diagonal; + } + + // The first two columns contain the missing third and fourth roots. The + // last two are exact high-energy eigenvectors and would otherwise lock + // before their lower replacements can enter the active subspace. + std::vector psi(ld * nband, T(0)); + psi[0] = T(inv_sqrt_two, 0); + psi[2] = T(inv_sqrt_two, 0); + psi[1 + ld] = T(inv_sqrt_two, 0); + psi[3 + ld] = T(inv_sqrt_two, 0); + psi[6 + 2 * ld] = T(1, 0); + psi[7 + 3 * ld] = T(1, 0); + + std::vector eval(nband, Real(0)); + std::vector ethr(nband, 1e-12); + hsolver::DiagoPPCG solver( + /* diag_thr = */ 1e-12, + /* max_iter = */ 50, + /* sbsize = */ 2, + /* rr_step = */ 4); + + auto h_op = [&h_mat, n_dim](T* in, T* out, int ld_in, int ncol) { + dense_h_multiply(h_mat.data(), n_dim, in, out, ld_in, ncol); + }; + + const double avg_iter = solver.diag( + h_op, nullptr, ld, nband, n_dim, psi.data(), eval.data(), ethr, prec.data()); + for (int i = 0; i < nband; ++i) + { + EXPECT_NEAR(eval[i], Real(i + 1), 1e-10) << "root[" << i << "] mismatch"; + } + EXPECT_TRUE(solver.converged()); + EXPECT_LE(avg_iter, 50.0); +} + +TEST(DiagoPPCGRootSelectionTest, StartupCorrectionKeepsLowestRootsWithoutInitialGap) +{ + const int n_dim = 8; + const int nband = 4; + const int ld = n_dim; + const Real inv_sqrt_two = Real(1) / std::sqrt(Real(2)); + + std::vector h_mat(n_dim * n_dim, T(0)); + std::vector prec(n_dim); + for (int i = 0; i < n_dim; ++i) + { + const Real diagonal = Real(i + 1); + h_mat[i + i * n_dim] = T(diagonal, 0); + prec[i] = diagonal; + } + + // Every initial vector has a nonzero residual. Independent two-column + // blocks would nevertheless retain roots 1, 2, 5, and 6, discarding the + // third and fourth roots before a global comparison is possible. + std::vector psi(ld * nband, T(0)); + psi[0] = T(inv_sqrt_two, 0); + psi[2] = T(inv_sqrt_two, 0); + psi[1 + ld] = T(inv_sqrt_two, 0); + psi[3 + ld] = T(inv_sqrt_two, 0); + psi[4 + 2 * ld] = T(inv_sqrt_two, 0); + psi[6 + 2 * ld] = T(inv_sqrt_two, 0); + psi[5 + 3 * ld] = T(inv_sqrt_two, 0); + psi[7 + 3 * ld] = T(inv_sqrt_two, 0); + + std::vector eval(nband, Real(0)); + std::vector ethr(nband, 1e-12); + hsolver::DiagoPPCG solver( + /* diag_thr = */ 1e-12, + /* max_iter = */ 50, + /* sbsize = */ 2, + /* rr_step = */ 4); + + auto h_op = [&h_mat, n_dim](T* in, T* out, int ld_in, int ncol) { + dense_h_multiply(h_mat.data(), n_dim, in, out, ld_in, ncol); + }; + + const double avg_iter = solver.diag( + h_op, nullptr, ld, nband, n_dim, psi.data(), eval.data(), ethr, prec.data()); + for (int i = 0; i < nband; ++i) + { + EXPECT_NEAR(eval[i], Real(i + 1), 1e-10) << "root[" << i << "] mismatch"; + } + EXPECT_TRUE(solver.converged()); + EXPECT_LE(avg_iter, 50.0); +} + +TEST(DiagoPPCGRootSelectionTest, StartupCorrectionKeepsConvergedPrefixInactive) +{ + const int n_dim = 10; + const int nband = 6; + const int ld = n_dim; + const Real inv_sqrt_two = Real(1) / std::sqrt(Real(2)); + + std::vector h_mat(n_dim * n_dim, T(0)); + std::vector prec(n_dim); + for (int i = 0; i < n_dim; ++i) + { + const Real diagonal = Real(i + 1); + h_mat[i + i * n_dim] = T(diagonal, 0); + prec[i] = diagonal; + } + + // The first two roots are already converged. The active tail still needs + // one global comparison: its first two vectors contain roots 3--6, while + // its last two vectors are exact high-energy states. + std::vector psi(ld * nband, T(0)); + psi[0] = T(1, 0); + psi[1 + ld] = T(1, 0); + psi[2 + 2 * ld] = T(inv_sqrt_two, 0); + psi[4 + 2 * ld] = T(inv_sqrt_two, 0); + psi[3 + 3 * ld] = T(inv_sqrt_two, 0); + psi[5 + 3 * ld] = T(inv_sqrt_two, 0); + psi[8 + 4 * ld] = T(1, 0); + psi[9 + 5 * ld] = T(1, 0); + + std::vector eval(nband, Real(0)); + const std::vector ethr(nband, 1e-12); + hsolver::DiagoPPCG solver( + /* diag_thr = */ 1e-12, + /* max_iter = */ 1, + /* sbsize = */ 2, + /* rr_step = */ 4); + const auto h_op = [&h_mat, n_dim](T* in, T* out, int ld_in, int ncol) { + dense_h_multiply(h_mat.data(), n_dim, in, out, ld_in, ncol); + }; + + const double avg_iter = solver.diag( + h_op, nullptr, ld, nband, n_dim, psi.data(), eval.data(), ethr, prec.data()); + + for (int i = 0; i < nband; ++i) + { + EXPECT_NEAR(eval[i], Real(i + 1), 1e-10) << "root[" << i << "] mismatch"; + } + const double expected_avg_iter = 1.0 + 4.0 / 6.0; + EXPECT_TRUE(solver.converged()); + EXPECT_NEAR(avg_iter, expected_avg_iter, 1e-12); +} + +TEST(DiagoPPCGRootSelectionTest, WarmStartAfterHamiltonianChangeKeepsLowestRoots) +{ + const int n_dim = 8; + const int nband = 4; + const int ld = n_dim; + const Real inv_sqrt_two = Real(1) / std::sqrt(Real(2)); + + std::vector h_mat(n_dim * n_dim, T(0)); + std::vector prec(n_dim, Real(1)); + std::vector psi(ld * nband, T(0)); + const int first_indices[nband] = {0, 1, 4, 5}; + const int second_indices[nband] = {2, 3, 6, 7}; + for (int state = 0; state < nband; ++state) + { + const int first = first_indices[state]; + const int second = second_indices[state]; + const Real low_eigenvalue = Real(state + 1); + const Real high_eigenvalue = Real(state + 5); + const Real diagonal = Real(0.5) * (low_eigenvalue + high_eigenvalue); + const Real off_diagonal = Real(0.5) * (low_eigenvalue - high_eigenvalue); + h_mat[first + first * n_dim] = T(diagonal, 0); + h_mat[second + second * n_dim] = T(diagonal, 0); + h_mat[first + second * n_dim] = T(off_diagonal, 0); + h_mat[second + first * n_dim] = T(off_diagonal, 0); + psi[first + state * ld] = T(inv_sqrt_two, 0); + psi[second + state * ld] = T(inv_sqrt_two, 0); + prec[first] = diagonal; + prec[second] = diagonal; + } + + std::vector eval(nband, Real(0)); + std::vector ethr(nband, 1e-12); + hsolver::DiagoPPCG solver( + /* diag_thr = */ 1e-12, + /* max_iter = */ 50, + /* sbsize = */ 2, + /* rr_step = */ 4); + auto h_op = [&h_mat, n_dim](T* in, T* out, int ld_in, int ncol) { + dense_h_multiply(h_mat.data(), n_dim, in, out, ld_in, ncol); + }; + + solver.diag( + h_op, nullptr, ld, nband, n_dim, psi.data(), eval.data(), ethr, prec.data()); + ASSERT_TRUE(solver.converged()); + for (int i = 0; i < nband; ++i) + { + ASSERT_NEAR(eval[i], Real(i + 1), 1e-10); + } + + // Model a later SCF call: the previous eigenvectors are a valid warm start, + // but the changed Hamiltonian makes each one mix a low and a high root. + std::fill(h_mat.begin(), h_mat.end(), T(0)); + for (int i = 0; i < n_dim; ++i) + { + const Real diagonal = Real(i + 1); + h_mat[i + i * n_dim] = T(diagonal, 0); + prec[i] = diagonal; + } + + solver.diag( + h_op, nullptr, ld, nband, n_dim, psi.data(), eval.data(), ethr, prec.data()); + EXPECT_TRUE(solver.converged()); + for (int i = 0; i < nband; ++i) + { + EXPECT_NEAR(eval[i], Real(i + 1), 1e-10) << "root[" << i << "] mismatch"; + } +} + +TEST(DiagoPPCGRootSelectionTest, StartupCorrectionKeepsLowestGeneralizedRoots) +{ + const int n_dim = 8; + const int nband = 4; + const int ld = n_dim; + const Real inv_sqrt_two = Real(1) / std::sqrt(Real(2)); + + std::vector h_mat(n_dim * n_dim, T(0)); + std::vector s_mat(n_dim * n_dim, T(0)); + std::vector s_diagonal(n_dim); + std::vector prec(n_dim); + for (int i = 0; i < n_dim; ++i) + { + const Real eigenvalue = Real(i + 1); + const Real overlap_value = Real(1) + Real(0.1) * Real(i % 3); + const Real hamiltonian_value = eigenvalue * overlap_value; + h_mat[i + i * n_dim] = T(hamiltonian_value, 0); + s_mat[i + i * n_dim] = T(overlap_value, 0); + s_diagonal[i] = overlap_value; + prec[i] = hamiltonian_value; + } + + // H = S diag(1,...,8), so the generalized eigenvectors are the + // coordinate vectors. The initial columns create the same cross-block + // root-selection hazard as the identity-overlap regression above. + std::vector psi(ld * nband, T(0)); + std::vector normalized_coefficients(n_dim); + for (int i = 0; i < n_dim; ++i) + { + normalized_coefficients[i] = inv_sqrt_two / std::sqrt(s_diagonal[i]); + } + psi[0] = T(normalized_coefficients[0], 0); + psi[2] = T(normalized_coefficients[2], 0); + psi[1 + ld] = T(normalized_coefficients[1], 0); + psi[3 + ld] = T(normalized_coefficients[3], 0); + psi[4 + 2 * ld] = T(normalized_coefficients[4], 0); + psi[6 + 2 * ld] = T(normalized_coefficients[6], 0); + psi[5 + 3 * ld] = T(normalized_coefficients[5], 0); + psi[7 + 3 * ld] = T(normalized_coefficients[7], 0); + + std::vector eval(nband, Real(0)); + std::vector ethr(nband, 1e-12); + hsolver::DiagoPPCG solver( + /* diag_thr = */ 1e-12, + /* max_iter = */ 50, + /* sbsize = */ 2, + /* rr_step = */ 4); + + auto h_op = [&h_mat, n_dim](T* in, T* out, int ld_in, int ncol) { + dense_h_multiply(h_mat.data(), n_dim, in, out, ld_in, ncol); + }; + auto s_op = [&s_mat, n_dim](T* in, T* out, int ld_in, int ncol) { + dense_h_multiply(s_mat.data(), n_dim, in, out, ld_in, ncol); + }; + + const double avg_iter = solver.diag( + h_op, s_op, ld, nband, n_dim, psi.data(), eval.data(), ethr, prec.data()); + for (int i = 0; i < nband; ++i) + { + EXPECT_NEAR(eval[i], Real(i + 1), 1e-10) << "root[" << i << "] mismatch"; + } + EXPECT_TRUE(solver.converged()); + EXPECT_LE(avg_iter, 50.0); +} + +TEST(DiagoPPCGRootSelectionTest, LockingKeepsRequestedLowRoots) +{ + const int n_dim = 4; + const int nband = 2; + const int ld = n_dim; + std::vector h_mat(n_dim * n_dim, T(0)); + h_mat[0] = T(1, 0); + h_mat[1 + ld] = T(2, 0); + h_mat[2 + 2 * ld] = T(4, 0); + h_mat[3 + 3 * ld] = T(3, 0); + const std::vector prec = {1, 2, 4, 3}; + const std::vector ethr(nband, 1e-12); + const Real epsilons[] = {1e-4, 1e-3, 1e-2, 1e-1}; + + const auto h_op = [&h_mat, n_dim](T* in, T* out, int ld_in, int ncol) { + dense_h_multiply(h_mat.data(), n_dim, in, out, ld_in, ncol); + }; + + for (const Real epsilon : epsilons) + { + SCOPED_TRACE("epsilon = " + std::to_string(epsilon)); + std::vector psi(ld * nband, T(0)); + const Real norm = std::sqrt(Real(2) + epsilon * epsilon); + psi[0] = T(Real(1) / norm, 0); + psi[1] = T(epsilon / norm, 0); + psi[2] = T(Real(1) / norm, 0); + psi[3 + ld] = T(1, 0); + std::vector eval(nband, Real(0)); + + hsolver::DiagoPPCG solver( + /* diag_thr = */ 1e-12, + /* max_iter = */ 80, + /* sbsize = */ 2, + /* rr_step = */ 1); + const double avg_iter = solver.diag( + h_op, nullptr, ld, nband, n_dim, psi.data(), eval.data(), ethr, prec.data()); + EXPECT_NEAR(eval[0], 1.0, 1e-8); + EXPECT_NEAR(eval[1], 2.0, 1e-8); + EXPECT_TRUE(solver.converged()); + EXPECT_LE(avg_iter, 80.0); + } +} + +TEST(DiagoPPCGRootSelectionTest, LaterCorrectionKeepsLowestRootsAcrossBlocks) +{ + const int n_dim = 4; + const int nband = 2; + const int ld = n_dim; + const std::vector diagonal = {1.0, 2.0, 3.0, 100.0}; + const std::vector prec(n_dim, 1.0); + const std::vector ethr(nband, 1.0e-10); + const Real coefficient = 1.0 / std::sqrt(3.0); + const int block_sizes[] = {1, 2}; + const auto h_op = [&diagonal, n_dim](T* in, T* out, int ld_in, int ncol) { + for (int column = 0; column < ncol; ++column) + { + for (int row = 0; row < n_dim; ++row) + { + out[row + column * ld_in] = diagonal[row] * in[row + column * ld_in]; + } + } + }; + + for (const int block_size : block_sizes) + { + SCOPED_TRACE("block_size = " + std::to_string(block_size)); + // The first expanded subspace keeps a low mixture and the exact + // third root. The missing second root needs a later correction; + // convergence of the retained vectors alone does not certify it. + std::vector psi(ld * nband, T(0)); + psi[0] = T(coefficient); + psi[1] = T(coefficient); + psi[3] = T(coefficient); + psi[2 + ld] = T(1); + std::vector eval(nband, 0.0); + hsolver::DiagoPPCG solver( + 1.0e-10, 100, block_size, 1); + solver.diag(h_op, nullptr, ld, nband, n_dim, + psi.data(), eval.data(), ethr, prec.data()); + + double max_residual = 0.0; + for (int column = 0; column < nband; ++column) + { + double residual_squared = 0.0; + for (int row = 0; row < n_dim; ++row) + { + const T residual = (diagonal[row] - eval[column]) * psi[row + column * ld]; + residual_squared += std::norm(residual); + } + const double residual_norm = std::sqrt(residual_squared); + max_residual = std::max(max_residual, residual_norm); + for (int other = 0; other < nband; ++other) + { + T overlap(0); + for (int row = 0; row < n_dim; ++row) + { + overlap += std::conj(psi[row + column * ld]) * psi[row + other * ld]; + } + const T expected_overlap = column == other ? T(1) : T(0); + const double overlap_error = std::abs(overlap - expected_overlap); + EXPECT_LE(overlap_error, 1.0e-12); + } + } + EXPECT_TRUE(solver.converged()); + EXPECT_LE(max_residual, 1.0e-10); + EXPECT_NEAR(solver.max_residual(), max_residual, 1.0e-13); + EXPECT_NEAR(eval[0], 1.0, 1.0e-9); + EXPECT_NEAR(eval[1], 2.0, 1.0e-9) + << "returned roots = " << eval[0] << ", " << eval[1] + << "; fresh max residual = " << max_residual; + if (block_size == 1) + { + EXPECT_GT(solver.root_recovery_count(), 0); + } + else + { + EXPECT_EQ(solver.root_recovery_count(), 0); + } + } +} + +TEST(DiagoPPCGRootSelectionTest, LaterCorrectionsRecoverMultipleComplexGeneralizedRoots) +{ + const int dimension = 6; + const int bands = 3; + const Real spectrum[] = {1.0, 2.0, 3.0, 4.0, 5.0, 100.0}; + const Real overlap_spectrum[] = {1.0, 2.0, 4.0, 7.0, 3.0, 5.0}; + std::vector unitary(dimension * dimension); + const Real normalization = 1.0 / std::sqrt(Real(dimension)); + for (int column = 0; column < dimension; ++column) + { + for (int row = 0; row < dimension; ++row) + { + const Real angle = 2.0 * PPCG_PI * row * column / dimension; + const T phase(std::cos(angle), std::sin(angle)); + unitary[row + column * dimension] = normalization * phase; + } + } + std::vector hamiltonian(dimension * dimension, T(0)); + std::vector overlap(dimension * dimension, T(0)); + for (int column = 0; column < dimension; ++column) + { + for (int row = 0; row < dimension; ++row) + { + for (int state = 0; state < dimension; ++state) + { + const T projector = unitary[row + state * dimension] + * std::conj(unitary[column + state * dimension]); + hamiltonian[row + column * dimension] += projector * overlap_spectrum[state] * spectrum[state]; + overlap[row + column * dimension] += projector * overlap_spectrum[state]; + } + } + } + const auto h_op = [&hamiltonian, dimension](T* in, T* out, int ld, int count) { + dense_h_multiply(hamiltonian.data(), dimension, in, out, ld, count); + }; + const auto s_op = [&overlap, dimension](T* in, T* out, int ld, int count) { + dense_h_multiply(overlap.data(), dimension, in, out, ld, count); + }; + const std::vector preconditioner(dimension, 1.0); + const std::vector thresholds(bands, 1.0e-9); + const int block_sizes[] = {1, 2, 3}; + const int ritz_intervals[] = {1, 4}; + for (const int block_size : block_sizes) + { + for (const int ritz_interval : ritz_intervals) + { + SCOPED_TRACE("block_size = " + std::to_string(block_size) + + ", ritz_interval = " + std::to_string(ritz_interval)); + std::vector psi(dimension * bands, T(0)); + // Four spectral components in the first vector require more + // than the first [X,W] expansion to resolve all three low roots. + const int mixed_states[] = {0, 1, 2, 5}; + for (int row = 0; row < dimension; ++row) + { + for (const int state : mixed_states) + { + psi[row] += 0.5 * unitary[row + state * dimension] + / std::sqrt(overlap_spectrum[state]); + } + psi[row + dimension] = unitary[row + 3 * dimension] / std::sqrt(overlap_spectrum[3]); + psi[row + 2 * dimension] = unitary[row + 4 * dimension] / std::sqrt(overlap_spectrum[4]); + } + std::vector eigenvalues(bands, 0.0); + hsolver::DiagoPPCG solver( + 1.0e-9, 120, block_size, ritz_interval); + solver.diag(h_op, s_op, dimension, bands, dimension, + psi.data(), eigenvalues.data(), thresholds, preconditioner.data()); + std::vector hpsi(psi.size(), T(0)); + std::vector spsi(psi.size(), T(0)); + h_op(psi.data(), hpsi.data(), dimension, bands); + s_op(psi.data(), spsi.data(), dimension, bands); + EXPECT_TRUE(solver.converged()); + for (int column = 0; column < bands; ++column) + { + EXPECT_NEAR(eigenvalues[column], spectrum[column], 1.0e-8); + double residual_squared = 0.0; + for (int row = 0; row < dimension; ++row) + { + const int offset = row + column * dimension; + const T residual = hpsi[offset] - eigenvalues[column] * spsi[offset]; + residual_squared += std::norm(residual); + } + const double residual_norm = std::sqrt(residual_squared); + EXPECT_LE(residual_norm, 1.0e-9); + for (int other = 0; other < bands; ++other) + { + T inner_product(0); + for (int row = 0; row < dimension; ++row) + { + inner_product += std::conj(psi[row + column * dimension]) * spsi[row + other * dimension]; + } + const T expected = column == other ? T(1) : T(0); + const double orthogonality_error = std::abs(inner_product - expected); + EXPECT_LE(orthogonality_error, 1.0e-11); + } + } + } + } +} + +TEST(DiagoPPCGRankTruncationTest, ComplexHermitianStartupSubspace) +{ + const int n_dim = 6; + const int nband = 3; + const int ld = n_dim; + const Real inv_sqrt_two = Real(1) / std::sqrt(Real(2)); + + // Build H = U diag(1,...,6) U^H. The first two columns of U form a + // complex rotation, so the reduced Hamiltonian needs a conjugate + // transpose when it is projected into the retained overlap subspace. + std::vector unitary(n_dim * n_dim, T(0)); + unitary[0] = T(inv_sqrt_two, 0); + unitary[1] = T(0, inv_sqrt_two); + unitary[n_dim] = T(0, inv_sqrt_two); + unitary[1 + n_dim] = T(inv_sqrt_two, 0); + for (int i = 2; i < n_dim; ++i) + { + unitary[i + i * n_dim] = T(1, 0); + } + + std::vector h_mat(n_dim * n_dim, T(0)); + std::vector eigenvalues(n_dim); + for (int state = 0; state < n_dim; ++state) + { + eigenvalues[state] = Real(state + 1); + for (int col = 0; col < n_dim; ++col) + { + for (int row = 0; row < n_dim; ++row) + { + h_mat[row + col * n_dim] += unitary[row + state * n_dim] + * eigenvalues[state] + * std::conj(unitary[col + state * n_dim]); + } + } + } + EXPECT_GT(std::abs(std::imag(h_mat[n_dim])), Real(0.1)); + + std::vector prec(n_dim); + for (int i = 0; i < n_dim; ++i) + { + const Real diagonal_value = std::real(h_mat[i + i * n_dim]); + prec[i] = std::max(Real(1), diagonal_value); + } + + // The third column is an exact high-energy state. Its zero residual + // makes [psi, w] rank deficient during the first global correction, + // while the first two columns and their residuals span the three lowest + // roots. This exercises Q^H K Q for a complex Hermitian matrix. + std::vector psi(ld * nband, T(0)); + for (int i = 0; i < n_dim; ++i) + { + psi[i] = (unitary[i] + unitary[i + 3 * n_dim]) * inv_sqrt_two; + psi[i + ld] = (unitary[i + n_dim] + unitary[i + 2 * n_dim]) * inv_sqrt_two; + psi[i + 2 * ld] = unitary[i + 5 * n_dim]; + } + + std::vector eval(nband, Real(0)); + std::vector ethr(nband, 1e-12); + hsolver::DiagoPPCG solver( + /* diag_thr = */ 1e-12, + /* max_iter = */ 30, + /* sbsize = */ 2, + /* rr_step = */ 4); + + auto h_op = [&h_mat, n_dim](T* in, T* out, int ld_in, int ncol) { + dense_h_multiply(h_mat.data(), n_dim, in, out, ld_in, ncol); + }; + + const double avg_iter = solver.diag( + h_op, nullptr, ld, nband, n_dim, psi.data(), eval.data(), ethr, prec.data()); + for (int i = 0; i < nband; ++i) + { + EXPECT_NEAR(eval[i], Real(i + 1), 1e-10) << "root[" << i << "] mismatch"; + } + EXPECT_TRUE(solver.converged()); + EXPECT_LE(avg_iter, 30.0); +} + +TEST_F(DiagoPPCGDiagonalTest, EmptyHOperatorThrows) +{ + std::vector psi_run = psi; + std::vector eval(nband, 0.0); + + hsolver::DiagoPPCG solver( + /* diag_thr = */ 1e-12, + /* max_iter = */ 50, + /* sbsize = */ 3, + /* rr_step = */ 3); + + hsolver::DiagoPPCG::HPsiFunc h_op; + EXPECT_THROW(solver.diag(h_op, nullptr, ld, nband, n_dim, psi_run.data(), eval.data(), ethr, prec.data()), + std::invalid_argument); +} + +TEST_F(DiagoPPCGDiagonalTest, NonFiniteInputThrows) +{ + std::vector psi_run = psi; + std::vector eval(nband, 0.0); + + hsolver::DiagoPPCG solver( + /* diag_thr = */ 1e-12, + /* max_iter = */ 50, + /* sbsize = */ 3, + /* rr_step = */ 3); + + auto h_op = [this](T* in, T* out, int ld_in, int ncol) { + dense_h_multiply(H_mat.data(), n_dim, in, out, ld_in, ncol); + }; + + std::vector bad_ethr = ethr; + bad_ethr[0] = std::numeric_limits::quiet_NaN(); + EXPECT_THROW(solver.diag(h_op, nullptr, ld, nband, n_dim, psi_run.data(), eval.data(), bad_ethr, prec.data()), + std::invalid_argument); + + std::vector bad_prec = prec; + bad_prec[0] = std::numeric_limits::infinity(); + EXPECT_THROW(solver.diag(h_op, nullptr, ld, nband, n_dim, psi_run.data(), eval.data(), ethr, bad_prec.data()), + std::invalid_argument); +} + +TEST(DiagoPPCGWorkspaceTest, RejectsOversizedWorkspacesBeforeUsingInputBuffers) +{ + struct Shape + { + const char* description; + int leading_dimension; + int bands; + }; + const int int_max = std::numeric_limits::max(); + const Shape shapes[] = { + {"wavefunction storage", int_max / 2 + 1, 2}, + {"local basis storage", int_max / 3 + 1, 2}, + {"projected matrix storage", 1, 32768}, + {"LAPACK workspace", 1, 20000}, + {"root recovery workspace", 1, 12000}, + }; + + for (const Shape& shape : shapes) + { + SCOPED_TRACE(shape.description); + T psi_value(1.0, 0.0); + Real eigenvalue = -1.0; + const Real preconditioner = 1.0; + // The size diagnostic must precede threshold/buffer validation. An + // empty threshold vector also keeps this regression bounded if that + // diagnostic is accidentally removed: no huge allocations are needed. + const std::vector thresholds; + int h_calls = 0; + int s_calls = 0; + const auto h_op = [&h_calls](T*, T*, int, int) { ++h_calls; }; + const auto s_op = [&s_calls](T*, T*, int, int) { ++s_calls; }; + hsolver::DiagoPPCG solver(1e-12, 1, 1, 1); + + try + { + solver.diag(h_op, + s_op, + shape.leading_dimension, + shape.bands, + 1, + &psi_value, + &eigenvalue, + thresholds, + &preconditioner); + FAIL() << "Oversized PPCG workspace was accepted"; + } + catch (const std::invalid_argument& error) + { + EXPECT_STREQ(error.what(), "PPCG: workspace dimensions exceed integer indexing or LAPACK limits."); + } + EXPECT_EQ(h_calls, 0); + EXPECT_EQ(s_calls, 0); + EXPECT_EQ(psi_value, T(1.0, 0.0)); + EXPECT_DOUBLE_EQ(eigenvalue, -1.0); + } +} + +TEST_F(DiagoPPCGDiagonalTest, HandlesSubblockSizeAboveBandCount) +{ + std::vector psi_run = psi; + std::vector eval(nband, 0.0); + const int subblock_size = std::numeric_limits::max(); + hsolver::DiagoPPCG solver(1e-12, 100, subblock_size, 3); + const auto h_op = [this](T* in, T* out, int ld_in, int ncol) { + dense_h_multiply(H_mat.data(), n_dim, in, out, ld_in, ncol); + }; + + solver.diag(h_op, nullptr, ld, nband, n_dim, psi_run.data(), eval.data(), ethr, prec.data()); + + EXPECT_TRUE(solver.converged()); + for (int band = 0; band < nband; ++band) + { + EXPECT_NEAR(eval[band], exact[band], ethr[band]); + } +} + +TEST(DiagoPPCGLeadingDimensionTest, BlockSubspaceWithPadding) +{ + const int n_dim = 5; + const int nband = 3; + const int ld = 8; + + std::vector H_mat(n_dim * n_dim, T(0)); + for (int i = 0; i < n_dim; ++i) + { + H_mat[i + i * n_dim] = T(Real(i + 1), 0); + } + + std::vector prec(n_dim); + for (int i = 0; i < n_dim; ++i) + { + prec[i] = Real(i + 1); + } + + std::vector psi(ld * nband, T(17.0, -3.0)); + std::mt19937 rng(7); + std::uniform_real_distribution dist(-1.0, 1.0); + for (int j = 0; j < nband; ++j) + { + for (int i = 0; i < n_dim; ++i) + { + psi[i + j * ld] = T(dist(rng), 0.0); + } + } + + for (int j = 0; j < nband; ++j) + { + for (int k = 0; k < j; ++k) + { + T dot = 0; + for (int i = 0; i < n_dim; ++i) + { + dot += std::conj(psi[i + k * ld]) * psi[i + j * ld]; + } + for (int i = 0; i < n_dim; ++i) + { + psi[i + j * ld] -= dot * psi[i + k * ld]; + } + } + Real nrm = 0; + for (int i = 0; i < n_dim; ++i) + { + nrm += std::norm(psi[i + j * ld]); + } + nrm = std::sqrt(nrm); + for (int i = 0; i < n_dim; ++i) + { + psi[i + j * ld] /= nrm; + } + } + + std::vector eval(nband, 0.0); + std::vector ethr(nband, 1e-10); + hsolver::DiagoPPCG solver( + /* diag_thr = */ 1e-12, + /* max_iter = */ 80, + /* sbsize = */ 2, + /* rr_step = */ 3); + + auto h_op = [&H_mat, n_dim](T* in, T* out, int ld_in, int ncol) { + dense_h_multiply(H_mat.data(), n_dim, in, out, ld_in, ncol); + }; + + double avg_iter = solver.diag(h_op, nullptr, ld, nband, n_dim, psi.data(), eval.data(), ethr, prec.data()); + + const Real exact[] = {1.0, 2.0, 3.0}; + for (int i = 0; i < nband; ++i) + { + EXPECT_NEAR(eval[i], exact[i], 1e-8) << "Padded ld BLOCK: eigenvalue[" << i << "] mismatch"; + } + EXPECT_LE(avg_iter, double(80)) << "Padded ld BLOCK: too many iterations"; +} + +// ============================================================================= +// Test fixture: 2×2 matrix — smallest non-trivial case +// H = [[2, 1], [1, 2]], eigenvalues: 1, 3 +// ============================================================================= +TEST_F(DiagoPPCGMatrixTest, TwoByTwo) +{ + initialize(2, 2, 123, 1e-10); + // H = [[2, 1], [1, 2]] + H_mat.assign(n_dim * n_dim, T(0)); + H_mat[0 + 0 * n_dim] = T(2.0, 0); + H_mat[1 + 1 * n_dim] = T(2.0, 0); + H_mat[0 + 1 * n_dim] = T(1.0, 0); + H_mat[1 + 0 * n_dim] = T(1.0, 0); + + prec.assign(n_dim, 2.0); + + // λ₁ = 1, λ₂ = 3 + exact = {1.0, 3.0}; + check_eigenvalues(1e-12, 50, 2, 2, 1e-8); +} + +TEST(DiagoPPCGComplexHermitianTest, DefaultKeepsImaginaryProjection) +{ + const int n_dim = 2; + const int nband = 2; + const int ld = n_dim; + + // H = [[2, i], [-i, 3]]. Dropping Im() would incorrectly + // diagonalize diag(2, 3); the Hermitian eigenvalues are 2.5 +/- sqrt(1.25). + std::vector H_mat(n_dim * n_dim, T(0)); + H_mat[0 + 0 * n_dim] = T(2.0, 0.0); + H_mat[1 + 1 * n_dim] = T(3.0, 0.0); + H_mat[0 + 1 * n_dim] = T(0.0, 1.0); + H_mat[1 + 0 * n_dim] = T(0.0, -1.0); + + std::vector psi(ld * nband, T(0)); + psi[0 + 0 * ld] = T(1.0, 0.0); + psi[1 + 1 * ld] = T(1.0, 0.0); + + std::vector prec(n_dim, 2.0); + std::vector ethr(nband, 1e-12); + std::vector eval(nband, 0.0); + + hsolver::DiagoPPCG solver( + /* diag_thr = */ 1e-12, + /* max_iter = */ 10, + /* sbsize = */ 2, + /* rr_step = */ 1); + + auto h_op = [&H_mat, n_dim](T* in, T* out, int ld_in, int ncol) { + dense_h_multiply(H_mat.data(), n_dim, in, out, ld_in, ncol); + }; + + solver.diag(h_op, nullptr, ld, nband, n_dim, psi.data(), eval.data(), ethr, prec.data()); + + EXPECT_TRUE(std::isfinite(eval[0])); + EXPECT_TRUE(std::isfinite(eval[1])); + const Real delta = std::sqrt(1.25); + EXPECT_NEAR(eval[0], 2.5 - delta, 1e-10); + EXPECT_NEAR(eval[1], 2.5 + delta, 1e-10); +} + +TEST(DiagoPPCGComplexHermitianTest, ExactImaginaryLeadingCoefficientNeedsNoProjection) +{ + const int n_dim = 3; + const int nband = 1; + const int ld = n_dim; + const int matrix_size = n_dim * n_dim; + const Real diagonal[n_dim] = {Real(1), Real(4), Real(9)}; + std::vector h_mat(matrix_size, T(0)); + for (int i = 0; i < n_dim; ++i) + { + h_mat[i + i * n_dim] = T(diagonal[i], 0); + } + + // A valid full-complex eigenvector may have an imaginary first + // coefficient. Clearing that coefficient would destroy its norm and + // rank. A Ritz solve is free to choose another overall phase, so check + // the eigenpair and normalization rather than a particular phase. + std::vector psi(ld, T(0)); + psi[0] = T(0, 1); + std::vector eval(nband, Real(0)); + const std::vector ethr(nband, 1e-12); + const std::vector prec(n_dim, Real(1)); + const auto h_op = [&h_mat, n_dim](T* in, T* out, int ld_in, int ncol) { + dense_h_multiply(h_mat.data(), n_dim, in, out, ld_in, ncol); + }; + hsolver::DiagoPPCG solver( + /* diag_thr = */ 1e-12, + /* max_iter = */ 0, + /* sbsize = */ 1, + /* rr_step = */ 4); + + solver.diag( + h_op, nullptr, ld, nband, n_dim, psi.data(), eval.data(), ethr, prec.data()); + + Real norm_squared = Real(0); + Real residual_squared = Real(0); + for (int i = 0; i < n_dim; ++i) + { + const T residual = diagonal[i] * psi[i] - eval[0] * psi[i]; + norm_squared += std::norm(psi[i]); + residual_squared += std::norm(residual); + } + const Real residual_norm = std::sqrt(residual_squared); + EXPECT_TRUE(solver.converged()); + EXPECT_FALSE(solver.iteration_limit_reached()); + EXPECT_NEAR(eval[0], diagonal[0], 1e-12); + EXPECT_NEAR(norm_squared, Real(1), 1e-12); + EXPECT_LE(residual_norm, Real(1e-12)); +} + +// ============================================================================= +// Test fixture: degenerate eigenvalues +// H = I + J (identity plus all-ones), 4×4. +// J has eigenvector [1,1,1,1]^T with eigenvalue 4. +// All vectors orthogonal to [1,1,1,1]^T are eigenvectors with eigenvalue 0. +// So H = I + J has: λ₁ = 1 (multiplicity 3), λ₄ = 5. +// ============================================================================= +TEST_F(DiagoPPCGMatrixTest, DegenerateEigenvalues) +{ + initialize(4, 4, 456, 1e-10); + // H = I + J where J is the all-ones matrix + H_mat.assign(n_dim * n_dim, T(0)); + for (int i = 0; i < n_dim; ++i) + { + for (int j = 0; j < n_dim; ++j) + { + H_mat[i + j * n_dim] = T(1.0, 0); // all-ones J + } + } + for (int i = 0; i < n_dim; ++i) + { + H_mat[i + i * n_dim] += T(1.0, 0); // J → I+J + } + + // Preconditioner: diagonal = 2 + prec.assign(n_dim, 2.0); + + // λ = {1, 1, 1, 5} + exact = {1.0, 1.0, 1.0, 5.0}; + check_eigenvalues(1e-12, 100, 4, 4, 1e-8); +} + +// ============================================================================= +// Test fixture: larger tridiagonal, more bands +// 20×20 tridiagonal Laplacian, nband=5. +// ============================================================================= +TEST_F(DiagoPPCGMatrixTest, Laplacian20) +{ + initialize(20, 5, 789, 1e-10); + set_laplacian(1.0); + check_eigenvalues(1e-12, 150, 5, 5, 1e-8); +} + +// ============================================================================= +// Test fixture: dense matrix with known eigenvalues +// H = Q^T * D * Q where Q is a known orthogonal matrix (a Givens rotation +// repeated on different index pairs) and D is diagonal. +// For an 8×8 case: D = diag(1, 2, 3, 4, 5, 6, 7, 8), then apply several +// Givens rotations to mix all rows/cols. The exact eigenvalues remain 1..8. +// ============================================================================= +TEST_F(DiagoPPCGMatrixTest, DenseKnownSpectrum) +{ + initialize(8, 4, 222, 1e-10); + // Start with diagonal matrix + std::vector dense(n_dim * n_dim, Real(0)); + for (int i = 0; i < n_dim; ++i) + { + dense[i + i * n_dim] = Real(i + 1); + } + + // Apply several Givens rotations to make it dense while preserving + // eigenvalues. Each rotation: A' = G(i,j,θ)^T * A * G(i,j,θ) + auto apply_givens = [&](int p, int q, Real theta) { + Real c = std::cos(theta); + Real s = std::sin(theta); + // Apply to columns + for (int i = 0; i < n_dim; ++i) + { + Real aip = dense[i + p * n_dim]; + Real aiq = dense[i + q * n_dim]; + dense[i + p * n_dim] = c * aip + s * aiq; + dense[i + q * n_dim] = -s * aip + c * aiq; + } + // Apply to rows + for (int j = 0; j < n_dim; ++j) + { + Real apj = dense[p + j * n_dim]; + Real aqj = dense[q + j * n_dim]; + dense[p + j * n_dim] = c * apj + s * aqj; + dense[q + j * n_dim] = -s * apj + c * aqj; + } + }; + + // Several rotations with different angles to create a genuinely + // dense matrix (all off-diagonals become non-zero) + std::mt19937 rng_dense(111); + std::uniform_real_distribution angle_dist(Real(0.2), Real(1.3)); + for (int k = 0; k < 20; ++k) + { + int p = k % (n_dim - 1); + int q = p + 1 + (k / (n_dim - 1)) % (n_dim - 1 - p); + if (q >= n_dim) + { + q = n_dim - 1; + } + if (p == q) + { + continue; + } + const Real angle = angle_dist(rng_dense); + apply_givens(p, q, angle); + } + + // Copy to complex H_mat + H_mat.assign(n_dim * n_dim, T(0)); + for (int i = 0; i < n_dim * n_dim; ++i) + { + H_mat[i] = T(dense[i], 0); + } + + // Preconditioner: use diagonal of the rotated H + prec.resize(n_dim); + for (int i = 0; i < n_dim; ++i) + { + prec[i] = std::real(H_mat[i + i * n_dim]); + } + + // Lowest 4 eigenvalues: 1, 2, 3, 4 + exact = {1.0, 2.0, 3.0, 4.0}; + check_eigenvalues(1e-12, 200, 4, 4, 1e-8); +} + +// ============================================================================= +// Euclidean column norm for independent residual checks +// ============================================================================= +static Real column_norm(const T* x, int n_dim) +{ + Real nrm = 0; + for (int i = 0; i < n_dim; ++i) + { + nrm += std::norm(x[i]); + } + return std::sqrt(nrm); +} + +// ============================================================================= +// Test fixture: non-trivial S matrix (diagonal overlap, S ≠ I) +// H = tridiag 6×6 Laplacian, S = diag(1.1, 1.0, 0.9, 1.0, 1.1, 1.0) +// Tests that the solver correctly handles a non-identity overlap matrix. +// ============================================================================= +class DiagoPPCGWithSTest : public ::testing::Test +{ + protected: + void SetUp() override + { + n_dim = 6; + nband = 3; + ld = n_dim + 2; // exercise custom S with padded leading dimension + + // Tridiagonal H + H_mat.assign(n_dim * n_dim, T(0)); + for (int i = 0; i < n_dim; ++i) + { + H_mat[i + i * n_dim] = T(2.0, 0); + if (i > 0) + { + H_mat[i + (i - 1) * n_dim] = T(-1.0, 0); + } + if (i < n_dim - 1) + { + H_mat[i + (i + 1) * n_dim] = T(-1.0, 0); + } + } + + // S = diag(1.1, 1.0, 0.9, 1.0, 1.1, 1.0) + s_diag = {1.1, 1.0, 0.9, 1.0, 1.1, 1.0}; + S_mat.assign(n_dim * n_dim, T(0)); + for (int i = 0; i < n_dim; ++i) + { + S_mat[i + i * n_dim] = T(s_diag[i], 0); + } + + prec.assign(n_dim, 2.0); + + ethr.assign(nband, 1e-8); + + std::mt19937 rng(333); + std::uniform_real_distribution dist(-1.0, 1.0); + + psi.assign(ld * nband, T(0)); + for (int j = 0; j < nband; ++j) + { + for (int i = 0; i < n_dim; ++i) + { + psi[i + j * ld] = T(dist(rng), 0.0); + } + } + + // S-orthonormalize initial guess + for (int j = 0; j < nband; ++j) + { + for (int k = 0; k < j; ++k) + { + T dot = 0; + for (int i = 0; i < n_dim; ++i) + { + dot += std::conj(psi[i + k * ld]) * T(s_diag[i], 0) * psi[i + j * ld]; + } + for (int i = 0; i < n_dim; ++i) + { + psi[i + j * ld] -= dot * psi[i + k * ld]; + } + } + Real nrm = 0; + for (int i = 0; i < n_dim; ++i) + { + nrm += s_diag[i] * std::norm(psi[i + j * ld]); + } + nrm = std::sqrt(nrm); + for (int i = 0; i < n_dim; ++i) + { + psi[i + j * ld] /= nrm; + } + } + } + + int n_dim; + int nband; + int ld; + std::vector H_mat; + std::vector S_mat; + std::vector s_diag; + std::vector prec; + std::vector ethr; + std::vector psi; +}; + +TEST_F(DiagoPPCGWithSTest, BlockSubspace) +{ + std::vector psi_run = psi; + std::vector eval(nband, 0.0); + + // S-apply function: S * psi = diag(s_diag) * psi (element-wise) + auto spsi_func = [this](T* in, T* out, int ld_in, int ncol) { + for (int j = 0; j < ncol; ++j) + { + for (int i = 0; i < n_dim; ++i) + { + out[i + j * ld_in] = T(s_diag[i], 0) * in[i + j * ld_in]; + } + } + }; + + hsolver::DiagoPPCG solver( + /* diag_thr = */ 1e-10, + /* max_iter = */ 100, + /* sbsize = */ 3, + /* rr_step = */ 3); + + auto h_op = [this](T* in, T* out, int ld_in, int ncol) { + dense_h_multiply(H_mat.data(), n_dim, in, out, ld_in, ncol); + }; + + double avg_iter = solver.diag(h_op, spsi_func, ld, nband, n_dim, psi_run.data(), eval.data(), ethr, prec.data()); + + // Eigenvalue check: skip absolute comparison (exact values not + // analytically known for non-trivial S). Instead verify via residual. + // Just check eigenvalues are reasonable (not NaN, not negative for + // this positive-definite problem). + for (int i = 0; i < nband; ++i) + { + EXPECT_GT(eval[i], 0.0) << "WithS BLOCK: eigenvalue[" << i << "] should be positive"; + EXPECT_LT(eval[i], 10.0) << "WithS BLOCK: eigenvalue[" << i << "] unreasonably large"; + } + + // Independently check the absolute residual ||Hψ_i - ε_i S ψ_i||. + std::vector hpsi(n_dim); + std::vector spsi(n_dim); + std::vector res(n_dim); + for (int i = 0; i < nband; ++i) + { + T* psi_column = psi_run.data() + i * ld; + dense_h_multiply(H_mat.data(), n_dim, psi_column, hpsi.data(), n_dim, 1); + spsi_func(psi_column, spsi.data(), n_dim, 1); + for (int j = 0; j < n_dim; ++j) + { + res[j] = hpsi[j] - T(eval[i], 0) * spsi[j]; + } + Real res_nrm = column_norm(res.data(), n_dim); + const double scaled_threshold = 1e2 * ethr[i]; + const double residual_threshold = std::max(1e-6, scaled_threshold); + EXPECT_LE(res_nrm, residual_threshold) + << "WithS BLOCK: residual[" << i << "] too large, r=" << res_nrm; + } + + EXPECT_LE(avg_iter, double(100)) << "WithS BLOCK: too many iterations"; +} + +// ============================================================================= +// Test fixture: single-band (nband = 1) +// Minimal test — extract only the lowest eigenvalue of a 5×5 tridiagonal +// Laplacian. This exercises the degenerate code paths for a single band. +// ============================================================================= +TEST_F(DiagoPPCGMatrixTest, SingleBand) +{ + initialize(5, 1, 42, 1e-10); + set_laplacian(1.0); + check_eigenvalues(1e-12, 50, 1, 1, 1e-8); +} + +// ============================================================================= +// Test fixture: eigenvector quality — verify Hψ ≈ εψ and ψ^H ψ = I +// Uses the 10×10 tridiagonal Laplacian. After convergence, check: +// 1. ||Hψ_i - ε_i ψ_i|| < tol for each band +// 2. |ψ_i^H ψ_j - δ_ij| < tol for all i,j +// ============================================================================= +class DiagoPPCGEigenvectorTest : public DiagoPPCGMatrixTest +{ + protected: + void SetUp() override + { + initialize(10, 3, 888, 1e-08); + set_laplacian(1.0); + } +}; + +TEST_F(DiagoPPCGEigenvectorTest, BlockSubspace) +{ + std::vector psi_run = psi; + std::vector eval(nband, 0.0); + + hsolver::DiagoPPCG solver( + /* diag_thr = */ 1e-12, + /* max_iter = */ 100, + /* sbsize = */ 3, + /* rr_step = */ 3); + + auto h_op = [this](T* in, T* out, int ld_in, int ncol) { + dense_h_multiply(H_mat.data(), n_dim, in, out, ld_in, ncol); + }; + + double avg_iter = solver.diag(h_op, nullptr, ld, nband, n_dim, psi_run.data(), eval.data(), ethr, prec.data()); + + // --- Eigenvalue check --- + for (int i = 0; i < nband; ++i) + { + EXPECT_NEAR(eval[i], exact[i], 1e-8) << "Eigenvec BLOCK: eigenvalue[" << i << "] mismatch"; + } + + // Independently check the returned eigenvectors using a fresh H application. + std::vector hpsi(n_dim); + std::vector res(n_dim); + for (int i = 0; i < nband; ++i) + { + const T* psi_column = psi_run.data() + i * ld; + dense_h_multiply(H_mat.data(), n_dim, psi_column, hpsi.data(), n_dim, 1); + for (int j = 0; j < n_dim; ++j) + { + res[j] = hpsi[j] - eval[i] * psi_run[j + i * ld]; + } + Real res_nrm = column_norm(res.data(), n_dim); + EXPECT_LT(res_nrm, 1e-4) << "Eigenvec BLOCK: residual[" << i << "] too large: " << res_nrm; + } + + // --- Orthogonality check: |ψ_i^H ψ_j - δ_ij| < 1e-8 --- + for (int i = 0; i < nband; ++i) + { + for (int j = 0; j < nband; ++j) + { + T dot = 0; + for (int k = 0; k < n_dim; ++k) + { + dot += std::conj(psi_run[k + i * ld]) * psi_run[k + j * ld]; + } + if (i == j) + { + EXPECT_NEAR(std::abs(dot), 1.0, 1e-8) + << "Eigenvec BLOCK: ψ[" << i << "] not normalized, |dot|=" << std::abs(dot); + } + else + { + EXPECT_LT(std::abs(dot), 1e-8) + << "Eigenvec BLOCK: ψ[" << i << "] not orthogonal to ψ[" << j << "], |dot|=" << std::abs(dot); + } + } + } + + EXPECT_LE(avg_iter, double(100)) << "Eigenvec BLOCK: too many iterations"; +} + +// ============================================================================= +// Test fixture: all eigenvalues of a small matrix (nband = n_dim) +// 3×3 tridiagonal Laplacian, compute all 3 eigenvalues. +// Exercises the degenerate case where every band is requested. +// ============================================================================= +TEST_F(DiagoPPCGMatrixTest, AllBands) +{ + initialize(3, 3, 101, 1e-10); + set_laplacian(1.0); + check_eigenvalues(1e-12, 100, 3, 3, 1e-8); +} + +// ============================================================================= +// Test fixture: medium-sized tridiagonal (15×15, nband=4) +// Bridges the gap between the 10×10 and 20×20 tests. +// ============================================================================= +TEST_F(DiagoPPCGMatrixTest, Laplacian15) +{ + initialize(15, 4, 202, 1e-10); + set_laplacian(1.0); + check_eigenvalues(1e-12, 120, 4, 4, 1e-8); +} + +// ============================================================================= +// Test fixture: pentadiagonal Toeplitz (discrete biharmonic operator) +// H[i,i]=6, H[i,i±1]=-4, H[i,i±2]=1. Eigenvalues: +// λ_k = 16·sin⁴(k·π / (2·(n+1))), k = 1,...,n +// Wider bandwidth (5 vs 3) tests the solver with more off-diagonal coupling. +// ============================================================================= +TEST_F(DiagoPPCGMatrixTest, SquaredLaplacian) +{ + initialize(8, 4, 505, 1e-10); + // H = T² where T is the tridiagonal Laplacian (2 on diag, -1 on off-diag). + // The corners of T² have diag=5 (not 6) since (T²)[0,0] = 2² + (-1)² = 5. + // Interior: (T²)[i,i] = (-1)² + 2² + (-1)² = 6. + H_mat.assign(n_dim * n_dim, T(0)); + for (int i = 0; i < n_dim; ++i) + { + const Real diagonal_value = (i == 0 || i == n_dim - 1) ? 5.0 : 6.0; + H_mat[i + i * n_dim] = T(diagonal_value, 0); + if (i >= 1) + { + H_mat[i + (i - 1) * n_dim] = T(-4.0, 0); + } + if (i < n_dim - 1) + { + H_mat[i + (i + 1) * n_dim] = T(-4.0, 0); + } + if (i >= 2) + { + H_mat[i + (i - 2) * n_dim] = T(1.0, 0); + } + if (i < n_dim - 2) + { + H_mat[i + (i + 2) * n_dim] = T(1.0, 0); + } + } + + prec.assign(n_dim, 6.0); + prec[0] = 5.0; + prec[n_dim - 1] = 5.0; + + exact.resize(nband); + for (int k = 0; k < nband; ++k) + { + Real theta = Real(k + 1) * PPCG_PI / Real(2 * (n_dim + 1)); + Real s = std::sin(theta); + exact[k] = Real(16) * s * s * s * s; + } + check_eigenvalues(1e-12, 150, 4, 4, 1e-8); +} + +// ============================================================================= +// Test fixture: gapped spectrum +// H = diag(0.1, 0.5, 5.0, 6.0, 10.0), nband=3. +// Large gaps between eigenvalue groups test the solver's band separation. +// ============================================================================= +TEST_F(DiagoPPCGMatrixTest, GappedSpectrum) +{ + initialize(5, 3, 606, 1e-10); + H_mat.assign(n_dim * n_dim, T(0)); + H_mat[0 + 0 * n_dim] = T(0.1, 0); + H_mat[1 + 1 * n_dim] = T(0.5, 0); + H_mat[2 + 2 * n_dim] = T(5.0, 0); + H_mat[3 + 3 * n_dim] = T(6.0, 0); + H_mat[4 + 4 * n_dim] = T(10.0, 0); + + prec.assign(n_dim, 1.0); + + exact = {0.1, 0.5, 5.0}; + check_eigenvalues(1e-12, 100, 3, 3, 1e-8); +} + +// ============================================================================= +// Test fixture: preconditioner stress test +// Uses the 10×10 tridiagonal Laplacian with a fixed, nonuniform positive +// preconditioner that differs from the exact diagonal. The test checks that +// the solver still converges and returns a small independently computed residual. +// ============================================================================= +class DiagoPPCGBadPrecTest : public DiagoPPCGMatrixTest +{ + protected: + void SetUp() override + { + initialize(10, 3, 707, 1e-10); + set_laplacian(1.0); + for (int i = 0; i < n_dim; ++i) + { + prec[i] = 0.6 + 0.1 * Real(i); + } + } +}; + +TEST_F(DiagoPPCGBadPrecTest, BlockSubspace) +{ + std::vector psi_run = psi; + std::vector eval(nband, 0.0); + + hsolver::DiagoPPCG solver( + /* diag_thr = */ 1e-12, + /* max_iter = */ 200, + /* sbsize = */ 4, + /* rr_step = */ 4); + + auto h_op = [this](T* in, T* out, int ld_in, int ncol) { + dense_h_multiply(H_mat.data(), n_dim, in, out, ld_in, ncol); + }; + + double avg_iter = solver.diag(h_op, nullptr, ld, nband, n_dim, psi_run.data(), eval.data(), ethr, prec.data()); + + for (int i = 0; i < nband; ++i) + { + EXPECT_NEAR(eval[i], exact[i], 1e-8) << "BadPrec BLOCK: eigenvalue[" << i << "] mismatch"; + } + EXPECT_TRUE(solver.converged()) << "BadPrec BLOCK: solver did not converge"; + EXPECT_LE(avg_iter, double(200)) << "BadPrec BLOCK: too many iterations"; + + for (int j = 0; j < nband; ++j) + { + std::vector hpsi(n_dim, T(0)); + const T* psi_column = psi_run.data() + j * ld; + dense_h_multiply(H_mat.data(), n_dim, psi_column, hpsi.data(), n_dim, 1); + Real residual_squared = 0; + for (int i = 0; i < n_dim; ++i) + { + const T residual = hpsi[i] - T(eval[j]) * psi_run[i + j * ld]; + residual_squared += std::norm(residual); + } + const Real residual = std::sqrt(residual_squared); + EXPECT_TRUE(std::isfinite(residual)); + EXPECT_LT(residual, 1e-8) << "BadPrec BLOCK: residual[" << j << "] too large"; + } +} + +// ============================================================================= +// Test fixture: n_dim = 1, nband = 1 — absolute minimum +// H is a 1×1 matrix [5.0], eigenvalue = 5.0 +// ============================================================================= +TEST_F(DiagoPPCGMatrixTest, ScalarMatrix) +{ + initialize(1, 1, 0, 1e-10); + H_mat = {T(5.0, 0)}; + prec = {5.0}; + exact = {5.0}; + psi = {T(1.0, 0)}; // already normalized + check_eigenvalues(1e-12, 10, 1, 1, 1e-8); +} + +// ============================================================================= +// Test fixture: scaled tridiagonal (eigenvalues × 100) +// H = 100 × tridiag(2, -1, -1). Tests convergence with large eigenvalues. +// ============================================================================= +TEST_F(DiagoPPCGMatrixTest, ScaledLaplacian) +{ + initialize(8, 3, 808, 1e-08); + set_laplacian(100.0); + check_eigenvalues(1e-10, 120, 4, 4, 1e-6); +} + +// ============================================================================= +// Test fixture: many bands (n_dim=12, nband=4) +// Moderate band-to-dimension ratio (1:3). +// ============================================================================= +TEST_F(DiagoPPCGMatrixTest, ManyBands) +{ + initialize(12, 4, 909, 1e-10); + set_laplacian(1.0); + check_eigenvalues(1e-12, 150, 4, 4, 1e-8); +} + +// ============================================================================= +// Test fixture: rr_step = 1 (most frequent Rayleigh-Ritz refinement) +// 8×8 tridiagonal, nband=3. Tests aggressive subspace diagonalization. +// ============================================================================= +TEST_F(DiagoPPCGMatrixTest, RitzEveryIteration) +{ + initialize(8, 3, 111, 1e-10); + set_laplacian(1.0); + check_eigenvalues(1e-12, 100, 3, 1, 1e-8); +} + +// ============================================================================= +// Test fixture: Neumann boundary tridiagonal Laplacian +// H[0,0]=H[n-1,n-1]=1, interior diag=2, off-diag=-1. +// Eigenvalues: λ_k = 2 - 2·cos(k·π/n), k=0,1,...,n-1 +// ============================================================================= +TEST_F(DiagoPPCGMatrixTest, NeumannZeroEigenvalue) +{ + initialize(8, 4, 222, 1e-10); + H_mat.assign(n_dim * n_dim, T(0)); + for (int i = 0; i < n_dim; ++i) + { + const Real diagonal_value = (i == 0 || i == n_dim - 1) ? 1.0 : 2.0; + H_mat[i + i * n_dim] = T(diagonal_value, 0); + if (i > 0) + { + H_mat[i + (i - 1) * n_dim] = T(-1.0, 0); + } + if (i < n_dim - 1) + { + H_mat[i + (i + 1) * n_dim] = T(-1.0, 0); + } + } + prec.assign(n_dim, 2.0); + prec[0] = 1.0; + prec[n_dim - 1] = 1.0; + exact.resize(nband); + for (int k = 0; k < nband; ++k) + { + exact[k] = 2.0 - 2.0 * std::cos(Real(k) * PPCG_PI / Real(n_dim)); + } + check_eigenvalues(1e-12, 100, 4, 4, 1e-8); +} + +// ============================================================================= +// Test fixture: tight convergence threshold (ethr = 1e-14) +// 6×6 tridiagonal, nband=2. Tests deep convergence. +// ============================================================================= +TEST_F(DiagoPPCGMatrixTest, TightThreshold) +{ + initialize(6, 2, 333, 1e-14); + set_laplacian(1.0); + check_eigenvalues(1e-14, 200, 3, 3, 1e-8); +} + +// ============================================================================= +// Test fixture: non-diagonal S matrix (tridiagonal overlap) +// H = tridiag(2,-1,-1), S = tridiag(1.0, 0.2, 0.2) — both tridiagonal. +// Verifies S-orthogonalization with a non-diagonal S. +// ============================================================================= +class DiagoPPCGTridiagSTest : public ::testing::Test +{ + protected: + void SetUp() override + { + n_dim = 6; + nband = 2; + ld = n_dim; + H_mat.assign(n_dim * n_dim, T(0)); + for (int i = 0; i < n_dim; ++i) + { + H_mat[i + i * n_dim] = T(2.0, 0); + if (i > 0) + { + H_mat[i + (i - 1) * n_dim] = T(-1.0, 0); + } + if (i < n_dim - 1) + { + H_mat[i + (i + 1) * n_dim] = T(-1.0, 0); + } + } + S_mat.assign(n_dim * n_dim, T(0)); + for (int i = 0; i < n_dim; ++i) + { + S_mat[i + i * n_dim] = T(1.0, 0); + if (i > 0) + { + S_mat[i + (i - 1) * n_dim] = T(0.2, 0); + S_mat[(i - 1) + i * n_dim] = T(0.2, 0); + } + } + prec.assign(n_dim, 2.0); + // H and S share sine eigenvectors. Their generalized eigenvalues are + // (2 - 2 cos(theta)) / (1 + 0.4 cos(theta)), theta = k*pi/(n_dim+1). + exact.resize(nband); + for (int k = 0; k < nband; ++k) + { + const Real theta = Real(k + 1) * PPCG_PI / Real(n_dim + 1); + const Real cosine = std::cos(theta); + exact[k] = (2.0 - 2.0 * cosine) / (1.0 + 0.4 * cosine); + } + ethr.assign(nband, 1e-8); + init_psi(444); + } + void init_psi(int seed) + { + std::mt19937 rng(seed); + std::uniform_real_distribution dist(-1.0, 1.0); + psi.assign(ld * nband, T(0)); + for (int j = 0; j < nband; ++j) + { + for (int i = 0; i < n_dim; ++i) + { + psi[i + j * ld] = T(dist(rng), 0.0); + } + } + // S-orthonormalize: S = tridiag(1.0, 0.2, 0.2) + for (int j = 0; j < nband; ++j) + { + for (int k = 0; k < j; ++k) + { + T d = 0; + for (int i = 0; i < n_dim; ++i) + { + T si = 0; + si += T(1.0, 0) * psi[i + k * ld]; + if (i > 0) + { + si += T(0.2, 0) * psi[(i - 1) + k * ld]; + } + if (i < n_dim - 1) + { + si += T(0.2, 0) * psi[(i + 1) + k * ld]; + } + d += std::conj(si) * psi[i + j * ld]; + } + for (int i = 0; i < n_dim; ++i) + { + psi[i + j * ld] -= d * psi[i + k * ld]; + } + } + Real nr = 0; + for (int i = 0; i < n_dim; ++i) + { + T si = 0; + si += T(1.0, 0) * psi[i + j * ld]; + if (i > 0) + { + si += T(0.2, 0) * psi[(i - 1) + j * ld]; + } + if (i < n_dim - 1) + { + si += T(0.2, 0) * psi[(i + 1) + j * ld]; + } + nr += std::real(std::conj(psi[i + j * ld]) * si); + } + nr = std::sqrt(nr); + for (int i = 0; i < n_dim; ++i) + { + psi[i + j * ld] /= nr; + } + } + } + int n_dim; + int nband; + int ld; + std::vector H_mat; + std::vector S_mat; + std::vector prec; + std::vector exact; + std::vector ethr; + std::vector psi; +}; + +TEST_F(DiagoPPCGTridiagSTest, BlockSubspace) +{ + std::vector psi_run = psi; + std::vector eval(nband, 0.0); + auto spsi_func = [this](T* in, T* out, int ldi, int nc) { + for (int j = 0; j < nc; ++j) + { + for (int i = 0; i < n_dim; ++i) + { + out[i + j * ldi] = T(1.0, 0) * in[i + j * ldi]; + if (i > 0) + { + out[i + j * ldi] += T(0.2, 0) * in[(i - 1) + j * ldi]; + } + if (i < n_dim - 1) + { + out[i + j * ldi] += T(0.2, 0) * in[(i + 1) + j * ldi]; + } + } + } + }; + hsolver::DiagoPPCG solver(1e-10, 150, 3, 3); + auto h_op = [this](T* in, T* out, int ldi, int nc) { dense_h_multiply(H_mat.data(), n_dim, in, out, ldi, nc); }; + const double avg_iter = solver.diag(h_op, spsi_func, ld, nband, n_dim, psi_run.data(), eval.data(), ethr, prec.data()); + + EXPECT_TRUE(solver.converged()); + EXPECT_EQ(solver.active_band_count(), 0); + EXPECT_FALSE(solver.iteration_limit_reached()); + for (int i = 0; i < nband; ++i) + { + EXPECT_NEAR(eval[i], exact[i], ethr[i]) << "TridiagS BLOCK: wrong eigenvalue at band " << i; + } + + // Reapply the dense matrices independently of the solver's operator images. + const int image_size = ld * nband; + std::vector hpsi(image_size); + std::vector spsi(image_size); + std::vector res(n_dim); + dense_h_multiply(H_mat.data(), n_dim, psi_run.data(), hpsi.data(), ld, nband); + dense_h_multiply(S_mat.data(), n_dim, psi_run.data(), spsi.data(), ld, nband); + for (int i = 0; i < nband; ++i) + { + for (int j = 0; j < n_dim; ++j) + { + const int offset = j + i * ld; + res[j] = hpsi[offset] - eval[i] * spsi[offset]; + } + const Real residual_norm = column_norm(res.data(), n_dim); + EXPECT_LE(residual_norm, ethr[i]) << "TridiagS BLOCK: residual at band " << i; + + for (int k = 0; k < nband; ++k) + { + T overlap = 0; + for (int j = 0; j < n_dim; ++j) + { + overlap += std::conj(psi_run[j + k * ld]) * spsi[j + i * ld]; + } + const T expected_overlap = (i == k) ? T(1) : T(0); + const T overlap_difference = overlap - expected_overlap; + const Real overlap_error = std::abs(overlap_difference); + EXPECT_LE(overlap_error, 1e-10) << "TridiagS BLOCK: S-overlap at (" << k << ", " << i << ")"; + } + } + EXPECT_LE(avg_iter, 150.0) << "TridiagS BLOCK: too many iterations"; +} + +// A deterministic sparse matrix exercises convergence beyond analytical fixtures. +TEST(DiagoPPCGRandomMatrixTest, ConvergesWithDiagonalPreconditioner) +{ + const int n_dim = 80; + const int nband = 8; + const int sparsity_pct = 60; + const int matrix_size = n_dim * n_dim; + const unsigned seed = n_dim * 100 + sparsity_pct; + std::vector h(matrix_size, T(0)); + std::mt19937 rng(seed); + std::uniform_real_distribution dist(-1.0, 1.0); + for (int i = 0; i < n_dim; ++i) + { + for (int j = i; j < n_dim; ++j) + { + if (i != j && (rng() % 100) < sparsity_pct) + { + continue; + } + const Real value = (i == j) ? std::abs(dist(rng)) * n_dim + 1.0 : dist(rng) * 0.5; + h[i + j * n_dim] = T(value, 0); + h[j + i * n_dim] = T(value, 0); + } + } + std::vector prec(n_dim); + for (int i = 0; i < n_dim; ++i) + { + const Real diagonal = std::real(h[i + i * n_dim]); + prec[i] = std::max(diagonal, 1e-6); + } + std::vector psi; + initialize_vectors(psi, n_dim, nband, n_dim, 42); + std::vector eval(nband, 0.0); + const std::vector ethr(nband, 1e-4); + const auto h_op = [&h, n_dim](T* in, T* out, int ld, int ncol) { + dense_h_multiply(h.data(), n_dim, in, out, ld, ncol); + }; + hsolver::DiagoPPCG solver(1e-8, 500, nband, 4); + const double iterations = solver.diag(h_op, nullptr, n_dim, nband, n_dim, + psi.data(), eval.data(), ethr, prec.data()); + EXPECT_TRUE(solver.converged()); + EXPECT_TRUE(std::isfinite(solver.max_residual())); + EXPECT_LE(solver.max_residual(), 1e-4); + EXPECT_LE(iterations, 500.0); +} + +int main(int argc, char** argv) +{ +#ifdef __MPI + int nproc = 1; + int myrank = 0; + setupmpi(argc, argv, nproc, myrank); + int nproc_in_pool = 1; + int kpar = nproc; + int mypool = 0; + int rank_in_pool = 0; + divide_pools(nproc, myrank, nproc_in_pool, kpar, mypool, rank_in_pool); +#endif + testing::InitGoogleTest(&argc, argv); + const int result = RUN_ALL_TESTS(); +#ifdef __MPI + finishmpi(); +#endif + return result; +} diff --git a/source/source_hsolver/unittests/test_diago_ppcg_float.cpp b/source/source_hsolver/unittests/test_diago_ppcg_float.cpp new file mode 100644 index 00000000000..4bf6c61f8ce --- /dev/null +++ b/source/source_hsolver/unittests/test_diago_ppcg_float.cpp @@ -0,0 +1,304 @@ +/** + * diago_ppcg_float_test.cpp — single-precision unit test for DiagoPPCG. + * + * Exercises the std::complex instantiation of the BLOCK_SUBSPACE and + * BLOCK_SUBSPACE strategy on dense matrices with analytical reference + * eigenvalues. Tolerances are looser than the double-precision suite because + * single precision has roughly 7 significant digits. + */ + +#include "../diago_ppcg.h" +#ifdef __MPI +#include "source_basis/module_pw/test/test_tool.h" +#endif + +#include +#include +#include +#include +#include +#include + +static constexpr float PPCG_PI = 3.14159265358979323846F; + +using T = std::complex; +using Real = float; + +// ----------------------------------------------------------------------------- +// Helper: dense H-matrix times a set of column vectors (column-major H). +// ----------------------------------------------------------------------------- +static void dense_h_multiply(const T* H_data, int n_dim, + const T* in, T* out, int ld, int ncol) +{ + for (int j = 0; j < ncol; ++j) + { + for (int i = 0; i < n_dim; ++i) + { + T sum = T(0.0f, 0.0f); + for (int k = 0; k < n_dim; ++k) + { + sum += H_data[i + k * n_dim] * in[k + j * ld]; + } + out[i + j * ld] = sum; + } + } +} + +// Orthonormalize columns of psi in-place (S = I). +static void gram_schmidt(std::vector& psi, int ld, int n_dim, int nband) +{ + for (int j = 0; j < nband; ++j) + { + for (int k = 0; k < j; ++k) + { + T dot = T(0.0f, 0.0f); + for (int i = 0; i < n_dim; ++i) + { + dot += std::conj(psi[i + k * ld]) * psi[i + j * ld]; + } + for (int i = 0; i < n_dim; ++i) + { + psi[i + j * ld] -= dot * psi[i + k * ld]; + } + } + Real nrm = 0.0f; + for (int i = 0; i < n_dim; ++i) + { + nrm += std::norm(psi[i + j * ld]); + } + nrm = std::sqrt(nrm); + for (int i = 0; i < n_dim; ++i) + { + psi[i + j * ld] /= nrm; + } + } +} + +// ----------------------------------------------------------------------------- +// Diagonal matrix: H = diag(1, 2, 3). All eigenvalues are computed (nband == +// n_dim), so there is no ambiguity about which end of the spectrum to converge +// to; single-precision Rayleigh-Ritz can otherwise drift toward the upper +// eigenvalues on some platforms. +// ----------------------------------------------------------------------------- +TEST(DiagoPPCGFloatTest, DiagonalBlockSubspace) +{ + const int n_dim = 3; + const int nband = 3; + const int ld = n_dim; + + std::vector H_mat(n_dim * n_dim, T(0.0f, 0.0f)); + for (int i = 0; i < n_dim; ++i) + { + H_mat[i + i * n_dim] = T(Real(i + 1), 0.0f); + } + + std::vector prec(n_dim); + for (int i = 0; i < n_dim; ++i) + { + prec[i] = Real(i + 1); + } + + const Real exact[3] = {1.0f, 2.0f, 3.0f}; + std::vector ethr(nband, 1e-4); + + std::mt19937 rng(42); + std::uniform_real_distribution dist(-1.0f, 1.0f); + std::vector psi(ld * nband, T(0.0f, 0.0f)); + for (int j = 0; j < nband; ++j) + { + for (int i = 0; i < n_dim; ++i) + { + psi[i + j * ld] = T(dist(rng), 0.0f); + } + } + gram_schmidt(psi, ld, n_dim, nband); + + std::vector psi_run = psi; + std::vector eval(nband, 0.0f); + + hsolver::DiagoPPCG solver( + /* diag_thr = */ 1e-5f, + /* max_iter = */ 100, + /* sbsize = */ 3, + /* rr_step = */ 3); + + auto h_op = [&](T* in, T* out, int ld_in, int ncol) { + dense_h_multiply(H_mat.data(), n_dim, in, out, ld_in, ncol); + }; + + double avg_iter = solver.diag(h_op, nullptr, ld, nband, n_dim, + psi_run.data(), eval.data(), ethr, prec.data()); + + for (int i = 0; i < nband; ++i) + { + EXPECT_NEAR(double(eval[i]), double(exact[i]), 1e-4) + << "Diagonal float BLOCK: eigenvalue[" << i << "] mismatch"; + } + EXPECT_LE(avg_iter, 100.0) << "Diagonal float BLOCK: too many iterations"; +} + +// ----------------------------------------------------------------------------- +// Tridiagonal Laplacian: H[i,i]=2, H[i,i±1]=-1, exact λ_k = 2 - 2cos(kπ/(n+1)) +// ----------------------------------------------------------------------------- +TEST(DiagoPPCGFloatTest, TridiagonalBlockSubspace) +{ + const int n_dim = 10; + const int nband = 3; + const int ld = n_dim; + + std::vector H_mat(n_dim * n_dim, T(0.0f, 0.0f)); + for (int i = 0; i < n_dim; ++i) + { + H_mat[i + i * n_dim] = T(2.0f, 0.0f); + if (i > 0) + { + H_mat[i + (i - 1) * n_dim] = T(-1.0f, 0.0f); + } + if (i < n_dim - 1) + { + H_mat[i + (i + 1) * n_dim] = T(-1.0f, 0.0f); + } + } + + std::vector prec(n_dim, 2.0f); + std::vector exact(nband); + for (int k = 0; k < nband; ++k) + { + exact[k] = 2.0f - 2.0f * std::cos(Real(k + 1) * PPCG_PI + / Real(n_dim + 1)); + } + std::vector ethr(nband, 1e-4); + + std::mt19937 rng(42); + std::uniform_real_distribution dist(-1.0f, 1.0f); + std::vector psi(ld * nband, T(0.0f, 0.0f)); + for (int j = 0; j < nband; ++j) + { + for (int i = 0; i < n_dim; ++i) + { + psi[i + j * ld] = T(dist(rng), 0.0f); + } + } + gram_schmidt(psi, ld, n_dim, nband); + + std::vector psi_run = psi; + std::vector eval(nband, 0.0f); + + hsolver::DiagoPPCG solver( + /* diag_thr = */ 1e-5f, + /* max_iter = */ 100, + /* sbsize = */ 4, + /* rr_step = */ 4); + + auto h_op = [&](T* in, T* out, int ld_in, int ncol) { + dense_h_multiply(H_mat.data(), n_dim, in, out, ld_in, ncol); + }; + + double avg_iter = solver.diag(h_op, nullptr, ld, nband, n_dim, + psi_run.data(), eval.data(), ethr, prec.data()); + + for (int i = 0; i < nband; ++i) + { + EXPECT_NEAR(double(eval[i]), double(exact[i]), 1e-4) + << "Tridiagonal float BLOCK: eigenvalue[" << i << "] mismatch"; + } + EXPECT_LE(avg_iter, 100.0) << "Tridiagonal float BLOCK: too many iterations"; +} + +// ----------------------------------------------------------------------------- +// Non-finite input validation (throws). +// ----------------------------------------------------------------------------- +TEST(DiagoPPCGFloatTest, NonFiniteInputThrows) +{ + const int n_dim = 5; + const int nband = 3; + const int ld = n_dim; + + std::vector H_mat(n_dim * n_dim, T(0.0f, 0.0f)); + for (int i = 0; i < n_dim; ++i) + { + H_mat[i + i * n_dim] = T(Real(i + 1), 0.0f); + } + + std::vector prec(n_dim, 1.0f); + std::vector psi(ld * nband, T(1.0f, 0.0f)); + std::vector eval(nband, 0.0f); + std::vector ethr(nband, 1e-4); + + hsolver::DiagoPPCG solver( + /* diag_thr = */ 1e-5f, 100, 3, 3); + + auto h_op = [&](T* in, T* out, int ld_in, int ncol) { + dense_h_multiply(H_mat.data(), n_dim, in, out, ld_in, ncol); + }; + + std::vector bad_ethr = ethr; + bad_ethr[0] = std::numeric_limits::quiet_NaN(); + EXPECT_THROW(solver.diag(h_op, nullptr, ld, nband, n_dim, + psi.data(), eval.data(), bad_ethr, prec.data()), + std::invalid_argument); + + std::vector bad_prec = prec; + bad_prec[0] = std::numeric_limits::quiet_NaN(); + EXPECT_THROW(solver.diag(h_op, nullptr, ld, nband, n_dim, + psi.data(), eval.data(), ethr, bad_prec.data()), + std::invalid_argument); +} + +TEST(DiagoPPCGFloatTest, RankDeficientStartupSubspace) +{ + const int n_dim = 4; + const int nband = 2; + const int ld = n_dim; + const Real inv_sqrt_two = Real(1) / std::sqrt(Real(2)); + std::vector h_mat(n_dim * n_dim, T(0)); + std::vector prec(n_dim); + for (int i = 0; i < n_dim; ++i) + { + const Real value = Real(i + 1); + h_mat[i + i * n_dim] = T(value, 0); + prec[i] = value; + } + + std::vector psi(ld * nband, T(0)); + psi[0] = T(inv_sqrt_two, 0); + psi[1] = T(inv_sqrt_two, 0); + psi[3 + ld] = T(1, 0); + std::vector eval(nband, Real(0)); + const std::vector ethr(nband, 1e-4); + const auto h_op = [&h_mat, n_dim](T* in, T* out, int ld_in, int ncol) { + dense_h_multiply(h_mat.data(), n_dim, in, out, ld_in, ncol); + }; + hsolver::DiagoPPCG solver( + /* diag_thr = */ 1e-5f, + /* max_iter = */ 20, + /* sbsize = */ 1, + /* rr_step = */ 2); + + solver.diag(h_op, nullptr, ld, nband, n_dim, + psi.data(), eval.data(), ethr, prec.data()); + + EXPECT_NEAR(eval[0], Real(1), Real(2e-4)); + EXPECT_NEAR(eval[1], Real(2), Real(2e-4)); + EXPECT_TRUE(solver.converged()); +} + +int main(int argc, char** argv) +{ +#ifdef __MPI + int nproc = 1; + int myrank = 0; + setupmpi(argc, argv, nproc, myrank); + int nproc_in_pool = 1; + int kpar = nproc; + int mypool = 0; + int rank_in_pool = 0; + divide_pools(nproc, myrank, nproc_in_pool, kpar, mypool, rank_in_pool); +#endif + testing::InitGoogleTest(&argc, argv); + const int result = RUN_ALL_TESTS(); +#ifdef __MPI + finishmpi(); +#endif + return result; +} diff --git a/source/source_hsolver/unittests/test_diago_ppcg_parallel.cpp b/source/source_hsolver/unittests/test_diago_ppcg_parallel.cpp new file mode 100644 index 00000000000..f34fe49651b --- /dev/null +++ b/source/source_hsolver/unittests/test_diago_ppcg_parallel.cpp @@ -0,0 +1,415 @@ +/** + * diago_ppcg_parallel_test.cpp — MPI parallel test for DiagoPPCG. + * + * Distributes the rows of a diagonal matrix across MPI processes: each process + * owns a slice of the diagonal and the corresponding rows of psi, computes the + * partial Gram matrix / residual locally, and relies on the solver's pooled + * MPI reductions (reduce_pool) to sum the partial results. The eigenvalues of + * the global diagonal matrix must be recovered identically on every process. + * + * Run with: mpirun -np ./MODULE_HSOLVER_ppcg_parallel + */ + +#include "../diago_ppcg.h" + +#include "source_base/parallel_comm.h" +#include "source_base/parallel_common.h" +#include "source_base/parallel_reduce.h" +#include "source_basis/module_pw/test/test_tool.h" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// A local block can discover two low roots while another block already holds +// an exact higher root. The discarded low direction must survive global root +// selection, including when its rows and projection norms are distributed. +static bool check_competing_roots(const int rank_in_pool, const int nproc_in_pool) +{ + using Scalar = std::complex; + const int bands = 2; + int local_rows = 2; + if (nproc_in_pool == 1) + { + local_rows = 4; + } + const int size = local_rows * bands; + std::vector diagonal(local_rows); + std::vector preconditioner(local_rows, 1.0); + std::vector thresholds(bands, 1.0e-10); + std::vector vectors(size, Scalar(0)); + const double scale = 1.0 / std::sqrt(3.0); + for (int row = 0; row < local_rows; ++row) + { + const int global_row = rank_in_pool * local_rows + row; + if (global_row < 3) + { + diagonal[row] = global_row + 1.0; + } + else + { + diagonal[row] = 97.0 + global_row; + } + if (global_row == 0 || global_row == 1 || global_row == 3) + { + vectors[row] = Scalar(scale); + } + if (global_row == 2) + { + vectors[row + local_rows] = Scalar(1); + } + } + bool operator_capacity_ok = true; + const auto h_op = [&](Scalar* in, Scalar* out, int ld, int ncol) { + operator_capacity_ok = operator_capacity_ok && ncol <= bands; + for (int column = 0; column < ncol; ++column) + { + for (int row = 0; row < local_rows; ++row) + { + out[row + column * ld] = diagonal[row] * in[row + column * ld]; + } + } + }; + hsolver::DiagoPPCG solver(1.0e-10, 80, 1, 4); + std::vector eigenvalues(bands); + solver.diag(h_op, nullptr, local_rows, bands, local_rows, + vectors.data(), eigenvalues.data(), thresholds, preconditioner.data()); + + std::vector residuals(bands, 0.0); + std::vector overlap(bands * bands, Scalar(0)); + for (int column = 0; column < bands; ++column) + { + for (int row = 0; row < local_rows; ++row) + { + const Scalar value = vectors[row + column * local_rows]; + const Scalar residual = (diagonal[row] - eigenvalues[column]) * value; + residuals[column] += std::norm(residual); + for (int other = 0; other < bands; ++other) + { + overlap[other + column * bands] += std::conj(vectors[row + other * local_rows]) * value; + } + } + } + const int residual_count = bands; + Parallel_Reduce::reduce_pool(residuals.data(), residual_count); + const int overlap_size = bands * bands; + Parallel_Reduce::reduce_pool(overlap.data(), overlap_size); + bool correct = solver.converged() && operator_capacity_ok && solver.root_recovery_count() > 0; + for (int column = 0; column < bands; ++column) + { + const double expected = column + 1.0; + correct = correct && std::isfinite(eigenvalues[column]); + correct = correct && std::abs(eigenvalues[column] - expected) < 1.0e-8; + correct = correct && std::isfinite(residuals[column]) && residuals[column] < 1.0e-16; + for (int other = 0; other < bands; ++other) + { + Scalar expected_overlap = Scalar(0); + if (other == column) + { + expected_overlap = Scalar(1); + } + const double error = std::abs(overlap[other + column * bands] - expected_overlap); + correct = correct && std::isfinite(error) && error < 1.0e-8; + } + } + if (!correct) + { + std::fprintf(stderr, "pool rank %d: competing-root regression returned %.12g, %.12g\n", + rank_in_pool, eigenvalues[0], eigenvalues[1]); + } + return correct; +} + +int main(int argc, char** argv) +{ + int nproc = 1; + int myrank = 0; + setupmpi(argc, argv, nproc, myrank); + + int nproc_in_pool = 0; + int kpar = 1; + if (argc > 1 && argv[1][0] != '-') + { + char* end = nullptr; + const long parsed_kpar = std::strtol(argv[1], &end, 10); + if (end != argv[1] && *end == '\0' && parsed_kpar > 0 && parsed_kpar <= INT_MAX) + { + kpar = static_cast(parsed_kpar); + } + } + int mypool = 0; + int rank_in_pool = 0; + divide_pools(nproc, myrank, nproc_in_pool, kpar, mypool, rank_in_pool); + const bool world_root_is_correct = Parallel_Reduce::is_world_root() == (myrank == 0); + BP_WORLD = MPI_COMM_WORLD; + + using T = std::complex; + using Real = double; + + const int nband = 3; + const int n_local = 2 * nband; // each process owns this many rows + const int n_dim_total = nproc_in_pool * n_local; + + // With one pool, only rank zero requests a trace to exercise inconsistent + // process environments. With multiple pools, every rank requests the same + // path to verify that only the MPI world root writes it. + const char* trace_environment = "ABACUS_PPCG_RESIDUAL_TRACE"; + std::string trace_path; + bool trace_path_ok = true; + if (myrank == 0) + { + std::string trace_template = testing::TempDir() + "ppcg_residual_trace_parallel_test_XXXXXX"; + char* trace_template_data = &trace_template[0]; + const int trace_fd = ::mkstemp(trace_template_data); + trace_path_ok = trace_fd >= 0; + if (trace_path_ok) + { + const int close_status = ::close(trace_fd); + char* resolved_path = ::realpath(trace_template_data, nullptr); + trace_path_ok = close_status == 0 && resolved_path != nullptr; + if (trace_path_ok) + { + trace_path = resolved_path; + } + else + { + std::remove(trace_template_data); + } + std::free(resolved_path); + } + } + Parallel_Common::bcast_bool(trace_path_ok); + if (!trace_path_ok) + { + if (myrank == 0) + { + std::fprintf(stderr, "rank 0: cannot create a unique residual trace file\n"); + } + finishmpi(); + return 1; + } + // One root reserves a unique file per launch; all ranks still request the + // same absolute path when the two-pool single-writer check is enabled. + Parallel_Common::bcast_string(trace_path); + int trace_environment_ok = 1; + if (myrank == 0) + { + trace_environment_ok = ::setenv(trace_environment, trace_path.c_str(), 1) == 0; + } + else if (kpar > 1) + { + trace_environment_ok = ::setenv(trace_environment, trace_path.c_str(), 1) == 0; + } + else + { + trace_environment_ok = ::unsetenv(trace_environment) == 0 ? 1 : 0; + } + + // Diagonal H with entries 1, 2, ..., n_dim_total; process owns a slice. + std::vector diag_local(n_local); + std::vector prec(n_local); + for (int i = 0; i < n_local; ++i) + { + diag_local[i] = Real(rank_in_pool * n_local + i + 1); + prec[i] = diag_local[i]; + } + + std::vector ethr(nband, 1e-8); + + // Random initial guess (fixed seed for reproducibility). + std::mt19937 rng(42); + std::uniform_real_distribution dist(-1.0, 1.0); + std::vector psi(n_local * nband, T(0.0, 0.0)); + for (int j = 0; j < nband; ++j) + { + for (int i = 0; i < n_local; ++i) + { + psi[i + j * n_local] = T(dist(rng), 0.0); + } + } + + hsolver::DiagoPPCG solver( + /* diag_thr = */ 1e-12, + /* max_iter = */ 100, + /* sbsize = */ 2, + /* rr_step = */ nband); + + int h_calls = 0; + auto h_op = [&](T* in, T* out, int ld, int ncol) { + ++h_calls; + for (int j = 0; j < ncol; ++j) + { + for (int i = 0; i < n_local; ++i) + { + out[i + j * ld] = diag_local[i] * in[i + j * ld]; + } + } + }; + + std::vector eval(nband, 0.0); + // One invalid rank must stop every rank in its pool before H/S collectives. + // Otherwise the valid ranks would hang while the invalid rank unwinds. + std::vector invalid_prec = prec; + if (rank_in_pool == 0) + { + invalid_prec[0] = std::numeric_limits::quiet_NaN(); + } + bool collective_validation_failed = false; + try + { + solver.diag(h_op, nullptr, n_local, nband, n_local, + psi.data(), eval.data(), ethr, invalid_prec.data()); + } + catch (const std::invalid_argument&) + { + collective_validation_failed = true; + } + const bool validation_before_operator = collective_validation_failed && h_calls == 0; + + // Only one rank advertises an unrepresentable local [X,W,P] workspace. + // The wavefunction product still fits, and the tiny real buffers must not + // be read on that rank or passed to H on any of its peers. + int checked_ld = n_local; + if (rank_in_pool == 0) + { + checked_ld = std::numeric_limits::max() / nband; + } + bool collective_size_failed = false; + try + { + solver.diag(h_op, nullptr, checked_ld, nband, n_local, + psi.data(), eval.data(), ethr, prec.data()); + } + catch (const std::invalid_argument&) + { + collective_size_failed = true; + } + const bool size_validation_before_operator = collective_size_failed && h_calls == 0; + + solver.diag(h_op, nullptr, n_local, nband, n_local, + psi.data(), eval.data(), ethr, prec.data()); + + // The lowest nband eigenvalues of the global diagonal matrix are 1..nband. + int ok = 1; + if (!world_root_is_correct + || !validation_before_operator + || !size_validation_before_operator + || trace_environment_ok != 1) + { + ok = 0; + } + if (!world_root_is_correct) + { + std::fprintf(stderr, "rank %d: MPI world-root detection is inconsistent\n", myrank); + } + if (!validation_before_operator) + { + std::fprintf(stderr, "rank %d: invalid input was not rejected before operator calls\n", myrank); + } + if (!size_validation_before_operator) + { + std::fprintf(stderr, "rank %d: workspace overflow was not rejected before operator calls\n", myrank); + } + if (!solver.converged() + || solver.active_band_count() != 0 + || solver.iteration_limit_reached() + || !std::isfinite(solver.max_residual()) + || solver.max_residual() > 1.0e-8) + { + const int converged = solver.converged() ? 1 : 0; + const int active_band_count = solver.active_band_count(); + const Real max_residual = solver.max_residual(); + std::fprintf(stderr, + "rank %d: solver status is inconsistent (converged=%d, active=%d, max_residual=%.12e)\n", + myrank, + converged, + active_band_count, + max_residual); + ok = 0; + } + for (int i = 0; i < nband; ++i) + { + if (std::abs(eval[i] - Real(i + 1)) > 1e-6) + { + std::printf("rank %d: eval[%d] = %.12f != %d\n", myrank, i, eval[i], i + 1); + ok = 0; + } + } + + if (myrank == 0) + { + std::ifstream trace(trace_path.c_str()); + std::string line; + int header_count = 0; + int record_count = 0; + bool has_final_record = false; + while (std::getline(trace, line)) + { + if (line.find("iteration,stage,max_residual") == 0) + { + ++header_count; + } + else if (!line.empty()) + { + ++record_count; + if (line.find(",final,") != std::string::npos) + { + has_final_record = true; + } + } + } + if (!trace.eof() || header_count != 1 || record_count == 0 || !has_final_record) + { + const int final_record = static_cast(has_final_record); + std::fprintf(stderr, + "rank 0: residual trace is incomplete (headers=%d, records=%d, final=%d)\n", + header_count, + record_count, + final_record); + ok = 0; + } + std::remove(trace_path.c_str()); + } + + // Complete the single-trace assertions before this additional solve writes + // its own trace. Each pool must independently recover both lowest roots. + const int trace_unset_status = ::unsetenv(trace_environment); + const bool competing_roots_ok = check_competing_roots(rank_in_pool, nproc_in_pool); + if (trace_unset_status != 0 || !competing_roots_ok) + { + ok = 0; + } + + int global_ok = ok; + Parallel_Reduce::reduce_min(global_ok); + + if (myrank == 0) + { + if (global_ok == 1) + { + std::printf("PPCG MPI parallel test PASSED\n"); + } + else + { + std::printf("PPCG MPI parallel test FAILED\n"); + } + } + int exit_code = 1; + if (global_ok == 1) + { + exit_code = 0; + } + finishmpi(); + return exit_code; +} diff --git a/source/source_hsolver/unittests/test_hsolver_pw.cpp b/source/source_hsolver/unittests/test_hsolver_pw.cpp index 2086ab34862..67d4d57d936 100644 --- a/source/source_hsolver/unittests/test_hsolver_pw.cpp +++ b/source/source_hsolver/unittests/test_hsolver_pw.cpp @@ -6,11 +6,16 @@ #include "hsolver_pw_sup.h" #include "hsolver_supplementary_mock.h" #include "source_hsolver/diag_comm_info.h" +#include "source_hsolver/diago_params.h" #include "source_hsolver/hs_operator.h" #include "source_hsolver/hsolver_lcaopw.h" #include "source_hsolver/hsolver_pw.h" +#ifdef __MPI +#include "source_basis/module_pw/test/test_tool.h" +#endif #include +#include #include /// H = S = identity: the simplest operator the solvers can be handed @@ -36,6 +41,41 @@ template class IdentityHSOperator : public hsolver::HSOperator, base_device::DEVICE_CPU> +{ + public: + explicit DiagonalHSOperator(const std::vector& diagonal) : diagonal_(diagonal) {} + + void update_k(const int ik) override {} + + void hpsi(const std::complex* x, + std::complex* hx, + const int ld, + const int nvec) const override + { + for (int column = 0; column < nvec; ++column) + { + for (int row = 0; row < ld; ++row) + { + hx[row + column * ld] = diagonal_[row] * x[row + column * ld]; + } + } + } + + void spsi(const std::complex* x, + std::complex* sx, + const int ld, + const int nvec) const override + { + const size_t count = static_cast(ld) * nvec; + const std::complex* end = x + count; + std::copy(x, end, sx); + } + + private: + std::vector diagonal_; +}; + // Mock implementations for the template functions causing linking errors namespace ModulePW { @@ -184,6 +224,24 @@ template void diago_hs_para>( class TestHSolverPW : public ::testing::Test { public: + using DiagoAssist = hsolver::DiagoIterAssist, base_device::DEVICE_CPU>; + + void SetUp() override + { + saved_diag_iter_max = DiagoAssist::PW_DIAG_NMAX; + saved_diag_thr = DiagoAssist::PW_DIAG_THR; + saved_scf_iter = DiagoAssist::SCF_ITER; + saved_need_subspace = DiagoAssist::need_subspace; + } + + void TearDown() override + { + DiagoAssist::PW_DIAG_NMAX = saved_diag_iter_max; + DiagoAssist::PW_DIAG_THR = saved_diag_thr; + DiagoAssist::SCF_ITER = saved_scf_iter; + DiagoAssist::need_subspace = saved_need_subspace; + } + // HSolverPW declares this fixture a friend, but a TEST_F body lives in a // class derived from it and friendship is not inherited, so the call into // the protected hamiltSolvePsiK() is routed through here. @@ -194,14 +252,71 @@ class TestHSolverPW : public ::testing::Test psi::Psi &ps, std::vector::type> &pre, typename GetTypeReal::type *eig, - const int ntry) + const int ntry, + std::ostream& log) { hs.hamiltSolvePsiK( op, ps, pre, eig, - ntry); + ntry, + log); + } + + template + static void set_ethr_band(hsolver::HSolverPW& hs, const int count, const double value) + { + hs.ethr_band.assign(count, value); + } + + std::vector run_ppcg_dispatch(const std::string& calculation, + const bool exact_vectors, + const int iteration_budget, + std::ostream& log) + { + Input_para inp; + inp.calculation = calculation; + inp.ks_solver = "ppcg"; + inp.pw_diag_nmax = iteration_budget; + hsolver::setup_diago_params_pw, base_device::DEVICE_CPU>(0, 1, 1.0e-12, inp); + hsolver::HSolverPW, base_device::DEVICE_CPU> hs( + &pwbk, + calculation, + "pw", + "ppcg", + false, + 1, + DiagoAssist::SCF_ITER, + DiagoAssist::PW_DIAG_NMAX, + DiagoAssist::PW_DIAG_THR, + DiagoAssist::need_subspace, + 2, + false, + 1, + 0, + 0); + EXPECT_EQ(hs.diag_iter_max, iteration_budget); + const std::vector diagonal = {1.0, 2.0, 4.0}; + DiagonalHSOperator op(diagonal); + psi::Psi> ps; + ps.resize(1, 2, 3); + std::complex* vectors = ps.get_pointer(); + std::complex* end = vectors + 6; + const std::complex zero(0.0, 0.0); + std::fill(vectors, end, zero); + vectors[0] = std::complex(1.0, 0.0); + vectors[4] = std::complex(1.0, 0.0); + if (!exact_vectors) + { + vectors[1] = std::complex(1.0, 0.0); + vectors[5] = std::complex(1.0, 0.0); + } + std::vector preconditioner(3, 1.0); + std::vector eigenvalues(2, 0.0); + set_ethr_band(hs, 2, 1.0e-12); + hamiltSolvePsiK(hs, op, ps, preconditioner, eigenvalues.data(), 1, log); + return eigenvalues; } ModulePW::PW_Basis_K pwbk; @@ -253,6 +368,12 @@ class TestHSolverPW : public ::testing::Test std::vector ekb_f; std::ofstream temp_ofs; + + private: + int saved_diag_iter_max = 0; + double saved_diag_thr = 0.0; + int saved_scf_iter = 0; + bool saved_need_subspace = false; }; // TEST_F(TestHSolverPW, solve) { @@ -502,6 +623,7 @@ TEST_F(TestHSolverPW, NpwxLessThanNbandsDeath) psi_test_cd.resize(1, 5, 3); std::vector precond(3, 0.0); std::vector eigenvalues(5, 0.0); + std::ostringstream log; // Expect death from WARNING_QUIT due to npwx < nbands EXPECT_EXIT( hamiltSolvePsiK( @@ -510,6 +632,99 @@ TEST_F(TestHSolverPW, NpwxLessThanNbandsDeath) psi_test_cd, precond, eigenvalues.data(), - 1), + 1, + log), ::testing::ExitedWithCode(1), ".*"); } + +TEST_F(TestHSolverPW, PPCGScfDispatchReportsExhaustedBudget) +{ + std::ostringstream log; + const std::vector eigenvalues = run_ppcg_dispatch("scf", false, 0, log); + EXPECT_NE(log.str().find("PPCG did not converge"), std::string::npos); + EXPECT_NE(log.str().find("active bands=2"), std::string::npos); + EXPECT_NE(log.str().find("max residual="), std::string::npos); + EXPECT_NE(log.str().find("iteration limit reached"), std::string::npos); + for (const double eigenvalue : eigenvalues) + { + EXPECT_TRUE(std::isfinite(eigenvalue)); + } +} + +TEST_F(TestHSolverPW, PPCGNscfDispatchRejectsExhaustedBudget) +{ + std::ostringstream log; + EXPECT_EXIT( + run_ppcg_dispatch("nscf", false, 0, log), + ::testing::ExitedWithCode(1), + "PPCG did not converge.*iteration limit reached"); +} + +TEST_F(TestHSolverPW, PPCGConvergedDispatchHasNoWarning) +{ + std::ostringstream log; + const std::vector eigenvalues = run_ppcg_dispatch("scf", true, 0, log); + EXPECT_EQ(log.str().find("PPCG did not converge"), std::string::npos); + EXPECT_NEAR(eigenvalues[0], 1.0, 1.0e-12); + EXPECT_NEAR(eigenvalues[1], 2.0, 1.0e-12); +} + +TEST_F(TestHSolverPW, PPCGNscfDispatchUsesConfiguredBudget) +{ + DiagoAssist::PW_DIAG_NMAX = 30; + std::ostringstream log; + const std::vector eigenvalues = run_ppcg_dispatch("nscf", false, 77, log); + EXPECT_EQ(log.str().find("PPCG did not converge"), std::string::npos); + EXPECT_NEAR(eigenvalues[0], 1.0, 1.0e-12); + EXPECT_NEAR(eigenvalues[1], 2.0, 1.0e-12); +} + +TEST_F(TestHSolverPW, DiagoParamsPreservesOtherNscfSolverBudgets) +{ + Input_para inp; + inp.calculation = "nscf"; + inp.pw_diag_nmax = 77; + for (const std::string solver : {"cg", "dav", "dav_subspace", "bpcg"}) + { + SCOPED_TRACE(solver); + inp.ks_solver = solver; + DiagoAssist::PW_DIAG_NMAX = 30; + hsolver::setup_diago_params_pw, base_device::DEVICE_CPU>(0, 1, 1.0e-12, inp); + EXPECT_EQ(DiagoAssist::PW_DIAG_NMAX, 30); + } +} + +TEST_F(TestHSolverPW, DiagoParamsSetsScfSolverBudgets) +{ + Input_para inp; + inp.calculation = "scf"; + inp.pw_diag_nmax = 77; + for (const std::string solver : {"cg", "dav", "dav_subspace", "bpcg", "ppcg"}) + { + SCOPED_TRACE(solver); + inp.ks_solver = solver; + DiagoAssist::PW_DIAG_NMAX = 30; + hsolver::setup_diago_params_pw, base_device::DEVICE_CPU>(0, 1, 1.0e-12, inp); + EXPECT_EQ(DiagoAssist::PW_DIAG_NMAX, 77); + } +} + +int main(int argc, char** argv) +{ + ::testing::InitGoogleTest(&argc, argv); +#ifdef __MPI + int nproc = 1; + int myrank = 0; + setupmpi(argc, argv, nproc, myrank); + int nproc_in_pool = 1; + int kpar = nproc; + int mypool = 0; + int rank_in_pool = 0; + divide_pools(nproc, myrank, nproc_in_pool, kpar, mypool, rank_in_pool); +#endif + const int result = RUN_ALL_TESTS(); +#ifdef __MPI + finishmpi(); +#endif + return result; +} diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index 1fa8f482855..b9b89c19571 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -94,7 +94,7 @@ struct Input_para ///< use mesh, which is used in QE. int nspin = 1; ///< LDA ; LSDA ; non-linear spin int pw_diag_nmax = 50; - double pw_diag_thr = 0.01; ///< used in cg method + double pw_diag_thr = 0.01; ///< solver-dependent threshold for iterative PW solvers bool diago_smooth_ethr = false; ///< smooth ethr for iter methods int pw_diag_ndim = 4; ///< dimension of workspace for Davidson diagonalization int diago_cg_prec = 1; ///< mohan add 2012-03-31 diff --git a/source/source_io/module_parameter/read_inp_estruc.cpp b/source/source_io/module_parameter/read_inp_estruc.cpp index 777dbdbfe78..b760c775552 100644 --- a/source/source_io/module_parameter/read_inp_estruc.cpp +++ b/source/source_io/module_parameter/read_inp_estruc.cpp @@ -37,11 +37,12 @@ void ReadInput::item_elec_stru() // Electronic Structure { Input_Item item("ks_solver"); - item.domain.set_choices({"cg", "dav", "bpcg", "dav_subspace", "genelpa", "elpa", "lapack", "scalapack_gvx", "cusolver", "cusolvermp", "pexsi"}); - item.domain.add_requirement("basis_type==pw", "ks_solver in [cg, dav, bpcg, dav_subspace]"); + item.domain.set_choices({"cg", "dav", "bpcg", "dav_subspace", "ppcg", "genelpa", "elpa", "lapack", "scalapack_gvx", "cusolver", "cusolvermp", "pexsi"}); + item.domain.add_requirement("basis_type==pw", "ks_solver in [cg, dav, bpcg, dav_subspace, ppcg]"); item.domain.add_requirement("basis_type==lcao", "ks_solver in [genelpa, elpa, lapack, scalapack_gvx, cusolver, cusolvermp, pexsi]"); + item.domain.add_requirement("ks_solver==ppcg", "esolver_type!=sdft"); item.domain.add_requirement("ks_solver==genelpa", "device==cpu"); - item.annotation = "cg; dav; lapack; genelpa; elpa; scalapack_gvx; cusolver"; + item.annotation = "cg; dav; bpcg; dav_subspace; ppcg; lapack; genelpa; elpa; scalapack_gvx; cusolver"; item.category = "Electronic structure"; item.type = "String"; item.description = R"(Choose the diagonalization methods for the Hamiltonian matrix expanded in a certain basis set. @@ -51,6 +52,7 @@ For plane-wave basis, * cg: The conjugate-gradient (CG) method. * dav: The Davidson algorithm. * dav_subspace: The Davidson algorithm without orthogonalization operation, this method is the most recommended for efficiency. `pw_diag_ndim` can be set to 2 for this method. +* ppcg: The projection preconditioned conjugate-gradient method. It is optimized and validated for CPU plane-wave calculations; non-CPU devices use a transitional host/device bridge. It is not supported for esolver_type=sdft, including mixed stochastic-deterministic calculations. * bpcg: The BPCG method, which is a block-parallel Conjugate Gradient (CG) method, typically exhibits higher acceleration in a GPU environment. The BPCG method is currently under testing and is not recommended for use. For numerical atomic orbitals basis, @@ -1000,10 +1002,10 @@ Use case: When experimental or high-level theoretical results suggest that the S } { Input_Item item("pw_diag_thr"); - item.annotation = "threshold for eigenvalues is cg electron iterations"; + item.annotation = "threshold for iterative PW diagonalization"; item.category = "Plane wave related variables"; item.type = "Real"; - item.description = "Only used when you use ks_solver = cg/dav/dav_subspace/bpcg. It indicates the threshold for the first electronic iteration, from the second iteration the pw_diag_thr will be updated automatically. For nscf calculations with planewave basis set, pw_diag_thr should be <= 1e-3."; + item.description = "Only used when you use ks_solver = cg/dav/dav_subspace/bpcg/ppcg. It indicates the threshold for the first electronic iteration, from the second iteration the pw_diag_thr will be updated automatically. The stopping criterion is solver-dependent: CG and Davidson methods compare eigenvalue changes against the threshold, whereas PPCG requires the norm of H psi - epsilon S psi to be at most max(ethr_band, ethr, 1e-14). The same numerical threshold therefore does not imply the same accuracy or cost across these methods. PPCG keeps the direct residual tolerance because an energy-only or square-root residual tolerance can leave forces and stresses insufficiently converged. For nscf calculations with planewave basis set, pw_diag_thr should be <= 1e-3."; item.default_value = "0.01"; item.unit = ""; read_sync_double(input.pw_diag_thr); @@ -1044,24 +1046,25 @@ Use case: When experimental or high-level theoretical results suggest that the S } { Input_Item item("pw_diag_nmax"); - item.annotation = "max iteration number for cg"; + item.annotation = "maximum iteration count for iterative PW diagonalization"; item.category = "Plane wave related variables"; item.type = "Integer"; - item.description = "Only useful when you use ks_solver = cg/dav/dav_subspace/bpcg. It indicates the maximal iteration number for cg/david/dav_subspace/bpcg method."; + item.description = "Sets the maximum number of iterations for the selected iterative PW diagonalization method. It is used by CG, Davidson, Davidson-subspace, and PPCG; the value is a solver-specific iteration limit and does not define a common number of operator applications or a convergence tolerance. Increasing it does not guarantee convergence. BPCG currently does not use this parameter."; item.default_value = "50"; item.unit = ""; - item.set_availability("basis_type==pw and ks_solver in [cg, dav, dav_subspace, bpcg]"); + item.set_availability("basis_type==pw and ks_solver in [cg, dav, dav_subspace, bpcg, ppcg]"); read_sync_int(input.pw_diag_nmax); this->add_item(item); } { Input_Item item("pw_diag_ndim"); - item.annotation = "dimension of workspace for Davidson diagonalization"; + item.annotation = "dimension of workspace for iterative PW diagonalization"; item.category = "Plane wave related variables"; item.type = "Integer"; - item.description = "Only useful when you use ks_solver = dav or ks_solver = dav_subspace. It indicates dimension of workspace(number of wavefunction packets, at least 2 needed) for the Davidson method. A larger value may yield a smaller number of iterations in the algorithm but uses more memory and more CPU time in subspace diagonalization."; + item.description = "Only useful when you use ks_solver = dav, dav_subspace, or ppcg. It indicates dimension of workspace(number of wavefunction packets, at least 2 needed) for the Davidson method, and the block size for the PPCG method. A larger value may yield a smaller number of iterations in the algorithm but uses more memory and more CPU time in subspace diagonalization."; item.default_value = "4"; item.unit = ""; + item.set_availability("basis_type==pw and ks_solver in [dav, dav_subspace, ppcg]"); read_sync_int(input.pw_diag_ndim); this->add_item(item); } diff --git a/source/source_pw/module_stodft/sto_hsolver_pw.cpp b/source/source_pw/module_stodft/sto_hsolver_pw.cpp index dcb53ec9ddc..e24096ef30c 100644 --- a/source/source_pw/module_stodft/sto_hsolver_pw.cpp +++ b/source/source_pw/module_stodft/sto_hsolver_pw.cpp @@ -67,7 +67,7 @@ void StoHSolverPW::solve(const UnitCell& ucell, this->update_precondition(precondition, ik, this->wfc_basis->npwk[ik], pes->pot->get_vl_of_0()); /// solve eigenvector and eigenvalue for H(k) double* p_eigenvalues = &(pes->ekb(ik, 0)); - this->hamiltSolvePsiK(op, psi, precondition, p_eigenvalues, nks); + this->hamiltSolvePsiK(op, psi, precondition, p_eigenvalues, nks, log); } #ifdef __MPI diff --git a/source/source_pw/module_stodft/test/CMakeLists.txt b/source/source_pw/module_stodft/test/CMakeLists.txt index 534b5143c4c..66a6ac20e88 100644 --- a/source/source_pw/module_stodft/test/CMakeLists.txt +++ b/source/source_pw/module_stodft/test/CMakeLists.txt @@ -28,6 +28,7 @@ addtest( ../sto_hsolver_pw.cpp ../../../source_hsolver/hsolver_pw.cpp ../../../source_hsolver/diago_bpcg.cpp + ../../../source_hsolver/diago_ppcg.cpp ../../../source_hsolver/diago_dav_subspace.cpp ../../../source_hsolver/diag_const_nums.cpp ../../../source_hsolver/diago_iter_assist.cpp diff --git a/tests/01_PW/817_PW_PPCG/INPUT b/tests/01_PW/817_PW_PPCG/INPUT new file mode 100644 index 00000000000..49cf788d32d --- /dev/null +++ b/tests/01_PW/817_PW_PPCG/INPUT @@ -0,0 +1,32 @@ +INPUT_PARAMETERS +#Parameters (General) +suffix autotest +pseudo_dir ../../PP_ORB +pw_seed 1 + +gamma_only 0 +calculation scf +symmetry 1 +out_level ie +smearing_method gaussian +smearing_sigma 0.02 + +#Parameters (3.PW) +ecutwfc 40 +scf_thr 1e-6 +scf_nmax 50 + +#Parameters (LCAO) +basis_type pw +ks_solver ppcg +device cpu +chg_extrap second-order +pw_diag_thr 0.00001 +pw_diag_ndim 4 + +cal_force 1 +cal_stress 1 + +mixing_type broyden +mixing_beta 0.4 +mixing_gg0 1.5 diff --git a/tests/01_PW/817_PW_PPCG/KPT b/tests/01_PW/817_PW_PPCG/KPT new file mode 100644 index 00000000000..b5b3bdb1ae2 --- /dev/null +++ b/tests/01_PW/817_PW_PPCG/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 2 0 0 0 diff --git a/tests/01_PW/817_PW_PPCG/README b/tests/01_PW/817_PW_PPCG/README new file mode 100644 index 00000000000..c3a6155993b --- /dev/null +++ b/tests/01_PW/817_PW_PPCG/README @@ -0,0 +1 @@ +pw basis for GaAs with ks_solver ppcg (projection preconditioned conjugate-gradient), multi k diff --git a/tests/01_PW/817_PW_PPCG/STRU b/tests/01_PW/817_PW_PPCG/STRU new file mode 100644 index 00000000000..b03baadd25e --- /dev/null +++ b/tests/01_PW/817_PW_PPCG/STRU @@ -0,0 +1,23 @@ +ATOMIC_SPECIES +As 1 As_dojo.upf upf201 +Ga 1 Ga_dojo.upf upf201 + +LATTICE_CONSTANT +1 // add lattice constant, 10.58 ang + +LATTICE_VECTORS +5.33 5.33 0.0 +0.0 5.33 5.33 +5.33 0.0 5.33 +ATOMIC_POSITIONS +Direct //Cartesian or Direct coordinate. + +As +0 +1 +0.300000 0.3300000 0.27000000 0 0 0 + +Ga //Element Label +0 +1 //number of atom +0.00000 0.00000 0.000000 0 0 0 diff --git a/tests/01_PW/817_PW_PPCG/result.ref b/tests/01_PW/817_PW_PPCG/result.ref new file mode 100644 index 00000000000..2bb5340777f --- /dev/null +++ b/tests/01_PW/817_PW_PPCG/result.ref @@ -0,0 +1,8 @@ +etotref -4862.3309704772291298 +etotperatomref -2431.1654852386 +totalforceref 9.105146 +totalstressref 37223.049159 +pointgroupref C_1 +spacegroupref C_1 +nksibzref 2 +totaltimeref 4.64 diff --git a/tests/01_PW/818_PW_PPCG_Si/INPUT b/tests/01_PW/818_PW_PPCG_Si/INPUT new file mode 100644 index 00000000000..38436bc1e41 --- /dev/null +++ b/tests/01_PW/818_PW_PPCG_Si/INPUT @@ -0,0 +1,33 @@ +INPUT_PARAMETERS +#Parameters (General) +suffix autotest +pseudo_dir ../../PP_ORB +pw_seed 1 + +gamma_only 0 +calculation scf +symmetry 1 +out_level ie +smearing_method gaussian +smearing_sigma 0.02 + +#Parameters (3.PW) +ecutwfc 40 +scf_thr 1e-6 +scf_nmax 50 + +#Parameters (LCAO) +basis_type pw +ks_solver ppcg +device cpu +nbands 6 +chg_extrap second-order +pw_diag_thr 0.00001 +pw_diag_ndim 4 + +cal_force 1 +cal_stress 1 + +mixing_type broyden +mixing_beta 0.4 +mixing_gg0 1.5 diff --git a/tests/01_PW/818_PW_PPCG_Si/KPT b/tests/01_PW/818_PW_PPCG_Si/KPT new file mode 100644 index 00000000000..28006d5e2df --- /dev/null +++ b/tests/01_PW/818_PW_PPCG_Si/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +2 2 2 0 0 0 diff --git a/tests/01_PW/818_PW_PPCG_Si/README b/tests/01_PW/818_PW_PPCG_Si/README new file mode 100644 index 00000000000..5f6c7b60ade --- /dev/null +++ b/tests/01_PW/818_PW_PPCG_Si/README @@ -0,0 +1 @@ +pw basis for Si with ks_solver ppcg (projection preconditioned conjugate-gradient), multi k diff --git a/tests/01_PW/818_PW_PPCG_Si/STRU b/tests/01_PW/818_PW_PPCG_Si/STRU new file mode 100644 index 00000000000..f3bade0c3b3 --- /dev/null +++ b/tests/01_PW/818_PW_PPCG_Si/STRU @@ -0,0 +1,19 @@ +ATOMIC_SPECIES +Si 1.000 Si_dojo_nsoc.upf + +LATTICE_CONSTANT +10.2 + +LATTICE_VECTORS +0.0 0.5 0.5 +0.5 0.0 0.5 +0.5 0.5 0.0 + +ATOMIC_POSITIONS +Direct + +Si +0.0 +2 +0.00 0.00 0.00 0 0 0 +0.25 0.25 0.25 0 0 0 diff --git a/tests/01_PW/818_PW_PPCG_Si/result.ref b/tests/01_PW/818_PW_PPCG_Si/result.ref new file mode 100644 index 00000000000..bd11216ce9f --- /dev/null +++ b/tests/01_PW/818_PW_PPCG_Si/result.ref @@ -0,0 +1,8 @@ +etotref -227.6796325586651903 +etotperatomref -113.8398162793 +totalforceref 0.000000 +totalstressref 342.892641 +pointgroupref T_d +spacegroupref O_h +nksibzref 3 +totaltimeref 1.45 diff --git a/tests/01_PW/819_PW_PPCG_Al/INPUT b/tests/01_PW/819_PW_PPCG_Al/INPUT new file mode 100644 index 00000000000..9e6b92cc088 --- /dev/null +++ b/tests/01_PW/819_PW_PPCG_Al/INPUT @@ -0,0 +1,36 @@ +INPUT_PARAMETERS +#Parameters (General) +suffix autotest +pseudo_dir ../../PP_ORB +pw_seed 1 + +gamma_only 0 +calculation scf +symmetry 1 +out_level ie +smearing_method gaussian +smearing_sigma 0.02 + +#Parameters (3.PW) +ecutwfc 40 +# Resolve stress with a real-space density criterion; see README. +scf_thr_type 2 +scf_thr 1e-10 +scf_nmax 50 + +#Parameters (LCAO) +basis_type pw +ks_solver ppcg +device cpu +nbands 8 +chg_extrap second-order +pw_diag_thr 0.00001 +pw_diag_nmax 200 +pw_diag_ndim 4 + +cal_force 1 +cal_stress 1 + +mixing_type broyden +mixing_beta 0.4 +mixing_gg0 1.5 diff --git a/tests/01_PW/819_PW_PPCG_Al/KPT b/tests/01_PW/819_PW_PPCG_Al/KPT new file mode 100644 index 00000000000..2df8740b6fa --- /dev/null +++ b/tests/01_PW/819_PW_PPCG_Al/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +4 4 4 0 0 0 diff --git a/tests/01_PW/819_PW_PPCG_Al/README b/tests/01_PW/819_PW_PPCG_Al/README new file mode 100644 index 00000000000..31f82ae2382 --- /dev/null +++ b/tests/01_PW/819_PW_PPCG_Al/README @@ -0,0 +1,7 @@ +Plane-wave Al metal case using `ks_solver=ppcg`. + +The case checks the final energy, forces, stress, symmetry and irreducible +k-point count with the standard integration-test reference. The real-space +SCF density criterion and the 200-step inner budget keep the stress check +converged while leaving the PPCG stopping rule unchanged. The case is a +correctness regression; its recorded wall time is not a performance claim. diff --git a/tests/01_PW/819_PW_PPCG_Al/STRU b/tests/01_PW/819_PW_PPCG_Al/STRU new file mode 100644 index 00000000000..1db1a774bdf --- /dev/null +++ b/tests/01_PW/819_PW_PPCG_Al/STRU @@ -0,0 +1,18 @@ +ATOMIC_SPECIES +Al 26.98 Al_ONCV_PBE-1.0.upf + +LATTICE_CONSTANT +7.63 + +LATTICE_VECTORS +0.0 0.5 0.5 +0.5 0.0 0.5 +0.5 0.5 0.0 + +ATOMIC_POSITIONS +Direct + +Al +0.0 +1 +0.00 0.00 0.00 1 1 1 diff --git a/tests/01_PW/819_PW_PPCG_Al/result.ref b/tests/01_PW/819_PW_PPCG_Al/result.ref new file mode 100644 index 00000000000..5183a1facb7 --- /dev/null +++ b/tests/01_PW/819_PW_PPCG_Al/result.ref @@ -0,0 +1,8 @@ +etotref -1882.1779472340826942 +etotperatomref -1882.1779472341 +totalforceref 0.000000 +totalstressref 1107.675822 +pointgroupref O_h +spacegroupref O_h +nksibzref 8 +totaltimeref 1.84 diff --git a/tests/01_PW/CASES_CPU.txt b/tests/01_PW/CASES_CPU.txt index 1b8261b7a8a..b263a49b980 100644 --- a/tests/01_PW/CASES_CPU.txt +++ b/tests/01_PW/CASES_CPU.txt @@ -122,6 +122,9 @@ scf_brav_bco_smear_mv 814_PW_LT_triclinic relax_dftu_spin2 relax_dftu_spin4 +817_PW_PPCG +818_PW_PPCG_Si +819_PW_PPCG_Al scf_deltaspin2 scf_deltaspin4 nscf_hse_same diff --git a/tests/11_PW_GPU/CASES_GPU.txt b/tests/11_PW_GPU/CASES_GPU.txt index a7a12920b4d..68969d84535 100644 --- a/tests/11_PW_GPU/CASES_GPU.txt +++ b/tests/11_PW_GPU/CASES_GPU.txt @@ -3,6 +3,7 @@ scf_cg scf_cg_single scf_dav scf_dav_sub +scf_ppcg scf_out_wf scf_out_wf_norm scf_out_wf_spinor diff --git a/tests/11_PW_GPU/scf_ppcg/INPUT b/tests/11_PW_GPU/scf_ppcg/INPUT new file mode 100644 index 00000000000..9af711ac742 --- /dev/null +++ b/tests/11_PW_GPU/scf_ppcg/INPUT @@ -0,0 +1,34 @@ +INPUT_PARAMETERS +#Parameters (General) +suffix autotest +pseudo_dir ../../PP_ORB + +gamma_only 0 +calculation scf +symmetry 1 +relax_nmax 1 +out_level ie +smearing_method gaussian +smearing_sigma 0.02 + +#Parameters (3.PW) +ecutwfc 40 +scf_thr 1e-7 +scf_nmax 100 + +#Parameters (LCAO) +basis_type pw +ks_solver ppcg +device gpu +chg_extrap second-order +pw_diag_thr 0.00001 +pw_diag_ndim 4 + +cal_force 1 +cal_stress 1 + +mixing_type broyden +mixing_beta 0.4 +mixing_gg0 1.5 + +pw_seed 1 diff --git a/tests/11_PW_GPU/scf_ppcg/KPT b/tests/11_PW_GPU/scf_ppcg/KPT new file mode 100644 index 00000000000..28006d5e2df --- /dev/null +++ b/tests/11_PW_GPU/scf_ppcg/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +2 2 2 0 0 0 diff --git a/tests/11_PW_GPU/scf_ppcg/README b/tests/11_PW_GPU/scf_ppcg/README new file mode 100644 index 00000000000..1457417f2cd --- /dev/null +++ b/tests/11_PW_GPU/scf_ppcg/README @@ -0,0 +1,7 @@ +GPU smoke case for the PPCG host/device bridge. + +It uses the same GaAs structure as the CPU PPCG case, but samples a 2x2x2 +Gamma-centered k-point mesh instead of the CPU case's 1x1x2 mesh. It checks the +standard energy, force and stress reference with the existing GPU acceptance +limits. Run it on a CUDA-enabled build with one OpenMP/BLAS thread; no CPU +timing or GPU speedup is inferred from this fixture. diff --git a/tests/11_PW_GPU/scf_ppcg/STRU b/tests/11_PW_GPU/scf_ppcg/STRU new file mode 100644 index 00000000000..b03baadd25e --- /dev/null +++ b/tests/11_PW_GPU/scf_ppcg/STRU @@ -0,0 +1,23 @@ +ATOMIC_SPECIES +As 1 As_dojo.upf upf201 +Ga 1 Ga_dojo.upf upf201 + +LATTICE_CONSTANT +1 // add lattice constant, 10.58 ang + +LATTICE_VECTORS +5.33 5.33 0.0 +0.0 5.33 5.33 +5.33 0.0 5.33 +ATOMIC_POSITIONS +Direct //Cartesian or Direct coordinate. + +As +0 +1 +0.300000 0.3300000 0.27000000 0 0 0 + +Ga //Element Label +0 +1 //number of atom +0.00000 0.00000 0.000000 0 0 0 diff --git a/tests/11_PW_GPU/scf_ppcg/result.ref b/tests/11_PW_GPU/scf_ppcg/result.ref new file mode 100644 index 00000000000..d851d7baa46 --- /dev/null +++ b/tests/11_PW_GPU/scf_ppcg/result.ref @@ -0,0 +1,8 @@ +etotref -4869.7470517268475305 +etotperatomref -2434.8735258634 +totalforceref 5.206812 +totalstressref 37241.543010 +pointgroupref C_1 +spacegroupref C_1 +nksibzref 8 +totaltimeref 4.30 diff --git a/tests/11_PW_GPU/scf_ppcg/threshold b/tests/11_PW_GPU/scf_ppcg/threshold new file mode 100644 index 00000000000..a044ada1bc4 --- /dev/null +++ b/tests/11_PW_GPU/scf_ppcg/threshold @@ -0,0 +1,4 @@ +threshold 0.000001 +force_threshold 0.02 +stress_threshold 0.5 +fatal_threshold 1