PARALLEL.RU

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

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

Одним из основных результатов разработанной технологии является численный метод решения разреженных СЛАУ по функциональным блокам. Основные черты метода опишем на примере моделирования некоторой системы, состоящей из набора связанных между собой функциональных блоков.

Итак, пусть имеется система, которую можно представить в виде набора функциональных блоков, соединённых друг с другом в N узлах (причём N<<R, где R общее количество узлов в системе). Поведение каждого из функциональных блоков описывается системой дифференциальных уравнений. При моделировании от систем дифференциальных уравнений методами численного интегрирования переходят к дискретным моделям (к СЛАУ) блоков. Далее формируют ?большую системную матрицу? (порядка R) путём включения в неё СЛАУ блоков. Причём если блоки функционально связаны, то они в матрице частично перекрываются (складываются). Полученную ?большую системную матрицу? решают прямыми или чаще итерационными методами. Это стандартный путь.

При моделировании системы по частям выполняются следующие шаги:
1. вместо формирования ?большой системной матрицы? формируется так называемая ?матрица сшивания?, ранг которой не превышает N. С помощью матрицы сшивания мы учитываем способ связи функциональных блоков между собой и их влияние друг на друга;
2. СЛАУ функциональных блоков решаются (любыми известными методами) параллельно друг другу относительно заданных значений векторов b в правой части;
3. с помощью решений, полученных на предыдущем шаге, формируем правую часть матрицы сшивания и решаем её. При формировании правой части матрицы сшивания из векторов решений СЛАУ блоков используется только небольшая часть отсчётов, т.е. процессоры на 2- 3 и 3-4 шагах обмениваются небольшим количеством данных;
4. далее опять параллельно друг другу решаем СЛАУ блоков, но с учётом решения, полученного на предыдущем шаге. На этом всё. Решения ?большой? СЛАУ для заданного значения общего вектора b* найдены. Если необходимо модифицируем вектора b блоков и ищем новые решения, перейдя к шагу 2.

Естественно результаты решения ?большой? СЛАУ стандартными методами и по функциональным блокам с учётом влияния ошибки округления совпадают. Как видно из описания, скорость решения СЛАУ по блокам на многопроцессорной системе зависит от максимального размера СЛАУ блока, размера матрицы сшивания и скорости обмена данными между процессорами, и практически не зависит от количества СЛАУ блоков и их способа соединения между собой, т.е. слабо зависит от размерности решаемой задачи.

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

Основные отличия и преимущества нового метода от известных состоят в следующем:
1. Он не требует формирования и решения целиком ?большой? СЛАУ;
2. Не фрагменты ?большой? СЛАУ, а СЛАУ функциональных блоков системы решаются отдельно и параллельно друг другу.
3. Если в системе нет явно выраженных функциональных блоков, например при использовании метода конечных элементов, то их размер можно задать самостоятельно, руководствуясь, например, вычислительной эффективностью.
4. При решении сверх больших задач возможна вложенная реализация метода (многоярусная реализация).
5. Матрица сшивания формируется путём включения в неё обратных матриц размерности менее N, которые на первом шаге рассчитываются отдельно и параллельно друг другу каждым из функциональных блоков. Следовательно, функциональный блок легко можно исключить или включить в исследуемую систему.
6. При наличии нелинейности в одном из функциональных блоков для нахождения решения можно, либо быстро откорректировать ?матрицу сшивания?, либо выполнить несколько итераций, воспользовавшись итерационной устойчивостью метода.

Надеюсь, что это будет кому-нибудь интересно.

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

Некоторые сведения о реализации "классической" диакоптики можно найти по адресу:
http://software.intel.com/ru-ru/article ... y-fortran/

Безусловно, всё есть. Поэтому решил написать маленькую библиотеку, в которой реализуется мой метод. Думаю, что к марту сделаю бета версию и вывешу её для тестирования в "облаке".

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

Yuri

С "классической" диакоптикой в России у меня прочно ассоциируются фамилии: д.т.н. Шакиров М.А., д.т.н. Курганов С.А., д.т.н. Филаретов В.В., их работы легко можно найти в интернете. Причём я полностью согласен с их принципиальными выводами, но не с методами реализации. Проблема этих работ заключается в том, что они сформулированы в терминах теории цепей, и поэтому воспринимаются как чисто прикладные узкоспециализированные ( т.е. математики их не читают и не воспринимают, а жаль). Во всяком случае, ни в одном из академических учебников по численным методам диакоптика не упоминается, хотя методы диакоптики - это фактически и есть прямые численные методы решения систем по частям.

В книге Г.Крона ?Исследование сложных систем по частям ? диакоптика?, можно найти
послесловие редактора А.В.Баранова (МГУ), где приводится ?перевод? диакоптики на математический язык блочных матриц. На этом математическом аппарате приводятся варианты метода расчета больших СЛАУ по частям. Возможно, это то, что Вам надо.
В настоящее время подобный подход развивается в рамках мультилинейной (тензорной) алгебры, материалы соответствующих конференций можно найти в Интернете.

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

Я тут посчитал, что если фрагменты системы имеют одинаковый порядок, и их количество равно N, и они моделируются по моему алгоритму на N процессорах, то в соответствии с законом Амдала (но без учёта накладных расходов на обмен данными) ускорение составит N/2, т.е. ускорение линейно зависит от N.

Поставил численный эксперимент на кластере. Предположил, что имеется система, состоящая из 22 функциональных блоков последовательно связанных друг с другом, причём последний связан с первым. Количество связей между блоками 80, порядок СЛАУ блоков 2000 (если сформировать ?большую? матрицу, то получается матрица порядка 42240) . Решил систему целиком (100с) и по частям (9.3 с), причём и там и там использовал функцию PARDISO из пакета MKL. Модуль вектора невязки составил в первом случае 1.3е-12, во втором - 4.5 е-12, что в принципе не плохо. В результате использования метода моделирования по частям, получил ускорение в 11 раз, что хорошо согласуется с результатом N/2.

Безусловно, можно возразить, что решая систему целиком, я использовал 1 процессор, а при решении по частям 22, поэтому выигрыш должен быть в любом случае. На это можно ответить, что если я могу эффективно распараллелить решение системы целиком с помощь общеизвестных алгоритмов, например, на 8 процессорах, то и при решении по частям, решая каждый из фрагментов системы тем же самым известным методом, мне ничто не мешает также использовать по 8 процессоров в каждом фрагменте. В результате при моделировании целиком мы сможем эффективно использовать 8 процессоров, а при моделировании по частям 22*8=176. Отсюда и выигрыш.

В послесловии А.В. Баранова приводится вариант расчета по частям системы, заданной отдельными блоками СЛАУ.
Взаимодействие подсистем выражается вспомогательными переменными. Приводятся условия эквивалентности
метода Крона и блочно-матричного метода (с.522-527).
Следует заметить, что в методе Г.Крона "большая системная матрица" не составляется. Рассчитываются
отдельные подсистемы и цепь пересечений, размер которой, как правило, не превосходит размера отдельной подсистемы.

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

Кроме того, я не утверждаю, что я первый изобрёл метод решения систем по частям. Как заметил Шакиров М.А., все методы диакоптики в чём-то схожи друг с другом. Просто, мой способ проще и эффективней других, как с алгоритмической, так и с вычислительной точки зрения. Как я могу это доказать? На данном этапе, поскольку я не хочу показывать математику, то только с помощью численных экспериментов с моей библиотекой, которую я планирую вывести на рынок.
.

Будем считать, что приведенных сведений достаточно, чтобы получить представление о Вашем методе.

Ещё раз проанализировал свой метод и пришёл к выводу, что действительно можно считать, что каждый функциональный блок системы вычисляет нечто напоминающее дополнение Шура, из которых потом формируется матрица сшивания. Это дополнение легко может быть вычислено после того, как была выполнена операция факторизации СЛАУ функционального блока.

Если блоки функционально связаны, то это значит, что их СЛАУ будут частично перекрываться (складываться) при формировании ?большой матрицы?. В численном примере, который я описал выше, 22 функциональных блока связаны друг с другом каскадно (последний связан с первым) количество связей 80. Это значит, что если сформировать "большую матрицу" то её портрет будет выглядеть следующим образом. 22 блока размещены на главной диагонали "большой матрицы" причём они перекрываются сверху и снизу друг с другом в 80 вершинах. Кроме того блоки 40 на 40 будут размещены в верхнем правом и нижнем левом углу матрицы (это обратная связь последнего функционального блока с первым). Хотелось бы отметить ещё раз, что мой метод позволяет решать систему с произвольным количеством функциональных блоков, которые связаны друг с другом любым способом (не только каскадным).
Надеюсь - этого будет достаточно, чтобы получить представление о методе.

Я буду участвовать в работе конференции "Параллельные вычислительные технологии" (ПаВТ'2011) Там же я планирую ознакомить желающих с демоверсией моей библиотеки. Кому интересна моя работа может встретиться со мной на конференции.

Сразу скажу, что совсем не силен в диакоптике (почти не знаком с ней). Но метод, о котором здесь идет речь, лично для меня весьма интересен, хотя он и не совсем нов.
Нечто подобное делается в Методе Конечных Элементов (МКЭ). Он очень широко распространен и, думаю, большинству известен, но, на всякий случай в вкратце о не скажу. Используется для решения сложных задач матфизики, т.е. позволяет решать системы диф. уравнений в ЧП в сложных областях при сложных начальных/граничных условиях.
Система представляется в виде совокупности конечных элементов (КЭ) малого размера. В каждом элементе есть некоторое конечное число узловых точек (чем их больше внутри КЭ, тем точнее результат). Значения неизвестных параметров в этих точках (узловые параметры) в совокупности во всех КЭ принимаются за глобальные неизвестные. Внутри же каждого элемента неизвестные находят по заданной аппроксимации
Далее система диф. ур-ий сводится к вариационной задаче, что уже снижает порядок исходных уравнений вдвое, т.е. находят такой набор глобальных неизвестных, чтобы некоторый функционал достигал стационарного значения (функционал может иметь либо физический смысл: потенциальная энергия, длина пути, время и т.д., либо создается искусственно по формуле Эйлера или, используя процедуры вариационных методов).
В процессе этого глобальная матрица системы в стандартном случае как раз формируется из блоков (которые располагаются на главной диагонали и идут как бы "внахлест" друг другу), каждая из которых - матрица каждого отдельного КЭ. Т.е. в процессе метода на разных этапах происходит то разбиение системы, то синтез.
Так вот как раз Ваша задача (т.е., как я понял, именно разбиение глобальной матрицы на несколько других меньшего порядка, их решение, а затем сшитие решений) решается в рамках т.н. Метода Суперэлементов - это развитие МКЭ, предназначенное для экономии вычислительных ресурсов. Система разбивается на подсистемы, те в свою очередь на КЭ (что весьма удобно - для разных подсистем можно задавать разную дискретизацию), для каждой из подсистем находят свою глобальную матрицу, находят решение каждой подсистемы, а затем склеивают. На практике это выглядит так. Например, если необходимо найти распределения поля скоростей, давления и т.д. воздуха вокруг самолета, то он разбивается на конструктивные составляющие: крыло, фюзеляж, хвост и т.д., для каждой такой составляющей независимо находят решение, используя МКЭ, потом все решения объединяют и получают глобальное распределение скоростей, давлений и т.д.
Но, насколько мне известно, последний метод не очень то оптимизировался под параллельные вычисления, так что я сомневаюсь в его эффективности, скажем, на кластерных системах. В этом "фишка" вашего алгоритма - он специализирован именно под параллельные вычисления, и это весьма интересно.
P.S. Надеюсь Вам какая-нибудь информация из вышеизложенного поможет. Успехов Вам на конференциях и конкурсах.
P.S.S. Уважаю Ваши коммерческие интересы, но все же очень бы хотелось, чтоб Вы хотя бы в двух словах (детали не столь интересны) познакомили нас с математической концепцией Вашего метода.

Во-первых, спасибо за комментарий. Я уже честно сказать думал, что этот форум мёртвый.

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

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

Если говорить о моём методе, то на вопрос о том, как Я разбиваю систему на части, я отвечаю по функциональным блокам, т.е. если функциональные блоки связаны, то их матрицы и вектора свободных членов частично перекрываются (складываются). В моём методе они не складываются, а решаются по отдельности, т.е. по функциональным блокам. Если говорить об устойчивости моего метода ? то он является прямым численным методом (не итерационным). Таким образом, его устойчивость определяется свойствами моделируемой системы (её обусловленностью и т.д.). Как видно из описания, количество частей, на которое надо разбивать систему, напрямую связано с количеством фрагментов, из которых она состоит, и их размера.

Что касается метода конечных суперэлементов (МКСЭ) Вы писали ?находят решение каждой подсистемы, а затем склеивают? весь вопрос как склеивают? Я думаю, с помощью моего метода это было бы гораздо быстрее.

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

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

Всё, тема закрыта. Предложенный метод является модификацией метода подструктур.

Работая над задачей hardware in the loop, т.е. когда одна часть системы D1 моделируется численно на компьютере, а другая часть D2 является реальным железом, пытался найти алгоритм сшивания. В результате нашёл матричный оператор S, причём в континуальной области, который по своим свойствам при правильном выборе параметров (коэффициенты при экспонентах) является идеальной линей задержки без потерь в системе D1-S-D2. При переходе в дискретную область оператор вырождается в два дополнения Шура, выполняя простые преобразования над которыми получаем модификацию метода подструктур. Оператор скорей всего ещё одна ипостасия оператора Пуанкаре-Стеклова, хотя это надо ещё доказать. Был бы благодарен за хорошую ссылку на литературу по свойствам операторов Пуанкаре-Стеклова (с работами Агошкова В.И. и Лебедева В.И. уже знаком).