AVX: матричный вектор точек, но игнорировать диагональ

Question

Мне нужно выполнить скалярное произведение между квадратной матрицей и вектором. Однако диагональ всегда должна игнорироваться во время этой конкретной операции. Я делаю это с AVX.

Как я могу изменить свой существующий код, чтобы эффективно игнорировать все [i, i] места, чтобы он оставался кеш-памяти?

template<bool add_to_result=false>
inline void dot( const float *const f_vecStart,  size_t fVecCount,  float *result ) const {
    m2d_assert(fVecCount == _y);
    assert_is_16_aligned(f_vecStart);
    assert_is_16_aligned(result);

    float *f_rowStart = _matArray;
    float *f_rowEnd = _matArray + _y;

    for (size_t i=0; i<_x; ++i){ //for every row:
        const __m256 *a = (__m256*)f_rowStart;
        const __m256 *rowEnd = (__m256*)(f_rowEnd);
        const __m256 *val_vec = (__m256*)(f_vecStart);

        __m256 mRowSum = _mm256_set1_ps(0.0f);

        while(a<rowEnd){
            __m256 mul = _mm256_mul_ps(*a, *val_vec);
            mRowSum = _mm256_add_ps(mRowSum, mul);

            ++val_vec;  ++a;
        }

        // finally, horizontally add the gathered sum (m256 vector), completing the
        // computation for this entire row:
        if (add_to_result){//<--known at compile time.
            result[i] += Mathf::fast_hAdd_ps(mRowSum);
        }
        else {
            result[i] = Mathf::fast_hAdd_ps(mRowSum);//<--notice EQUALS (overwrite garbage values)
        }

        f_rowStart += _y;
        f_rowEnd += _y;
    }//end for every row

    check_isNan(result, _x);
}

Answer 1

Если ошибка округления FP не является проблемой, вы, конечно, можете просто вычесть vec[i] * matrow[i] из этого точечного произведения после l oop. Но это терпит неудачу, если это приводит к NaN или Inf, или к огромному числу, которое приводит к огромной ошибке округления для остального скалярного произведения.

---

Для небольших матриц, вы можете рассмотреть возможность добавления некоторого вида безответственное маскирование на внутреннюю l oop например, вектор {0,1,2,3,4,5,6,7}, который вы увеличиваете с +=8, и vpcmpeqd против set1_epi32(i). Используйте это как маску с ANDN.

Для больших матриц вы хотите что-то, что только добавляет накладные расходы к внешнему l oop, а не к внутреннему l oop.

You может временно сделать соответствующий исходный векторный элемент 0.0, если это не const. Это «только» приведет к остановке пересылки магазина в первых рядах; последующие должны иметь время для скалярного хранилища, чтобы покинуть буфер хранилища до того, как перезагрузка вектора достигнет этой точки. (Особенно, если вы сделаете это как можно раньше, например, перед суммой предыдущей строки.)

Конечно, если вы можете обнулить диагональ матрицы, это тоже работает. Но гораздо хуже пространственное расположение для этих магазинов. (Делать это прямо перед чтением этой строки было бы наименее вредно для локальности на больших матрицах)

В противном случае вы могли бы разделить внутренний l oop на два диапазона, используя векторную маску для одного вектора, который охватывает диагональ. Это может ввести в заблуждение некоторые отрасли. Если вы выполняете распараллеливание и работаете на многопоточном компьютере, SMT может поддерживать производительность. В противном случае вам следует надеяться, что он хорошо предсказывает.

Кстати, вы можете собирать от 2 до 8 строк параллельно, развернув l oop, чтобы скрыть задержку FP и дать вам место для перемешивания и -хадд в уборке. Транспонирование и суммирование 4 или 8 векторов дешевле, чем 8x хсумов до одного вектора. (Но это проблема для изменения фактического исходного вектора).

---

Я не проверял ни одного из них; это только с моей головы.

Answer 2

Вы можете применить vblendps к результату умножения для диагональных блоков. Это очень легко при развертывании 8 л oop итераций. Развернув, вы также сохраняете перезагрузку вашего вектора каждый раз. После этого просто умножьте + добавьте к нему оставшиеся блоки (с FMA, если доступно). И, наконец, уменьшите векторные результаты до одного вектора (это один из немногих случаев, когда haddps может быть полезным).

Вот версия, которая предполагает, что ваш размер кратен 8 (а ваш матрица квадратная). Он либо добавляет к существующему вектору result, либо «добавляет» к неявному нулевому вектору. Я добавил _matArray в качестве параметра функции, и я предполагаю, что _x, _y и fVecCount в вашем примере одинаковы. Не проверено / протестировано ( ", но он компилируется" ) и достаточно много кода для копирования + вставки (вы можете попытаться автоматически развернуть части этого компилятора). Что касается локальности кэша, может быть лучше развернуть 4 строки (8 *1024* 4B = 32 КБ, что больше, чем кэш L1 Piledriver). Или может помочь добавление некоторых предварительных выборок.

void dot(float* result, float const* _matArray, float const* f_vecStart) {
  bool const add_to_result = true; // can also be a template parameter
  int const rowlength = 1024;      // can also be a function parameter
  for (int row = 0; row < rowlength; row += 8) {
    float const* blockStart = _matArray + rowlength * row;
    // registers for accumulation
    __m256 s0, s1, s2, s3, s4, s5, s6, s7;

    // calculate diagonal-block:
    {
      // register to blend out or replace value by pre-existing value
      __m256 res0 = add_to_result ? _mm256_loadu_ps(result + row) : _mm256_setzero_ps();

      __m256 vec_r = _mm256_loadu_ps(f_vecStart + row);
      s0 = _mm256_mul_ps(vec_r, _mm256_loadu_ps(blockStart + row + 0 * rowlength));
      s1 = _mm256_mul_ps(vec_r, _mm256_loadu_ps(blockStart + row + 1 * rowlength));
      s2 = _mm256_mul_ps(vec_r, _mm256_loadu_ps(blockStart + row + 2 * rowlength));
      s3 = _mm256_mul_ps(vec_r, _mm256_loadu_ps(blockStart + row + 3 * rowlength));
      s4 = _mm256_mul_ps(vec_r, _mm256_loadu_ps(blockStart + row + 4 * rowlength));
      s5 = _mm256_mul_ps(vec_r, _mm256_loadu_ps(blockStart + row + 5 * rowlength));
      s7 = _mm256_mul_ps(vec_r, _mm256_loadu_ps(blockStart + row + 7 * rowlength));
      s6 = _mm256_mul_ps(vec_r, _mm256_loadu_ps(blockStart + row + 6 * rowlength));

      // replace diagonal product by zero or pre-existing result value
      s0 = _mm256_blend_ps(s0, res0, 1 << 0);
      s1 = _mm256_blend_ps(s1, res0, 1 << 1);
      s2 = _mm256_blend_ps(s2, res0, 1 << 2);
      s3 = _mm256_blend_ps(s3, res0, 1 << 3);
      s4 = _mm256_blend_ps(s4, res0, 1 << 4);
      s5 = _mm256_blend_ps(s5, res0, 1 << 5);
      s6 = _mm256_blend_ps(s6, res0, 1 << 6);
      s7 = _mm256_blend_ps(s7, res0, 1 << 7);
    }

    // add remaining elements
    for (int col = 0; col < rowlength; col += 8) {
      // skip diagonal block. Maybe it is worth splitting the loop into two halves
      if (col == row) continue;

      __m256 vec_r = _mm256_loadu_ps(f_vecStart + col);
      s0 = _mm256_fmadd_ps(vec_r, _mm256_loadu_ps(blockStart + col + 0 * rowlength), s0);
      s1 = _mm256_fmadd_ps(vec_r, _mm256_loadu_ps(blockStart + col + 1 * rowlength), s1);
      s2 = _mm256_fmadd_ps(vec_r, _mm256_loadu_ps(blockStart + col + 2 * rowlength), s2);
      s3 = _mm256_fmadd_ps(vec_r, _mm256_loadu_ps(blockStart + col + 3 * rowlength), s3);
      s4 = _mm256_fmadd_ps(vec_r, _mm256_loadu_ps(blockStart + col + 4 * rowlength), s4);
      s5 = _mm256_fmadd_ps(vec_r, _mm256_loadu_ps(blockStart + col + 5 * rowlength), s5);
      s6 = _mm256_fmadd_ps(vec_r, _mm256_loadu_ps(blockStart + col + 6 * rowlength), s6);
      s7 = _mm256_fmadd_ps(vec_r, _mm256_loadu_ps(blockStart + col + 7 * rowlength), s7);
    }

    // reduce s0-s7 horizontally and store
    {
      // Perhaps on Piledriver doing vshufps+blend is more efficient?
      __m256 s01 = _mm256_hadd_ps(s0, s1);
      __m256 s23 = _mm256_hadd_ps(s2, s3);
      __m256 s45 = _mm256_hadd_ps(s4, s5);
      __m256 s67 = _mm256_hadd_ps(s6, s7);

      __m256 s0123 = _mm256_hadd_ps(s01, s23);
      __m256 s4567 = _mm256_hadd_ps(s45, s67);

      // inter-lane reduction
      // combine upper half of s0123 with lower half of s4567:
      __m256 res = _mm256_permute2f128_ps(s0123, s4567, 0x21);
      // blend lower half of s0123 with upper half of s4567 and add:
      res = _mm256_add_ps(res, _mm256_blend_ps(s0123, s4567, 0xF0));

      // store result. Ideally, replace by store_ps, if you can guarantee 32byte alignment
      _mm256_storeu_ps(result + row, res);
    }
  }
}