Ответ 1
Вот пример с помощью Eigen library:
#include <Eigen/Core>
std::complex<float> f(const std::complex<float> *x, int n)
{
return Eigen::VectorXcf::Map(x, n).prod();
}
Если вы скомпилируете это с помощью clang или g++ и sse или avx enabled (и -O2), вы должны получить довольно приличный машинный код. Он также работает для некоторых других архитектур, таких как Altivec или NEON. Если вы знаете, что первая запись x выровнена, вы можете использовать MapAligned вместо Map.
Вы получаете еще лучший код, если вы случайно знаете размер вашего вектора при компиляции:
template<int n>
std::complex<float> f(const std::complex<float> *x)
{
return Eigen::Matrix<std::complex<float>, n, 1> >::MapAligned(x).prod();
}
Примечание. Вышеуказанные функции непосредственно соответствуют функции f OP.
Однако, как указывал @PeterCordes, обычно сложно хранить сложные числа, чередующиеся, так как это потребует много перетасовки для умножения. Вместо этого нужно хранить реальные и мнимые части так, чтобы они могли сразу загружать один пакет за один раз.
Изменить/добавить. Чтобы реализовать структуру массивов, например комплексное умножение, вы можете написать что-то вроде:
typedef Eigen::Array<float, 8, 1> v8sf; // Eigen::Array allows element-wise standard operations
typedef std::complex<v8sf> complex8;
complex8 prod(const complex8& a, const complex8& b)
{
return a*b;
}
Или более общий (с использованием С++ 11):
template<int size, typename Scalar = float> using complexX = std::complex<Eigen::Array<Scalar, size, 1> >;
template<int size>
complexX<size> prod(const complexX<size>& a, const complexX<size>& b)
{
return a*b;
}
При компиляции с помощью -mavx -O2 это скомпилируется примерно так (с помощью g++ - 5.4):
vmovaps 32(%rsi), %ymm1
movq %rdi, %rax
vmovaps (%rsi), %ymm0
vmovaps 32(%rdi), %ymm3
vmovaps (%rdi), %ymm4
vmulps %ymm0, %ymm3, %ymm2
vmulps %ymm4, %ymm1, %ymm5
vmulps %ymm4, %ymm0, %ymm0
vmulps %ymm3, %ymm1, %ymm1
vaddps %ymm5, %ymm2, %ymm2
vsubps %ymm1, %ymm0, %ymm0
vmovaps %ymm2, 32(%rdi)
vmovaps %ymm0, (%rdi)
vzeroupper
ret
По причинам, не очевидным для меня, это фактически скрыто в методе, который вызывается фактическим методом, который просто перемещается вокруг некоторой памяти - я не знаю, почему Eigen/gcc не предполагает, что аргументы уже правильно выровнен. Если я скомпилирую то же самое с clang 3.8.0 (и теми же аргументами), он скомпилируется только:
vmovaps (%rsi), %ymm0
vmovaps %ymm0, (%rdi)
vmovaps 32(%rsi), %ymm0
vmovaps %ymm0, 32(%rdi)
vmovaps (%rdi), %ymm1
vmovaps (%rdx), %ymm2
vmovaps 32(%rdx), %ymm3
vmulps %ymm2, %ymm1, %ymm4
vmulps %ymm3, %ymm0, %ymm5
vsubps %ymm5, %ymm4, %ymm4
vmulps %ymm3, %ymm1, %ymm1
vmulps %ymm0, %ymm2, %ymm0
vaddps %ymm1, %ymm0, %ymm0
vmovaps %ymm0, 32(%rdi)
vmovaps %ymm4, (%rdi)
movq %rdi, %rax
vzeroupper
retq
Опять же, движение памяти в начале является странным, но, по крайней мере, оно векторизовано. Однако для gcc и clang это оптимизируется при вызове в цикле:
complex8 f8(complex8 x[], int n) {
if(n==0)
return complex8(v8sf::Ones(),v8sf::Zero()); // I guess you want p = 1 + 0*i at the beginning?
complex8 p = x[0];
for (int i = 1; i < n; i++) p = prod(p, x[i]);
return p;
}
Различие заключается в том, что clang будет разворачивать этот внешний цикл на 2 умножения на цикл. С другой стороны, при компиляции с -mfma gcc будет использовать команды с плавным перемножением-добавлением.
Конечно, функция f8 также может быть обобщена на произвольные размеры:
template<int size>
complexX<size> fX(complexX<size> x[], int n) {
using S= typename complexX<size>::value_type;
if(n==0)
return complexX<size>(S::Ones(),S::Zero());
complexX<size> p = x[0];
for (int i = 1; i < n; i++) p *=x[i];
return p;
}
И для уменьшения complexX<N> до одного std::complex можно использовать следующую функцию:
// only works for powers of two
template<int size> EIGEN_ALWAYS_INLINE
std::complex<float> redux(const complexX<size>& var) {
complexX<size/2> a(var.real().template head<size/2>(), var.imag().template head<size/2>());
complexX<size/2> b(var.real().template tail<size/2>(), var.imag().template tail<size/2>());
return redux(a*b);
}
template<> EIGEN_ALWAYS_INLINE
std::complex<float> redux(const complexX<1>& var) {
return std::complex<float>(var.real()[0], var.imag()[0]);
}
Однако, в зависимости от того, использую ли я clang или g++, я получаю совсем другой выход ассемблера. В целом, g++ имеет тенденцию к сбою встроенной загрузки входных аргументов, а clang не использует операции FMA (YMMV...) По сути, вам все равно нужно проверить сгенерированный код ассемблера. И что еще более важно, вам следует сравнить код (не уверен, насколько сильно эта процедура работает в вашей общей проблеме).
Кроме того, я хотел бы отметить, что Eigen фактически является библиотекой линейной алгебры. Использование его для создания чистого портативного кода SIMD на самом деле не предназначено для.