SGEMM на видеокарте и CPU, серия 6

Уже становится какой-то дурной традицией переделывать тест умножения матриц каждые несколько месяцев. Однако с прошлого забега случилось несколько событий:

  • В форумах по CUDA появилась очень быстрая реализация SGEMM, сделанная Василием Волковым.
  • Вышел CUDA Toolkit 2.0 beta, использующий реализацию Волкова в библиотеке CuBLAS.
  • 4-ядерный Penryn стал не просто доступен, а появился у меня прямо под столом.

Участники забега

Hardware

  • Intel Core2Quad 9300, на частоте 3GHz (работает отлично, температура 49-59 градусов), память DDR2 PC6400 4-4-4-12
  • NVidia GeForce 8800GTX на штатной частоте
Software
  • Windows XP Professional x64
  • Intel MKL 10.0.3 (использовалась совместно с MS Visual C 2005 в 32-битном режиме)
  • Goto BLAS 1.26 (использовалась под Cygwin, 32-битный режим)
  • NVidia CUDA 2.0 beta (32-битная версия)

Методика тестирования

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

Для видеокарты жизнь чуть сложнее, там есть время потраченное на вычисления и дополнительное время, нужное для загрузки-выгрузки данных, поэтому для CUDA приводятся два значения производительности.

Результаты

Производительность умножения матриц (SGEMM), GFLOP/s
Описание теста размер задачи
  1024 2048 4096
Intel MKL 10.0.3 54.8 59.5 60.5
Goto BLAS 1.26 76.7 83.4 87.1
NVidia CUDA 2.0 beta, чистое время вычислений 202 204 205
NVidia CUDA 2.0 beta, вычисления+ввод-вывод 151 177 190

Обсуждение результатов

Сначала проговорим сами результаты:

  1. Penryn - феноменально быстрый. Надо смотреть не только на абсолютные цифры производительности (которые тоже большие), но и в количество реальных операций на такт. Предельная (теоретическая) производительность у этого процессора - 8 операций с одинарной точностью на ядро * 4 ядра * 3Ghz = 96 GFLOP/s, а мы получили 87 или больше 90%. И это все - за $350
  2. Несмотря на анонсированные в 10-й версии Intel MKL ускорения, Goto BLAS так и остается раза в полтора быстрее, чем библиотека от производителя.
  3. CUDA 2.0 ускорилась относительно первой версии в 1.7 раза на вычислениях и в 1.9 раза на вычислениях+I/O (вклад в скорость I/O могла дать и замена hardware хост-компьютера, хотя по тестам bandwidth я принципиального ускорения не вижу).

Удивительно, но произошло следущее

  • Стоимость гигафлопса на CPU упала с октября прошлого года примерно втрое (для десктопных процессоров) и достигла стоимости гигафлопса на GPU полугодовалой давности.
  • При этом, стоимость гигафлопса (только по одному тесту) на GPU тоже упала почти вдвое при том, что сам GPU остался прежним и подешевел несильно. Падение стоимости исключительно софтовое.
  • Конечно, 8800GTX - далеко не лидер по price-performance сейчас, новые карты могут на некоторых задачах и обгонять (но надо проверять).
  • Впрочем, соотношение в очередной раз изменится буквально в ближайшие недели, ходят упорные слухи, что новый чип от NVidia всех опять порвет.

Tags: 

Comments

GTX 280 будет примерно в 1.5

GTX 280 будет примерно в 1.5 раза быстрее, чем 9800GTX. печально, что это достигается банальным увеличением количества вычислительных блоков. чип становится больше, дороже и горячее. рано или поздно этот подход упрётся в ограничения техпроцесса. неспроста стали появляться двухчиповые видюхи, значит у производителя пока нет других идей по увеличению производительности.

По блокам там 224/128 = 1.75

По блокам там 224/128 = 1.75 раза. Частоту забыл, может и чуть ниже.
Но сам путь - наращиванием блоков (и ядер)- мне кажется куда менее тупиковым, чем наращивание частоты.

Всё идёт по кругу

Повторят историю с 8800 -- сначала сделают прожорливый чип (8800GTX), потом освоят новый техпроцесс и выпустят сравнимую по производительности, но более холодную и дешёвую версию (8800GTS на G92, он же 9800GTX).

Тоже G92 (он на этом чипе был

Тоже G92 (он на этом чипе был изначально), но с урезанным числом исполнительных блоков.

Тоже G92 (он на нём был

Тоже G92 (он на нём был изначально), но с урезанным числом исполнительных блоков.

Умножение

Не знаю, затрагивалась ли раньше эта тема. Думаю, ни для кого не секрет, что в численных методах умножение матрицы на матрицу в чистом виде встречается крайне редко. По сути, результаты такого умножения - это те же попугаи. Поэтому было бы очень интересно увидеть результаты умножения матрицы на вектор, а также умножение матрицы на вектор в разреженном формате, например, CSR. Не знаю, позволяет ли CUDA работать с разреженными форматами. У Гото Sparce BLAS нет, в MKL есть.

"крайне редко" --- сильное преувеличение

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

не сказал бы

(предыдущий анонимус - это я)
Неточно выразился. Я имел в виду, что в реальных задачах такое умножение (заполненных матриц, да и вообще матриц) очень редко встречается.

По методу Гаусса, QR и прочая. А где это вы там умножение заполненных матриц нашли? Или будем умножение на матрицу поворота реализовывать по полной программе, как в тестах? http://blog.lexa.ru/2007/01/04/o_peremnozhenii_matric_i_prochix_arxitekt...

Когда решаем задачи мат физики, умножения матрицы на матрицу нет. Только умножение матрицы на вектор. Притом - разреженной матрицы, в том же формате CSR, например. Потому и сказал о попугаях.

На практике для реализации

На практике для реализации метода Гаусса (LU и Cholesky факторизации) используются блочные алгоритмы, смотрите, например, их реализацию в библиотеке LAPACK --- процедуры sgetrf и sportf. В блочных алгоритмах большинство операций совершается в умножении матриц и других, аналогичных процедурах из BLAS3. Блочные алгоритмы потребляют гораздо меньше пропускной способности памяти, поэтому работают значительно быстрее на современных машинах при больших размерах матриц.

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

new

Хорошо, разобрались, в Гауссе это нужно.

Тем не менее, очень интересен пример задачи матфизики (имеющий реальное применение!) с получающейся плотной матрицей - раз. С используемым умножением матрица*матрица - два. Слова "разреженные матрицы тоже встречаются очень часто" видятся мне не совсем корректными. Продолжаю настаивать, что при решении уравнений матфизики применительно к реальным задачам подавляющее большинство получающихся систем имеют разреженные матрицы.

Интегральные уравнения,

Интегральные уравнения, по-моему, приводят к плотным системам. Ещё метод Хартри-Фока, задачи рассеяния.

new

Копаю дальше, расширяю кругозор. А какие задачи приводят к интегральным уравнениям?

например, дифур по объёму

например, дифур по объёму иногда сводят к интегральному уравнению по поверхности. Если не ошибаюсь, в механике сплошных сред встречается. Ещё вроде в элекромагнетизме. Не моя область, плохо разбираюсь.

Так это классика. Только что

Так это классика. Только что там потом? При получении системы линейных уравнений будут использоваться базисные функции, которые, конечно, можно так подобрать, что будет заполненная матриц, но зачем? Берется функция-шапочка, все довольны.

это в методе конечных

это в методе конечных элементов так получается. Но там и интеграл берётся по всему объёму, и степеней свободы O(N^3). А если свести к интегралу по поверхности, то степеней свободы будет O(N^2), но и зависимость будет не локальной. Насколько я понимаю, в это уравнение будет входить свёртка с функцией Грина.

Я гораздо ближе знаком с аналогичной задачей в компьютерной графике --- уравнения рендеринга и радиосити. Есть источники света, есть непрозрачная поверхность, возможно неодносвязная, рассеивает свет по известному закону. Пространство однородное, свет не рассеивает и не поглощает. Нужно рассчитать освещённость каждой точки поверхности. Получается, что она зависит освещённости каждой другой видимой точки поверхности. Выходит интегральное уравнение для решения на поверхности с интегралами по поверхности. Как вариант, можно оперировать не с освещённостью в точке, а с распределением по полусфере входящего и исходящего излучения в точке. Причём, если среда неоднородная и рассеивающая, то возможно придётся решать уравнение переноса излучения по всему пространству, которое уже дифференциальное (точнее, интегро-дифференциальное).

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

Где можно почитать сведение

Где можно почитать сведение метода Хартри-Фока к интегральным уравнениям?

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

Почему сразу матфизика ? Вот

Почему сразу матфизика ?

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

химия химия

Почему матфизика? Да как обычно, у кого что болит...

Химия разная бывает. Горение с помощью МКЭ считается. Кстати, тот же Ansys МКЭ пользует. А МКЭ - это разреженные матрицы, притом умножения матрицы на матрицу при решении систем нет. А Ansys - это Ansys.

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

Я вот ещё почему об этом говорю: в третьей серии про умножении матриц была получена превосходная масштабируемость в Double precision (Intel MKL 9.0). А вот у меня на C2Q 6600, Intel MKL версии 10.0.3.020 и Intel Compiler версии 10.1.015 на разреженной матрице самый лучший результат - это увеличение скорости расчета в 1.4 раза на 4 ядрах по сравнению с 1. Ключ компиляции -O3. С ключем -O0 там ситуация заметно получше.

Я немножко не по этой части,

Я немножко не по этой части, поэтому стандартных тестов на разреженные матрицы не знаю, буду искать.

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

Про МКЭ - если правильно элементы нумеровать, то часто удается свести к матрицам, где ненулевые элементы вдоль диагонали, что жизнь сильно упрощает.

Да, я это и имел в виду,

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

Перенумерация элементов в МКЭ дело хорошее, но достаточно затратное. А если задача параболическая и сетка изменяется каждый шаг по времени?

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

В задачах, которыми я

В задачах, которыми я занимался (лет 10 тому как) генерация сетки была, конечно, затратной, но менее затратной чем все остальное.

Но занимался я 2D и осесимметричными задачами (течение смеси сред в пористой матрице, геология), на 3D тогда ресурсов не хватало. Сейчас - хватит, но уже не занимаюсь.

лучшее что я слышал --- 3-5

лучшее что я слышал --- 3-5 Гфлоп/c в SpMV на GPU. Недавно Университет Мериленда объявил конкурс на самую быструю реализацию, может там что лучше появится --- http://scriptroute.cs.umd.edu/gpucompete.

есть неплохая статья по изучению производительности SpMV на разных платформах, правда GPU там нет --- Williams et al. "Optimization of Sparse Matrix-Vector Multiplication on Emerging Multicore Platforms", SC07. Там тоже цитируется плохая масштабируемость по количеству ядер как на Intel так и на AMD процессорах. Дело в том, что производительность SpMV на больших матрицах ограничена пропускной способностью памяти. Поэтому количество задействованных ядер не столь важно --- они все на одном контроллере/шине сидят. Зато хорошо масштабировалось на двух процессорах. Лучше всего SpMV работало на Cell процессоре.

благодарю за ссылки

Спасибо за ссылки, пошукаю.

"...одинарной точности на GeForce 8800 получилось в 6 раз быстрее чем в двойной на CPU/Intel MKL..." - это некорректное сравнение, все-таки. Single и double сравнивать?

Насчет ПСП я в курсе. Инсайдерская, так сказать, информация. В MKL изначально оптимизация так сделана, что в ПСП всё упирается. Это накладывается на то, что у процов штеуда как раз там узкое место. Думаю, в тех тестах процы AMD заметно лучше масштабируемость показали.

если интересны какие-то

если интересны какие-то результаты умножения разреженная матрица - вектор на GPU, то можно посмотреть He et al. "Efficient Gather and Scatter Operations on Graphics Processors", SC07. У них в одинарной точности на GeForce 8800 получилось в 6 раз быстрее чем в двойной на CPU/Intel MKL. А в сопряжённых градиентах 1.6-1.8 Гфлоп/c на GPU.

результаты для ATI

Не плохо бы еще добавить результаты для ATI карт.
К сожалению, новых 48*0 карт у меня нет, но вот результат, полученный мною 3 года назад, для старенькой ATI1950XT, с использованием CTM (sgemm):
1024x1024 - 77.5 Gflops
2048x2048 - 80.4 Gflops

Интересно, сколько может выдать на гора ATI4870 ?

А разве три года назад был

А разве три года назад был CTM ?

Но это придирки, конечно. Давайте результаты свежих карт, опубликуем. У меня карт от АТИ нет

СTM родился в 2005. Только

СTM родился в 2005. Только всех кому его давали заставили подписать NDA. Дошли до 5 беты, а потом ATI всех кинула, заявив, что дальше будет только brook+ и об CTM можно забыть.
ххх...... Канадская компания ....ххх..
Я вот думаю, NDA уже кончилась или нет? Могу я уже продемонстрировать ATIшный алгоритм перемножения матриц ( весьма крутой) народу, например на cude, или у меня будут юр.последствия? Ребята с http://sourceforge.net/projects/amdctm вроде написали, что для CTM NDA закончилась. А то сижу на исходниках CTM BLAS как собака на сене, обидно.

Так можно же спросить

Так можно же спросить напрямую. Хотя, конечно, на такой прямой вопрос лично я бы ответил "нет, все продолжается".

Ну или можно завести анонима в ЖЖ, выложить все там и сделать вид что оно стало общеизвестным. (но это не я посоветовал , от всего откажусь, это хакеры!)

Потоковый BLAS это на самом деле жутко интересно, я бы почитал.

а что в нём крутого? Имеется

а что в нём крутого? Имеется в виду именно BLAS или просто умножение матриц? ATI публиковало детали реализации умножения матриц на CTM (тогда CTM назывался DVPM) в 2006. Потом, насколько я знаю, CTM не просто умер, а был развит и переименован в CAL, который сейчас свободно доступен для скачивания. CAL включает исходники умножения матриц и главу в программинг-гиде, расписывающую на три страницы как они работают.

Да крутого может и ничего,

Да крутого может и ничего, просто делать такие вещи эффективно без local store - не очень удобно.

1.крутой алгоритм

1.крутой алгоритм перемножения матриц для CTM нигде не публиковался. Публиковался стандартный алгоритм, демонстрирующий общую идею матричных операций на GPU. А это и так все знают.
2. имеется ввиду именно BLAS3
3. CTM и CAL - абсолютно разные вещи
4. в 2006 уже была бета. Называлась ,правильно, Data Parallel Virtual Machine.

странно, я что-то не видел

странно, я что-то не видел BLAS3 в CTM, может плохо глядел. Увы, так же не слышал про >120Gflop/s на X1900, которая была как раз актуальна в то время. Реализация, которая публиковалась (http://ati.amd.com/developer/siggraph06/dpvm_e.pdf) тоже была не совсем простой, она работала через fetch-4 расширение доступное только на картах ATI.

Intel MKL vs GotoBLAS

по моему опыту, SGEMM в Intel MKL (какая-то из 10-х версий) работает гораздо быстрее в 64-битном режиме. Разница примерно такая: 24 Gflop/s при запуске в 32-битных виндах и 38 Gflop/s в 64-битных на 2.67 GHz Core2 Duo E6700. То есть 64-битный Intel MKL тоже может работать под 90% пика, как и Гото. Наблюдал аналогичную картину и на четырёх ядерных процессорах.

Одно объяснений этому явлению что я слышал --- в 2 раза больший размер доступного файла регистров в 64-битном режиме по сравнению с 32-битным. Соответственно, можно использовать большие блоки в регистрах чтобы разрешить нехватку скорости кеша.

Все-таки 880 (980 - это

Все-таки 880 (980 - это теоретический расчет), но безумно круто в любом случае.

А между тем ATI Stream SDK

А между тем ATI Stream SDK 2.0 Beta2:
* First public beta release of ATI Stream SDK with OpenCL CPU support.
* The OpenCL 1.0 conformance logs from this release have been submitted to the Khronos OpenCL Working Group.

http://developer.amd.com/GPU/ATISTREAMSDKBETAPROGRAM/Pages/default.aspx

Надо, пожалуй, купить ATI-шную карточку, ибо есть тут меня идея-фикус. И на маке уже (вроде) все есть или с минуты на минуту будет....

Я думал, что Core2Quad

Я думал, что Core2Quad способен выполнять только 4 операции за такт. Где упоминается 8 операций? Можно ссылочку?

А это просто креативный

А это просто креативный подсчет, такой же как у производителей видеокарт.

SSE делает 4 32-битных MAD за такт. А MAD считается за две операции.

подскажите, пожалуйста, по

подскажите, пожалуйста, по какой формуле вы считаете гигафлопсы, для программы перемножения матриц на GPU? Допустим размер матрицы N*N, программа считает за m секунд.

Pages