Оптимизация для быстрого умножения, но медленное добавление: FMA и doubleedouble

Когда я впервые получил процессор Haswell, я попробовал реализовать FMA для определения набора Мандельброта. Основной алгоритм таков:

Это определяет, находятся ли n пиксели в наборе Мандельброта. Таким образом, для двойной с плавающей запятой она работает на 4 пикселя (floatn = __m256d, intn = __m256i). Это требует 4 SIMD-операций с плавающей запятой и четырех добавок с плавающей запятой SIMD.

Затем я изменил это, чтобы работать с FMA, как это

где mul_add вызывает _mm256_fmad_pd и mul_sub вызывает _mm256_fmsub_pd. Этот метод использует 4 операции SIM-карты FMA и два SIMD-умножения, которые являются двумя менее арифметическими операциями, а затем без FMA. Кроме того, FMA и умножение могут использовать два порта и добавить только один.

Чтобы мои тесты были менее предвзятыми, я увеличил масштаб до области, которая полностью находится в наборе Мандельброта, поэтому все значения maxiter. В этом случае метод, использующий FMA, примерно на 27% быстрее.. Конечно, улучшение, но переход от SSE к AVX удвоил мою производительность, поэтому я надеялся, что, возможно, еще один фактор из двух с FMA.

Но затем я нашел этот ответ в отношении FMA, где он говорит

а затем приводит пример double * double to double-double умножение

Из этого я пришел к выводу, что я внедряю FMA не оптимально, поэтому решил реализовать SIMD double-double. Я реализовал double-double на основе Числа с плавающей запятой с расширенной точностью для вычисления графических процессоров. Бумага предназначена для двойного плавания, поэтому я изменил ее для двойного двойного. Кроме того, вместо того, чтобы упаковывать одно двойное значение в регистры SIMD, я упаковываю 4 двойных значения в один высокий регистр AVX и один низкий регистр AVX.

Для набора Мандельброта то, что мне действительно нужно, это двойное двойное умножение и добавление. В этой статье это функции df64_add и df64_mult. На рисунке ниже показана сборка для моей функции df64_mult для программного обеспечения FMA (слева) и аппаратного FMA (справа). Это ясно показывает, что аппаратное FMA является большим улучшением для двойного двойного умножения.

Итак, как аппаратное FMA выполняет в двойном двойном вычислении Мандельброта? Ответ заключается в том, что только на 15% быстрее, чем с программным обеспечением FMA. Это намного меньше, чем я надеялся. Для вычисления двойного двойника Мандельброта требуется 4 двойных двойных дополнения и четыре двойных двойных умножения (x*x, y*y, x*y и 2*(x*y)). Тем не менее, 2*(x*y) умножение тривиально для double-double, поэтому это умножение можно игнорировать в стоимости. Поэтому причина, по которой я думаю, что улучшение с использованием аппаратного FMA настолько мало, заключается в том, что в расчете преобладает медленное двойное двойное добавление (см. Сборку ниже).

Раньше было, что умножение было медленнее, чем добавление (и программисты использовали несколько трюков, чтобы избежать умножения), но с Haswell кажется, что это наоборот. Не только из-за FMA, но и потому, что умножение может использовать два порта, но только одно.

В случае, если кто-то задается вопросом, что двойной двойной метод примерно в десять раз медленнее, чем двойной. Это не так плохо, я думаю, что, если бы существовал аппаратный четырехступенчатый тип, он, вероятно, был бы как минимум в два раза медленнее, чем двойной, поэтому мой программный метод примерно в пять раз медленнее, чем я ожидал бы для аппаратного обеспечения, если бы он существовал.

Ответы

Ответ 1

Чтобы ответить на третий вопрос, я нашел более быстрое решение для двойного двойного добавления. Я нашел альтернативное определение в статье Реализация Операторы float-float на графике аппаратное обеспечение.

Theorem 5 (Add22 theorem) Let be ah+al and bh+bl the float-float arguments of the following
algorithm:
Add22 (ah ,al ,bh ,bl)
1 r = ah ⊕ bh
2 if | ah | ≥ | bh | then
3     s = ((( ah ⊖ r ) ⊕ bh ) ⊕ b l ) ⊕ a l
4 e l s e
5     s = ((( bh ⊖ r ) ⊕ ah ) ⊕ a l ) ⊕ b l
6 ( rh , r l ) = add12 ( r , s )
7 return (rh , r l)

Вот как я реализовал этот (псевдокод):

static inline doubledoublen add22(doubledoublen const &a, doubledouble const &b) {
    doublen aa,ab,ah,bh,al,bl;
    booln mask;
    aa = abs(a.hi);                //_mm256_and_pd
    ab = abs(b.hi); 
    mask = aa >= ab;               //_mm256_cmple_pd
    // z = select(cut,x,y) is a SIMD version of z = cut ? x : y;
    ah = select(mask,a.hi,b.hi);   //_mm256_blendv_pd
    bh = select(mask,b.hi,a.hi);
    al = select(mask,a.lo,b.lo);
    bl = select(mask,b.lo,a.lo);

    doublen r, s;
    r = ah + bh;
    s = (((ah - r) + bh) + bl ) + al;
    return two_sum(r,s);
}

Это определение Add22 использует 11 дополнений вместо 20, но для определения того, < <22 > , требуется некоторый дополнительный код. Ниже приведено обсуждение того, как реализовать функции SIMM minmag и maxmag. К счастью, большая часть дополнительного кода не использует порт 1. Теперь только 12 команд отправляются в порт 1 вместо 20.

Вот форма анализа пропускной способности IACA для нового Add22

Throughput Analysis Report
--------------------------
Block Throughput: 12.05 Cycles       Throughput Bottleneck: Port1

Port Binding In Cycles Per Iteration:
---------------------------------------------------------------------------------------
|  Port  |  0   -  DV  |  1   |  2   -  D   |  3   -  D   |  4   |  5   |  6   |  7   |
---------------------------------------------------------------------------------------
| Cycles | 0.0    0.0  | 12.0 | 2.5    2.5  | 2.5    2.5  | 2.0  | 10.0 | 0.0  | 2.0  |
---------------------------------------------------------------------------------------


| Num Of |                    Ports pressure in cycles                     |    |
|  Uops  |  0  - DV  |  1  |  2  -  D  |  3  -  D  |  4  |  5  |  6  |  7  |    |
---------------------------------------------------------------------------------
|   1    |           |     | 0.5   0.5 | 0.5   0.5 |     |     |     |     |    | vmovapd ymm3, ymmword ptr [rip]
|   1    |           |     | 0.5   0.5 | 0.5   0.5 |     |     |     |     |    | vmovapd ymm0, ymmword ptr [rdx]
|   1    |           |     | 0.5   0.5 | 0.5   0.5 |     |     |     |     |    | vmovapd ymm4, ymmword ptr [rsi]
|   1    |           |     |           |           |     | 1.0 |     |     |    | vandpd ymm2, ymm4, ymm3
|   1    |           |     |           |           |     | 1.0 |     |     |    | vandpd ymm3, ymm0, ymm3
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vcmppd ymm2, ymm3, ymm2, 0x2
|   1    |           |     | 0.5   0.5 | 0.5   0.5 |     |     |     |     |    | vmovapd ymm3, ymmword ptr [rsi+0x20]
|   2    |           |     |           |           |     | 2.0 |     |     |    | vblendvpd ymm1, ymm0, ymm4, ymm2
|   2    |           |     |           |           |     | 2.0 |     |     |    | vblendvpd ymm4, ymm4, ymm0, ymm2
|   1    |           |     | 0.5   0.5 | 0.5   0.5 |     |     |     |     |    | vmovapd ymm0, ymmword ptr [rdx+0x20]
|   2    |           |     |           |           |     | 2.0 |     |     |    | vblendvpd ymm5, ymm0, ymm3, ymm2
|   2    |           |     |           |           |     | 2.0 |     |     |    | vblendvpd ymm0, ymm3, ymm0, ymm2
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vaddpd ymm3, ymm1, ymm4
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vsubpd ymm2, ymm1, ymm3
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vaddpd ymm1, ymm2, ymm4
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vaddpd ymm1, ymm1, ymm0
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vaddpd ymm0, ymm1, ymm5
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vaddpd ymm2, ymm3, ymm0
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vsubpd ymm1, ymm2, ymm3
|   2^   |           |     |           |           | 1.0 |     |     | 1.0 |    | vmovapd ymmword ptr [rdi], ymm2
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vsubpd ymm0, ymm0, ymm1
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vsubpd ymm1, ymm2, ymm1
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vsubpd ymm3, ymm3, ymm1
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vaddpd ymm0, ymm3, ymm0
|   2^   |           |     |           |           | 1.0 |     |     | 1.0 |    | vmovapd ymmword ptr [rdi+0x20], ymm0

и вот анализ пропускной способности из старого

Throughput Analysis Report
--------------------------
Block Throughput: 20.00 Cycles       Throughput Bottleneck: Port1

Port Binding In Cycles Per Iteration:
---------------------------------------------------------------------------------------
|  Port  |  0   -  DV  |  1   |  2   -  D   |  3   -  D   |  4   |  5   |  6   |  7   |
---------------------------------------------------------------------------------------
| Cycles | 0.0    0.0  | 20.0 | 2.0    2.0  | 2.0    2.0  | 2.0  | 0.0  | 0.0  | 2.0  |
---------------------------------------------------------------------------------------

| Num Of |                    Ports pressure in cycles                     |    |
|  Uops  |  0  - DV  |  1  |  2  -  D  |  3  -  D  |  4  |  5  |  6  |  7  |    |
---------------------------------------------------------------------------------
|   1    |           |     | 1.0   1.0 |           |     |     |     |     |    | vmovapd ymm0, ymmword ptr [rsi]
|   1    |           |     |           | 1.0   1.0 |     |     |     |     |    | vmovapd ymm1, ymmword ptr [rdx]
|   1    |           |     | 1.0   1.0 |           |     |     |     |     |    | vmovapd ymm3, ymmword ptr [rsi+0x20]
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vaddpd ymm4, ymm0, ymm1
|   1    |           |     |           | 1.0   1.0 |     |     |     |     |    | vmovapd ymm5, ymmword ptr [rdx+0x20]
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vsubpd ymm2, ymm4, ymm0
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vsubpd ymm1, ymm1, ymm2
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vsubpd ymm2, ymm4, ymm2
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vsubpd ymm0, ymm0, ymm2
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vaddpd ymm2, ymm0, ymm1
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vaddpd ymm1, ymm3, ymm5
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vsubpd ymm6, ymm1, ymm3
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vsubpd ymm5, ymm5, ymm6
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vsubpd ymm6, ymm1, ymm6
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vaddpd ymm1, ymm2, ymm1
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vsubpd ymm3, ymm3, ymm6
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vaddpd ymm2, ymm4, ymm1
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vaddpd ymm3, ymm3, ymm5
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vsubpd ymm4, ymm2, ymm4
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vsubpd ymm1, ymm1, ymm4
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vaddpd ymm0, ymm1, ymm3
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vaddpd ymm1, ymm2, ymm0
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vsubpd ymm2, ymm1, ymm2
|   2^   |           |     |           |           | 1.0 |     |     | 1.0 |    | vmovapd ymmword ptr [rdi], ymm1
|   1    |           | 1.0 |           |           |     |     |     |     | CP | vsubpd ymm0, ymm0, ymm2
|   2^   |           |     |           |           | 1.0 |     |     | 1.0 |    | vmovapd ymmword ptr [rdi+0x20], ymm0

Лучшее решение было бы, если бы кроме FMA было три режима одиночного округления с одним операндом. Мне кажется, что для

a + b + c
a * b + c //FMA - this is the only one in x86 so far
a * b * c