LINUX.ORG.RU

[openmp]помогите составить алгоритм

 


0

1

хочу распараллелить последовательную подпрограмму с помощью openmp. подпрограмма представляет собой обход дефрагментированого разреженного массива F(N1,N2), где N1~500, N2~100000. в служебном массиве V(N1)<N2 хранится значение последнего значимого элемента подмассива F(i,:).

значимые элементы содержат в себе фазовые вектора частиц p=(x,v,i) где x — координата, v — скорость, i — номер узла, причём в F(i,:) хранятся лишь те вектора, чей (x) лежит в заданном интервале (т.е. все вектора локализованы)

в процессе счёта значения x меняются так, что частица может или перескочить в соседний узел или остаться в том же.

сейчас я просто обхожу последовательно все узлы i=0..N1, все вектора j=0..V(i), провожу вычисления, и в случае, если x не соотвествует узлу i, то переношу его значение в соседний узел i2, а на его место записываю нижнее значение из F(i,V(i)) (и делаю V(i)-1, V(i2)+1).

как организовать такой обход массива при параллельном расчёте? пока на ум приходит только завести массив такого же ранга F2, все сложные вычисления проводить параллельно в F, а затем обойти F последовательно, копируя значения в уже в нужные узлы F2, но это явно непотимально. какие есть ещё варианты? с учётом того, что порядок в котором вектора записаны внутри F(i,:) не имеет значения. ссылки на статьи приветствуются.

★★★★★
Ответ на: комментарий от VLev

Нет и не будет.

жаль.

От сетки как минимум требуется структурированность.

в регулярной сетки 2D внутренний 1 узел сетки = 4 ячейки. Для нерегулярной сетки 1 узел = N-ячеек. LRnLA генерирует иерархию наборов данных, на одном уровне шаблон наборов данных один. нельзя ли для нерегулярной сетки генерировать иерахию уровней, которые содержат различные шаблоны наборов данных?

marat2345
()
Ответ на: комментарий от marat2345

нельзя ли для нерегулярной сетки генерировать иерахию уровней, которые содержат различные шаблоны наборов данных?

можно, но проблема неструктурированных сеток не в этом, а в отсутствии структуры (КО).
Потому я и не хочу ими заниматься.
Что касается LRnLA, то суть в том, что в настоящем виде эти алгоритмы работают с «сетками» вообще одной и той же универсальной структуры (т.н. обобщенный граф зависимостей LRnLA). Причем количество интересных мне задач, которые описываются этим самым графом столь велико, что я не уверен, что до каких-то принципиально других структур у меня дойдут когда-нибудь руки.

С другой стороны, всегда есть возможность согласования любимой Вами сетки с обобщенным графом зависимостей LRnLA.
Например, неструктурированную сетку можно взвесить на структурированной с ячейками бОльшего объема, и свести тем самым задачу к предыдущей.

VLev
()
Ответ на: комментарий от AIv

Если поставленная задача на 100% ясна программисту с самого верха до самого низа её реализации, то конечно, вручную он может реализовать её эффективней всего. Но столь всеобъемлющее понимание программистом поставленной задачи — это, имхо, фантастика. IRL, как известно, имеет место специализация. Иначе не было всей этой туевой хучи ЯП.

blexey ★★★★★
()
Ответ на: комментарий от blexey

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

blexey ★★★★★
()
Ответ на: комментарий от blexey

Если поставленная задача на 100% ясна программисту с самого верха до самого низа её реализации,

В числ моделировании надо знать всю цепочку: физ-мат модель->числ схема->код. Иначе просто не сможешь отладить программу. Ес-но используются всякие библиотеки, но узкая специализация и хорошая оптимизация - вещи слабосовместимые;-) Оптимизация подразумевает понимание неск уровней.

AIv ★★★★★
()
Ответ на: комментарий от qnikst

а почему #pragma omp не перед первым циклом?

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

VLev
()
Ответ на: комментарий от AIv

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

Rakot ★★
()
Ответ на: комментарий от VLev

Это я под анонимом писал.
У меня для каждой ячейки есть ссылка на первую частицу динамического списка частиц данной ячейке. Тип частица содержит координаты частицы, компоненты скорости (трёхмерные векторы), целое число для указания компоненты газа(плазмы), и указатель на следующее число в списке. Таким образом, у каждой ячейки свой список. При движении, соответственно, указатели перебрасываются от одной ячейки к другой.

Rakot ★★
()
Ответ на: комментарий от VLev

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

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

#pragma omp parrallel for
for (i=0;i<N;i+=N/cores) {
   for (j=i; j<i+N/cores;j++) {
      //actions
   }
}
qnikst ★★★★★
()
Ответ на: комментарий от Rakot

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

Вот о чем и речь...

LRnLA на неструктурированные сетки (тераэдроны) не перенесли?

[openmp]помогите составить алгоритм (комментарий)

что бы вы посоветовали для этого?

LRnLA на прямоугольной сетке;-) М.б. Вадим озвучит тут параметры нового 3D3V PIC-кода... хотя было два доклада на последней звенигородской конференции как раз по этому коду.

AIv ★★★★★
()
Ответ на: комментарий от Rakot

При движении, соответственно, указатели перебрасываются от одной ячейки к другой.

угу. это как раз и есть главный источник проблем с эффективностью в методе PIC.
В начале расчёта, когда все частицы упорядочены, доступ в память последовательный. Такой тип доступа ограничивается темпом (~>10GB/сек)
А если в процессе расчёта частицы полностью перемешаются --- станет случайным. А такой тип ограничивается латентностью ~>40нсек/частицу.
В результате замедление может составить сотни раз.

Кстати что бы вы посоветовали для этого?

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

VLev
()
Ответ на: комментарий от qnikst

Не обращая внимание на ошибки в применении omp...
Ваше решение не решает основной проблемы «распаралелливания» таких задач --- гонки данных (http://www.viva64.com/ru/t/0042/).
Он не дает никаких гарантий, что два потока не попытаются изменять данные «граничных» ячеек i=(N/cores)*icore «одновременно»

Кстати, для Rakot:
несмотря на весь примитивизм решения организации хранения частиц в массиве, описанный автором треда --- он будет весьма эффективным, особенно если второй его предел (N2=100000) сделать простым числом.

VLev
()
Ответ на: комментарий от VLev

Не обращая внимание на ошибки в применении omp...

?!

Ваше решение не решает основной проблемы «распаралелливания» таких задач --- гонки данных (http://www.viva64.com/ru/t/0042/).

спасибо за ссылку на 2+2=4, если уж и учить необразованных, то надо уж на что-нить адекватное давать. спасибо что не на википедию

Он не дает никаких гарантий, что два потока не попытаются изменять данные «граничных» ячеек i=(N/cores)*icore «одновременно»

оно даёт ровно такие же гарантии, как и приведённое Вами. Разница в том, что «домены» в массиве в Вашем случае будут: xovxovxov..., моём xxxxooovvvvkkkk (где кол-во доменов должно быть равно кол-ву процессоров^W процессов используемых openmp). В вашем случае по реальным процессам данные разойдутся как-то так (в случае static) 1??1??1??2??2??2??... в моём 1111222233334444..... В обоих случаях наличие/отсутсвие гонок определятся тем действием, которое производится над массивом, я _честно_ не увидел в описании ТС-а того, как именно производится подсчёт, поэтому я не хочу _угадывать_ где будут гонки, а где нет.

...особенно если второй его предел (N2=100000) сделать простым числом.

а можно ссылочку на работу показывающую, что простые числа (не кратные размерам страниц памяти/кеша и т.п.) используются где-то кроме криптографии?

P.S. я лично не до конца понял задачу (по описанию ТС), поэтому не уверен, что решал бы её правильно

P.P.S. перечитал задачу и понял почему было желание локализовать по 3 элемента, если есть хоть какая-нибудь гарантия, что одна и та же частица не пересчитается 2жды за один проход

qnikst ★★★★★
()
Ответ на: комментарий от qnikst

плохое лекарство от гонки данных в OPENMP директива !DEC$OMP CRITICAL (для Фортрана)

marat2345
()
Ответ на: комментарий от qnikst

лучше сделать декомпозицию данных по потокам, найти общие узлы и для них сделать отдельный код с OMP CRITICAL. либо на каждый общий узел зависти массив размерностью NUM_THREADS, каждый поток в OPENMP будет писать в свой элемент вспомогательного массива. После OPENMP секции reduce вспомогательного массива

marat2345
()
Ответ на: комментарий от qnikst

Разница в том, что «домены» в массиве в Вашем случае будут: xovxovxov..., моём xxxxooovvvvkkkk

а почему не xxxx........oooo..........vvvv....... ?

marat2345
()
Ответ на: комментарий от marat2345

мне кажется, что OMP_CRITICAL можно заменить атомарными операциями (хотя тут массивы и может не прокатить). собственно когда приходилось помогать оптимизировать данные для расчёта какой-то страшной биофизики, я предложил именно такой же способ :) но там считает жирный кластер и ограничений на оперативку (и кеши, чтобы команда VLev и Alv меня не съела) не было.

qnikst ★★★★★
()
Ответ на: комментарий от qnikst

Тут вот какая петрушка получается

Честно говоря, я не собирался много времени тратить на обсуждение этой задачи --- сама по себе она тривиальная и может быть решена миллионом разных способов, а своё решение дал лишь потому, что оно показалось мне достаточно элегантным и очень в стиле «дешёвого параллелизма» --- в последовательный код добавляется всего 2 строчки, в третьей меняется один символ --- и все.

Но, такое впечатление, что никто кроме меня этой элегантности не почувствовал... :(
В любом случае, нижеидущее объяснение уже никак нельзя назвать «элегантным».

я _честно_ не увидел в описании ТС-а того, как именно производится подсчёт, поэтому я не хочу _угадывать_ где будут гонки, а где нет.

thunar> в процессе счёта значения x меняются так, что частица может или перескочить в соседний узел или остаться в том же.

соседей у «узлов» два --- левый и правый. То есть узел с номером i на шаге k может изменить данные только узлов i-1 и i+1.
Для предотвращения гонки данных надо запретить в этот момент производить вычисления в этих узлах, а также в их соседях --- i+2 и i-2. А вот вычисления в i+3 и i-3 уже можно производить параллельно (точнее, асинхронно).
В результате выполнения этих асинхронных действий вычисления будут произведены в каждом третьем узле, конца чего надо дождаться (это барьерная синхронизация), после чего надо сдвинуться на единичку по i и сделать то же самое, а потом ещё раз.

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

VLev
()
Ответ на: комментарий от marat2345

SoA лучше AoS

эээ, откуда такое умозаключение?
Я нас, кстати, все наоборот. Данные LRnLA организованы в стиле массива структур.

VLev
()
Ответ на: комментарий от VLev

я рад за Вас. думаю лучше чтобы каждый остался при своем мнении. для Вас хорошо AoS, для меня SoA.

marat2345
()
Ответ на: комментарий от qnikst

по моему скромному опыту работы с OMP --- он сам разберётся как наиболее эффективно разбить цикл на домены.

VLev
()
Ответ на: Тут вот какая петрушка получается от VLev

при четвёртом перечитывании задачи я уже понял откуда было разбиение по 3, что я отметил в PPS. И данное решение действительно будет очень хорошим при условиях: того в формуле пересчёта не используются данные из соседних столбцов, ничего страшного, если мы одну частицу пересчитаем 3 раза за 1 проход.

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

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

qnikst ★★★★★
()
Ответ на: комментарий от qnikst

последовательный алгоритм, вобще формально позволяет частице за 1 проход перейти в последний столбец,

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

	type p1d3v !тип для электронов
	real(8)::x
	real(8),dimension(3)::v !v1=vx; v2=vy; v3=vz
	real(8)::ev
	logical::s !отмечает обработанную/вторичную частицу
	end type p1d3v
thunar ★★★★★
() автор топика
Ответ на: комментарий от VLev

Спасибо за замечания.

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

Rakot ★★
()

Чемто знакомым потянуло

N1~500, N2~100000

N1 - аминокислота, N2 - длина белка? Folding, не?

anonymous
()
Вы не можете добавлять комментарии в эту тему. Тема перемещена в архив.