Indexing arrays of matrices within operations results in Eigen compilation errors

[andrjohns]

Debugging a build failure for the ctsem R package against the development rstan, and narrowed it down to some Eigen-expression-issues with Stan's indexing.

Given model:

data {
  real y_mean;
}

parameters {
  real y;
}

model {
  matrix[2, 2] m1 = [[0,0],[0,0]];
  array[1,1] matrix[2, 2] m2 = {{ m1 }};
  matrix[2,2] res;

  // No compilation error
  res = m1 * m1';

  // Compilation error
  res = m2[1,1] * m2[1,1]';

  y ~ normal(y_mean, 1);
}

Compilation fails with:

--- Compiling C++ code ---
g++ -std=c++17 -pthread -D_REENTRANT -Wno-sign-compare -Wno-ignored-attributes -Wno-class-memaccess      -isystem stan/lib/stan_math/lib/tbb_2020.3/include    -O3 -I src -I stan/src -I stan/lib/rapidjson_1.1.0/ -I lib/CLI11-1.9.1/ -I stan/lib/stan_math/ -isystem stan/lib/stan_math/lib/eigen_5.0.1 -isystem stan/lib/stan_math/lib/boost_1.87.0 -isystem stan/lib/stan_math/lib/sundials_6.1.1/include -isystem stan/lib/stan_math/lib/sundials_6.1.1/src/sundials    -DBOOST_DISABLE_ASSERTS          -c -Wno-ignored-attributes   -x c++ -o examples/bernoulli/model.o examples/bernoulli/model.hpp
In file included from stan/lib/stan_math/lib/eigen_5.0.1/Eigen/Core:418,
                 from stan/lib/stan_math/lib/eigen_5.0.1/Eigen/Dense:1,
                 from stan/lib/stan_math/stan/math/prim/fun/Eigen.hpp:32,
                 from stan/lib/stan_math/stan/math/rev/core/Eigen_NumTraits.hpp:4,
                 from stan/lib/stan_math/stan/math/rev/core.hpp:4,
                 from stan/src/stan/model/model_base.hpp:8,
                 from stan/src/stan/model/model_header.hpp:4:
stan/lib/stan_math/lib/eigen_5.0.1/Eigen/src/Core/products/GeneralMatrixMatrix.h: In instantiation of ‘void Eigen::internal::gemm_functor<Scalar, Index, Gemm, Lhs, Rhs, Dest, BlockingType>::operator()(Index, Index, Index, Index, Eigen::internal::GemmParallelInfo<Index>*) const [with Scalar = double; Index = long int; Gemm = Eigen::internal::general_matrix_matrix_product<long int, double, 0, false, double, 1, false, 0, 1>; Lhs = Eigen::Matrix<double, -1, -1>; Rhs = stan::math::Holder<Eigen::Transpose<Eigen::Matrix<double, -1, -1> >, Eigen::Matrix<double, -1, -1, 0, -1, -1> >; Dest = Eigen::Matrix<double, -1, -1>; BlockingType = Eigen::internal::gemm_blocking_space<0, double, double, -1, -1, -1, 1, false>]’:
stan/lib/stan_math/lib/eigen_5.0.1/Eigen/src/Core/products/Parallelizer.h:110:7:   required from ‘void Eigen::internal::parallelize_gemm(const Functor&, Index, Index, Index, bool) [with bool Condition = true; Functor = gemm_functor<double, long int, general_matrix_matrix_product<long int, double, 0, false, double, 1, false, 0, 1>, Eigen::Matrix<double, -1, -1>, stan::math::Holder<Eigen::Transpose<Eigen::Matrix<double, -1, -1> >, Eigen::Matrix<double, -1, -1, 0, -1, -1> >, Eigen::Matrix<double, -1, -1>, gemm_blocking_space<0, double, double, -1, -1, -1, 1, false> >; Index = long int]’
  110 |   func(0, rows, 0, cols);
      |   ~~~~^~~~~~~~~~~~~~~~~~
stan/lib/stan_math/lib/eigen_5.0.1/Eigen/src/Core/products/GeneralMatrixMatrix.h:449:107:   required from ‘static void Eigen::internal::generic_product_impl<Lhs, Rhs, Eigen::DenseShape, Eigen::DenseShape, 8>::scaleAndAddTo(Dest&, const Lhs&, const Rhs&, const Scalar&) [with Dest = Eigen::Matrix<double, -1, -1>; Lhs = Eigen::Matrix<double, -1, -1>; Rhs = stan::math::Holder<Eigen::Transpose<Eigen::Matrix<double, -1, -1> >, Eigen::Matrix<double, -1, -1, 0, -1, -1> >; Scalar = double]’
  449 |     internal::parallelize_gemm<(Dest::MaxRowsAtCompileTime > 32 || Dest::MaxRowsAtCompileTime == Dynamic)>(
      |     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^
  450 |         GemmFunctor(lhs, rhs, dst, actualAlpha, blocking), a_lhs.rows(), a_rhs.cols(), a_lhs.cols(),
      |         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~       
  451 |         Dest::Flags & RowMajorBit);
      |         ~~~~~~~~~~~~~~~~~~~~~~~~~~                                                                         
stan/lib/stan_math/lib/eigen_5.0.1/Eigen/src/Core/products/GeneralMatrixMatrix.h:392:20:   required from ‘static void Eigen::internal::generic_product_impl<Lhs, Rhs, Eigen::DenseShape, Eigen::DenseShape, 8>::evalTo(Dst&, const Lhs&, const Rhs&) [with Dst = Eigen::Matrix<double, -1, -1>; Lhs = Eigen::Matrix<double, -1, -1>; Rhs = stan::math::Holder<Eigen::Transpose<Eigen::Matrix<double, -1, -1> >, Eigen::Matrix<double, -1, -1, 0, -1, -1> >]’
  392 |       scaleAndAddTo(dst, lhs, rhs, Scalar(1));
      |       ~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~
stan/lib/stan_math/lib/eigen_5.0.1/Eigen/src/Core/ProductEvaluators.h:112:75:   required from ‘Eigen::internal::product_evaluator<Eigen::Product<Lhs, Rhs, Option>, ProductTag, LhsShape, RhsShape>::product_evaluator(const XprType&) [with Lhs = Eigen::Matrix<double, -1, -1>; Rhs = stan::math::Holder<Eigen::Transpose<Eigen::Matrix<double, -1, -1> >, Eigen::Matrix<double, -1, -1, 0, -1, -1> >; int Options = 0; int ProductTag = 8; LhsShape = Eigen::DenseShape; RhsShape = Eigen::DenseShape; typename Eigen::internal::traits<typename Eigen::Product<Lhs, Rhs, Option>::Lhs>::Scalar = double; typename Eigen::Product<Lhs, Rhs, Option>::Lhs = Eigen::Matrix<double, -1, -1>; typename Eigen::internal::traits<typename Eigen::Product<Lhs, Rhs, Option>::Rhs>::Scalar = double; typename Eigen::Product<Lhs, Rhs, Option>::Rhs = stan::math::Holder<Eigen::Transpose<Eigen::Matrix<double, -1, -1> >, Eigen::Matrix<double, -1, -1, 0, -1, -1> >; XprType = Eigen::Product<Eigen::Matrix<double, -1, -1>, stan::math::Holder<Eigen::Transpose<Eigen::Matrix<double, -1, -1> >, Eigen::Matrix<double, -1, -1, 0, -1, -1> >, 0>]’
  112 |     generic_product_impl<Lhs, Rhs, LhsShape, RhsShape, ProductTag>::evalTo(m_result, xpr.lhs(), xpr.rhs());
      |     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
stan/lib/stan_math/lib/eigen_5.0.1/Eigen/src/Core/ProductEvaluators.h:35:90:   required from ‘Eigen::internal::evaluator<Eigen::Product<Lhs, Rhs, Option> >::evaluator(const XprType&) [with Lhs = Eigen::Matrix<double, -1, -1>; Rhs = stan::math::Holder<Eigen::Transpose<Eigen::Matrix<double, -1, -1> >, Eigen::Matrix<double, -1, -1, 0, -1, -1> >; int Options = 0; XprType = Eigen::Product<Eigen::Matrix<double, -1, -1>, stan::math::Holder<Eigen::Transpose<Eigen::Matrix<double, -1, -1> >, Eigen::Matrix<double, -1, -1, 0, -1, -1> >, 0>]’
   35 |   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE explicit evaluator(const XprType& xpr) : Base(xpr) {}
      |                                                                                          ^
stan/lib/stan_math/stan/math/prim/meta/holder.hpp:277:44:   [ skipping 8 instantiation contexts, use -ftemplate-backtrace-limit=0 to disable ]
stan/src/stan/model/indexing/access_helpers.hpp:92:5:   required from ‘void stan::model::internal::assign_impl(T1&&, T2&&, const char*) [with T1 = Eigen::Matrix<double, -1, -1>&; T2 = stan::math::Holder<Eigen::Product<Eigen::Matrix<double, -1, -1>, stan::math::Holder<Eigen::Transpose<Eigen::Matrix<double, -1, -1> >, Eigen::Matrix<double, -1, -1, 0, -1, -1> >, 0>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, stan::math::Holder<Eigen::Transpose<Eigen::Matrix<double, -1, -1, 0, -1, -1> >, Eigen::Matrix<double, -1, -1, 0, -1, -1> > >; stan::require_all_eigen_t<T_desired, T_actual>* <anonymous> = 0]’
   92 |   x = std::forward<T2>(y);
      |   ~~^~~~~~~~~~~~~~~~~~~~~
stan/src/stan/model/indexing/assign.hpp:65:24:   required from ‘void stan::model::assign(T&&, U&&, const char*) [with T = Eigen::Matrix<double, -1, -1>&; U = stan::math::Holder<Eigen::Product<Eigen::Matrix<double, -1, -1>, stan::math::Holder<Eigen::Transpose<Eigen::Matrix<double, -1, -1> >, Eigen::Matrix<double, -1, -1, 0, -1, -1> >, 0>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, stan::math::Holder<Eigen::Transpose<Eigen::Matrix<double, -1, -1, 0, -1, -1> >, Eigen::Matrix<double, -1, -1, 0, -1, -1> > >; stan::require_t<std::is_assignable<typename std::decay<_Tp>::type&, typename std::decay<_Tp2>::type> >* <anonymous> = 0; stan::require_all_not_std_vector_t<T1, T2>* <anonymous> = 0; stan::math::require_all_not_tuple_t<T, U>* <anonymous> = 0]’
   65 |   internal::assign_impl(x, std::forward<U>(y), name);
      |   ~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~
examples/bernoulli/model.hpp:113:28:   required from ‘stan::scalar_type_t<T2> model_model_namespace::model_model::log_prob_impl(VecR&, VecI&, std::ostream*) const [with bool propto__ = false; bool jacobian__ = false; VecR = Eigen::Matrix<double, -1, 1>; VecI = Eigen::Matrix<int, -1, 1>; stan::require_vector_like_t<VecR>* <anonymous> = 0; stan::require_vector_like_vt<std::is_integral, VecI>* <anonymous> = 0; stan::require_not_st_var<VecR>* <anonymous> = 0; stan::scalar_type_t<T2> = double; std::ostream = std::basic_ostream<char>]’
  113 |         stan::model::assign(res,
      |         ~~~~~~~~~~~~~~~~~~~^~~~~
  114 |           stan::math::multiply(
      |           ~~~~~~~~~~~~~~~~~~~~~
  115 |             stan::model::rvalue(m2, "m2", stan::model::index_uni(1),
      |             ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  116 |               stan::model::index_uni(1)),
      |               ~~~~~~~~~~~~~~~~~~~~~~~~~~~
  117 |             stan::math::transpose(
      |             ~~~~~~~~~~~~~~~~~~~~~~
  118 |               stan::model::rvalue(m2, "m2", stan::model::index_uni(1),
      |               ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  119 |                 stan::model::index_uni(1)))), "assigning variable res");
      |                 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
examples/bernoulli/model.hpp:375:47:   required from ‘T_ model_model_namespace::model_model::log_prob(Eigen::Matrix<T, -1, 1>&, std::ostream*) const [with bool propto__ = false; bool jacobian__ = false; T_ = double; std::ostream = std::basic_ostream<char>]’
  375 |     return log_prob_impl<propto__, jacobian__>(params_r, params_i, pstream);
      |            ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~
stan/src/stan/model/model_base_crtp.hpp:93:80:   required from ‘double stan::model::model_base_crtp<M>::log_prob(Eigen::VectorXd&, std::ostream*) const [with M = model_model_namespace::model_model; Eigen::VectorXd = Eigen::Matrix<double, -1, 1>; std::ostream = std::basic_ostream<char>]’
   93 |     return static_cast<const M*>(this)->template log_prob<false, false, double>(
      |            ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^
   94 |         theta, msgs);
      |         ~~~~~~~~~~~~                                                            
stan/src/stan/model/model_base_crtp.hpp:91:17:   required from here
   91 |   inline double log_prob(Eigen::VectorXd& theta,
      |                 ^~~~~~~~
stan/lib/stan_math/lib/eigen_5.0.1/Eigen/src/Core/products/GeneralMatrixMatrix.h:206:102: error: passing ‘const stan::math::Holder<Eigen::Transpose<Eigen::Matrix<double, -1, -1> >, Eigen::Matrix<double, -1, -1, 0, -1, -1> >’ as ‘this’ argument discards qualifiers [-fpermissive]
  206 |     Gemm::run(rows, cols, m_lhs.cols(), &m_lhs.coeffRef(row, 0), m_lhs.outerStride(), &m_rhs.coeffRef(0, col),
      |                                                                                        ~~~~~~~~~~~~~~^~~~~~~~
In file included from stan/lib/stan_math/lib/eigen_5.0.1/Eigen/Core:350:
stan/lib/stan_math/lib/eigen_5.0.1/Eigen/src/Core/DenseCoeffsBase.h:316:59: note:   in call to ‘constexpr Eigen::DenseCoeffsBase<Derived, 1>::Scalar& Eigen::DenseCoeffsBase<Derived, 1>::coeffRef(Eigen::Index, Eigen::Index) [with Derived = stan::math::Holder<Eigen::Transpose<Eigen::Matrix<double, -1, -1> >, Eigen::Matrix<double, -1, -1, 0, -1, -1> >; Scalar = double; Eigen::Index = long int]’
  316 |   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Scalar& coeffRef(Index row, Index col) {
      |                                                           ^~~~~~~~
make: *** [make/program:77: examples/bernoulli/model.o] Error 1