From fd2fe30c6f324d661e8dd734a23846f51f88487c Mon Sep 17 00:00:00 2001 From: Chen Nuo <49788094+Cstandardlib@users.noreply.github.com> Date: Thu, 18 Dec 2025 23:51:15 +0800 Subject: [PATCH 1/8] Call output_iterInfo in hsolver --- source/source_hsolver/hsolver_pw.cpp | 41 ++++++++++++++++------------ 1 file changed, 23 insertions(+), 18 deletions(-) diff --git a/source/source_hsolver/hsolver_pw.cpp b/source/source_hsolver/hsolver_pw.cpp index f3e5a7ca60..5577907018 100644 --- a/source/source_hsolver/hsolver_pw.cpp +++ b/source/source_hsolver/hsolver_pw.cpp @@ -105,7 +105,7 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, for (int i = 0; i < this->wfc_basis->nks; ++i) { const int ik = k_order[i]; - + // update H(k) for each k point pHamilt->updateHk(ik); @@ -140,6 +140,7 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, // solve eigenvector and eigenvalue for H(k) this->hamiltSolvePsiK(pHamilt, psi, precondition, eigenvalues.data() + ik * psi.get_nbands(), this->wfc_basis->nks); + // output iteration information and reset avg_iter if (skip_charge) { GlobalV::ofs_running << "Average iterative diagonalization steps for k-points " << ik @@ -147,6 +148,10 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, << " ; where current threshold is: " << this->diag_thr << " . " << std::endl; DiagoIterAssist::avg_iter = 0.0; } + else + { + this->output_iterInfo(); + } } } else { @@ -192,7 +197,7 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, /// calculate the contribution of Psi for charge density rho } } - + count++; // END Loop over k points @@ -533,39 +538,39 @@ void HSolverPW::build_k_neighbors() { kvecs_c.resize(nk); k_order.clear(); k_order.reserve(nk); - + // Store k-points and corresponding indices struct KPoint { ModuleBase::Vector3 kvec; int index; double norm; - - KPoint(const ModuleBase::Vector3& v, int i) : + + KPoint(const ModuleBase::Vector3& v, int i) : kvec(v), index(i), norm(v.norm()) {} }; - + // Build k-point list std::vector klist; for (int ik = 0; ik < nk; ++ik) { kvecs_c[ik] = this->wfc_basis->kvec_c[ik]; klist.push_back(KPoint(kvecs_c[ik], ik)); } - + // Sort k-points by distance from origin std::sort(klist.begin(), klist.end(), [](const KPoint& a, const KPoint& b) { return a.norm < b.norm; }); - + // Build parent-child relationships k_order.push_back(klist[0].index); - + // Find nearest processed k-point as parent for each k-point for (int i = 1; i < nk; ++i) { int current_k = klist[i].index; double min_dist = 1e10; int parent = -1; - + // find the nearest k-point as parent for (int j = 0; j < k_order.size(); ++j) { int processed_k = k_order[j]; @@ -575,7 +580,7 @@ void HSolverPW::build_k_neighbors() { parent = processed_k; } } - + k_parent[current_k] = parent; k_order.push_back(current_k); } @@ -585,34 +590,34 @@ template void HSolverPW::propagate_psi(psi::Psi& psi, const int from_ik, const int to_ik) { const int nbands = psi.get_nbands(); const int npwk = this->wfc_basis->npwk[to_ik]; - + // Get k-point difference ModuleBase::Vector3 dk = kvecs_c[to_ik] - kvecs_c[from_ik]; - + // Allocate porter locally T* porter = nullptr; resmem_complex_op()(porter, this->wfc_basis->nmaxgr, "HSolverPW::porter"); - + // Process each band for (int ib = 0; ib < nbands; ib++) { // Fix current k-point and band // psi.fix_k(from_ik); - + // FFT to real space // this->wfc_basis->recip_to_real(this->ctx, psi.get_pointer(ib), porter, from_ik); this->wfc_basis->recip_to_real(this->ctx, &psi(from_ik, ib, 0), porter, from_ik); - + // Apply phase factor // // TODO: Check how to get the r vector // ModuleBase::Vector3 r = this->wfc_basis->get_ir2r(ir); // double phase = this->wfc_basis->tpiba * (dk.x * r.x + dk.y * r.y + dk.z * r.z); // psi_real[ir] *= std::exp(std::complex(0.0, phase)); // } - + // Fix k-point for target // psi.fix_k(to_ik); - + // FFT back to reciprocal space // this->wfc_basis->real_to_recip(this->ctx, porter, psi.get_pointer(ib), to_ik, true); this->wfc_basis->real_to_recip(this->ctx, porter, &psi(to_ik, ib, 0), to_ik); From 4b4d689984094f472772c40d54bb6bea9a630301 Mon Sep 17 00:00:00 2001 From: Chen Nuo <49788094+Cstandardlib@users.noreply.github.com> Date: Fri, 19 Dec 2025 00:06:17 +0800 Subject: [PATCH 2/8] Add average iter output for cg --- source/source_hsolver/diago_cg.cpp | 6 ++++-- source/source_hsolver/diago_cg.h | 21 +++++++++++---------- source/source_hsolver/hsolver_pw.cpp | 4 +++- source/source_hsolver/test/hsolver_pw_sup.h | 4 ++-- 4 files changed, 20 insertions(+), 15 deletions(-) diff --git a/source/source_hsolver/diago_cg.cpp b/source/source_hsolver/diago_cg.cpp index 1a72774303..dbb3806762 100644 --- a/source/source_hsolver/diago_cg.cpp +++ b/source/source_hsolver/diago_cg.cpp @@ -575,7 +575,7 @@ bool DiagoCG::test_exit_cond(const int& ntry, const int& notconv) con } template -void DiagoCG::diag(const Func& hpsi_func, +double DiagoCG::diag(const Func& hpsi_func, const Func& spsi_func, ct::Tensor& psi, ct::Tensor& eigen, @@ -626,6 +626,8 @@ void DiagoCG::diag(const Func& hpsi_func, psi.zero(); // copy psi_temp to psi for 0 to npw. psi.sync(psi_temp); + + return avg_iter_; } namespace hsolver @@ -644,4 +646,4 @@ template class DiagoCG; template class DiagoCG; #endif #endif -} // namespace hsolver \ No newline at end of file +} // namespace hsolver diff --git a/source/source_hsolver/diago_cg.h b/source/source_hsolver/diago_cg.h index 20a5890552..ebf984c113 100644 --- a/source/source_hsolver/diago_cg.h +++ b/source/source_hsolver/diago_cg.h @@ -35,13 +35,14 @@ class DiagoCG final const Real& pw_diag_thr, const int& pw_diag_nmax, const int& nproc_in_pool); - + ~DiagoCG(); // virtual void init(){}; // refactor hpsi_info // this is the diag() function for CG method - void diag(const Func& hpsi_func, + // returns avg_iter + double diag(const Func& hpsi_func, const Func& spsi_func, ct::Tensor& psi, ct::Tensor& eigen, @@ -59,7 +60,7 @@ class DiagoCG final /// col size for input psi matrix int n_basis_ = 0; /// average iteration steps for cg diagonalization - int avg_iter_ = 0; + double avg_iter_ = 0; /// threshold for cg diagonalization Real pw_diag_thr_ = 1e-5; /// maximum iteration steps for cg diagonalization @@ -87,15 +88,15 @@ class DiagoCG final ct::Tensor& pphi); void orth_grad( - const ct::Tensor& psi, - const int& m, - ct::Tensor& grad, + const ct::Tensor& psi, + const int& m, + ct::Tensor& grad, ct::Tensor& scg, ct::Tensor& lagrange); void calc_gamma_cg( const int& iter, - const Real& cg_norm, + const Real& cg_norm, const Real& theta, const ct::Tensor& prec, const ct::Tensor& scg, @@ -110,8 +111,8 @@ class DiagoCG final const ct::Tensor& cg, const ct::Tensor& scg, const double& ethreshold, - Real &cg_norm, - Real &theta, + Real &cg_norm, + Real &theta, Real &eigen, ct::Tensor& phi_m, ct::Tensor& sphi, @@ -133,4 +134,4 @@ class DiagoCG final } // namespace hsolver -#endif // MODULE_HSOLVER_DIAGO_CG_H_ \ No newline at end of file +#endif // MODULE_HSOLVER_DIAGO_CG_H_ diff --git a/source/source_hsolver/hsolver_pw.cpp b/source/source_hsolver/hsolver_pw.cpp index 5577907018..e3f23f9c07 100644 --- a/source/source_hsolver/hsolver_pw.cpp +++ b/source/source_hsolver/hsolver_pw.cpp @@ -346,7 +346,9 @@ void HSolverPW::hamiltSolvePsiK(hamilt::Hamilt* hm, .to_device() .slice({0}, {psi.get_current_ngk()}); - cg.diag(hpsi_func, spsi_func, psi_tensor, eigen_tensor, this->ethr_band, prec_tensor); + DiagoIterAssist::avg_iter += static_cast( + cg.diag(hpsi_func, spsi_func, psi_tensor, eigen_tensor, this->ethr_band, prec_tensor) + ); // TODO: Double check tensormap's potential problem // ct::TensorMap(psi.get_pointer(), psi_tensor, {psi.get_nbands(), psi.get_nbasis()}).sync(psi_tensor); } diff --git a/source/source_hsolver/test/hsolver_pw_sup.h b/source/source_hsolver/test/hsolver_pw_sup.h index 07cc650013..fb3757a08b 100644 --- a/source/source_hsolver/test/hsolver_pw_sup.h +++ b/source/source_hsolver/test/hsolver_pw_sup.h @@ -92,7 +92,7 @@ DiagoCG::~DiagoCG() { } template -void DiagoCG::diag(const Func& hpsi_func, +double DiagoCG::diag(const Func& hpsi_func, const Func& spsi_func, ct::Tensor& psi, ct::Tensor& eigen, @@ -112,7 +112,7 @@ void DiagoCG::diag(const Func& hpsi_func, eigen_pack[ib] /= n_basis; } DiagoIterAssist::avg_iter += 1.0; - return; + return avg_iter_; } template class DiagoCG, base_device::DEVICE_CPU>; From 9da27335c12fe0cba0ef95972306defb154da5dd Mon Sep 17 00:00:00 2001 From: Chen Nuo <49788094+Cstandardlib@users.noreply.github.com> Date: Fri, 19 Dec 2025 01:01:02 +0800 Subject: [PATCH 3/8] Fix average iter output --- source/source_hsolver/diago_cg.cpp | 13 +++++++++++++ source/source_hsolver/diago_cg.h | 3 +++ source/source_hsolver/hsolver_pw.cpp | 16 ++++++++-------- 3 files changed, 24 insertions(+), 8 deletions(-) diff --git a/source/source_hsolver/diago_cg.cpp b/source/source_hsolver/diago_cg.cpp index dbb3806762..7e1623ff66 100644 --- a/source/source_hsolver/diago_cg.cpp +++ b/source/source_hsolver/diago_cg.cpp @@ -171,6 +171,7 @@ void DiagoCG::diag_once(const ct::Tensor& prec_in, { ++this->notconv_; } + iter_band.push_back(iter); avg += static_cast(iter) + 1.00; // reorder eigenvalue if they are not in the right order @@ -627,6 +628,18 @@ double DiagoCG::diag(const Func& hpsi_func, // copy psi_temp to psi for 0 to npw. psi.sync(psi_temp); +#ifdef __DEBUG +// only output iter count for each band if DEBUG! +// this should not be output in production log + std::cout << "\n DiagoCG::diag' avg_iter_ = " << avg_iter_; + std::cout << "\n DiagoCG::diag' iter_band = "; + for (auto iter : iter_band) + { + std::cout << iter << " "; + } + std::cout << "\n"; +#endif + return avg_iter_; } diff --git a/source/source_hsolver/diago_cg.h b/source/source_hsolver/diago_cg.h index ebf984c113..bf03cb5850 100644 --- a/source/source_hsolver/diago_cg.h +++ b/source/source_hsolver/diago_cg.h @@ -2,6 +2,7 @@ #define MODULE_HSOLVER_DIAGO_CG_H_ #include +#include #include #include @@ -61,6 +62,8 @@ class DiagoCG final int n_basis_ = 0; /// average iteration steps for cg diagonalization double avg_iter_ = 0; + /// std::vector for iter count of each band + std::vector iter_band; /// threshold for cg diagonalization Real pw_diag_thr_ = 1e-5; /// maximum iteration steps for cg diagonalization diff --git a/source/source_hsolver/hsolver_pw.cpp b/source/source_hsolver/hsolver_pw.cpp index e3f23f9c07..39f152c5f5 100644 --- a/source/source_hsolver/hsolver_pw.cpp +++ b/source/source_hsolver/hsolver_pw.cpp @@ -140,7 +140,6 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, // solve eigenvector and eigenvalue for H(k) this->hamiltSolvePsiK(pHamilt, psi, precondition, eigenvalues.data() + ik * psi.get_nbands(), this->wfc_basis->nks); - // output iteration information and reset avg_iter if (skip_charge) { GlobalV::ofs_running << "Average iterative diagonalization steps for k-points " << ik @@ -148,12 +147,8 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, << " ; where current threshold is: " << this->diag_thr << " . " << std::endl; DiagoIterAssist::avg_iter = 0.0; } - else - { - this->output_iterInfo(); - } } - } + } // if (use_k_continuity) else { // Original code without k-point continuity for (int ik = 0; ik < this->wfc_basis->nks; ++ik) @@ -187,6 +182,7 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, // solve eigenvector and eigenvalue for H(k) this->hamiltSolvePsiK(pHamilt, psi, precondition, eigenvalues.data() + ik * psi.get_nbands(), this->wfc_basis->nks); + // output iteration information and reset avg_iter if (skip_charge) { GlobalV::ofs_running << " k(" << ik+1 << "/" << pes->klist->get_nkstot() @@ -194,9 +190,13 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, << " threshold=" << this->diag_thr << std::endl; DiagoIterAssist::avg_iter = 0.0; } + /// calculate the contribution of Psi for charge density rho } - } + } // else (use_k_continuity) + + // output iteration information and reset avg_iter + this->output_iterInfo(); count++; // END Loop over k points @@ -528,7 +528,7 @@ void HSolverPW::output_iterInfo() { GlobalV::ofs_running << "Average iterative diagonalization steps: " << DiagoIterAssist::avg_iter / this->wfc_basis->nks - << " ; where current threshold is: " << this->diag_thr << " . " << std::endl; + << "; where current threshold is: " << this->diag_thr << ". " << std::endl; // reset avg_iter DiagoIterAssist::avg_iter = 0.0; } From 391ad61619696350bbb84c8d273dcc06f92f5a92 Mon Sep 17 00:00:00 2001 From: Chen Nuo <49788094+Cstandardlib@users.noreply.github.com> Date: Fri, 19 Dec 2025 01:28:28 +0800 Subject: [PATCH 4/8] Add comment --- source/source_hsolver/hsolver_pw.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/source_hsolver/hsolver_pw.cpp b/source/source_hsolver/hsolver_pw.cpp index 39f152c5f5..e524c85401 100644 --- a/source/source_hsolver/hsolver_pw.cpp +++ b/source/source_hsolver/hsolver_pw.cpp @@ -144,7 +144,7 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, { GlobalV::ofs_running << "Average iterative diagonalization steps for k-points " << ik << " is: " << DiagoIterAssist::avg_iter - << " ; where current threshold is: " << this->diag_thr << " . " << std::endl; + << "; where current threshold is: " << this->diag_thr << ". " << std::endl; DiagoIterAssist::avg_iter = 0.0; } } @@ -195,7 +195,7 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, } } // else (use_k_continuity) - // output iteration information and reset avg_iter + // output average iteration information and reset avg_iter this->output_iterInfo(); count++; From ff7e3a24cdb59d9c3a1d6bbf89b47c0e3faf7dd8 Mon Sep 17 00:00:00 2001 From: Chen Nuo <49788094+Cstandardlib@users.noreply.github.com> Date: Fri, 19 Dec 2025 10:26:34 +0800 Subject: [PATCH 5/8] Align format --- source/source_hsolver/diago_cg.cpp | 4 ++-- source/source_hsolver/hsolver_pw.cpp | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/source/source_hsolver/diago_cg.cpp b/source/source_hsolver/diago_cg.cpp index 7e1623ff66..564c36e74b 100644 --- a/source/source_hsolver/diago_cg.cpp +++ b/source/source_hsolver/diago_cg.cpp @@ -633,9 +633,9 @@ double DiagoCG::diag(const Func& hpsi_func, // this should not be output in production log std::cout << "\n DiagoCG::diag' avg_iter_ = " << avg_iter_; std::cout << "\n DiagoCG::diag' iter_band = "; - for (auto iter : iter_band) + for (auto iter_in_band : iter_band) { - std::cout << iter << " "; + std::cout << iter_in_band << " "; } std::cout << "\n"; #endif diff --git a/source/source_hsolver/hsolver_pw.cpp b/source/source_hsolver/hsolver_pw.cpp index e524c85401..6f828a908c 100644 --- a/source/source_hsolver/hsolver_pw.cpp +++ b/source/source_hsolver/hsolver_pw.cpp @@ -142,7 +142,7 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, if (skip_charge) { - GlobalV::ofs_running << "Average iterative diagonalization steps for k-points " << ik + GlobalV::ofs_running << " Average iterative diagonalization steps for k-points " << ik << " is: " << DiagoIterAssist::avg_iter << "; where current threshold is: " << this->diag_thr << ". " << std::endl; DiagoIterAssist::avg_iter = 0.0; @@ -526,9 +526,9 @@ void HSolverPW::output_iterInfo() // in PW base, average iteration steps for each band and k-point should be printing if (DiagoIterAssist::avg_iter > 0.0) { - GlobalV::ofs_running << "Average iterative diagonalization steps: " + GlobalV::ofs_running << " Average iterative diagonalization steps: " << DiagoIterAssist::avg_iter / this->wfc_basis->nks - << "; where current threshold is: " << this->diag_thr << ". " << std::endl; + << ";\n where current threshold is: " << this->diag_thr << ". " << std::endl; // reset avg_iter DiagoIterAssist::avg_iter = 0.0; } From 882685a0cdefe8139bdaff84cddc6973a6a27457 Mon Sep 17 00:00:00 2001 From: Chen Nuo <49788094+Cstandardlib@users.noreply.github.com> Date: Fri, 19 Dec 2025 10:33:17 +0800 Subject: [PATCH 6/8] Update output format --- source/source_hsolver/hsolver_pw.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/source_hsolver/hsolver_pw.cpp b/source/source_hsolver/hsolver_pw.cpp index 6f828a908c..026165ccfc 100644 --- a/source/source_hsolver/hsolver_pw.cpp +++ b/source/source_hsolver/hsolver_pw.cpp @@ -526,7 +526,7 @@ void HSolverPW::output_iterInfo() // in PW base, average iteration steps for each band and k-point should be printing if (DiagoIterAssist::avg_iter > 0.0) { - GlobalV::ofs_running << " Average iterative diagonalization steps: " + GlobalV::ofs_running << " Average iterative diagonalization steps for k-points: " << DiagoIterAssist::avg_iter / this->wfc_basis->nks << ";\n where current threshold is: " << this->diag_thr << ". " << std::endl; // reset avg_iter From ea8216f5c6480138f51bddb31262bdb3c6e7775d Mon Sep 17 00:00:00 2001 From: Chen Nuo <49788094+Cstandardlib@users.noreply.github.com> Date: Fri, 19 Dec 2025 10:49:44 +0800 Subject: [PATCH 7/8] Make clear output format --- source/source_hsolver/hsolver_pw.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/source/source_hsolver/hsolver_pw.cpp b/source/source_hsolver/hsolver_pw.cpp index 026165ccfc..43f027ed40 100644 --- a/source/source_hsolver/hsolver_pw.cpp +++ b/source/source_hsolver/hsolver_pw.cpp @@ -143,8 +143,8 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, if (skip_charge) { GlobalV::ofs_running << " Average iterative diagonalization steps for k-points " << ik - << " is: " << DiagoIterAssist::avg_iter - << "; where current threshold is: " << this->diag_thr << ". " << std::endl; + << " = " << DiagoIterAssist::avg_iter + << "; current threshold of diagonalization = " << this->diag_thr << ". " << std::endl; DiagoIterAssist::avg_iter = 0.0; } } @@ -526,9 +526,9 @@ void HSolverPW::output_iterInfo() // in PW base, average iteration steps for each band and k-point should be printing if (DiagoIterAssist::avg_iter > 0.0) { - GlobalV::ofs_running << " Average iterative diagonalization steps for k-points: " + GlobalV::ofs_running << " Average iterative diagonalization steps for k-points = " << DiagoIterAssist::avg_iter / this->wfc_basis->nks - << ";\n where current threshold is: " << this->diag_thr << ". " << std::endl; + << ";\n current threshold of diagonalizaiton = " << this->diag_thr << ". " << std::endl; // reset avg_iter DiagoIterAssist::avg_iter = 0.0; } From 5c8e95c20eb5acbeb573b72b31be887c37123320 Mon Sep 17 00:00:00 2001 From: Chen Nuo <49788094+Cstandardlib@users.noreply.github.com> Date: Fri, 19 Dec 2025 11:31:49 +0800 Subject: [PATCH 8/8] Make clear output format --- source/source_hsolver/hsolver_pw.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/source/source_hsolver/hsolver_pw.cpp b/source/source_hsolver/hsolver_pw.cpp index 43f027ed40..7d178f699b 100644 --- a/source/source_hsolver/hsolver_pw.cpp +++ b/source/source_hsolver/hsolver_pw.cpp @@ -143,8 +143,8 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, if (skip_charge) { GlobalV::ofs_running << " Average iterative diagonalization steps for k-points " << ik - << " = " << DiagoIterAssist::avg_iter - << "; current threshold of diagonalization = " << this->diag_thr << ". " << std::endl; + << " is " << DiagoIterAssist::avg_iter + << "\n current threshold of diagonalization is " << this->diag_thr << std::endl; DiagoIterAssist::avg_iter = 0.0; } } @@ -526,9 +526,9 @@ void HSolverPW::output_iterInfo() // in PW base, average iteration steps for each band and k-point should be printing if (DiagoIterAssist::avg_iter > 0.0) { - GlobalV::ofs_running << " Average iterative diagonalization steps for k-points = " + GlobalV::ofs_running << " Average iterative diagonalization steps for k-points is " << DiagoIterAssist::avg_iter / this->wfc_basis->nks - << ";\n current threshold of diagonalizaiton = " << this->diag_thr << ". " << std::endl; + << "\n current threshold of diagonalizaiton is " << this->diag_thr << std::endl; // reset avg_iter DiagoIterAssist::avg_iter = 0.0; }