From 5ce812393c310059f5b9631f37b0af437714c119 Mon Sep 17 00:00:00 2001 From: Pratik Nayak Date: Fri, 12 Jan 2024 21:28:38 +0100 Subject: [PATCH] Update multigrid-precond example --- .../doc/results.dox | 14 +- .../multigrid-preconditioned-solver.cpp | 150 ++++++++++++++---- 2 files changed, 128 insertions(+), 36 deletions(-) diff --git a/examples/multigrid-preconditioned-solver/doc/results.dox b/examples/multigrid-preconditioned-solver/doc/results.dox index dccd3ccad93..9049debd1a9 100644 --- a/examples/multigrid-preconditioned-solver/doc/results.dox +++ b/examples/multigrid-preconditioned-solver/doc/results.dox @@ -3,18 +3,20 @@ This is the expected output: @code{.cpp} +Using PGM coarsening +Problem size: (8000, 8000) Initial residual norm sqrt(r^T r): %%MatrixMarket matrix array real general 1 1 -4.3589 +163.578 Final residual norm sqrt(r^T r): %%MatrixMarket matrix array real general 1 1 -1.69858e-09 -CG iteration count: 39 -CG generation time [ms]: 2.04293 -CG execution time [ms]: 22.3874 -CG execution time per iteration[ms]: 0.574036 +7.61038e-09 +CG iteration count: 1 +CG generation time [ms]: 14.4458 +CG execution time [ms]: 10.975 +CG execution time per iteration[ms]: 10.975 @endcode diff --git a/examples/multigrid-preconditioned-solver/multigrid-preconditioned-solver.cpp b/examples/multigrid-preconditioned-solver/multigrid-preconditioned-solver.cpp index 51f17b7821c..130e72569d2 100644 --- a/examples/multigrid-preconditioned-solver/multigrid-preconditioned-solver.cpp +++ b/examples/multigrid-preconditioned-solver/multigrid-preconditioned-solver.cpp @@ -2,9 +2,6 @@ // // SPDX-License-Identifier: BSD-3-Clause -#include - - #include #include #include @@ -12,6 +9,9 @@ #include +#include + + int main(int argc, char* argv[]) { // Some shortcuts @@ -20,13 +20,21 @@ int main(int argc, char* argv[]) using vec = gko::matrix::Dense; using mtx = gko::matrix::Csr; using cg = gko::solver::Cg; + using ir = gko::solver::Ir; using mg = gko::solver::Multigrid; using pgm = gko::multigrid::Pgm; + using bj = gko::preconditioner::Jacobi; + using uniform_coarsening = + gko::multigrid::UniformCoarsening; // Print version information std::cout << gko::version_info::get() << std::endl; const auto executor_string = argc >= 2 ? argv[1] : "reference"; + const auto coarse_type = argc >= 3 ? argv[2] : "pgm"; + const unsigned num_jumps = argc >= 4 ? std::atoi(argv[3]) : 2u; + const unsigned grid_dim = argc >= 5 ? std::atoi(argv[4]) : 20u; + // Figure out where to run the code std::map()>> exec_map{ @@ -49,21 +57,52 @@ 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 num_rows = grid_dim * grid_dim * grid_dim; // Read data - auto A = share(gko::read(std::ifstream("data/A.mtx"), exec)); + gko::matrix_data A_data; + gko::matrix_data b_data; + gko::matrix_data x_data; + A_data.size = {num_rows, num_rows}; + b_data.size = {num_rows, 1}; + x_data.size = {num_rows, 1}; + for (int i = 0; i < grid_dim; i++) { + for (int j = 0; j < grid_dim; j++) { + for (int k = 0; k < grid_dim; k++) { + auto idx = i * grid_dim * grid_dim + j * grid_dim + k; + if (i > 0) + A_data.nonzeros.emplace_back(idx, idx - grid_dim * grid_dim, + -1); + if (j > 0) + A_data.nonzeros.emplace_back(idx, idx - grid_dim, -1); + if (k > 0) A_data.nonzeros.emplace_back(idx, idx - 1, -1); + A_data.nonzeros.emplace_back(idx, idx, 8); + if (k < grid_dim - 1) + A_data.nonzeros.emplace_back(idx, idx + 1, -1); + if (j < grid_dim - 1) + A_data.nonzeros.emplace_back(idx, idx + grid_dim, -1); + if (i < grid_dim - 1) + A_data.nonzeros.emplace_back(idx, idx + grid_dim * grid_dim, + -1); + b_data.nonzeros.emplace_back( + idx, 0, std::sin(i * 0.01 + j * 0.14 + k * 0.056)); + x_data.nonzeros.emplace_back(idx, 0, 1.0); + } + } + } + + auto A = share(mtx::create(exec, A_data.size)); + A->read(A_data); + // auto A = share(gko::read(std::ifstream("data/A.mtx"), exec)); + + A->set_strategy(std::make_shared()); // 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)); - for (auto i = 0; i < size; i++) { - host_x->at(i, 0) = 0.; - host_b->at(i, 0) = 1.; - } auto x = vec::create(exec); auto b = vec::create(exec); - x->copy_from(host_x); - b->copy_from(host_b); + x->read(x_data); + b->read(b_data); + auto b_clone = gko::clone(b); // Calculate initial residual by overwriting b auto one = gko::initialize({1.0}, exec); @@ -73,24 +112,76 @@ int main(int argc, char* argv[]) b->compute_norm2(initres); // copy b again - b->copy_from(host_b); + b->copy_from(b_clone); - // Create multigrid factory - std::shared_ptr multigrid_gen; - multigrid_gen = - mg::build() - .with_mg_level(pgm::build().with_deterministic(true)) - .with_criteria(gko::stop::Iteration::build().with_max_iters(1u)) - .on(exec); + // Prepare the stopping criteria const gko::remove_complex tolerance = 1e-8; - auto solver_gen = - cg::build() - .with_criteria(gko::stop::Iteration::build().with_max_iters(100u), - gko::stop::ResidualNorm::build() - .with_baseline(gko::stop::mode::absolute) - .with_reduction_factor(tolerance)) - .with_preconditioner(multigrid_gen) - .on(exec); + auto iter_stop = gko::share( + gko::stop::Iteration::build().with_max_iters(1000u).on(exec)); + auto tol_stop = gko::share(gko::stop::ResidualNorm::build() + .with_reduction_factor(tolerance) + .with_baseline(gko::stop::mode::absolute) + .on(exec)); + + std::shared_ptr> logger = + gko::log::Convergence::create(); + iter_stop->add_logger(logger); + tol_stop->add_logger(logger); + + // Create smoother factory (ir with bj) + auto inner_solver_gen = + gko::share(bj::build().with_max_block_size(1u).on(exec)); + auto smoother_gen = gko::share( + ir::build() + .with_solver(inner_solver_gen) + .with_relaxation_factor(0.9) + .with_criteria( + gko::stop::Iteration::build().with_max_iters(1u).on(exec)) + .on(exec)); + // Create MultigridLevel factory + auto mg_level_gen = + gko::share(pgm::build().with_deterministic(true).on(exec)); + // Create MultigridLevel factory + auto coarse_unif_gen = gko::share( + uniform_coarsening::build().with_num_jumps(num_jumps).on(exec)); + + // Create CoarsestSolver factory + auto coarsest_gen = gko::share( + ir::build() + .with_solver(inner_solver_gen) + .with_relaxation_factor(0.9) + .with_criteria( + gko::stop::Iteration::build().with_max_iters(4u).on(exec)) + .on(exec)); + // Create multigrid factory + auto multigrid_gen = gko::share(mg::build() + .with_max_levels(10u) + .with_min_coarse_rows(10u) + .with_pre_smoother(smoother_gen) + .with_post_uses_pre(true) + .with_mg_level(mg_level_gen) + .with_coarsest_solver(coarsest_gen) + .with_criteria(iter_stop, tol_stop) + .on(exec)); + if (coarse_type == std::string("uniform")) { + std::cout << "Using Uniform coarsening" << std::endl; + multigrid_gen = gko::share(mg::build() + .with_max_levels(10u) + .with_min_coarse_rows(10u) + .with_pre_smoother(smoother_gen) + .with_post_uses_pre(true) + .with_mg_level(coarse_unif_gen) + .with_coarsest_solver(coarsest_gen) + .with_criteria(iter_stop, tol_stop) + .on(exec)); + } else { + std::cout << "Using PGM coarsening" << std::endl; + } + // Create solver factory + auto solver_gen = cg::build() + .with_criteria(iter_stop, tol_stop) + .with_preconditioner(multigrid_gen) + .on(exec); // Create solver std::chrono::nanoseconds gen_time(0); auto gen_tic = std::chrono::steady_clock::now(); @@ -101,8 +192,6 @@ int main(int argc, char* argv[]) std::chrono::duration_cast(gen_toc - gen_tic); // Add logger - std::shared_ptr> logger = - gko::log::Convergence::create(); solver->add_logger(logger); // Solve system @@ -117,6 +206,7 @@ int main(int argc, char* argv[]) // Calculate residual auto res = gko::as(logger->get_residual_norm()); + std::cout << "Problem size: " << A->get_size() << std::endl; std::cout << "Initial residual norm sqrt(r^T r): \n"; write(std::cout, initres); std::cout << "Final residual norm sqrt(r^T r): \n";