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?