amgcl.cpp

#include <vector>
#include <tuple>

#include <amgcl/backend/vexcl.hpp>
#include <amgcl/make_solver.hpp>
#include <amgcl/preconditioner/runtime.hpp>
#include <amgcl/solver/runtime.hpp>
#include <amgcl/adapter/crs_tuple.hpp>
#include <amgcl/adapter/zero_copy.hpp>
#include <amgcl/profiler.hpp>

// Generates matrix for poisson problem in a unit cube.
template <typename PointerType, typename ColumnType, typename ValueType, typename RhsType>
int sample_problem(
        ptrdiff_t                n,
        std::vector<PointerType> &ptr,
        std::vector<ColumnType>  &col,
        std::vector<ValueType>   &val,
        std::vector<RhsType>     &rhs,
        double anisotropy = 1.0
        )
{
    ptrdiff_t n3  = n * n * n;

    ptr.clear();
    col.clear();
    val.clear();
    rhs.clear();

    ptr.reserve(n3 + 1);
    col.reserve(n3 * 7);
    val.reserve(n3 * 7);
    rhs.reserve(n3);

    ValueType one = amgcl::math::identity<ValueType>();

    double hx = 1;
    double hy = hx * anisotropy;
    double hz = hy * anisotropy;

    ptr.push_back(0);
    for(ptrdiff_t k = 0, idx = 0; k < n; ++k) {
        for(ptrdiff_t j = 0; j < n; ++j) {
            for (ptrdiff_t i = 0; i < n; ++i, ++idx) {
                if (k > 0) {
                    col.push_back(idx - n * n);
                    val.push_back(-1.0/(hz * hz) * one);
                }

                if (j > 0) {
                    col.push_back(idx - n);
                    val.push_back(-1.0/(hy * hy) * one);
                }

                if (i > 0) {
                    col.push_back(idx - 1);
                    val.push_back(-1.0/(hx * hx) * one);
                }

                col.push_back(idx);
                val.push_back((2 / (hx * hx) + 2 / (hy * hy) + 2 / (hz * hz)) * one);

                if (i + 1 < n) {
                    col.push_back(idx + 1);
                    val.push_back(-1.0/(hx * hx) * one);
                }

                if (j + 1 < n) {
                    col.push_back(idx + n);
                    val.push_back(-1.0/(hy * hy) * one);
                }

                if (k + 1 < n) {
                    col.push_back(idx + n * n);
                    val.push_back(-1.0/(hz * hz) * one);
                }

                rhs.push_back( amgcl::math::constant<RhsType>(1.0) );
                ptr.push_back( static_cast<PointerType>(col.size()) );
            }
        }
    }

    return n3;
}

//---------------------------------------------------------------------------
template <typename PointerType, typename ColumnType, typename ValueType, typename RhsType>
void solve(
        // Matrix in CRS format:
        const std::vector<PointerType> &ptr,
        const std::vector<ColumnType>  &col,
        const std::vector<ValueType>   &val,
        // RHS:
        const std::vector<RhsType>     &rhs,
        // Initial solution on input, solution on output:
        std::vector<RhsType>           &x
        )
{
    typedef amgcl::backend::vexcl<ValueType, ColumnType, PointerType> Backend;
    typename Backend::params bprm;

    vex::Context ctx(vex::Filter::Env && vex::Filter::Count(1));
    std::cout << ctx << std::endl;

    bprm.q = ctx;

    vex::vector<RhsType> X(ctx, x);
    vex::vector<RhsType> F(ctx, rhs);

    typedef amgcl::make_solver<
        amgcl::runtime::preconditioner<Backend>,
        amgcl::runtime::solver::wrapper<Backend>
        > Solver;

    boost::property_tree::ptree prm;
    prm.put("solver.tol", 1e-6);

    PointerType n = ptr.size() - 1;

    amgcl::profiler<> prof;

    // setup the solver:
    prof.tic("setup");
    auto A = amgcl::adapter::zero_copy_direct(n, ptr.data(), col.data(), val.data());
    amgcl::backend::sort_rows(*A);
    Solver solve(A, prm, bprm);
    prof.toc("setup");
    std::cout << solve << std::endl;

    // solve the problem:
    int    iters;
    double error;
    prof.tic("solve");
    std::tie(iters, error) = solve(F, X);
    prof.toc("solve");

    vex::copy(X, x);

    std::cout << "Iterations: " << iters << std::endl;
    std::cout << "Error     : " << error << std::endl;
    std::cout << prof << std::endl;
}

//---------------------------------------------------------------------------
int main(int argc, char *argv[]) {
    int n = argc > 1 ? std::stoi(argv[1]) : 32;

    std::vector<int> ptr;
    std::vector<int> col;
    std::vector<double> val, rhs;

    int n3 = sample_problem(n, ptr, col, val, rhs);

    std::vector<double> x(n3, 0.0);

    solve(ptr, col, val, rhs, x);
}

CMakeLists.txt

cmake_minimum_required(VERSION 3.17) 
project(hello)

find_package(amgcl)
find_package(VexCL)

add_executable(hello amgcl.cpp)
target_link_libraries(hello amgcl::amgcl)
target_link_libraries(hello VexCL::OpenCL)