Оптимизация 2D вращения

Дана классическая формула для вращения точки в двумерном пространстве: 

cv::Point pt[NPOINTS];
cv::Point rotated[NPOINTS];
float angle = WHATEVER;
float cosine = cos(angle);
float sine = sin(angle);

for (int i = 0; i < NPOINTS; i++)
{
    rotated[i].x = pt[i].x * cosine - pt[i].y * sine;
    rotated[i].y = pt[i].x * sine   + pt[i].y * cosine;
}

Учитывая, что N точек равно 32 и массивы выровнены, как можно оптимизировать код для SSE или AVX? Поиск здесь и в других местах не нашел ничего полезного, и я заблудился здесь: 

__m128i onePoint = _mm_set_epi32(pt[i].x, pt[i].y, pt[i].x, pt[i].y);
__m128 onefPoint = _m128_cvtepi32_ps(onePoint);
__m128 sinCos = _mm_set_ps(cosine, -sine, sine, cosine);
__m128 rotated = _mm_mul_ps(onefPoint, sinCos);
Но как перейти от [y*cosine, -x*sine, x*sine, y*cosine] к [y*cosine + -x*sine, x*sine + y*cosine]? Является ли это лучшим подходом? Легко ли он масштабируется до __m512?

UPDATE : я провел еще немного исследований, и теперь у меня есть приблизительно: 

__m128i onePoint = _mm_set_epi32(pt[i].x, pt[i].y, pt[i].x, pt[i].y);
__m128 onefPoint = _m128_cvtepi32_ps(onePoint);
__m128i twoPoint = _mm_set_epi32(pt[i+1].x, pt[i+1].y, pt[i+1].x, pt[i+1].y);
__m128 twofPoint = _m128_cvtepi32_ps(twoPoint);
__m128 sinCos = _mm_set_ps(cosine, -sine, sine, cosine);
__m128 rotated1 = _mm_mul_ps(onefPoint, sinCos);
__m128 rotated2 = _mm_mul_ps(twofPoint, sinCos);
__m128 added = _mm_hadd_ps(rotated1, rotated2);
__m128i intResult = _mm_cvtps_epi32(added);
int results[4];
_mm_storeu_si128((__m128i*)results, intResult);
Это дает 50% - ную скорость-от 11% процессорного времени до примерно 6%. Расширение до __m256 и выполнение четырех точек одновременно дает еще одно ускорение. Это выглядит довольно ужасно, но я иду в правильном направлении?
1 3
c++ opencv optimization avx

1 ответ:

Используйте массив struct of arrays (AoSoA) и обрабатывайте восемь точек одновременно. В приведенном ниже коде point8 есть структура массивов, содержащих восемь точек. Функция rotate_point8 вращает восемь точек и имеет ту же алгебраическую структуру, что и функция rotate_point, которая вращает одну точку. Функция rotate_all8 вращает 32 точки, используя AoSoA point8*.

Одноточечный код вращения выполняет 4 умножения, одно сложение и одно вычитание.

Если мы посмотрим на сборку для rotate_point8 мы видим, что GCC развернул цикл и делает 4 умножения SIMD, одно сложение SIMD и одно вычитание SIMD за размотку. Это лучшее, что вы можете сделать: восемь по цене одного. 

#include <x86intrin.h>
#include <stdio.h>
#include <math.h>

struct point8 {
  __m256 x;
  __m256 y;
};

struct point {
  float x;
  float y;
};

static point rotate_point(point p, float a, float b) {
  point r;
  r.x = p.x*a - p.y*b;
  r.y = p.x*b + p.y*a;
  return r;
}

static point8 rotate_point8(point8 p, float a, float b) {
  __m256 va = _mm256_set1_ps(a), vb = _mm256_set1_ps(b);
  point8 r;
  r.x = _mm256_sub_ps(_mm256_mul_ps(p.x,va), _mm256_mul_ps(p.y,vb));
  r.y = _mm256_add_ps(_mm256_mul_ps(p.x,vb), _mm256_mul_ps(p.y,va));
  return r;
}

void rotate_all(point* points, point* r, float angle) {
  float a = cos(angle), b = sin(angle);
  for(int i=0; i<32; i++) r[i] = rotate_point(points[i], a, b);
}

void rotate_all8(point8* points, point8* r8, float angle) {
  float a = cos(angle), b = sin(angle);
  for(int i=0; i<4; i++) r8[i] = rotate_point8(points[i], a, b);
}

int main(void) {
  float x[32], y[32];
  point p[32], r[32];
  point8 p8[4], r8[4];
  float angle = 3.14159f/4;

  for(int i=0; i<32; i++) y[i] = 1.0*i/31, x[i] = sqrt(1-y[i]*y[i]);
  for(int i=0; i<32; i++) p[i].x = x[i], p[i].y = y[i];
  for(int i=0; i<4; i++) p8[i].x = _mm256_load_ps(&x[8*i]), p8[i].y = _mm256_load_ps(&y[8*i]); 

  for(int i=0; i<32; i++) printf("%f %f\n", p[i].x, p[i].y); puts("");

  rotate_all(p, r, angle);
  for(int i=0; i<32; i++) printf("%f %f\n", r[i].x, r[i].y); puts("");

  rotate_all8(p8, r8, angle);
  for(int i=0; i<4; i++) {
    _mm256_storeu_ps(x, r8[i].x), _mm256_storeu_ps(y, r8[i].y);
    for(int j=0; j<8; j++) printf("%f %f\n", x[j], y[j]);
  }
}
1