From 716de160db7099c38105b287e922fcbfd5ae27f1 Mon Sep 17 00:00:00 2001 From: "Yuhsiang M. Tsai" Date: Mon, 7 Nov 2022 21:19:32 +0800 Subject: [PATCH] use cycle_info and clean the example --- core/solver/multigrid.cpp | 147 +++++++++++------- .../mixed-multigrid-solver.cpp | 16 +- .../mixed-precision-ir/mixed-precision-ir.cpp | 4 +- .../multigrid-preconditioned-solver.cpp | 20 ++- 4 files changed, 109 insertions(+), 78 deletions(-) diff --git a/core/solver/multigrid.cpp b/core/solver/multigrid.cpp index e843bc393c1..7cf5b92ea0a 100644 --- a/core/solver/multigrid.cpp +++ b/core/solver/multigrid.cpp @@ -172,6 +172,42 @@ void clear_and_reserve(Vec& vec, size_type size) namespace multigrid { + + +/** + * The enum class is to combine the cycle infomation It's legal to use a binary + * or(|) operation to combine several properties. + */ +enum class cycle_info { + /** + * indicate input is zero + */ + x_is_zero = 1, + + /** + * current process is the first one of the cycle + */ + first_of_cycle = 2, + + /** + * current procees is the end one of the cycle + */ + end_of_cycle = 4 +}; + + +GKO_ATTRIBUTES GKO_INLINE cycle_info operator|(cycle_info a, cycle_info b) +{ + return static_cast(static_cast(a) | static_cast(b)); +} + + +GKO_ATTRIBUTES GKO_INLINE bool operator&(cycle_info a, cycle_info b) +{ + return static_cast(static_cast(a) & static_cast(b)); +} + + /** * MultigridState is to store the necessary cache and run the operation of all * levels. @@ -213,11 +249,13 @@ struct MultigridState { * @param matrix the system matrix of current level * @param b the right hand side * @param x the input vectors + * @param info the info of cycle */ void run_cycle(multigrid::cycle cycle, size_type level, const std::shared_ptr& matrix, const LinOp* b, - LinOp* x, bool x_is_zero = false, bool is_first = true, - bool is_end = true); + LinOp* x, + cycle_info info = cycle_info::first_of_cycle | + cycle_info::end_of_cycle); /** * @copydoc run_cycle @@ -229,7 +267,7 @@ struct MultigridState { template void run_cycle(multigrid::cycle cycle, size_type level, const std::shared_ptr& matrix, const LinOp* b, - LinOp* x, bool x_is_zero, bool is_first, bool is_end); + LinOp* x, cycle_info info); // current level's nrows x nrhs std::vector> r_list; @@ -318,8 +356,7 @@ void MultigridState::allocate_memory(int level, multigrid::cycle cycle, void MultigridState::run_cycle(multigrid::cycle cycle, size_type level, const std::shared_ptr& matrix, - const LinOp* b, LinOp* x, bool x_is_zero, - bool is_first, bool is_end) + const LinOp* b, LinOp* x, cycle_info info) { if (level == multigrid->get_mg_level_list().size()) { multigrid->get_coarsest_solver()->apply(b, x); @@ -331,8 +368,7 @@ void MultigridState::run_cycle(multigrid::cycle cycle, size_type level, mg_level, [&, this](auto mg_level) { using value_type = typename std::decay_t< gko::detail::pointee>::value_type; - this->run_cycle(cycle, level, matrix, b, x, x_is_zero, - is_first, is_end); + this->run_cycle(cycle, level, matrix, b, x, info); }); } @@ -340,16 +376,13 @@ void MultigridState::run_cycle(multigrid::cycle cycle, size_type level, template void MultigridState::run_cycle(multigrid::cycle cycle, size_type level, const std::shared_ptr& matrix, - const LinOp* b, LinOp* x, bool x_is_zero, - bool is_first, bool is_end) + const LinOp* b, LinOp* x, cycle_info info) { auto total_level = multigrid->get_mg_level_list().size(); auto r = r_list.at(level); auto g = g_list.at(level); auto e = e_list.at(level); - LinOp* x_ptr = x; - const LinOp* b_ptr = b; // get mg_level auto mg_level = multigrid->get_mg_level_list().at(level); // get the pre_smoother @@ -366,37 +399,36 @@ void MultigridState::run_cycle(multigrid::cycle cycle, size_type level, auto next_one = next_one_list.at(level).get(); auto neg_one = neg_one_list.at(level).get(); // origin or next or first - bool use_pre = is_first || mid_case == multigrid::mid_smooth_type::both || + bool use_pre = (info & cycle_info::first_of_cycle) || + mid_case == multigrid::mid_smooth_type::both || mid_case == multigrid::mid_smooth_type::pre_smoother; if (use_pre && pre_smoother) { - if (x_is_zero) { - // when level is zero, the x_ptr is already filled by zero - if (level != 0) { - dynamic_cast*>(x_ptr)->fill( - zero()); - } + if (info & cycle_info::x_is_zero) { if (auto pre_allow_zero_input = std::dynamic_pointer_cast( pre_smoother)) { pre_allow_zero_input->apply_with_initial_guess( - b_ptr, x_ptr, initial_guess_mode::zero); + b, x, initial_guess_mode::zero); } else { - pre_smoother->apply(b_ptr, x_ptr); + // x in first level is already filled by zero outside. + if (level != 0) { + dynamic_cast*>(x)->fill( + zero()); + } + pre_smoother->apply(b, x); } } else { - pre_smoother->apply(b_ptr, x_ptr); + pre_smoother->apply(b, x); } - // split the check - // Thus, when the IR only contains iter limit, there's no additional - // residual computation - r->copy_from(b); // n * b - matrix->apply(neg_one, x_ptr, one, r.get()); - } else if (level != 0) { - // move the residual computation at level 0 to out-of-cycle if there - // is no pre-smoother at level 0 - r->copy_from(b); - matrix->apply(neg_one, x_ptr, one, r.get()); } + // The common smoother is wrapped by IR and IR already split the iter and + // residual check. Thus, when the IR only contains iter limit, there's no + // additional residual computation + // TODO: if already computes the residual outside, the first level may not + // need this residual computation when no presmoother in the first level. + r->copy_from(b); // n * b + matrix->apply(neg_one, x, one, r.get()); + // first cycle mg_level->get_restrict_op()->apply(r.get(), g.get()); // next level @@ -408,38 +440,45 @@ void MultigridState::run_cycle(multigrid::cycle cycle, size_type level, (level + 1 < total_level) ? multigrid->get_mg_level_list().at(level + 1)->get_fine_op() : mg_level->get_coarse_op(); - this->run_cycle(cycle, level + 1, next_level_matrix, g.get(), e.get(), true, - true, cycle == multigrid::cycle::v); + auto next_info = cycle_info::x_is_zero | cycle_info::first_of_cycle; + if (cycle == multigrid::cycle::v) { + // v cycle only contains one step + next_info = next_info | cycle_info::end_of_cycle; + } + this->run_cycle(cycle, level + 1, next_level_matrix, g.get(), e.get(), + next_info); if (level < multigrid->get_mg_level_list().size() - 1) { // additional work for non-v_cycle // next level if (cycle == multigrid::cycle::f) { // f_cycle call v_cycle in the second cycle this->run_cycle(multigrid::cycle::v, level + 1, next_level_matrix, - g.get(), e.get(), false, false, true); + g.get(), e.get(), cycle_info::end_of_cycle); } else if (cycle == multigrid::cycle::w) { this->run_cycle(cycle, level + 1, next_level_matrix, g.get(), - e.get(), false, false, true); + e.get(), cycle_info::end_of_cycle); } } // prolong - mg_level->get_prolong_op()->apply(next_one, e.get(), next_one, x_ptr); + mg_level->get_prolong_op()->apply(next_one, e.get(), next_one, x); // end or origin previous - bool use_post = is_end || mid_case == multigrid::mid_smooth_type::both || + bool use_post = (info & cycle_info::end_of_cycle) || + mid_case == multigrid::mid_smooth_type::both || mid_case == multigrid::mid_smooth_type::post_smoother; // post-smooth if (use_post && post_smoother) { - post_smoother->apply(b_ptr, x_ptr); + post_smoother->apply(b, x); } // put the mid smoother into the end of previous cycle // only W/F cycle bool use_mid = (cycle == multigrid::cycle::w || cycle == multigrid::cycle::f) && - !is_end && mid_case == multigrid::mid_smooth_type::standalone; + !(info & cycle_info::end_of_cycle) && + mid_case == multigrid::mid_smooth_type::standalone; if (use_mid && mid_smoother) { - mid_smoother->apply(b_ptr, x_ptr); + mid_smoother->apply(b, x); } } @@ -617,44 +656,42 @@ void Multigrid::apply_dense_impl(const VectorType* b, VectorType* x, using value_type = typename std::decay_t< gko::detail::pointee>::value_type; auto exec = this->get_executor(); - static auto neg_one_op = - initialize>({-one()}, exec); - static auto one_op = - initialize>({one()}, exec); + auto neg_one_op = cache_.state->neg_one_list.at(0); + auto one_op = cache_.state->one_list.at(0); constexpr uint8 RelativeStoppingId{1}; auto& stop_status = this->template create_workspace_array( ws::stop, b->get_size()[1]); bool one_changed{}; exec->run(multigrid::make_initialize(&stop_status)); - // compute the residual at the r_list(0); - // auto r = state->r_list.at(0); - // r->copy_from(b) - // system_matrix->apply(lend(neg_one_op), x, lend(one_op), r.get()); auto stop_criterion = this->get_stop_criterion_factory()->generate( this->get_system_matrix(), std::shared_ptr(b, null_deleter{}), x, nullptr); int iter = -1; + while (true) { ++iter; this->template log(this, iter, nullptr, x); if (stop_criterion->update() .num_iterations(iter) - // .residual(r.get()) + // TODO: combine the out-of-cycle residual computation + // currently, the residual will computed additionally in + // stop_criterion when users require the corresponding + // residual check. .solution(x) .check(RelativeStoppingId, true, &stop_status, &one_changed)) { break; } - + auto info = multigrid::cycle_info::first_of_cycle | + multigrid::cycle_info::end_of_cycle; + if (iter == 0 && guess == initial_guess_mode::zero) { + info = info | multigrid::cycle_info::x_is_zero; + } cache_.state->run_cycle(this->get_parameters().cycle, 0, - this->get_system_matrix(), b, x, - guess == initial_guess_mode::zero); - // r->copy_from(b); - // system_matrix->apply(lend(neg_one_op), x, lend(one_op), - // r.get()); + this->get_system_matrix(), b, x, info); } }; diff --git a/examples/mixed-multigrid-solver/mixed-multigrid-solver.cpp b/examples/mixed-multigrid-solver/mixed-multigrid-solver.cpp index 3b5f67f950e..13cd230e44c 100644 --- a/examples/mixed-multigrid-solver/mixed-multigrid-solver.cpp +++ b/examples/mixed-multigrid-solver/mixed-multigrid-solver.cpp @@ -86,16 +86,15 @@ int main(int argc, char* argv[]) // executor where Ginkgo will perform the computation const auto exec = exec_map.at(executor_string)(); // throws if not valid - const bool use_mixed = argc >= 5; + const int mixed_int = argc >= 3 ? std::atoi(argv[2]) : 1; + const bool use_mixed = mixed_int != 0; // nonzero uses mixed std::cout << "Using mixed precision? " << use_mixed << std::endl; // Read data - auto A = share(gko::read(std::ifstream(argv[2]), exec)); + auto A = share(gko::read(std::ifstream("data/A.mtx"), exec)); // Create RHS as 1 and initial guess as 0 gko::size_type size = A->get_size()[0]; auto host_x = vec::create(exec->get_master(), gko::dim<2>(size, 1)); auto host_b = vec::create(exec->get_master(), gko::dim<2>(size, 1)); - // auto host_b = - // share(gko::read(std::ifstream(argv[3]), exec->get_master())); for (auto i = 0; i < size; i++) { host_x->at(i, 0) = 0.; host_b->at(i, 0) = 1.; @@ -227,15 +226,6 @@ int main(int argc, char* argv[]) std::cout << "Final residual norm sqrt(r^T r): \n"; write(std::cout, lend(res)); - // auto mg_level_list = solver->get_mg_level_list(); - // auto smoother_list = solver->get_pre_smoother_list(); - // // Check the MultigridLevel and smoother. - // // throw error if there is mismatch - // auto level0 = gko::as(mg_level_list.at(0)); - // auto level1 = gko::as(mg_level_list.at(1)); - // auto smoother0 = gko::as(smoother_list.at(0)); - // auto smoother1 = gko::as(smoother_list.at(1)); - // Print solver statistics std::cout << "Multigrid iteration count: " << logger->get_num_iterations() << std::endl; diff --git a/examples/mixed-precision-ir/mixed-precision-ir.cpp b/examples/mixed-precision-ir/mixed-precision-ir.cpp index a17d8d161e0..07a1118a065 100644 --- a/examples/mixed-precision-ir/mixed-precision-ir.cpp +++ b/examples/mixed-precision-ir/mixed-precision-ir.cpp @@ -86,8 +86,8 @@ int main(int argc, char* argv[]) }}, {"dpcpp", [] { - return gko::DpcppExecutor::create(0, - gko::OmpExecutor::create()); + return gko::DpcppExecutor::create( + 0, gko::ReferenceExecutor::create()); }}, {"reference", [] { return gko::ReferenceExecutor::create(); }}}; diff --git a/examples/multigrid-preconditioned-solver/multigrid-preconditioned-solver.cpp b/examples/multigrid-preconditioned-solver/multigrid-preconditioned-solver.cpp index 1370b0f1659..9a6a230bdd1 100644 --- a/examples/multigrid-preconditioned-solver/multigrid-preconditioned-solver.cpp +++ b/examples/multigrid-preconditioned-solver/multigrid-preconditioned-solver.cpp @@ -89,20 +89,18 @@ int main(int argc, char* argv[]) // executor where Ginkgo will perform the computation const auto exec = exec_map.at(executor_string)(); // throws if not valid - // const auto fileA = argc >= 3 ? argv[2] : "data/A.mtx"; - const bool use_mixed = argc >= 4; - std::cout << "Use Mixed: " << use_mixed << std::endl; + const int mixed_int = argc >= 3 ? std::atoi(argv[2]) : 1; + const bool use_mixed = mixed_int != 0; // nonzero uses mixed + std::cout << "Using mixed precision? " << use_mixed << std::endl; // Read data - auto A = share(gko::read(std::ifstream(argv[2]), exec)); + auto A = share(gko::read(std::ifstream("data/A.mtx"), exec)); // Create RHS as 1 and initial guess as 0 gko::size_type size = A->get_size()[0]; auto host_x = vec::create(exec->get_master(), gko::dim<2>(size, 1)); - // auto host_b = vec::create(exec->get_master(), gko::dim<2>(size, 1)); - auto host_b = - share(gko::read(std::ifstream(argv[3]), exec->get_master())); + auto host_b = vec::create(exec->get_master(), gko::dim<2>(size, 1)); for (auto i = 0; i < size; i++) { host_x->at(i, 0) = 0.; - // host_b->at(i, 0) = 1.; + host_b->at(i, 0) = 1.; } auto x = vec::create(exec); auto b = vec::create(exec); @@ -184,6 +182,12 @@ int main(int argc, char* argv[]) .with_mg_level(mg_level_gen, mg_level_gen_f) .with_level_selector([](const gko::size_type level, const gko::LinOp*) -> gko::size_type { + // The first (index 0) level will use the first + // mg_level_gen, smoother_gen which are the factories with + // ValueType. The rest of levels (>= 1) will use the second + // (index 1) mg_level_gen2 and smoother_gen2 which use the + // MixedType. The rest of levels will use different type + // than the normal multigrid. return level >= 1 ? 1 : 0; }) .with_coarsest_solver(coarsest_gen_f)