На тему векторизации написано немало интересного. Вот скажем, отличный пост, который много полезного объясняет по работе автовекторизации, очень рекомендовал бы его к прочтению. Мне интересен другой вопрос. Сейчас в руках у разработчиков большое количество способов, чтобы создать «векторный» код – от чистого ассемблера до того же автовекторизатора. На каком же способе остановиться? Как найти баланс между необходимым и достаточным? Об этом и поговорим.
Итак, получить «заветные» векторные инструкции можно несколькими способами. Представим схематично в виде следующей таблицы:
Если мы опытные гуру и можем себе позволить писать на чистом ассемблере, то, пожалуй, этот путь даст нам 100% уверенности в использовании максимума в нашем коде. Ещё бы, мы сразу напишем на нужных инструкциях и используем все возможности процессора. Вот только «заточен» он будет под конкретный набор инструкций, а значит, и под конкретное «железо». Выход очередных новых инструкций (а прогресс не стоит на месте), потребует глобальной переработки и новых трудозатрат. Очевидно, стоит подумать о чём-то более легком в использовании. И на следующей «ступеньке» появляются intrinsic функции.
Это уже не чистый ассемблер, тем не менее, времени на переписывание кода всё равно потребуется немало. Скажем, простой цикл, в котором складывают два массива, будет выглядеть так:
#include <immintrin.h>
double A[100], B[100], C[100];
for (int i = 0; i < 100; i += 4) {
__m256d a = _mm256_load_pd(&A[i]);
__m256d b = _mm256_load_pd(&B[i]);
__m256d c = _mm256_add_pd(a, b);
_mm256_store_pd(&C[i], c);
}
В данном случае мы используем AVX intrinsic функции. Тем самым мы гарантируем генерацию соответствующих AVX инструкций, то есть опять привязаны к конкретному «железу». Трудозатраты уменьшились, но использовать этот код в будущем мы не сможем – рано или поздно его придётся снова переписывать. И так будет всегда, пока мы явно выбираем инструкции, «прописывая» их в исходном коде. Будь то чистый ассемблер, intrinsic функции или SIMD intrinsic классы. Тоже, кстати, интересная вещь, представляющая следующий уровень абстракции.
Тот же самый пример перепишется следующим образом:
#include <dvec.h>
// 4 elements per vector * 25 = 100 elements
F64vec4 A[25], B[25], C[25];
for(int i = 0; i < 25; i++)
C[i] = A[i] + B[i];
В этом случае нам уже не нужно знать, какие функции использовать. Код сам по себе выглядит вполне элегантно, а разработчику достаточно создавать данные нужного класса. В этом примере, F64 означает тип float размера 64 бита, а vec4 говорит об использовании Intel AVX (vec2 для SSE).
Я думаю, всем понятно, почему и этот способ нельзя назвать лучшим по соотношению «цена/качество». Правильно, портируемость всё ещё неидеальна. Поэтому разумным решением является использование компилятора для решения подобных задач. С ним, пересобирая наш код, мы сможем создавать бинарники под нужную нам архитектуру, какой бы она не была, и использовать последние наборы инструкций. При этом нам нужно убедиться, что компилятор в состоянии векторизовать код.
Пока мы шли по табличке снизу вверх, обсуждая «сложные» пути векторизации кода. Поговорим о более простых способах.
Очевидно, самый простой – переложить всю ответственность на компилятор и наслаждаться жизнью. Но не всё так просто. Каким бы умным не был компилятор, до сих пор существует множество случаев, когда он бессилен что-либо сделать с циклом без дополнительных данных или подсказок. Кроме того, в ряде случаев, код, удачно векторизованный с одной версией компилятора, уже не векторизуется с другой. Всё дело в хитрых компиляторных эвристиках, поэтому полагаться на автовекторизацию на 100% нельзя, хотя штука, безусловно, полезная. Например, современный компилятор может векторизовать такой код:
double A[1000], B[1000], C[1000], D[1000], E[1000];
for (int i = 0; i < 1000; i++)
E[i] = (A[i] < B[i]) ? C[i] : D[i];
Если бы мы пытались создать аналог кода на intrinsic функциях, гарантируя векторизацию, то получилось бы что-то подобное:
double A[1000], B[1000], C[1000], D[1000], E[1000];
for (int i = 0; i < 1000; i += 2) {
__m128d a = _mm_load_pd(&A[i]);
__m128d b = _mm_load_pd(&B[i]);
__m128d c = _mm_load_pd(&C[i]);
__m128d d = _mm_load_pd(&D[i]);
__m128d e;
__m128d mask = _mm_cmplt_pd(a, b);
e = _mm_or_pd(
_mm_and_pd (mask, c),
_mm_andnot_pd(mask, d));
_mm_store_pd(&E[i], e);
}
Хорошо, когда компилятор умеет делать это за нас! Жаль, что не всегда… и в тех случаях, когда компилятор не справляется, разработчик может помочь ему сам. Для этого можно использовать специальные «хитрости» в виде директив. Например, #pragma ivdep подскажет, что в цикле нет зависимостей, а #pragma vector always позволит не обращать внимания на «политику эффективности» векторизации (часто, если компилятор считает, что векторизовать цикл неэффективно, скажем, из-за доступа к памяти, он этого не делает). Но эти директивы из разряда «возможно помогут». Если компилятор уверен, что зависимости есть, то цикл он так и не будет векторизовывать, даже если имеется pragma ivdep.
Поэтому я выделил ещё один способ, который основан на директивах, но несколько других принципах работы. Это директивы из нового стандарта OpenMP 4.0 и Inte Cilk Plus #pragma omp simd и #pragma simd соответственно. Они позволяют совсем «забывать» компилятору о собственных проверках, и целиком полагаться на то, что говорит разработчик. Ответственность, в этом случае, естественно, перекладывается на его плечи и голову, так что действовать нужно аккуратно. Отсюда и появляется потребность в ещё одном способе.
Как бы сделать так, чтобы проверки всё же оставались, но и код гарантированно был векторизован? К сожалению, с существующим в С/С++ синтаксисе, пока никак. А вот с использованием возможностей специального синтаксиса для работы с массивами (array notation), являющегося часть Cilk Plus’a (см. предыдущий пост для того чтобы понять, как много там всего есть), это возможно. Причём синтаксис весьма простой, чем то напоминает Фортран, и имеет следующий вид:
base[first:length:stride]
Задаём имя, начальный индекс, число элементов, шаг(опционален) и вперёд. Предыдущий пример перепишется с ним так:
double A[1000], B[1000], C[1000], D[1000], E[1000];
E[:] = (A[:] < B[:]) ? C[:] : D[:];
Двоеточие означает, что мы обращаемся ко всем элементам. Можно так же осуществлять более сложные манипуляции. Скажем этот код
for (i = 0; i < 5; i++)
A[(i*2) + 1] = B[(i*1) + 1];
перепишется более компактно, и главное, гарантирует векторизацию:
int A[10], *B;
A[1:5:2] = B[1:5];
Таким образом, мы видим, что способов добиться векторизации действительно много. Если говорить о балансе, то раз уж мы перечислили все способы в табличке «от простого к сложному», то с точки зрения требуемых затрат и результата на выходе, золотая середина сходится на Cilk Plus’е. Но это не значит, что всё так очевидно. Если пациент болен, ему же не всегда сразу выписывают антибиотики, верно? Так же и здесь. Для кого-то может оказаться достаточно автовектризации, для кого-то директивы ivdep и vector always будут вполне разумным решением. Важнее вовлекать компилятор, чтобы не болела голова при выходе новых инструкций и «железа», а здесь у Intel всегда есть что предложить. Так что до новых постов, друзья!
