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

Ответ 1

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

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

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

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

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

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

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

В IMHO реализация BLAS не использует алгоритмы, такие как алгоритм Coppersmith-Winograd или алгоритм Штрассена. Я не совсем уверен в причине, но это мое предположение:

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

Edit/Update:

Новый и основополагающий документ для этой темы - это статьи 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% времени. В других вариантах я использовал код intrinsics и ассемблера для повышения производительности. Вы можете увидеть учебное пособие, в котором перечислены все варианты:

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

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

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

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

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

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

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

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

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

Итак, вот вы: Кодирование оптимизированной параллельной LU-оптимизации кеша

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

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

Здесь аналогичный проект с BLAZE:

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

Ответ 2

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

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

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

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

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

• Инструкции в стиле SIMD
• Параллельность уровня инструкций
• Кэш-информирование

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

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

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

Ответ 3

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

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

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

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

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

Ответ 4

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

Ответ 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. Я думаю, что эта процедура на самом деле не Fortran, а код ассемблера (я не знаю, что есть в библиотеке, у меня нет источников). Для меня совершенно непонятно, почему cblas_dgemm не так быстро, поскольку, насколько мне известно, это всего лишь обертка для dgemm _.

Ответ 6

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

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

Ответ 7

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

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

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

Ответ 8

Мы собираемся запустить MOOC по этому вопросу на edX, мы надеемся, в начале июня 2019 года. Мы все еще собираем материалы вместе, но, если вы нетерпеливы, вы можете посмотреть записи (которые обновляются несколько раз в день). ): ulaff.net (Программирование LAFF-On для высокой производительности).

Эти материалы раскрывают методы, лежащие в основе реализации матричного умножения в матрице BLIS (и являются рефакторингом алгоритма Гото).

Пожалуйста, не связывайтесь с нами по поводу материалов (подождите, пока курс не выйдет в сеть, затем вы можете оставлять вопросы на форуме). (Будут рассмотрены исключения для миллиардеров, которые хотят финансировать усилия...)

Если вы хотите получать информацию об этих и других наших усилиях, вы можете присоединиться к группе ULAFF-On, упомянутой на этой странице.

Наслаждайтесь!

Ответ 9

Так как мой последний ответ был заблокирован, позвольте мне попробовать еще раз.

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

Чтобы сделать результирующий код переносимым на разные архитектуры, вычисления приводятся в терминах "микроядра", которое обновляет небольшой блок C (который хранится в векторных регистрах) с последовательностью обновлений ранга 1 (внешние продукты), Для этого нужно разбить задачу на умножения матрицы на матрицу с блоками, чтобы сохранить данные во всех трех уровнях кэша (от L1 до L3). Важно отметить, что нужно "упаковать" данные в непрерывную память, чтобы вычисления обращались к данным с "единичным шагом" (непрерывно).

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

• Анатомия высокопроизводительного матричного умножения К. Гото, Р. А. ван де Гейн. Транзакции ACM по математическому программному обеспечению (TOMS) 34 (3), 12. (Упоминается в предыдущих статьях и часто используется в классах)

• BLIS: платформа для быстрой реализации функциональности BLAS Ф. Г. Ван Зее, Р. А. ван де Гейн. Транзакции ACM по математическому программному обеспечению 41 (3), статья 14.

мы преодолели проблему создания 4-недельного курса под названием "Программирование для высокой производительности", который тщательно подбирает упражнения, которые ведут как новичка, так и эксперта в этом понимании.

Надеемся, что цензура Qaru поймет, что, поскольку этот курс будет предлагаться бесплатно для тех, кто хочет его проверять, а заметки будут доступны бесплатно, и разрешите указатель на то, где эти заметки можно найти: http://ulaff.net. Курс, основанный на этих заметках, будет запущен на edX в начале июня.

Это сообщение явно относится к этой дискуссии.

Ответ 10

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

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

Другое дело, что 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), если вы используете тривиальный алгоритм. Существуют алгоритмы, которые масштабируются намного лучше.