г. СLANG- дает вполне компактный и эффективный код, но все равно может пропустить очевидные вещи (как с повторным divsd на одних и тех же данных)
div — это так, досадное недоразумение. Не в div-ах суть. Я даже не буду спрашивать, сколько div-ов на таком коде генерирует Дельфи :) Нормальный разработчик делает замену в коде множественного деления на умножение «на автопилоте».
А суть вот в чём: я переделал код на обработку массива матриц, и теперь векторизатор отработал лучше. Он считает по 2 матрицы за итерацию, за счёт этого все команды цикла векторные.
Причём можно сделать загрузку матриц в регистры эффективнее, если хранить их по-другому: не как массив структур, а как раздельные массивы для каждого элемента матрицы. Смотрите сами, насколько это допустимо в вашей ситуации.
Раздельные массивы хороши тем, что лучше масштабируются на любую систему команд, можно считать по 4 матрицы за итерацию с AVX.
И всё это в 20-30 строчках максимально простого кода, никаких ассемблерных портянок вручную.
периодически выходят скромные статьи с разбором профайлингов и сказками про то, как какие-то одинокие самоучки делают умножение матриц на CUDE лучше, чем отдел разработки NVidia. На хабре было несколько примеров
Мне кажется, конкретно операции с матрицами — очень ходовая и часто используемая штука, код стандартных библиотек должен быть вылизан до блестящего состояния и соревноваться с ними сложно.
Даже Ермолаев, при всей продвинутости его оптимизации, пишет в комментариях, что предположительно уступает MKL 5-10%.
Возможно, в вашем случае действительно будут влиять накладные расходы на вызов, у вас матрицы мелкие. В MKL для них предусмотрен какой-то хитрый инлайн, но при вызове из Дельфи он конечно работать не будет. Хотя можно написать промежуточную dll на Си, которая реализует функцию обработки массива матриц, с заинлайненной ф-ей MKL.
Мне также попадался одиночка, который заявляет, что у него быстрее, но для матриц общего вида — незначительно, те же 10% (transform inverse это, как я понимаю, упрощённый алгоритм для матриц трансформации в 3D-графике)
про читабельность на Си- я смотрел всякие портянки из интринсиков- честно- не вижу я там читабельности
Интринсики не особо читабельны, согласен.
Я обычно либо пытаюсь подтюнить код под автовекторизацию (это возможно, если понимать её логику и ограничения), либо использую векторные расширения Clang/GCC, они позволяют писать в «шейдерном» стиле. Ну знаете, в шейдерах/OpenCL/CUDA есть типы вроде float4, с которыми можно делать математику как с обычными float. Здесь тот же принцип.
А почему нельзя без сырцов?
Вы неявно используете массу dll из Винды, и для них тоже нет сырцов.
Предлагать перейти на компиляторы Си- это, простите, моветон
Ну да, дельфистам сразу слышатся отголоски холиваров «Дельфи vs С++».
Нет, я предлагаю менять не Дельфи, а ассемблер на Си.
А в случае с обращением матрицы- я в принципе не могу родить код (ни на С, ни на Fortran, ни на Pascal), который можно упаковать в одни регистры
В регистры AVX (или скажем AVX-512) влезает гораздо больше.
хочу. меня интересует скорость обращения массива из 1млн матриц 4*4*FP64 в 1 поток на 1 ядре
Поскольку это ваша задача, то давайте вы напишете тестовую программу на Дельфи, а я к ней прикручу сишный вызов. Вам виднее, какие матрицы должны быть, какие примеры данных. Это и как приложение к статье будет полезно.
Но главное- у меня только один divsd, а по Вашей ссылке- два дива
С делением нехорошо получилось, да. Но это же легко правится вручную, дописать одну строку D = 1 / D и поменять деление на умножение. Или переключиться на систему команд AVX (-mavx) или AVX2/FMA (-march=haswell), будет одно деление.
В целом я не исключаю, что может быть медленнее, в конце концов, я на эту «оптимизацию» потратил 2-3 минуты. А вы на свою сколько дней?
Даже если придётся векторизовать вручную, на Си это будет компактнее, читабельнее и более гибко, можно легко переключаться между 32/64 битами и относительно легко — между наборами команд.
Во-первых, операции с матрицами реализованы уже, наверное, тысячи раз и странно, что не нашлось ни одной подходящей библиотеки.
Я сначала подумал про Intel IPP, но пишут, что для операций с мелкими матрицами надо использовать MKL. Эти библиотеки сейчас бесплатные. Можно даже найти не очень старые заголовки для Дельфи (2019 года).
Во-вторых, если всё-таки хотите писать самостоятельно, рекомендую обратить внимание на компиляторы Си. Они могут генерировать эффективный ассемблер (в т.ч. с векторизацией) даже из самого примитивного кода безо всяких зачатков оптимизации. Вот например ваши формулы для расчёта обратной матрицы 3*3 обычным методом, просто «в лоб» скопированные в Си-код: gcc.godbolt.org/z/8c4b1G
Получилось с виду неплохо, около половины команд в ассемблере векторная (c окончанием на «pd»), по длине примерно так же, как у вас. Если хотите, можем сравнить реальную скорость.
Значит, в версии из статьи ничто не мешает увеличить кол-во вызовов? Я пробовал replace по всему вашему исходнику 1024 на 2048 или другое значение, получалось медленнее. Или так делать неправильно?
Ещё на будущее — инициализация OpenCL у вас сильно упрощённая, у меня установлены драйвера OpenCL для CPU, и GPU оказывается не в первой платформе, пришлось править.
Код
Const
MaxPlatforms = 4;
var
platform_id : array [0..MaxPlatforms-1] of CL_platform_id;
begin
// Получаем доступные платформы
ret := clGetPlatformIDs(MaxPlatforms, @platform_id, @ret_num_platforms);
if CL_SUCCESS<>ret then writeln('clGetPlatformIDs Error ',ret);
// Получаем доступные устройства
For i:=0 to ret_num_platforms-1 do begin
ret := clGetDeviceIDs(platform_id[i], CL_DEVICE_TYPE_GPU, 1, @device_id, @ret_num_devices);
If ret = 0 then Break;
end;
if CL_SUCCESS<>ret then writeln('clGetDeviceIDs Error ',ret);
И это ещё не учитывает системы с несколькими GPU (напр. ноутбуки с гибридной графикой).
Нет, CL_DEVICE_MAX_WORK_GROUP_SIZE это ограничение для размера local work group (у вас не задаётся), а для общего кол-ва вызовов ограничений нет.
Векторизация на современных GPU вроде бы не очень актуальна. Помню, что какое-то поколение AMD было на VLIW, там да.
Основная фишка примера NVidia — более эффективный доступ к памяти, что-то вручную кэшируется в локальную память, что-то переупорядочивается для более эффективного доступа.
Алгоритм в целом неудобен для GPU именно из-за произвольного доступа к памяти и малого кол-ва вычислений. Все эти безумные терафлопсы просто некуда приложить.
Если интересно, OpenCL на GF1070 всё-таки быстрее CPU, 0.097 секунды против 0.143.
Самый медленный этап это собственно сортировка (sort4), 0.016-0.013 * 4 прохода. Пересылка 0.010-0.013 * 2, туда и обратно.
В теории можно ещё ускорить, если увеличить кол-во вызовов, сейчас 1024, а у карты 1920 ядер.
Собрал для сравнения пример NVidia OpenCL Radix Sort, почему-то он не работает более чем с 4 млн. элементов — на 4 млн. получается 19 мс против 40 ваших с пересылкой, или 9 против 30 без пересылки.
У меня примерно так же.
Clang 10
RSort_Node1: 0.14313
RSort_Node2: 0.15709
RSort_Node3: 0.14532
GCC 9.3
RSort_Node1: 0.14448
RSort_Node2: 0.15864
RSort_Node3: 0.14716
если параллелить в данном месте, то в зависимости от входных данных нагрузка неравномерно распределиться по ядрам.
Для компенсации неравномерной нагрузки есть режимы с динамическим распределением потоков. В том же OpenMP для параллельного цикла это #pragma omp parallel for schedule(dynamic)
Итерации цикла могут занимать разное время, но их должно быть много, значительно больше кол-ва потоков/ядер.
Там же цикл не по корзинам, а по массиву данных. Я пробовал менять в Паскале порядок обхода на прямой — начинает сортировать в обратном порядке. Надо что-то править в инициализации корзин?
(пардон, в алгоритм не особо вникал — но вы ведь его и не описываете)
Паскаль
procedure Sort_step(offset: PIntArr; num: Integer; dest, source : PNodeArr; len : Longword);
var b : PByte;
i, n : Integer;
begin
// b:=@source[len]; // обратный
b:=@source[0]; // прямой
inc(b,num);
For i:=0 to len do begin
n := offset[b^]-1;
dest[n] := PNode(Integer(b)-num)^;
offset[b^] := n;
// Dec(b, SizeOf(TNode)); // обратный
Inc(b, SizeOf(TNode)); // прямой
end;
end;
Посмотрел внимательнее сишный вариант — вы не совсем правильно скопировали код picul, добавили фактически копирование структуры
//TNode const& src = source[i]; // < — здесь оперсанд опечатка?
TNode const src = source[i];
Я сам сишник-любитель и плохо знаю синтаксис ссылок (вот это &), но можно сделать обычный указатель:
TNode * src = &source[i];
unsigned int off = (src->key >> (sortable_bit * 8)) & 0xFF;
dest[offset[off]++] = *src;
тогда получаем:
RSort_Node1: 0.14094 seconds
RSort_Node2: 0.14432 seconds
разница уже близка к погрешности измерений
Ну так-то да — разные компиляторы, разный результат.
Clang, как и GCC, предпочитает вариант «задом наперёд»:
RSort_Node1: 0.14152 seconds
RSort_Node2: 0.15728 seconds
А Паскалевский вариант с прямым направлением есть?
Вообще порядок обхода элементов в RSort_step важен? Если нет, то можно в этом месте распараллелить.
Насколько я разглядел, на более ярких поверхностях размер белых точек в текстурах увеличивается, что в принципе правильно.
Фактически получился дизеринг по методу Halftone: en.wikipedia.org/wiki/Dither
Метод этот очень примитивный, и если бы он применялся к финальной картинке (а не «оклейкой стен»), то вообще ничего не было бы видно. Со стенами — приемлемо, только нужно специально объяснять, в чём фишка и зачем это нужно.
А в Obra Dinn сразу понятно, что это стилизация под монохромную/малоцветную графику, потому что полутона делаются именно таким дизерингом, как делались тогда.
(Вспомнился лайфхак, как быстро получить дизеринг картинки по методу Байера: это делает функция Windows/GDI StretchBlt при преобразовании из 8-битного битмапа в 1-битный с флагом HALFTONE. Во всяком случае, результат очень похож на картинку из Вики)
Идея из той статьи вполне применима, и вообще более «правильна» с точки зрения SIMD (делать одну операцию над несколькими данными, а не пытаться ускорять работу единственной операции).
Ну да, это удобно ещё и в том смысле, что логически векторный код не отличается от скалярного. В реализации ErmIg основная сортирующая функция вообще почти не меняется, благо интринсики спрятаны на более низком уровне.
Я добавил вариант того же «параллельного» алгоритма ErmIg для AVX2, и он всего в 1.36x раз медленнее.
Понятно. ErmIg получил почти такую же разницу между AVX2 и SSE2 — 1.33.
Похоже это общая закономерность SIMD-оптимизаций, наибольший прирост даёт SSE (хотя тоже не всегда пропорционально размеру вектора), а c AVX+ темп снижается.
А, пардон, до меня дошло, что у вас окно 7*1 и напрямую сравнивать с кодом по приведённой ссылке (там 3*3 для картинок) нельзя.
Но вопрос про железо остаётся в силе.
Спасибо за статью.
Не сравнивали AVX-512 с другими SIMD — SSE/AVX2? Вот например: habr.com/ru/post/204682
То есть можно догадаться, что в идеале AVX-512 должен быть в 2 раза быстрее AVX2, поскольку вектор в 2 раза шире. Но все эти сдвиги и перетасовки данных не добавляют ли нагрузки?
Также интересно, на каком железе вы тестируете. И будет ли оно работать IceLake-U — вроде как там поддерживаются не все подмножества AVX-512.
часть алгоритмов хорошо «ложится» на GPU, а многие совсем нет.
Согласен, но сложные алгоритмы часто появляются при попытках оптимизировать под CPU — а если закинуть ту же задачу на GPU, то может оказаться, что и брут-форс вывозит.
Отрисовка КТ на GPU, кстати, совсем не похожа на полигональную графику в играх, так что «для этого созданы» не совсем верно.
Я понимаю, что это две отдельные задачи, реконструкция и визуализация. Я хотел сказать (и подтвердить предположение Civil), что GPU здесь в целом подходят лучше, так как могут делать и реконструкцию, и визуализацию в реальном времени.
Сравнить скорость отрисовки на GPU и CPU можно например в этой софтине: ngavrilov.ru/invols/index.php?id=Download
Тогда непонятно, как можно хвалиться визуализацией, которая работает несколько секунд (второе видео в материале) вместо реального времени на GPU.
Для сравнения, 2015-й год и слабая ноутбучная видеокарта: www.youtube.com/watch?v=d9I9OoqHK2o
А суть вот в чём: я переделал код на обработку массива матриц, и теперь векторизатор отработал лучше. Он считает по 2 матрицы за итерацию, за счёт этого все команды цикла векторные.
Причём можно сделать загрузку матриц в регистры эффективнее, если хранить их по-другому: не как массив структур, а как раздельные массивы для каждого элемента матрицы. Смотрите сами, насколько это допустимо в вашей ситуации.
Раздельные массивы хороши тем, что лучше масштабируются на любую систему команд, можно считать по 4 матрицы за итерацию с AVX.
И всё это в 20-30 строчках максимально простого кода, никаких ассемблерных портянок вручную.
Можно и просто архивом выложить, без Гитхабов.
Даже Ермолаев, при всей продвинутости его оптимизации, пишет в комментариях, что предположительно уступает MKL 5-10%.
Возможно, в вашем случае действительно будут влиять накладные расходы на вызов, у вас матрицы мелкие. В MKL для них предусмотрен какой-то хитрый инлайн, но при вызове из Дельфи он конечно работать не будет. Хотя можно написать промежуточную dll на Си, которая реализует функцию обработки массива матриц, с заинлайненной ф-ей MKL.
Мне также попадался одиночка, который заявляет, что у него быстрее, но для матриц общего вида — незначительно, те же 10% (transform inverse это, как я понимаю, упрощённый алгоритм для матриц трансформации в 3D-графике)
Интринсики не особо читабельны, согласен.
Я обычно либо пытаюсь подтюнить код под автовекторизацию (это возможно, если понимать её логику и ограничения), либо использую векторные расширения Clang/GCC, они позволяют писать в «шейдерном» стиле. Ну знаете, в шейдерах/OpenCL/CUDA есть типы вроде float4, с которыми можно делать математику как с обычными float. Здесь тот же принцип.
Вы неявно используете массу dll из Винды, и для них тоже нет сырцов.
Ну да, дельфистам сразу слышатся отголоски холиваров «Дельфи vs С++».
Нет, я предлагаю менять не Дельфи, а ассемблер на Си.
В регистры AVX (или скажем AVX-512) влезает гораздо больше.
Поскольку это ваша задача, то давайте вы напишете тестовую программу на Дельфи, а я к ней прикручу сишный вызов. Вам виднее, какие матрицы должны быть, какие примеры данных. Это и как приложение к статье будет полезно.
С делением нехорошо получилось, да. Но это же легко правится вручную, дописать одну строку D = 1 / D и поменять деление на умножение. Или переключиться на систему команд AVX (-mavx) или AVX2/FMA (-march=haswell), будет одно деление.
В целом я не исключаю, что может быть медленнее, в конце концов, я на эту «оптимизацию» потратил 2-3 минуты. А вы на свою сколько дней?
Даже если придётся векторизовать вручную, на Си это будет компактнее, читабельнее и более гибко, можно легко переключаться между 32/64 битами и относительно легко — между наборами команд.
Я сначала подумал про Intel IPP, но пишут, что для операций с мелкими матрицами надо использовать MKL. Эти библиотеки сейчас бесплатные. Можно даже найти не очень старые заголовки для Дельфи (2019 года).
Во-вторых, если всё-таки хотите писать самостоятельно, рекомендую обратить внимание на компиляторы Си. Они могут генерировать эффективный ассемблер (в т.ч. с векторизацией) даже из самого примитивного кода безо всяких зачатков оптимизации. Вот например ваши формулы для расчёта обратной матрицы 3*3 обычным методом, просто «в лоб» скопированные в Си-код: gcc.godbolt.org/z/8c4b1G
Получилось с виду неплохо, около половины команд в ассемблере векторная (c окончанием на «pd»), по длине примерно так же, как у вас. Если хотите, можем сравнить реальную скорость.
Ещё на будущее — инициализация OpenCL у вас сильно упрощённая, у меня установлены драйвера OpenCL для CPU, и GPU оказывается не в первой платформе, пришлось править.
Векторизация на современных GPU вроде бы не очень актуальна. Помню, что какое-то поколение AMD было на VLIW, там да.
Основная фишка примера NVidia — более эффективный доступ к памяти, что-то вручную кэшируется в локальную память, что-то переупорядочивается для более эффективного доступа.
Алгоритм в целом неудобен для GPU именно из-за произвольного доступа к памяти и малого кол-ва вычислений. Все эти безумные терафлопсы просто некуда приложить.
Самый медленный этап это собственно сортировка (sort4), 0.016-0.013 * 4 прохода. Пересылка 0.010-0.013 * 2, туда и обратно.
В теории можно ещё ускорить, если увеличить кол-во вызовов, сейчас 1024, а у карты 1920 ядер.
Собрал для сравнения пример NVidia OpenCL Radix Sort, почему-то он не работает более чем с 4 млн. элементов — на 4 млн. получается 19 мс против 40 ваших с пересылкой, или 9 против 30 без пересылки.
Clang 10
RSort_Node1: 0.14313
RSort_Node2: 0.15709
RSort_Node3: 0.14532
GCC 9.3
RSort_Node1: 0.14448
RSort_Node2: 0.15864
RSort_Node3: 0.14716
Для компенсации неравномерной нагрузки есть режимы с динамическим распределением потоков. В том же OpenMP для параллельного цикла это #pragma omp parallel for schedule(dynamic)
Итерации цикла могут занимать разное время, но их должно быть много, значительно больше кол-ва потоков/ядер.
(пардон, в алгоритм не особо вникал — но вы ведь его и не описываете)
//TNode const& src = source[i]; // < — здесь оперсанд опечатка?
TNode const src = source[i];
Я сам сишник-любитель и плохо знаю синтаксис ссылок (вот это &), но можно сделать обычный указатель:
тогда получаем:
RSort_Node1: 0.14094 seconds
RSort_Node2: 0.14432 seconds
разница уже близка к погрешности измерений
Clang, как и GCC, предпочитает вариант «задом наперёд»:
RSort_Node1: 0.14152 seconds
RSort_Node2: 0.15728 seconds
А Паскалевский вариант с прямым направлением есть?
Вообще порядок обхода элементов в RSort_step важен? Если нет, то можно в этом месте распараллелить.
Фактически получился дизеринг по методу Halftone:
en.wikipedia.org/wiki/Dither
Метод этот очень примитивный, и если бы он применялся к финальной картинке (а не «оклейкой стен»), то вообще ничего не было бы видно. Со стенами — приемлемо, только нужно специально объяснять, в чём фишка и зачем это нужно.
А в Obra Dinn сразу понятно, что это стилизация под монохромную/малоцветную графику, потому что полутона делаются именно таким дизерингом, как делались тогда.
(Вспомнился лайфхак, как быстро получить дизеринг картинки по методу Байера: это делает функция Windows/GDI StretchBlt при преобразовании из 8-битного битмапа в 1-битный с флагом HALFTONE. Во всяком случае, результат очень похож на картинку из Вики)
Понятно. ErmIg получил почти такую же разницу между AVX2 и SSE2 — 1.33.
Похоже это общая закономерность SIMD-оптимизаций, наибольший прирост даёт SSE (хотя тоже не всегда пропорционально размеру вектора), а c AVX+ темп снижается.
Но вопрос про железо остаётся в силе.
Не сравнивали AVX-512 с другими SIMD — SSE/AVX2? Вот например:
habr.com/ru/post/204682
То есть можно догадаться, что в идеале AVX-512 должен быть в 2 раза быстрее AVX2, поскольку вектор в 2 раза шире. Но все эти сдвиги и перетасовки данных не добавляют ли нагрузки?
Также интересно, на каком железе вы тестируете. И будет ли оно работать IceLake-U — вроде как там поддерживаются не все подмножества AVX-512.
Отрисовка КТ на GPU, кстати, совсем не похожа на полигональную графику в играх, так что «для этого созданы» не совсем верно.
Сравнить скорость отрисовки на GPU и CPU можно например в этой софтине:
ngavrilov.ru/invols/index.php?id=Download
Тогда непонятно, как можно хвалиться визуализацией, которая работает несколько секунд (второе видео в материале) вместо реального времени на GPU.
Для сравнения, 2015-й год и слабая ноутбучная видеокарта:
www.youtube.com/watch?v=d9I9OoqHK2o