distributed.cpp

#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include <string>

#include <amgcl/backend/builtin.hpp>
#include <amgcl/adapter/crs_tuple.hpp>

#include <amgcl/mpi/util.hpp>
#include <amgcl/mpi/make_solver.hpp>
#include <amgcl/mpi/coarsening/aggregation.hpp>
#include <amgcl/mpi/relaxation/spai0.hpp>
#include <amgcl/mpi/solver/cg.hpp>
#include <amgcl/mpi/amg.hpp>

typedef amgcl::backend::builtin<double> Backend;
typedef amgcl::mpi::make_solver<
    amgcl::mpi::amg<
        Backend,
        amgcl::mpi::coarsening::aggregation<Backend>,
        amgcl::mpi::relaxation::spai0<Backend>>,
        amgcl::mpi::solver::cg<Backend>
    > Solver;

typedef Backend::matrix Matrix;
typedef amgcl::mpi::distributed_matrix<Backend> DMatrix;

std::shared_ptr<Matrix> read_matrix(int step, int rank) {
    std::ostringstream fname;
    fname << "step" << step << "/A_" << rank << ".dat";

    std::vector<ptrdiff_t> ptr, col;
    std::vector<double> val;

    ptr.push_back(0);
    ptrdiff_t last_row = 0;

    std::ifstream f(fname.str());
    while (true) {
        ptrdiff_t r, c;
        double v;

        if (!(f >> r >> c >> v)) break;

        if (r != last_row) ptr.push_back(col.size());

        last_row = r;
        col.push_back(c);
        val.push_back(v);
    }

    ptr.push_back(col.size());
    ptrdiff_t n = ptr.size() - 1;

    return std::make_shared<Matrix>(std::tie(n, ptr, col, val));
}

std::vector<double> read_vector(int step, int rank) {
    std::ostringstream fname;
    fname << "step" << step << "/b_" << rank << ".dat";

    std::vector<double> rhs;

    std::ifstream f(fname.str());
    while (true) {
        double v;
        if (!(f >> v)) break;
        rhs.push_back(v);
    }

    return rhs;
}

int main(int argc, char *argv[]) {
    int nsteps = 2;

    amgcl::mpi::init_thread mpi(&argc, &argv);
    amgcl::mpi::communicator comm(MPI_COMM_WORLD);

    amgcl::precondition(comm.size == 4, "mpi world should have 4 processes!");
    
    Backend::params bprm;
    Solver::params prm;

    prm.precond.allow_rebuild = true;
    prm.solver.tol = 1e-7;
    prm.solver.abstol = 1e-2;
    prm.solver.maxiter = 1000;

    std::shared_ptr<Solver> solver;

    for (int step = 0; step < nsteps; ++step) {
        if (comm.rank == 0) {
            std::cout << "step " << step << std::endl;
        }

        auto A = read_matrix(step, comm.rank);
        auto b = read_vector(step, comm.rank);
        std::vector<double> x(b.size());

        comm.check(A->nrows == b.size(), "matrix and vector have inconsistent sizes");

        auto Ad = std::make_shared<DMatrix>(comm, *A);

        if (solver) {
            if (comm.rank == 0) {
                std::cout << "rebuilding solver..." << std::endl;
            }
            solver->precond().rebuild(Ad);
        } else {
            solver = std::make_shared<Solver>(comm, Ad, prm, bprm);
        }
        if (comm.rank == 0) {
            std::cout << *solver << std::endl;
        }

        int iters;
        double error;
        std::tie(iters, error) = (*solver)(b, x);

        if (comm.rank == 0) {
            std::cout << "iters=" << iters << " error=" << error << std::endl;
        }
    }
}

Makefile

amgcl_root=$(HOME)/work/amgcl

distributed: distributed.cpp
	mpic++ -o $@ -g $^ -I $(amgcl_root)

serial: serial.cpp
	g++ -o $@ -g $^ -I $(amgcl_root)

run_distr: distributed
	mpirun -np 4 ./$^

run_serial: serial
	./$^

serial.cpp

#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include <string>

#include <amgcl/backend/builtin.hpp>
#include <amgcl/adapter/crs_tuple.hpp>

#include <amgcl/make_solver.hpp>
#include <amgcl/coarsening/aggregation.hpp>
#include <amgcl/relaxation/spai0.hpp>
#include <amgcl/solver/cg.hpp>
#include <amgcl/amg.hpp>

typedef amgcl::backend::builtin<double> Backend;
typedef amgcl::make_solver<
    amgcl::amg<
        Backend,
        amgcl::coarsening::aggregation,
        amgcl::relaxation::spai0
    >,
    amgcl::solver::cg<Backend>
> Solver;

typedef Backend::matrix Matrix;

const int np = 4;

std::shared_ptr<Matrix> read_matrix(int step) {
    std::vector<ptrdiff_t> ptr, col;
    std::vector<double> val;
    ptr.push_back(0);

    ptrdiff_t last_row = 0;

    for (int rank = 0; rank < np; ++rank) {
        std::ostringstream fname;
        fname << "step" << step << "/A_" << rank << ".dat";

        std::ifstream f(fname.str());
        while (true) {
            ptrdiff_t r, c;
            double v;

            if (!(f >> r >> c >> v)) break;

            if (r != last_row) ptr.push_back(col.size());

            last_row = r;
            col.push_back(c);
            val.push_back(v);
        }
    }

    ptr.push_back(col.size());
    ptrdiff_t n = ptr.size() - 1;

    return std::make_shared<Matrix>(std::tie(n, ptr, col, val));
}

std::vector<double> read_vector(int step) {
    std::vector<double> rhs;

    for (int rank = 0; rank < np; ++rank) {
        std::ostringstream fname;
        fname << "step" << step << "/b_" << rank << ".dat";

        std::ifstream f(fname.str());
        while (true) {
            double v;
            if (!(f >> v)) break;
            rhs.push_back(v);
        }
    }

    return rhs;
}

int main() {
    int nsteps = 2;

    Backend::params bprm;
    Solver::params prm;

    prm.precond.allow_rebuild = true;
    prm.solver.tol = 1e-7;
    prm.solver.abstol = 1e-2;
    prm.solver.maxiter = 1000;

    std::shared_ptr<Solver> solver;

    for (int step = 0; step < nsteps; ++step) {
        std::cout << "step " << step << std::endl;

        auto A = read_matrix(step);
        auto b = read_vector(step);
        std::vector<double> x(b.size());

        amgcl::precondition(A->nrows == b.size(), "matrix and vector have inconsistent sizes");

        if (solver) {
            std::cout << "rebuilding solver..." << std::endl;
            solver->precond().rebuild(*A);
        } else {
            solver = std::make_shared<Solver>(*A, prm, bprm);
        }
        std::cout << *solver << std::endl;

        int iters;
        double error;
        std::tie(iters, error) = (*solver)(b, x);

        std::cout << "iters=" << iters << " error=" << error << std::endl;
    }
}