Работа с непрерывным вектором матриц фиксированного размера для обеих схем хранения в Eigen

Question

Внешняя библиотека дает мне необработанный указатель double s, который я хочу отобразить на собственный тип. Необработанный массив логически представляет собой большой упорядоченный набор небольших плотных матриц фиксированного размера, одинакового размера. Основная проблема заключается в том, что маленькие плотные матрицы могут располагаться в мажорном порядке по строкам или по столбцам, и я хочу разместить их обоих.

Мой текущий подход заключается в следующем. Обратите внимание, что все записи небольшого блока фиксированного размера (в массиве блоков) должны быть смежными в памяти.

template<int bs, class Mattype>
void block_operation(double *const vals, const int numblocks)
{
    Eigen::Map<Mattype> mappedvals(vals, 
                Mattype::IsRowMajor ? numblocks*bs : bs,
                Mattype::IsRowMajor ? bs : numblocks*bs
            );

    for(int i = 0; i < numblocks; i++)
      if(Mattype::isRowMajor)
          mappedvals.template block<bs,bs>(i*bs,0) = block_operation_rowmajor(mappedvals);
      else
          mappedvals.template block<bs,bs>(0,i*bs) = block_operation_colmajor(mappedvals);
}

Вызывающая функция сначала определяет Mattype (из двух вариантов), а затем вызывает вышеуказанную функцию с правильным параметром шаблона.

Таким образом, все мои алгоритмы должны быть написаны дважды, и мой код чередуется с этими проверками компоновки. Есть ли способ сделать это независимым от макета способом? Помните, что этот код должен быть максимально быстрым.

В идеале я хотел бы Map данных только один раз и использовать их для всех необходимых операций. Однако единственное решение, которое я мог придумать, - это вызывать конструктор Map один раз для каждого небольшого блока, когда мне нужно получить доступ к блоку.

template<int bs, StorageOptions layout>
inline Map<Matrix<double,bs,bs,layout>> extractBlock(double *const vals, 
                                                     const int bindex)
{
  return Map<Matrix<double,bs,bs,layout>>(vals+bindex*bs*bs);
}

Будет ли эта функция полностью оптимизирована (с помощью современного компилятора, такого как GCC 7.3 или Intel 2017 под -std=c++14 -O3), или я буду платить небольшой штраф каждый раз, когда я вызываю эту функцию (один раз для каждого блока, и там много маленьких блоков)? Есть ли лучший способ сделать это?

Answer 1

Ваш extractBlock в порядке, более простое, но несколько уродливое решение - использовать переосмысление приведения в начале block_operation:

using BlockType = Matrix<double,bs,bs,layout|DontAlign>;
BlockType* blocks = reinterpret_cast<BlockType*>(vals);

for(int i...)
  block[i] = ...;

Это будет работать только для матриц фиксированных размеров. Также обратите внимание на DontAlign, который важен, если только вы не можете гарантировать, что vals выровнен на 16 или даже 32 байта в зависимости от наличия AVX и bs .... так что просто используйте DontAlign!