From e20d5cb13cb17531620a00295d772d7527effa8c Mon Sep 17 00:00:00 2001 From: zzlinpku <2601110378@stu.pku.edu.cn> Date: Thu, 17 Sep 2026 17:15:02 +0800 Subject: [PATCH 1/2] Refactor OFDFT line-search task buffer to std::array --- source/source_esolver/esolver_of.cpp | 2 -- source/source_esolver/esolver_of.h | 4 +++- .../source_esolver/esolver_of_interface.cpp | 20 +++++++++---------- 3 files changed, 13 insertions(+), 13 deletions(-) diff --git a/source/source_esolver/esolver_of.cpp b/source/source_esolver/esolver_of.cpp index 50e94b2e558..cfdf1a56e25 100644 --- a/source/source_esolver/esolver_of.cpp +++ b/source/source_esolver/esolver_of.cpp @@ -18,7 +18,6 @@ namespace ModuleESolver ESolver_OF::ESolver_OF() { this->classname = "ESolver_OF"; - this->task_ = new char[60]; } ESolver_OF::~ESolver_OF() @@ -43,7 +42,6 @@ ESolver_OF::~ESolver_OF() delete[] this->nelec_; delete[] this->theta_; - delete[] this->task_; delete this->ptemp_rho_; delete this->kedf_manager_; diff --git a/source/source_esolver/esolver_of.h b/source/source_esolver/esolver_of.h index 940a26b7e80..f252adfa10b 100644 --- a/source/source_esolver/esolver_of.h +++ b/source/source_esolver/esolver_of.h @@ -1,6 +1,8 @@ #ifndef ESOLVER_OF_H #define ESOLVER_OF_H +#include + #include "esolver_fp.h" #include "source_base/opt_dcsrch.h" #include "source_base/opt_tn.hpp" @@ -58,7 +60,7 @@ class ESolver_OF : public ESolver_FP double** pdEdphi_ = nullptr; // dE/dphi double** pdLdphi_ = nullptr; // dL/dphi double** pphi_ = nullptr; // pphi[i] = ppsi.get_pointer(i), which will be freed in ~Psi(). - char* task_ = nullptr; // used in line search + std::array task_{}; // used in line search int tn_spin_flag_ = -1; // spin flag used in cal_potential, which will be called by opt_tn int max_dcsrch_ = 200; // max no. of line search int flag_ = -1; // flag of TN diff --git a/source/source_esolver/esolver_of_interface.cpp b/source/source_esolver/esolver_of_interface.cpp index 4240db9b94c..90c60e0809f 100644 --- a/source/source_esolver/esolver_of_interface.cpp +++ b/source/source_esolver/esolver_of_interface.cpp @@ -107,7 +107,7 @@ void ESolver_OF::get_step_length(double* dEdtheta, double** ptemp_phi, UnitCell& if (this->inp_->nspin == 1) { int numDC = 0; // iteration number of line search - strcpy(this->task_, "START"); + strcpy(this->task_.data(), "START"); while (true) { // update energy @@ -122,11 +122,11 @@ void ESolver_OF::get_step_length(double* dEdtheta, double** ptemp_phi, UnitCell& temp_energy += kinetic_energy + pseudopot_energy; // line search to update theta[0] - this->opt_dcsrch_->dcSrch(temp_energy, dEdtheta[0], this->theta_[0], this->task_); + this->opt_dcsrch_->dcSrch(temp_energy, dEdtheta[0], this->theta_[0], this->task_.data()); numDC++; // decide what to do next according to the output of line search - if (strncmp(this->task_, "FG", 2) == 0) // continue line search + if (strncmp(this->task_.data(), "FG", 2) == 0) // continue line search { // update tempPhi and tempRho for (int i = 0; i < this->pw_rho->nrxx; ++i) @@ -146,20 +146,20 @@ void ESolver_OF::get_step_length(double* dEdtheta, double** ptemp_phi, UnitCell& break; } } - else if (strncmp(this->task_, "CO", 2) == 0) // convergence achieved + else if (strncmp(this->task_.data(), "CO", 2) == 0) // convergence achieved { break; } - else if (strncmp(this->task_, "WA", 2) == 0) // warning of line search + else if (strncmp(this->task_.data(), "WA", 2) == 0) // warning of line search { - GlobalV::ofs_warning << "ESolver_OF linesearch: WARNING " << this->task_ << std::endl; - std::cout << this->task_ << std::endl; + GlobalV::ofs_warning << "ESolver_OF linesearch: WARNING " << this->task_.data() << std::endl; + std::cout << this->task_.data() << std::endl; break; } - else if (strncmp(this->task_, "ER", 2) == 0) // ERROR in line search + else if (strncmp(this->task_.data(), "ER", 2) == 0) // ERROR in line search { - GlobalV::ofs_warning << "ESolver_OF linesearch: ERROR " << this->task_ << std::endl; - std::cout << this->task_ << std::endl; + GlobalV::ofs_warning << "ESolver_OF linesearch: ERROR " << this->task_.data() << std::endl; + std::cout << this->task_.data() << std::endl; break; } } From 66a05560a7aa24b52745e3e41c1a95b3844af463 Mon Sep 17 00:00:00 2001 From: zzlinpku <2601110378@stu.pku.edu.cn> Date: Fri, 18 Sep 2026 11:49:22 +0800 Subject: [PATCH 2/2] Use std::string throughout OFDFT line-search status handling --- source/source_base/opt_dcsrch.cpp | 47 +++++++------- source/source_base/opt_dcsrch.h | 7 ++- source/source_base/test/CMakeLists.txt | 5 ++ source/source_base/test/opt_cg_test.cpp | 14 ++--- source/source_base/test/opt_tn_test.cpp | 14 ++--- source/source_base/test/test_opt_dcsrch.cpp | 61 +++++++++++++++++++ source/source_esolver/esolver_of.h | 4 +- .../source_esolver/esolver_of_interface.cpp | 30 ++++----- 8 files changed, 122 insertions(+), 60 deletions(-) create mode 100644 source/source_base/test/test_opt_dcsrch.cpp diff --git a/source/source_base/opt_dcsrch.cpp b/source/source_base/opt_dcsrch.cpp index 303ac658181..189ed554cdd 100644 --- a/source/source_base/opt_dcsrch.cpp +++ b/source/source_base/opt_dcsrch.cpp @@ -1,10 +1,10 @@ #include "opt_dcsrch.h" #include -#include // This file is translated from fortran codes dcstep.f of scipy. -// The structure and all annotation of the original file have been retained. +// The numerical algorithm and original Fortran annotations are retained; +// status storage uses std::string. // See original source at https://github.com/scipy/scipy/blob/main/scipy/optimize/minpack2/dcstep.f. // sunliang 2022-05-30 @@ -16,7 +16,7 @@ int dcsrch(double& stp, double& ftol, double& gtol, double& xtol, - char* task, + std::string& task, double& stpmin, double& stpmax, int* isave, @@ -115,7 +115,7 @@ int dcsrch(double& stp, // c is less than xtol. // c On exit xtol is unchanged. // c - // c task is a character variable of length at least 60. + // c task is a status string updated in place. // c On initial entry task must be set to 'START'. // c On exit task indicates the required action: // c @@ -184,45 +184,45 @@ int dcsrch(double& stp, double&, double&); // c Initialization block. - if (strncmp(task, "START", 5) == 0) + if (task.compare(0, 5, "START") == 0) { // c Check the input arguments for errors. if (stp < stpmin) { - strcpy(task, "ERROR: STP .LT. STPMIN"); + task = "ERROR: STP .LT. STPMIN"; } if (stp > stpmax) { - strcpy(task, "ERROR: STP .GT. STPMAX"); + task = "ERROR: STP .GT. STPMAX"; } if (g >= 0.) { - strcpy(task, "ERROR: INITIAL G .GE. ZERO"); + task = "ERROR: INITIAL G .GE. ZERO"; } if (ftol < 0.) { - strcpy(task, "ERROR: FTOL .LT. ZERO"); + task = "ERROR: FTOL .LT. ZERO"; } if (gtol < 0.) { - strcpy(task, "ERROR: GTOL .LT. ZERO"); + task = "ERROR: GTOL .LT. ZERO"; } if (xtol < 0.) { - strcpy(task, "ERROR: XTOL .LT. ZERO"); + task = "ERROR: XTOL .LT. ZERO"; } if (stpmin < 0.) { - strcpy(task, "ERROR: STPMIN .LT. ZERO"); + task = "ERROR: STPMIN .LT. ZERO"; } if (stpmax < stpmin) { - strcpy(task, "ERROR: STPMAX .LT. STPMIN"); + task = "ERROR: STPMAX .LT. STPMIN"; } // c Exit if there are errors on input. - if (strncmp(task, "ERROR", 5) == 0) + if (task.compare(0, 5, "ERROR") == 0) { return 0; } @@ -250,7 +250,7 @@ int dcsrch(double& stp, gy = ginit; stmin = zero; stmax = stp + stp * xtrapu; - strcpy(task, "FG"); + task = "FG"; goto L10; } else @@ -293,32 +293,31 @@ int dcsrch(double& stp, if (brackt && (stp <= stmin || stp >= stmax)) { - strcpy(task, "WARNING: ROUNDING ERRORS PREVENT PROGRESS"); + task = "WARNING: ROUNDING ERRORS PREVENT PROGRESS"; } if (brackt && stmax - stmin <= xtol * stmax) { - strcpy(task, "WARNING: XTOL TEST SATISFIED"); + task = "WARNING: XTOL TEST SATISFIED"; } if (stp == stpmax && f <= ftest && g <= gtest) { - strcpy(task, "WARNING: STP = STPMAX"); + task = "WARNING: STP = STPMAX"; } if (stp == stpmin && (f > ftest || g >= gtest)) { - strcpy(task, "WARNING: STP = STPMIN"); + task = "WARNING: STP = STPMIN"; } // c Test for convergence. if (f <= ftest && std::abs(g) <= gtol * (-ginit)) { - strcpy(task, "CONVERGENCE"); - // strcpy(task, "CONVERGENCE", 11); + task = "CONVERGENCE"; } // c Test for termination. - if (strncmp(task, "WARN", 4) == 0 || strncmp(task, "CONV", 4) == 0) + if (task.compare(0, 4, "WARN") == 0 || task.compare(0, 4, "CONV") == 0) { goto L10; } @@ -389,7 +388,7 @@ int dcsrch(double& stp, } // c Obtain another function and derivative. - strcpy(task, "FG"); + task = "FG"; L10: // c Save local variables. if (brackt) @@ -715,7 +714,7 @@ int dcsrch(double& stp, stp = stpf; } -void Opt_DCsrch::dcSrch(double& f, double& g, double& rstp, char* rtask) +void Opt_DCsrch::dcSrch(double& f, double& g, double& rstp, std::string& rtask) { dcsrch(rstp, f, diff --git a/source/source_base/opt_dcsrch.h b/source/source_base/opt_dcsrch.h index 6432833350c..cbfbfb4bbc3 100644 --- a/source/source_base/opt_dcsrch.h +++ b/source/source_base/opt_dcsrch.h @@ -2,6 +2,7 @@ #define OPT_DCSRCH_H #include +#include #include "constants.h" @@ -64,7 +65,7 @@ class Opt_DCsrch * @param g the derivative of the function at 0 on initial entry. * On subsequent entries g is the derivative of the function at x + stp * d. * @param rstp the optimized step length, assert the initial value is larger than zero. - * @param rtask a character variable of length at least 60. + * @param rtask the line-search status string, updated in place. * On initial entry task must be set to 'START'. * On exit task indicates the required action: * If task(1:2) = 'FG' then evaluate the function and derivative at stp and call dcsrch again. @@ -73,7 +74,7 @@ class Opt_DCsrch * The exit value of stp contains the best point found during the search. * If task(1:5) = 'ERROR' then there is an error in the input arguments. */ - void dcSrch(double& f, double& g, double& rstp, char* rtask); + void dcSrch(double& f, double& g, double& rstp, std::string& rtask); private: double ftol_ = 1e-4; // nonnegative tolerance for the sufficient decrease condition. @@ -87,4 +88,4 @@ class Opt_DCsrch }; } // namespace ModuleBase -#endif \ No newline at end of file +#endif diff --git a/source/source_base/test/CMakeLists.txt b/source/source_base/test/CMakeLists.txt index 3b8fb64deee..ce9df8138ba 100644 --- a/source/source_base/test/CMakeLists.txt +++ b/source/source_base/test/CMakeLists.txt @@ -153,6 +153,11 @@ AddTest( SOURCES opt_cg_test.cpp opt_test_tools.cpp mpi_test_main.cpp ) +AddTest( + TARGET MODULE_BASE_opt_dcsrch + SOURCES test_opt_dcsrch.cpp ../opt_dcsrch.cpp +) + AddTest( TARGET MODULE_BASE_opt_tn LIBS parameter base device diff --git a/source/source_base/test/opt_cg_test.cpp b/source/source_base/test/opt_cg_test.cpp index c2793f292d5..d6b1d38a72c 100644 --- a/source/source_base/test/opt_cg_test.cpp +++ b/source/source_base/test/opt_cg_test.cpp @@ -21,7 +21,7 @@ class CG_test : public testing::Test double residual = 10.; double tol = 1e-5; int final_iter = 0; - char *task = nullptr; + std::string task; double *Ap = nullptr; double *p = nullptr; double *x = nullptr; @@ -31,7 +31,6 @@ class CG_test : public testing::Test cg.set_para(1.); cg.allocate(tools.nx); cg.init_b(tools.le.b); - task = new char[60]; Ap = new double[tools.nx]; p = new double[tools.nx]; x = new double[tools.nx]; @@ -39,7 +38,6 @@ class CG_test : public testing::Test void TearDown() { - delete[] task; delete[] Ap; delete[] p; delete[] x; @@ -116,27 +114,27 @@ class CG_test : public testing::Test cg.next_direct(gradient, cg_label, p); for (int i = 0; i < 3; ++i) { temp_x[i] = x[i]; } - task[0] = 'S'; task[1] = 'T'; task[2] = 'A'; task[3] = 'R'; task[4] = 'T'; + task = "START"; while (true) { f = tools.func(temp_x, func_label); g = tools.dfuncdstp(temp_x, p, func_label); ds.dcSrch(f, g, step, task); - if (task[0] == 'F' && task[1] == 'G') + if (task.compare(0, 2, "FG") == 0) { for (int j = 0; j < 3; ++j) { temp_x[j] = x[j] + step * p[j]; } continue; } - else if (task[0] == 'C' && task[1] == 'O') + else if (task.compare(0, 2, "CO") == 0) { break; } - else if (task[0] == 'W' && task[1] == 'A') + else if (task.compare(0, 2, "WA") == 0) { break; } - else if (task[0] == 'E' && task[1] == 'R') + else if (task.compare(0, 2, "ER") == 0) { break; } diff --git a/source/source_base/test/opt_tn_test.cpp b/source/source_base/test/opt_tn_test.cpp index 8eaa7229bb5..df921d4a7ad 100644 --- a/source/source_base/test/opt_tn_test.cpp +++ b/source/source_base/test/opt_tn_test.cpp @@ -17,7 +17,7 @@ class TN_test : public testing::Test double tol = 1e-5; int final_iter = 0; int flag = 0; - char *task = nullptr; + std::string task; double *p = nullptr; double *x = nullptr; @@ -25,14 +25,12 @@ class TN_test : public testing::Test { tn.set_para(1.); tn.allocate(tools.nx); - task = new char[60]; p = new double[tools.nx]; x = new double[tools.nx]; } void TearDown() { - delete[] task; delete[] p; delete[] x; } @@ -78,27 +76,27 @@ class TN_test : public testing::Test } for (int i = 0; i < 3; ++i) { temp_x[i] = x[i]; } - task[0] = 'S'; task[1] = 'T'; task[2] = 'A'; task[3] = 'R'; task[4] = 'T'; + task = "START"; while (true) { f = tools.func(temp_x, func_label); g = tools.dfuncdstp(temp_x, p, func_label); ds.dcSrch(f, g, step, task); - if (task[0] == 'F' && task[1] == 'G') + if (task.compare(0, 2, "FG") == 0) { for (int j = 0; j < 3; ++j) { temp_x[j] = x[j] + step * p[j]; } continue; } - else if (task[0] == 'C' && task[1] == 'O') + else if (task.compare(0, 2, "CO") == 0) { break; } - else if (task[0] == 'W' && task[1] == 'A') + else if (task.compare(0, 2, "WA") == 0) { break; } - else if (task[0] == 'E' && task[1] == 'R') + else if (task.compare(0, 2, "ER") == 0) { break; } diff --git a/source/source_base/test/test_opt_dcsrch.cpp b/source/source_base/test/test_opt_dcsrch.cpp new file mode 100644 index 00000000000..5fc786fbba2 --- /dev/null +++ b/source/source_base/test/test_opt_dcsrch.cpp @@ -0,0 +1,61 @@ +#include "../opt_dcsrch.h" +#include "gtest/gtest.h" + +TEST(OptDCsrch, ConvergenceAndRestart) +{ + ModuleBase::Opt_DCsrch search; + search.set_paras(); + std::string task; + for (int run = 0; run < 2; ++run) + { + double f = 1.0; + double g = -2.0; + double step = 1.0; + task = "START"; + search.dcSrch(f, g, step, task); + ASSERT_EQ(task, "FG"); + EXPECT_EQ(task.size(), 2u); + f = (step - 1.0) * (step - 1.0); + g = 2.0 * (step - 1.0); + search.dcSrch(f, g, step, task); + EXPECT_EQ(task, "CONVERGENCE"); + EXPECT_EQ(task.size(), 11u); + EXPECT_DOUBLE_EQ(step, 1.0); + } +} + +TEST(OptDCsrch, ErrorReplacesStatus) +{ + ModuleBase::Opt_DCsrch search; + search.set_paras(); + double f = 1.0; + double g = 1.0; + double step = 1.0; + std::string task = "START"; + search.dcSrch(f, g, step, task); + EXPECT_EQ(task, "ERROR: INITIAL G .GE. ZERO"); + EXPECT_EQ(task.size(), std::string("ERROR: INITIAL G .GE. ZERO").size()); + + task = "START"; + g = -2.0; + search.dcSrch(f, g, step, task); + EXPECT_EQ(task, "FG"); + EXPECT_EQ(task.size(), 2u); +} + +TEST(OptDCsrch, WarningAtMaximumStep) +{ + ModuleBase::Opt_DCsrch search; + search.set_paras(1e-4, 0.2, 1e-12, 0.0, 0.1); + double f = 1.0; + double g = -2.0; + double step = 0.1; + std::string task = "START"; + search.dcSrch(f, g, step, task); + ASSERT_EQ(task, "FG"); + f = (step - 1.0) * (step - 1.0); + g = 2.0 * (step - 1.0); + search.dcSrch(f, g, step, task); + EXPECT_EQ(task, "WARNING: STP = STPMAX"); + EXPECT_DOUBLE_EQ(step, 0.1); +} diff --git a/source/source_esolver/esolver_of.h b/source/source_esolver/esolver_of.h index f252adfa10b..d57abb87379 100644 --- a/source/source_esolver/esolver_of.h +++ b/source/source_esolver/esolver_of.h @@ -1,7 +1,7 @@ #ifndef ESOLVER_OF_H #define ESOLVER_OF_H -#include +#include #include "esolver_fp.h" #include "source_base/opt_dcsrch.h" @@ -60,7 +60,7 @@ class ESolver_OF : public ESolver_FP double** pdEdphi_ = nullptr; // dE/dphi double** pdLdphi_ = nullptr; // dL/dphi double** pphi_ = nullptr; // pphi[i] = ppsi.get_pointer(i), which will be freed in ~Psi(). - std::array task_{}; // used in line search + std::string task_; // used in line search int tn_spin_flag_ = -1; // spin flag used in cal_potential, which will be called by opt_tn int max_dcsrch_ = 200; // max no. of line search int flag_ = -1; // flag of TN diff --git a/source/source_esolver/esolver_of_interface.cpp b/source/source_esolver/esolver_of_interface.cpp index 90c60e0809f..adc73e33981 100644 --- a/source/source_esolver/esolver_of_interface.cpp +++ b/source/source_esolver/esolver_of_interface.cpp @@ -107,7 +107,7 @@ void ESolver_OF::get_step_length(double* dEdtheta, double** ptemp_phi, UnitCell& if (this->inp_->nspin == 1) { int numDC = 0; // iteration number of line search - strcpy(this->task_.data(), "START"); + this->task_ = "START"; while (true) { // update energy @@ -122,11 +122,11 @@ void ESolver_OF::get_step_length(double* dEdtheta, double** ptemp_phi, UnitCell& temp_energy += kinetic_energy + pseudopot_energy; // line search to update theta[0] - this->opt_dcsrch_->dcSrch(temp_energy, dEdtheta[0], this->theta_[0], this->task_.data()); + this->opt_dcsrch_->dcSrch(temp_energy, dEdtheta[0], this->theta_[0], this->task_); numDC++; // decide what to do next according to the output of line search - if (strncmp(this->task_.data(), "FG", 2) == 0) // continue line search + if (this->task_.compare(0, 2, "FG") == 0) // continue line search { // update tempPhi and tempRho for (int i = 0; i < this->pw_rho->nrxx; ++i) @@ -146,20 +146,20 @@ void ESolver_OF::get_step_length(double* dEdtheta, double** ptemp_phi, UnitCell& break; } } - else if (strncmp(this->task_.data(), "CO", 2) == 0) // convergence achieved + else if (this->task_.compare(0, 2, "CO") == 0) // convergence achieved { break; } - else if (strncmp(this->task_.data(), "WA", 2) == 0) // warning of line search + else if (this->task_.compare(0, 2, "WA") == 0) // warning of line search { - GlobalV::ofs_warning << "ESolver_OF linesearch: WARNING " << this->task_.data() << std::endl; - std::cout << this->task_.data() << std::endl; + GlobalV::ofs_warning << "ESolver_OF linesearch: WARNING " << this->task_ << std::endl; + std::cout << this->task_ << std::endl; break; } - else if (strncmp(this->task_.data(), "ER", 2) == 0) // ERROR in line search + else if (this->task_.compare(0, 2, "ER") == 0) // ERROR in line search { - GlobalV::ofs_warning << "ESolver_OF linesearch: ERROR " << this->task_.data() << std::endl; - std::cout << this->task_.data() << std::endl; + GlobalV::ofs_warning << "ESolver_OF linesearch: ERROR " << this->task_ << std::endl; + std::cout << this->task_ << std::endl; break; } } @@ -202,12 +202,12 @@ void ESolver_OF::get_step_length(double* dEdtheta, double** ptemp_phi, UnitCell& // // line search along thetaDir to find thetaAlpha // this->opt_dcsrch_->set_paras(1e-4, 1e-2, 1e-12, 0., - // ModuleBase::PI/maxThetaDir); strcpy(this->task_, "START"); + // ModuleBase::PI/maxThetaDir); this->task_ = "START"; // numDC = 0; // while(true) // { // this->pelec->f_en.calculate_etot(this->pw_rho->nrxx, - // if (strncmp(this->task_, "FG", 2) == 0) + // if (this->task_.compare(0, 2, "FG") == 0) // { // for (int is = 0; is < this->inp_->nspin; ++is) // { @@ -233,17 +233,17 @@ void ESolver_OF::get_step_length(double* dEdtheta, double** ptemp_phi, UnitCell& // endl; break; // } // } - // else if (strncmp(this->task_, "CO", 2) == 0) + // else if (this->task_.compare(0, 2, "CO") == 0) // { // break; // } - // else if (strncmp(this->task_, "WA", 2) == 0) + // else if (this->task_.compare(0, 2, "WA") == 0) // { // GlobalV::ofs_warning << "ESolver_OF linesearch: // WARNING " << this->task_ << std::endl; cout << // this->task_ << endl; break; // } - // else if (strncmp(this->task_, "ER", 2) == 0) + // else if (this->task_.compare(0, 2, "ER") == 0) // { // GlobalV::ofs_warning << "ESolver_OF linesearch: ERROR // " << this->task_ << std::endl; cout << this->task_ <<