Привет, Хабр! Меня зовут Евгений Буевич, я работаю в Рунити. Как-то раз у меня возникла непреодолимая потребность умножать матрицы определенного размера, смотреть, что получится и умножать опять до тех пор, пока что-нибудь не получится. =) Дело, как я понял после знакомства с литературой, любимое многими.

Остановился на BLIS, скомпилировал, подключил, и было мне счастье. Матрицы стали подрастать в числе и размере, скорость процесса, как ей и положено, падала в кубе от размера и кратно от числа. В конце концов стало ощущаться, что на ЦПУ 486,4 GFLOPS и ни флопсом больше, а замеры показывали, что на самом деле их около 350. Оно бы и ладно, но и на одном треде их было не больше 107. Стало интересно, куда пропадают остальные. Наибольшую ясность в вопрос внесла статья «Умножение матриц: эффективная реализация шаг за шагом» , после которой стало более понятно, что происходит внутри библиотеки. Так совпало, что мой процессор похож на используемый в ней, и пальцы зачесались что-нибудь улучшить.

Вкратце о том, как это работает

Матрицы умножаются блочно, площадь блока определяется числом векторных регистров. В данном случае их 16, из которых 4 заняты, размеры блока взяты 6x16 (6 строк, 2 регистра в ширину).

Данные B собираются в буфер Bbuf, предположительно находящийся в кеше L1 так, чтобы они позже читались последовательно:

Аналогично данные A собираются в другой буфер Abuf, предположительно находящийся в L2, порциями, соответствующими размеру Bbuf (назовем такую порцию микропанель Abuf):

Алгоритм состоит из трех (по числу измерений матриц) + число задействованных уровней кеша вложенных циклов. Самый внутренний цикл умножает mr строк матрицы A из Abuf на nr столбцов матрицы B из Bbuf. Он называется микроядром и по идеологии BLIS должен быть написан на ассемблере для конкретной архитектуры. В таблице — краткая характеристика всех 6 циклов (удобнее в порядке снизу вверх):

Уровень размещения данных определяется частотой доступа. В таком подходе Bbuf при умножении на микропанели Abuf будет переиспользован чаще всего и удерживается в L1. Реже будет использоваться Abuf, умножаемый на разные Bbuf. C L3 теоретически все не очень хорошо, поскольку при определенном размере матриц обращения к C могут вытеснить данные ранее использованного Bbuf до того, как алгоритм к нему вернется. Это решаемо, но такой размер матриц нам и не нужен, остановимся на скромном и кратном всем мыслимым блокам 1152х1152

Все это намного более подробно описано здесь, откуда я частично заимствовал стартовую точку для экспериментов. Также можно почитать статьи из документации к BLIS.

В поисках гигафлопсов

После более-менее близкого воспроизведения кода получился примерно тот же результат: 99.2 GFLOPS.

void micro_core_6x16(matrixtype_t alpha,unsigned B_height, unsigned C_width, 
                     const matrixtype_t *A, const matrixtype_t *B, matrixtype_t *C)
   {
   unsigned k;
   __m256 c00 = _mm256_setzero_ps(), c01 = _mm256_setzero_ps();
   __m256 c10 = c00, c11 = c01, c20 = c00, c21 = c01;
   __m256 c30 = c00, c31 = c01, c40 = c00, c41 = c01, c50 = c00, c51 = c01;
   __m256 b0, b1, a0, a1;

   for (k = 0; k < B_height; k++)
      {
      b0 = _mm256_loadu_ps(B + 0); b1 = _mm256_loadu_ps(B + ITEMS_PER_REGISTRY);   
      B += MICROCORE_WIDTH * ITEMS_PER_REGISTRY; 

      a0 = _mm256_set1_ps(A[0]); a1 = _mm256_set1_ps(A[1]);                  
      c00 = _mm256_fmadd_ps(a0, b0, c00); c01 = _mm256_fmadd_ps(a0, b1, c01);
      c10 = _mm256_fmadd_ps(a1, b0, c10); c11 = _mm256_fmadd_ps(a1, b1, c11);
      a0 = _mm256_set1_ps(A[2]); a1 = _mm256_set1_ps(A[3]);                  
      c20 = _mm256_fmadd_ps(a0, b0, c20); c21 = _mm256_fmadd_ps(a0, b1, c21);
      c30 = _mm256_fmadd_ps(a1, b0, c30); c31 = _mm256_fmadd_ps(a1, b1, c31);
      a0 = _mm256_set1_ps(A[4]); a1 = _mm256_set1_ps(A[5]);                  
      c40 = _mm256_fmadd_ps(a0, b0, c40); c41 = _mm256_fmadd_ps(a0, b1, c41);
      c50 = _mm256_fmadd_ps(a1, b0, c50); c51 = _mm256_fmadd_ps(a1, b1, c51);
      A += MICROCORE_HEIGHT;
      }
   a0 = _mm256_set1_ps(alpha);
   b0 = _mm256_fmadd_ps(a0,c00,_mm256_loadu_ps(C + 0)); 
   _mm256_storeu_ps(C + 0, b0);
   b1 = _mm256_fmadd_ps(a0,c01,_mm256_loadu_ps(C + ITEMS_PER_REGISTRY)); 
   _mm256_storeu_ps(C + ITEMS_PER_REGISTRY, b1);
   C += C_width;
   b0 = _mm256_fmadd_ps(a0,c10,_mm256_loadu_ps(C + 0)); 
   _mm256_storeu_ps(C + 0, b0);
   b1 = _mm256_fmadd_ps(a0,c11,_mm256_loadu_ps(C + ITEMS_PER_REGISTRY)); 
   _mm256_storeu_ps(C + ITEMS_PER_REGISTRY, b1);
   C += C_width;
   b0 = _mm256_fmadd_ps(a0,c20,_mm256_loadu_ps(C + 0)); 
   _mm256_storeu_ps(C + 0, b0);
   b1 = _mm256_fmadd_ps(a0,c21,_mm256_loadu_ps(C + ITEMS_PER_REGISTRY)); 
   _mm256_storeu_ps(C + ITEMS_PER_REGISTRY, b1);
   C += C_width;
   b0 = _mm256_fmadd_ps(a0,c30,_mm256_loadu_ps(C + 0)); 
   _mm256_storeu_ps(C + 0, b0);
   b1 = _mm256_fmadd_ps(a0,c31,_mm256_loadu_ps(C + ITEMS_PER_REGISTRY)); 
   _mm256_storeu_ps(C + ITEMS_PER_REGISTRY, b1);
   C += C_width;
   b0 = _mm256_fmadd_ps(a0,c40,_mm256_loadu_ps(C + 0)); 
   _mm256_storeu_ps(C + 0, b0);
   b1 = _mm256_fmadd_ps(a0,c41,_mm256_loadu_ps(C + ITEMS_PER_REGISTRY)); 
   _mm256_storeu_ps(C + ITEMS_PER_REGISTRY, b1);
   C += C_width;
   b0 = _mm256_fmadd_ps(a0,c50,_mm256_loadu_ps(C + 0)); 
   _mm256_storeu_ps(C + 0, b0);
   b1 = _mm256_fmadd_ps(a0,c51,_mm256_loadu_ps(C + ITEMS_PER_REGISTRY)); 
   _mm256_storeu_ps(C + ITEMS_PER_REGISTRY, b1); 
   }

Код целиком (и последующие оптимизации тоже) можно посмотреть тут. 

Попробуем поискать, куда пропадают гигафлопсы.

На первый взгляд, наибольшие задержки дают обращения к данным C. Исходя из этого, микроядра меньше используют данные B и больше данные A, чтобы цикл, определяемый размером Bbuf в L1 кеше, был длиннее. Добавим предвыборки для C:

...
   __m256 b0, b1, a0, a1;
   _mm_prefetch((const char *)&C[0 * C_width],PT_LEVEL);
   _mm_prefetch((const char *)&C[1 * C_width],PT_LEVEL);
   _mm_prefetch((const char *)&C[2 * C_width],PT_LEVEL);
   _mm_prefetch((const char *)&C[3 * C_width],PT_LEVEL);
   _mm_prefetch((const char *)&C[4 * C_width],PT_LEVEL);
   _mm_prefetch((const char *)&C[5 * C_width],PT_LEVEL);

   for (k = 0; k < B_height; k++)
...

Результат подрос до 101,8 GFLOPS, L3 промахов стало меньше.

Если поглядеть микроядра в BLIS, видно что причина не в этом: предвыборки там есть, но результат далек от максимального. Единственное, что вызывает некоторое сомнение — запрошенный уровень L1 перед проходящим по L1 циклом, хотя этому есть теоретическое объяснение (ниже в блоке про используемый объем кеша).

Второе очевидное место поиска — заполнение буферов упорядоченными данными. Несомненно, оно отнимает какое-то полезное время, и, судя по этой статье, авторов библиотеки это тоже беспокоит. Вопрос в том, какое. При размере 1152х1152 матрица B поместится в L3 целиком, и таким образом будет переупорядочена один раз. Матрица A также будет переупорядочена 1 раз. Полное умножение-сложение всех элементов требует 95М тактов. Примем эту величину за 100%. На переупорядочивание матрицы B потребуется 165К тактов, но переупорядочивание А содержит вычисления, вызванные необходимостью собрать данные из 6 разных строк матрицы в векторном регистре (10 unpckxps). Они выполняются на единственном 5-м порту, поэтому потребуют 552К тактов.

void fill_micro_A_6x16(int A_width, int micro_A_height, int micro_A_width, 
                       int micro_A_step, const matrixtype_t *A, 
                       matrixtype_t *micro_A)
   {
   int i,k;
   for (i = 0; i < micro_A_height; i += MICROCORE_HEIGHT)
      {
      matrixtype_t *mA = micro_A;
      for (k = 0; k < micro_A_width; k += 4)
         {
         const float * pA = A + k;
         __m128 a0 = _mm_loadu_ps(pA + 0 * A_width);
         __m128 a1 = _mm_loadu_ps(pA + 1 * A_width);
         __m128 a2 = _mm_loadu_ps(pA + 2 * A_width); 
         __m128 a3 = _mm_loadu_ps(pA + 3 * A_width);
         __m128 a4 = _mm_loadu_ps(pA + 4 * A_width); 
         __m128 a5 = _mm_loadu_ps(pA + 5 * A_width);
         __m128 a00 = _mm_unpacklo_ps(a0, a2), a01 = _mm_unpacklo_ps(a1, a3);
         __m128 a10 = _mm_unpackhi_ps(a0, a2), a11 = _mm_unpackhi_ps(a1, a3);
         __m128 a20 = _mm_unpacklo_ps(a4, a5), a21 = _mm_unpackhi_ps(a4, a5);
         _mm_storeu_ps(mA + 0 * MICROCORE_HEIGHT, _mm_unpacklo_ps(a00, a01));
         _mm_storel_pi((__m64*)(mA + 0 * MICROCORE_HEIGHT + 4), a20);
         _mm_storeu_ps(mA + 1 * MICROCORE_HEIGHT, _mm_unpackhi_ps(a00, a01));
         _mm_storeh_pi((__m64*)(mA + 1 * MICROCORE_HEIGHT + 4), a20);
         _mm_storeu_ps(mA + 2 * MICROCORE_HEIGHT, _mm_unpacklo_ps(a10, a11));
         _mm_storel_pi((__m64*)(mA + 2 * MICROCORE_HEIGHT + 4), a21);
         _mm_storeu_ps(mA + 3 * MICROCORE_HEIGHT, _mm_unpackhi_ps(a10, a11));
         _mm_storeh_pi((__m64*)(mA + 3 * MICROCORE_HEIGHT + 4), a21);
         mA += MICROCORE_HEIGHT * 4;
         }
      micro_A += micro_A_step;
      A += MICROCORE_HEIGHT * A_width;
      }
   }

Нашли ~0,8% простоя. Ниже будет много отсылок к номерам портов, поэтому напомню важные для нас: 0,1 – FMA инструкции, 2,3 – загрузка из памяти, 4 – сохранение в память, 5 – векторные перестановки.

Чтение происходит из основной памяти, поэтому можно добавить предвыборки для A и здесь, поскольку нас ограничивает порт 5, на порты 2 и 3 есть возможность что-нибудь послать бесплатно. Аналогично для B — ограничивает порт 4, есть возможность для предвыборок (для небольших матриц у B должна включиться аппаратная предвыборка, но согласно Intel Optimization Manual она не умеет шаги больше 2Кб). Однако для первых циклов заполнения это не так просто: предвыборки должны быть где-то внутри микроядра. Это достаточно неудачно для концепции BLIS разделения кода на си фреймворк и специфичное для архитектуры микроядро. По крайней мере остальные итерации циклов можно ускорить. Получилось примерно 104 GFLOPS.

Теперь посмотрим зависимости по данным. Микроядро считает 12 переменных, каждая из которых зависит от своего предыдущего значения [в регистре], одного из двух значений одной матрицы и одного из 6 значений другой матрицы. Два векторных порта выполняют 12 операций за 6 тактов. Единственная неустранимая зависимость между итерациями —это предыдущее значение регистра аккумулятора, а смещения в массивах можно увеличивать один раз в несколько итераций через разворот цикла. Загрузка данных для умножения может начаться еще на предыдущей итерации (если считать временем итерации фактическое выполнение умножений, а не получение микроопераций из очереди) для всех итераций, кроме первой. На первой мы получим простой векторных портов в ожидании данных. Достаточно интересный вопрос — латентность инструкций vbroadcastss и vmovups, в который транслируются интрисики из кода. По таблицам Agner Fog она равна минимум 3 циклам, по Intel Optmization Manual — 1+4 циклам, по uops.info — <=5;<=8. Судя по экспериментам в последнем источнике, под меньшими значениями подразумевается латентность самой инструкции без доступа в память, потому что данные получаются через store forwarding. Значит, если мы озаботимся предвыборкой первого значения Abuf в L1 и найдем, где ее поместить, в начале цикла потеряются 8+1 тактов.

Для этого участка у меня получилась следующая таблица (из предположения что микрооперации из очереди поступают по 4 за такт, как в предыдущих процессорах):

Здесь RS — ожидание аргументов или портов, L — выполнение загрузки из памяти, V — выполнение FMA, RT — получение результата. Красные линии показывают готовность аккумулятора С_xy на следующем цикле. Видно, что порты V полностью загружены на каждом такте, кроме первой итерации.

Еще минимум 6 тактов потеряются в конце микроядра, поскольку ограничивающим будет единственный порт сохранения данных.

Таким образом, на каждый вызов микроядра есть некоторые потери. Число этих вызовов обратно пропорционально числу итераций в микроядре или размеру Bbuf, или доле кеша, которую мы под него выделяем. Про это стоит рассказать более подробно.

Нам нужно добиться, чтобы данные Abuf не вытесняли данные Bbuf, поэтому использование полного объема L1 для Bbuf — не лучшая идея. Также где-то надо размещать и данные C, и если взаимным расположением упорядоченных Abuf и Bbuf мы можем управлять, то адреса данных C произвольны. Тут на выручку приходит ассоциативность кеша, равная 8 на Skylake. Каждая кешлиния может попасть в 1 из 8 мест в кеше (ассоциативный набор), и это место выбирается по алгоритму LRU (если точно, то Tree-PLRU, и это не совсем одно и тоже). По данным Bbuf микроядро пробегает последовательно, таким образом нам нужно выделить под него столько кеш-линий, чтобы за цикл микроядра первое значение (Bbuf[0]) не было вытеснено. На итерацию приходится 6 значений Abuf и 16 значений Bbuf. Примерно в таком соотношении и надо распределять кеш. Для лучшей управляемости этим процессом микропанели Abuf выравниваются по странице памяти (страница памяти соответствует условному «слою» L1 — данные одной страницы не могут вытеснить друг друга). Вопрос возникает про данные C: стоит ли выделять им кеш-линию в каждом наборе. Теоретически да, на практике результаты варианта 0,75 иногда лучше, чем 0,625, так как некоторое число L1 промахов может окупаться уменьшением количества вызовов.

Для определенности, остановимся на 0,625. При этом в нулевом наборе L1 к моменту второго вызова микроядра должно получиться что-то вроде бутерброда на картинке.

Вызовов, если использовать эту часть объема L1, будет 55К, 15 тактов на вызов — итого теряются 825K тактов. Также можно предположить, что 6 последовательных предвыборок (порты 2 и 3) перед 4 загружающими инструкциями в начале микроядра дадут еще 3 такта простоя, что добавит 165К. Всего оверхеда таким образом насчиталось еще около процента. Эти выкладки можно проверить, выкинув лишнее из вызывающей функции — заполнение буферов и сдвиги по ним. Код получится примерно таким:

void one_6x16_empty(const FCalcParams *params)
   {
   unsigned i,j,k,n,m,ni;

   for(n = 0; n < params->n_size; n += params->n_step)
      { 
      unsigned npart = min(params->n_size, n + params->n_step) - n;
      for(k = 0; k < params->k_size; k += params->k_step)
         { 
         unsigned micro_B_height = min(params->k_size, k + params->k_step) - k;
         for (m = 0; m < params->m_size; m += params->m_step)
            { 
            unsigned micro_A_height = min(params->m_size, m + params->m_step) - m;
            for (j = 0; j < npart; j += MICROCORE_WIDTH * ITEMS_PER_REGISTRY)
               for (i = 0, ni = 0; i < micro_A_height; i += MICROCORE_HEIGHT, ni++) 
                  micro_core_6x16(params->alpha,micro_B_height,params->n_size,
                        &params->micro_A[0], &params->micro_B[0],&params->C[0]);
            }
         }
      }
   }

Лучшее, что мне из него удалось извлечь — 120GFLOPS, что меньше ожидаемого (98,6% а не 99%). Если убрать предвыборки (так как данные C во всех вызовах одинаковы, они помогают только в первом цикле), результат еще чуть-чуть улучшится.

Дайте две

Очевидно, что остальные 11% оверхеда — это кеш-промахи. И тут можно поcтупить, как авторы BLIS и других библиотек, отлавливая их и подбирая место для предвыборки. Возможно, по дороге придется разбираться, почему Tree-PLRU не равно LRU. Но есть альтернативный путь — Simultaneous Multithreading (SMT), более известный под именем HyperThreading (его первой реализации Intel).

Идея в том, что, когда одному потоку нечем занять векторные (и вообще любые) порты, это делает второй поток. Интуитивно кажется оптимальным сделать Bbuf общим, а Abuf разделить на части и выделить каждую своему потоку, что даст максимальную длину цикла микроядра (ну и использовать пошаренный L1 кэш — воплощенная мечта). Однако придется обойти ряд подводных камней.

Первый из них — синхронизации. Признаться, я сильно переживал из-за них, поскольку мутекс слишком дорого, а спинлок тратит ресурсы второго потока. Но все оказалось не так страшно, синхронизироваться надо только возле заполнения буферов. В конечной реализации 1 цикл спинлока приходится примерно на 40 вызовов микроядра. И из-за большой латентности PAUSE это не простой, а почти полноценная работа в 1 поток.

Второй подводный камень: если использовать 0,675 L1 — получим кучу кеш-промахов. Потоков стало 2 и каждый поднимает свои данные Abuf. Таким образом нужно 5 кеш-линий на B, 4 на А и еще одну, а всего их 8. Если же распределить пропорционально, используемая доля кеша ощутимо уменьшится. Решительно переворачиваем микроядро на бок — используем 6 значений B и всего 2 из A. Соотношение A к B стало 2:48, теперь можно использовать 0,75 или даже 0,875 L1, оставив единственную кеш-линию в наборе для A. Приложив некоторые усилия (а именно выкинув выравнивание по страницам), можно обеспечить отсутствие коллизий между потоками. Упрощенно код такого микроядра выглядит так:

void micro_core_2x48(matrixtype_t alpha,unsigned B_height, unsigned C_width, 
                     const matrixtype_t *A, const matrixtype_t *B, matrixtype_t *C)
   {
   unsigned k;
   __m256 c00 = _mm256_setzero_ps(), c01 = _mm256_setzero_ps();
   __m256 c02 = _mm256_setzero_ps(), c03 = _mm256_setzero_ps();
   __m256 c04 = _mm256_setzero_ps(), c05 = _mm256_setzero_ps();
   __m256 c10 = c00, c11 = c01, c12 = c02, c13 = c03, c14 = c04, c15 = c05; 
   __m256 b0, b1, a0, a1;

   for (k = 0; k < B_height; k++)
      {
      a0 = _mm256_set1_ps(A[0]); a1 = _mm256_set1_ps(A[1]);   
      A += MICROCORE_HEIGHT_2;
      b0 = _mm256_loadu_ps(B + ITEMS_PER_REGISTRY * 0);      
      b1 = _mm256_loadu_ps(B + ITEMS_PER_REGISTRY * 1);   
      c00 = _mm256_fmadd_ps(a0, b0, c00); c10 = _mm256_fmadd_ps(a1, b0, c10);
      c01 = _mm256_fmadd_ps(a0, b1, c01); c11 = _mm256_fmadd_ps(a1, b1, c11);
      b0 = _mm256_loadu_ps(B + ITEMS_PER_REGISTRY * 2);                     
      b1 = _mm256_loadu_ps(B + ITEMS_PER_REGISTRY * 3);   
      c02 = _mm256_fmadd_ps(a0, b0, c02); c12 = _mm256_fmadd_ps(a1, b0, c12);
      c03 = _mm256_fmadd_ps(a0, b1, c03); c13 = _mm256_fmadd_ps(a1, b1, c13);
      b0 = _mm256_loadu_ps(B + ITEMS_PER_REGISTRY * 4);                     
      b1 = _mm256_loadu_ps(B + ITEMS_PER_REGISTRY * 5);   
      c04 = _mm256_fmadd_ps(a0, b0, c04); c14 = _mm256_fmadd_ps(a1, b0, c14);
      c05 = _mm256_fmadd_ps(a0, b1, c05); c15 = _mm256_fmadd_ps(a1, b1, c15);
      B += MICROCORE_WIDTH_2 * ITEMS_PER_REGISTRY; 
      }
   b0 = _mm256_set1_ps(alpha);
   b1 = _mm256_fmadd_ps(b0,c00, _mm256_loadu_ps(C + ITEMS_PER_REGISTRY * 0));
   _mm256_storeu_ps(C + ITEMS_PER_REGISTRY * 0, b1);
   a1 = _mm256_fmadd_ps(b0,c01, _mm256_loadu_ps(C + ITEMS_PER_REGISTRY * 1));
   _mm256_storeu_ps(C + ITEMS_PER_REGISTRY * 1, a1);
   b1 = _mm256_fmadd_ps(b0,c02, _mm256_loadu_ps(C + ITEMS_PER_REGISTRY * 2));
   _mm256_storeu_ps(C + ITEMS_PER_REGISTRY * 2, b1);
   a1 = _mm256_fmadd_ps(b0,c03, _mm256_loadu_ps(C + ITEMS_PER_REGISTRY * 3));
   _mm256_storeu_ps(C + ITEMS_PER_REGISTRY * 3, a1);
   b1 = _mm256_fmadd_ps(b0,c04, _mm256_loadu_ps(C + ITEMS_PER_REGISTRY * 4));
   _mm256_storeu_ps(C + ITEMS_PER_REGISTRY * 4, b1);
   a1 = _mm256_fmadd_ps(b0,c05, _mm256_loadu_ps(C + ITEMS_PER_REGISTRY * 5));
   _mm256_storeu_ps(C + ITEMS_PER_REGISTRY * 5, a1);
   C += C_width;
   c00 = _mm256_fmadd_ps(b0,c10, _mm256_loadu_ps(C + ITEMS_PER_REGISTRY * 0));
   _mm256_storeu_ps(C + ITEMS_PER_REGISTRY * 0, c00);
   c01 = _mm256_fmadd_ps(b0,c11, _mm256_loadu_ps(C + ITEMS_PER_REGISTRY * 1));
   _mm256_storeu_ps(C + ITEMS_PER_REGISTRY * 1, c00);
   c02 = _mm256_fmadd_ps(b0,c12, _mm256_loadu_ps(C + ITEMS_PER_REGISTRY * 2));
   _mm256_storeu_ps(C + ITEMS_PER_REGISTRY * 2, c02);
   c03 = _mm256_fmadd_ps(b0,c13, _mm256_loadu_ps(C + ITEMS_PER_REGISTRY * 3));
   _mm256_storeu_ps(C + ITEMS_PER_REGISTRY * 3, c03);
   c04 = _mm256_fmadd_ps(b0,c14, _mm256_loadu_ps(C + ITEMS_PER_REGISTRY * 4));
   _mm256_storeu_ps(C + ITEMS_PER_REGISTRY * 4, c04);
   c05 = _mm256_fmadd_ps(b0,c15, _mm256_loadu_ps(C + ITEMS_PER_REGISTRY * 5));
   _mm256_storeu_ps(C + ITEMS_PER_REGISTRY * 5, c05);
   }

Это порождает следующую проблему — теперь данные Bbuf «расходуются» в три раза быстрее, а значит вызовов микроядра и соответствующих простоев тоже стало в три раза больше. Но действительно ли это проблема? В таблице ниже показано, что происходит, когда два идентичных потока выполняют команды на ограничивающем порту (у нас это FMA порт, для упрощения можно представить что он один), при этом вероятность использования этого порта одним потоком p, а сдвиг времени старта одного потока относительно другого равномерно распределен.

Тогда загрузку порта двумя потоками можно выразить как

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

У нас p = 106,7 / 121,6 = 0,877, отсюда p' ≥ 0,991. Или утилизация (0,991 – 0,877) / (1 – 0,877)  ≈ 92,6% простоев. Если число потерь при вызове цикла в однопоточной версии обозначить как L, то с двумя потоками и ядром 2x48 оно будет L' = L * 3 * 0,074 = 0.222L

Получилась формула «универсального выпрямления рук» — про точно расположенные предвыборки и многие другие сложности однопоточной реализации, требующие глубокого погружения в архитектуру процессора, можно спокойно забыть. Достаточно обеспечить более-менее равномерный сдвиг итераций циклов в потоках. Схематично это выглядит так:

Интересно, что это не такая простая задача, nop, кроме дырки в очереди микроинструкций, ни к чему не приводит, а рассчитать, сколько нужно подряд таких дырок, чтобы заданное количество тактов ничего не делать, — довольно сложно. Я добавил после вызова микроядра 4 бестолковые предвыборки данных, которые уже заведомо в L1, и это дало некоторый видимый эффект. Хотя синхронизация на спинлоке сама по себе рандомизирует смещение потоков в пределах минимум 140 тактов (латентность современной pause довольно загадочна, но по интеловскому руководству не менее 140 тактов, так что можно ориентироваться на это значение).

Размеры микроядра менялись, чтобы уменьшить число вызовов, но в результате число вызовов только увеличилось. Можно съехидничать, что в этом суть современного мира, но если число вызовов для SMT реализации не так важно, может, просто использовать микроядро 6x16? Однако у 2x48 есть пара побочных достоинств. Во-первых, при заполнении Abuf мы избавляемся от громоздкой конструкции на 5-м порту — теперь чередовать надо данные всего двух строк, что снижает нагрузку на него (по формуле «выпрямления рук» это мелочь, но все же). Во-вторых, код vmovups короче на байт, чем vbroadcastss. И это уже не мелочь.

Потому что главный подводный камень приходится прямо в киль: при наличии двух потоков в ядре часть ресурсов фронтенда режется пополам. В Skylake это и выборка инструкций (16 байт за такт, по 8 на поток) и кеш микроопераций, который должен спасти, если она не справляется. Loopback-буфер, насколько я понял, в нем поломали и не починили. В сочетании с длинными векторными инструкциями (4–7 байт) все это приводит к периодическим простоям конвейера. Вручную решать пазл максимального процента укладки двух инструкций в 8 байт не хотелось, заставить это сделать компилятор не удалось. Более того, он, видимо, разворачивает циклы под полный кеш микроопераций и после прохода кода заполнения буферов повторно собирает простои, вытеснив микроядро из обрезанного кеша. Но в версии 2x48 это происходит значительно реже, именно потому что vmovups (которых в ней шесть штук, а не две), на 1 байт короче чем vbroadcastss (которых в ней две штуки, а не шесть).

Еще одним моментом остается необходимость как-то разнести по времени заполнение буферов и умножение. Для Abuf это решилось само собой через рандомизацию сдвига потоков. Поскольку каждый поток заполняет свою часть Abuf, тут не требуется синхронизация, и они через несколько итераций разъезжаются достаточно, чтобы фазы заполнения микропанелей Abuf, нужных для следующего цикла, не перекрывались. С Bbuf все сложнее. Если его заполнять пополам каждым потоком, то это будет впритык к точке синхронизации, то есть перекрытие практически неизбежно. Но можно схитрить и делать это одним потоком, выделив ему чуть меньшую долю Abuf и таким образом сократив в нем число вызовов микроядра. Получится какое-то количество L1 промахов, так как в L1 одновременно будут и данные B для нового Bbuf и данные старого Bbuf. Однако в это время заполняющий поток не будет поднимать данные Abuf, а без предвыборок собрать Bbuf — не такое быстрое дело: опытным путем установлено что второй поток успевает за это время сделать 9–10 вызовов микроядра. Данные B по большей части заменят в L1 отсутствующую половину данных Abuf. Для оставшихся промахов есть формула выше.

Поскольку проект стал из рабочего домашним и в какой-то момент плавно мигрировал под Windows, вызов микроядра пришлось заинлайнить, чтобы избавиться от массивных и бесполезных в этом месте виндовых пролога и эпилога. Итоговый результат 118,1 из 121,6 GFLOPS на ядре. При этом предвыборки C в L2, предвыборок A и B нет (чтобы попусту не занимать порты), ассемблера нет. Пасьянс из регистров в последней итерации цикла им не считается.

Бочка дегтя в ложке меда

Вроде бы взяли детально изученный алгоритм умножения матриц и разогнали его процентов на 10. Выглядит неплохо, осталось решить пару небольших проблем:

1. Для совместного использования L1 потоки должны работать на одном ядре. Это решается через привязку к ядрам. Проблема в том, что они должны это делать одновременно. К потерям можно добавить ожидание второго потока на синхронизации, пока то, что работает вместо него, сносит наш кеш, и, собственно, восстановление кеша. В принципе, однопоточная реализация тоже чувствительна к засорению кеша, разница в масштабах.

1². тоже самое под гипервизором. Гипервизоры стали умнее и не делят одно ядро между гостевыми системами (подозреваю, что им помогли преодолеть этот рубеж различные уязвимости, а не забота о производительности). Но как получить маппинг видимых ядер на физические из гостевой системы? Под WSL2 после долгого копания в CPUID так и не удалось конвенционным путем понять, на чем работают треды (бит наличия гипервизора при этом погашен). Можно пойти "проверенным" путем: через тайминг последовательного доступа к одной ячейке памяти из обоих потоков, но выглядит это печально.

2. В режиме SMT значительная доля процессоров что-то режет пополам. Компиляторы в такие тонкости не вникают и оптимизируют под полные ресурсы (в частности, разворот циклов под кеш микроопераций). Я не нашел способа подсказать MSVC или GCC, что в определенных функциях он будет уполовинен.

В процессе работы сложилось ощущение что SMT — жертва подхода SMT=SMP, или «в два раза больше ядер за ту же цену». Если бы была возможность сообщить ОС, что такие-то потоки предназначены для совместной работы на одном физическом ядре, все было бы проще.

Вывод

Одновременной многопоточностью можно пользоваться в контролируемой среде, где она может дать статистически заметное преимущество и даже сохранит код относительно читаемым. Умножение матриц, возможно, не лучший пример, поскольку основная фаза, полностью загружающая ограничивающий порт, за много лет эволюции алгоритма выросла до 90+% от общего времени выполнения (хотя AMX опять меняет ограничивающий фактор этого процесса на скорость памяти, и тут есть над чем подумать).