如何在SSE优化C++代码时避免if条件分支？

针对你这段带双层粒子循环+位移循环的代码，核心痛点是if (nx != dx || ny != dy)这个条件分支会打断SIMD流水线，除了用SSE逻辑比较生成掩码的方法外，还有以下几种实用方案：

1. 掩码扩展：直接清零自交互贡献 #

不用分支判断，先把所有粒子对（包括自交互）的贡献都算出来，再用掩码把自交互的项置零。这种方法和你提到的逻辑比较思路一致，但更直接适配你的计算逻辑：

// 以SIMD打包处理多组(dx, dy)为例
__m128d nx_vec = _mm_set1_pd((double)nx);
__m128d dx_vec = _mm_loadu_pd(dx_batch); // 打包当前要处理的dx值
__m128d ny_vec = _mm_set1_pd((double)ny);
__m128d dy_vec = _mm_loadu_pd(dy_batch);

// 生成掩码：非自交互的位置为1.0，自交互为0.0
__m128d nx_neq_dx = _mm_cmpneq_pd(nx_vec, dx_vec);
__m128d ny_neq_dy = _mm_cmpneq_pd(ny_vec, dy_vec);
__m128d valid_mask = _mm_or_pd(nx_neq_dx, ny_neq_dy);

// 打包计算rlen、vmag、delta等
__m128d rlen_vec = _mm_set1_pd(sqrt((double)((nx-dx)*(nx-dx)+(ny-dy)*(ny-dy)))*sep);
__m128d xdiff_vec = _mm_sub_pd(_mm_loadu_pd(x + dx*n + dy), _mm_set1_pd(x[index]));
__m128d ydiff_vec = _mm_sub_pd(_mm_loadu_pd(y + dx*n + dy), _mm_set1_pd(y[index]));
__m128d zdiff_vec = _mm_sub_pd(_mm_loadu_pd(z + dx*n + dy), _mm_set1_pd(z[index]));
__m128d vmag_vec = _mm_sqrt_pd(_mm_add_pd(_mm_add_pd(_mm_mul_pd(xdiff_vec, xdiff_vec), _mm_mul_pd(ydiff_vec, ydiff_vec)), _mm_mul_pd(zdiff_vec, zdiff_vec)));
__m128d delta_vec = _mm_sub_pd(vmag_vec, rlen_vec);

// 计算势能和力的增量，并用掩码过滤自交互项
__m128d pe_inc = _mm_mul_pd(_mm_mul_pd(delta_vec, delta_vec), _mm_set1_pd(fcon));
pe_inc = _mm_mul_pd(pe_inc, valid_mask);

__m128d coeff = _mm_mul_pd(_mm_set1_pd(fcon), _mm_div_pd(delta_vec, vmag_vec));
__m128d fx_inc = _mm_mul_pd(_mm_mul_pd(coeff, xdiff_vec), valid_mask);
__m128d fy_inc = _mm_mul_pd(_mm_mul_pd(coeff, ydiff_vec), valid_mask);
__m128d fz_inc = _mm_mul_pd(_mm_mul_pd(coeff, zdiff_vec), valid_mask);

// 累加结果到pe和fx/fy/fz
pe += _mm_cvtsd_f64(_mm_hadd_pd(pe_inc, pe_inc)); // 水平加法汇总SIMD结果
_mm_storeu_pd(fx + index, _mm_add_pd(_mm_loadu_pd(fx + index), fx_inc));
// fy、fz同理

2. 拆分循环：从根源消除分支 #

直接调整dx、dy的循环范围，让循环里只处理非自交互的粒子对，彻底去掉if判断：

// 处理dx < nx的所有dy
for (dy = MAX(ny - delta, 0); dy < MIN(ny + delta + 1, n); dy++) {
    for (dx = MAX(nx - delta, 0); dx < nx; dx++) {
        calculate_interaction(nx, ny, dx, dy, pe, fx, fy, fz, ...);
    }
}

// 处理dx == nx时，dy != ny的情况
if (nx >= MAX(nx - delta, 0) && nx <= MIN(nx + delta + 1, n)) {
    // dy < ny
    for (dy = MAX(ny - delta, 0); dy < ny; dy++) {
        calculate_interaction(nx, ny, nx, dy, pe, fx, fy, fz, ...);
    }
    // dy > ny
    for (dy = ny + 1; dy < MIN(ny + delta + 1, n); dy++) {
        calculate_interaction(nx, ny, nx, dy, pe, fx, fy, fz, ...);
    }
}

// 处理dx > nx的所有dy
for (dy = MAX(ny - delta, 0); dy < MIN(ny + delta + 1, n); dy++) {
    for (dx = nx + 1; dx < MIN(nx + delta + 1, n); dx++) {
        calculate_interaction(nx, ny, dx, dy, pe, fx, fy, fz, ...);
    }
}

把循环拆成三个部分后，每个子循环里的粒子对都满足非自交互，完全不需要分支判断，CPU流水线能持续高效运行，SIMD也能充分发挥作用。

3. 预计算rlen表：减少重复计算+适配SIMD #

rlen只和(nx-dx)、(ny-dy)的相对位置有关，和粒子坐标无关，可以提前预计算所有可能的rlen值，循环内直接查表：

// 提前预计算rlen表（假设delta是固定值，或每次函数调用前生成）
double rlen_table[2*delta+1][2*delta+1];
for (int ddy = -delta; ddy <= delta; ddy++) {
    for (int ddx = -delta; ddx <= delta; ddx++) {
        if (ddx != 0 || ddy != 0) {
            rlen_table[ddy + delta][ddx + delta] = sqrt((double)(ddx*ddx + ddy*ddy)) * sep;
        } else {
            rlen_table[ddy + delta][ddx + delta] = 0.0; // 自交互用不到
        }
    }
}

// 循环内查表获取rlen
int ddx = dx - nx;
int ddy = dy - ny;
double rlen = rlen_table[ddy + delta][ddx + delta];

预计算后，循环内少了sqrt和乘法运算，配合前面的掩码或循环拆分方法，能进一步提升SIMD的执行效率。

4. 融合乘加（FMA）指令优化计算 #

虽然不是直接消除分支，但FMA可以把fcon * xdiff * (vmag - rlen) / vmag这类复合运算合并成单条指令，减少延迟，配合掩码使用效果更好：

__m128d delta_over_vmag = _mm_div_pd(delta_vec, vmag_vec);
__m128d coeff = _mm_mul_pd(_mm_set1_pd(fcon), delta_over_vmag);
// FMA指令：fx_inc = coeff * xdiff_vec + 0.0
__m128d fx_inc = _mm_fmadd_pd(coeff, xdiff_vec, _mm_setzero_pd());
// 应用掩码过滤自交互
fx_inc = _mm_mul_pd(fx_inc, valid_mask);

FMA指令能同时完成乘法和加法，减少指令数和流水线延迟，让SIMD单元的利用率更高。

内容的提问来源于stack exchange，提问作者Beute