Как BLAS получает такую ​​экстремальную производительность?

Question

Из любопытства я решил сравнить свою функцию умножения матриц с реализацией BLAS ... Я был, мягко говоря, удивлен результатом:

Индивидуальная реализация, 10 испытаний 1000x1000 матричное умножение:

Took: 15.76542 seconds.

Реализация BLAS, 10 испытаний 1000x1000 матричное умножение:

Took: 1.32432 seconds.

Используются числа с плавающей запятой одинарной точности.

Моя реализация:

template<class ValT>
void mmult(const ValT* A, int ADim1, int ADim2, const ValT* B, int BDim1, int BDim2, ValT* C)
{
    if ( ADim2!=BDim1 )
        throw std::runtime_error("Error sizes off");

    memset((void*)C,0,sizeof(ValT)*ADim1*BDim2);
    int cc2,cc1,cr1;
    for ( cc2=0 ; cc2<BDim2 ; ++cc2 )
        for ( cc1=0 ; cc1<ADim2 ; ++cc1 )
            for ( cr1=0 ; cr1<ADim1 ; ++cr1 )
                C[cc2*ADim2+cr1] += A[cc1*ADim1+cr1]*B[cc2*BDim1+cc1];
}

У меня есть два вопроса:

1. Учитывая, что матрично-матричное умножение говорит: nxm * mxn требует n * n * m умножений, поэтому в случае выше 1000 ^ 3 или 1e9 операций. Как на моем 2,6 ГГц процессоре для BLAS можно выполнить 10 * 1e9 операций за 1,32 секунды? Даже если умножение было единственной операцией, и больше ничего не было сделано, это должно занять ~ 4 секунды.
2. Почему моя реализация намного медленнее?

Answer 1

Хорошей отправной точкой является великая книга Наука программирования матричных вычислений Роберта А. ван де Гейна и Энрике С. Кинтана-Орти. Они предоставляют бесплатную версию для загрузки.

BLAS делится на три уровня:

• Уровень 1 определяет набор функций линейной алгебры, которые работают только с векторами. Эти функции выигрывают от векторизации (например, от использования SSE).

• Функции уровня 2 - это матрично-векторные операции, например, некоторое матрично-векторное произведение. Эти функции могут быть реализованы в терминах функций уровня 1. Однако вы можете повысить производительность этих функций, если сможете обеспечить выделенную реализацию, которая использует некоторую многопроцессорную архитектуру с разделяемой памятью.

• Функции уровня 3 - это операции, подобные матрично-матричному произведению Опять же, вы можете реализовать их в терминах функций уровня 2. Но функции уровня 3 выполняют O (N ^ 3) операций с данными O (N ^ 2). Поэтому, если ваша платформа имеет иерархию кеша, вы можете повысить производительность, если предоставите выделенную реализацию, оптимизированную для кеша / дружественную для кеша . Это хорошо описано в книге. Основной импульс функций уровня 3 связан с оптимизацией кэша. Это усиление значительно превышает второе усиление от параллелизма и других аппаратных оптимизаций.

Кстати, большинство (или даже все) высокопроизводительных реализаций BLAS НЕ реализованы на Фортране. ATLAS реализован на C. GotoBLAS / OpenBLAS реализован на C, а его критически важные компоненты - на Ассемблере. Только эталонная реализация BLAS реализована на Фортране. Однако все эти реализации BLAS предоставляют интерфейс Fortran, так что он может быть связан с LAPACK (LAPACK получает всю свою производительность от BLAS).

Оптимизированные компиляторы играют незначительную роль в этом отношении (и для GotoBLAS / OpenBLAS компилятор не имеет никакого значения).

ИМХО в реализации BLAS не используются такие алгоритмы, как алгоритм Копперсмита-Винограда или алгоритм Штрассена. Я не совсем уверен в причине, но это мое предположение:

• Возможно, невозможно обеспечить оптимизированную кэш-памятью реализацию этих алгоритмов (т. Е. Вы потеряете больше, чем выиграете)
• Эти алгоритмы численно нестабильны. Поскольку BLAS является вычислительным ядром LAPACK, это не нужно.

Редактировать / Update:

Новая и новаторская статья по этой теме - BLIS paper . Они исключительно хорошо написаны. Для своей лекции «Основы программного обеспечения для высокопроизводительных вычислений» я реализовал матрично-матричный продукт, следуя их статье. На самом деле я реализовал несколько вариантов матрично-матричного продукта. Простейшие варианты полностью написаны на простом C и содержат менее 450 строк кода. Все остальные варианты просто оптимизируют циклы

    for (l=0; l<MR*NR; ++l) {
        AB[l] = 0;
    }
    for (l=0; l<kc; ++l) {
        for (j=0; j<NR; ++j) {
            for (i=0; i<MR; ++i) {
                AB[i+j*MR] += A[i]*B[j];
            }
        }
        A += MR;
        B += NR;
    }

Общая производительность матрично-матричного продукта * только 1042 * зависит от этих циклов. Около 99,9% времени проводится здесь. В других вариантах я использовал встроенный код и ассемблерный код для повышения производительности. Вы можете увидеть учебник, который рассматривает все варианты здесь:

ulmBLAS: Учебное пособие по GEMM (Matrix-Matrix Product)

Вместе с документами BLIS становится довольно легко понять, как такие библиотеки, как Intel MKL, могут добиться такой производительности. И почему не важно, используете ли вы главное хранилище строк или столбцов!

Финальные тесты здесь (мы назвали наш проект ulmBLAS):

Тесты для ulmBLAS, BLIS, MKL, openBLAS и Eigen

Другое редактирование / обновление:

Я также написал учебник о том, как BLAS используется для решения задач численной линейной алгебры, таких как решение системы линейных уравнений:

Высокая производительность LU Факторизация

(Эта факторизация LU, например, используется Matlab для решения системы линейных уравнений.)

Я надеюсь найти время для расширения учебного руководства, чтобы описать и продемонстрировать, как реализовать масштабируемую параллельную реализацию факторизации LU, как в PLASMA .

Хорошо, вот и все: Кодирование кэш-оптимизированной параллельной LU-факторизации

П.С .: Я также провел несколько экспериментов по улучшению производительности uBLAS. На самом деле довольно просто повысить (да, играть словами) производительность uBLAS:

Эксперименты на uBLAS .

Здесь похожий проект с BLAZE :

Эксперименты на BLAZE .

Answer 2

Итак, прежде всего, BLAS - это всего лишь интерфейс из 50 функций. Есть много конкурирующих реализаций интерфейса.

Во-первых, я упомяну вещи, которые в основном не связаны:

• Фортран против С, без разницы
• Продвинутые матричные алгоритмы, такие как Strassen, реализации не используют их, так как они не помогают на практике

Большинство реализаций разбивают каждую операцию на матричные или векторные операции небольшого размера более или менее очевидным способом. Например, большое умножение матрицы 1000x1000 может быть разбито на последовательность умножений матрицы 50x50.

Эти операции небольшого размера фиксированного размера (называемые ядрами) жестко запрограммированы в специфическом для процессора коде сборки с использованием нескольких функций процессора их цели:

• Инструкции в стиле SIMD
• Параллелизм уровня команд
• Cache-осведомленность

Кроме того, эти ядра могут выполняться параллельно друг другу с использованием нескольких потоков (ядер ЦП) в типичном шаблоне проектирования с уменьшением карты.

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

(Совет: вот почему, если вы используете ATLAS, вам лучше собирать и настраивать библиотеку вручную для вашей конкретной машины, чем использовать предварительно собранную.)

Answer 3

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

Во-вторых, ваш процессор может выполнять одновременно более одной инструкции.

Ваш ЦП выполняет 3-4 инструкции за цикл, и если используются блоки SIMD, каждая команда обрабатывает 4 числа с плавающей запятой или 2 двойных. (конечно, эта цифра также не точна, поскольку ЦП обычно может обрабатывать только одну SIMD-инструкцию за цикл)

В-третьих, ваш код далек от оптимального:

• Вы используете необработанные указатели, что означает, что компилятор должен предположить, что они могут иметь псевдоним. Существуют специфичные для компилятора ключевые слова или флаги, которые вы можете указать компилятору, что они не являются псевдонимами. В качестве альтернативы, вы должны использовать другие типы, кроме необработанных указателей, которые решают проблему.
• Вы очищаете кэш, выполняя наивный обход каждой строки / столбца входных матриц. Вы можете использовать блокировку, чтобы выполнить как можно больше работы с меньшим блоком матрицы, который помещается в кэш ЦП, прежде чем переходить к следующему блоку.
• Для чисто числовых задач Фортран в значительной степени непобедим, и С ++ требует много уговоров, чтобы набрать аналогичную скорость. Это можно сделать, и есть несколько библиотек, демонстрирующих это (обычно с использованием шаблонов выражений), но это не тривиально, и не происходит , просто .

Answer 4

Я не знаю конкретно о реализации BLAS, но есть более эффективные алгоритмы для умножения матриц, которые имеют сложность лучше, чем O (n3). Хорошо известным является алгоритм Штрассена

Answer 5

Большинство аргументов по второму вопросу - ассемблер, разбиение на блоки и т. Д. (Но не менее чем N ^ 3 алгоритмов, они действительно слишком развиты) - играют роль. Но низкая скорость вашего алгоритма в основном обусловлена ​​размером матрицы и неудачным расположением трех вложенных циклов. Ваши матрицы настолько велики, что не помещаются сразу в кеш-память. Вы можете переставить циклы так, чтобы как можно больше было выполнено для строки в кэше, таким образом значительно уменьшая обновления кэша (BTW-разбиение на маленькие блоки имеет аналоговый эффект, лучше всего, если циклы над блоками расположены одинаково). Реализация модели для квадратных матриц приведена ниже. На моем компьютере его время было около 1:10 по сравнению со стандартной реализацией (как у вас). Другими словами: никогда не программируйте умножение матриц по схеме «столбец времен строки», которую мы изучали в школе. После перестановки циклов можно получить больше улучшений, развернув циклы, код ассемблера и т. Д.

    void vector(int m, double ** a, double ** b, double ** c) {
      int i, j, k;
      for (i=0; i<m; i++) {
        double * ci = c[i];
        for (k=0; k<m; k++) ci[k] = 0.;
        for (j=0; j<m; j++) {
          double aij = a[i][j];
          double * bj = b[j];
          for (k=0; k<m; k++)  ci[k] += aij*bj[k];
        }
      }
    }

Еще одно замечание: эта реализация даже лучше на моем компьютере, чем замена всего на подпрограмму BLAS cblas_dgemm (попробуйте на своем компьютере!). Но гораздо быстрее (1: 4) вызывает dgemm_ из библиотеки Фортрана напрямую. Я думаю, что эта процедура на самом деле не Fortran, а ассемблерный код (я не знаю, что находится в библиотеке, у меня нет исходников). Мне совершенно непонятно, почему cblas_dgemm не такой быстрый, поскольку, насколько мне известно, это просто оболочка для dgemm _.

Answer 6

Это реалистичное ускорение. Пример того, что можно сделать с помощью ассемблера SIMD через код C ++, см. В некоторых примерах матричные функции iPhone - они были в 8 раз быстрее, чем версия C, и даже не "оптимизированы" для сборки - нет конвейерная прокладка еще и есть ненужные операции стека.

Кроме того, ваш код не является " Ограничить правильно " - как компилятор узнает, что когда он изменяет C, он не изменяет A и B?

Answer 7

Что касается исходного кода в множителе MM, ссылка на память для большинства операций является основной причиной плохой производительности. Память работает в 100-1000 раз медленнее, чем кеш.

Большая часть ускорения достигается за счет использования методов оптимизации цикла для этой функции тройного цикла в умножении MM. Используются два метода оптимизации основного цикла; Развертывание и блокировка. Что касается развертывания, мы развернем два самых внешних цикла и заблокируем их для повторного использования данных в кеше. Развертывание внешнего цикла помогает временно оптимизировать доступ к данным, уменьшая количество ссылок на память для одних и тех же данных в разное время в течение всей операции. Блокировка индекса цикла по определенному номеру помогает сохранить данные в кеше. Вы можете выбрать оптимизацию для кэша L2 или кэша L3.

https://en.wikipedia.org/wiki/Loop_nest_optimization

Answer 8

По многим причинам.

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

Другое дело, что Fortran хранит вещи по столбцам, а C - по строкам. Я не проверил ваш код, но будьте осторожны с тем, как вы выполняете продукт. В C вы должны сканировать по строкам: таким образом вы сканируете свой массив по непрерывной памяти, уменьшая количество кеш-ошибок. Промах кэша - первый источник неэффективности.

В-третьих, это зависит от используемой вами реализации blas. Некоторые реализации могут быть написаны на ассемблере и оптимизированы для конкретного процессора, который вы используете. Версия netlib написана на фортране 77.

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

Например, вы могли бы переработать свой код таким образом

template<class ValT>
void mmult(const ValT* A, int ADim1, int ADim2, const ValT* B, int BDim1, int BDim2, ValT* C)
{
if ( ADim2!=BDim1 ) throw std::runtime_error("Error sizes off");

memset((void*)C,0,sizeof(ValT)*ADim1*BDim2);
int cc2,cc1,cr1, a1,a2,a3;
for ( cc2=0 ; cc2<BDim2 ; ++cc2 ) {
    a1 = cc2*ADim2;
    a3 = cc2*BDim1
    for ( cc1=0 ; cc1<ADim2 ; ++cc1 ) {
          a2=cc1*ADim1;
          ValT b = B[a3+cc1];
          for ( cr1=0 ; cr1<ADim1 ; ++cr1 ) {
                    C[a1+cr1] += A[a2+cr1]*b;
           }
     }
  }
} 

Попробуйте, я уверен, что вы что-нибудь сохраните.

По вашему вопросу № 1 причина в том, что умножение матриц масштабируется как O (n ^ 3), если вы используете тривиальный алгоритм. Существуют алгоритмы, которые гораздо лучше масштабируются .