Шаблоны выражений и рукописный код

Question

В настоящее время я пишу библиотеку выражений шаблонов C++ и сравниваю некоторые экземпляры с рукописным кодом на уровне сборки. Рукописная функция заключается в следующем:

spinor multiply(vector const& a, vector const& b)
{
        spinor result = {
                a.at<1>() * b.at<1>() - a.at<2>() * b.at<2>()
                          - a.at<4>() * b.at<4>() - a.at<8>() * b.at<8>(),
                a.at<1>() * b.at<2>() - a.at<2>() * b.at<1>(),
                a.at<1>() * b.at<4>() - a.at<4>() * b.at<1>(),
                a.at<1>() * b.at<8>() - a.at<8>() * b.at<1>(),
                a.at<2>() * b.at<4>() - a.at<4>() * b.at<2>(),
                a.at<2>() * b.at<8>() - a.at<8>() * b.at<2>(),
                a.at<4>() * b.at<8>() - a.at<8>() * b.at<4>()
        };

        return result;
}

vector класс это просто обертка над четырьмя двойными числами, которые можно прочитать с помощью at<index>() функция-член. Из-за проектных решений индексы для четырех компонентов 1, 2, 4, 8 которые доступны с at<index>() вместо обычного 0, 1, 2, 3,

Цель этой функции - вернуть результат умножения двух векторов (в пространстве Минковского). Если вы знакомы с геометрической алгеброй, вы увидите точечный продукт (первый компонент resultсимметричный при обмене a а также b) и клин-продукт (остальные компоненты, антисимметричные при обмене a а также b). Если вы не знакомы с геометрической алгеброй, просто используйте эту функцию как рецепт для умножения векторов.

Если я скомпилирую вышеупомянутую функцию с GCC 4.7 и посмотрю на разборку, objdump -SC a.out это дает мне следующий вывод:

400bc0: movsd  0x8(%rsi),%xmm6
400bc5: mov    %rdi,%rax
400bc8: movsd  (%rsi),%xmm8
400bcd: movsd  0x8(%rdx),%xmm5
400bd2: movapd %xmm6,%xmm9
400bd7: movsd  (%rdx),%xmm7
400bdb: movapd %xmm8,%xmm0
400be0: mulsd  %xmm5,%xmm9
400be5: movsd  0x10(%rsi),%xmm4
400bea: mulsd  %xmm7,%xmm0
400bee: movsd  0x10(%rdx),%xmm1
400bf3: movsd  0x18(%rdx),%xmm3
400bf8: movsd  0x18(%rsi),%xmm2
400bfd: subsd  %xmm9,%xmm0
400c02: movapd %xmm4,%xmm9
400c07: mulsd  %xmm1,%xmm9
400c0c: subsd  %xmm9,%xmm0
400c11: movapd %xmm3,%xmm9
400c16: mulsd  %xmm2,%xmm9
400c1b: subsd  %xmm9,%xmm0
400c20: movapd %xmm6,%xmm9
400c25: mulsd  %xmm7,%xmm9
400c2a: movsd  %xmm0,(%rdi)
400c2e: movapd %xmm5,%xmm0
400c32: mulsd  %xmm8,%xmm0
400c37: subsd  %xmm9,%xmm0
400c3c: movapd %xmm4,%xmm9
400c41: mulsd  %xmm7,%xmm9
400c46: mulsd  %xmm2,%xmm7
400c4a: movsd  %xmm0,0x8(%rdi)
400c4f: movapd %xmm1,%xmm0
400c53: mulsd  %xmm8,%xmm0
400c58: mulsd  %xmm3,%xmm8
400c5d: subsd  %xmm9,%xmm0
400c62: subsd  %xmm7,%xmm8
400c67: movapd %xmm4,%xmm7
400c6b: mulsd  %xmm5,%xmm7
400c6f: movsd  %xmm0,0x10(%rdi)
400c74: mulsd  %xmm2,%xmm5
400c78: movapd %xmm1,%xmm0
400c7c: mulsd  %xmm6,%xmm0
400c80: movsd  %xmm8,0x18(%rdi)
400c86: mulsd  %xmm3,%xmm6
400c8a: mulsd  %xmm2,%xmm1
400c8e: mulsd  %xmm4,%xmm3
400c92: subsd  %xmm7,%xmm0
400c96: subsd  %xmm5,%xmm6
400c9a: subsd  %xmm1,%xmm3
400c9e: movsd  %xmm0,0x20(%rdi)
400ca3: movsd  %xmm6,0x28(%rdi)
400ca8: movsd  %xmm3,0x30(%rdi)
400cad: retq   
400cae: nop
400caf: nop

Это выглядит довольно хорошо для меня - компоненты первого (%rsi) и второе (%rdx) векторы доступны только один раз, а фактические вычисления выполняются только в регистрах. В конце результат записывается по адресу в реестре %rdi, Поскольку это первый регистр аргументов, я думаю, здесь используется оптимизация возвращаемого значения.

Сравните это со следующим листингом для версии шаблона выражения для функции выше:

400cb0: mov    (%rsi),%rdx
400cb3: mov    0x8(%rsi),%rax
400cb7: movsd  0x1f1(%rip),%xmm4        # 400eb0 <_IO_stdin_used+0x10>
400cbe: 
400cbf: movsd  0x10(%rdx),%xmm3
400cc4: movsd  0x18(%rdx),%xmm0
400cc9: mulsd  0x10(%rax),%xmm3
400cce: xorpd  %xmm4,%xmm0
400cd2: mulsd  0x18(%rax),%xmm0
400cd7: movsd  0x8(%rdx),%xmm2
400cdc: movsd  (%rdx),%xmm1
400ce0: mulsd  0x8(%rax),%xmm2
400ce5: mulsd  (%rax),%xmm1
400ce9: subsd  %xmm3,%xmm0
400ced: subsd  %xmm2,%xmm0
400cf1: addsd  %xmm0,%xmm1
400cf5: movsd  %xmm1,(%rdi)
400cf9: movsd  (%rdx),%xmm0
400cfd: movsd  0x8(%rdx),%xmm1
400d02: mulsd  0x8(%rax),%xmm0
400d07: mulsd  (%rax),%xmm1
400d0b: subsd  %xmm1,%xmm0
400d0f: movsd  %xmm0,0x8(%rdi)
400d14: movsd  (%rdx),%xmm0
400d18: movsd  0x10(%rdx),%xmm1
400d1d: mulsd  0x10(%rax),%xmm0
400d22: mulsd  (%rax),%xmm1
400d26: subsd  %xmm1,%xmm0
400d2a: movsd  %xmm0,0x10(%rdi)
400d2f: movsd  0x8(%rdx),%xmm0
400d34: movsd  0x10(%rdx),%xmm1
400d39: mulsd  0x10(%rax),%xmm0
400d3e: mulsd  0x8(%rax),%xmm1
400d43: subsd  %xmm1,%xmm0
400d47: movsd  %xmm0,0x18(%rdi)
400d4c: movsd  (%rdx),%xmm0
400d50: movsd  0x18(%rdx),%xmm1
400d55: mulsd  0x18(%rax),%xmm0
400d5a: mulsd  (%rax),%xmm1
400d5e: subsd  %xmm1,%xmm0
400d62: movsd  %xmm0,0x20(%rdi)
400d67: movsd  0x8(%rdx),%xmm0
400d6c: movsd  0x18(%rdx),%xmm1
400d71: mulsd  0x18(%rax),%xmm0
400d76: mulsd  0x8(%rax),%xmm1
400d7b: subsd  %xmm1,%xmm0
400d7f: movsd  %xmm0,0x28(%rdi)
400d84: movsd  0x10(%rdx),%xmm0
400d89: movsd  0x18(%rdx),%xmm1
400d8e: mulsd  0x18(%rax),%xmm0
400d93: mulsd  0x10(%rax),%xmm1
400d98: subsd  %xmm1,%xmm0
400d9c: movsd  %xmm0,0x30(%rdi)
400da1: retq   

Сигнатура этой функции

spinor<product<vector, vector>>(product<vector, vector> const&)

Надеюсь, вы мне доверяете, что обе версии дают одинаковый результат. Первые две строки извлекают первый и второй вектор, которые хранятся в виде ссылок в product, Я задавался вопросом о следующих вещах:

• Что значит movsd 0x1f1(%rip),%xmm4 в комбинации с xorpd %xmm4,%xmm0 делать? Я уже выяснил, что это называется "относительная адресация RIP", см. http://www.x86-64.org/documentation/assembly.html
• Почему GCC не использует больше регистров, например, для кэширования 0x10(%rax) который читается четыре раза?

Я также протестировал обе функции, сгенерировав 100000000 случайных векторов и потратив время на обе функции:

ET: 7.5 sec
HW: 6.8 sec

Рукописная функция примерно на 10% быстрее. Кто-нибудь имеет опыт работы с шаблонами выражений и знает, как заставить их работать ближе к своему рукописному аналогу?

Answer 1

Было бы понятно, если бы мы точно знали содержание адреса 0x400eb0, но я подозреваю, что это 0x8000 0000 0000 0000 8000 0000 0000 0000или аналогичный (возможно, с начальным 0, потому что код не векторизован), записанный как 128-битное целое число.

В этом случае xorpd меняет знак второго операнда

Почему чтение реестра не кэшируется - лучше спросите об этом в списке рассылки gcc-help. Возможно, компилятор не может доказать, что два вектора или промежуточный результат не являются псевдонимами.

Но вопреки общему мнению, компиляторы не всегда полностью оптимизируют, а только лучше, чем 90% (или 99%?) Всех программистов (если они пытаются писать на ассемблере), а иногда (редко) они создают действительно медленный код.

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

PS: код может быть ускорен с помощью векторных инструкций (mulpd вместо mulsd), который умножает два или четыре двойных за один раз), он же SSE или AVX. Но для перестановки значений в нужные места в регистрах необходимы некоторые инструкции, поэтому усиление всегда медленнее, чем в два или четыре раза.