Все длиннобородые HPC-профессионалы , пожалуйста, разрешите немного научный подход, который может (по моему честному мнению) стать интересным, если не понравится нашим членам Сообщества, которые чувствуют сами они немного младше, чем вы профессионально себя чувствуете, и которые могут заинтересоваться более глубоким взглядом на проектирование кода, мотивируемое производительностью, настройкой производительности и другими рисками и преимуществами параллельного кода, которые вы знаете на своем собственном жестком ядре HPC- Компьютерный опыт так хорошо и так глубоко. Спасибо.
а) АЛГОРИТМ (как есть) может получить ~ 2-кратное ускорение для низко висящих фруктов + еще больше сюрпризов2в
б) ДРУГОЙ АЛГОРИТМ может получить ~ 40 ~ 80-кратное ускорение из-за геометрии
в) СОВЕТЫ ПО ЛУЧШЕМУ ПАРАЛЛЕЛЬНОМУ КОДУ + ОТЛИЧНАЯ ЭФФЕКТИВНОСТЬ
ЦЕЛЬ : целевое время выполнения для 10k точек в 300k треугольников будет 2-3 мин на компьютере с i5 8500, 3GHz, 6core, NVIDIA Quadro P400 (нужно попробовать вычисления на GPU, даже не уверен, стоит ли это того)
Хотя это может показаться долгим путешествием, проблема приятна и заслуживает более пристального внимания, поэтому, пожалуйста, потерпите меня во время потока мотивированных мыслей с максимальной эффективностью.
а) АЛГОРИТМ (как есть) АНАЛИЗ:
Барицентрическая система координат «как есть» является хорошим трюком, прямая реализация которого стоит чуть больше, чем примерно (20 FLOPs + 16 MEM / REG-I / O-ops) в лучшем случае и чуть выше ( 30 FLOPs + 30 MEM / REG-I / O-ops).
Есть несколько штрихов, которые могут снизить эти затраты на исполнение, избегая дорогостоящих и даже несущественных операций:
--------------------------------------------------------------------------------------
double denominator = ( ( Tri2D( 1, 1 )
- Tri2D( 2, 1 ) // -------------------------- 2x MEM + OP-1.SUB
) * ( Tri2D( 0, 0 ) //--------------------- + OP-3.MUL
- Tri2D( 2, 0 ) //--------------------- 2x MEM + OP-2.SUB
) + ( Tri2D( 2, 0 ) //--------------- + OP-7.ADD
- Tri2D( 1, 0 ) //--------------- 2x MEM + OP-4.SUB
) * ( Tri2D( 0, 1 ) //--------- + OP-6.MUL
- Tri2D( 2, 1 ) //--------- 2x MEM + OP-5.SUB
)
);
// Computing the first coefficient ------------------------------------------------------------------------------------------------------
double a = ( ( Tri2D( 1, 1 )
- Tri2D( 2, 1 ) //-------------------------- 2x MEM + OP-8.SUB
) * ( Point(0) //------------------------ + OP-A.MUL
- Tri2D( 2, 0 ) //--------------------- 2x MEM + OP-9.SUB
) + ( Tri2D( 2, 0 ) //--------------- + OP-E.ADD
- Tri2D( 1, 0 ) //--------------- 2x MEM + OP-B.SUB
) * ( Point(1) //-------------- + OP-D.MUL
- Tri2D( 2, 1 ) //--------- 2x MEM + OP-C.MUL
)
) / denominator; //-------------------------- 1x REG + OP-F.DIV //----------- MAY DEFER THE MOST EXPENSIVE DIVISION UNTIL a third coeff is indeed first needed, if ever------------[3]
if (a < 0 || a>1) { // ----------------------------------------------------------------------------- a < 0 ~~ ( sign( a ) * sign( denominator ) ) < 0
return(false); // ------------------------------------------------------------------------------ a > 1 ~~ || a > denominator
}
// Computing the second coefficient
double b = ( ( Tri2D( 2, 1 ) - Tri2D( 0, 1 ) ) //--------- 2x MEM + OP-16.SUB
* ( Point(0) - Tri2D( 2, 0 ) ) //--------- 2x MEM + OP-17.SUB + OP-18.MUL
+ ( Tri2D( 0, 0 ) - Tri2D( 2, 0 ) ) //--------- 2x MEM + OP-19.SUB + OP-22.ADD
* ( Point(1) - Tri2D( 2, 1 ) ) //--------- 2x MEM + OP-20.SUB + OP-21.MUL
) / denominator; //-------------------------- 1x REG + OP-23.DIV //---------- MAY DEFER THE MOST EXPENSIVE DIVISION UNTIL a third coeff is indeed first needed, if ever -----------[3]
if (b < 0 || b>1) { // ----------------------------------------------------------------------------- b < 0 ~~ ( sign( b ) * sign( denominator ) ) < 0
return(false); // ------------------------------------------------------------------------------ b > 1 ~~ || b > denominator
}
// Computing the third coefficient
double c = 1 - a - b; // ------------------------------------------- 2x REG + OP-24.SUB + OP-25.SUB
// 1 -(a - b)/denominator; //--------------------------------------------------------------- MAY DEFER THE MOST EXPENSIVE DIVISION EXECUTED BUT HERE, IFF INDEED FIRST NEEDED <---HERE <----------[3]
повторные переоценки, которые появляются в оригинале, могут быть явно обработаны путем ручного присвоения / повторного использования, однако, есть шанс, что хороший оптимизирующий компилятор может получить их при использовании -O3 флаг принудительной оптимизации.
не стесняйтесь профилировать даже этот низко висящий фрукт для полировки самых дорогих деталей.
//------------------------------------------------------------------
double Tri2D_11_sub_21 = ( Tri2D( 1, 1 )
- Tri2D( 2, 1 )
), //====================================================== 2x MEM + OP-a.SUB (REG re-used 2x)
Tri2D_20_sub_10 = ( Tri2D( 2, 0 )
- Tri2D( 1, 0 )
), //====================================================== 2x MEM + OP-b.SUB (REG re-used 2x)
Tri2D_00_sub_20 = ( Tri2D( 0, 0 )
- Tri2D( 2, 0 )
); //====================================================== 2x MEM + OP-c.SUB (REG re-used 1~2x)
//-----------------------
double denominator = ( ( /*
Tri2D( 1, 1 )
- Tri2D( 2, 1 ) // -------------------------- 2x MEM + OP-1.SUB (avoided by re-use) */
Tri2D_11_sub_21 //=========================================== 1x REG + OP-d.MUL
) * ( /*
Tri2D( 0, 0 ) //--------------------- + OP-3.MUL
- Tri2D( 2, 0 ) //--------------------- 2x MEM + OP-2.SUB (avoided by re-use) */
Tri2D_00_sub_20 //===================================== 1x REG + OP-f.ADD
) + ( /*
Tri2D( 2, 0 ) //--------------- + OP-7.ADD
- Tri2D( 1, 0 ) //--------------- 2x MEM + OP-4.SUB (avoided by re-use) */
Tri2D_20_sub_10 //=============================== 1x REG + OP-e.MUL
) * ( Tri2D( 0, 1 ) //--------- + OP-6.MUL
- Tri2D( 2, 1 ) //--------- 2x MEM + OP-5.SUB
)
);
//\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/
//
// Computing the first coefficient ------------------------------------------------------------------------------------------------------
//
double enumer_of_a = ( ( /*
Tri2D( 1, 1 )
- Tri2D( 2, 1 ) //-------------------------- 2x MEM + OP-8.SUB (avoided by re-use) */
Tri2D_11_sub_21 //=========================================== 1x REG + OP-g.MUL
) * ( Point(0) //------------------------------------------ + OP-i.MUL
- Tri2D( 2, 0 ) //--------------------------------------- 2x MEM + OP-h.SUB
) + ( /*
Tri2D( 2, 0 ) //--------------- + OP-E.ADD
- Tri2D( 1, 0 ) //--------------- 2x MEM + OP-B.SUB (avoided by re-use) */
Tri2D_20_sub_10 //=============================== 1x REG + OP-l.ADD
) * ( Point(1) //-------------------------------- + OP-k.MUL
- Tri2D( 2, 1 ) //--------------------------- 2x MEM + OP-j.MUL
)
);/*denominator; *///------------------------ 1x REG + OP-F.DIV (avoided by DEFERRAL THE MOST EXPENSIVE DIVISION UNTIL a third coeff is indeed first needed, if ever-----------[3]
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~J.I.T./non-MISRA-C-RET-->
// TEST CONDITIONS FOR A CHEAPEST EVER J.I.T./non-MISRA-C-RET-->
//
if ( enumer_of_a > denominator // in a > 1, THE SIZE DECIDES, the a / denominator > 1, iff enumer_of_a > denominator a rather expensive .FDIV is avoided at all
|| enumer_of_a * denominator < 0 ) return(false); // in a < 0, THE SIGN DECIDES, not the VALUE matters, so will use a cheaper .FMUL, instead of a rather expensive .FDIV ~~ ( sign( a ) * sign( denominator ) ) < 0
//\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/
//
// Computing the second coefficient
//
double enumer_of_b = ( ( Tri2D( 2, 1 ) - Tri2D( 0, 1 ) ) //---------------------------------------- 2x MEM + OP-m.SUB
* ( Point(0) - Tri2D( 2, 0 ) ) //---------------------------------------- 2x MEM + OP-n.SUB + OP-o.MUL
+ ( /*
Tri2D( 0, 0 ) - Tri2D( 2, 0 ) //--------- 2x MEM + OP-19.SUB + OP-22.ADD (avoided by re-use) */
Tri2D_00_sub_20 //======================================================== 1x REG + OP-p.ADD
)
* ( Point(1) - Tri2D( 2, 1 ) ) //---------------------------------------- 2x MEM + OP-r.SUB + OP-q.MUL
);/*denominator; *///------------------------ 1x REG + OP-23.DIV (avoided by DEFERRAL THE MOST EXPENSIVE DIVISION UNTIL a third coeff is indeed first needed, if ever-----------[3]
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~J.I.T./non-MISRA-C-RET-->
// TEST CONDITIONS FOR A 2nd CHEAPEST J.I.T./non-MISRA-C-RET-->
//
if ( enumer_of_b > denominator // in b > 1, THE SIZE DECIDES, the a / denominator > 1, iff enumer_of_a > denominator a rather expensive .FDIV is avoided at all
|| enumer_of_b * denominator < 0 ) return(false); // in b < 0, THE SIGN DECIDES, not the VALUE matters, so will use a cheaper .FMUL, instead of a rather expensive .FDIV ~~ ( sign( a ) * sign( denominator ) ) < 0
//\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/
//
// Computing the third coefficient
//
double c = 1 - ( ( enumer_of_a
- enumer_of_b
)
/ denominator
); // --------------------------------------------- 3x REG + OP-s.SUB + OP-t.FDIC + OP-u.SUB <----THE MOST EXPENSIVE .FDIV BUT HERE, IFF INDEED FIRST NEEDED <---HERE <------------[3]
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~J.I.T./non-MISRA-C-RET-->
// TEST CONDITIONS FOR PRE-FINAL RET J.I.T./non-MISRA-C-RET-->
//
if ( c < 0 || c > 1 ) return( false );
return( true ); //~~~~~~~~~~~~~~~~~~~ "the-last-resort" RET-->
б) ДРУГОЙ ПОДХОД К АЛГОРИТМУ:
Давайте рассмотрим другой подход, который кажется и более быстрым, и более дешевым, из-за меньшего количества данных и наименее инструкций, а также обещает иметь высокую плотность, когда разумное использование векторизованных векторных инструкций AVX-2 или более эффективных AVX-512 приведет к получить упряжку для каждого ядра: АНИМИРОВАННОЕ, полностью ИНТЕРАКТИВНОЕ объяснение в совокупности с переформулировкой аналитической проблемы здесь .
Тройной тест с расстоянием точка-линия (каждая строка ax + by + c = 0) обходится по низкой цене: достаточно ~ 2 FMA3 + тест знака (и даже лучше, если вектор-4-в-1-компактный AVX-2/8-в-1 AVX-512 VFMADD -s)
Хотя возможно «быстрое» решение о том, имеет ли смысл проверять точку на соответствующем треугольнике, возможно со статическим предварительным «кадрированием» каждого треугольника в polar(R,Theta) координатном пространстве с помощьюстатический, предварительно вычисленный кортеж ( R_min, R_max, Theta_min, Theta_max ) и «быстрый» различают каждую точку, если она не помещается внутри такого полярного сегмента.Тем не менее, затраты на это (затраты на доступ к данным + стоимость этих «быстрых» инструкций) превысят любые потенциально «сохраненные» пути инструкций, которые не должны иметь место (если точка находится за пределами полярного сегмента).).Достигнув диапазона производительности в 24 теста «точки в треугольнике» при стоимости ~ 9 инструкций ЦП на 6 ядер ЦП при 3,0+ ГГц, «предварительное тестирование» в полярном сегменте внезапно станет чрезмерно дорогим, если не говоритьо негативном эффекте второго порядка (вызванном гораздо худшим соотношением попаданий в кэш / пропаданий в кэше, учитывая, что больше данных нужно хранить и считывать в «быстрый» предварительный тест ~ + 16B на треугольник, «обрамляющий» полярныйсегмент сегмента + 8B на точку (с наихудшим влиянием на соотношение попаданий и промахов в кэше).
Это явно не хорошее направление для дальнейшего движения, так как производительность будет снижаться, а не увеличиваться, что является нашейглобальная стратегия здесь.
Процессор Intel i5 8500 можно использовать, но не AVX-2, поэтому наиболее компактное использование 8-треугольников на такт такта процессора на ядрооставлено для достижения даже в 2 раза более высокой производительности, если необходимо.
TRIPLE-"point-above-line"-TEST per POINT per TRIANGLE:
---------------------------------------------------------------------------------
PRE-COMPUTE STATIC per TRIANGLE CONSTANTS:
LINE_1: C0__L1, C1__L1, C2__L1, bool_L1DistanceMustBePOSITIVE
LINE_2: C0__L2, C1__L2, C2__L2, bool_L2DistanceMustBePOSITIVE
LINE_3: C0__L3, C1__L3, C2__L3, bool_L3DistanceMustBePOSITIVE
TEST per TRIANGLE per POINT (Px,Py) - best executed in an AVX-vectorised fashion
LINE_1_______________________________________________________________
C0__L1 ( == pre-calc'd CONST = c1 / sqrt( a1^2 + b1^2 ) ) //
Px * C1__L1 ( == pre-calc'd CONST = a1 / sqrt( a1^2 + b1^2 ) ) // OP-1.FMA REG-Px,C1__L1,C0__L1
Py * C2__L1 ( == pre-calc'd CONST = b1 / sqrt( a1^2 + b1^2 ) ) // OP-2.FMA REG-Py,C2__L1, +
.GT./.LT. 0 // OP-3.SIG
LINE_2_______________________________________________________________
C0__L2 ( == pre-calc'd CONST = c2 / sqrt( a2^2 + b2^2 ) ) //
Px * C1__L2 ( == pre-calc'd CONST = a2 / sqrt( a2^2 + b2^2 ) ) // OP-4.FMA REG-Px,C1__L2,C0__L2
Py * C2__L2 ( == pre-calc'd CONST = b2 / sqrt( a2^2 + b2^2 ) ) // OP-5.FMA REG-Py,C2__L2, +
.GT./.LT. 0 // OP-6.SIG
LINE_3_______________________________________________________________
C0__L3 ( == pre-calc'd CONST = c3 / sqrt( a3^2 + b3^2 ) ) //
Px * C1__L3 ( == pre-calc'd CONST = a3 / sqrt( a3^2 + b3^2 ) ) // OP-7.FMA REG-Px,C1__L3,C0__L3
Py * C2__L3 ( == pre-calc'd CONST = b3 / sqrt( a3^2 + b3^2 ) ) // OP-8.FMA REG-Py,C2__L3, +
.GT./.LT. 0 // OP-9.SIG
( using AVX-2 intrinsics or inlined assembler will deliver highest performance due to COMPACT 4-in-1 VFMADDs )
____________________________________________
| __________________________________________triangle A: C1__L1
| | ________________________________________triangle B: C1__L1
| | | ______________________________________triangle C: C1__L1
| | | | ____________________________________triandle D: C1__L1
| | | | |
| | | | | ______________________________
| | | | | | ____________________________triangle A: Px
| | | | | | | __________________________triangle B: Px
| | | | | | | | ________________________triangle C: Px
| | | | | | | | | ______________________triandle D: Px
| | | | | | | | | |
|1|2|3|4| | | | | |
| | | | | |1|2|3|4| ________________
| | | | | | | | | | | ______________triangle A: C0__L1
| | | | | | | | | | | | ____________triangle B: C0__L1
| | | | | | | | | | | | | __________triangle C: C0__L1
| | | | | | | | | | | | | | ________triandle D: C0__L1
| | | | | | | | | | | | | | |
|1|2|3|4| | | | | | | | | | |
| | | | | |1|2|3|4| | | | | |
| | | | | | | | | | |1|2|3|4|
(__m256d) __builtin_ia32_vfmaddpd256 ( (__v4df )__A, (__v4df )__B, (__v4df )__C ) ~/ per CPU-core @ 3.0 GHz ( for actual uops durations check Agner or Intel CPU documentation )
can
perform 4-( point-in-triangle ) PARALLEL-test in just about ~ 9 ASSEMBLY INSTRUCTIONS / per CPU-core @ 3.0 GHz
24-( point-in-triangle ) PARALLEL-test in just about ~ 9 ASSEMBLY INSTRUCTIONS / per CPU
using AVX-512 empowered CPU, can use 8-in-1 VFMADDs
could
perform 8-( point-in-triangle ) PARALLEL-test in just about ~ 9 ASSEMBLY INSTRUCTIONS / per CPU-core @ 3.0 GHz
48-( point-in-triangle ) PARALLEL-test in just about ~ 9 ASSEMBLY INSTRUCTIONS / per CPU
c) СОВЕТЫ ДЛЯ ЛУЧШЕГО ПАРАЛЛЕЛЬНОГО КОДА + ОТЛИЧНАЯ ЭФФЕКТИВНОСТЬ:
Шаг -1: GPU / CUDA стоит v / sПреимущества
Если ваш PhD-наставник, профессор, начальник или руководитель проекта действительно настаивает на том, чтобы вы разработали решение для кода на языке C ++ / CUDA на GPU-вычислениях для этой самой проблемы, лучшим следующим шагом будет поиск лучшего решения.GPU-карта для такой задачи, чем та, которую вы выложили.
На указанной вами карте, Q400 GPU есть только 2 SMX (48 КБ L1-кеш каждый) не совсем подходит для серьезно подразумеваемых параллельных CUDA-вычислений, имея около 30-кратных вычислительных устройств для выполнения любых SIMD-потоковых вычислений, не говоря уже о его малой памяти и крошечных L1 / L2-кешах на SMX.Таким образом, после всех усилий по проектированию и оптимизации, связанных с CUDA, в SIMD-потоках будет не больше, а одна (!) Пара потоковых блоков GPU-SMX warp32 в масштабе всего потока, поэтому большихЦирк следует ожидать от этого устройства на базе GP107, и для некоторых действительно высокопроизводительных параллельных обработок COTS имеется более 30 устройств с улучшенной X-скоростью
Шаг 0: предварительное вычисление иПредварительно упорядочить данные ( MAXIMIZE строка кэша COHERENCY ):
Здесь имеет смысл оптимизировать наилучший макроскопический цикл алгоритма, чтобы вы "извлекали выгоду" больше всего из попаданий в кэш (т. е. лучшее повторное использование «быстрых» данных, уже предварительно извлеченных)
Итак, может проверить, является ли кэш-память быстрее распределять работу среди N-одновременных рабочих, которые обрабатывают разрозненную работу-блоки треугольников, где каждый из них проходит по наименьшей области памяти --- все точки (~ 10k * 2 * 4B ~ 80 кБ), прежде чем перейти к следующему треугольнику в рабочем блоке.Очень важно обеспечить выравнивание массива по строке в области памяти (ребята из FORTRAN могут многое рассказать о затратах и преимуществах только этого трюка для высокопроизводительной компактной / векторизованной векторной алгебры и матричной обработки HPC)
Преимущество?
Кэшированные коэффициенты будут повторно использоваться ~ 10 000 раз (по стоимости ~ 0.5~1.0 [ns] вместо затрат на повторное извлечение + 100 ~ 300 [ns], если они должны быть перечитаны при доступе к памяти RAM).Разница в том, что около ~ 200X ~ 600X стоит усилий, чтобы лучше выровнять рабочий процесс, подчиненный шаблонам доступа к данным и ресурсам строки кэша.
Результаты являются атомарными - любая точка будет принадлежатьк одному и только одному (непересекающемуся) треугольнику.
Это упрощает априорные неколлизионные и относительно редкие записи в результирующий вектор, в который любая обнаруженная точка внутри треугольника может свободно сообщать индекс найденного такого треугольника.
Использование этого вектора результатовпотенциально избегать повторного тестирования в точках, где уже было выполнено совпадение (индекс не отрицательный), не очень эффективно, поскольку затраты на повторное считывание такой индикации и реорганизацию компактного выравнивания из 4 вТесты «точка-в-треугольнике» 1 или 8 в 1 могут стать слишком дорогими для потенциальной «экономии», если не проводить повторное тестирование уже нанесенной на карту точки.
Таким образом, 10k указывает на блок 300k треугольники могут дать рабочий процесс:
разделение некоторых 300k / 6 ядер
~ 50k треугольники / 1 ядро
~ 50k * 10 k тестирование точек на треугольник на ядро
~ 500M тесты «точка-треугольник» на ядро @ 3.0+ GHz
~ 125M AVX-2 компактные тесты «вектор-4-в-1» на ядро
~ 125M тесты x 10 uops -инструкции @ 3.0 GHz ...
что такое 1.25G uops @ 3.0+ GHz ... секунд?Да, здесь есть довольно сильная мотивация, чтобы идти к этой конечной производительности и направлять дальнейшую работу таким образом.
ТАК, НАМ МЫ:
Принципиально достижимая цель HPC - этов диапазоне нескольких секунд для 300 + k треугольников / 10 + k очков на 6-ядерном AVX-2 i5 8500 @ 3.0+ ГГц
Стоит усилий, не так ли?