test_amgcl_mpi.cpp

#include <vector>
#include <iostream>

#include <amgcl/backend/builtin.hpp>
#include <amgcl/value_type/static_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>

// 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;
    std::vector<ptrdiff_t> ptr, col;
    std::vector<double> val;
    ptr.push_back(0);
    if (world.rank == 0)
    {
        col.push_back(0); val.push_back(1.0); ptr.push_back(col.size());
        col.push_back(1); val.push_back(1.0); ptr.push_back(col.size());
        col.push_back(2); val.push_back(1.0); ptr.push_back(col.size());
    }
    else
    {
        col.push_back(3); val.push_back(1.0); ptr.push_back(col.size());
        col.push_back(4); val.push_back(1.0); ptr.push_back(col.size());
    }

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

    // 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;

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

    // Initialize the solver:
    auto A = std::make_shared<amgcl::mpi::distributed_matrix<DBackend>>(
        world, std::tie(localsize, ptr, col, val));

    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(localsize, 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;
    }
}