Ginkgo  Generated from pipelines/1589998975 branch based on develop. Ginkgo version 1.10.0
A numerical linear algebra library targeting many-core architectures
The poisson-solver program

The poisson solver example..

This example depends on simple-solver.

Table of contents
  1. Introduction
  2. The commented program
  1. Results
  2. The plain program

Introduction

This example solves a 1D Poisson equation:

$ u : [0, 1] \rightarrow R\\ u'' = f\\ u(0) = u0\\ u(1) = u1 $

using a finite difference method on an equidistant grid with K discretization points (K can be controlled with a command line parameter). The discretization is done via the second order Taylor polynomial:

$ u(x + h) = u(x) - u'(x)h + 1/2 u''(x)h^2 + O(h^3)\\ u(x - h) = u(x) + u'(x)h + 1/2 u''(x)h^2 + O(h^3) / +\\ ---------------------- \\ -u(x - h) + 2u(x) + -u(x + h) = -f(x)h^2 + O(h^3) $

For an equidistant grid with K "inner" discretization points $x1, ..., xk, $and step size $ h = 1 / (K + 1)$, the formula produces a system of linear equations

$ 2u_1 - u_2 = -f_1 h^2 + u0\\ -u_(k-1) + 2u_k - u_(k+1) = -f_k h^2, k = 2, ..., K - 1\\ -u_(K-1) + 2u_K = -f_K h^2 + u1\\ $

which is then solved using Ginkgo's implementation of the CG method preconditioned with block-Jacobi. It is also possible to specify on which executor Ginkgo will solve the system via the command line. The function $`f` $is set to $`f(x) = 6x`$ (making the solution $`u(x) = x^3`$), but that can be changed in the main function.

The intention of the example is to show how Ginkgo can be used to build an application solving a real-world problem, which includes a solution of a large, sparse linear system as a component.

About the example

The commented program

row_ptrs[i + 1] = pos;
}
}

Generates the RHS vector given f and the boundary conditions.

template <typename Closure, typename ValueType>
void generate_rhs(Closure f, ValueType u0, ValueType u1,
gko::matrix::Dense<ValueType>* rhs)
{
const auto discretization_points = rhs->get_size()[0];
auto values = rhs->get_values();
const ValueType h = 1.0 / static_cast<ValueType>(discretization_points + 1);
for (gko::size_type i = 0; i < discretization_points; ++i) {
const auto xi = static_cast<ValueType>(i + 1) * h;
values[i] = -f(xi) * h * h;
}
values[0] += u0;
values[discretization_points - 1] += u1;
}

Prints the solution u.

template <typename Closure, typename ValueType>
void print_solution(ValueType u0, ValueType u1,
const gko::matrix::Dense<ValueType>* u)
{
std::cout << u0 << '\n';
for (int i = 0; i < u->get_size()[0]; ++i) {
std::cout << u->get_const_values()[i] << '\n';
}
std::cout << u1 << std::endl;
}

Computes the 1-norm of the error given the computed u and the correct solution function correct_u.

template <typename Closure, typename ValueType>
gko::remove_complex<ValueType> calculate_error(
int discretization_points, const gko::matrix::Dense<ValueType>* u,
Closure correct_u)
{
const ValueType h = 1.0 / static_cast<ValueType>(discretization_points + 1);
gko::remove_complex<ValueType> error = 0.0;
for (int i = 0; i < discretization_points; ++i) {
using std::abs;
const auto xi = static_cast<ValueType>(i + 1) * h;
error +=
abs(u->get_const_values()[i] - correct_u(xi)) / abs(correct_u(xi));
}
return error;
}
int main(int argc, char* argv[])
{

Some shortcuts

using ValueType = double;
using IndexType = int;
using vec = gko::matrix::Dense<ValueType>;
using mtx = gko::matrix::Csr<ValueType, IndexType>;
using cg = gko::solver::Cg<ValueType>;
using bj = gko::preconditioner::Jacobi<ValueType, IndexType>;

Print version information

std::cout << gko::version_info::get() << std::endl;
if (argc == 2 && (std::string(argv[1]) == "--help")) {
std::cerr << "Usage: " << argv[0]
<< " [executor] [DISCRETIZATION_POINTS]" << std::endl;
std::exit(-1);
}

Get number of discretization points

const unsigned int discretization_points =
argc >= 3 ? std::atoi(argv[2]) : 100;

Get the executor string

const auto executor_string = argc >= 2 ? argv[1] : "reference";

Figure out where to run the code

std::map<std::string, std::function<std::shared_ptr<gko::Executor>()>>
exec_map{
{"omp", [] { return gko::OmpExecutor::create(); }},
{"cuda",
[] {
return gko::CudaExecutor::create(0,
gko::OmpExecutor::create());
}},
{"hip",
[] {
return gko::HipExecutor::create(0, gko::OmpExecutor::create());
}},
{"dpcpp",
[] {
return gko::DpcppExecutor::create(0,
gko::OmpExecutor::create());
}},
{"reference", [] { return gko::ReferenceExecutor::create(); }}};

executor where Ginkgo will perform the computation

const auto exec = exec_map.at(executor_string)(); // throws if not valid

executor used by the application

const auto app_exec = exec->get_master();

Set up the problem: define the exact solution, the right hand side and the Dirichlet boundary condition.

auto correct_u = [](ValueType x) { return x * x * x; };
auto f = [](ValueType x) { return ValueType(6) * x; };
auto u0 = correct_u(0);
auto u1 = correct_u(1);

initialize matrix and vectors

auto matrix = mtx::create(app_exec, gko::dim<2>(discretization_points),
3 * discretization_points - 2);
generate_stencil_matrix(matrix.get());
auto rhs = vec::create(app_exec, gko::dim<2>(discretization_points, 1));
generate_rhs(f, u0, u1, rhs.get());
auto u = vec::create(app_exec, gko::dim<2>(discretization_points, 1));
for (int i = 0; i < u->get_size()[0]; ++i) {
u->get_values()[i] = 0.0;
}
const gko::remove_complex<ValueType> reduction_factor = 1e-7;

Generate solver and solve the system

cg::build()
.with_criteria(
gko::stop::Iteration::build().with_max_iters(discretization_points),
gko::stop::ResidualNorm<ValueType>::build().with_reduction_factor(
reduction_factor))
.with_preconditioner(bj::build())
.on(exec)
->generate(clone(exec, matrix)) // copy the matrix to the executor
->apply(rhs, u);

Uncomment to print the solution print_solution<ValueType>(u0, u1, u.get());

std::cout << "Solve complete.\nThe average relative error is "
<< calculate_error(discretization_points, u.get(), correct_u) /
static_cast<gko::remove_complex<ValueType>>(
discretization_points)
<< std::endl;
}

Results

This is the expected output:

Solve complete.
The average relative error is 2.52236e-11

Comments about programming and debugging

The plain program

#include <iostream>
#include <map>
#include <string>
#include <vector>
#include <ginkgo/ginkgo.hpp>
template <typename ValueType, typename IndexType>
void generate_stencil_matrix(gko::matrix::Csr<ValueType, IndexType>* matrix)
{
const auto discretization_points = matrix->get_size()[0];
auto row_ptrs = matrix->get_row_ptrs();
auto col_idxs = matrix->get_col_idxs();
auto values = matrix->get_values();
int pos = 0;
const ValueType coefs[] = {-1, 2, -1};
row_ptrs[0] = pos;
for (int i = 0; i < discretization_points; ++i) {
for (auto ofs : {-1, 0, 1}) {
if (0 <= i + ofs && i + ofs < discretization_points) {
values[pos] = coefs[ofs + 1];
col_idxs[pos] = i + ofs;
++pos;
}
}
row_ptrs[i + 1] = pos;
}
}
template <typename Closure, typename ValueType>
void generate_rhs(Closure f, ValueType u0, ValueType u1,
gko::matrix::Dense<ValueType>* rhs)
{
const auto discretization_points = rhs->get_size()[0];
auto values = rhs->get_values();
const ValueType h = 1.0 / static_cast<ValueType>(discretization_points + 1);
for (gko::size_type i = 0; i < discretization_points; ++i) {
const auto xi = static_cast<ValueType>(i + 1) * h;
values[i] = -f(xi) * h * h;
}
values[0] += u0;
values[discretization_points - 1] += u1;
}
template <typename Closure, typename ValueType>
void print_solution(ValueType u0, ValueType u1,
const gko::matrix::Dense<ValueType>* u)
{
std::cout << u0 << '\n';
for (int i = 0; i < u->get_size()[0]; ++i) {
std::cout << u->get_const_values()[i] << '\n';
}
std::cout << u1 << std::endl;
}
template <typename Closure, typename ValueType>
gko::remove_complex<ValueType> calculate_error(
int discretization_points, const gko::matrix::Dense<ValueType>* u,
Closure correct_u)
{
const ValueType h = 1.0 / static_cast<ValueType>(discretization_points + 1);
gko::remove_complex<ValueType> error = 0.0;
for (int i = 0; i < discretization_points; ++i) {
using std::abs;
const auto xi = static_cast<ValueType>(i + 1) * h;
error +=
abs(u->get_const_values()[i] - correct_u(xi)) / abs(correct_u(xi));
}
return error;
}
int main(int argc, char* argv[])
{
using ValueType = double;
using IndexType = int;
using vec = gko::matrix::Dense<ValueType>;
using mtx = gko::matrix::Csr<ValueType, IndexType>;
using cg = gko::solver::Cg<ValueType>;
using bj = gko::preconditioner::Jacobi<ValueType, IndexType>;
std::cout << gko::version_info::get() << std::endl;
if (argc == 2 && (std::string(argv[1]) == "--help")) {
std::cerr << "Usage: " << argv[0]
<< " [executor] [DISCRETIZATION_POINTS]" << std::endl;
std::exit(-1);
}
const unsigned int discretization_points =
argc >= 3 ? std::atoi(argv[2]) : 100;
const auto executor_string = argc >= 2 ? argv[1] : "reference";
std::map<std::string, std::function<std::shared_ptr<gko::Executor>()>>
exec_map{
{"omp", [] { return gko::OmpExecutor::create(); }},
{"cuda",
[] {
return gko::CudaExecutor::create(0,
gko::OmpExecutor::create());
}},
{"hip",
[] {
return gko::HipExecutor::create(0, gko::OmpExecutor::create());
}},
{"dpcpp",
[] {
return gko::DpcppExecutor::create(0,
gko::OmpExecutor::create());
}},
{"reference", [] { return gko::ReferenceExecutor::create(); }}};
const auto exec = exec_map.at(executor_string)(); // throws if not valid
const auto app_exec = exec->get_master();
auto correct_u = [](ValueType x) { return x * x * x; };
auto f = [](ValueType x) { return ValueType(6) * x; };
auto u0 = correct_u(0);
auto u1 = correct_u(1);
auto matrix = mtx::create(app_exec, gko::dim<2>(discretization_points),
3 * discretization_points - 2);
generate_stencil_matrix(matrix.get());
auto rhs = vec::create(app_exec, gko::dim<2>(discretization_points, 1));
generate_rhs(f, u0, u1, rhs.get());
auto u = vec::create(app_exec, gko::dim<2>(discretization_points, 1));
for (int i = 0; i < u->get_size()[0]; ++i) {
u->get_values()[i] = 0.0;
}
const gko::remove_complex<ValueType> reduction_factor = 1e-7;
cg::build()
.with_criteria(
gko::stop::Iteration::build().with_max_iters(discretization_points),
gko::stop::ResidualNorm<ValueType>::build().with_reduction_factor(
reduction_factor))
.with_preconditioner(bj::build())
.on(exec)
->generate(clone(exec, matrix)) // copy the matrix to the executor
->apply(rhs, u);
std::cout << "Solve complete.\nThe average relative error is "
<< calculate_error(discretization_points, u.get(), correct_u) /
static_cast<gko::remove_complex<ValueType>>(
discretization_points)
<< std::endl;
}
gko::matrix::Csr::get_col_idxs
index_type * get_col_idxs() noexcept
Returns the column indexes of the matrix.
Definition: csr.hpp:886
gko::matrix::Csr
CSR is a matrix format which stores only the nonzero coefficients by compressing each row of the matr...
Definition: matrix.hpp:28
gko::matrix::Dense
Dense is a matrix format which explicitly stores all values of the matrix.
Definition: dense_cache.hpp:19
gko::abs
constexpr std::enable_if_t<!is_complex_s< T >::value, T > abs(const T &x)
Returns the absolute value of the object.
Definition: math.hpp:931
gko::size_type
std::size_t size_type
Integral type used for allocation quantities.
Definition: types.hpp:89
gko::Executor::get_master
virtual std::shared_ptr< Executor > get_master() noexcept=0
Returns the master OmpExecutor of this Executor.
gko::clone
detail::cloned_type< Pointer > clone(const Pointer &p)
Creates a unique clone of the object pointed to by p.
Definition: utils_helper.hpp:173
gko::HipExecutor::create
static std::shared_ptr< HipExecutor > create(int device_id, std::shared_ptr< Executor > master, bool device_reset, allocation_mode alloc_mode=default_hip_alloc_mode, CUstream_st *stream=nullptr)
Creates a new HipExecutor.
gko::version_info::get
static const version_info & get()
Returns an instance of version_info.
Definition: version.hpp:139
gko::preconditioner::Jacobi
A block-Jacobi preconditioner is a block-diagonal linear operator, obtained by inverting the diagonal...
Definition: jacobi.hpp:187
gko::stop::ResidualNorm
The ResidualNorm class is a stopping criterion which stops the iteration process when the actual resi...
Definition: residual_norm.hpp:113
gko::matrix::Csr::get_row_ptrs
index_type * get_row_ptrs() noexcept
Returns the row pointers of the matrix.
Definition: csr.hpp:905
gko::dim< 2 >
gko::solver::Cg
CG or the conjugate gradient method is an iterative type Krylov subspace method which is suitable for...
Definition: cg.hpp:48
gko::CudaExecutor::create
static std::shared_ptr< CudaExecutor > create(int device_id, std::shared_ptr< Executor > master, bool device_reset, allocation_mode alloc_mode=default_cuda_alloc_mode, CUstream_st *stream=nullptr)
Creates a new CudaExecutor.
gko::solver::initial_guess_mode::rhs
the input is right hand side
gko::OmpExecutor::create
static std::shared_ptr< OmpExecutor > create(std::shared_ptr< CpuAllocatorBase > alloc=std::make_shared< CpuAllocator >())
Creates a new OmpExecutor.
Definition: executor.hpp:1396
gko::matrix::Dense::get_const_values
const value_type * get_const_values() const noexcept
Returns a pointer to the array of values of the matrix.
Definition: dense.hpp:860
gko::LinOp::get_size
const dim< 2 > & get_size() const noexcept
Returns the size of the operator.
Definition: lin_op.hpp:210
gko::remove_complex
typename detail::remove_complex_s< T >::type remove_complex
Obtain the type which removed the complex of complex/scalar type or the template parameter of class b...
Definition: math.hpp:260
gko::DpcppExecutor::create
static std::shared_ptr< DpcppExecutor > create(int device_id, std::shared_ptr< Executor > master, std::string device_type="all", dpcpp_queue_property property=dpcpp_queue_property::in_order)
Creates a new DpcppExecutor.
gko::matrix::Csr::get_values
value_type * get_values() noexcept
Returns the values of the matrix.
Definition: csr.hpp:867
Generated by   doxygen 1.8.16