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

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

Пользовательская реализация, 10 испытаний умножения матрицы 1000×1000:

Took: 15.76542 seconds. 

BLAS, 10 испытаний умножения матрицы 1000×1000:

 Took: 1.32432 seconds. 

Это использование чисел с плавающей точкой с одинарной точностью.

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

 template 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]; } : template 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.6Ghz процессоре для BLAS выполнять 10 * 1e9 операций за 1,32 секунды? Даже если множественность была единственной операцией, и ничего больше не было сделано, она должна занять около 4 секунд.
  2. Почему моя реализация намного медленнее?

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

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

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

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

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

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

 template 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 : template 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 

Попробуй, я уверен, что ты что-то сбережешь.

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

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

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

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

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

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

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

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

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

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

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

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

  for (l=0; l 

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

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 функций. Существует много конкурирующих реализаций интерфейса.

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

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

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

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

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

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

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

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

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

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

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

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

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

Я не знаю подробно о реализации BLAS, но есть более эффективные алогорифмы для Matrix Multiplication, которые лучше, чем 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 

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

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

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

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

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

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

  • Как вы используете данные модуля Fortran 90
  • gfortran для макетов: что делает mcmodel = medium?
  • Форматированный формат файла Fortran
  • Заявление Fortran SAVE
  • Есть ли лучшее назначение двойной точности в Fortran 90?
  • Хороший профилер для Fortran и MPI
  • Как передать распределяемые массивы подпрограмм в Fortran
  • Линия усечена, Синтаксическая ошибка в списке аргументов
  • Отправка 2D-массивов в Fortran с помощью MPI_Gather
  • Как псевдоним имени функции в Fortran
  • Вызов 32-битного кода из 64-битного процесса
  • Interesting Posts

    Как запустить OpenCV-код без OpenCv Manager

    Как виртуальное наследование разрешает двусмысленность «бриллиантов» (множественное наследование)?

    Можно ли выразить «тип» lambda-выражения?

    Ошибка TypeScript в коде Angular2: не удается найти имя ‘module’

    WinForms Глобальная обработка исключений?

    Добавление локального принтера через создание стандартного порта TCP / IP и добавление сетевого принтера

    Как настроить Chrome для открытия новых экземпляров браузера в новых окнах, а не на вкладке?

    Каким образом значения хеша MD5 не являются обратимыми?

    Пользовательские предупреждения компилятора

    Класс Java Comparator для сортировки массивов

    Applescript; открытие приложения в пробеле N

    Firefox: не-Vimperator способ сделать mouseless просмотра?

    Как получить аргументы vm из java-приложения?

    Что такое идеальное соглашение об именах переменных для переменных цикла?

    Чрезмерное торрентирование и здоровье жесткого диска

    Давайте будем гением компьютера.