amgcl.cpp

#include <vector>
#include <tuple>

#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::builtin<ValueType, ColumnType, PointerType> Backend;

    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);
    prof.toc("setup");
    std::cout << solve << std::endl;

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

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

Makefile

amgcl: amgcl.cpp
	g++ -O3 -fopenmp -DNDEBUG -I$(HOME)/work/amgcl -o amgcl amgcl.cpp

./amgcl 64
Solver
======
Type:             BiCGStab
Unknowns:         262144
Memory footprint: 14.00 M

Preconditioner
==============
Number of levels:    3
Operator complexity: 1.61
Grid complexity:     1.13
Memory footprint:    103.90 M

level     unknowns       nonzeros      memory
---------------------------------------------
    0       262144        1810432     68.52 M (62.20%)
    1        31868         955994     20.34 M (32.84%)
    2         2279         144245     15.04 M ( 4.96%)

Iterations: 6
Error     : 6.59412e-07

[Profile:     1.010 s] (100.00%)
[  setup:     0.729 s] ( 72.23%)
[  solve:     0.280 s] ( 27.75%)