我实际上在优化算法时遇到了一些麻烦:
我有一个磁盘(以0为中心，半径为1)，上面充满了三角形(不一定具有相同的面积/长度)。可能会有大量的三角形(比如说从1k到 300k 三角形)
我的目标是尽快找到一个点所属的三角形。
必须重复操作多次(大约是 10k乘以)。
现在，我使用的算法是:我正在计算每个三角形中点的重心坐标。如果第一个系数在0到1之间，我继续。如果不是，我停下来。然后，我用相同的想法计算第二个系数，第三个计算，然后对每个三角形进行计算。
我想不出一种方法来利用我在光盘上工作的事实(以及我有一个欧几里得距离来帮助我直接“瞄准”好三角形的事实):如果我尝试计算我指向三角形的每个“中心”:
1)它比我用重心坐标蛮力地进行的操作要多
2)我将必须订购一个 vector ，该 vector 包含所有三角形到我的点的欧几里得距离。
3)我绝对不能保证最接近我的三角形就是好三角形。
我觉得自己缺少了一些东西，可以预先计算，这可以帮助我在开始使用蛮力部分之前发现好“区域”。
该算法已经并行化(使用OpenMP):我正在并行调用以下函数。
bool Triangle2D::is_in_triangle(Vector2d Point) {
double denominator = ((Tri2D(1, 1) - Tri2D(2, 1))*(Tri2D(0, 0) - Tri2D(2, 0)) + (Tri2D(2, 0) - Tri2D(1, 0))*(Tri2D(0, 1) - Tri2D(2, 1)));
// Computing the first coefficient
double a = ((Tri2D(1, 1) - Tri2D(2, 1))*(Point(0) - Tri2D(2, 0)) + (Tri2D(2, 0) - Tri2D(1, 0))*(Point(1) - Tri2D(2, 1))) / denominator;
if (a < 0 || a>1) {
return(false);
}
// Computing the second coefficient
double b = ((Tri2D(2, 1) - Tri2D(0, 1))*(Point(0) - Tri2D(2, 0)) + (Tri2D(0, 0) - Tri2D(2, 0))*(Point(1) - Tri2D(2, 1))) / denominator;
if (b < 0 || b>1) {
return(false);
}
// Computing the third coefficient
double c = 1 - a - b;
if (c < 0 || c>1) {
return(false);
}
return(true);
}
下一步可能是研究GPU并行化，但是我需要确保代码背后的想法足够好。
目前，对于
2min30三角形和75k点，大约需要 10k ，但这还不够快。编辑:
Triangle2D使用特征矩阵存储坐标
最佳答案
所有长 mustache 的 HPC专业人士都允许在这里进行一些精心设计的方法，(对于我来说，这可能会使我们的社区成员感到有趣，即使不被我们的社区成员所享用，他们也会觉得自己比您年轻一些)专业的自我感觉，可能会对对性能驱动的代码设计，性能调整以及其他并行代码的风险和 yield 有更深入了解的人感兴趣，而您对硬核HPC计算经验的了解是如此之深。谢谢。
a)算法(原样)可以使〜2倍的速度提速-令人垂涎的水果+更多惊喜
b)其他算法可能会获得〜40〜80X的加速提升
c)最佳并行代码+最佳性能的提示
GOAL :在具有i5 8500、3GHz，6core，NVIDIA Quadro P400的计算机上，即使在 10k三角形300k三角形中的的目标运行时间为 2-3分钟，即使没有尝试这很值得)
尽管这似乎是一段漫长的旅程，但问题仍然存在，值得我们仔细观察，因此，请在表现为最高性能的思考过程中耐心等待。
a)算法(原样)分析:
重心坐标系的原样使用是一个不错的技巧，在最佳情况下，直接实现该方法的成本大约是(20 FLOP + 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-->
b)其他方法到算法:
让我们回顾一下另一种方法，它由于数据量少和指令少而显得既快又便宜，并且还有望具有高密度，一旦巧妙地使用矢量化AVX-2或更好的AVX-512 vector 指令，每个核心:ANIMATED, fully-INTERACTIVE explanation连同分析问题重新制定here。
对点到线距离(每条线为
ax + by + c = 0)进行三重测试只需2 ozt_code的便宜成本就可以了+进行符号测试(如果vector-4-in-1-compact AVX- 2/8合1 AVX-512 FMA3 -s)尽管可以“快速”决定是否要针对相应的三角形测试一个点，但可以通过静态的，预先计算的
VFMADD元组和“”将 polar(R,Theta) 坐标空间中的每个三角形静态地预先“定帧”，快速”区分每个点(如果它不适合这种极段)。但是，执行此操作的成本(数据访问模式成本+这些“快速”指令的成本)将超出不需要的任何可能的“已保存”指令路径(如果在极坐标段之外找到了一个点) )。在3.0+ GHz的频率下，每6个CPU内核以大约9条CPU指令的成本实现了24个三角点测试的性能，极地段的“预测试”将突然变得昂贵得令人无法接受。大约是二阶负面影响(通过更糟糕的缓存命中率/缓存未命中率引入，因为要存储更多数据并将其读入“快速”预测试〜每个三角形“成帧”极性为+ 16B，分段元组每点+ 8B(对高速缓存命中率/未命中率的影响最大)。显然，这不是进一步行动的好方向，因为绩效将下降而不是提高，这是我们的全局战略。
Intel i5 8500 CPU只能使用AVX-2，因此，如果需要，可以最紧凑地使用每个内核每个CPU时钟滴答8个三角形，以实现甚至两倍的更高性能。
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带来 yield 的优势
如果您的博士导师，教授，老板或项目经理确实坚持要您针对此问题开发GPU计算c++ / CUDA代码解决方案，那么最好的下一步就是要求获得用于此类任务的任何更合适的GPU卡，而不是您发布的那个。
您所显示的Q400 GPU has just 2 SMX (48KB L1-cache each)卡不太适合进行严重的并行CUDA计算，它几乎没有30X的处理设备可以实际执行任何SIMD线程块计算，更不用说它的小内存和SMX上的小L1 / L2了。 -快取。因此，在完成所有与CUDA相关的设计和优化工作之后，SIMD线程中将不会有更多(但只有一对)GPU-SMX
( R_min, R_max, Theta_min, Theta_max )范围内的线程块执行，因此，可以从基于GP107的设备中获得期望，并且有30台X更好的设备可用于某些确实高性能的并行处理可用COTS)步骤0:预先计算和预先安排数据( MAXIMIZE高速缓存行COHERENCY ):
在这里，优化算法的最佳宏范围循环是有意义的，这样您就可以从缓存命中中“受益”最多(即最好地重用已经预先提取的“快速”数据)
因此，可以测试一下，将工作分配给N个并行工作的人是否更快，而N个并行工作的人处理离散的三角形工作块，它们每个都在最小的内存区域上循环-所有点( 〜10k * 2 * 4B〜80 kB)，然后再移入工作块中的下一个三角形。确保行优先阵列对齐到内存区域是至关重要的(FORTRAN的人可以用HPC快速紧凑/矢量化 vector 代数和矩阵处理的技巧告诉很多费用/ yield )
有好处吗？
缓存的系数将被重复使用约1万次(以
warp32 的成本为代价，而不是重新获取~ 0.5~1.0 [ns] 的成本，如果这些必须通过RAM存储器访问来重新读取)。区别在于 + 100 ~ 300 [ns]是值得努力的，以使从属于数据访问模式和缓存行资源的工作流最佳地对齐。 结果是原子的-任何点将属于一个且仅属于一个(不重叠)三角形。
这简化了对结果 vector 的先验非冲突且相对稀疏的写入，在该 vector 中，任何检测到的三角形内点都可以自由报告找到的这种三角形的索引。
使用此结果 vector 来潜在地避免对已经进行过匹配(索引为非负数)的点进行任何重新测试不是很有效，因为重新读取此类指示并重新排列紧凑对齐方式的成本对于不重新测试已经映射的点的潜在“节省”，进行4合1或8合1三角形点测试将变得不利。
因此，
~ 200X ~ 600X指向一块 10k 三角形上的点可能会产生以下工作流:每个
300k核心300k / 6三角形的分割/每个核心1个核心~ 50k三角测试 ~ 50k * 10 k每个核心三角测试~ 500M @ 3.0+ GHz vector 4合1紧凑型测试执行每个核心 ~ 125M AVX-2 测试x ~ 125M-指令10 uops ...这是@ 3.0 GHz ...秒？是的，这里有很大的动力去实现这一最终性能，并以此方式指导进一步的工作。所以，我们在这里:
仅在6核AVX-2 i5 8500 @ 3.0+ GHz上，可实现的HPC目标在几秒钟内即可形成300 + k三角形/ 10 + k点的范围
值得努力，不是吗？
关于c++ - 查找点属于哪个三角形的优化技巧，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/56937577/