Программа Simd Matmul дает разные численные результаты

Question

Я пытаюсь запрограммировать умножение матриц в C, используя встроенные функции simd. Я был почти уверен в своей реализации, но когда я выполняю, я получаю некоторые числовые ошибки, начиная с 5-й цифры полученных коэффициентов матрицы.

REAL_T - это просто float с typedef

/* This is my matmul Version with simd, using floating simple precision*/
void matmul(int n, REAL_T *A, REAL_T *B, REAL_T *C){
  int i,j,k;
  __m256 vA, vB, vC, vRes;
  for (i=0; i<n; i++){
    for (j=0; j<n; j++){  
      for (k=0; k<n; k= k+8){
        vA = _mm256_load_ps(&A[i*n+k]);
        vB = _mm256_loadu_ps(&B[k*n+j]);
        vC = _mm256_mul_ps(vA, vB);
        vC = _mm256_hadd_ps(vC, vC);
        vC = _mm256_hadd_ps(vC, vC);
        /*To get the resulting coefficient, after doing 2 hadds,
        I have to get the first and the last element of the resulting
        Vector vC*/
        C[i*n+j] += ((float )(vC[0])) + ((float )(vC[7]));
      } /* for k */
    } /* for j */
  } /* for i */
}
*/End of program

/*And this is the sequential Version*/
void matmul(int n, REAL_T *A, REAL_T *B, REAL_T *C){
  int i,j,k;
  for (i=0; i<n; i++){ 
    for (j=0; j<n; j++){
      for (k=0; k<n; k++){
        C[i*n+j] +=  A[i*n+k] *  B[k*n+j];  
      } /* for k */
    } /* for j */
  } /* for i */  
}
/*End of program*/

/*The matrix are initialized as follows*/
  for (i = 0; i < n; i++)
    for (j = 0; j < n; j++){
      *(A+i*n+j) = 1 / ((REAL_T) (i+j+1));
      *(B+i*n+j) = 1.0;
      *(C+i*n+j) = 1.0;
    }
/*End of initialization*/

Протестированная матрица имеет размер 512 * 512. Для последовательной версии левый верхний квадрат полученной матрицы дает:

+6.916512e+01  +6.916512e+01  
+5.918460e+01  +5.918460e+01  

+7.946186e+00  +7.946186e+00  
+7.936391e+00  +7.936391e+00  

Однако, для версии simd, квадрат:

+6.916510e+01  +6.916510e+01  
+5.918463e+01  +5.918463e+01  

+7.946147e+00  +7.946147e+00  
+7.936355e+00  +7.936355e+00 

Как показано, числовая ошибка между двумя версиями. Любая помощь будет принята с благодарностью!

Answer 1

Это выглядит нормально; добавление чисел в другом порядке приводит к разным округлениям во временных значениях.

FP математика не ассоциативна; оптимизация, как будто это изменит результаты. ¹ Являются ли сложения с плавающей запятой и умножение ассоциативными? / Являются ли операции с плавающей запятой в C ассоциативными?

Количество изменений зависит от данных. Различия только в 5-м десятичном знаке кажутся разумными для float.

---

Если вы не предприняли специальных численных мер предосторожности, таких как сложение первых чисел, результат в последовательном порядке не является «более правильным», они просто имеют разные ошибки.

Фактически, при использовании нескольких аккумуляторов обычно увеличивает точность для больших списков, предполагая, что все ваши числа имеют одинаковую величину. (В идеале несколько векторов SIMD, каждый из которых состоит из нескольких элементов, чтобы скрыть задержку FP-add или FMA). https://en.wikipedia.org/wiki/Pairwise_summation - числовая техника, которая выводит это на следующий уровень: суммирование подмножеств списка в дереве, чтобы избежать добавления отдельных элементов массива к гораздо большему значению. См., Например, Как избежать менее точной суммы для массивов с несколькими столбцами

Использование фиксированного количества аккумуляторов (например, 8x __m256 = 64 float аккумуляторов) может уменьшить ожидаемую ошибку в 64 раза вместо N до log N для полного попарного суммирования.

---

Сноска 1: Ассоциативность необходима для распараллеливания, SIMD и нескольких аккумуляторов. Ассоциативность дает нам возможность распараллеливания. Но что дает коммутативность?

На машине, например, с 4-тактовой задержкой FMA с пропускной способностью 2 на такт, с шириной SIMD 8 поплавков, т. Е. Системой Skylake с AVX2, потенциальное ускорение 4 * 2 = 8 от нескольких аккумуляторов, * 8 от ширины SIMD, умноженного на количество ядер, по сравнению с чисто последовательной версией, даже для задач, где может быть менее точным, чем просто другим.

Большинство людей считают, что фактор 8*8 = 64 стоит того! (И теоретически вы можете распараллелить для еще одного множителя, равного 4, на четырехъядерном процессоре, предполагая идеальное масштабирование для больших матриц).

Вы уже используете float вместо double для производительности.

См. Также Почему mulss занимает только 3 цикла в Haswell, в отличие от таблиц инструкций Агнера? для получения дополнительной информации об использовании нескольких аккумуляторов, чтобы скрыть задержку FMA в сокращении, выставляя этот другой фактор ускорения в 8 раз.

Кроме того, не используйте hadd внутри самого внутреннего цикла. Суммируйте по вертикали и используйте эффективное сокращение в конце цикла. ( Самый быстрый способ сделать горизонтальную векторную сумму с плавающей запятой на x86 ). Вы действительно хотите избежать того, чтобы компилятор извлекал ваши векторы для скалярного вычисления на каждом шагу, что лишает большинство преимуществ SIMD! Помимо того, что hadd не стоит использовать для горизонтальных сумм 1 вектора; стоит 2 шаффла + обычный add на всех существующих процессорах.