Комментарии 17
А пробовали SageMath для символьных выражений?
Первым делом, я сделал так, чтобы повторяющиеся части в выражениях вычислялись один раз и хранились во временных переменных.
Это должен уметь делать компилятор. https://en.wikipedia.org/wiki/Common_subexpression_elimination
Далее я заметил, что часть подвыражений вообще не зависит от i,j,k, то есть является константами! Добавил автоматический вынос их за пределы цикла.
А это https://en.wikipedia.org/wiki/Loop-invariant_code_motion
Я заметил, что некоторые выражения повторяются, но со смещением на 1 по одному из индексов (индексу самого внутреннего цикла). В принципе, это ожидаемо – мы вычисляем поток через правую грань ячейки, потом сдвигаемся на 1 вправо и снова вычисляем то же значение для уже левой грани текущей ячейки. Почему бы не вычислить его один раз и не сохранить для следующей итерации цикла?
Не знаю, есть ли специальное название для такой оптимизации, но некоторые компиляторы это делать умеют. Правда тут требуются некоторые моменты учесть, дать компилятору определенные гарантии https://godbolt.org/z/Eo4dxz3x5
Компилятор должен понимать, что вызываемая в теле цикла функция является функционально чистой, и что входной массив не изменяется. Clang эту оптимизацию сделать не осиливает
Дело в том, что кусочки выражений, которые подставляются в зависимости от условий, в численных схемах не такие уж и большие. В принципе, дешевле будет вычислять сразу обе ветки ветвления, а потом код вида ...
Проблема тут в том, что компилятор может на эти "подсказки" наплевать и всё равно сгенерировать бранч. Вот например https://godbolt.org/z/xx5vf3jrK тут branchless получился только для функции test4. Инструкция bltu это "branch if less than unsigned"
Лучше было бы для таких вещей специальные директивы для компилятора придумать, чтобы пометить конкретную область кода или функцию, в которой бранчи не использовать. Может быть такое уже есть где-то.
Да, безусловно, многое может и сам компилятор. Но мне нужно было всё в комплексе -- распаралливание, векторизация, генерация кода для GPU и т.д. Я пробовал давать компилятору "сырые" выражения и он просто зависал при попытке их скомпилировать) Мой "компилятор" как раз знал всю подоплеку кода и мог делать эти оптимизации без проверок, что ускоряет и гарантирует.
Делать на уровне DSL то, что называется common subexpression elimination (CSE) чисто для того, чтобы компилятор быстрее компилировал - это еще имеет смысл, если куски подвыражений много раз повторяются и компилятор при оптимизации их долго сам находит и оптимизирует. А вот loop-invariant code motion это достаточно простая и понятная вещь, не думаю что компилятор будет сильно дольше компилировать (или хуже оптимизировать), если это не делать за него.
Чтобы максимально избежать простоев железа необходимо запускать счет на каждом ядре и GPU сразу же, как только для него готовы входные данные. К моменту, когда ядро или GPU завершило обработку, для него должен уже быть готов новый входной массив, чтобы следующий цикл счета запускался без задержек.
Почитайте про LRnLA алгоритмы, вроде это как раз из этой области
Спасибо! Можете дать ссылку? А то гуглится что-то не то, что нужно.
https://stackoverflow.com/questions/77846522/how-can-i-generalize-diamond-tiling-to-higher-dimensions вот тут например, в самом низу есть ссылки на ряд статей. Сразу скажу, я не специалист в этой области, так что напишу в рамках своего понимания. Насколько я понял, симулируются некие процессы в пространстве, пространство разбивается на некие кирпичики определенной формы, притом разбивается рекурсивно. Для симуляции некоей одной клетки нужно учитывать состояние соседних клеток вокруг нее, время в клетках вокруг допустим будет t, просимулировали процессы в клетке и получили ее состояние для времени t+n. И пространство можно каким-то хитрым образом разбивать на такие куски, чтобы простоя было меньше. Собственно, вот цитата из той ссылки:
The only hint I can find in the literature on how it may be done comes from a Russian research group at Keldysh Institute of Applied Mathematics. This group developed a series of different parallel spacetime decomposition algorithms in the last 20 years using a methodology called Locally-Recursive non-Locally Asynchronous (LRnLA) algorithms. LRnLa is not a single algorithm but a general guideline of how to design them. Basically: (1) The tessellation must should be a recursive process to utilize the entire memory hierarchy (it's not necessarily automatic, manual parameter tuning for different machines is allow as long as tiling algorithm itself can be generalized). (2) Parallelism should exist between different tiles, the dependency conflict problem should be solvable in some natural way. Using both requirements as starting point, researchers would manually examine the stencil dependencies and use their intuition in solid geometry to design custom algorithms to satisfy these goals. Unlike polyhedron compilers, these are custom solutions designed by human domain experts for human use, with geometric interpretations that ease implementations (but only from the perspective of mathematicians and physicists...).
А, иерархические сетки. Для задач газодинамики, где есть ударные волны, с ними нужно обращаться крайне осторожно. Одно неверное движение, индикатор не сработал вовремя, сетка не перестроилась и ударная волна пропала) Мы пробовали такие сетки, но для наших задач они оказались не очень хороши.
Но у вас ведь тоже на задачах газодинамики используется некий алгоритм разбиения пространства на некие куски, чтобы симуляцию этих кусков разбросать по отдельным вычислительным узлам. Каким образом вы решаете, на какие куски разделить пространство и как оптимальней его разбросать на разные ядра с учетом особенностей оперативной и кеш памяти? Этот LRnLA как раз об этом. Может тут еше подойдет полиэдральная модель оптимизации циклов https://xeno-by.livejournal.com/28122.html https://en.wikipedia.org/wiki/Frameworks_supporting_the_polyhedral_model
Может для вашей задачи получается, что некие области требуют более точного моделирования (с использованием более мелкой сетки) т.к. в них происходят более сложные процессы, и эту сетку можно динамически перестраивать, когда этот сложномоделируемый процесс (ударная волна или может турбулентность какая-то) будет смещаться в пространстве (т.е. придумать эвристику, которая отслеживала границы процесса и уплотняла под него сетку)? Применяется ли у вас что-то подобное?
Это как раз простая часть задачи. У нас все узлы одинаковы, так что просто делим сетку на куски равного объема (в числе ячеек). Там в статье есть картинка с разбиением. Конечно, так мы можем использовать не совершенно любое число узлов, а только такое, которое является произведением Npx, Npy и Npz, но это не проблема, так как на суперкомпьютере мы заказываем нужное нам число узлов для счета (и задача стоит в очереди, пока нужное число узлов не освободится).
Что касается использования более мелких сеток в отдельных частях -- это у нас решается более простым способом -- сетка остается топологически декартовой, но она может быть деформирована произвольным образом, сгущаясь к нужным местам. При этом у нас отсутствуют резкие переходы между сетками с разным разрешением, что большой плюс, так как на этих переходах всегда возникают артефакты.
У вас не возникала мысль попробовать Haskell вместо Lisp'a? Он, кажется,
доставляет больше удобств для работы с выражениями. Хотя я понимаю, что
порог вхождения там высок, а ваши основные задачи не в этом...
Есть еще такой проект https://herbie.uwplse.org/ для оптимизации вычислений с плавающей запятой. Там можно находить некие компромиссные оптимизации (например теряем в точности 5% для такого-то диапазона входных значений, но зато считаем в 5 раз быстрее). Написан на Racket (диалект лиспа)
О, спасибо, не знал про этот проект. Обычно такие места встречаются в численных схемах и авторы схем их вручную оптимизируют.
А до какого уровня авторы численных схем их оптимизируют, какие нюансы учитываются? Например, этот herbie "знает" о том, что в процессорах есть FMA инструкции (fused multiply-add) и может их использовать, в числовых схемах такое учитывают? Учитывается ли наличие SIMD инструкций (в подобных схемах могут явно указывать, что вот тут это уравнение должно считаться через SIMD?) ?
Сложный способ писать программы