Векторизация вычислений (SIMD): аппаратное ускорение численных алгоритмов
Когда закон Мура перестает ускорять код
На протяжении десятилетий программисты и математики привыкли к простому правилу: чтобы программа для численного интегрирования или умножения матриц заработала быстрее, нужно просто подождать пару лет, пока Intel или AMD выпустят процессор с более высокой тактовой частотой. Этот феномен, подпитываемый законом Мура, позволял писать неоптимальный код и все равно получать рост производительности. Однако в начале 2000-х годов рост тактовых частот уперся в непреодолимый физический барьер — тепловой предел кремниевых чипов (частота застыла на уровне 3-5 ГГц). Чтобы продолжать наращивать вычислительную мощь, создателям процессоров пришлось менять архитектуру, внедряя многоядерность и, что еще более важно для математики, векторные инструкции.
Концепция SIMD (Single Instruction, Multiple Data — одна инструкция, множество данных) является фундаментом высокопроизводительных вычислений (HPC). В классической, скалярной модели процессора (SISD), чтобы сложить два массива из четырех чисел (A1+B1, A2+B2, A3+B3, A4+B4), процессору нужно выполнить цикл: загрузить два числа, сложить их, записать результат, и повторить это 4 раза. Архитектура SIMD внедряет в процессор сверхширокие регистры. Загрузив сразу 4 числа из массива A в один 256-битный регистр, и 4 числа из массива B во второй регистр, процессор за один единственный такт времени выполняет сложение сразу всех четырех пар чисел одновременно! Это обеспечивает мгновенный 4-кратный (или даже 8-кратный) прирост чистой математической производительности.
Наборы инструкций AVX и выравнивание памяти
Современные процессоры оснащены мощнейшими блоками векторной математики. Инструкции SSE (128 бит) позволяют складывать 4 числа типа float за такт. Инструкции AVX и AVX2 (256 бит) удваивают этот показатель до 8 чисел. Новейшие серверные процессоры с инструкциями AVX-512 оперируют колоссальными 512-битными регистрами, выполняя по 16 операций сложения или умножения 32-битных чисел с плавающей запятой за один микроскопический такт процессора. Более того, существует специальная инструкция FMA (Fused Multiply-Add), которая за один такт выполняет сразу две операции (A * B + C), удваивая пиковую производительность алгоритма матричного умножения.
Но эта фантастическая аппаратная мощь недоступна программам автоматически. Чтобы векторные блоки заработали, данные в оперативной памяти компьютера должны быть расположены (упакованы) строго последовательно, без разрывов, и желательно быть выровненными (aligned) по границам 32-байтовых или 64-байтовых кеш-линий. Если вы написали алгоритм метода конечных разностей, который скачет по памяти (собирая данные из разных классов или объектов), процессор просто не сможет загрузить их широким вектором. Он откатится к медленным скалярным вычислениям (потеряв до 90% скорости).
Компиляторы, раскрутка циклов и библиотеки BLAS
Чтобы заставить численный код летать, разработчики применяют агрессивную раскрутку циклов (Loop Unrolling) и прагмы для компиляторов (например, #pragma omp simd в OpenMP). Важнейшим требованием для векторизации является отсутствие логических ветвлений (операторов if/else) внутри тяжелого математического цикла, так как ветвление разрушает векторный конвейер.
В мире профессиональной вычислительной математики никто не пишет циклы умножения матриц вручную. Вся базовая алгебра (нахождение скалярных произведений, умножение матрицы на вектор, разложение матриц) делегируется стандартизированным библиотекам BLAS (Basic Linear Algebra Subprograms) и LAPACK. Такие реализации, как Intel MKL (Math Kernel Library) или OpenBLAS, написаны на глубоком ассемблерном уровне. Они досконально учитывают размер L1/L2 кэша процессора и микроархитектуру AVX регистров, нарезая матрицы на микро-блоки. Именно благодаря идеальной SIMD-векторизации внутри библиотек BLAS, скрипты на языке Python (использующие NumPy) умудряются работать со скоростью скомпилированного языка C++ при решении гигантских систем дифференциальных уравнений.