Undefined behavior in ODE solver

CRAN is telling me that there is some undefined behavior in the new StanHeaders that I need help tracking down. The code in question comes from

https://cran.r-project.org/web/packages/StanHeaders/vignettes/stanmath.html

specifically this (trimmed) C++ chunk with a simple harmonic oscillator:

// [[Rcpp::depends(BH)]]
// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppParallel)]]
// [[Rcpp::depends(StanHeaders)]]
#include <stan/math.hpp>                         // pulls in everything from rev/ and prim/
#include <Rcpp.h>
#include <RcppEigen.h>                           // do this AFTER including stan/math

// [[Rcpp::plugins(cpp14)]]

/* Ordinary Differential Equations */

// [[Rcpp::export]]
auto nonstiff(Eigen::MatrixXd A, Eigen::VectorXd y0) {
  using stan::math::integrate_ode_rk45;
  using stan::math::to_vector;
  using stan::math::to_array_1d;
  std::vector<double> theta;
  std::vector<double> times = {1, 2};
  auto y = integrate_ode_rk45([&A](auto t, auto y, 
                                   auto theta, auto x_r, auto x_i, std::ostream *msgs) {
    return to_array_1d( (A * to_vector(y)).eval() );
  }, to_array_1d(y0), 0, times, theta, {}, {});
  return to_vector(y[0]) - (stan::math::matrix_exp(A) * y0).eval(); // should be "zero"
}

The nonstiff function is then compiled into R’s global environment via Rcpp::sourceCpp (the same mechanism we use for rstan::expose_stan_functions) and called from R like

A <- matrix(c(0, -1, 1, -runif(1)), nrow = 2, ncol = 2)
y0 <- rexp(2)
all.equal(nonstiff(A, y0), c(0, 0), tol = 1e-5)

If it is compiled by clang++-9 with -fsanitize=undefined, then you get

SUMMARY: UndefinedBehaviorSanitizer: undefined-behavior /usr/local/lib/clang/9.0.1/include/emmintrin.h:1581:10 in
*** caught segfault ***
address (nil), cause ‘memory not mapped’
Traceback:
1: .Call(<pointer: 0x7fb289b17920>, A, y0)
2: nonstiff(A, y0)

and a traceback with a bunch of Eigen function calls

/usr/local/lib/clang/9.0.1/include/emmintrin.h:1581:10: runtime error: load of null pointer of type ‘__m128d’ (vector of 2 ‘double’ values)
#0 0x7fb289c5bef8 in double vector[2] Eigen::internal::pload<double vector[2]>(Eigen::internal::unpacket_traits<double vector[2]>::type const*) /data/gannet/ripley/R/test-clang/RcppEigen/include/Eigen/src/Core/arch/SSE/PacketMath.h:307:112
#1 0x7fb289c5bef8 in double vector[2] Eigen::internal::ploadt<double vector[2], 16>(Eigen::internal::unpacket_traits<double vector[2]>::type const*) /data/gannet/ripley/R/test-clang/RcppEigen/include/Eigen/src/Core/GenericPacketMath.h:463:12
#2 0x7fb289c5bef8 in double vector[2] Eigen::internal::evaluator<Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > >::packet<16, double vector[2]>(long) const /data/gannet/ripley/R/test-clang/RcppEigen/include/Eigen/src/Core/CoreEvaluators.h:204:12
#3 0x7fb289c5bef8 in double vector[2] Eigen::internal::binary_evaluator<Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const>, Eigen::internal::IndexBased, Eigen::internal::IndexBased, double, double>::packet<16, double vector[2]>(long) const /data/gannet/ripley/R/test-clang/RcppEigen/include/Eigen/src/Core/CoreEvaluators.h:734:50
#4 0x7fb289c5bc23 in void Eigen::internal::generic_dense_assignment_kernel<Eigen::internal::evaluator<Eigen::Matrix<double, -1, 1, 0, -1, 1> >, Eigen::internal::evaluator<Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> >, Eigen::internal::assign_op<double, double>, 0>::assignPacket<16, 16, double vector[2]>(long) /data/gannet/ripley/R/test-clang/RcppEigen/include/Eigen/src/Core/AssignEvaluator.h:658:87
#5 0x7fb289c5b91a in Eigen::internal::dense_assignment_loop<Eigen::internal::generic_dense_assignment_kernel<Eigen::internal::evaluator<Eigen::Matrix<double, -1, 1, 0, -1, 1> >, Eigen::internal::evaluator<Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> >, Eigen::internal::assign_op<double, double>, 0>, 3, 0>::run(Eigen::internal::generic_dense_assignment_kernel<Eigen::internal::evaluator<Eigen::Matrix<double, -1, 1, 0, -1, 1> >, Eigen::internal::evaluator<Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> >, Eigen::internal::assign_op<double, double>, 0>&) /data/gannet/ripley/R/test-clang/RcppEigen/include/Eigen/src/Core/AssignEvaluator.h:416:23
#6 0x7fb289c5b5db in void Eigen::internal::call_dense_assignment_loop<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const>, Eigen::internal::assign_op<double, double> >(Eigen::Matrix<double, -1, 1, 0, -1, 1>&, Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> const&, Eigen::internal::assign_op<double, double> const&) /data/gannet/ripley/R/test-clang/RcppEigen/include/Eigen/src/Core/AssignEvaluator.h:741:3
#7 0x7fb289c5b31b in Eigen::Matrix<double, -1, 1, 0, -1, 1>& Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >::_set_noalias<Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> >(Eigen::DenseBase<Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> > const&) /data/gannet/ripley/R/test-clang/RcppEigen/include/Eigen/src/Core/PlainObjectBase.h:732:7
#8 0x7fb289c5b146 in void Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >::_init1<Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const>, Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> >(Eigen::DenseBase<Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> > const&) /data/gannet/ripley/R/test-clang/RcppEigen/include/Eigen/src/Core/PlainObjectBase.h:816:13
#9 0x7fb289c5b146 in Eigen::Matrix<double, -1, 1, 0, -1, 1>::Matrix<Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> >(Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> const&) /data/gannet/ripley/R/test-clang/RcppEigen/include/Eigen/src/Core/Matrix.h:294:22
#10 0x7fb289c59f72 in Eigen::DenseBase<Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> >::eval() const /data/gannet/ripley/R/test-clang/RcppEigen/include/Eigen/src/Core/DenseBase.h:406:14
#11 0x7fb289c59f72 in SEXPREC* Rcpp::RcppEigen::eigen_wrap_is_plain<Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> >(Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> const&, Rcpp::traits::integral_constant<bool, false>) /data/gannet/ripley/R/test-clang/RcppEigen/include/RcppEigenWrap.h:138:45
#12 0x7fb289b17b7d in SEXPREC* Rcpp::RcppEigen::eigen_wrap<Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> >(Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> const&) /data/gannet/ripley/R/test-clang/RcppEigen/include/RcppEigenWrap.h:149:20
#13 0x7fb289b17b7d in SEXPREC* Rcpp::internal::wrap_dispatch_eigen<Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> >(Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> const&, Rcpp::traits::integral_constant<bool, true>) /data/gannet/ripley/R/test-clang/Rcpp/include/Rcpp/internal/wrap.h:776:20
#14 0x7fb289b17b7d in SEXPREC* Rcpp::internal::wrap_dispatch_unknown_importable<Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> >(Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> const&, Rcpp::traits::integral_constant<bool, false>) /data/gannet/ripley/R/test-clang/Rcpp/include/Rcpp/internal/wrap.h:787:20
#15 0x7fb289b17b7d in SEXPREC* Rcpp::internal::wrap_dispatch<Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> >(Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> const&, Rcpp::traits::wrap_type_unknown_tag) /data/gannet/ripley/R/test-clang/Rcpp/include/Rcpp/internal/wrap.h:807:20
#16 0x7fb289b17b7d in SEXPREC* Rcpp::wrap<Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> >(Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> const&) /data/gannet/ripley/R/test-clang/Rcpp/include/Rcpp/internal/wrap_end.h:30:15
#17 0x7fb289b17b7d in sourceCpp_3_nonstiff /tmp/RtmpzeUiMS/working_dir/RtmpNSAfLH/sourceCpp-x86_64-pc-linux-gnu-1.0.3/sourcecpp_bf3246adfc6b/filebf325c97c225.cpp:74:23

For the full log, see

https://www.stats.ox.ac.uk/pub/bdr/memtests/clang-UBSAN/StanHeaders/build_vignettes.log

Does anyone see what the problem is?

1 Like

There is a little bit more information from g++:

0x6030000515f0 is located 16 bytes inside of 32-byte region [0x6030000515e0,0x603000051600)
freed by thread T0 here:
#0 0x7f54499db85f in __interceptor_free (/lib64/libasan.so.5+0x10d85f)
#1 0x7f54327bb008 in void Eigen::internal::conditional_aligned_delete_auto<double, true>(double*, unsigned long) /data/gannet/ripley/R/test-4.0/RcppEigen/include/Eigen/src/Core/util/Memory.h:98
#2 0x7f54327bb008 in Eigen::DenseStorage<double, -1, -1, 1, 0>::~DenseStorage() /data/gannet/ripley/R/test-4.0/RcppEigen/include/Eigen/src/Core/DenseStorage.h:542
#3 0x7f54327bb008 in Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >::~PlainObjectBase() /data/gannet/ripley/R/test-4.0/RcppEigen/include/Eigen/src/Core/PlainObjectBase.h:98
#4 0x7f54327bb008 in Eigen::Matrix<double, -1, 1, 0, -1, 1>::~Matrix() /data/gannet/ripley/R/test-4.0/RcppEigen/include/Eigen/src/Core/Matrix.h:178
#5 0x7f54327bb008 in nonstiff(Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>) /tmp/RtmpaumA1z/working_dir/RtmpDuSnMH/sourceCpp-x86_64-pc-linux-gnu-1.0.3/sourcecpp_7c7b40a7fef1/file7c7b76deac8b.cpp:36
#6 0x7f54327c37f4 in sourceCpp_3_nonstiff /tmp/RtmpaumA1z/working_dir/RtmpDuSnMH/sourceCpp-x86_64-pc-linux-gnu-1.0.3/sourcecpp_7c7b40a7fef1/file7c7b76deac8b.cpp:74
#7 0x5787eb in R_doDotCall /data/gannet/ripley/R/svn/R-devel/src/main/dotcode.c:601
#8 0xfffb04d70c1 ()
previously allocated by thread T0 here:
#0 0x7f54499dbc58 in __interceptor_malloc (/lib64/libasan.so.5+0x10dc58)
#1 0x7f543285bca8 in Eigen::internal::handmade_aligned_malloc(unsigned long) /data/gannet/ripley/R/test-4.0/RcppEigen/include/Eigen/src/Core/util/Memory.h:88
#2 0x7f543285bca8 in Eigen::internal::aligned_malloc(unsigned long) /data/gannet/ripley/R/test-4.0/RcppEigen/include/Eigen/src/Core/util/Memory.h:164
#3 0x7f54327bd6df in void* Eigen::internal::conditional_aligned_malloc(unsigned long) /data/gannet/ripley/R/test-4.0/RcppEigen/include/Eigen/src/Core/util/Memory.h:214
#4 0x7f54327bd6df in double* Eigen::internal::conditional_aligned_new_auto<double, true>(unsigned long) /data/gannet/ripley/R/test-4.0/RcppEigen/include/Eigen/src/Core/util/Memory.h:374
#5 0x7f54327bd6df in Eigen::DenseStorage<double, -1, -1, 1, 0>::resize(long, long, long) /data/gannet/ripley/R/test-4.0/RcppEigen/include/Eigen/src/Core/DenseStorage.h:557
#6 0x7f54327bd6df in Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >::resize(long, long) /data/gannet/ripley/R/test-4.0/RcppEigen/include/Eigen/src/Core/PlainObjectBase.h:293
#7 0x7f54327bd6df in Eigen::internal::Assignment<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Product<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, 0>, Eigen::internal::assign_op<double, double>, Eigen::internal::Dense2Dense, void>::run(Eigen::Matrix<double, -1, 1, 0, -1, 1>&, Eigen::Product<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, 0> const&, Eigen::internal::assign_op<double, double> const&) /data/gannet/ripley/R/test-4.0/RcppEigen/include/Eigen/src/Core/ProductEvaluators.h:146
#8 0x7f54327bd6df in void Eigen::internal::call_assignment_no_alias<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Product<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, 0>, Eigen::internal::assign_op<double, double> >(Eigen::Matrix<double, -1, 1, 0, -1, 1>&, Eigen::Product<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, 0> const&, Eigen::internal::assign_op<double, double> const&) /data/gannet/ripley/R/test-4.0/RcppEigen/include/Eigen/src/Core/AssignEvaluator.h:836
#9 0x7f54327bd6df in Eigen::Matrix<double, -1, 1, 0, -1, 1>& Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >::_set_noalias<Eigen::Product<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, 0> >(Eigen::DenseBase<Eigen::Product<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, 0> > const&) /data/gannet/ripley/R/test-4.0/RcppEigen/include/Eigen/src/Core/PlainObjectBase.h:732
#10 0x7f54327bd6df in void Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >::_init1<Eigen::Product<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, 0>, Eigen::Product<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, 0> >(Eigen::DenseBase<Eigen::Product<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, 0> > const&) /data/gannet/ripley/R/test-4.0/RcppEigen/include/Eigen/src/Core/PlainObjectBase.h:816
#11 0x7f54327bd6df in Eigen::Matrix<double, -1, 1, 0, -1, 1>::Matrix<Eigen::Product<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, 0> >(Eigen::Product<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, 0> const&) /data/gannet/ripley/R/test-4.0/RcppEigen/include/Eigen/src/Core/Matrix.h:294
#12 0x7f54327bd6df in Eigen::DenseBase<Eigen::Product<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, 0> >::eval() const /data/gannet/ripley/R/test-4.0/RcppEigen/include/Eigen/src/Core/DenseBase.h:406
#13 0x7f54327bd6df in nonstiff(Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>) /tmp/RtmpaumA1z/working_dir/RtmpDuSnMH/sourceCpp-x86_64-pc-linux-gnu-1.0.3/sourcecpp_7c7b40a7fef1/file7c7b76deac8b.cpp:36
#14 0x7f54327c37f4 in sourceCpp_3_nonstiff /tmp/RtmpaumA1z/working_dir/RtmpDuSnMH/sourceCpp-x86_64-pc-linux-gnu-1.0.3/sourcecpp_7c7b40a7fef1/file7c7b76deac8b.cpp:74
#15 0x5787eb in R_doDotCall /data/gannet/ripley/R/svn/R-devel/src/main/dotcode.c:601
#16 0xfffb04d70c1 ()
SUMMARY: AddressSanitizer: heap-use-after-free /usr/lib/gcc/x86_64-redhat-linux/9/include/emmintrin.h:124 in _mm_load_pd
Shadow bytes around the buggy address:
0x0c0680002260: fa fa fd fd fd fd fa fa fd fd fd fd fa fa fd fd
0x0c0680002270: fd fd fa fa fd fd fd fd fa fa fd fd fd fd fa fa
0x0c0680002280: fd fd fd fd fa fa fd fd fd fd fa fa fd fd fd fd
0x0c0680002290: fa fa fd fd fd fd fa fa fd fd fd fd fa fa fd fd
0x0c06800022a0: fd fd fa fa fd fd fd fd fa fa fd fd fd fd fa fa
=>0x0c06800022b0: fd fd fd fa fa fa fd fd fd fa fa fa fd fd[fd]fd
0x0c06800022c0: fa fa fd fd fd fd fa fa 00 00 00 00 fa fa fa fa
0x0c06800022d0: fa fa fa fa fa fa fa fa fa fa fa fa fa fa fa fa
0x0c06800022e0: fa fa fa fa fa fa fa fa fa fa fa fa fa fa fa fa
0x0c06800022f0: fa fa fa fa fa fa fa fa fa fa fa fa fa fa fa fa
0x0c0680002300: fa fa fa fa fa fa fa fa fa fa fa fa fa fa fa fa
Shadow byte legend (one shadow byte represents 8 application bytes):
Addressable: 00
Partially addressable: 01 02 03 04 05 06 07
Heap left redzone: fa
Freed heap region: fd
Stack left redzone: f1
Stack mid redzone: f2
Stack right redzone: f3
Stack after return: f5
Stack use after scope: f8
Global redzone: f9
Global init order: f6
Poisoned by user: f7
Container overflow: fc
Array cookie: ac
Intra object redzone: bb
ASan internal: fe
Left alloca redzone: ca
Right alloca redzone: cb
Shadow gap: cc
==31867==ABORTING

If it’s something with intrinsics you can get segfaults when you pass unaligned memory addresses to special instructions that require aligned addresses.

I don’t know why this would have to do with our ODE stuff.

I don’t know either. I do know it didn’t happen with StanHeaders 2.19, but I don’t know offhand what changed with ODEs in 2.20 or 2.21.

Where can I install this development StanHeaders from?

It is on CRAN now.

Oh yeah this is one of the instructions that needs to be 16 byte aligned. Not that it’s guaranteed to be this, but it could be.

I just got 2.21 off Cran and things compiled and ran for me. Is there something special about this system?

Compiled with -fsanitize=undefined in CXX14FLAGS. See also

https://cran.r-project.org/doc/manuals/r-release/R-exts.html#Using-Undefined-Behaviour-Sanitizer

I’ve not been able to get this to trigger for me yet.

Did you need to build RcppEigen with -fsanitize=undefined? Or just the c++ you posted?

Since I haven’t been able to dig in, are you pretty convinced this is something inside the ODE solver? There isn’t an easy way to trigger it with regular Eigen stuff outside the ODE stuff?

Maybe Rcpp and RcppEigen need to be rebuilt with -fsanitize=undefined. The last time one of these things happened @yizhang pointed out that I didn’t use .eval() but I think I got it right here.

Have you been able to reproduce this not on the Cran computers?

(Gave up getting it to blow up with g++, trying to get clang++ to work but getting weird link errors)

It is necessary to compile Rcpp, RcppEigen, and the StanHeaders vignette with the same compiler and settings. The laptop I am on right now does not have clang++-9, but I will try to reproduce it later.

I wasn’t ever able to get the example to compile with Rcpp (linking errors with the sanitizer), but I did get it to compile as a test model (edit: in a clone of stan-math) with clang+±9.0.0 with the sanitizer and no errors. I didn’t manage to ever reproduce the error though, so I guess that doesn’t tell us much at all.

Would valgrind be of help here? Perhaps it could tell something if the program compiles as a standalone executable outside of R, unless the problem is in RcppEigen or somewhere in R 4.0.

What makes this all suspect to me is that the complaint is about Eigen while the non-stiff ODE solver does not use any Eigen stuff at all.

Could you rewrite the example without the use of auto and please also pass into integrate_ode_rk45 a proper instance of a functor instead of lambda magic. Inside the ODE RHS I would also start with declaring things one by one. So create an Eigen temporary for the vector, for the dy/dt and then return to_array_1d of that.

Maybe/hopefully that helps already.

I will try digging into this today. It may actually be with the matrix exponential that the numerical solution is being compared to.

I think I have fixed this.

What was it?