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