Skip to content

Commit

Permalink
use cycle_info and clean the example
Browse files Browse the repository at this point in the history
  • Loading branch information
@yhmtsai
yhmtsai committed Nov 7, 2022
1 parent 00a4e75 commit 716de16
Show file tree
Hide file tree
Showing 4 changed files with 109 additions and 78 deletions.
147 changes: 92 additions & 55 deletions core/solver/multigrid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<cycle_info>(static_cast<int>(a) | static_cast<int>(b));
}


GKO_ATTRIBUTES GKO_INLINE bool operator&(cycle_info a, cycle_info b)
{
return static_cast<bool>(static_cast<int>(a) & static_cast<int>(b));
}


/**
* MultigridState is to store the necessary cache and run the operation of all
* levels.
Expand Down Expand Up @@ -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<const LinOp>& 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
Expand All @@ -229,7 +267,7 @@ struct MultigridState {
template <typename ValueType>
void run_cycle(multigrid::cycle cycle, size_type level,
const std::shared_ptr<const LinOp>& 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<std::shared_ptr<LinOp>> r_list;
Expand Down Expand Up @@ -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<const LinOp>& 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);
Expand All @@ -331,25 +368,21 @@ 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<decltype(mg_level)>>::value_type;
this->run_cycle<value_type>(cycle, level, matrix, b, x, x_is_zero,
is_first, is_end);
this->run_cycle<value_type>(cycle, level, matrix, b, x, info);
});
}


template <typename ValueType>
void MultigridState::run_cycle(multigrid::cycle cycle, size_type level,
const std::shared_ptr<const LinOp>& 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
Expand All @@ -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<matrix::Dense<ValueType>*>(x_ptr)->fill(
zero<ValueType>());
}
if (info & cycle_info::x_is_zero) {
if (auto pre_allow_zero_input =
std::dynamic_pointer_cast<const ApplyWithInitialGuess>(
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<matrix::Dense<ValueType>*>(x)->fill(
zero<ValueType>());
}
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
Expand All @@ -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);
}
}

Expand Down Expand Up @@ -617,44 +656,42 @@ void Multigrid::apply_dense_impl(const VectorType* b, VectorType* x,
using value_type = typename std::decay_t<
gko::detail::pointee<decltype(mg_level)>>::value_type;
auto exec = this->get_executor();
static auto neg_one_op =
initialize<matrix::Dense<value_type>>({-one<value_type>()}, exec);
static auto one_op =
initialize<matrix::Dense<value_type>>({one<value_type>()}, 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<stopping_status>(
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<const LinOp>(b, null_deleter<const LinOp>{}), x,
nullptr);
int iter = -1;

while (true) {
++iter;
this->template log<log::Logger::iteration_complete>(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);
}
};

Expand Down
16 changes: 3 additions & 13 deletions examples/mixed-multigrid-solver/mixed-multigrid-solver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<mtx>(std::ifstream(argv[2]), exec));
auto A = share(gko::read<mtx>(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<vec>(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.;
Expand Down Expand Up @@ -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<pgm>(mg_level_list.at(0));
// auto level1 = gko::as<pgm2>(mg_level_list.at(1));
// auto smoother0 = gko::as<ir>(smoother_list.at(0));
// auto smoother1 = gko::as<ir2>(smoother_list.at(1));

// Print solver statistics
std::cout << "Multigrid iteration count: "
<< logger->get_num_iterations() << std::endl;
Expand Down
4 changes: 2 additions & 2 deletions examples/mixed-precision-ir/mixed-precision-ir.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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(); }}};

Expand Down
20 changes: 12 additions & 8 deletions examples/multigrid-preconditioned-solver/multigrid-preconditioned-solver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<mtx>(std::ifstream(argv[2]), exec));
auto A = share(gko::read<mtx>(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<vec>(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);
Expand Down Expand Up @@ -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)
Expand Down

0 comments on commit 716de16

Please sign in to comment.