From 19df99629478fe1fa0c3ff8526b129332db7481b Mon Sep 17 00:00:00 2001 From: Yi Zhang Date: Sat, 21 Sep 2024 13:01:10 +0800 Subject: [PATCH] tmp update --- example/ex7.cpp | 4 ++ lib/optimization/clcg.cpp | 129 ++++++++++++++++++++++++++++++++------ lib/optimization/lcg.cpp | 12 ++-- 3 files changed, 120 insertions(+), 25 deletions(-) diff --git a/example/ex7.cpp b/example/ex7.cpp index e34752b..ebd4380 100644 --- a/example/ex7.cpp +++ b/example/ex7.cpp @@ -119,6 +119,10 @@ int main(int argc, char const *argv[]) test.set_clcg_para(my_para); test.CLCG_Minimize(m, B, gctl::CLCG_BICG_SYM); + std::clog << "maximal difference: " << max_diff(fm, m) << std::endl; + + m.assign_all(std::complex(0.0, 0.0)); + test.CLCG_Minimize(m, B, gctl::CLCG_BICG); std::clog << "maximal difference: " << max_diff(fm, m) << std::endl; return 0; } \ No newline at end of file diff --git a/lib/optimization/clcg.cpp b/lib/optimization/clcg.cpp index 275d623..134d779 100644 --- a/lib/optimization/clcg.cpp +++ b/lib/optimization/clcg.cpp @@ -232,23 +232,23 @@ void gctl::clcg_solver::CLCG_Minimize(array > &m, const arr #ifdef GCTL_OPENMP double start = omp_get_wtime(); - if (solver_id == CLCG_BICG) return clbicg(m, B, ss); - else if (solver_id == CLCG_BICG_SYM) return clbicg_symmetric(m, B, ss); - else if (solver_id == CLCG_CGS) return clcgs(m, B, ss); - else if (solver_id == CLCG_BICGSTAB) return clbicgstab(m, B, ss); - else if (solver_id == CLCG_TFQMR) return cltfqmr(m, B, ss); - else throw std::invalid_argument("Invalid solver type. gctl::clcg_solver::Minimize(...)"); + if (solver_id == CLCG_BICG) clbicg(m, B, ss); + else if (solver_id == CLCG_BICG_SYM) clbicg_symmetric(m, B, ss); + else if (solver_id == CLCG_CGS) clcgs(m, B, ss); + else if (solver_id == CLCG_BICGSTAB) clbicgstab(m, B, ss); + else if (solver_id == CLCG_TFQMR) cltfqmr(m, B, ss); + else throw std::invalid_argument("Invalid solver type."); double end = omp_get_wtime(); double costime = 1000*(end-start); #else clock_t start = clock(); - if (solver_id == CLCG_BICG) return clbicg(m, B, ss); - else if (solver_id == CLCG_BICG_SYM) return clbicg_symmetric(m, B, ss); - else if (solver_id == CLCG_CGS) return clcgs(m, B, ss); - else if (solver_id == CLCG_BICGSTAB) return clbicgstab(m, B, ss); - else if (solver_id == CLCG_TFQMR) return cltfqmr(m, B, ss); - else throw std::invalid_argument("Invalid solver type. gctl::clcg_solver::Minimize(...)"); + if (solver_id == CLCG_BICG) clbicg(m, B, ss); + else if (solver_id == CLCG_BICG_SYM) clbicg_symmetric(m, B, ss); + else if (solver_id == CLCG_CGS) clcgs(m, B, ss); + else if (solver_id == CLCG_BICGSTAB) clbicgstab(m, B, ss); + else if (solver_id == CLCG_TFQMR) cltfqmr(m, B, ss); + else throw std::invalid_argument("Invalid solver type."); clock_t end = clock(); double costime = 1000*(end-start)/(double)CLOCKS_PER_SEC; @@ -259,17 +259,17 @@ void gctl::clcg_solver::CLCG_Minimize(array > &m, const arr switch (solver_id) { case CLCG_BICG: - std::clog << "Solver: " << std::setw(9) << "Bi-CG, Times cost: " << costime << " ms" << std::endl; break; + ss << "Solver: " << std::setw(9) << "Bi-CG, Times cost: " << costime << " ms" << std::endl; break; case CLCG_BICG_SYM: - std::clog << "Solver: " << std::setw(9) << "Bi-CG (symmetrically accelerated), Times cost: " << costime << " ms" << std::endl; break; + ss << "Solver: " << std::setw(9) << "Bi-CG (symmetrically accelerated), Times cost: " << costime << " ms" << std::endl; break; case CLCG_CGS: - std::clog << "Solver: " << std::setw(9) << "CGS, Times cost: " << costime << " ms" << std::endl; break; + ss << "Solver: " << std::setw(9) << "CGS, Times cost: " << costime << " ms" << std::endl; break; case CLCG_BICGSTAB: - std::clog << "Solver: " << std::setw(9) << "CGS, Times cost: " << costime << " ms" << std::endl; break; + ss << "Solver: " << std::setw(9) << "CGS, Times cost: " << costime << " ms" << std::endl; break; case CLCG_TFQMR: - std::clog << "Solver: " << std::setw(9) << "TFQMR, Times cost: " << costime << " ms" << std::endl; break; + ss << "Solver: " << std::setw(9) << "TFQMR, Times cost: " << costime << " ms" << std::endl; break; default: - std::clog << "Solver: " << std::setw(9) << "Unknown, Times cost: " << costime << " ms" << std::endl; break; + ss << "Solver: " << std::setw(9) << "Unknown, Times cost: " << costime << " ms" << std::endl; break; } } return; @@ -277,7 +277,98 @@ void gctl::clcg_solver::CLCG_Minimize(array > &m, const arr void gctl::clcg_solver::clbicg(array > &m, const array > &B, std::ostream &ss) { - return; + size_t n_size = B.size(); + //check parameters + if (n_size <= 0) return clcg_error_str(CLCG_INVILAD_VARIABLE_SIZE, ss); + if (clcg_param_.max_iterations < 0) return clcg_error_str(CLCG_INVILAD_MAX_ITERATIONS, ss); + if (clcg_param_.epsilon <= 0.0 || clcg_param_.epsilon >= 1.0) return clcg_error_str(CLCG_INVILAD_EPSILON, ss); + + r1k.resize(n_size); r2k.resize(n_size); + d1k.resize(n_size); d2k.resize(n_size); + Ax.resize(n_size); + + std::complex ak, Ad1d2, r1r2_next, betak; + + CLCG_Ax(m, Ax, gctl::NoTrans, gctl::NoConj); + + std::complex one_z(1.0, 0.0); + vecdiff(r1k, B, Ax, one_z, one_z); + veccpy(d1k, r1k, one_z); + vecconj(r2k, r1k); + veccpy(d2k, r2k, one_z); + + std::complex r1r2 = vecinner(r2k, r1k); + + double r0_square, rk_square; + std::complex r0_mod; + std::complex rk_mod = vecinner(r1k, r1k); + + r0_square = rk_square = std::norm(rk_mod); + if (r0_square < 1.0) r0_square = 1.0; + + if (clcg_param_.abs_diff && sqrt(rk_square)/n_size <= clcg_param_.epsilon) + { + CLCG_Progress(m, sqrt(rk_square)/n_size, clcg_param_, 0, ss); + return clcg_error_str(CLCG_ALREADY_OPTIMIZIED, ss); + } + else if (rk_square/r0_square <= clcg_param_.epsilon) + { + CLCG_Progress(m, rk_square/r0_square, clcg_param_, 0, ss); + return clcg_error_str(CLCG_ALREADY_OPTIMIZIED, ss); + } + + double residual; + size_t t = 0; + while(1) + { + if (clcg_param_.abs_diff) residual = sqrt(rk_square)/n_size; + else residual = rk_square/r0_square; + + if (CLCG_Progress(m, residual, clcg_param_, t, ss)) + { + return clcg_error_str(CLCG_STOP, ss); + } + + if (residual <= clcg_param_.epsilon) + { + return clcg_error_str(CLCG_CONVERGENCE, ss); + } + + if (clcg_param_.max_iterations > 0 && t+1 > clcg_param_.max_iterations) + { + return clcg_error_str(CLCG_REACHED_MAX_ITERATIONS, ss); + } + + t++; + + CLCG_Ax(d1k, Ax, gctl::NoTrans, gctl::NoConj); + Ad1d2 = vecinner(d2k, Ax); + ak = r1r2/Ad1d2; + + vecapp(m, d1k, ak); + vecsub(r1k, Ax, ak); + + if (!vecvalid(m)) + { + return clcg_error_str(CLCG_NAN_VALUE, ss); + } + + rk_mod = vecinner(r1k, r1k); + rk_square = std::norm(rk_mod); + + CLCG_Ax(d2k, Ax, gctl::Trans, gctl::Conj); + + vecsub(r2k, Ax, std::conj(ak)); + + r1r2_next = vecinner(r2k, r1k); + betak = r1r2_next/r1r2; + r1r2 = r1r2_next; + + vecadd(d1k, d1k, r1k, betak, one_z); + vecadd(d2k, d2k, r2k, std::conj(betak), one_z); + } + + return clcg_error_str(CLCG_UNKNOWN_ERROR, ss); } void gctl::clcg_solver::clbicg_symmetric(array > &m, const array > &B, std::ostream &ss) diff --git a/lib/optimization/lcg.cpp b/lib/optimization/lcg.cpp index 697564b..990909a 100644 --- a/lib/optimization/lcg.cpp +++ b/lib/optimization/lcg.cpp @@ -280,17 +280,17 @@ void gctl::lcg_solver::LCG_Minimize(array &m, const array &b, lc switch (solver_id) { case LCG_CG: - ss << "Solver: " << std::setw(9) << "CG" << ", Time cost: " << costime << " ms\n"; break; + ss << "Solver: " << std::setw(9) << "CG, Time cost: " << costime << " ms\n"; break; case LCG_CGS: - ss << "Solver: " << std::setw(9) << "CGS" << ", Time cost: " << costime << " ms\n"; break; + ss << "Solver: " << std::setw(9) << "CGS, Time cost: " << costime << " ms\n"; break; case LCG_PCG: - ss << "Solver: " << std::setw(9) << "PCG" << ", Time cost: " << costime << " ms\n"; break; + ss << "Solver: " << std::setw(9) << "PCG, Time cost: " << costime << " ms\n"; break; case LCG_BICGSTAB: - ss << "Solver: " << std::setw(9) << "BICGSTAB" << ", Times cost: " << costime << " ms\n"; break; + ss << "Solver: " << std::setw(9) << "BICGSTAB, Times cost: " << costime << " ms\n"; break; case LCG_BICGSTAB2: - ss << "Solver: " << std::setw(9) << "BICGSTAB2" << ", Time cost: " << costime << " ms\n"; break; + ss << "Solver: " << std::setw(9) << "BICGSTAB2, Time cost: " << costime << " ms\n"; break; default: - ss << "Solver: " << std::setw(9) << "Unknown" << ", Time cost: " << costime << " ms\n"; break; + ss << "Solver: " << std::setw(9) << "Unknown, Time cost: " << costime << " ms\n"; break; } } return;