diff --git a/example/ex7.cpp b/example/ex7.cpp index ebd4380..c5981cd 100644 --- a/example/ex7.cpp +++ b/example/ex7.cpp @@ -123,6 +123,14 @@ int main(int argc, char const *argv[]) 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; + + m.assign_all(std::complex(0.0, 0.0)); + test.CLCG_Minimize(m, B, gctl::CLCG_CGS); + 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_BICGSTAB); 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 134d779..c1a1bdc 100644 --- a/lib/optimization/clcg.cpp +++ b/lib/optimization/clcg.cpp @@ -30,7 +30,7 @@ /** * Default parameter for conjugate gradient methods */ -static const gctl::clcg_para clcg_defparam = {0, 1e-8, 0}; +static const gctl::clcg_para clcg_defparam = {0, 1e-8, 0, 1e-8}; gctl::clcg_solver::clcg_solver() { @@ -259,13 +259,13 @@ void gctl::clcg_solver::CLCG_Minimize(array > &m, const arr switch (solver_id) { case CLCG_BICG: - ss << "Solver: " << std::setw(9) << "Bi-CG, Times cost: " << costime << " ms" << std::endl; break; + ss << "Solver: " << std::setw(9) << "BiCG, Times cost: " << costime << " ms" << std::endl; break; case CLCG_BICG_SYM: - ss << "Solver: " << std::setw(9) << "Bi-CG (symmetrically accelerated), Times cost: " << costime << " ms" << std::endl; break; + ss << "Solver: " << std::setw(9) << "BiCG-Sym, Times cost: " << costime << " ms" << std::endl; break; case CLCG_CGS: ss << "Solver: " << std::setw(9) << "CGS, Times cost: " << costime << " ms" << std::endl; break; case CLCG_BICGSTAB: - ss << "Solver: " << std::setw(9) << "CGS, Times cost: " << costime << " ms" << std::endl; break; + ss << "Solver: " << std::setw(9) << "BiCG-Stab, Times cost: " << costime << " ms" << std::endl; break; case CLCG_TFQMR: ss << "Solver: " << std::setw(9) << "TFQMR, Times cost: " << costime << " ms" << std::endl; break; default: @@ -462,12 +462,219 @@ void gctl::clcg_solver::clbicg_symmetric(array > &m, const void gctl::clcg_solver::clcgs(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); + uk.resize(n_size); qk.resize(n_size); wk.resize(n_size); + Ax.resize(n_size); + + CLCG_Ax(m, Ax, gctl::NoTrans, gctl::NoConj); + + std::complex one_z(1.0, 0.0); + std::complex one_one(1.0, 1.0); + vecdiff(r1k, B, Ax, one_z, one_z); + veccpy(uk, r1k, one_z); + veccpy(d1k, r1k, one_z); + + srand(time(0)); + std::complex rhok; + do + { + for (size_t i = 0; i < n_size; i++) + { + r2k[i] = random(one_z, 2.0*one_one); + } + + rhok = vecinner(r2k, r1k); + } while (std::norm(rhok) < clcg_param_.lamda); + + double r0_square, rk_square; + std::complex r0_mod, rk_mod; + + 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; + std::complex ak, rhok2, sigma, betak; + + 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); // vk = Apk + sigma = vecinner(r2k, Ax); + ak = rhok/sigma; + + vecdiff(qk, uk, Ax, one_z, ak); + vecadd(wk, uk, qk, one_z, one_z); + + CLCG_Ax(wk, Ax, gctl::NoTrans, gctl::NoConj); + + vecapp(m, wk, ak); + vecsub(r1k, Ax, ak); + + rk_mod = vecinner(r1k, r1k); + rk_square = std::norm(rk_mod); + + if (!vecvalid(m)) + { + return clcg_error_str(CLCG_NAN_VALUE, ss); + } + + rhok2 = vecinner(r2k, r1k); + betak = rhok2/rhok; + rhok = rhok2; + + vecadd(uk, r1k, qk, one_z, betak); + vecadd(d1k, qk, d1k, one_z, betak); + vecadd(d1k, uk, d1k, one_z, betak); + } + + return clcg_error_str(CLCG_UNKNOWN_ERROR, ss); } void gctl::clcg_solver::clbicgstab(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); Ad.resize(n_size); + + CLCG_Ax(m, Ax, gctl::NoTrans, gctl::NoConj); + + std::complex one_z(1.0, 0.0); + std::complex one_one(1.0, 1.0); + vecdiff(r1k, B, Ax, one_z, one_z); + veccpy(d1k, r1k, one_z); + + srand(time(0)); + std::complex rhok; + do + { + for (size_t i = 0; i < n_size; i++) + { + r2k[i] = random(one_z, 2.0*one_one); + } + + rhok = vecinner(r2k, r1k); + } while (std::norm(rhok) < clcg_param_.lamda); + + double r0_square, rk_square; + std::complex r0_mod, rk_mod; + 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; + std::complex ak, rhok2, sigma, omega, betak, Ass, AsAs; + + 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); + sigma = vecinner(r2k, Ax); + ak = rhok/sigma; + + vecdiff(d2k, r1k, Ax, one_z, ak); + + CLCG_Ax(d2k, Ad, gctl::NoTrans, gctl::NoConj); + Ass = vecinner(Ad, d2k); + AsAs = vecinner(Ad, Ad); + omega = Ass/AsAs; + + vecdiff(r1k, d2k, Ad, one_z, omega); + vecadd(d2k, d1k, d2k, ak, omega); + vecapp(m, d2k, one_z); + + rk_mod = vecinner(r1k, r1k); + rk_square = std::norm(rk_mod); + + if (!vecvalid(m)) + { + return clcg_error_str(CLCG_NAN_VALUE, ss); + } + + rhok2 = vecinner(r2k, r1k); + betak = rhok2*ak/(rhok*omega); + rhok = rhok2; + + vecdiff(d1k, d1k, Ax, one_z, omega); + vecadd(d1k, r1k, d1k, one_z, betak); + } + + return clcg_error_str(CLCG_UNKNOWN_ERROR, ss); } void gctl::clcg_solver::cltfqmr(array > &m, const array > &B, std::ostream &ss) diff --git a/lib/optimization/clcg.h b/lib/optimization/clcg.h index 229c4c9..87c4c87 100644 --- a/lib/optimization/clcg.h +++ b/lib/optimization/clcg.h @@ -132,6 +132,13 @@ namespace gctl * applied to the non-constrained methods. */ int abs_diff; + + /** + * Minimal value for testing if a vector's module is bigger than the threshold. + * The default value is 1e-8 + * + */ + double lamda; }; class clcg_solver @@ -143,7 +150,8 @@ namespace gctl // make them class variables are more suitable for repetitively usages array > r1k, r2k, d1k, d2k; - array > Ax; + array > uk, qk, wk; + array > Ax, Ad; void clcg_error_str(clcg_return_code err_code, std::ostream &ss);