Матрица Eigen :: Map'd из необработанного буфера дает OBJECT_ALLOCATED_ON_STACK_IS_TOO_BIG

Question

Недавно я работал с матрицами Eigen, полученными из необработанных буферов, и заметил этот любопытный случай:

#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;
}

Если у меня большие размеры матрицы, я не могу выделить в стеке или получуOBJECT_ALLOCATED_ON_STACK_IS_TOO_BIG, поэтому я использую собственный динамический распределитель.

Однако, похоже, что если у меня есть необработанный буфер, и я сопоставляю его с матрицей, я не могу выполнить BLAS 3-подобное умножение gemm (C+= (alpha*A)*(beta*B)) из-за ошибки компиляции: OBJECT_ALLOCATED_ON_STACK_IS_TOO_BIG.Если я сделаю простой C += A*B, все будет работать, как и ожидалось.

В этом примере я отображаю необработанный буфер из матрицы, выделенной Eigen, но в принципе это может быть необработанный буфер из чего угодно (например,std::vector).

Есть идеи, что здесь происходит?Насколько я могу сказать, здесь все должно быть выделено кучей, и даже если бы этого не было, почему C += A*B будет работать с отображенными матрицами памяти, а C+= (alpha*A)*(beta*B) не будет?

Приветствия,

Ник

Answer 1

Для таких больших матриц лучше использовать размеры времени выполнения, как в ответе Ави Гинзбурга.При этом я сейчас объясню, что происходит внутри Эйгена.Проблема в том, что в реализации матричного продукта у нас есть такая ветвь (упрощенная):

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

Если продукт слишком маленький, вместо того, чтобы вызывать тяжелую подпрограмму «гем», мы возвращаемся креализация на основе коэффициентов, в вашем случае:

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

Этот путь не переписывает выражение как (alpha*beta)*(map_a*map_b) и, следовательно, во избежание повторного вычисления map_a*alpha и map_b*beta много раз, стратегия заключается в том, чтобы поддерживать их во временных файлах ... отсюда и ошибка компиляции.

Конечно, в вашем случае этот путь никогда не будет выбран, и он даже будет полностью удален компилятором, есливы увеличиваете EIGEN_STACK_ALLOCATION_LIMIT, потому что условие if(<too small>) известно во время компиляции.Как грустно.

Answer 2

Ваши Map s оборачивают матрицы статического размера, например:

Eigen::Map<Eigen::Matrix<int32_t, M, K, Eigen::ColMajor> > 
                                  ^  ^

Используйте взамен Map s динамического размера:

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

Это не означает, что выможет изменить размер Map s, просто указывает, как выделяются временные значения.