From b17e7b529e81ca091477257e6344fa63af2ab32d Mon Sep 17 00:00:00 2001 From: Yi Zhang Date: Mon, 24 Feb 2025 23:58:35 +0800 Subject: [PATCH] tmp --- example/ex1.cpp | 2 ++ lib/optimization/lcg.cpp | 29 +++++++++++++++++++++++++++++ lib/optimization/lcg.h | 12 ++++++++++-- 3 files changed, 41 insertions(+), 2 deletions(-) diff --git a/example/ex1.cpp b/example/ex1.cpp index 7fbbe06..e8d8eae 100644 --- a/example/ex1.cpp +++ b/example/ex1.cpp @@ -125,6 +125,8 @@ int main(int argc, char const *argv[]) test.LCG_Minimize(m, B, gctl::LCG_CG, ofile); ofile << "maximal difference: " << max_diff(fm, m) << std::endl; + test.save_convergence("convergence"); + m.assign(0.0); test.LCG_Minimize(m, B, gctl::LCG_PCG, ofile); diff --git a/lib/optimization/lcg.cpp b/lib/optimization/lcg.cpp index 990909a..c889b91 100644 --- a/lib/optimization/lcg.cpp +++ b/lib/optimization/lcg.cpp @@ -93,6 +93,19 @@ void gctl::lcg_solver::set_lcg_para(const lcg_para &in_param) return; } +void gctl::lcg_solver::save_convergence(std::string filename) +{ + std::ofstream ofs; + open_outfile(ofs, filename, ".txt"); + + for (size_t i = 0; i < rcd.size(); ++i) + { + ofs << i << " " << rcd[i] << std::endl; + } + ofs.close(); + return; +} + gctl::lcg_para gctl::lcg_solver::default_lcg_para() { lcg_para dp = lcg_defparam; @@ -345,6 +358,7 @@ void gctl::lcg_solver::lcg(array &m, const array &B, std::ostrea gk.resize(n_size); dk.resize(n_size); Adk.resize(n_size); + if (!rcd.empty()) rcd.clear(); LCG_Ax(m, Adk); @@ -361,6 +375,7 @@ void gctl::lcg_solver::lcg(array &m, const array &B, std::ostrea { if (lcg_param_.abs_diff) residual = sqrt(gk_mod)/n_size; else residual = gk_mod/g0_mod; + rcd.push_back(residual); if (LCG_Progress(m, residual, lcg_param_, t, ss)) { @@ -415,6 +430,7 @@ void gctl::lcg_solver::lpcg(array &m, const array &B, std::ostre // locate memory rk.resize(n_size); zk.resize(n_size); dk.resize(n_size); Adk.resize(n_size); + if (!rcd.empty()) rcd.clear(); LCG_Ax(m, Adk); @@ -436,6 +452,7 @@ void gctl::lcg_solver::lpcg(array &m, const array &B, std::ostre { if (lcg_param_.abs_diff) residual = sqrt(rk_mod)/n_size; else residual = rk_mod/r0_mod; + rcd.push_back(residual); if (LCG_Progress(m, residual, lcg_param_, t, ss)) { @@ -499,6 +516,7 @@ void gctl::lcg_solver::lcgs(array &m, const array &B, std::ostre uk.resize(n_size); qk.resize(n_size); wk.resize(n_size); + if (!rcd.empty()) rcd.clear(); LCG_Ax(m, Apx); @@ -520,6 +538,7 @@ void gctl::lcg_solver::lcgs(array &m, const array &B, std::ostre { if (lcg_param_.abs_diff) residual = sqrt(rk_mod)/n_size; else residual = rk_mod/r0_mod; + rcd.push_back(residual); if (LCG_Progress(m, residual, lcg_param_, t, ss)) { @@ -584,6 +603,7 @@ void gctl::lcg_solver::lbicgstab(array &m, const array &B, std:: pk.resize(n_size); Adk.resize(n_size); sk.resize(n_size); Apx.resize(n_size); yk.resize(n_size); + if (!rcd.empty()) rcd.clear(); LCG_Ax(m, Adk); @@ -602,6 +622,7 @@ void gctl::lcg_solver::lbicgstab(array &m, const array &B, std:: { if (lcg_param_.abs_diff) residual = sqrt(rk_mod)/n_size; else residual = rk_mod/r0_mod; + rcd.push_back(residual); if (LCG_Progress(m, residual, lcg_param_, t, ss)) { @@ -669,6 +690,7 @@ void gctl::lcg_solver::lbicgstab2(array &m, const array &B, std: pk.resize(n_size); Adk.resize(n_size); sk.resize(n_size); Apx.resize(n_size); yk.resize(n_size); + if (!rcd.empty()) rcd.clear(); LCG_Ax(m, Adk); @@ -688,6 +710,7 @@ void gctl::lcg_solver::lbicgstab2(array &m, const array &B, std: { if (lcg_param_.abs_diff) residual = sqrt(rk_mod)/n_size; else residual = rk_mod/r0_mod; + rcd.push_back(residual); if (LCG_Progress(m, residual, lcg_param_, t, ss)) { @@ -804,6 +827,8 @@ void gctl::lcg_solver::lpg(array &m, const array &B, gk_new.resize(n_size); sk.resize(n_size); yk.resize(n_size); + if (!rcd.empty()) rcd.clear(); + double alpha_k = lcg_param_.step; vectop(m, hig); @@ -823,6 +848,7 @@ void gctl::lcg_solver::lpg(array &m, const array &B, { if (lcg_param_.abs_diff) residual = sqrt(gk_mod)/n_size; else residual = gk_mod/g0_mod; + rcd.push_back(residual); if (LCG_Progress(m, residual, lcg_param_, t, ss)) { @@ -889,6 +915,8 @@ void gctl::lcg_solver::lspg(array &m, const array &B, yk.resize(n_size); dk.resize(n_size); qk_m.resize(lcg_param_.maxi_m); + if (!rcd.empty()) rcd.clear(); + double lambda_k = lcg_param_.step; double qk = 0; @@ -923,6 +951,7 @@ void gctl::lcg_solver::lspg(array &m, const array &B, { if (lcg_param_.abs_diff) residual = sqrt(gk_mod)/n_size; else residual = gk_mod/g0_mod; + rcd.push_back(residual); if (LCG_Progress(m, residual, lcg_param_, t, ss)) { diff --git a/lib/optimization/lcg.h b/lib/optimization/lcg.h index 20212f6..fe5b51e 100644 --- a/lib/optimization/lcg.h +++ b/lib/optimization/lcg.h @@ -28,7 +28,7 @@ #ifndef _GCTL_LCG_H #define _GCTL_LCG_H -#include "iostream" +#include #include "gctl_optimization_config.h" #ifdef GCTL_OPTIMIZATION_TOML @@ -40,7 +40,7 @@ #endif // _WINDOWS || __WIN32__ #include "gctl/utility/stream_t.h" -#include "gctl/core/array.h" +#include "gctl/utility/stream.h" #include "gctl/maths/linear_algebra.h" namespace gctl @@ -201,6 +201,7 @@ namespace gctl array Apx, uk, qk, qk_m, wk; array m_new, gk_new; array sk, yk; + std::vector rcd; /** * @brief Display info of a given return code. This is a private function @@ -266,6 +267,13 @@ namespace gctl * @param param Input lcg parameters. */ void set_lcg_para(const lcg_para ¶m); + + /** + * @brief Save the convergence history to a file. + * + * @param filename Filename of the output file. + */ + void save_convergence(std::string filename); /** * @brief Return a lcg_para object with default values