How to keep grid hierarchy once and modify solver parameter when solve many right hand side

[qzxuhui]

Sir, I want to solve many right-hand sides. I want the grid hierarchy only to be created once. What I want can be expressed using follow pseudocode.

SparseMatrix A;
Vector b1,b2,b3;
Solver solver;
solver.compute(A);
solver.setTolerence(1e-5);
Vector x1 = solver.solve(b1);
solver.setTolerence(1e-8);
Vector x2 = solver.solve(b2);
solver.setTolerence(1e-10);
Vector x3 = solver.solve(b3);

I have tried something like follow but seems not to work.

        typedef amgcl::make_solver<
            amgcl::amg<Backend, amgcl::coarsening::smoothed_aggregation, amgcl::relaxation::damped_jacobi>,
            amgcl::solver::cg<Backend> >
            Solver;

        Solver::params prm;
        prm.solver.maxiter = 500;
        prm.solver.tol = 1e-14;

        // we first set tol to 1e-14
        Solver solve(amgcl::adapter::block_matrix<mat_type>(A), prm);

        // now I want to 1e-10
        solve.prm.solver.tol = 1e-10;

         solve(x,f);

How can I implement that in amgcl ? Thanks for your time.

[ddemidov]

That is not directly possible. However, you could split the make_solver class into the solver and the preconditioner instances. The major time during the setup is spent constructing the preconditioner, so recreating the solver should not be a problem:

amgcl::amg<Backend, amgcl::coarsening::smoothed_aggregation, amgcl::relaxation::damped_jacobi> P(A, p_prm);
amgcl::solver::cg<Backend> S(n, s_prm);
std::tie(iters, error) = S(A, P, f, x);

In fact, the solver params are a public member of the solver class, so you could update the tolerance directly before each solve: https://github.com/ddemidov/amgcl/blob/e82f164864d62648cda26b3f0e13a80b1edd682a/amgcl/solver/cg.hpp#L213-L214

[qzxuhui]

Sir, I tried your suggestions. But it leads to Segmentation fault in my case. I do not know why this happen.

Could you please help me?

My code shown as follows, it is modified from https://github.com/ddemidov/amgcl/issues/175

void usingAMGCL()
{
    amgcl::profiler<> prof;

    // Read sparse matrix from MatrixMarket format.
    // In general this should come pre-assembled.
    const std::string A_file = "matrix.txt";
    SparseMatrix<double, RowMajor> A = readSparseMatrixFromFile(A_file);

    const std::string f_file = "prd_15.txt";
    Eigen::VectorXd f = readVectorFromFile(f_file);

    // Zero initial approximation:
    Eigen::VectorXd x = Eigen::VectorXd::Constant(A.rows(), 0.0);

    // Setup the solver:
    typedef amgcl::static_matrix<double, 3, 3> mat_type;
    typedef amgcl::static_matrix<double, 3, 1> rhs_type;
    typedef amgcl::backend::builtin<mat_type> Backend;

    typedef amgcl::make_solver<
        amgcl::amg<
            Backend,
            amgcl::coarsening::smoothed_aggregation,
            amgcl::relaxation::damped_jacobi>,
        amgcl::solver::cg<Backend>>
        AMGSolver;

    AMGSolver::params prm;
    prm.solver.maxiter = 500;
    prm.solver.tol = 1e-14;

    // first try -----------------------------------

    prof.tic("setup");
    AMGSolver solve(amgcl::adapter::block_matrix<mat_type>(A), prm);
    prof.toc("setup");
    std::cout << solve << std::endl;

    // Solve the system for the given RHS:
    int iters;
    double error;

    prof.tic("first solve");
    auto f_ptr = reinterpret_cast<rhs_type *>(f.data());
    auto x_ptr = reinterpret_cast<rhs_type *>(x.data());
    std::tie(iters, error) = solve(
        amgcl::make_iterator_range(f_ptr, f_ptr + A.rows() / 3),
        amgcl::make_iterator_range(x_ptr, x_ptr + A.rows() / 3));
    prof.toc("first solve");

    std::cout << "first try success..." << std::endl;
    std::cout << iters << " " << error << std::endl;

    // second try -----------------------------------
    amgcl::amg<Backend, amgcl::coarsening::smoothed_aggregation, amgcl::relaxation::damped_jacobi> P(amgcl::adapter::block_matrix<mat_type>(A));
    amgcl::solver::cg<Backend> S(A.rows());
    for (int i = 0;i<A.rows();i++)
    {
        x(i) = 0;
    }
    std::tie(iters, error) = S(amgcl::adapter::block_matrix<mat_type>(A),
                               P,
                               amgcl::make_iterator_range(f_ptr, f_ptr + A.rows() / 3),
                               amgcl::make_iterator_range(x_ptr, x_ptr + A.rows() / 3));

    std::cout << "second try success..." << std::endl;
    std::cout << iters << " " << error << std::endl;

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

The first try can solve success, as you told me yesterday. But the second try is failed. The log is shown as follows.

Solver
======
Type:             CG
Unknowns:         73284
Memory footprint: 6.71 M

Preconditioner
==============
Number of levels:    3
Operator complexity: 1.23
Grid complexity:     1.07
Memory footprint:    146.40 M

level     unknowns       nonzeros      memory
---------------------------------------------
    0        73284         998886    124.28 M (81.41%)
    1         5063         217965     19.94 M (17.76%)
    2          242          10102      2.18 M ( 0.82%)

first try success...
222 8.59069e-15
./run.sh: line 6:  5770 Segmentation fault

[ddemidov]

I think the problem is with this line:

S(A.rows());

The block-valued matrix has A.rows() / 3 rows. This version works: https://gist.github.com/ddemidov/5ebb67a5a79e370f53f325dd4655a031

[qzxuhui]

You are so kind sir. I am really sorry for my stupid! Really thanks for your time, sir. Everything is working well now~