Simple solver#
The shortest possible end-to-end Ginkgo program: pick an executor, read a matrix and two vectors from disk, generate a CG solver with a stopping criterion, apply it, and print the residual.
Kind: basic
Upstream source: examples/simple-solver/simple-solver.cpp in the Ginkgo repository.
Introduction#
If you’ve never written a Ginkgo program, start here. The walkthrough below covers the four ingredients that every solve needs:
An
Executor: where the computation runs (CPU, GPU, …).A
LinOp: the system matrix, usuallygko::matrix::Csrfor sparse problems.A right-hand-side vector and a solution vector, both
gko::matrix::Dense.A solver factory and a stopping criterion.
Once you have those four, the apply is one line. The same pattern scales up to multigrid, distributed, batched, and everything else in the library.
The commented program#
Includes and type aliases#
The main Ginkgo header ginkgo/ginkgo.hpp pulls in the core library headers:
executors, matrices, solvers, preconditioners and stopping criteria. The
other headers are only used for command-line parsing and I/O.
#include <ginkgo/ginkgo.hpp>
#include <fstream>
#include <iostream>
#include <map>
#include <string>
Ginkgo’s class names are heavily templated on value type / index type, so local aliases pay for themselves quickly:
using ValueType = double;
using RealValueType = gko::remove_complex<ValueType>;
using IndexType = int;
using vec = gko::matrix::Dense<ValueType>;
using real_vec = gko::matrix::Dense<RealValueType>;
using mtx = gko::matrix::Csr<ValueType, IndexType>;
using cg = gko::solver::Cg<ValueType>;
gko::matrix::Csr is the recommended sparse format for SPD problems like
this one. Swap to Coo, Ell, Sellp, Hybrid, Fbcsr, etc. and the
rest of the program is unchanged: every matrix format implements the same
LinOp interface. See
Matrix formats.
Pick an executor#
Map the command-line [executor] argument to a factory. The
ReferenceExecutor (single-threaded, fully portable) is always available;
the others depend on which backends were enabled at Ginkgo build time.
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)();
The point of the executor abstraction is exactly this: choose exec here,
and every kernel below dispatches automatically to that backend. Nothing
else in the program changes. See
Execution backends for the
broader picture.
Read the system from disk#
gko::read<T>(stream, exec) reads a Matrix Market file and produces a
T on exec. We share A with a shared_ptr because the solver factory
takes shared ownership; b and x0 are read once and not shared.
auto A = gko::share(gko::read<mtx>(std::ifstream("data/A.mtx"), exec));
auto b = gko::read<vec>(std::ifstream("data/b.mtx"), exec);
auto x = gko::read<vec>(std::ifstream("data/x0.mtx"), exec);
If you want to build the system in-process instead of reading it from disk, the Poisson on CPU tutorial walks through it.
Build the CG solver factory#
The factory describes how to build a solver; calling generate(A) on it
actually builds one bound to A. Two stopping criteria are combined: an
absolute iteration cap so a divergent run can’t loop forever, and a
relative residual reduction by \(10^{-7}\):
const RealValueType reduction_factor{1e-7};
auto solver_gen =
cg::build()
.with_criteria(
gko::stop::Iteration::build().with_max_iters(20u),
gko::stop::ResidualNorm<ValueType>::build()
.with_reduction_factor(reduction_factor))
.on(exec);
For other Krylov methods (Bicgstab, Gmres, Fcg, …) just swap the
factory; the criteria don’t change. See
Solver taxonomy.
Solve and verify#
generate(A) returns a unique_ptr<cg> (the bound solver); apply(b, x)
runs the iteration and writes the result into x. The residual check uses
Ginkgo’s own SpMV via A->apply(one, x, neg_one, b), which computes
\(b \leftarrow A x - b\) in place; then compute_norm2 gives us
\(\lVert A x - b \rVert_2\).
auto solver = solver_gen->generate(A);
solver->apply(b, x);
std::cout << "Solution (x):\n";
write(std::cout, x);
auto one = gko::initialize<vec>({1.0}, exec);
auto neg_one = gko::initialize<vec>({-1.0}, exec);
auto res = gko::initialize<real_vec>({0.0}, exec);
A->apply(one, x, neg_one, b);
b->compute_norm2(res);
std::cout << "Residual norm sqrt(r^T r):\n";
write(std::cout, res);
The one / neg_one scalars are themselves stored on exec so the
expression doesn’t trigger a host round-trip on a GPU executor. See
Cross-executor data movement
for the broader pattern.
Build and run#
CMakeLists.txt:
cmake_minimum_required(VERSION 3.16)
project(simple_solver CXX)
find_package(Ginkgo REQUIRED)
add_executable(simple-solver simple-solver.cpp)
target_link_libraries(simple-solver PRIVATE Ginkgo::ginkgo)
The example’s three Matrix Market files live next to the source in the
upstream repo (examples/simple-solver/data/). Copy them to a data/
directory next to the executable, then:
cmake -B build && cmake --build build
./build/simple-solver # reference (default)
./build/simple-solver omp # multi-threaded CPU
./build/simple-solver cuda # NVIDIA GPU
Results#
The Matrix Market files ship with a tiny 19-row SPD system. Expected console output (truncated):
Solution (x):
%%MatrixMarket matrix array real general
19 1
0.252218
0.108645
...
Residual norm sqrt(r^T r):
%%MatrixMarket matrix array real general
1 1
2.10788e-15
The exact final digit varies with compiler / backend / scheduling; any residual below ~\(10^{-13}\) on this problem means CG converged to round-off.
Variations to try#
Different solver. Replace
gko::solver::CgwithFcg, another method for SPD systems, or withBicgstab,Gmres, orCgs, which also handle non-symmetric and indefinite systems.Different matrix format. Change
using mtx = gko::matrix::Csr<...>;toCoo,Ell,Sellp,Hybrid, orFbcsr; everything else stays the same.Add a preconditioner. Drop in a
.with_preconditioner(...)call; the Preconditioned solver example does exactly that.Larger problem. Read your own Matrix Market file. The SuiteSparse Matrix Collection is a good source of standard test problems.
The plain program#
// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors
//
// SPDX-License-Identifier: BSD-3-Clause
// @sect3{Include files}
// This is the main ginkgo header file.
#include <ginkgo/ginkgo.hpp>
// Add the fstream header to read from data from files.
#include <fstream>
// Add the C++ iostream header to output information to the console.
#include <iostream>
// Add the STL map header for the executor selection
#include <map>
// Add the string manipulation header to handle strings.
#include <string>
int main(int argc, char* argv[])
{
// Use some shortcuts. In Ginkgo, vectors are seen as a gko::matrix::Dense
// with one column/one row. The advantage of this concept is that using
// multiple vectors is a now a natural extension of adding columns/rows are
// necessary.
using ValueType = double;
using RealValueType = gko::remove_complex<ValueType>;
using IndexType = int;
using vec = gko::matrix::Dense<ValueType>;
using real_vec = gko::matrix::Dense<RealValueType>;
// The gko::matrix::Csr class is used here, but any other matrix class such
// as gko::matrix::Coo, gko::matrix::Hybrid, gko::matrix::Ell or
// gko::matrix::Sellp could also be used.
using mtx = gko::matrix::Csr<ValueType, IndexType>;
// The gko::solver::Cg is used here, but any other solver class can also be
// used.
using cg = gko::solver::Cg<ValueType>;
// Print the ginkgo version information.
std::cout << gko::version_info::get() << std::endl;
// Print help on how to execute this example.
if (argc == 2 && (std::string(argv[1]) == "--help")) {
std::cerr << "Usage: " << argv[0] << " [executor] " << std::endl;
std::exit(-1);
}
// @sect3{Where do you want to run your solver ?}
// The gko::Executor class is one of the cornerstones of Ginkgo. Currently,
// we have support for
// an gko::OmpExecutor, which uses OpenMP multi-threading in most of its
// kernels, a gko::ReferenceExecutor, a single threaded specialization of
// the OpenMP executor and a gko::CudaExecutor which runs the code on a
// NVIDIA GPU if available.
// @note With the help of C++, you see that you only ever need to change the
// executor and all the other functions/ routines within Ginkgo should
// automatically work and run on the executor with any other changes.
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(); }}};
// executor where Ginkgo will perform the computation
const auto exec = exec_map.at(executor_string)(); // throws if not valid
// @sect3{Reading your data and transfer to the proper device.}
// Read the matrix, right hand side and the initial solution using the @ref
// read function.
// @note Ginkgo uses C++ smart pointers to automatically manage memory. To
// this end, we use our own object ownership transfer functions that under
// the hood call the required smart pointer functions to manage object
// ownership. gko::share and gko::give are the functions that you would need
// to use.
auto A = gko::share(gko::read<mtx>(std::ifstream("data/A.mtx"), exec));
auto b = gko::read<vec>(std::ifstream("data/b.mtx"), exec);
auto x = gko::read<vec>(std::ifstream("data/x0.mtx"), exec);
// @sect3{Creating the solver}
// Generate the gko::solver factory. Ginkgo uses the concept of Factories to
// build solvers with certain
// properties. Observe the Fluent interface used here. Here a cg solver is
// generated with a stopping criteria of maximum iterations of 20 and a
// residual norm reduction of 1e-7. You also observe that the stopping
// criteria(gko::stop) are also generated from factories using their build
// methods. You need to specify the executors which each of the object needs
// to be built on.
const RealValueType reduction_factor{1e-7};
auto solver_gen =
cg::build()
.with_criteria(gko::stop::Iteration::build().with_max_iters(20u),
gko::stop::ResidualNorm<ValueType>::build()
.with_reduction_factor(reduction_factor))
.on(exec);
// Generate the solver from the matrix. The solver factory built in the
// previous step takes a "matrix"(a gko::LinOp to be more general) as an
// input. In this case we provide it with a full matrix that we previously
// read, but as the solver only effectively uses the apply() method within
// the provided "matrix" object, you can effectively create a gko::LinOp
// class with your own apply implementation to accomplish more tasks. We
// will see an example of how this can be done in the custom-matrix-format
// example
auto solver = solver_gen->generate(A);
// Finally, solve the system. The solver, being a gko::LinOp, can be applied
// to a right hand side, b to
// obtain the solution, x.
solver->apply(b, x);
// Print the solution to the command line.
std::cout << "Solution (x):\n";
write(std::cout, x);
// To measure if your solution has actually converged, you can measure the
// error of the solution.
// one, neg_one are objects that represent the numbers which allow for a
// uniform interface when computing on any device. To compute the residual,
// all you need to do is call the apply method, which in this case is an
// spmv and equivalent to the LAPACK z_spmv routine. Finally, you compute
// the euclidean 2-norm with the compute_norm2 function.
auto one = gko::initialize<vec>({1.0}, exec);
auto neg_one = gko::initialize<vec>({-1.0}, exec);
auto res = gko::initialize<real_vec>({0.0}, exec);
A->apply(one, x, neg_one, b);
b->compute_norm2(res);
std::cout << "Residual norm sqrt(r^T r):\n";
write(std::cout, res);
}