Convergence issue with gmres + Eigen backend

Hi all,

I have been playing with AMGCL + Eigen backend for complex matrices and had issues getting convergence, while everything was fine for the builtin backend (or Eigen adapter w/ builtin backend). TL;DR is : I think the inner product of the eigen backend is at fault:

template < class V1, class V2 >
struct inner_product_impl<
    V1, V2, 
    typename std::enable_if<
        is_eigen_type<V1>::value &&
        is_eigen_type<V2>::value
        >::type
    >   
{
    typedef typename value_type<V1>::type real;
    static real get(const V1 &x, const V2 &y) 
    {   
        //return y.dot(x);
        return x.dot(y);
    }   
};

With the original x.dot(y) gmres stagnates, while with y.dot(x) I see the same results as with the builtin backend. Any idea what could cause this assymetric behavior? (FYI the conjugation is done on the second parameter in struct inner_product_impl< std::complex<T> > while eigen will conjugate the first one...)

Here is the MWE I have (run with ./a.out mat.txt rhs.txt), you can switch between the backends by switching the USE_EIGEN_BACKEND macro rhs.txt mat.txt

I can create PR with the 'fix' but I would like to have a second opinion on why this should happen...

#include <iostream>
#include <complex>

#include <Eigen/Sparse>
#include <Eigen/Dense>

#include <unsupported/Eigen/SparseExtra> // For reading MatrixMarket files

#define USE_EIGEN_BACKEND 1

#include <amgcl/value_type/complex.hpp>
#include <amgcl/adapter/complex.hpp>
#include <amgcl/adapter/eigen.hpp>
#if USE_EIGEN_BACKEND
    #include <amgcl/backend/eigen.hpp>
#else
    AMGCL_USE_EIGEN_VECTORS_WITH_BUILTIN_BACKEND()
#endif
#include <amgcl/backend/builtin.hpp>
#include <amgcl/adapter/crs_tuple.hpp>
#include <amgcl/make_solver.hpp>
#include <amgcl/solver/bicgstab.hpp>
#include <amgcl/solver/gmres.hpp>
#include <amgcl/amg.hpp>
#include <amgcl/coarsening/smoothed_aggregation.hpp>
#include <amgcl/relaxation/spai0.hpp>
#include <amgcl/io/mm.hpp>

using complex_t = std::complex<double>;

int main(int argc, char *argv[]) {
    if (argc < 3) {
        std::cerr << "Usage: " << argv[0] << " <matrix.mm>" << std::endl;
        return 1;
    }

#if USE_EIGEN_BACKEND
    {
        Eigen::SparseMatrix<complex_t, Eigen::RowMajor> A;
        Eigen::loadMarket(A, argv[1]);

        Eigen::VectorXcd f = Eigen::VectorXcd::Constant(A.rows(), 1.0);
        Eigen::loadMarketVector(f, argv[2]);

        Eigen::VectorXcd x = Eigen::VectorXcd::Zero(A.rows());

        typedef amgcl::make_solver<
            amgcl::amg<
            amgcl::backend::eigen<complex_t>,
            amgcl::coarsening::smoothed_aggregation,
            amgcl::relaxation::spai0
                >,
            amgcl::solver::gmres<amgcl::backend::eigen<complex_t> >
                > Solver;

        Solver solve(A);
        std::cout << solve << std::endl;

        int    iters;
        double error;
        std::tie(iters, error) = solve(f, x);

        std::cout << iters << " " << error << std::endl;
    }
#else
    {
        Eigen::SparseMatrix<complex_t, Eigen::RowMajor> A;
        Eigen::loadMarket(A, argv[1]);

        Eigen::VectorXcd f = Eigen::VectorXcd::Constant(A.rows(), 1.0);
        Eigen::loadMarketVector(f, argv[2]);

        Eigen::VectorXcd x = Eigen::VectorXcd::Zero(A.rows());

        // Setup the solver:
        typedef amgcl::make_solver<
            amgcl::amg<
            amgcl::backend::builtin<complex_t>,
            amgcl::coarsening::smoothed_aggregation,
            amgcl::relaxation::spai0
                >,
            amgcl::solver::gmres<amgcl::backend::builtin<complex_t> >
                > Solver;

        Solver solve(A);
        std::cout << solve << std::endl;

        // Solve the system for the given RHS:
        int    iters;
        double error;
        std::tie(iters, error) = solve(f, x);

        std::cout << iters << " " << error << std::endl;
    }
#endif

    {
        ptrdiff_t rows, cols;
        std::vector<ptrdiff_t> ptr, col;
        std::vector<complex_t> val, f;

        std::tie(rows, cols) = amgcl::io::mm_reader(argv[1])(ptr, col, val);
        std::tie(rows, cols) = amgcl::io::mm_reader(argv[2])(f);
        std::vector<complex_t> x(rows, 0.0);

        // Setup the solver:
        typedef amgcl::make_solver<
            amgcl::amg<
            amgcl::backend::builtin<complex_t>,
            amgcl::coarsening::smoothed_aggregation,
            amgcl::relaxation::spai0
                >,
            amgcl::solver::gmres<amgcl::backend::builtin<complex_t> >
                > Solver;

        auto A = std::tie(rows, ptr, col, val);

        Solver solve(A);
        std::cout << solve << std::endl;

        // Solve the system for the given RHS:
        int    iters;
        double error;
        std::tie(iters, error) = solve(f, x);

        std::cout << iters << " " << error << std::endl;
    }
    return EXIT_SUCCESS;
}

ddemidov / amgcl

Convergence issue with gmres + Eigen backend #258