From 2e11533444f8860f361301da7fef0a6add61590e Mon Sep 17 00:00:00 2001 From: dyzheng Date: Sat, 30 May 2026 22:37:15 +0800 Subject: [PATCH] add generalized eigenvalue test for diago_cg_float_test --- .../test/diago_cg_float_test.cpp | 57 ++++++++++++++++++- source/source_hsolver/test/diago_mock.h | 42 ++++++++++---- 2 files changed, 85 insertions(+), 14 deletions(-) diff --git a/source/source_hsolver/test/diago_cg_float_test.cpp b/source/source_hsolver/test/diago_cg_float_test.cpp index 32514c92e31..2380431af9b 100644 --- a/source/source_hsolver/test/diago_cg_float_test.cpp +++ b/source/source_hsolver/test/diago_cg_float_test.cpp @@ -58,6 +58,33 @@ void lapackEigen(int &npw, std::vector> &hm, float *e, bool delete[] work2; } +// LAPACK reference for generalized eigenproblem when S is diagonal (complex) +void lapackGeneralEigen(int &npw, std::vector> &hm, const std::vector> &sdiag, float *e, bool outtime = false) +{ + // build transformed matrix tmp = S^{-1/2} H S^{-1/2} + std::vector> tmp(npw * npw); + for (int i = 0; i < npw; ++i) { + std::complex si = std::sqrt(sdiag[i]); + for (int j = 0; j < npw; ++j) { + std::complex sj = std::sqrt(sdiag[j]); + tmp[i * npw + j] = hm[i * npw + j] / (si * sj); + } + } + + // call cheev_ on transformed matrix + clock_t start = clock(), end; + int lwork = 2 * npw; + std::complex *work2 = new std::complex[lwork]; + float *rwork = new float[3 * npw - 2]; + int info = 0; + char tmp_c1 = 'V', tmp_c2 = 'U'; + cheev_(&tmp_c1, &tmp_c2, &npw, tmp.data(), &npw, e, work2, &lwork, rwork, &info); + end = clock(); + if (outtime) std::cout << "Lapack General Run time: " << (float)(end - start) / CLOCKS_PER_SEC << " S" << std::endl; + delete[] rwork; + delete[] work2; +} + class DiagoCGPrepare { public: @@ -85,8 +112,14 @@ class DiagoCGPrepare float *e_lapack = new float[npw]; auto ev = DIAGOTEST::hmatrix_f; - if(mypnum == 0) { lapackEigen(npw, ev, e_lapack, false); -} + if(mypnum == 0) { + if (DIAGOTEST::sdiag_f.empty()) { + lapackEigen(npw, ev, e_lapack); + } else { + auto hm_copy = ev; // operate on a copy + lapackGeneralEigen(npw, hm_copy, DIAGOTEST::sdiag_f, e_lapack); + } + } // initial guess of psi by perturbing lapack psi std::vector> psiguess(nband * npw); @@ -228,12 +261,32 @@ TEST_P(DiagoCGFloatTest, RandomHamilt) //std::cout<<"eps "<>::PW_DIAG_THR<> hpsi(dcp.nband, dcp.npw, dcp.sparsity); DIAGOTEST::hmatrix_f = hpsi.hamilt(); + DIAGOTEST::sdiag_f.clear(); // ensure sdiag is empty for this test DIAGOTEST::npw = dcp.npw; // ModuleBase::ComplexMatrix psi = hpsi.psi(); dcp.CompareEigen(hpsi.precond()); } +TEST_P(DiagoCGFloatTest, RandomHamiltAndS) +{ + DiagoCGPrepare dcp = GetParam(); + hsolver::DiagoIterAssist>::PW_DIAG_NMAX = dcp.maxiter; + hsolver::DiagoIterAssist>::PW_DIAG_THR = dcp.eps; + HPsi> hpsi(dcp.nband, dcp.npw, dcp.sparsity); + DIAGOTEST::hmatrix_f = hpsi.hamilt(); + + DIAGOTEST::npw = dcp.npw; + + // 生成正定对角 S（范围 0.5..1.5） +DIAGOTEST::sdiag_f.resize(dcp.npw); +std::default_random_engine eng(123); +std::uniform_real_distribution ud(0.5f, 1.5f); +for (int i = 0; i < dcp.npw; ++i) DIAGOTEST::sdiag_f[i] = std::complex(ud(eng), 0.0f); + + dcp.CompareEigen(hpsi.precond()); +} + INSTANTIATE_TEST_SUITE_P(VerifyCG, DiagoCGFloatTest, ::testing::Values( diff --git a/source/source_hsolver/test/diago_mock.h b/source/source_hsolver/test/diago_mock.h index 75cced8409a..81fa6bc62a7 100644 --- a/source/source_hsolver/test/diago_mock.h +++ b/source/source_hsolver/test/diago_mock.h @@ -11,6 +11,12 @@ namespace DIAGOTEST std::vector> hmatrix_local; std::vector> hmatrix_f; std::vector> hmatrix_local_f; + + // diagonal representation of overlap (S) for simple mock of generalized eigenproblem + // if empty, sPsi will treat S as identity + std::vector sdiag_d; // for double / complex + std::vector > sdiag_f; + std::vector> sdiag; int h_nr; int h_nc; int npw; @@ -409,11 +415,15 @@ void hamilt::HamiltPW::sPsi(const double* psi_i const int npw, const int nbands) const { - for (size_t i = 0; i < static_cast(nbands * nrow); i++) - { - spsi[i] = psi_in[i]; + if (DIAGOTEST::sdiag_d.size() < static_cast(nrow)) { + DIAGOTEST::sdiag_d.assign(nrow, 1.0); // 默认单位 S + } + for (int v = 0; v < nbands; ++v) { + for (int i = 0; i < nrow; ++i) { + size_t idx = static_cast(v) * nrow + i; + spsi[idx] = psi_in[idx] * DIAGOTEST::sdiag_d[i]; + } } - return; } template <> void hamilt::HamiltPW, base_device::DEVICE_CPU>::sPsi(const std::complex* psi_in, @@ -422,11 +432,15 @@ void hamilt::HamiltPW, base_device::DEVICE_CPU>::sPsi(const const int npw, const int nbands) const { - for (size_t i = 0; i < static_cast(nbands * nrow); i++) - { - spsi[i] = psi_in[i]; + if (DIAGOTEST::sdiag_d.size() < static_cast(nrow)) { + DIAGOTEST::sdiag_d.assign(nrow, 1.0); + } + for (int v = 0; v < nbands; ++v) { + for (int i = 0; i < nrow; ++i) { + size_t idx = static_cast(v) * nrow + i; + spsi[idx] = psi_in[idx] * DIAGOTEST::sdiag_d[i]; + } } - return; } template <> void hamilt::HamiltPW, base_device::DEVICE_CPU>::sPsi(const std::complex* psi_in, @@ -435,11 +449,15 @@ void hamilt::HamiltPW, base_device::DEVICE_CPU>::sPsi(const const int npw, const int nbands) const { - for (size_t i = 0; i < static_cast(nbands * nrow); i++) - { - spsi[i] = psi_in[i]; + if (DIAGOTEST::sdiag_f.size() < static_cast(nrow)) { + DIAGOTEST::sdiag_f.assign(nrow, 1.0f); + } + for (int v = 0; v < nbands; ++v) { + for (int i = 0; i < nrow; ++i) { + size_t idx = static_cast(v) * nrow + i; + spsi[idx] = psi_in[idx] * DIAGOTEST::sdiag_f[i]; + } } - return; } //Mock function h_psi