使用 AVX2 对嵌套循环进行矢量化

Question

我正在尝试将内部循环矢量化为以下嵌套循环。首先，这是一种好的做法，还是应该避免尝试矢量化嵌套循环？

下面的工作，它已经有一些基本的循环展开。

int sparsemv(struct mesh *A, const double * const x, double * const y) {

  const int nrow = (const int) A->local_nrow;
  int j = 0;
  double sum = 0.0;
  #pragma omp parallel for private(j, sum)
  for (int i=0; i< nrow; i++) {
    sum = 0.0;
    const double * const cur_vals = (const double * const) A->ptr_to_vals_in_row[i];
    const int * const cur_inds = (const int * const) A->ptr_to_inds_in_row[i];
    const int cur_nnz = (const int) A->nnz_in_row[i];
    int unroll = (cur_nnz/4)*4;
    for (j=0; j< unroll; j+=4) {
      sum += cur_vals[j] * x[cur_inds[j]];
      sum += cur_vals[j+1] * x[cur_inds[j+1]];
      sum += cur_vals[j+2] * x[cur_inds[j+2]];
      sum += cur_vals[j+3] * x[cur_inds[j+3]];
    }
    for (; j < cur_nnz; j++) {
      sum += cur_vals[j] * x[cur_inds[j]];
    }
    y[i] = sum;
  }
  return 0;
}

然而，当我尝试在 AVX2 中使用 256 位矢量寄存器进行矢量化时，我得到的要么是不正确的答案，要么是段错误。 x 和 y 对齐，但 A 未对齐，但目前，所有加载和存储都是使用未对齐操作完成的，因为这是我唯一一次没有遇到 seg 错误：

int sparsemv(struct mesh *A, const double * const x, double * const y) {

  const int nrow = (const int) A->local_nrow;
  int j = 0;
  double sum = 0.0;
  #pragma omp parallel for private(j, sum)
  for (int i=0; i< nrow; i++) {
    sum = 0.0;
    const double * const cur_vals = (const double * const) A->ptr_to_vals_in_row[i];
    const int * const cur_inds = (const int * const) A->ptr_to_inds_in_row[i];
    const int cur_nnz = (const int) A->nnz_in_row[i];
    int unroll = (cur_nnz/4)*4;
    __m256d sumVec = _mm256_set1_pd(sum);
    for (j=0; j< unroll; j+=4) {
      __m256d cur_valsVec = _mm256_loadu_pd(cur_vals + j);
      __m256d xVec = _mm256_loadu_pd(x + cur_inds[j]);
      sumVec = _mm256_add_pd(sumVec, _mm256_mul_pd(cur_valsVec, xVec));
    }
    _mm256_storeu_pd(y + i, sumVec);  // Is this storing in y + i + 1, 2 and 3 aswell?
    for (; j < cur_nnz; j++) {
      sum += cur_vals[j] * x[cur_inds[j]];
    }
    y[i] += sum;
  }
  return 0;
}