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

Из любопытства я решил сравнить свою функцию умножения матриц с реализацией 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. Почему моя реализация намного медленнее?

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

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

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

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

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

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

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

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

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

Изменить / обновить:

Новым и новаторским документом по этой теме являются документы BLIS . Они исключительно хорошо написаны. Для моей лекции «Основы программного обеспечения для высокопроизводительных вычислений» я реализовал матрично-матричное произведение после их статьи. Фактически я реализовал несколько вариантов матричного произведения. Простейшие варианты полностью написаны на простом 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;
    }

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

ulmBLAS: Учебник по GEMM (матрично-матричный продукт)

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

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

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

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

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

Факторизация LU с высокими характеристиками

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

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

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

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

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

Вот похожий проект с BLAZE :

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

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

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

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

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

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

• Инструкции в стиле SIMD
• Параллелизм на уровне инструкций
• Осведомленность о кэше

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

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

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

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

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

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

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

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

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

Большинство аргументов ко второму вопросу - ассемблер, разбиение на блоки и т.д. (но не менее 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_.

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

Также ваш код не является « ограничивающим правильным » - как компилятор узнает, что когда он изменяет C, он не изменяет A и B?

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

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

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

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

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

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

В-третьих, это зависит от используемой вами реализации blas. Некоторые реализации могут быть написаны на ассемблере и оптимизированы для конкретного процессора, который вы используете. Версия netlib написана на fortran 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), если вы используете тривиальный алгоритм. Есть алгоритмы, которые масштабируются намного лучше .