Makefile

solve: solve.cpp
	g++ -o solve -std=c++11 -O3 solve.cpp          \
	    -I${HOME}/work/amgcl -I/usr/include/eigen3 \
	    -l boost_system

solve.cpp

#include <iostream>
#include <fstream>
#include <iterator>
#include <chrono>
#include <cassert>

#include <Eigen/Dense>
#include <Eigen/SparseCore>

#include <amgcl/amgcl.hpp>
#include <amgcl/interp_smoothed_aggr.hpp>
#include <amgcl/aggr_plain.hpp>
#include <amgcl/level_cpu.hpp>
#include <amgcl/operations_eigen.hpp>
#include <amgcl/cg.hpp>
#include <amgcl/profiler.hpp>

typedef Eigen::Triplet<double> triplet;

namespace Eigen {
    std::istream& operator>>(std::istream &is, triplet &t) {
        int i, j;
        double v;
        is >> i >> j >> v;
        t = triplet(i, j, v);
        return is;
    }
}

int main() {
    amgcl::profiler<std::chrono::steady_clock> prof("Total");

    // Read RHS and A.
    prof.tic("read");
    Eigen::VectorXd rhs;
    {
        std::ifstream f("sample_rhs.txt");
        std::istream_iterator<double> bf(f);
        std::istream_iterator<double> ef;
        std::vector<double> v(bf, ef);
        rhs = Eigen::Map<Eigen::VectorXd>(v.data(), v.size());
    }
    std::cout << "#unknowns: " << rhs.size() << std::endl;

    const int n = rhs.size();

    Eigen::SparseMatrix<double> A(n, n);
    {
        std::vector<triplet> triplets;
        std::ifstream f("sample_matrix.txt");
        std::istream_iterator<triplet> bf(f);
        std::istream_iterator<triplet> ef;
        triplets.assign(bf, ef);

        A.setFromTriplets(triplets.begin(), triplets.end());
    }
    std::cout << "#nonzeros: " << A.nonZeros() << std::endl;
    prof.toc("read");

#if 0
    // Scale the diagonal
    for(int i = 0; i < n; ++i)
        A.coeffRef(i,i) *= 1.1;
#endif

    // Build the preconditioner:
    prof.tic("setup");
    typedef amgcl::solver<
        double, int,
        amgcl::interp::smoothed_aggregation<amgcl::aggr::plain>,
        amgcl::level::cpu<amgcl::relax::ilu0>
        > AMG;

    AMG amg( amgcl::sparse::map(A) );
    prof.toc("setup");

    std::cout << amg << std::endl;

    prof.tic("solve");
    Eigen::VectorXd x = Eigen::VectorXd::Ones(n);
    auto cnv = amgcl::solve(A, rhs, amg, x, amgcl::cg_tag());
    prof.toc("solve");

    std::cout << "Iterations: " << cnv.first  << std::endl
              << "Error:      " << cnv.second << std::endl
              << std::endl;

    std::cout << prof << std::endl;
}