amgcl icon indicating copy to clipboard operation
amgcl copied to clipboard

Difference between using make_shared on distributed matrix vs. solver

Open allaffa opened this issue 1 year ago • 1 comments

I see that there are different options proposed to solve a linear system with a distributed matrix in different examples.

Option 1 is described in this example: https://github.com/ddemidov/amgcl/blob/master/examples/mpi/mpi_solver.cpp

  1. Reading chunks of a matrix across different processes
  2. Use make_shared on the solver
solve = std::make_shared<Solver>(comm,
                amgcl::adapter::block_matrix<val_type>(std::tie(chunk, ptr, col, val)),
                prm, bprm);

Option 2 is described in this example: https://amgcl.readthedocs.io/en/latest/tutorial/poisson3DbMPI.html

  1. Reading chunks of a matrix across different processes
  2. Create the distributed matrix from the local parts.
auto A = std::make_shared<amgcl::mpi::distributed_matrix<DBackend>>(
            world, std::tie(chunk, ptr, col, val));
  1. Initialize the solver without calling make_shared anymore:
prof.tic("setup");
Solver solve(world, A);
prof.toc("setup");

I have two questions to ask:

Question 1: Are the two options equivalent, or should one be preferable over the other depending on the situation (e.g., scale, sparsity pattern, ...)? Question 2: In the example https://amgcl.readthedocs.io/en/latest/tutorial/poisson3DbMPI.html, the matrix is re-partitioned before being passed as argument to the solver's constructor. Is this repartition needed, or can one use the original partition as well?

if (world.size > 1) {
        prof.tic("partition");
        Partition part;

        // part(A) returns the distributed permutation matrix:
        auto P = part(*A);
        auto R = transpose(*P);

        // Reorder the matrix:
        A = product(*R, *product(*A, *P));

        // and the RHS vector:
        std::vector<double> new_rhs(R->loc_rows());
        R->move_to_backend(typename DBackend::params());
        amgcl::backend::spmv(1, *R, rhs, 0, new_rhs);
        rhs.swap(new_rhs);

        // Update the number of the local rows
        // (it may have changed as a result of permutation):
        chunk = A->loc_rows();
        prof.toc("partition");
    }

Thank you again for your support and timely response to the questions :-)

allaffa avatar Jan 25 '23 19:01 allaffa