Ginkgo  Generated from pipelines/1556235455 branch based on develop. Ginkgo version 1.9.0
A numerical linear algebra library targeting many-core architectures
The inverse-iteration program

The inverse iteration 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 shows how components available in Ginkgo can be used to implement higher-level numerical methods. The method used here will be the shifted inverse iteration method for eigenvalue computation which find the eigenvalue and eigenvector of A closest to z, for some scalar z. The method requires repeatedly solving the shifted linear system (A - zI)x = b, as well as performing matrix-vector products with the matrix A. Here is the complete pseudocode of the method:

x_0 = initial guess
for i = 0 .. max_iterations:
solve (A - zI) y_i = x_i for y_i+1
x_(i+1) = y_i / || y_i || # compute next eigenvector approximation
g_(i+1) = x_(i+1)^* A x_(i+1) # approximate eigenvalue (Rayleigh quotient)
if ||A x_(i+1) - g_(i+1)x_(i+1)|| < tol * g_(i+1): # check convergence
break

About the example

The commented program

Figure out where to run the code

if (argc == 2 && (std::string(argv[1]) == "--help")) {
std::cerr << "Usage: " << argv[0] << " [executor]" << std::endl;
std::exit(-1);
}
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
auto this_exec = exec->get_master();

linear system solver parameters

auto system_max_iterations = 100u;
auto system_residual_goal = real_precision{1e-16};

eigensolver parameters

auto max_iterations = 20u;
auto residual_goal = real_precision{1e-8};
auto z = precision{20.0, 2.0};

Read data

auto A = share(gko::read<mtx>(std::ifstream("data/A.mtx"), exec));

Generate shifted matrix A - zI

  • we avoid duplicating memory by not storing both A and A - zI, but compute A - zI on the fly by using Ginkgo's utilities for creating linear combinations of operators
auto one = share(gko::initialize<vec>({precision{1.0}}, exec));
auto neg_one = share(gko::initialize<vec>({-precision{1.0}}, exec));
auto neg_z = gko::initialize<vec>({-z}, exec);
auto system_matrix = share(gko::Combination<precision>::create(
one, A, gko::initialize<vec>({-z}, exec),
gko::matrix::Identity<precision>::create(exec, A->get_size()[0])));

Generate solver operator (A - zI)^-1

auto solver =
solver_type::build()
.with_criteria(gko::stop::Iteration::build().with_max_iters(
system_max_iterations),
gko::stop::ResidualNorm<precision>::build()
.with_reduction_factor(system_residual_goal))
.on(exec)
->generate(system_matrix);

inverse iterations

start with guess [1, 1, ..., 1]

auto x = [&] {
auto work = vec::create(this_exec, gko::dim<2>{A->get_size()[0], 1});
const auto n = work->get_size()[0];
for (int i = 0; i < n; ++i) {
work->get_values()[i] = precision{1.0} / sqrt(n);
}
return clone(exec, work);
}();
auto y = clone(x);
auto tmp = clone(x);
auto norm = gko::initialize<real_vec>({1.0}, exec);
auto inv_norm = clone(this_exec, one);
auto g = clone(one);
for (auto i = 0u; i < max_iterations; ++i) {
std::cout << "{ ";

(A - zI)y = x

solver->apply(x, y);
system_matrix->apply(one, y, neg_one, x);
x->compute_norm2(norm);
std::cout << "\"system_residual\": "
<< clone(this_exec, norm)->get_values()[0] << ", ";
x->copy_from(y);

x = y / || y ||

x->compute_norm2(norm);
inv_norm->get_values()[0] =
real_precision{1.0} / clone(this_exec, norm)->get_values()[0];
x->scale(clone(exec, inv_norm));

g = x^* A x

A->apply(x, tmp);
x->compute_dot(tmp, g);
auto g_val = clone(this_exec, g)->get_values()[0];
std::cout << "\"eigenvalue\": " << g_val << ", ";

||Ax - gx|| < tol * g

auto v = gko::initialize<vec>({-g_val}, exec);
tmp->add_scaled(v, x);
tmp->compute_norm2(norm);
auto res_val = clone(exec->get_master(), norm)->get_values()[0];
std::cout << "\"residual\": " << res_val / g_val << " }," << std::endl;
if (abs(res_val) < residual_goal * abs(g_val)) {
break;
}
}
}

Results

This is the expected output:

{ "system_residual": +1.61736920e-14, "eigenvalue": (+2.03741410e+01,-1.17744356e-16), "residual": (+2.92231055e-01,+1.68883476e-18) },
{ "system_residual": +4.98014795e-15, "eigenvalue": (+1.94878474e+01,+1.25948378e-15), "residual": (+7.94370276e-02,-5.13395071e-18) },
{ "system_residual": +3.39296916e-15, "eigenvalue": (+1.93282121e+01,-1.19329332e-15), "residual": (+4.11149623e-02,+2.53837290e-18) },
{ "system_residual": +3.35953656e-15, "eigenvalue": (+1.92638912e+01,+3.28657016e-16), "residual": (+2.34717040e-02,-4.00445585e-19) },
{ "system_residual": +2.91474009e-15, "eigenvalue": (+1.92409166e+01,+3.65597737e-16), "residual": (+1.34709547e-02,-2.55962367e-19) },
{ "system_residual": +3.09863953e-15, "eigenvalue": (+1.92331106e+01,-1.07919176e-15), "residual": (+7.72060707e-03,+4.33212063e-19) },
{ "system_residual": +2.31198069e-15, "eigenvalue": (+1.92305014e+01,-2.89755360e-16), "residual": (+4.42106625e-03,+6.66143651e-20) },
{ "system_residual": +3.02771202e-15, "eigenvalue": (+1.92296339e+01,+8.04259901e-16), "residual": (+2.53081312e-03,-1.05848687e-19) },
{ "system_residual": +2.02954523e-15, "eigenvalue": (+1.92293461e+01,+7.81834016e-16), "residual": (+1.44862114e-03,-5.88985854e-20) },
{ "system_residual": +2.31762332e-15, "eigenvalue": (+1.92292506e+01,-1.11718775e-16), "residual": (+8.29183451e-04,+4.81741912e-21) },
{ "system_residual": +8.12541038e-15, "eigenvalue": (+1.92292190e+01,-6.55606254e-16), "residual": (+4.74636702e-04,+1.61823936e-20) },
{ "system_residual": +2.77259926e-15, "eigenvalue": (+1.92292085e+01,+4.30588140e-16), "residual": (+2.71701077e-04,-6.08403935e-21) },
{ "system_residual": +8.87888675e-14, "eigenvalue": (+1.92292051e+01,+9.67936313e-18), "residual": (+1.55539937e-04,-7.82937998e-23) },
{ "system_residual": +2.85077117e-15, "eigenvalue": (+1.92292039e+01,-4.52923128e-16), "residual": (+8.90457139e-05,+2.09737561e-21) },
{ "system_residual": +6.46865302e-14, "eigenvalue": (+1.92292035e+01,+1.58710681e-17), "residual": (+5.09805252e-05,-4.20774259e-23) },
{ "system_residual": +4.18913713e-15, "eigenvalue": (+1.92292034e+01,+1.06839590e-15), "residual": (+2.91887365e-05,-1.62175862e-21) },
{ "system_residual": +1.06421578e-11, "eigenvalue": (+1.92292034e+01,+3.26089685e-17), "residual": (+1.67126561e-05,-2.83413965e-23) },
{ "system_residual": +2.97434420e-13, "eigenvalue": (+1.92292034e+01,-7.85427712e-16), "residual": (+9.56961199e-06,+3.90876227e-22) },
{ "system_residual": +1.63230281e-11, "eigenvalue": (+1.92292033e+01,+3.69307000e-16), "residual": (+5.47975753e-06,-1.05241636e-22) },
{ "system_residual": +6.14939758e-14, "eigenvalue": (+1.92292033e+01,+1.36057865e-15), "residual": (+3.13794996e-06,-2.22028320e-22) },

Comments about programming and debugging

The plain program

#include <cmath>
#include <complex>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <map>
#include <string>
#include <ginkgo/ginkgo.hpp>
int main(int argc, char* argv[])
{
using precision = std::complex<double>;
using real_precision = gko::remove_complex<precision>;
using vec = gko::matrix::Dense<precision>;
using real_vec = gko::matrix::Dense<real_precision>;
using mtx = gko::matrix::Csr<precision>;
using solver_type = gko::solver::Bicgstab<precision>;
using std::abs;
using std::sqrt;
std::cout << gko::version_info::get() << std::endl;
std::cout << std::scientific << std::setprecision(8) << std::showpos;
if (argc == 2 && (std::string(argv[1]) == "--help")) {
std::cerr << "Usage: " << argv[0] << " [executor]" << std::endl;
std::exit(-1);
}
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
auto this_exec = exec->get_master();
auto system_max_iterations = 100u;
auto system_residual_goal = real_precision{1e-16};
auto max_iterations = 20u;
auto residual_goal = real_precision{1e-8};
auto z = precision{20.0, 2.0};
auto A = share(gko::read<mtx>(std::ifstream("data/A.mtx"), exec));
auto one = share(gko::initialize<vec>({precision{1.0}}, exec));
auto neg_one = share(gko::initialize<vec>({-precision{1.0}}, exec));
auto neg_z = gko::initialize<vec>({-z}, exec);
auto system_matrix = share(gko::Combination<precision>::create(
one, A, gko::initialize<vec>({-z}, exec),
gko::matrix::Identity<precision>::create(exec, A->get_size()[0])));
auto solver =
solver_type::build()
.with_criteria(gko::stop::Iteration::build().with_max_iters(
system_max_iterations),
gko::stop::ResidualNorm<precision>::build()
.with_reduction_factor(system_residual_goal))
.on(exec)
->generate(system_matrix);
auto x = [&] {
auto work = vec::create(this_exec, gko::dim<2>{A->get_size()[0], 1});
const auto n = work->get_size()[0];
for (int i = 0; i < n; ++i) {
work->get_values()[i] = precision{1.0} / sqrt(n);
}
return clone(exec, work);
}();
auto y = clone(x);
auto tmp = clone(x);
auto norm = gko::initialize<real_vec>({1.0}, exec);
auto inv_norm = clone(this_exec, one);
auto g = clone(one);
for (auto i = 0u; i < max_iterations; ++i) {
std::cout << "{ ";
solver->apply(x, y);
system_matrix->apply(one, y, neg_one, x);
x->compute_norm2(norm);
std::cout << "\"system_residual\": "
<< clone(this_exec, norm)->get_values()[0] << ", ";
x->copy_from(y);
x->compute_norm2(norm);
inv_norm->get_values()[0] =
real_precision{1.0} / clone(this_exec, norm)->get_values()[0];
x->scale(clone(exec, inv_norm));
A->apply(x, tmp);
x->compute_dot(tmp, g);
auto g_val = clone(this_exec, g)->get_values()[0];
std::cout << "\"eigenvalue\": " << g_val << ", ";
auto v = gko::initialize<vec>({-g_val}, exec);
tmp->add_scaled(v, x);
tmp->compute_norm2(norm);
auto res_val = clone(exec->get_master(), norm)->get_values()[0];
std::cout << "\"residual\": " << res_val / g_val << " }," << std::endl;
if (abs(res_val) < residual_goal * abs(g_val)) {
break;
}
}
}
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::Combination
The Combination class can be used to construct a linear combination of multiple linear operators c1 *...
Definition: combination.hpp:31
gko::log::profile_event_category::solver
Solver events.
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:945
gko::solver::Bicgstab
BiCGSTAB or the Bi-Conjugate Gradient-Stabilized is a Krylov subspace solver.
Definition: bicgstab.hpp:51
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::stop::ResidualNorm
The ResidualNorm class is a stopping criterion which stops the iteration process when the actual resi...
Definition: residual_norm.hpp:111
gko::dim< 2 >
gko::share
detail::shared_type< OwningPointer > share(OwningPointer &&p)
Marks the object pointed to by p as shared.
Definition: utils_helper.hpp:224
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::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::Identity
This class is a utility which efficiently implements the identity matrix (a linear operator which map...
Definition: identity.hpp:35
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:325
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::one
constexpr T one()
Returns the multiplicative identity for T.
Definition: math.hpp:652
Generated by   doxygen 1.8.16