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

#include <ginkgo/ginkgo.hpp>


#include <fstream>
#include <iostream>
#include <string>
#include <thread>


/ **
 * The ByInteraction class is a criterion which asks for user input to stop
 * the iteration process. Using this criterion is slightly more complex than the
 * other ones, because it is asynchronous therefore requires the use of threads.
 * /
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)
    {
        / **
         * Boolean set by the user to stop the iteration process
         * /
        std::add_pointer<volatile bool>::type GKO_FACTORY_PARAMETER(
            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)
{

Some shortcuts

using mtx = gko::matrix::Csr<>;
using vec = gko::matrix::Dense<>;
using bicg = gko::solver::Bicgstab<>;

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))
                  .on(exec)
                  ->generate(A);
solver->add_logger(gko::log::Stream<>::create(
    exec, gko::log::Logger::iteration_complete_mask, std::cout, true));
solver->apply(lend(b), lend(x));

std::cout << "Solver stopped" << std::endl;

Print solution

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

Calculate residual

    auto one = gko::initialize<vec>({1.0}, exec);
    auto neg_one = gko::initialize<vec>({-1.0}, exec);
    auto res = gko::initialize<vec>({0.0}, exec);
    A->apply(lend(one), lend(x), lend(neg_one), lend(b));
    b->compute_norm2(lend(res));

    std::cout << "Residual norm sqrt(r^T r): \n";
    write(std::cout, lend(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

std::shared_ptr<gko::Executor> exec;
if (argc == 1 || std::string(argv[1]) == "reference") {
    exec = gko::ReferenceExecutor::create();
} else if (argc == 2 && std::string(argv[1]) == "omp") {
    exec = gko::OmpExecutor::create();
} else if (argc == 2 && std::string(argv[1]) == "cuda" &&
           gko::CudaExecutor::get_num_devices() > 0) {
    exec = gko::CudaExecutor::create(0, gko::OmpExecutor::create());
} else {
    std::cerr << "Usage: " << argv[0] << " [executor]" << std::endl;
    std::exit(-1);
}

Declare a user controled 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 15331 completed with solver LinOp[gko::solver::Bicgstab<double>,0x7f2f38003c10] with residual LinOp[gko::matrix::Dense<double>,0x7f2f380048e0], solution LinOp[gko::LinOp const*,0] and residual_norm LinOp[gko::LinOp const*,0]
LinOp[gko::matrix::Dense<double>,0x7f2f380048e0][
        6.21705e-164
        -1.18919e-164
        7.89129e-165
        -6.78013e-165
        -2.42405e-164
        -4.29503e-165
        6.16567e-166
        -3.34064e-164
        6.38335e-165
        7.86768e-165
        -1.80969e-165
        -4.17609e-166                   
        2.5395e-165
        -5.34283e-166
        -4.10476e-166
        -1.50132e-166
        -1.25732e-165
        -1.82819e-166
        -2.0927e-165
]        

// Typing 'stop' stops the solver.

User input command 'stop' - The solver will stop

LinOp[gko::matrix::Dense<double>,0x7f2f38004730][
        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
1.06135e-15

Comments about programming and debugging

The plain program

/*******************************<GINKGO LICENSE>******************************
Copyright (c) 2017-2019, the Ginkgo authors
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:

1. Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.

3. Neither the name of the copyright holder nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
******************************<GINKGO LICENSE>*******************************/

#include <ginkgo/ginkgo.hpp>


#include <fstream>
#include <iostream>
#include <string>
#include <thread>


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(
            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 mtx = gko::matrix::Csr<>;
    using vec = gko::matrix::Dense<>;
    using bicg = gko::solver::Bicgstab<>;

    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))
                      .on(exec)
                      ->generate(A);
    solver->add_logger(gko::log::Stream<>::create(
        exec, gko::log::Logger::iteration_complete_mask, std::cout, true));
    solver->apply(lend(b), lend(x));

    std::cout << "Solver stopped" << std::endl;

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

    auto one = gko::initialize<vec>({1.0}, exec);
    auto neg_one = gko::initialize<vec>({-1.0}, exec);
    auto res = gko::initialize<vec>({0.0}, exec);
    A->apply(lend(one), lend(x), lend(neg_one), lend(b));
    b->compute_norm2(lend(res));

    std::cout << "Residual norm sqrt(r^T r): \n";
    write(std::cout, lend(res));
}


int main(int argc, char *argv[])
{
    std::cout << gko::version_info::get() << std::endl;

    std::shared_ptr<gko::Executor> exec;
    if (argc == 1 || std::string(argv[1]) == "reference") {
        exec = gko::ReferenceExecutor::create();
    } else if (argc == 2 && std::string(argv[1]) == "omp") {
        exec = gko::OmpExecutor::create();
    } else if (argc == 2 && std::string(argv[1]) == "cuda" &&
               gko::CudaExecutor::get_num_devices() > 0) {
        exec = gko::CudaExecutor::create(0, gko::OmpExecutor::create());
    } else {
        std::cerr << "Usage: " << argv[0] << " [executor]" << std::endl;
        std::exit(-1);
    }

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