Hello

I get a "munmap_chunk: invalid pointer" error leading to a SIGABRT in my code when using the SparseLU solver. I suspect that I'm feeding the solver something fishy. However the error arrive after 1h+ of computation in release mode and most values are optimized away. So I don't really know what's going wrong. However I'd prefer if the sparseLU solver raises an error in such cases rather than abort the program.

Unfortunately, the long computation time means that the error may be hard to catch even if it's quite reproducible (the errors happened with Eigen 3.2.8 on 3 different machines). I run a debug session, and an assertion was raised instead, but I forgot to copy it :/. It was saying that something was wrong with the factorization.

The code (huge and ugly) is here :
https://bitbucket.org/specmicp/specmicp

The trouble occurs in the function "restart_timestep" (line 111) of this file : https://bitbucket.org/specmicp/specmicp ... ew-default

There is a small wrapper around the Eigen solvers that allows me to chose the one I want at runtime. The error does not occur with the SparseQR solver. The corresponding code is here : https://bitbucket.org/specmicp/specmicp ... ?at=master

I ran valgrind to get more information, the output is pasted below. Unfortunately I'm not sure what's it's telling me.

My questions are :
- Is there anything I'm doing obviously wrong ?
- IS there anything more I can do to track down the error ?
- Does the SparseLU solver assumes that the input is "correct" ? What is the definition of correct in this case ?
- Can the SparseLU solver be changed to catch the problem and abort the factorization rather the cause a memory error ?

Code: Select all

==3628== Invalid write of size 8
==3628==    at 0x4062F3F: copyCoeff<Eigen::Block<const Eigen::Block<Eigen::DiagonalProduct<Eigen::CwiseUnaryOp<Eigen::internal::scalar_opposite_op<double>, const Eigen::Matrix<double, -1, 1> >, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>, -1, 1, true>, 1, 1, false> > (DenseCoeffsBase.h:495)
==3628==    by 0x4062F3F: run (Assign.h:180)
==3628==    by 0x4062F3F: run (Assign.h:313)
==3628==    by 0x4062F3F: lazyAssign<Eigen::Block<const Eigen::Block<Eigen::DiagonalProduct<Eigen::CwiseUnaryOp<Eigen::internal::scalar_opposite_op<double>, const Eigen::Matrix<double, -1, 1> >, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>, -1, 1, true>, 1, 1, false> > (Assign.h:506)
==3628==    by 0x4062F3F: run (Assign.h:527)
==3628==    by 0x4062F3F: operator=<Eigen::Block<const Eigen::Block<Eigen::DiagonalProduct<Eigen::CwiseUnaryOp<Eigen::internal::scalar_opposite_op<double>, const Eigen::Matrix<double, -1, 1> >, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>, -1, 1, true>, 1, 1, false> > (Assign.h:571)
==3628==    by 0x4062F3F: evalTo<Eigen::Block<Eigen::Matrix<double, -1, 1>, -1, 1, true> > (PermutationMatrix.h:621)
==3628==    by 0x4062F3F: evalTo<Eigen::Block<Eigen::Matrix<double, -1, 1>, -1, 1, true> > (ReturnByValue.h:61)
==3628==    by 0x4062F3F: evalTo<Eigen::Block<Eigen::Matrix<double, -1, 1>, -1, 1, true>, Eigen::ReturnByValue<Eigen::internal::permut_matrix_product_retval<Eigen::PermutationMatrix<-1>, Eigen::Block<Eigen::DiagonalProduct<Eigen::CwiseUnaryOp<Eigen::internal::scalar_opposite_op<double>, const Eigen::Matrix<double, -1, 1> >, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>, -1, 1, true>, 1, false> > > (Assign.h:529)
==3628==    by 0x4062F3F: operator=<Eigen::internal::permut_matrix_product_retval<Eigen::PermutationMatrix<-1>, Eigen::Block<Eigen::DiagonalProduct<Eigen::CwiseUnaryOp<Eigen::internal::scalar_opposite_op<double>, const Eigen::Matrix<double, -1, 1> >, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>, -1, 1, true>, 1, false> > (Assign.h:585)
==3628==    by 0x4062F3F: _solve<Eigen::DiagonalProduct<Eigen::CwiseUnaryOp<Eigen::internal::scalar_opposite_op<double>, const Eigen::Matrix<double, -1, 1> >, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>, Eigen::Matrix<double, -1, 1> > (SparseLU.h:235)
==3628==    by 0x4062F3F: evalTo<Eigen::Matrix<double, -1, 1> > (SparseLU.h:786)
==3628==    by 0x4062F3F: evalTo<Eigen::Matrix<double, -1, 1> > (Solve.h:51)
==3628==    by 0x4062F3F: evalTo<Eigen::Matrix<double, -1, 1> > (ReturnByValue.h:61)
==3628==    by 0x4062F3F: evalTo<Eigen::Matrix<double, -1, 1>, Eigen::ReturnByValue<Eigen::internal::solve_retval_base<Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int> >, Eigen::DiagonalProduct<Eigen::CwiseUnaryOp<Eigen::internal::scalar_opposite_op<double>, const Eigen::Matrix<double, -1, 1> >, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1> > > > (Assign.h:529)
==3628==    by 0x4062F3F: operator=<Eigen::internal::solve_retval_base<Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int> >, Eigen::DiagonalProduct<Eigen::CwiseUnaryOp<Eigen::internal::scalar_opposite_op<double>, const Eigen::Matrix<double, -1, 1> >, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1> > > (Assign.h:585)
==3628==    by 0x4062F3F: operator=<Eigen::internal::solve_retval_base<Eigen::SparseLU<Eigen::SparseMatrix<double, 0, int>, Eigen::COLAMDOrdering<int> >, Eigen::DiagonalProduct<Eigen::CwiseUnaryOp<Eigen::internal::scalar_opposite_op<double>, const Eigen::Matrix<double, -1, 1> >, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1> > > (PlainObjectBase.h:421)
==3628==    by 0x4062F3F: operator=<Eigen::internal::solve_retval_base<Eigen::SparseLU<Eigen::SparseMatrix<double, 0, int>, Eigen::COLAMDOrdering<int> >, Eigen::DiagonalProduct<Eigen::CwiseUnaryOp<Eigen::internal::scalar_opposite_op<double>, const Eigen::Matrix<double, -1, 1> >, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1> > > (Matrix.h:190)
==3628==    by 0x4062F3F: specmicp::sparse_solvers::SparseSolverLU<Eigen::SparseMatrix<double, 0, int>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >::solve_scaling(Eigen::Matrix<double, -1, 1, 0, -1, 1> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1>&) (sparse_lu.hpp:93)
==3628==    by 0x409BABC: specmicp::dfpmsolver::ParabolicDriver<specmicp::reactmicp::systems::unsaturated::PressureEquation>::restart_timestep(Eigen::Matrix<double, -1, 1, 0, -1, 1>&) (parabolic_driver.inl:150)
==3628==    by 0x40AC81F: specmicp::reactmicp::systems::unsaturated::EquationTask<specmicp::reactmicp::systems::unsaturated::PressureTask, specmicp::reactmicp::systems::unsaturated::PressureTaskTraits>::restart_timestep(specmicp::reactmicp::systems::unsaturated::UnsaturatedVariables*) (transport_stagger.cpp:432)
==3628==    by 0x40A57C0: specmicp::reactmicp::systems::unsaturated::UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::restart_timestep(specmicp::reactmicp::systems::unsaturated::UnsaturatedVariables*) [clone ._omp_fn.4] (transport_stagger.cpp:1092)
==3628==    by 0x68FBCFE: GOMP_parallel (in /usr/lib64/gcc/x86_64-pc-linux-gnu/5.3.0/libgomp.so.1.0.0)
==3628==    by 0x40A541E: specmicp::reactmicp::systems::unsaturated::UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::restart_timestep(specmicp::reactmicp::systems::unsaturated::UnsaturatedVariables*) (transport_stagger.cpp:1085)
==3628==    by 0x407AAA6: specmicp::reactmicp::solver::ReactiveTransportSolver::one_iteration(specmicp::reactmicp::solver::VariablesBase*, specmicp::reactmicp::solver::internal::ReactiveTransportResiduals&) (reactive_transport_solver.cpp:193)
==3628==    by 0x407AE4D: specmicp::reactmicp::solver::ReactiveTransportSolver::solve_timestep(double, specmicp::reactmicp::solver::VariablesBase*) (reactive_transport_solver.cpp:120)
==3628==    by 0x407BB03: specmicp::reactmicp::solver::ReactiveTransportRunner::ReactiveTransportRunnerImpl::run_until(double, std::shared_ptr<specmicp::reactmicp::solver::VariablesBase>) (runner.cpp:167)
==3628==    by 0x407C037: specmicp::reactmicp::solver::ReactiveTransportRunner::run_until(double, std::shared_ptr<specmicp::reactmicp::solver::VariablesBase>) (runner.cpp:120)
==3628==    by 0x408D50: main (unsaturated.cpp:289)
==3628==  Address 0x7b6fb68 is 8 bytes before a block of size 400 alloc'd
==3628==    at 0x4C29FAF: malloc (vg_replace_malloc.c:299)
==3628==    by 0x40421A0: aligned_malloc (Memory.h:220)
==3628==    by 0x40421A0: conditional_aligned_malloc<true> (Memory.h:301)
==3628==    by 0x40421A0: double* Eigen::internal::conditional_aligned_new_auto<double, true>(unsigned long) (Memory.h:423)
==3628==    by 0x409B71B: DenseStorage (DenseStorage.h:386)
==3628==    by 0x409B71B: PlainObjectBase (PlainObjectBase.h:471)
==3628==    by 0x409B71B: Matrix (Matrix.h:236)
==3628==    by 0x409B71B: specmicp::dfpmsolver::ParabolicDriver<specmicp::reactmicp::systems::unsaturated::PressureEquation>::restart_timestep(Eigen::Matrix<double, -1, 1, 0, -1, 1>&) (parabolic_driver.inl:117)
==3628==    by 0x40AC81F: specmicp::reactmicp::systems::unsaturated::EquationTask<specmicp::reactmicp::systems::unsaturated::PressureTask, specmicp::reactmicp::systems::unsaturated::PressureTaskTraits>::restart_timestep(specmicp::reactmicp::systems::unsaturated::UnsaturatedVariables*) (transport_stagger.cpp:432)
==3628==    by 0x40A57C0: specmicp::reactmicp::systems::unsaturated::UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::restart_timestep(specmicp::reactmicp::systems::unsaturated::UnsaturatedVariables*) [clone ._omp_fn.4] (transport_stagger.cpp:1092)
==3628==    by 0x68FBCFE: GOMP_parallel (in /usr/lib64/gcc/x86_64-pc-linux-gnu/5.3.0/libgomp.so.1.0.0)
==3628==    by 0x40A541E: specmicp::reactmicp::systems::unsaturated::UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::restart_timestep(specmicp::reactmicp::systems::unsaturated::UnsaturatedVariables*) (transport_stagger.cpp:1085)
==3628==    by 0x407AAA6: specmicp::reactmicp::solver::ReactiveTransportSolver::one_iteration(specmicp::reactmicp::solver::VariablesBase*, specmicp::reactmicp::solver::internal::ReactiveTransportResiduals&) (reactive_transport_solver.cpp:193)
==3628==    by 0x407AE4D: specmicp::reactmicp::solver::ReactiveTransportSolver::solve_timestep(double, specmicp::reactmicp::solver::VariablesBase*) (reactive_transport_solver.cpp:120)
==3628==    by 0x407BB03: specmicp::reactmicp::solver::ReactiveTransportRunner::ReactiveTransportRunnerImpl::run_until(double, std::shared_ptr<specmicp::reactmicp::solver::VariablesBase>) (runner.cpp:167)
==3628==    by 0x407C037: specmicp::reactmicp::solver::ReactiveTransportRunner::run_until(double, std::shared_ptr<specmicp::reactmicp::solver::VariablesBase>) (runner.cpp:120)
==3628==    by 0x408D50: main (unsaturated.cpp:289)
==3628==
terminate called after throwing an instance of 'std::bad_alloc'
  what():  std::bad_alloc
==3628==
==3628== Process terminating with default action of signal 6 (SIGABRT)
==3628==    at 0x6F78538: raise (in /lib64/libc-2.22.so)
==3628==    by 0x6F79989: abort (in /lib64/libc-2.22.so)
==3628==    by 0x65FBB4C: __gnu_cxx::__verbose_terminate_handler() (in /usr/lib64/gcc/x86_64-pc-linux-gnu/5.3.0/libstdc++.so.6.0.21)
==3628==    by 0x65F9975: ??? (in /usr/lib64/gcc/x86_64-pc-linux-gnu/5.3.0/libstdc++.so.6.0.21)
==3628==    by 0x65F8968: ??? (in /usr/lib64/gcc/x86_64-pc-linux-gnu/5.3.0/libstdc++.so.6.0.21)
==3628==    by 0x65F92C4: __gxx_personality_v0 (in /usr/lib64/gcc/x86_64-pc-linux-gnu/5.3.0/libstdc++.so.6.0.21)
==3628==    by 0x6B21D22: ??? (in /lib64/libgcc_s.so.1)
==3628==    by 0x6B22238: _Unwind_Resume (in /lib64/libgcc_s.so.1)
==3628==    by 0x409BD10: aligned_free (Memory.h:243)
==3628==    by 0x409BD10: conditional_aligned_free<true> (Memory.h:317)
==3628==    by 0x409BD10: conditional_aligned_delete_auto<double, true> (Memory.h:445)
==3628==    by 0x409BD10: ~DenseStorage (DenseStorage.h:403)
==3628==    by 0x409BD10: ~PlainObjectBase (PlainObjectBase.h:88)
==3628==    by 0x409BD10: ~Matrix (Matrix.h:127)
==3628==    by 0x409BD10: specmicp::dfpmsolver::ParabolicDriver<specmicp::reactmicp::systems::unsaturated::PressureEquation>::restart_timestep(Eigen::Matrix<double, -1, 1, 0, -1, 1>&) (parabolic_driver.inl:117)
==3628==    by 0x40AC81F: specmicp::reactmicp::systems::unsaturated::EquationTask<specmicp::reactmicp::systems::unsaturated::PressureTask, specmicp::reactmicp::systems::unsaturated::PressureTaskTraits>::restart_timestep(specmicp::reactmicp::systems::unsaturated::UnsaturatedVariables*) (transport_stagger.cpp:432)
==3628==    by 0x40A57C0: specmicp::reactmicp::systems::unsaturated::UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::restart_timestep(specmicp::reactmicp::systems::unsaturated::UnsaturatedVariables*) [clone ._omp_fn.4] (transport_stagger.cpp:1092)
==3628==    by 0x68FBCFE: GOMP_parallel (in /usr/lib64/gcc/x86_64-pc-linux-gnu/5.3.0/libgomp.so.1.0.0)
==3628==
==3628== HEAP SUMMARY:
==3628==     in use at exit: 1,845,215 bytes in 13,247 blocks
==3628==   total heap usage: 333,868,041 allocs, 333,854,794 frees, 186,561,370,594 bytes allocated
==3628==
==3628== LEAK SUMMARY:
==3628==    definitely lost: 0 bytes in 0 blocks
==3628==    indirectly lost: 0 bytes in 0 blocks
==3628==      possibly lost: 1,048 bytes in 4 blocks
==3628==    still reachable: 1,844,167 bytes in 13,243 blocks
==3628==         suppressed: 0 bytes in 0 blocks
==3628== Rerun with --leak-check=full to see details of leaked memory
==3628==
==3628== For counts of detected and suppressed errors, rerun with: -v
==3628== ERROR SUMMARY: 49 errors from 1 contexts (suppressed: 0 from 0)
zsh: killed     OMP_NUM_THREADS=4 valgrind reactmicp_unsaturated -i acc_carbo.yml -w /tmp


oh, the first error is happening in the solve step and not during the factorization. Please paste the assertion with the associated back-trace. Also make sure the factorization succeeded before calling solve, like if(lu.info()==Eigen::Succeed) x=lu.solve(b); else ....

ggael wrote:Also make sure the factorization succeeded before calling solve, like if(lu.info()==Eigen::Succeed) x=lu.solve(b); else ....

I thought it was what I was doing.... but clearly not ! With this check I'm able to stop after a failed factorization and avoid trying to solve a badly formed system. In my case I guess that I run into some problems with numerical precision. The SparseQR solver is able to overcome them, not the SparseLU solver, but I guess it can be expected. Thanks. I'm going to inspect more closely the problem.

I'll post the exact assertion that was raised when I get it again, but it may take a while for the debug version to get to the point where the problems tart to arrive.

This is the assertion, which makes sense now :

Code: Select all

reactmicp_unsaturated: /usr/include/eigen3/Eigen/src/SparseLU/SparseLU.h:180:  const Eigen::internal::solve_retval<Eigen::SparseLU<_MatrixType, _OrderingType>, Rhs> Eigen::SparseLU<_MatrixType, _OrderingType>::solve(const  Eigen::MatrixBase<OtherDerived>&) const [with Rhs = Eigen::DiagonalProduct<Eigen::CwiseUnaryOp<Eigen::internal::scalar_opposite_op<double>, const Eigen::Matrix<double, -1, 1> >, Eigen::DiagonalWrapper<const
Eigen::Matrix<double, -1, 1> >, 1>; _MatrixType = Eigen::SparseMatrix<double, 0, int>; _OrderingType = Eigen::COLAMDOrdering<int>]: Assertion `m_factorizationIsOk && "SparseLU is not initialized."' failed.


And the associated backtrace :

Code: Select all

#0  0x00007f3c8a907348 in raise () from /lib64/libc.so.6
#1  0x00007f3c8a9087ca in abort () from /lib64/libc.so.6
#2  0x00007f3c8a9002f7 in __assert_fail_base () from /lib64/libc.so.6
#3  0x00007f3c8a9003a2 in __assert_fail () from /lib64/libc.so.6
#4  0x00000000005459b4 in Eigen::SparseLU<Eigen::SparseMatrix<double, 0, int>, Eigen::COLAMDOrdering<int> >::solve<Eigen::DiagonalProduct<Eigen::CwiseUnaryOp<Eigen::internal::scalar_opposite_op<double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const>, Eigen::DiagonalWrapper<Eigen::Matrix<double, -1, 1, 0, -1, 1> const>, 1> > (this=0x289f188, B=...) at /usr/include/eigen3/Eigen/src/SparseLU/SparseLU.h:180
#5  0x0000000000543314 in specmicp::sparse_solvers::SparseSolverLU<Eigen::SparseMatrix<double, 0, int>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >::solve_scaling (this=0x289f180, residuals=...,
    scaling=..., solution=...) at /home/manawy/princeton/progra/specmicp/src/specmicp_common/sparse_solvers/sparse_lu.hpp:93
#6  0x00000000005b401a in specmicp::dfpmsolver::ParabolicDriver<specmicp::reactmicp::systems::unsaturated::PressureEquation>::restart_timestep (this=0x2916908, displacement=...)
    at /home/manawy/princeton/progra/specmicp/src/dfpm/solver/parabolic_driver.inl:161
#7  0x000000000051c9f0 in specmicp::reactmicp::systems::unsaturated::EquationTask<specmicp::reactmicp::systems::unsaturated::PressureTask, specmicp::reactmicp::systems::unsaturated::PressureTaskTraits>::restart_timestep (
    this=0x29168c0, vars=0x29a2d50) at /home/manawy/princeton/progra/specmicp/src/reactmicp/systems/unsaturated/transport_stagger.cpp:446
#8  0x0000000000514011 in specmicp::reactmicp::systems::unsaturated::UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::restart_timestep (this=0x285fa80, vars=0x29a2d50)
    at /home/manawy/princeton/progra/specmicp/src/reactmicp/systems/unsaturated/transport_stagger.cpp:1108
#9  0x000000000051275b in specmicp::reactmicp::systems::unsaturated::UnsaturatedTransportStagger::restart_timestep (this=0x29897e0, var=0x29a2d50)
    at /home/manawy/princeton/progra/specmicp/src/reactmicp/systems/unsaturated/transport_stagger.cpp:226
#10 0x000000000052a6d5 in specmicp::reactmicp::solver::ReactiveTransportSolver::one_iteration (this=0x7ffcafa8be20, variables=0x29a2d50, residuals=...)
    at /home/manawy/princeton/progra/specmicp/src/reactmicp/solver/reactive_transport_solver.cpp:199
#11 0x000000000052a222 in specmicp::reactmicp::solver::ReactiveTransportSolver::solve_timestep (this=0x7ffcafa8be20, timestep=0.10000000000000001, variables=0x29a2d50)
    at /home/manawy/princeton/progra/specmicp/src/reactmicp/solver/reactive_transport_solver.cpp:120
#12 0x000000000050485f in specmicp::reactmicp::solver::ReactiveTransportRunner::ReactiveTransportRunnerImpl::run_until (this=0x291c7b0, target=2592000, variables=std::shared_ptr (count 4, weak 0) 0x29a2d50)
    at /home/manawy/princeton/progra/specmicp/src/reactmicp/solver/runner.cpp:167
#13 0x00000000005044c3 in specmicp::reactmicp::solver::ReactiveTransportRunner::run_until (this=0x7ffcafa8bd00, target=2592000, variables=std::shared_ptr (count 4, weak 0) 0x29a2d50)
    at /home/manawy/princeton/progra/specmicp/src/reactmicp/solver/runner.cpp:120
#14 0x00000000004da70e in main (argc=5, argv=0x7ffcafa8cbe8) at /home/manawy/princeton/progra/specmicp/src/bin/reactmicp/unsaturated.cpp:289
(gdb) frame 4
#4  0x00000000005459b4 in Eigen::SparseLU<Eigen::SparseMatrix<double, 0, int>, Eigen::COLAMDOrdering<int> >::solve<Eigen::DiagonalProduct<Eigen::CwiseUnaryOp<Eigen::internal::scalar_opposite_op<double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const>, Eigen::DiagonalWrapper<Eigen::Matrix<double, -1, 1, 0, -1, 1> const>, 1> > (this=0x289f188, B=...) at /usr/include/eigen3/Eigen/src/SparseLU/SparseLU.h:180
180           eigen_assert(m_factorizationIsOk && "SparseLU is not initialized.");
(gdb) frame 5
#5  0x0000000000543314 in specmicp::sparse_solvers::SparseSolverLU<Eigen::SparseMatrix<double, 0, int>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >::solve_scaling (this=0x289f180, residuals=...,
    scaling=..., solution=...) at /home/manawy/princeton/progra/specmicp/src/specmicp_common/sparse_solvers/sparse_lu.hpp:93
93              solution = m_solver.solve(scaling.asDiagonal()*(-residuals));


It was my fault for not checking the factorization step. Thank you for your answer.