test_amgcl_mpi.cpp

#include <vector>
#include <iostream>

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

#include <amgcl/mpi/distributed_matrix.hpp>
#include <amgcl/mpi/make_solver.hpp>
#include <amgcl/mpi/amg.hpp>
#include <amgcl/mpi/coarsening/smoothed_aggregation.hpp>
#include <amgcl/mpi/relaxation/spai0.hpp>
#include <amgcl/mpi/solver/bicgstab.hpp>

#include <amgcl/io/binary.hpp>
#include <amgcl/profiler.hpp>

#include <Eigen/Sparse>

// Block size
// const int B = 1;
// Get the global size of the matrix:
int rows = 5;
//---------------------------------------------------------------------------
int main(int argc, char *argv[])
{
    amgcl::mpi::init mpi(&argc, &argv);
    amgcl::mpi::communicator world(MPI_COMM_WORLD);

    int localsize = world.rank == 0 ? 3 : 2;
    Eigen::SparseMatrix<double> eA(localsize, rows);
    if (world.rank == 0)
    {
        eA.coeffRef(0, 0) = 1.0;
        eA.coeffRef(1, 1) = 1.0;
        eA.coeffRef(2, 2) = 1.0;
    }
    else
    {
        eA.coeffRef(0, 3) = 1.0;
        eA.coeffRef(1, 4) = 1.0;
    }

    eA.makeCompressed();
    // The profiler:
    amgcl::profiler<> prof("Serena MPI");

    prof.tic("read");

    int chunk = localsize;

    // Read our part of the system matrix.
    int nonzeros = eA.nonZeros();
    std::vector<int> ptr(eA.outerIndexPtr(), eA.outerIndexPtr() + localsize + 1);
    std::vector<int> col(eA.innerIndexPtr(), eA.innerIndexPtr() + nonzeros);
    std::vector<double> val(eA.valuePtr(), eA.valuePtr() + nonzeros);

    if (world.rank == 1) {
        std::cout << "nonzeros: " << nonzeros << std::endl;
        std::cout << "ptr: "; for(auto p : ptr) std::cout << p << " "; std::cout << std::endl;
        std::cout << "col: "; for(auto c : col) std::cout << c << " "; std::cout << std::endl;
        std::cout << "val: "; for(auto v : val) std::cout << v << " "; std::cout << std::endl;
    }

    prof.toc("read");

    // Declare the backend and the solver types
    typedef amgcl::backend::builtin<double> DBackend;
    typedef amgcl::backend::builtin<float> FBackend;

    typedef amgcl::mpi::make_solver<
        amgcl::mpi::amg<
            FBackend,
            amgcl::mpi::coarsening::smoothed_aggregation<FBackend>,
            amgcl::mpi::relaxation::spai0<FBackend>>,
        amgcl::mpi::solver::bicgstab<DBackend>>
        Solver;

    // Solver parameters
    Solver::params prm;
    prm.solver.maxiter = 200;

    auto A = std::make_shared<amgcl::mpi::distributed_matrix<DBackend>>(
        world, std::tie(chunk, ptr, col, val));

    std::vector<double> rhs(chunk, 1.0);

    // Initialize the solver:
    prof.tic("setup");
    Solver solve(world, A, prm);
    prof.toc("setup");

    // Show the mini-report on the constructed solver:
    if (world.rank == 0)
        std::cout << solve << std::endl;

    // Solve the system with the zero initial approximation:
    int iters;
    double error;
    std::vector<double> x(chunk, 0.0);

    prof.tic("solve");
    std::tie(iters, error) = solve(*A, rhs, x);
    prof.toc("solve");

    // Output the number of iterations, the relative error,
    // and the profiling data:
    if (world.rank == 0)
    {
        std::cout
            << "Iterations: " << iters << std::endl
            << "Error:      " << error << std::endl
            << prof << std::endl;

        std::cout << "sol: " << std::endl;
        for (auto &xi : x)
            std::cout << xi << std::endl;
    }
}