Eigen::Map'd matrix from raw buffer gives OBJECT_ALLOCATED_ON_STACK_IS_TOO_BIG

Question

Recently I have been working with Eigen matrices derived from raw buffers and I noticed this curious case:

#include <eigen3/Eigen/Dense>

int main(int argc, char const *argv[]) {
    /* code */
    const int M = 320;
    const int N = 640;
    const int K = 320;
    const int alpha = 2;
    const int beta = 1;
    Eigen::Matrix<int32_t, Eigen::Dynamic,Eigen::Dynamic> A = Eigen::Matrix<int32_t, Eigen::Dynamic,Eigen::Dynamic>::Random(M,K);
    Eigen::Matrix<int32_t, Eigen::Dynamic,Eigen::Dynamic> B = Eigen::Matrix<int32_t, Eigen::Dynamic,Eigen::Dynamic>::Random(K,N);
    Eigen::Matrix<int32_t, Eigen::Dynamic,Eigen::Dynamic> C = Eigen::Matrix<int32_t, Eigen::Dynamic,Eigen::Dynamic>::Random(M,N);

    //Following http://eigen.tuxfamily.org/dox/TopicWritingEfficientProductExpression.html

    C.noalias() += (A*alpha)*(B*beta); //WORKS

    C.noalias() += A*B;

    Eigen::Map<Eigen::Matrix<int32_t, M, K, Eigen::ColMajor> > map_a(A.data());
    Eigen::Map<Eigen::Matrix<int32_t, K, N, Eigen::ColMajor> > map_b(B.data());
    Eigen::Map<Eigen::Matrix<int32_t, M, N, Eigen::ColMajor> > map_c(C.data());

    map_c.noalias() += map_a*map_b; //WORKS

    map_c.noalias() += (map_a*alpha)*(map_b*beta); //COMPILE ERROR HERE

    return 0;
}

If I have big matrix dimensions, I can't allocate on the stack or I would get a OBJECT_ALLOCATED_ON_STACK_IS_TOO_BIG, therefore I use the Eigen dynamic allocator.

However, it seems that if I have a raw buffer and I map it to a matrix, i can not perform a BLAS 3 like gemm multiplication (C+= (alpha*A)*(beta*B)), due to compilation error: OBJECT_ALLOCATED_ON_STACK_IS_TOO_BIG. If I do a simple C += A*B everything works as expected.

In the example case, I map the raw buffer from a matrix allocated by Eigen, but in principle it could be the raw buffer from anything (such as std::vector).

Any idea what is happening here? As far as I can tell everything here should be heap allocated, and even if it weren't, why would C += A*B work with the mapped memory matrices and C+= (alpha*A)*(beta*B) would not?

Cheers,

Nick

Answer 1

For such large matrices better use runtime sizes as in Avi Ginsburg's answer. That being said, I'll now explain what's going on within Eigen. The problem is that within the matrix product implementation, we have a branch like that (simplified):

if(<too small>)
  lazyproduct::eval(dst, lhs, rhs);
else
  gemm::eval(dst,lhs, rhs);

If the product is too small, instead of calling the heavy "gemm" routine, we fall back to a coefficient-based implementation, in your case:

map_c.noalias() += (map_a*alpha).lazyProduct(map_b*beta);

This path does not rewrite the expression as (alpha*beta)*(map_a*map_b) and therefore, in order to avoid recomputing map_a*alpha and map_b*beta many times, the strategy is to backed them within temporaries... hence the compilation error.

Of course, in your case this path will never be taken, and it would even be completely removed by the compiler if you increase EIGEN_STACK_ALLOCATION_LIMIT because the condition if(<too small>) is known at compile-time. How sad.