Что происходит внутри Eigen на простом примере
Рассмотрим следующую программу-пример:
#include<Eigen/Core> int main() { int size = 50; // VectorXf is a vector of floats, with dynamic size. Eigen::VectorXf u(size), v(size), w(size); u = v + w; }
Цель этой страницы — понять, как Eigen компилирует её, предполагая, что включена векторная оптимизация SSE2 (опция компилятора GCC -msse2).
Почему это интересно
Возможно, вы думаете, что приведенная выше программа-пример настолько проста, что её компиляция не должна вызывать никаких интересных моментов. Поэтому перед началом объясним, что нетривиально в её корректной компиляции — то есть в генерации оптимизированного кода — так что сложность Eigen, которую мы здесь объясним, действительно полезна.
Посмотрите на строку кода
u = v + w; // (*)
Первое важное замечание по поводу компиляции состоит в том, что массивы должны быть пройдены только один раз, как в
for(int i = 0; i < size; i++) u[i] = v[i] + w[i];
Проблема в том, что если мы создадим простую библиотеку C++, где класс VectorXf имеет оператор +, возвращающий VectorXf, то строка кода (*) будет эквивалентна:
VectorXf tmp = v + w; VectorXf u = tmp;
Очевидно, введение временной переменной tmp здесь бесполезно. Это очень плохо сказывается на производительности, во-первых, потому что создание tmp требует динамического выделения памяти в данном контексте, а во-вторых, теперь есть два цикла for:
for(int i = 0; i < size; i++) tmp[i] = v[i] + w[i]; for(int i = 0; i < size; i++) u[i] = tmp[i];
Обход массивов дважды вместо одного раза ужасно сказывается на производительности, поскольку это означает, что мы выполняем много избыточных обращений к памяти.
Второе важное замечание по поводу компиляции вышеприведенной программы состоит в правильном использовании инструкций SSE2. Обратите внимание, что Eigen также поддерживает AltiVec, и все рассуждения, которые мы здесь приводим, также относятся к AltiVec.
SSE2, как и AltiVec, представляет собой набор инструкций, позволяющих выполнять вычисления над пакетами из 128 бит одновременно. Поскольку float занимает 32 бита, это означает, что инструкции SSE2 могут обрабатывать 4 float одновременно. Это означает, что при правильном использовании они могут ускорить наши вычисления в 4 раза.
Однако в вышеприведенной программе мы выбрали размер=50, поэтому наши векторы состоят из 50 float, а 50 не является кратным 4. Это означает, что мы не можем надеяться выполнить все эти вычисления с помощью инструкций SSE2. Лучший вариант — обработать первые 48 коэффициентов с помощью инструкций SSE2, так как 48 — это наибольшее кратное 4, меньшее 50, а затем отдельно, без SSE2, обработать 49-й и 50-й коэффициенты. Что-то вроде этого:
for(int i = 0; i < 4*(size/4); i+=4) u.packet(i) = v.packet(i) + w.packet(i); for(int i = 4*(size/4); i < size; i++) u[i] = v[i] + w[i];
Итак, давайте посмотрим на нашу программу-пример строка за строкой и проследим за тем, как Eigen её компилирует.
Создание векторов
Давайте проанализируем первую строку:
Eigen::VectorXf u(size), v(size), w(size);
Прежде всего, VectorXf — это следующий typedef:
typedef Matrix<float, Dynamic, 1> VectorXf;
Шаблон класса Matrix объявлен в src/Core/util/ForwardDeclarations.h с 6 параметрами шаблона, но последние 3 автоматически определяются первыми 3. Поэтому сейчас вы можете ими пренебречь. Здесь Matrix<float, Dynamic, 1> означает матрицу float с динамическим числом строк и 1 столбцом.
Класс Matrix наследует базовый класс MatrixBase. Не беспокойтесь об этом пока, достаточно сказать, что MatrixBase объединяет матрицы/векторы и все типы выражений — подробнее об этом ниже.
Когда мы делаем
Eigen::VectorXf u(size);
Вызывается конструктор Matrix::Matrix(int) в src/Core/Matrix.h. Помимо некоторых проверок, он только создаёт член m_storage, который имеет тип DenseStorage<float, Dynamic, Dynamic, 1>.
Вы можете спросить, неужели излишне иметь хранилище в отдельном классе? Причина в том, что шаблон класса Matrix охватывает все типы матриц и векторов: как фиксированного, так и динамического размера. Способ хранения не одинаков в этих двух случаях. Для матриц фиксированного размера коэффициенты матрицы хранятся как обычный массив-член. Для матриц динамического размера коэффициенты хранятся как указатель на динамически выделенный массив. Именно поэтому нам нужно абстрагировать хранилище от класса Matrix. Это и есть DenseStorage.
Давайте посмотрим на этот конструктор в src/Core/DenseStorage.h. Вы можете увидеть, что здесь есть много частных шаблонов специализации DenseStorages, отдельно рассматривающих случаи, когда размерности Dynamic или фиксированы во время компиляции. Частная специализация, которую мы ищем, это:
template<typename T, int _Cols> class DenseStorage<T, Dynamic, Dynamic, _Cols>
Здесь вызывается конструктор DenseStorage::DenseStorage(int size, int rows, int columns) с size=50, rows=50, columns=1.
Вот этот конструктор:
inline DenseStorage(int size, int rows, int) : m_data(internal::aligned_new<T>(size)), m_rows(rows) {}
Здесь член m_data — это фактический массив коэффициентов матрицы. Как вы видите, он динамически выделяется. Вместо вызова new[] или malloc() у нас есть своя функция internal::aligned_new, определённая в src/Core/util/Memory.h. Если включена векторная оптимизация, то она использует платформоспецифическую функцию для выделения 128-битно-выровненного массива, так как это очень полезно для векторной оптимизации как SSE2, так и AltiVec. Если векторная оптимизация отключена, она эквивалентна стандартному new[].
Как вы можете видеть, конструктор также устанавливает член m_rows в значение size. Обратите внимание, что члена m_columns нет: действительно, в этой частной специализации DenseStorage мы знаем количество столбцов во время компиляции, так как параметр шаблона _Cols отличается от Dynamic. То есть в нашем случае _Cols равно 1, что означает, что наш вектор — это просто матрица с 1 столбцом. Поэтому нет необходимости хранить количество столбцов как переменную во время выполнения.
При вызове VectorXf::data() для получения указателя на массив коэффициентов, возвращается DenseStorage::data(), которая возвращает член m_data.
При вызове VectorXf::size() для получения размера вектора, это фактически метод базового класса MatrixBase. Он определяет, что вектор является вектор-столбцом, поскольку ColsAtCompileTime==1 (это происходит из параметров шаблона в typedef VectorXf). Он вычисляет, что размер — это количество строк, поэтому возвращает VectorXf::rows(), которая возвращает DenseStorage::rows(), которая возвращает член m_rows, который был установлен в значение size конструктором.
Создание выражения суммы
Теперь, когда наши векторы созданы, перейдём к следующей строке:
u = v + w;
В итоге, оператор + возвращает выражение «сумма векторов», но фактически вычисления не выполняет. Вычисления выполняет оператор =, который вызывается позже.
Теперь давайте посмотрим, что делает Eigen, когда видит это:
v + w
Здесь v и w имеют тип VectorXf, который является typedef для специализации Matrix (как мы объяснили выше), который является подклассом MatrixBase. Таким образом, вызывается
MatrixBase::operator+(const MatrixBase&)
Тип возвращаемого значения этого оператора
CwiseBinaryOp<internal::scalar_sum_op<float>, VectorXf, VectorXf>
Класс CwiseBinaryOp — это наша первая встреча с шаблоном выражений. Как мы уже говорили, оператор + сам по себе не выполняет никаких вычислений, он просто возвращает абстрактное выражение «сумма векторов». Поскольку существуют также выражения «разность векторов» и «покомпонентное произведение векторов», мы объединяем их все как «покомпонентные бинарные операции», которые сокращаем до «CwiseBinaryOp». «Покомпонентное» означает, что операция выполняется по коэффициентам. «Бинарное» означает, что есть два операнда — мы складываем два вектора друг с другом.
Теперь вы можете спросить, что произойдёт, если мы сделаем что-то вроде
v + w + u;
Первый v + w вернёт CwiseBinaryOp, как и выше, поэтому для компиляции этого нам нужно определить оператор + также в классе CwiseBinaryOp… На этом этапе это начинает выглядеть как кошмар: нам придётся определять все операторы в каждом из классов выражений (как вы догадались, CwiseBinaryOp — только один из многих)? Это тупик!
Решение в том, что CwiseBinaryOp, а также Matrix и все другие типы выражений, являются подклассами MatrixBase. Поэтому достаточно один раз определить операторы в классе MatrixBase.
Поскольку MatrixBase — это общий базовый класс различных подклассов, аспекты, зависящие от подкласса, должны быть абстрагированы от MatrixBase. Это называется полиморфизмом.
Классический подход к полиморфизму в C++ — через виртуальные функции. Это динамический полиморфизм. Здесь нам не нужен динамический полиморфизм, потому что вся структура библиотеки Eigen основана на предположении, что вся сложность, вся абстракция, разрешается на этапе компиляции. Это крайне важно: если абстракция не может быть разрешена на этапе компиляции, механизмы оптимизации Eigen на этапе компиляции становятся бесполезными, не говоря уже о том, что разрешение этой абстракции во время выполнения само по себе будет накладными затратами.
Здесь мы хотим иметь единственный класс MatrixBase в качестве базового класса для многих подклассов таким образом, что каждый объект MatrixBase (будь то матрица, вектор или любой вид выражения) знает на этапе компиляции (в отличие от времени выполнения), каким именно подклассом он является (то есть, является ли это матрицей, выражением и какого типа выражение).
Решение — шаблонный паттерн Curiously Recurring Template Pattern. Давайте сделаем перерыв сейчас. Надеюсь, вы сможете прочитать эту страницу Википедии во время перерыва, если это необходимо, но во время экзамена это будет запрещено.
Короче говоря, MatrixBase принимает шаблонный параметр Derived. Всякий раз, когда мы определяем подкласс Subclass, мы фактически заставляем Subclass наследовать MatrixBase<Subclass>. Суть в том, что разные подклассы наследуют разные типы MatrixBase. Благодаря этому, когда у нас есть объект подкласса и мы вызываем метод MatrixBase, мы всё ещё помним, какой именно подкласс мы обсуждаем, даже внутри метода MatrixBase.
Это означает, что мы можем поместить практически все методы и операторы в базовый класс MatrixBase, а в подклассах оставить только самое необходимое. Если вы посмотрите на подклассы в Eigen, например, класс CwiseBinaryOp, они имеют очень мало методов. Есть методы coeff() и иногда coeffRef() для доступа к коэффициентам, есть методы rows() и cols() для возврата количества строк и столбцов, но больше ничего нет. Вся сложность находится в MatrixBase, поэтому его нужно программировать только один раз для всех типов выражений, матриц и векторов.
Итак, давайте закончим это отступление и вернёмся к куску кода из нашего примера программы, который мы в данный момент анализируем,
v + w
Теперь, когда MatrixBase — наш хороший друг, давайте полностью запишем прототип оператора +, который вызывается здесь (этот код из src/Core/MatrixBase.h):
template<typename Derived> class MatrixBase { // ... template<typename OtherDerived> const CwiseBinaryOp<internal::scalar_sum_op<typename internal::traits<Derived>::Scalar>, Derived, OtherDerived> operator+(const MatrixBase<OtherDerived> &other) const; // ... };
Здесь, конечно, Derived и OtherDerived — это VectorXf.
Как мы уже сказали, CwiseBinaryOp также используется для других операций, таких как вычитание, поэтому он принимает ещё один шаблонный параметр, определяющий операцию, которая будет применена к коэффициентам. Этот шаблонный параметр — функтор, то есть класс, в котором у нас есть оператор (), поэтому он ведёт себя как функция. Здесь используемый функтор — internal::scalar_sum_op. Он определён в src/Core/Functors.h.
Теперь давайте объясним internal::traits. Класс internal::scalar_sum_op принимает один шаблонный параметр: тип чисел для обработки. Здесь, конечно, мы хотим передать скалярный тип (также известный как числовой тип) VectorXf, который равен float. Как определить скалярный тип Derived? Во всей библиотеке Eigen все типы матриц и выражений определяют typedef Scalar, который даёт его скалярный тип. Например, VectorXf::Scalar — это typedef для float. Итак, здесь, если бы жизнь была лёгкой, мы могли бы найти числовой тип Derived просто как
typename Derived::Scalar
К сожалению, мы не можем сделать это здесь, так как компилятор будет жаловаться, что тип Derived ещё не определён. Поэтому мы используем обходной путь: в src/Core/util/ForwardDeclarations.h мы объявили (но не определили!) все наши подклассы, такие как Matrix, и также объявили следующий шаблонный класс:
template<typename T> struct internal::traits;
В src/Core/Matrix.h, непосредственно перед определением класса Matrix, мы определяем частичную специализацию internal::traits для T=Matrix<любые шаблонные параметры>. В этой специализации internal::traits мы определяем typedef Scalar. Таким образом, когда мы фактически определяем Matrix, допустимо ссылаться на "typename internal::traits<Matrix>::Scalar".
В любом случае, мы объявили наш оператор +. В нашем случае, где Derived и OtherDerived — это VectorXf, вышеприведённое объявление эквивалентно:
class MatrixBase<VectorXf> { // ... const CwiseBinaryOp<internal::scalar_sum_op<float>, VectorXf, VectorXf> operator+(const MatrixBase<VectorXf> &other) const; // ... };
Теперь давайте перейдём к src/Core/CwiseBinaryOp.h, чтобы увидеть, как он определён. Как вы можете увидеть там, всё, что он делает, — возвращает объект CwiseBinaryOp, и этот объект просто хранит ссылки на выражения левой и правой части — здесь это векторы v и w. Ну, объект CwiseBinaryOp также хранит экземпляр класса функтора (пустой), но вы не должны об этом беспокоиться, так как это незначительный деталь реализации.
Таким образом, оператор + не выполнил никакого фактического вычисления. Подводя итог, операция v + w просто вернула объект типа CwiseBinaryOp, который ничего не делал, кроме хранения ссылок на v и w.
Присваивание
На этом этапе выражение v + w завершило вычисление, поэтому в процессе компиляции строки кода
u = v + w;
мы теперь попадаем в оператор =.
Какой оператор = вызывается здесь? Вектор u — это объект класса VectorXf, т.е. Matrix. В src/Core/Matrix.h, внутри определения класса Matrix, мы видим это:
template<typename OtherDerived> inline Matrix& operator=(const MatrixBase<OtherDerived>& other) { eigen_assert(m_storage.data()!=0 && "you cannot use operator= with a non initialized matrix (instead use set()"); return Base::operator=(other.derived()); }
Здесь Base — это typedef для MatrixBase<Matrix>. Итак, вызывается оператор = класса MatrixBase. Давайте посмотрим его прототип в src/Core/MatrixBase.h:
template<typename OtherDerived> Derived& operator=(const MatrixBase<OtherDerived>& other);
Здесь Derived — это VectorXf (поскольку u — это VectorXf), а OtherDerived — это CwiseBinaryOp. Более конкретно, как было объяснено в предыдущем разделе, OtherDerived — это:
CwiseBinaryOp<internal::scalar_sum_op<float>, VectorXf, VectorXf>
Итак, полный прототип вызываемого оператора = выглядит так:
VectorXf& MatrixBase<VectorXf>::operator=(const MatrixBase<CwiseBinaryOp<internal::scalar_sum_op<float>, VectorXf, VectorXf> > & other);
Этот оператор = буквально означает "копирование суммы двух VectorXf в другой VectorXf".
Теперь давайте посмотрим на реализацию этого оператора =. Она находится в файле src/Core/Assign.h.
Что мы можем там увидеть:
template<typename Derived> template<typename OtherDerived> inline Derived& MatrixBase<Derived> ::operator=(const MatrixBase<OtherDerived>& other) { return internal::assign_selector<Derived,OtherDerived>::run(derived(), other.derived()); }
Хорошо, наша следующая задача — понять internal::assign_selector :)
Вот его объявление (всё это по-прежнему в том же файле src/Core/Assign.h):
template<typename Derived, typename OtherDerived, bool EvalBeforeAssigning = int(OtherDerived::Flags) & EvalBeforeAssigningBit, bool NeedToTranspose = Derived::IsVectorAtCompileTime && OtherDerived::IsVectorAtCompileTime && int(Derived::RowsAtCompileTime) == int(OtherDerived::ColsAtCompileTime) && int(Derived::ColsAtCompileTime) == int(OtherDerived::RowsAtCompileTime) && int(Derived::SizeAtCompileTime) != 1> struct internal::assign_selector;
Итак, internal::assign_selector принимает 4 шаблонных параметра, но 2 последних автоматически определяются по 2 первым.
EvalBeforeAssigning призван обеспечить соблюдение EvalBeforeAssigningBit. Как объясняется здесь, определённые выражения имеют этот флаг, что заставляет их автоматически вычисляться в временные переменные перед их присваиванием другому выражению. Это относится к выражению Product, чтобы избежать странных эффектов алиасинга при выполнении "m = m * m;". Однако, разумеется, наше выражение CwiseBinaryOp не имеет флага EvalBeforeAssigningBit: мы с самого начала решили, что не хотим, чтобы здесь вводилась временная переменная. Поэтому, если вы перейдёте по ссылке src/Core/CwiseBinaryOp.h, вы увидите, что флаги в internal::traits<CwiseBinaryOp> не включают EvalBeforeAssigningBit. Член Flags выражения CwiseBinaryOp затем импортируется из internal::traits с помощью макроса EIGEN_GENERIC_PUBLIC_INTERFACE. В любом случае, здесь шаблонный параметр EvalBeforeAssigning имеет значение false.
NeedToTranspose нужен в случае, когда пользователь хочет скопировать строковый вектор в столбцовый вектор. Мы разрешаем это как исключение из общего правила, что при присваивании требуется соответствие размерностей. В любом случае, здесь как левая, так и правая части являются столбцовыми векторами, в том смысле, что ColsAtCompileTime равно 1. Таким образом, NeedToTranspose также имеет значение false.
Итак, здесь мы находимся в частичной специализации:
internal::assign_selector<Derived, OtherDerived, false, false>
Вот как это определено:
template<typename Derived, typename OtherDerived> struct internal::assign_selector<Derived,OtherDerived,false,false> { static Derived& run(Derived& dst, const OtherDerived& other) { return dst.lazyAssign(other.derived()); } };
Хорошо, теперь наша следующая задача – понять, как работает lazyAssign :)
template<typename Derived> template<typename OtherDerived> inline Derived& MatrixBase<Derived> ::lazyAssign(const MatrixBase<OtherDerived>& other) { EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(Derived,OtherDerived) eigen_assert(rows() == other.rows() && cols() == other.cols()); internal::assign_impl<Derived, OtherDerived>::run(derived(),other.derived()); return derived(); }
Что мы видим здесь? Некоторые утверждения, и единственная интересная строка:
internal::assign_impl<Derived, OtherDerived>::run(derived(),other.derived());
Хорошо, теперь мы хотим узнать, что находится внутри internal::assign_impl.
Вот его объявление:
template<typename Derived1, typename Derived2, int Vectorization = internal::assign_traits<Derived1, Derived2>::Vectorization, int Unrolling = internal::assign_traits<Derived1, Derived2>::Unrolling> struct internal::assign_impl;
Опять же, internal::assign_selector принимает 4 шаблонных параметра, но 2 последних автоматически определяются первыми 2.
Эти два параметра Vectorization и Unrolling определяются вспомогательным классом internal::assign_traits. Его задача – определить, какую стратегию векторизации использовать (т. е. Vectorization) и какую стратегию развёртывания (т. е. Unrolling).
Мы не будем вдаваться в подробности о том, как выбираются эти стратегии (это реализовано в internal::assign_traits в верхней части того же файла). Скажем только, что здесь Vectorization имеет значение LinearVectorization, а Unrolling имеет значение NoUnrolling (последнее очевидно, так как наши векторы имеют динамический размер, поэтому развёртывание цикла во время компиляции невозможно).
Таким образом, частичная специализация internal::assign_impl, которую мы рассматриваем, выглядит так:
internal::assign_impl<Derived1, Derived2, LinearVectorization, NoUnrolling>
Вот как она определена:
template<typename Derived1, typename Derived2> struct internal::assign_impl<Derived1, Derived2, LinearVectorization, NoUnrolling> { static void run(Derived1 &dst, const Derived2 &src) { const int size = dst.size(); const int packetSize = internal::packet_traits<typename Derived1::Scalar>::size; const int alignedStart = internal::assign_traits<Derived1,Derived2>::DstIsAligned ? 0 : internal::first_aligned(&dst.coeffRef(0), size); const int alignedEnd = alignedStart + ((size-alignedStart)/packetSize)*packetSize; for(int index = 0; index < alignedStart; index++) dst.copyCoeff(index, src); for(int index = alignedStart; index < alignedEnd; index += packetSize) { dst.template copyPacket<Derived2, Aligned, internal::assign_traits<Derived1,Derived2>::SrcAlignment>(index, src); } for(int index = alignedEnd; index < size; index++) dst.copyCoeff(index, src); } };
Вот как это работает. LinearVectorization означает, что к выражениям левой и правой стороны можно обратиться линейно, т. е. вы можете ссылаться на их коэффициенты с помощью одного целого числа index, в отличие от необходимости ссылаться на их коэффициенты с помощью двух целых чисел row, column.
Как мы сказали в начале, векторизация работает с блоками из 4 чисел с плавающей точкой. Здесь PacketSize равно 4.
Существует две потенциальные проблемы, с которыми нам нужно справиться:
- во-первых, векторизация работает намного лучше, если пакеты выровнены по 128 битам. Это особенно важно для записи. Поэтому при записи в коэффициенты dst мы хотим сгруппировать эти коэффициенты по пакетам по 4 так, чтобы каждый из этих пакетов был выровнен по 128 битам. В общем случае это требует пропустить несколько коэффициентов в начале dst. Именно для этого предназначено alignedStart. Затем мы копируем эти первые несколько коэффициентов по одному, а не пакетами. Однако в нашем случае выражение dst – это VectorXf, и помните, что при создании векторов мы выделили выровненные массивы. Благодаря DstIsAligned, Eigen запоминает это, не выполняя никаких проверок во время выполнения, поэтому alignedStart равно нулю, и эта часть полностью избегается.
- во-вторых, количество копируемых коэффициентов не всегда является кратным packetSize. Здесь нужно скопировать 50 коэффициентов, а packetSize равно 4. Таким образом, нам придётся скопировать последние 2 коэффициента по одному, а не пакетами. Здесь alignedEnd равно 48.
Теперь приходят фактические циклы.
Сначала векторизованная часть: первые 48 коэффициентов из 50 будут скопированы пакетами по 4:
for(int index = alignedStart; index < alignedEnd; index += packetSize) { dst.template copyPacket<Derived2, Aligned, internal::assign_traits<Derived1,Derived2>::SrcAlignment>(index, src); }
Что такое copyPacket? Оно определено в src/Core/Coeffs.h:
template<typename Derived> template<typename OtherDerived, int StoreMode, int LoadMode> inline void MatrixBase<Derived>::copyPacket(int index, const MatrixBase<OtherDerived>& other) { eigen_internal_assert(index >= 0 && index < size()); derived().template writePacket<StoreMode>(index, other.derived().template packet<LoadMode>(index)); }
Хорошо, что такое writePacket() и packet() здесь?
Во-первых, writePacket() здесь – это метод левой стороны VectorXf. Поэтому мы переходим по ссылке src/Core/Matrix.h, чтобы посмотреть его определение:
template<int StoreMode> inline void writePacket(int index, const PacketScalar& x) { internal::pstoret<Scalar, PacketScalar, StoreMode>(m_storage.data() + index, x); }
Здесь StoreMode – это Aligned, указывающее, что мы выполняем запись с выравниванием по 128 битам, PacketScalar – это тип, представляющий "пакет SSE из 4 чисел с плавающей точкой", и internal::pstoret – это функция, записывающая такой пакет в память. Их определения зависят от архитектуры, мы найдём их в src/Core/arch/SSE/PacketMath.h:
Строка в src/Core/arch/SSE/PacketMath.h, определяющая тип PacketScalar (через typedef в Matrix.h), выглядит так:
template<> struct internal::packet_traits<float> { typedef __m128 type; enum {size=4}; };
Здесь __m128 – это специфичный для SSE тип. Обратите внимание, что перечисление size здесь использовалось для определения packetSize выше.
И вот реализация internal::pstoret:
template<> inline void internal::pstore(float* to, const __m128& from) { _mm_store_ps(to, from); }
Здесь __mm_store_ps – это специфичная для SSE встроенная функция, представляющая собой одну инструкцию SSE. Разница между internal::pstore и internal::pstoret заключается в том, что internal::pstoret – это диспетчер, обрабатывающий как выровненные, так и невыровненные случаи. Вы найдёте его определение в src/Core/GenericPacketMath.h:
template<typename Scalar, typename Packet, int LoadMode> inline void internal::pstoret(Scalar* to, const Packet& from) { if(LoadMode == Aligned) internal::pstore(to, from); else internal::pstoreu(to, from); }
Хорошо, это объясняет, как работает writePacket(). Теперь давайте посмотрим на вызов packet(). Помните, что мы анализируем эту строку кода внутри copyPacket():
derived().template writePacket<StoreMode>(index,
other.derived().template packet<LoadMode>(index));
Здесь other – это наше выражение суммы v + w. .derived() просто преобразует MatrixBase в подкласс, который здесь является CwiseBinaryOp. Поэтому давайте перейдём по ссылке src/Core/CwiseBinaryOp.h:
class CwiseBinaryOp { // ... template<int LoadMode> inline PacketScalar packet(int index) const { return m_functor.packetOp(m_lhs.template packet<LoadMode>(index), m_rhs.template packet<LoadMode>(index)); } };
Здесь m_lhs – это вектор v, а m_rhs – вектор w. Таким образом, функция packet() здесь – это Matrix::packet(). Шаблонный параметр LoadMode – это Aligned. Поэтому мы рассматриваем
class Matrix { // ... template<int LoadMode> inline PacketScalar packet(int index) const { return internal::ploadt<Scalar, LoadMode>(m_storage.data() + index); } };
Мы оставляем вам изучить определение internal::ploadt в GenericPacketMath.h и internal::pload в src/Core/arch/SSE/PacketMath.h. Оно очень похоже на предыдущее для internal::pstore.
Вернёмся к CwiseBinaryOp::packet(). После того, как пакеты из векторов v и w были возвращены, что делает эта функция? Она вызывает m_functor.packetOp() на них. Что такое m_functor? Здесь мы должны помнить, с какой конкретной шаблонной специализацией CwiseBinaryOp мы имеем дело:
CwiseBinaryOp<internal::scalar_sum_op<float>, VectorXf, VectorXf>
Таким образом, m_functor – это объект пустого класса internal::scalar_sum_op<float>. Как мы упоминали выше, не беспокойтесь о том, почему вообще был создан объект этого пустого класса – это деталь реализации, суть в том, что некоторые другие функторы должны хранить данные членов.
В любом случае, internal::scalar_sum_op определена в src/Core/Functors.h:
template<typename Scalar> struct internal::scalar_sum_op EIGEN_EMPTY_STRUCT { inline const Scalar operator() (const Scalar& a, const Scalar& b) const { return a + b; } template<typename PacketScalar> inline const PacketScalar packetOp(const PacketScalar& a, const PacketScalar& b) const { return internal::padd(a,b); } };
Как вы можете видеть, всё, что делает packetOp(), это вызов internal::padd на двух пакетах. Вот определение internal::padd из src/Core/arch/SSE/PacketMath.h:
template<> inline __m128 internal::padd(const __m128& a, const __m128& b) { return _mm_add_ps(a,b); }
Здесь _mm_add_ps – это специфичная для SSE встроенная функция, представляющая собой одну инструкцию SSE.
Подводя итог, цикл
for(int index = alignedStart; index < alignedEnd; index += packetSize) { dst.template copyPacket<Derived2, Aligned, internal::assign_traits<Derived1,Derived2>::SrcAlignment>(index, src); }
был скомпилирован в следующий код: для index от 0 до 11 ( = 48/4 - 1), считывайте i-й пакет (из 4 чисел с плавающей точкой) из вектора v и i-й пакет из вектора w с помощью двух инструкций SSE __mm_load_ps, затем складывайте их с помощью инструкции __mm_add_ps, затем записывайте результат с помощью инструкции __mm_store_ps.
Остаётся второй цикл, обрабатывающий последние несколько (здесь последние 2) коэффициента:
for(int index = alignedEnd; index < size; index++) dst.copyCoeff(index, src);
Однако он работает так же, как и тот, который мы только что объяснили, он просто проще, потому что здесь нет векторизации SSE. copyPacket() становится copyCoeff(), packet() – coeff(), writePacket() – coeffRef(). Если вы дошли до этого момента, вы, вероятно, сможете понять эту часть самостоятельно.
Мы видим, что вся абстракция C++ в Eigen исчезает во время компиляции, и что мы действительно точно контролируем, какие инструкции ассемблера мы генерируем. Такова прелесть C++! Поскольку у нас есть такой точный контроль над сгенерированными инструкциями ассемблера, но такая сложная логика для выбора правильных инструкций, мы можем сказать, что Eigen действительно ведёт себя как оптимизирующая компилятор. Если хотите, можете сказать, что Eigen ведёт себя как скрипт для компилятора. В некотором смысле, метапрограммирование шаблонов C++ – это скриптинг компилятора – и было показано, что этот язык сценариев является Turing-полным. См. Wikipedia.
© Eigen.
Licensed under the MPL2 License.
https://eigen.tuxfamily.org/dox/TopicInsideEigenExample.html
