LINUX.ORG.RU

[cuda] Решение систем ОДУ

 


0

2

Требуется решить численно систему N обыкновенных дифференциальных уравнений первого порядка (или одно уравнение N-го порядка, что, я так понял, сводится к первому случаю). Знаю, что есть много методов, в частности метод Рунге-Кутта.

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

Но проблема в том, что этот метод итерационный... для каждого последующего шага нам необходимо знать состояние предыдущего шага... то есть мы не можем выполнять шаги одновременно. А для эффективного ускорения это необходимо... Так как же это сделать? %) это вообще возможно?

Википедию и книжки по численным методам рыл - но ни одного параллельного метода решения диффуров не нашёл... они вообще есть?

думаю, в данном случае это параллелится достаточно просто. Перепиши систему в матричном виде, а потом уже параллель матричные операции. Вот полезная ссылка.

Rakot ★★
()

Думаю, что такие системы плохо параллелятся.

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

Понятно.. с матрицами проблем нет - всё довольно прозрачно. Осталось перевести систему в матрицы =)

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

По этой полезной ссылке, к сожалению, не нашёл, как решить систему ОДУ в матричном виде, используя матричные операции... может, ссылка не та? %)

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

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

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

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

про линейность Вы ничего в исходном посте не сказали.

Вообще то баян, есть 100 солверов для СЛАУ (к которым все сведется при неявной схеме) и 100500 библиотек линейной алгебры, к-е будут нужны при явных схемах (той же Рунге-Кутты). Из них наверняка есть скока то штук уже реализованных для GPU, берите и пользуйтесь.

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

И да, неявную схему не рекомендую.

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

стоп... СЛАУ - это системы линейных алгебраических уравнений? %) это, мягко говоря, совсем не то, что мне нужно... ну да... метод Крамера решения этих уравнений очень тривиальный. решается на раз два. и легко параллелится... %)

но мне нужна систему ОДУ... _дифференциальных) уравнений!

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

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

Возможно, я неточно сформулировал задачу.

Итак, есть система ОДУ (обыкновенных дифференциальных уравнений первого порядка). Например, N уравнений.

[code]dy0(x)/dx = f(x,y0,y1,...yN-1)

dy1(x)/dx = f(x,y0,y1,...yN-1)

...

dyN-1(x)/dx = f(x,y0,y1,...yN-1)[/code]

И дано, допустим, N краевых условий:

[code]y0(x0)=..., y1(x0)=..., ..., yN-1(x0)=.... [/code]

Короче, задача Коши. И её надо как-то решить... Вот.

Как это свести к СЛАУ? Возможно ли? Было бы здорово, конечно.

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

Решение системы ОДУ сводится к решению СЛАУ. Вот и всё. И да, Крамером СЛАУ никто не решает. Такие дела. Какой курс? 1 или 2?

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

Можно итерациями - раскладывая функции в ряд Тейлора.

Кстати, а функции заданы аналитически или численно?

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

Скорее всего численно... Тоже пока не знаю =) забываю о таких вещах сразу спросить.

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

Я никогда СЛАУ не решал, не приходилось, потому и назвал первое, что пришло в голову =) да, об этом на первом курсе рассказывали...

Сейчас на четвёртом как бы... если это как-то к делу относится.

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

Открываю «теплую ламповую» книжку Бахвалова. Численным методам решения ОДУ посвящена целая часть: 180 страниц.

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

Самый элементарный (но и неточный) - конечно-разностный метод, сводящийся к двум задачам, которые можно распараллелить: 1) вычисление коэффициентов конечных разностей (это можно делать отдельно для каждого уравнения - в параллели); 2) нахождение корней получившейся матрицы (cuBLAS в руки).

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

я читал какую-то книжку демидовича... численных методов решения ОДУ чуть более чем до хрена... например, Эйлера, Адамса, Рунге-Кутты, их я все знаю, более того, кодил их... приходилось.

Но нужен метод, который легко распараллелить ;) Рунге-Кутты, увы, не подходит...

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

Попробую... но мне такой метод тоже не нравится в силу не очень хорошей точности.

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

Ммм. вроде бы по предпоследней ссылке то, что нужно... Надо почитать =)

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

> Увы, вряд ли мне подойдёт что-то уже реализованное...

Почему??? Вы не первый в этой области, задача нетривиальная (я про качественную реализацию для GPU), но типовая. Вы посмотрите что уже сделано, потому как с Вашим уровнем это самому писать... неск. лет уйдет.

Мне для ВКР надо...

Что такое ВКР?

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

Заюзаете библиотеку, объясните как заюзали и как она работает.

провести анализ эффективности, ну и всё в том же духе.

Анализ эффективности проводится для уже работающего кода, требование «самостоятельной реализации» тут не при чем.

Вы определитесь для начала, сист ОДУ линейная или нет? И если хотите что б Вам помогли, все таки потрудитесь отвечать на вопросы: N ЧЕМУ РАВНО (ПОРЯДОК)?

И еще раз - СЛАУ появляется когда у Вас неявная схема. Для явной схемы (Рунге-Кутты той же) никаких СЛАУ и в помине нет, там все тупо, дали мяч - бери фигачь...параллелится по уравнениям системы на одном шаге.

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

> Но нужен метод, который легко распараллелить ;) Рунге-Кутты, увы, не подходит...

«Я тэбэ адын вещъ скажу, ты только нэ абижайся...»(c) ВСЕ ЯВНЫЕ ЧИСЛЕННЫЕ МЕТОДЫ ПАРАЛЛЕЛЯТСЯ ПРЕКРАСНО, ИМЕННО ПОТОМУ ЧТО ОНИ ЯВНЫЕ. Рунге-Кутты, Мунге-Путты, Вунге-Катты... все. Вот неявные методы в этом смысле куда хуже, но тут линал спасает, там свои финтифлюшки.

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

Что такое ВКР?

Дипломная работа.

параллелится по уравнениям системы на одном шаге

Синхронизацию придется либо через shm делать, либо через глобальную память. Т.е. уже немного тормоза будут появляться. Но все равно, сдается мне, быстрее будет распараллелить, несмотря на частое использование syncthreads.

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

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

Eddy_Em ☆☆☆☆☆
()
Ответ на: комментарий от registrant

Лучше по-английски гуглить. Материала намного больше получится.

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

> Не все параллелится: например, выделение изолиний

Я Вам приводил алгоритм построения изолиний со 100% внутренним параллелизмом.

А уш про явные методы для сист. лин. ОДУ и их параллелизацию... это ваще баян. Если N много больше числа ядер (а иначе о чем вообще разговор?) то накладные расходы на синхронизацию пренебрежимо малы. Про GPU не скажу, а так OpenMP в зубы (три строчки в код добавить) и все.

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

Я Вам приводил алгоритм построения изолиний со 100% внутренним параллелизмом.

Я в той простыне так и не разобрался. Shame me… Вот если бы Вы привели алгоритм…

явные методы для сист. лин. ОДУ и их параллелизацию...

Я этим пока не интересовался, но сдается мне, что в cuBLAS все нужное действительно уже есть.

OpenMP в зубы (три строчки в код добавить) и все

С GPU синхронизацию и прочее придется делать самому. Хотя, конечно, код для omp нуждается в минимальной доработке, чтобы стать кодом для GPU.

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

> Я в той простыне так и не разобрался. Shame me… Вот если бы Вы привели алгоритм…

Уф. Есть данные (чиселки) заданные в узлах(!) сетки.

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

2) Проходим контейнер(ы) по значениям изолиний и строим контуры (с-но изолинии) из отрезков. Понятно, что каждое значение обрабатывается отдельно и опять таки параллелится. Каждый контур - это список точек. Для каждого значения изолинии создаем контейнер с ключом - координатами на плоскости и значением - пойнтером на контур и флагом, определяющим куда точку добавлять (в конец или в начало списка). Фактически, мы отдельно храним концы контуров. Берем очередной отрезок, смотрим:

- если ни один из концов отрезка не соотв. концам контуров - это новый контур.

- если один из концов отрезка соотв концу контура - добавляем точку, меняем конец контура

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

- если оба конца отрезка соотв. концам одного контура - замыкаем контур, убираем его концы.

сдается мне, что в cuBLAS все нужное действительно уже есть.

Об чем и речь... наверное это первое, что вообще было реализованно;-)

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

Спасибо, положу себе в копилку - как соберусь, попробую реализовать. Правда, сложноватый алгоритм получается (с точки зрения сшивания контуров - в однопоточном-то режиме контуры сшивать не надо). Боюсь, что очень много перебирать придется - вычисления могут довольно сильно замедлиться.

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

Перебирать там ничего вообще не надо - для хранения концов контуров используется ассоциативный массив (я брал std::map). Такой алгоритм обеспечивает куда более высокую локальность данных, кроме того в трад алгоритме Вам все равно все ячейки походить...

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

> Что такое ВКР?

Выпускная квалификационная работа же, ну... Дипломная, проще говоря ;)

Почему??? Вы не первый в этой области, задача нетривиальная (я про качественную реализацию для GPU), но типовая. Вы посмотрите что уже сделано, потому как с Вашим уровнем это самому писать... неск. лет уйдет.

ok

«Я тэбэ адын вещъ скажу, ты только нэ абижайся...»(c) ВСЕ ЯВНЫЕ ЧИСЛЕННЫЕ МЕТОДЫ ПАРАЛЛЕЛЯТСЯ ПРЕКРАСНО, ИМЕННО ПОТОМУ ЧТО ОНИ ЯВНЫЕ. Рунге-Кутты, Мунге-Путты, Вунге-Катты... все. Вот неявные методы в этом смысле куда хуже, но тут линал спасает, там свои финтифлюшки.

В вот этого я что-то не догоняю %) как его параллелить, если каждый следующий шаг зависит от состояния, получаемого на предыдущем шаге? %) в том же методе рунге-кутты чтобы посчитать y(xn) надо знать y(xn-1)

Лучше по-английски гуглить. Материала намного больше получится.

Да, наверное с того и надо было начинать...

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

Про GPU не скажу, а так OpenMP в зубы (три строчки в код добавить) и все.

OpenMP мне не нужно. Мне нужно GPU. А это сотни и тысячи нитей... А какого порядка система уравнений - я пока не в курсе... надо уточнить. возможно, что меньше, чем количество нитей (потоков) ;)

BattleCoder ★★★★★
() автор топика

Тут недавно на конференции был. В общем китайцы на Тайвани написали DSMC на GPU. На 16 GF590 у них в некоторых тестах ускорение до 300 раз доходило по сравнению с 4-х ядерным зионом на 3-х ГГц. Так что тема интересная. Кстати поищи в инете по поводу МГУшного семинара по CUDA. На Физтехе тоже вроде этим занимаются.

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

Ну, у меня, к сожалению, нет 16 GD590, только одна GF9800 GT =) но хоть какое-нить ускорение по сравнению с обычным десктопным процем ведь будет, да?

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

> А уш про явные методы для сист. лин. ОДУ и их параллелизацию... это ваще баян.

Если верить формулировке выше (которая с ошибками), то система у товарища нелинейная. Как ты собрался параллелить?

Мне самому интересно. Если бы это было возможно, я бы сам, пожалуй, взялся за реализацию без всяких дипломных :) Но чувствую, что тут у вас ни хрена не получится. Производные то разные для каждого интеграла.

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

> В вот этого я что-то не догоняю %) как его параллелить, если каждый следующий шаг зависит от состояния, получаемого на предыдущем шаге? %) в том же методе рунге-кутты чтобы посчитать xn) надо знать y(xn-1)

Параллелится ОДИН шаг, в конце шага синхронизация. Поэтому ускорение может быть много меньше чем ожидается. Если в системе связаны только близжайшие уравнения (матрица скока-то-диагональная), то можно параллелить куда серьезнее, яндексите LRnLA. Но на GPU его реализовать.... не на диплом задача.

Ну, у меня, к сожалению, нет 16 GD590, только одна GF9800 GT =) но хоть какое-нить ускорение по сравнению с обычным десктопным процем ведь будет, да?

Не факт, заивисит от очень многих вещей. Мб и замедление...

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

А где на это можно посмотреть?

По-моему сам этап вычисления производных никак не распараллелить. Обновление после каждого шага еще куда ни шло, но производные?

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

Фиг знает. Потом синхронизировать это все. Может и сработает.

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

А какая разница? ну есть какие то правые части, их можно считать параллельно? А есть производные от правых частей, почему их нельзя считать параллельно? В чем разница то? И то и то некие алгебраические выражения, зависящие от вектора переменных известных с пред шага. Зависимости разрешены? Вперед.

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

Я еще не настолько знаком с CUDA. Как-то купил книжку. Попробовал почитать, но глубоко не влазил.

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

Я правильно представляю распараллеливание вычисления производных?

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