The custom-stopping-criterion program

The custom stopping criterion creation example..

This example depends on simple-solver, minimal-cuda-solver.

Table of contents
Introduction About the example The commented program	Results Comments about programming and debugging The plain program

Introduction

About the example

The commented program

    };
    GKO_ENABLE_CRITERION_FACTORY(ByInteraction, parameters, Factory);
    GKO_ENABLE_BUILD_METHOD(Factory);
 
protected:
    bool check_impl(gko::uint8 stoppingId, bool setFinalized,
                    gko::array<gko::stopping_status>* stop_status,
                    bool* one_changed, const Criterion::Updater&) override
    {
        bool result = *(parameters_.stop_iteration_process);
        if (result) {
            this->set_all_statuses(stoppingId, setFinalized, stop_status);
            *one_changed = true;
        }
        return result;
    }
 
    explicit ByInteraction(std::shared_ptr<const gko::Executor> exec)
        : EnablePolymorphicObject<ByInteraction, Criterion>(std::move(exec))
    {}
 
    explicit ByInteraction(const Factory* factory,
                           const gko::stop::CriterionArgs& args)
 
        : EnablePolymorphicObject<ByInteraction, Criterion>(
              factory->get_executor()),
          parameters_{factory->get_parameters()}
    {}
};
 
 
void run_solver(volatile bool* stop_iteration_process,
                std::shared_ptr<gko::Executor> exec)
{

Some shortcuts

using ValueType = double;
using RealValueType = gko::remove_complex<ValueType>;
using IndexType = int;
 
using mtx = gko::matrix::Csr<ValueType, IndexType>;
using vec = gko::matrix::Dense<ValueType>;
using real_vec = gko::matrix::Dense<RealValueType>;
using bicg = gko::solver::Bicgstab<ValueType>;

Read Data

auto A = 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);

Create solver factory and solve system

auto solver =
    bicg::build()
        .with_criteria(ByInteraction::build().with_stop_iteration_process(
            stop_iteration_process))
        .on(exec)
        ->generate(A);
solver->add_logger(gko::log::Stream<ValueType>::create(
    gko::log::Logger::iteration_complete_mask, std::cout, true));
solver->apply(b, x);
 
std::cout << "Solver stopped" << std::endl;

Print solution

std::cout << "Solution (x): \n";
write(std::cout, x);

Calculate residual

    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);
}
 
 
int main(int argc, char* argv[])
{

Print version information

std::cout << gko::version_info::get() << std::endl;

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);
}

Figure out where to run the code

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

Declare a user controlled boolean for the iteration process

volatile bool stop_iteration_process{};

Create a new a thread to launch the solver

std::thread t(run_solver, &stop_iteration_process, exec);

Look for an input command "stop" in the console, which sets the boolean to true

    std::cout << "Type 'stop' to stop the iteration process" << std::endl;
    std::string command;
    while (std::cin >> command) {
        if (command == "stop") {
            break;
        } else {
            std::cout << "Unknown command" << std::endl;
        }
    }
    std::cout << "User input command 'stop' - The solver will stop!"
              << std::endl;
    stop_iteration_process = true;
    t.join();
}

Results

This is the expected output:

.
.
.
.
.
.
[LOG] >>> iteration 22516 completed with solver LinOp[gko::solver::Bicgstab<double>,0x7fe6a4003710] with residual LinOp[gko::matrix::Dense<double>,0x7fe6a40050b0], solution LinOp[gko::matrix::Dense<double>,0x7fe6a40048e0] and residual_norm LinOp[gko::LinOp const*,0]
LinOp[gko::matrix::Dense<double>,0x7fe6a40050b0][
    5.17803e-164
    -7.6865e-165
    -2.06149e-164
    -4.84737e-165
    -3.36597e-164
    2.22353e-164
    1.47594e-165
    -1.78592e-165
    -6.17274e-166
    -3.02681e-166
    7.82009e-166
    8.57102e-165
    -1.28879e-164
    -2.62076e-165
    2.55329e-165
    -5.95988e-166
    -5.79273e-166
    -5.20172e-166
    -6.79458e-166
]
 
// Typing 'stop' stops the solver.
User input command 'stop' - The solver will stop
 
LinOp[gko::matrix::Dense<double>,0x7fe6a40048e0][
    0.252218
    0.108645
    0.0662811
    0.0630433
    0.0384088
    0.0396536
    0.0402648
    0.0338935
    0.0193098
    0.0234653
    0.0211499
    0.0196413
    0.0199151
    0.0181674
    0.0162722
    0.0150714
    0.0107016
    0.0121141
    0.0123025
]
 
Solver stopped
Solution (x):
%%MatrixMarket matrix array real general
19 1
0.252218
0.108645
0.0662811
0.0630433
0.0384088
0.0396536
0.0402648
0.0338935
0.0193098
0.0234653
0.0211499
0.0196413
0.0199151
0.0181674
0.0162722
0.0150714
0.0107016
0.0121141
0.0123025
Residual norm sqrt(r^T r):
%%MatrixMarket matrix array real general
1 1
6.50306e-16

Comments about programming and debugging

The plain program

#include <fstream>
#include <iostream>
#include <map>
#include <string>
#include <thread>
 
#include <ginkgo/ginkgo.hpp>
 
 
class ByInteraction
    : public gko::EnablePolymorphicObject<ByInteraction, gko::stop::Criterion> {
    friend class gko::EnablePolymorphicObject<ByInteraction,
                                              gko::stop::Criterion>;
    using Criterion = gko::stop::Criterion;
 
public:
    GKO_CREATE_FACTORY_PARAMETERS(parameters, Factory)
    {
        std::add_pointer<volatile bool>::type GKO_FACTORY_PARAMETER_SCALAR(
            stop_iteration_process, nullptr);
    };
    GKO_ENABLE_CRITERION_FACTORY(ByInteraction, parameters, Factory);
    GKO_ENABLE_BUILD_METHOD(Factory);
 
protected:
    bool check_impl(gko::uint8 stoppingId, bool setFinalized,
                    gko::array<gko::stopping_status>* stop_status,
                    bool* one_changed, const Criterion::Updater&) override
    {
        bool result = *(parameters_.stop_iteration_process);
        if (result) {
            this->set_all_statuses(stoppingId, setFinalized, stop_status);
            *one_changed = true;
        }
        return result;
    }
 
    explicit ByInteraction(std::shared_ptr<const gko::Executor> exec)
        : EnablePolymorphicObject<ByInteraction, Criterion>(std::move(exec))
    {}
 
    explicit ByInteraction(const Factory* factory,
                           const gko::stop::CriterionArgs& args)
 
        : EnablePolymorphicObject<ByInteraction, Criterion>(
              factory->get_executor()),
          parameters_{factory->get_parameters()}
    {}
};
 
 
void run_solver(volatile bool* stop_iteration_process,
                std::shared_ptr<gko::Executor> exec)
{
    using ValueType = double;
    using RealValueType = gko::remove_complex<ValueType>;
    using IndexType = int;
 
    using mtx = gko::matrix::Csr<ValueType, IndexType>;
    using vec = gko::matrix::Dense<ValueType>;
    using real_vec = gko::matrix::Dense<RealValueType>;
    using bicg = gko::solver::Bicgstab<ValueType>;
 
    auto A = 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);
 
    auto solver =
        bicg::build()
            .with_criteria(ByInteraction::build().with_stop_iteration_process(
                stop_iteration_process))
            .on(exec)
            ->generate(A);
    solver->add_logger(gko::log::Stream<ValueType>::create(
        gko::log::Logger::iteration_complete_mask, std::cout, true));
    solver->apply(b, x);
 
    std::cout << "Solver stopped" << std::endl;
 
    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);
}
 
 
int main(int argc, char* argv[])
{
    std::cout << gko::version_info::get() << std::endl;
 
    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
 
    volatile bool stop_iteration_process{};
 
    std::thread t(run_solver, &stop_iteration_process, exec);
 
    std::cout << "Type 'stop' to stop the iteration process" << std::endl;
    std::string command;
    while (std::cin >> command) {
        if (command == "stop") {
            break;
        } else {
            std::cout << "Unknown command" << std::endl;
        }
    }
    std::cout << "User input command 'stop' - The solver will stop!"
              << std::endl;
    stop_iteration_process = true;
    t.join();
}