手段（或其他优化）
注意：我希望更多地了解如何处理和提出这些解决方案，而不是解决方案本身。
在我的系统中，我有一个非常关键的性能功能，在特定环境下显示为头号性能分析热点。 它处于k-means迭代的中间（已经使用并行处理多线程来处理每个工作线程中的点的子范围）。
ClusterPoint& pt = points[j];
pt.min_index = -1;
pt.min_dist = numeric_limits<float>::max();
for (int i=0; i < num_centroids; ++i)
{
const ClusterCentroid& cent = centroids[i];
const float dist = ...;
if (dist < pt.min_dist) // <-- #1 hotspot
{
pt.min_dist = dist;
pt.min_index = i;
}
}
处理这部分代码所需的时间会大大减少，所以我经常一直在摆弄它。 例如，将质心循环置于外部可能是值得的，并且对于给定的质心并行地遍历点。 这里的聚类点的数量跨越数百万，而质心的数量跨越数千个。 该算法适用于少数迭代（通常在10以下）。 它不寻求完美的收敛/稳定性，只是一些“合理”的近似。
任何想法都是值得赞赏的，但是我真正渴望发现的是，如果这个代码可以做到无分支，因为它可以允许SIMD版本。 我并没有真正发展出那种能够轻易掌握如何提出无分支解决方案的思维能力：我的大脑在那里失败的时候很像我在早期第一次接触递归时所做的那样，所以有关如何编写无分支代码以及如何为它开发适当的思维方式也是有帮助的。
总之，我正在寻找关于如何微代码优化的指南和提示和建议（不一定是解决方案）。 它很可能有改进算法的空间，但是我的盲点一直在微型优化解决方案中（我很想知道如何更有效地应用它们而不会过度使用它）。 对于逻辑来说，它已经是紧凑的多线程并行处理器，所以我几乎被推入微型优化角落，作为一种更快速的尝试，而不需要更聪明的算法。 我们完全可以自由改变内存布局。
对算法建议的回应
关于在寻求在算法级别上明显改进的O（knm）算法的微观优化方面寻找这一切都是错误的，我完全同意。 这将这个具体问题推向了一个有点学术和不切实际的领域。 但是，如果我可以获得一个轶事，我来自高级程序设计的原始背景 - 重点强调广泛，大规模的观点，安全性，而对低级实施细节则很少。 我最近将项目切换到了一种非常不同的现代风格的项目，我正在学习来自我的同行的缓存效率，GPGPU，无分支技术，SIMD，专用mem分配器的各种新技巧，这些技巧实际上超过了malloc但针对特定场景）等
这正是我试图赶上最新性能趋势的地方，而且令人惊讶的是，我发现那些90年代我经常喜欢的那些通常是链接/树型结构的旧数据结构实际上正在大大超越更为天真，粗略的，微型优化的并行代码，在连续内存块上应用调优指令。 这让我感觉有些失望，因为我觉得现在我们正在为机器增加更多的算法，并通过这种方式缩小了可能性（特别是对于GPGPU）。
最有趣的是，我发现这种类型的微型优化快速阵列处理代码比我以前使用的复杂算法和数据结构更容易维护。 首先，他们更容易推广。 此外，我的同事经常会针对某个地区的特定减速情况向客户提出投诉，只是为了一个平行的SIMD和可能的某个SIMD，并称其加速完成。 算法的改进往往可以提供更多，但是这些微优化可以应用的速度和非侵入性让我想在这方面学到更多，因为阅读更好算法的论文可能需要一些时间（并且需要更多广泛的变化）。 因此，我最近一直在微调优化潮流，在这种情况下可能有点过分，但是我的好奇心更多的是扩大我在任何情况下可能的解决方案范围。
拆卸
注意：我在汇编时真的非常糟糕，所以我经常以一种反复的方式调整事情的方式，提出一些有教育意义的猜测，说明为什么vtune中显示的热点可能是瓶颈，然后尝试看看如果时间有所改善，假设如果时间确实有所改进，猜测会有一些真相，或者如果时间没有改变，则完全错过了标记。
000007FEEE3FB8A1 jl thread_partition+70h (7FEEE3FB780h)
{
ClusterPoint& pt = points[j];
pt.min_index = -1;
pt.min_dist = numeric_limits<float>::max();
for (int i = 0; i < num_centroids; ++i)
000007FEEE3FB8A7 cmp ecx,r10d
000007FEEE3FB8AA jge thread_partition+1F4h (7FEEE3FB904h)
000007FEEE3FB8AC lea rax,[rbx+rbx*2]
000007FEEE3FB8B0 add rax,rax
000007FEEE3FB8B3 lea r8,[rbp+rax*8+8]
{
const ClusterCentroid& cent = centroids[i];
const float x = pt.pos[0] - cent.pos[0];
const float y = pt.pos[1] - cent.pos[1];
000007FEEE3FB8B8 movss xmm0,dword ptr [rdx]
const float z = pt.pos[2] - cent.pos[2];
000007FEEE3FB8BC movss xmm2,dword ptr [rdx+4]
000007FEEE3FB8C1 movss xmm1,dword ptr [rdx-4]
000007FEEE3FB8C6 subss xmm2,dword ptr [r8]
000007FEEE3FB8CB subss xmm0,dword ptr [r8-4]
000007FEEE3FB8D1 subss xmm1,dword ptr [r8-8]
const float dist = x*x + y*y + z*z;
000007FEEE3FB8D7 mulss xmm2,xmm2
000007FEEE3FB8DB mulss xmm0,xmm0
000007FEEE3FB8DF mulss xmm1,xmm1
000007FEEE3FB8E3 addss xmm2,xmm0
000007FEEE3FB8E7 addss xmm2,xmm1
if (dist < pt.min_dist)
// VTUNE HOTSPOT
000007FEEE3FB8EB comiss xmm2,dword ptr [rdx-8]
000007FEEE3FB8EF jae thread_partition+1E9h (7FEEE3FB8F9h)
{
pt.min_dist = dist;
000007FEEE3FB8F1 movss dword ptr [rdx-8],xmm2
pt.min_index = i;
000007FEEE3FB8F6 mov dword ptr [rdx-10h],ecx
000007FEEE3FB8F9 inc ecx
000007FEEE3FB8FB add r8,30h
000007FEEE3FB8FF cmp ecx,r10d
000007FEEE3FB902 jl thread_partition+1A8h (7FEEE3FB8B8h)
for (int j = *irange.first; j < *irange.last; ++j)
000007FEEE3FB904 inc edi
000007FEEE3FB906 add rdx,20h
000007FEEE3FB90A cmp edi,dword ptr [rsi+4]
000007FEEE3FB90D jl thread_partition+31h (7FEEE3FB741h)
000007FEEE3FB913 mov rbx,qword ptr [irange]
}
}
}
}
我们被迫瞄准SSE 2--在我们这个时代有点落后，但是当我们假设即使SSE 4作为最低要求（用户有一些英特尔原型机）时，用户群实际上也会绊倒一次。
更新与独立测试：约5.6秒
我非常感谢所有提供的帮助！ 由于代码库相当广泛，触发代码的条件非常复杂（系统事件在多个线程间触发），因此每次进行实验性更改和配置文件都有点难以处理。 所以我已经做了一个表面测试，作为一个独立的应用程序，其他人也可以运行和尝试，以便我可以尝试所有这些优雅的解决方案。
#define _SECURE_SCL 0
#include <iostream>
#include <fstream>
#include <vector>
#include <limits>
#include <ctime>
#if defined(_MSC_VER)
#define ALIGN16 __declspec(align(16))
#else
#include <malloc.h>
#define ALIGN16 __attribute__((aligned(16)))
#endif
using namespace std;
// Aligned memory allocation (for SIMD).
static void* malloc16(size_t amount)
{
#ifdef _MSC_VER
return _aligned_malloc(amount, 16);
#else
void* mem = 0;
posix_memalign(&mem, 16, amount);
return mem;
#endif
}
template <class T>
static T* malloc16_t(size_t num_elements)
{
return static_cast<T*>(malloc16(num_elements * sizeof(T)));
}
// Aligned free.
static void free16(void* mem)
{
#ifdef _MSC_VER
return _aligned_free(mem);
#else
free(mem);
#endif
}
// Test parameters.
enum {num_centroids = 512};
enum {num_points = num_centroids * 2000};
enum {num_iterations = 5};
static const float range = 10.0f;
class Points
{
public:
Points(): data(malloc16_t<Point>(num_points))
{
for (int p=0; p < num_points; ++p)
{
const float xyz[3] =
{
range * static_cast<float>(rand()) / RAND_MAX,
range * static_cast<float>(rand()) / RAND_MAX,
range * static_cast<float>(rand()) / RAND_MAX
};
init(p, xyz);
}
}
~Points()
{
free16(data);
}
void init(int n, const float* xyz)
{
data[n].centroid = -1;
data[n].xyz[0] = xyz[0];
data[n].xyz[1] = xyz[1];
data[n].xyz[2] = xyz[2];
}
void associate(int n, int new_centroid)
{
data[n].centroid = new_centroid;
}
int centroid(int n) const
{
return data[n].centroid;
}
float* operator[](int n)
{
return data[n].xyz;
}
private:
Points(const Points&);
Points& operator=(const Points&);
struct Point
{
int centroid;
float xyz[3];
};
Point* data;
};
class Centroids
{
public:
Centroids(Points& points): data(malloc16_t<Centroid>(num_centroids))
{
// Naive initial selection algorithm, but outside the
// current area of interest.
for (int c=0; c < num_centroids; ++c)
init(c, points[c]);
}
~Centroids()
{
free16(data);
}
void init(int n, const float* xyz)
{
data[n].count = 0;
data[n].xyz[0] = xyz[0];
data[n].xyz[1] = xyz[1];
data[n].xyz[2] = xyz[2];
}
void reset(int n)
{
data[n].count = 0;
data[n].xyz[0] = 0.0f;
data[n].xyz[1] = 0.0f;
data[n].xyz[2] = 0.0f;
}
void sum(int n, const float* pt_xyz)
{
data[n].xyz[0] += pt_xyz[0];
data[n].xyz[1] += pt_xyz[1];
data[n].xyz[2] += pt_xyz[2];
++data[n].count;
}
void average(int n)
{
if (data[n].count > 0)
{
const float inv_count = 1.0f / data[n].count;
data[n].xyz[0] *= inv_count;
data[n].xyz[1] *= inv_count;
data[n].xyz[2] *= inv_count;
}
}
float* operator[](int n)
{
return data[n].xyz;
}
int find_nearest(const float* pt_xyz) const
{
float min_dist_squared = numeric_limits<float>::max();
int min_centroid = -1;
for (int c=0; c < num_centroids; ++c)
{
const float* cen_xyz = data[c].xyz;
const float x = pt_xyz[0] - cen_xyz[0];
const float y = pt_xyz[1] - cen_xyz[1];
const float z = pt_xyz[2] - cen_xyz[2];
const float dist_squared = x*x + y*y * z*z;
if (min_dist_squared > dist_squared)
{
min_dist_squared = dist_squared;
min_centroid = c;
}
}
return min_centroid;
}
private:
Centroids(const Centroids&);
Centroids& operator=(const Centroids&);
struct Centroid
{
int count;
float xyz[3];
};
Centroid* data;
};
// A high-precision real timer would be nice, but we lack C++11 and
// the coarseness of the testing here should allow this to suffice.
static double sys_time()
{
return static_cast<double>(clock()) / CLOCKS_PER_SEC;
}
static void k_means(Points& points, Centroids& centroids)
{
// Find the closest centroid for each point.
for (int p=0; p < num_points; ++p)
{
const float* pt_xyz = points[p];
points.associate(p, centroids.find_nearest(pt_xyz));
}
// Reset the data of each centroid.
for (int c=0; c < num_centroids; ++c)
centroids.reset(c);
// Compute new position sum of each centroid.
for (int p=0; p < num_points; ++p)
centroids.sum(points.centroid(p), points[p]);
// Compute average position of each centroid.
for (int c=0; c < num_centroids; ++c)
centroids.average(c);
}
int main()
{
Points points;
Centroids centroids(points);
cout << "Starting simulation..." << endl;
double start_time = sys_time();
for (int i=0; i < num_iterations; ++i)
k_means(points, centroids);
cout << "Time passed: " << (sys_time() - start_time) << " secs" << endl;
cout << "# Points: " << num_points << endl;
cout << "# Centroids: " << num_centroids << endl;
// Write the centroids to a file to give us some crude verification
// of consistency as we make changes.
ofstream out("centroids.txt");
for (int c=0; c < num_centroids; ++c)
out << "Centroid " << c << ": " << centroids[c][0] << "," << centroids[c][1] << "," << centroids[c][2] << endl;
}
我意识到表面测试的危险，但是由于它已经被认为是以前的真实世界会议的热点，所以我希望它是可以原谅的。 我也只是对与微代码优化相关的一般技术感兴趣。
在分析这一个时，我的确得到了稍微不同的结果。 时间在这里的循环中更加均匀地分散，我不知道为什么。 也许是因为数据较小（我省略了成员，并将min_dist成员挂起并将其作为本地变量）。 质心与点之间的确切比例也有点不同，但希望足够接近以将此处的改进转换为原始代码。 它在这个表面测试中也是单线程的，而且反汇编看起来完全不同，所以我可能冒着在没有原始的情况下优化这种表面测试的风险（我现在愿意采取的风险，因为我更感兴趣扩展我的知识可以优化这些情况的技术，而不是针对这种确切情况的解决方案）。
更新Yochai Timmer的建议 - 约12.5秒
哦，我不了解组装的情况，就面临微型优化的困境。 我取而代之：
-if (min_dist_squared > dist_squared)
-{
- min_dist_squared = dist_squared;
- pt.centroid = c;
-}
有了这个：
+const bool found_closer = min_dist_squared > dist_squared;
+pt.centroid = bitselect(found_closer, c, pt.centroid);
+min_dist_squared = bitselect(found_closer, dist_squared, min_dist_squared);
只是发现时间从〜5.6秒增加到〜12.5秒。 尽管如此，这不是他的错，也没有消除他解决方案的价值 - 这是我的理解，因为他没有理解机器层面发生的事情，并在黑暗中进行刺探。 那个显然错过了，显然我并不像我最初想的那样是分支预测失误的受害者。 尽管如此，他提出的解决方案在这种情况下仍然是一个非常棒的通用功能，我很感谢将它添加到我的提示和技巧工具箱中。 现在进行第二轮。
Harold的SIMD解决方案 - 2.496秒（参见注意事项）
这个解决方案可能很棒。 将集群代表转换为SoA之后，我得到这个时间约2.5秒的时间！ 不幸的是，似乎有某种问题。 对于最终输出，我得到的结果非常不同，它表明不仅仅是轻微的精度差异，还包括一些以0为结尾的质心（意味着它们在搜索中未找到）。 我一直试图通过调试器来检查SIMD逻辑，以查看可能会出现的情况 - 这可能仅仅是我的一个转录错误，但这里是代码，以防有人发现错误。
如果错误可以在不减慢结果的情况下得到纠正，那么这种速度的提升比我以前从纯微观优化中想象的要多！
// New version of Centroids::find_nearest (from harold's solution):
int find_nearest(const float* pt_xyz) const
{
__m128i min_index = _mm_set_epi32(3, 2, 1, 0);
__m128 xdif = _mm_sub_ps(_mm_set1_ps(pt_xyz[0]), _mm_load_ps(cen_x));
__m128 ydif = _mm_sub_ps(_mm_set1_ps(pt_xyz[1]), _mm_load_ps(cen_y));
__m128 zdif = _mm_sub_ps(_mm_set1_ps(pt_xyz[2]), _mm_load_ps(cen_z));
__m128 min_dist = _mm_add_ps(_mm_add_ps(_mm_mul_ps(xdif, xdif),
_mm_mul_ps(ydif, ydif)),
_mm_mul_ps(zdif, zdif));
__m128i index = min_index;
for (int i=4; i < num_centroids; i += 4)
{
xdif = _mm_sub_ps(_mm_set1_ps(pt_xyz[0]), _mm_load_ps(cen_x + i));
ydif = _mm_sub_ps(_mm_set1_ps(pt_xyz[1]), _mm_load_ps(cen_y + i));
zdif = _mm_sub_ps(_mm_set1_ps(pt_xyz[2]), _mm_load_ps(cen_z + i));
__m128 dist = _mm_add_ps(_mm_add_ps(_mm_mul_ps(xdif, xdif),
_mm_mul_ps(ydif, ydif)),
_mm_mul_ps(zdif, zdif));
__m128i mask = _mm_castps_si128(_mm_cmplt_ps(dist, min_dist));
min_dist = _mm_min_ps(min_dist, dist);
min_index = _mm_or_si128(_mm_and_si128(index, mask),
_mm_andnot_si128(mask, min_index));
index = _mm_add_epi32(index, _mm_set1_epi32(4));
}
ALIGN16 float mdist[4];
ALIGN16 uint32_t mindex[4];
_mm_store_ps(mdist, min_dist);
_mm_store_si128((__m128i*)mindex, min_index);
float closest = mdist[0];
int closest_i = mindex[0];
for (int i=1; i < 4; i++)
{
if (mdist[i] < closest)
{
closest = mdist[i];
closest_i = mindex[i];
}
}
return closest_i;
}
哈罗德的SIMD解决方案（更正） - 约2.5秒
在应用更正并测试它们之后，结果将保持不变并且可以正常运行，并对原始代码库进行类似的改进！
由于这打击了我希望更好地了解的知识（无分支SIMD），我将用一些额外的道具来奖励解决方案，使操作速度提高一倍以上。 我试着去理解它，因为我的目标不仅仅是为了缓解这个热点，而是为了扩展我个人对可能的解决方案的理解。
尽管如此，我非常感谢所有的贡献，从算法建议到非常酷的bitselect技巧！ 我希望我能接受所有的答案。 我可能最终会在某个时候尝试所有这些，但现在我已经在理解其中一些非算术SIMD操作的作业中删除了我的作业。
int find_nearest_simd(const float* pt_xyz) const
{
__m128i min_index = _mm_set_epi32(3, 2, 1, 0);
__m128 pt_xxxx = _mm_set1_ps(pt_xyz[0]);
__m128 pt_yyyy = _mm_set1_ps(pt_xyz[1]);
__m128 pt_zzzz = _mm_set1_ps(pt_xyz[2]);
__m128 xdif = _mm_sub_ps(pt_xxxx, _mm_load_ps(cen_x));
__m128 ydif = _mm_sub_ps(pt_yyyy, _mm_load_ps(cen_y));
__m128 zdif = _mm_sub_ps(pt_zzzz, _mm_load_ps(cen_z));
__m128 min_dist = _mm_add_ps(_mm_add_ps(_mm_mul_ps(xdif, xdif),
_mm_mul_ps(ydif, ydif)),
_mm_mul_ps(zdif, zdif));
__m128i index = min_index;
for (int i=4; i < num_centroids; i += 4)
{
xdif = _mm_sub_ps(pt_xxxx, _mm_load_ps(cen_x + i));
ydif = _mm_sub_ps(pt_yyyy, _mm_load_ps(cen_y + i));
zdif = _mm_sub_ps(pt_zzzz, _mm_load_ps(cen_z + i));
__m128 dist = _mm_add_ps(_mm_add_ps(_mm_mul_ps(xdif, xdif),
_mm_mul_ps(ydif, ydif)),
_mm_mul_ps(zdif, zdif));
index = _mm_add_epi32(index, _mm_set1_epi32(4));
__m128i mask = _mm_castps_si128(_mm_cmplt_ps(dist, min_dist));
min_dist = _mm_min_ps(min_dist, dist);
min_index = _mm_or_si128(_mm_and_si128(index, mask),
_mm_andnot_si128(mask, min_index));
}
ALIGN16 float mdist[4];
ALIGN16 uint32_t mindex[4];
_mm_store_ps(mdist, min_dist);
_mm_store_si128((__m128i*)mindex, min_index);
float closest = mdist[0];
int closest_i = mindex[0];
for (int i=1; i < 4; i++)
{
if (mdist[i] < closest)
{
closest = mdist[i];
closest_i = mindex[i];
}
}
return closest_i;
}
太糟糕了，我们不能使用SSE4.1，但是非常好，SSE2就是这样。 我没有测试过这个，只是编译它，看看是否有语法错误，看看程序集是否有意义（大多数情况都是正常的，尽管即使有一些xmm寄存器没有被使用，GCC也会泄漏min_index ，但不知道为什么会发生这种情况）
int find_closest(float *x, float *y, float *z,
float pt_x, float pt_y, float pt_z, int n) {
__m128i min_index = _mm_set_epi32(3, 2, 1, 0);
__m128 xdif = _mm_sub_ps(_mm_set1_ps(pt_x), _mm_load_ps(x));
__m128 ydif = _mm_sub_ps(_mm_set1_ps(pt_y), _mm_load_ps(y));
__m128 zdif = _mm_sub_ps(_mm_set1_ps(pt_z), _mm_load_ps(z));
__m128 min_dist = _mm_add_ps(_mm_add_ps(_mm_mul_ps(xdif, xdif),
_mm_mul_ps(ydif, ydif)),
_mm_mul_ps(zdif, zdif));
__m128i index = min_index;
for (int i = 4; i < n; i += 4) {
xdif = _mm_sub_ps(_mm_set1_ps(pt_x), _mm_load_ps(x + i));
ydif = _mm_sub_ps(_mm_set1_ps(pt_y), _mm_load_ps(y + i));
zdif = _mm_sub_ps(_mm_set1_ps(pt_z), _mm_load_ps(z + i));
__m128 dist = _mm_add_ps(_mm_add_ps(_mm_mul_ps(xdif, xdif),
_mm_mul_ps(ydif, ydif)),
_mm_mul_ps(zdif, zdif));
index = _mm_add_epi32(index, _mm_set1_epi32(4));
__m128i mask = _mm_castps_si128(_mm_cmplt_ps(dist, min_dist));
min_dist = _mm_min_ps(min_dist, dist);
min_index = _mm_or_si128(_mm_and_si128(index, mask),
_mm_andnot_si128(mask, min_index));
}
float mdist[4];
_mm_store_ps(mdist, min_dist);
uint32_t mindex[4];
_mm_store_si128((__m128i*)mindex, min_index);
float closest = mdist[0];
int closest_i = mindex[0];
for (int i = 1; i < 4; i++) {
if (mdist[i] < closest) {
closest = mdist[i];
closest_i = mindex[i];
}
}
return closest_i;
}
像往常一样，它期望指针是16对齐的。 此外，填充应该与无限点（所以他们永远不会最接近目标）。
SSE 4.1会让你取代这个
min_index = _mm_or_si128(_mm_and_si128(index, mask),
_mm_andnot_si128(mask, min_index));
这样
min_index = _mm_blendv_epi8(min_index, index, mask);
这是一个asm版本，为vsyasm制作，测试了一下（似乎工作）
bits 64
section .data
align 16
centroid_four:
dd 4, 4, 4, 4
centroid_index:
dd 0, 1, 2, 3
section .text
global find_closest
proc_frame find_closest
;
; arguments:
; ecx: number of points (multiple of 4 and at least 4)
; rdx -> array of 3 pointers to floats (x, y, z) (the points)
; r8 -> array of 3 floats (the reference point)
;
alloc_stack 0x58
save_xmm128 xmm6, 0
save_xmm128 xmm7, 16
save_xmm128 xmm8, 32
save_xmm128 xmm9, 48
[endprolog]
movss xmm0, [r8]
shufps xmm0, xmm0, 0
movss xmm1, [r8 + 4]
shufps xmm1, xmm1, 0
movss xmm2, [r8 + 8]
shufps xmm2, xmm2, 0
; pointers to x, y, z in r8, r9, r10
mov r8, [rdx]
mov r9, [rdx + 8]
mov r10, [rdx + 16]
; reference point is in xmm0, xmm1, xmm2 (x, y, z)
movdqa xmm3, [rel centroid_index] ; min_index
movdqa xmm4, xmm3 ; current index
movdqa xmm9, [rel centroid_four] ; index increment
paddd xmm4, xmm9
; calculate initial min_dist, xmm5
movaps xmm5, [r8]
subps xmm5, xmm0
movaps xmm7, [r9]
subps xmm7, xmm1
movaps xmm8, [r10]
subps xmm8, xmm2
mulps xmm5, xmm5
mulps xmm7, xmm7
mulps xmm8, xmm8
addps xmm5, xmm7
addps xmm5, xmm8
add r8, 16
add r9, 16
add r10, 16
sub ecx, 4
jna _tail
_loop:
movaps xmm6, [r8]
subps xmm6, xmm0
movaps xmm7, [r9]
subps xmm7, xmm1
movaps xmm8, [r10]
subps xmm8, xmm2
mulps xmm6, xmm6
mulps xmm7, xmm7
mulps xmm8, xmm8
addps xmm6, xmm7
addps xmm6, xmm8
add r8, 16
add r9, 16
add r10, 16
movaps xmm7, xmm6
cmpps xmm6, xmm5, 1
minps xmm5, xmm7
movdqa xmm7, xmm6
pand xmm6, xmm4
pandn xmm7, xmm3
por xmm6, xmm7
movdqa xmm3, xmm6
paddd xmm4, xmm9
sub ecx, 4
ja _loop
_tail:
; calculate horizontal minumum
pshufd xmm0, xmm5, 0xB1
minps xmm0, xmm5
pshufd xmm1, xmm0, 0x4E
minps xmm0, xmm1
; find index of the minimum
cmpps xmm0, xmm5, 0
movmskps eax, xmm0
bsf eax, eax
; index into xmm3, sort of
movaps [rsp + 64], xmm3
mov eax, [rsp + 64 + rax * 4]
movaps xmm9, [rsp + 48]
movaps xmm8, [rsp + 32]
movaps xmm7, [rsp + 16]
movaps xmm6, [rsp]
add rsp, 0x58
ret
endproc_frame
在C ++中：
extern "C" int find_closest(int n, float** points, float* reference_point);
您可以使用无分支三元运算符，有时称为bitselect（condition？true：false）。
只为2名成员使用，默认为无所事事。
不要担心额外的操作，它们与if语句分支相比没有任何意义。
bitselect实现：
inline static int bitselect(int condition, int truereturnvalue, int falsereturnvalue)
{
return (truereturnvalue & -condition) | (falsereturnvalue & ~(-condition)); //a when TRUE and b when FALSE
}
inline static float bitselect(int condition, float truereturnvalue, float falsereturnvalue)
{
//Reinterpret floats. Would work because it's just a bit select, no matter the actual value
int& at = reinterpret_cast<int&>(truereturnvalue);
int& af = reinterpret_cast<int&>(falsereturnvalue);
int res = (at & -condition) | (af & ~(-condition)); //a when TRUE and b when FALSE
return reinterpret_cast<float&>(res);
}
你的循环应该看起来像这样：
for (int i=0; i < num_centroids; ++i)
{
const ClusterCentroid& cent = centroids[i];
const float dist = ...;
bool isSmaeller = dist < pt.min_dist;
//use same value if not smaller
pt.min_index = bitselect(isSmaeller, i, pt.min_index);
pt.min_dist = bitselect(isSmaeller, dist, pt.min_dist);
}
C ++是一种高级语言。 假设C ++源代码中的控制流转换为分支指令是有缺陷的。 我没有你的例子中的某些类型的定义，所以我用一个类似的条件赋值做了一个简单的测试程序：
int g(int, int);
int f(const int *arr)
{
int min = 10000, minIndex = -1;
for ( int i = 0; i < 1000; ++i )
{
if ( arr[i] < min )
{
min = arr[i];
minIndex = i;
}
}
return g(min, minIndex);
}
请注意，使用未定义的“g”只是为了防止优化器删除所有内容。 我用G ++ 4.9.2用-O3和-S将它转换成x86_64程序集（甚至不需要更改-march的默认值），并且（不太令人意外的）结果是循环体不包含分支
movl (%rdi,%rax,4), %ecx
movl %edx, %r8d
cmpl %edx, %ecx
cmovle %ecx, %r8d
cmovl %eax, %esi
addq $1, %rax
除此之外，无分支必然更快的假设也可能是有缺陷的，因为新距离“击败”旧的概率减少了你所看到的更多元素。 这不是投掷硬币。 当编译器在生成“现在”组件时比现在更不积极，因此发明了“bitselect”技巧。 我宁愿建议查看一下编译器实际生成的汇编类型，然后再尝试重写代码，以便编译器能够更好地优化它，或者将结果作为手写汇编的基础。 如果你想研究SIMD，我会建议尝试一种减少数据依赖的“最小化”方法（在我的例子中，对“min”的依赖可能是一个瓶颈）。
链接地址: http://www.djcxy.com/p/15217.html