Введение в проблему
Я пытаюсь ускорить код пересечения (2d) трассировщика лучей, который я пишу. Я использую C# и библиотеку System.Numerics, чтобы увеличить скорость SIMD-инструкций.
Проблема в том, что я получаю странные результаты, с запредельным ускорением и довольно низким ускорением. У меня вопрос, почему один над крышей, а другой довольно низкий?
Контекст:
- Структура RayPack представляет собой серию (разных) лучей, упакованных в векторы System.Numerics.
- Структура BoundingBoxPack и CirclePack представляет собой единый bb/circle, упакованный в векторы System.Numerics.
- Используемый процессор — i7-4710HQ (Haswell) с инструкциями SSE 4.2, AVX(2) и FMA(3).
- Работает в режиме выпуска (64 бит). Проект работает в .Net Framework 472. Никаких дополнительных опций не установлено.
Попытки
Я уже пытался выяснить, могут ли некоторые операции поддерживаться должным образом или нет (обратите внимание: это для С++. https://fgiesen.wordpress.com/2016/04/03/sse-mind-the-gap/ или http://sci.tuomastonteri.fi/programming/sse), но, похоже, это не так, потому что ноутбук, на котором я работаю, поддерживает ССЕ 4.2.
В текущем коде применяются следующие изменения:
- Используя более правильные инструкции (например, упакованный мин).
- Неиспользование инструкции float * vector (вызывает множество дополнительных операций, см. сборку оригинала).
Код…фрагменты?
Извиняюсь за большой объем кода, но я не уверен, как мы можем обсудить это конкретно без такого количества кода.
Код луча -> BoundingBox
public bool CollidesWith(Ray ray, out float t)
{
// https://gamedev.stackexchange.com/questions/18436/most-efficient-aabb-vs-ray-collision-algorithms
// r.dir is unit direction vector of ray
float dirfracx = 1.0f / ray.direction.X;
float dirfracy = 1.0f / ray.direction.Y;
// lb is the corner of AABB with minimal coordinates - left bottom, rt is maximal corner
// r.org is origin of ray
float t1 = (this.rx.min - ray.origin.X) * dirfracx;
float t2 = (this.rx.max - ray.origin.X) * dirfracx;
float t3 = (this.ry.min - ray.origin.Y) * dirfracy;
float t4 = (this.ry.max - ray.origin.Y) * dirfracy;
float tmin = Math.Max(Math.Min(t1, t2), Math.Min(t3, t4));
float tmax = Math.Min(Math.Max(t1, t2), Math.Max(t3, t4));
// if tmax < 0, ray (line) is intersecting AABB, but the whole AABB is behind us
if (tmax < 0)
{
t = tmax;
return false;
}
// if tmin > tmax, ray doesn't intersect AABB
if (tmin > tmax)
{
t = tmax;
return false;
}
t = tmin;
return true;
}
Код RayPack -> BoundingBoxPack
public Vector<int> CollidesWith(ref RayPack ray, out Vector<float> t)
{
// ------------------------------------------------------- \\
// compute the collision. \\
Vector<float> dirfracx = Constants.ones / ray.direction.x;
Vector<float> dirfracy = Constants.ones / ray.direction.y;
Vector<float> t1 = (this.rx.min - ray.origin.x) * dirfracx;
Vector<float> t2 = (this.rx.max - ray.origin.x) * dirfracx;
Vector<float> t3 = (this.ry.min - ray.origin.y) * dirfracy;
Vector<float> t4 = (this.ry.max - ray.origin.y) * dirfracy;
Vector<float> tmin = Vector.Max(Vector.Min(t1, t2), Vector.Min(t3, t4));
Vector<float> tmax = Vector.Min(Vector.Max(t1, t2), Vector.Max(t3, t4));
Vector<int> lessThanZeroMask = Vector.GreaterThan(tmax, Constants.zeros);
Vector<int> greaterMask = Vector.LessThan(tmin, tmax);
Vector<int> combinedMask = Vector.BitwiseOr(lessThanZeroMask, greaterMask);
// ------------------------------------------------------- \\
// Keep track of the t's that collided. \\
t = Vector.ConditionalSelect(combinedMask, tmin, Constants.floatMax);
return combinedMask;
}
Код Луча -> Круг
public bool Intersect(Circle other)
{
// Step 0: Work everything out on paper!
// Step 1: Gather all the relevant data.
float ox = this.origin.X;
float dx = this.direction.X;
float oy = this.origin.Y;
float dy = this.direction.Y;
float x0 = other.origin.X;
float y0 = other.origin.Y;
float cr = other.radius;
// Step 2: compute the substitutions.
float p = ox - x0;
float q = oy - y0;
float r = 2 * p * dx;
float s = 2 * q * dy;
// Step 3: compute the substitutions, check if there is a collision.
float a = dx * dx + dy * dy;
float b = r + s;
float c = p * p + q * q - cr * cr;
float DSqrt = b * b - 4 * a * c;
// no collision possible! Commented out to make the benchmark more fair
//if (DSqrt < 0)
//{ return false; }
// Step 4: compute the substitutions.
float D = (float)Math.Sqrt(DSqrt);
float t0 = (-b + D) / (2 * a);
float t1 = (-b - D) / (2 * a);
float ti = Math.Min(t0, t1);
if(ti > 0 && ti < t)
{
t = ti;
return true;
}
return false;
}
Код RayPack -> CirclePack Оригинальный, неотредактированный код можно найти по адресу: https://pastebin.com/87nYgZrv
public Vector<int> Intersect(CirclePack other)
{
// ------------------------------------------------------- \\
// Put all the data on the stack. \\
Vector<float> zeros = Constants.zeros;
Vector<float> twos = Constants.twos;
Vector<float> fours = Constants.fours;
// ------------------------------------------------------- \\
// Compute whether the ray collides with the circle. This \\
// is stored in the 'mask' vector. \\
Vector<float> p = this.origin.x - other.origin.x; ;
Vector<float> q = this.origin.y - other.origin.y;
Vector<float> r = twos * p * this.direction.x;
Vector<float> s = twos * q * this.direction.y; ;
Vector<float> a = this.direction.x * this.direction.x + this.direction.y * this.direction.y;
Vector<float> b = r + s;
Vector<float> c = p * p + q * q - other.radius * other.radius;
Vector<float> DSqrt = b * b - fours * a * c;
Vector<int> maskCollision = Vector.GreaterThan(DSqrt, zeros);
// Commented out to make the benchmark more fair.
//if (Vector.Dot(maskCollision, maskCollision) == 0)
//{ return maskCollision; }
// ------------------------------------------------------- \\
// Update t if and only if there is a collision. Take \\
// note of the conditional where we compute t. \\
Vector<float> D = Vector.SquareRoot(DSqrt);
Vector<float> bMinus = Vector.Negate(b);
Vector<float> twoA = twos * a;
Vector<float> t0 = (bMinus + D) / twoA;
Vector<float> t1 = (bMinus - D) / twoA;
Vector<float> tm = Vector.ConditionalSelect(Vector.LessThan(t1, t0), t1, t0);
Vector<int> maskBiggerThanZero = Vector.GreaterThan(tm, zeros);
Vector<int> maskSmallerThanT = Vector.LessThan(tm, this.t);
Vector<int> mask = Vector.BitwiseAnd(
maskCollision,
Vector.BitwiseAnd(
maskBiggerThanZero,
maskSmallerThanT)
);
this.t = Vector.ConditionalSelect(
mask, // the bit mask that allows us to choose.
tm, // the smallest of the t's.
t); // if the bit mask is false (0), then we get our original t.
return mask;
}
Код сборки
Их можно найти на pastebin. Обратите внимание, что есть некоторая шаблонная сборка из инструмента для тестирования. Вам нужно посмотреть на вызовы функций.
- BoundingBox(Pack): https://pastebin.com/RYMQdZMh
- Круг (пакет) изменен: https://pastebin.com/YZHjc1vY
- Circle(Pack) Оригинал: https://pastebin.com/87nYgZrv
Сравнительный анализ
Я сравнивал ситуацию с BenchmarkDotNet.
Результаты для Circle / CirclePack (обновлено):
| Method | Mean | Error | StdDev |
|------------------- |---------:|----------:|----------:|
| Intersection | 9.710 ms | 0.0540 ms | 0.0505 ms |
| IntersectionPacked | 3.296 ms | 0.0055 ms | 0.0051 ms |
Результаты для BoundingBox/BoundingBoxPacked:
| Method | Mean | Error | StdDev |
|------------------- |----------:|----------:|----------:|
| Intersection | 24.269 ms | 0.2663 ms | 0.2491 ms |
| IntersectionPacked | 1.152 ms | 0.0229 ms | 0.0264 ms |x
Благодаря AVX ожидается ускорение примерно в 6-8 раз. Ускорение ограничивающей рамки значительно, тогда как ускорение круга довольно низкое.
Возвращаясь к вопросу вверху: почему одно ускорение выше крыши, а другое довольно низкое? И как нижний из двух (CirclePack) может стать быстрее?
Правки в отношении Питера Кордеса (комментарии)
- Сделали тест более справедливым: версия с одним лучом не имеет раннего ответвления, как только луч больше не может сталкиваться. Теперь ускорение примерно в 2,5 раза.
- Добавил ассемблерный код отдельным заголовком.
- Что касается квадратного корня: это влияет, но не так сильно, как кажется. Удаление квадратного корня из вектора сокращает общее время примерно на 0,3 мс. Код с одним лучом теперь также всегда выполняет извлечение квадратного корня.
- Вопрос о FMA (сложение с плавным умножением) в C #. Я думаю, что это подходит для скаляров (см. умножить-сложить?), но я не нашел подобной операции в структуре System.Numerics.Vector.
- Об инструкции C# для упакованного min: Да, это так. Я такой глупый. Я даже использовал его уже.
sqrt
?vsqrtps ymm
имеет гораздо худшую пропускную способность, чем FMA, в зависимости от аппаратного обеспечения. В худшем случае это будет Intel Haswell (2 256-битных FMA на такт, но без улучшений Broadwell и Skylake для делителя). Есть ли в C# опция быстрой математики, позволяющая сокращать операции mul + add в исходном коде в FMA в как м? И/или позволить использовать приближениеrsqrtps
? (Я знаю C++ и x86 asm/производительность.) - person Peter Cordes   schedule 09.07.2019vsqrtps ymm
пропускной способности 1 на 14 циклов; См. Деление с плавающей запятой и умножение с плавающей запятой. Пропускная способность и задержка FP sqrt очень похожи на деление FP на большинстве процессоров и идентичны на Haswell. Ваша SIMD-версия делает все это без ответвлений, вместо того, чтобы пытаться выйти раньше (что будет работать только в том случае, если все 8 векторных элементов имеют там отрицательное значение). - person Peter Cordes   schedule 09.07.2019Vector.ConditionalSelect
(vblendvps ymm
) в любом случае требует AVX, а не SSE4.1. - person Peter Cordes   schedule 09.07.2019Vector.ConditionalSelect( Vector.LessThan(...), ...)
Разве в C# нет встроенной функции для упакованных минут? SSE/AVX имеетminps
в аппаратном обеспечении, но я бы беспокоился, что этот источник заставит компилятор сделать что-то ужасное, напримерcmpps
->blendvps
- person Peter Cordes   schedule 09.07.2019if()
может избежать этого, в зависимости от того, насколько хорошо JIT справляется с использованием таких инструкций, какcmpss
вместо FP, сравнивает в целочисленные ФЛАГИ сucomiss
- person Peter Cordes   schedule 09.07.2019