Hi,

I run into a problem with the SparseQR solver (which does not occur with the SparseLU solver)

The problem is the following :

Code: Select all

Eigen::DenseCoeffsBase<Derived, 1>::Scalar& Eigen::DenseCoeffsBase<Derived, 1>::operator()(Eigen::DenseCoeffsBase<Derived, 1>::Index) [with Derived = Eigen::Matrix<int, -1, 1>; Eigen::DenseCoeffsBase<Derived, 1>::Scalar = int; Eigen::DenseCoeffsBase<Derived, 1>::Index = long int]: Assertion `index >= 0 && index < size()' failed.


The following assert is not verified :

Code: Select all

EIGEN_STRONG_INLINE Scalar&
    operator()(Index index)
    {
      eigen_assert(index >= 0 && index < size());
      return derived().coeffRef(index);
    }


Indeed the size of the matrix is 52*52 and index is equal to 52

The problem occurs on line 372 of SparseQR.h in the factorize function :

Code: Select all

   
       // Traverse the etree
      Index bi = nzcolR;
      for (; mark(st) != col; st = m_etree(st))
      {
        Ridx(nzcolR) = st;  // Add this row to the list,
        mark(st) = col;     // and mark this row as visited
        nzcolR++;
      }


I don't really know what I am doing wrong. I am not even sure that I am doing something wrong, since the SparseLU solver works. The only particularity is that the lhs matrix is row major. (I need it to reformulate some coefficients row-by-row, see viewtopic.php?f=74&t=122606 )

- Is it known that SparseQR may not work with row-major sparse matrix ? (I didn't found any indications in the documentation)
- The problem occurs at the third factorization, is it possible that it only occurs because my matrix is singular ?
- Is this a possible bug ?

Can you show the code where SParseQR is instantiated and called together with the exact types of the involved matrices. Thank you.

btw, it seems that you are using an old version of Eigen. Please give the 3.2.2 a try. Maybe you are just hitting a pretty old bug.

Updating to 3.2.2 didn't change the problem

It is part of a Semismooth Newton solver. This is the main code

Code: Select all

   
Eigen::SparseMatrix<double, Eigen::RowMajor> jacobian;
    setup_jacobian(x, dt, jacobian);
    Eigen::SparseQR<Eigen::SparseMatrix<double, Eigen::RowMajor>,
                    Eigen::NaturalOrdering<int>> solver;
    solver.compute(jacobian);


Changing the ordering method does not correct the problem.
setup_jacobian has two step, first I compute it :

Code: Select all

 
    list_triplet_t jacob;
    // filling jacob
    jacobian = SparseMatrix(get_neq(), get_neq());
    jacobian.setFromTriplets(jacob.begin(), jacob.end());


then I reformulate it using inner iterator :

Code: Select all

           jacobian.makeCompressed()
           // compute d and c
            for (SparseMatrix::InnerIterator it(jacobian,i); it; ++it)
            {
                it.valueRef() *= d;
                if (it.row() == it.col()) it.valueRef() += c;
            }


The solver doesn't converge even with the LU solver, so there is probably a problem somewhere in my code, it is possible that the jacobian is close to singular, leading to a problem with the QR algorithm. With the current parameters the problem only occurs at the 38th iteration. But I would still expect the solver to finish its jobs and indicate the failure when calling solver.info() not with a failed assertion. Is there a way for me to provide you the matrix so you can test ?

( unrelated note with Eigen-3.2.2 and gcc-4.8.3 I have the following non-critical warning

Code: Select all

In file included from eigen-3.2.2/Eigen/Cholesky:24:0,
                 from eigen-3.2.2/Eigen/Dense:3,
                 from common.hpp:37,
eigen-3.2.2/Eigen/src/Cholesky/LDLT.h: In member function 'void Eigen::internal::solve_retval<Eigen::LDLT<MatrixType, _UpLo>, Rhs>::evalTo(Dest&) const':
eigen-3.2.2/Eigen/src/Cholesky/LDLT.h:505:39: warning: typedef 'Scalar' locally defined but not used [-Wunused-local-typedefs]
     typedef typename LDLTType::Scalar Scalar;
                                       ^


)

using Eigen::AMDOrdering<int> as the ordering method with Eigen-3.2.2 leads to that compilation error :

Code: Select all

In file included from eigen-3.2.2/Eigen/OrderingMethods:63:0,
                 eigen-3.2.2/Eigen/Sparse:20,
                 from specmicp/src/reactmicp/common.hpp:38,
                 from specmicp/src/reactmicp/systems/diffusion/diffusion.hpp:42,
                 from specmicp/tests/reactmip/systems/diffusion/solving.cpp:4:
eigen-3.2.2/Eigen/src/OrderingMethods/Ordering.h: In instantiation of 'void Eigen::AMDOrdering<Index>::operator()(const MatrixType&, Eigen::AMDOrdering<Index>::PermutationType&) [with MatrixType = Eigen::SparseMatrix<double, 1>; Index = int; Eigen::AMDOrdering<Index>::PermutationType = Eigen::PermutationMatrix<-1>]':
eigen-3.2.2/Eigen/src/SparseQR/SparseQR.h:286:20:   required from 'void Eigen::SparseQR<MatrixType, OrderingType>::analyzePattern(const MatrixType&) [with _MatrixType = Eigen::SparseMatrix<double, 1>; _OrderingType = Eigen::AMDOrdering<int>; Eigen::SparseQR<MatrixType, OrderingType>::MatrixType = Eigen::SparseMatrix<double, 1>]'
eigen-3.2.2/Eigen/src/SparseQR/SparseQR.h:100:25:   required from 'void Eigen::SparseQR<MatrixType, OrderingType>::compute(const MatrixType&) [with _MatrixType = Eigen::SparseMatrix<double, 1>; _OrderingType = Eigen::AMDOrdering<int>; Eigen::SparseQR<MatrixType, OrderingType>::MatrixType = Eigen::SparseMatrix<double, 1>]'
eigen-3.2.2/Eigen/src/SparseQR/SparseQR.h:89:18:   required from 'Eigen::SparseQR<MatrixType, OrderingType>::SparseQR(const MatrixType&) [with _MatrixType = Eigen::SparseMatrix<double, 1>; _OrderingType = Eigen::AMDOrdering<int>; Eigen::SparseQR<MatrixType, OrderingType>::MatrixType = Eigen::SparseMatrix<double, 1>]'
specmicp/src/reactmicp/micpsolvers/micpsolver.inl:279:53:   required from 'specmicp::reactmicp::micpsolver::MiCPParabolicSolverReturnCode specmicp::reactmicp::micpsolver::MiCPParabolicSolver<program_t>::compute_search_direction(Eigen::VectorXd&, specmicp::reactmicp::ParabolicVariables&, double) [with program_t = specmicp::reactmicp::systems::DiffusionProgram; Eigen::VectorXd = Eigen::Matrix<double, -1, 1>]'
specmicp/src/reactmicp/micpsolvers/micpsolver.inl:42:47:   required from 'specmicp::reactmicp::micpsolver::MiCPParabolicSolverReturnCode specmicp::reactmicp::micpsolver::MiCPParabolicSolver<program_t>::solve_one_timestep(double, specmicp::reactmicp::ParabolicVariables&) [with program_t = specmicp::reactmicp::systems::DiffusionProgram]'
specmicp/tests/reactmip/systems/diffusion/solving.cpp:53:58:   required from here
eigen-3.2.2/Eigen/src/OrderingMethods/Ordering.h:64:51: error: no matching function for call to 'ordering_helper_at_plus_a(const Eigen::SparseMatrix<double, 1>&, Eigen::SparseMatrix<double, 0, int>&)'
       internal::ordering_helper_at_plus_a(mat,symm);
                                                   ^
eigen-3.2.2/Eigen/src/OrderingMethods/Ordering.h:64:51: note: candidate is:
eigen-3.2.2/Eigen/src/OrderingMethods/Ordering.h:26:6: note: template<class MatrixType> void Eigen::internal::ordering_helper_at_plus_a(const MatrixType&, MatrixType&)
 void ordering_helper_at_plus_a(const MatrixType& mat, MatrixType& symmat)
      ^
eigen-3.2.2/Eigen/src/OrderingMethods/Ordering.h:26:6: note:   template argument deduction/substitution failed:
eigen-3.2.2/Eigen/src/OrderingMethods/Ordering.h:64:51: note:   deduced conflicting types for parameter 'MatrixType' ('Eigen::SparseMatrix<double, 1>' and 'Eigen::SparseMatrix<double, 0, int>')
       internal::ordering_helper_at_plus_a(mat,symm);
                                                   ^


Also, the BiCGSTAB failed to solve the linear system

Yes, having access to the matrix will help to spot the issue. See http://eigen.tuxfamily.org/dox-devel/gr ... tml#title3 to export it in the matrix market format. There is no reason for SparseLU to be numerically more robust than SparseQR!

AMDOrdering is for symmetric problems. It is not compatible with QR.

Regarding BiCGSTAB, it is well known that BiCGSTAB does not converge on some problems, nonetheless, make sure you tried BiCGSTAB with the IncompleteLUT preconditioner, not the default Jacobi one.

The matrix : http://pastebin.com/6JeLCiab
The RHS : http://pastebin.com/LWxJNA19

This is not a "good" system, I still have some bugs in my code, and this is probably why the assertion is raised, but I would expect that the solver fails nicely and report the error with the solver.info() methods. The fact that SuperLU and BiCGSTAB manage to return an answer is not very relevant, since the answer is garbage, but at least they don't raise an error that terminates the program.

But I think this has something to do with the "row-major" thingy, since transforming the matrix to a column major leads to an answer (still garbage but the assertion is not raised) :

Code: Select all

    Eigen::SparseMatrix<double, Eigen::RowMajor> jacobian;
    setup_jacobian(x, dt, jacobian);
    Eigen::SparseMatrix<double, Eigen::ColMajor> new_jacobian(jacobian);
    Eigen::SparseQR<Eigen::SparseMatrix<double, Eigen::ColMajor>,
           Eigen::NaturalOrdering<int>> solver(new_jacobian);


The backtrace is :

Code: Select all

0   raise   /lib64/libc.so.6      0x7ffff703a385   
1   abort   /lib64/libc.so.6      0x7ffff703b808   
2   __assert_fail_base   /lib64/libc.so.6      0x7ffff7033452   
3   __assert_fail   /lib64/libc.so.6      0x7ffff7033502   
4   Eigen::DenseCoeffsBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>, 1>::operator()   DenseCoeffsBase.h   394   0x66c058   
5   Eigen::SparseQR<Eigen::SparseMatrix<double, 1, int>, Eigen::COLAMDOrdering<int> >::factorize   SparseQR.h   393   0x668ef6   
6   Eigen::SparseQR<Eigen::SparseMatrix<double, 1, int>, Eigen::COLAMDOrdering<int> >::compute   SparseQR.h   101   0x666ed4   
7   Eigen::SparseQR<Eigen::SparseMatrix<double, 1, int>, Eigen::COLAMDOrdering<int> >::SparseQR   SparseQR.h   89   0x664b0f   
8   specmicp::reactmicp::micpsolver::MiCPParabolicSolver<specmicp::reactmicp::systems::DiffusionProgram>::compute_search_direction   micpsolver.inl   284   0x6628af   
9   specmicp::reactmicp::micpsolver::MiCPParabolicSolver<specmicp::reactmicp::systems::DiffusionProgram>::solve_one_timestep   micpsolver.inl   44   0x66138d   
10   ____C_A_T_C_H____T_E_S_T____17   solving.cpp   54   0x65b834   
11   Catch::FreeFunctionTestCase::invoke   catch_test_case_registry_impl.hpp   88   0x62e78c   
12   Catch::TestCase::invoke   catch_test_case_info.hpp   168   0x61ce0b   
13   Catch::RunContext::runCurrentTest   catch_runner_impl.hpp   238   0x62c7e5   
14   Catch::RunContext::runTest   catch_runner_impl.hpp   107   0x62b7b7   
15   Catch::Runner::runTests   catch_runner.hpp   59   0x62d16b   
16   Catch::Session::run   catch_runner.hpp   179   0x62df9c   
17   Catch::Session::run   catch_runner.hpp   162   0x62de4c   
18   main   catch_default_main.hpp   15   0x621904   


In the factorize method :

Code: Select all

      // Traverse the etree
      Index bi = nzcolR;
      for (; mark(st) != col; st = m_etree(st))
      {
        Ridx(nzcolR) = st;  // Add this row to the list,
        mark(st) = col;     // and mark this row as visited
        nzcolR++;
      }


When the assertion is raised 'st' is equal to 52 (hence the error since mark is a Vector of size 52, m_etree is a vector of size 52 where m_etree[i] = i+1 (so st= m_etree[51] = 52 ))
This is the values of the local variables of the method : http://pastebin.com/LtFzrNNU

I tried to go deeper in the code, but I don't really know these algorithms so I am a bit lost...

alright, this is fixed in the 3.2 and default branch.

What do you mean by garbage? For me your example is working ok with both SparseQ, SparseLU and BiCGSTAB with a relative error of 1e-7. (your matrix is highly singular, so 1e-7 is pretty good)

Great thank you very much

I still have a question, your answer is similar to mine, you copy the matrix in a column-major matrix before you compute the column elimination tree. How expensive is it ?
This problem wasn't happening with every row-major matrix I encounter, so it may be too much...

and when I said garbage, I was just meaning, not at all the solution I expect (but that's my problem)

All the computation and analysis have to be performed column-wise, there is no alternative. The best would be to assemble a column-major matrix at the first place.

ok, so the fact that it sometimes worked was purely luck...
I have to reformulate the Jacobian row-by-row, so using a row-major seems required

Thank you for Eigen and your helpful support