Вызов Eigen GEMM поверх внешних данных

Question

Eigen имеет удивительно быструю реализацию GEMM, поэтому я хочу использовать ее в своей библиотеке тензоров проекта для домашних животных. Если я правильно понимаю, это возможно через Eigen :: Map. Я написал простой пример и определил EIGEN_NO_MALLOC, чтобы убедиться в отсутствии нежелательных выделений.

Хорошо работает с простым умножением матриц, например C += A * B. Но, к сожалению, он не может справиться с ситуацией C += alpha * A * B (как в GEMM).

#include <iostream>
#include <vector>

#define EIGEN_NO_MALLOC
#include "Eigen/Core"

int main()
{
    using Scalar = float;
    using namespace Eigen;
    std::vector<Scalar> aDat = {1, 2, 3, 4};
    std::vector<Scalar> bDat = {1, 2, 3, 4};
    std::vector<Scalar> cDat = {1, 2, 3, 4};
    Map<Matrix<Scalar, -1, -1, RowMajor>, Unaligned> a(aDat.data(), 2, 2);
    Map<Matrix<Scalar, -1, -1, RowMajor>, Unaligned> b(bDat.data(), 2, 2);
    Map<Matrix<Scalar, -1, -1, RowMajor>, Unaligned> c(cDat.data(), 2, 2);

    //Ok
    c.noalias() += a * b;

    //Assertion `false && "heap allocation is forbidden.....
    c.noalias() += 2 * a * b;

    return 0;
}

c.noalias() += 2 * a * b; дает мне следующую ошибку времени выполнения

a.out: path_to_eigen/Eigen/src/Core/util/Memory.h:129: void Eigen::internal::check_that_malloc_is_allowed(): Assertion `false && "heap allocation is forbidden (EIGEN_NO_MALLOC is defined)"' failed.

Можно ли позвонить c.noalias() += someScalar * a * b; без выделений?

PS Моя собственная версия 3.3.2, gcc версия 7.2.0

Извините за мой плохой английский

Answer 1

Это похоже на ошибку, вносимую оптимизацией небольшого размера. Для продуктов (KxM)*(MxN), где K+M+N < EIGEN_GEMM_TO_COEFFBASED_THRESHOLD (по умолчанию это 20), Eigen переключается на реализацию отложенного продукта, которая, очевидно, считает, что оценка во временную является хорошей идеей. Ваш пример действительно будет работать нормально, например для 5x5 * 5x10 продукта.

Я подал ошибку для этой проблемы: http://eigen.tuxfamily.org/bz/show_bug.cgi?id=1562

Если вы хотите быстро обойти, определите -D EIGEN_GEMM_TO_COEFFBASED_THRESHOLD=1, при этом всегда будет использоваться реализация big matrix -GEMM.

С другой стороны, если вы знаете свой размер матрицы во время компиляции, вам лучше использовать соответствующий тип фиксированного размера (Matrix<Scalar, 2, 2, RowMajor> в вашем случае).