Большая разреженная матрица

[roof]

Доброго времени суток! Я пишу программу для вычислений методом конечных элементов, сейчас возникла задача решения оооочень большой разреженной системы линейных уравнений. Если кто-то сталкивался с похожей задачей, помогите найти пример реализации на C++ матрицы с учетом ее разреженности, и еще подскажите, стоит здесь возиться с такой штукой как valarray? У меня есть работающая модель на фортране, но я его знаю плохо, да и хотелось бы использовать все доступные в языке средства оптимизации.

[Serguie]

Да, когда то я тоже занимался МКЭ! Дело в том что, нет смысла хранить всю разреженную матрицу целиком впамяти (т.е. нолики хранить не надо), есть несколько способов хранения разреженных матриц: разреженно-строчный формат, диагональный формат и их вариации. Причем сборка матрицы в программе должна вестить с учетом этих форматов, соответственно и решение СЛАУ тоже. Решение СЛАУ так же можно сделать несколькими способами.
Помню, на 3-4 курсах, я делал курсовую работу по МКЭ, у меня есть отчет и сами исходники на C++, если надо - могу выслать. Куда?

[roof]

Буду просто непредставимо благодарен!
Убрал адрес от греха подальше

[Serguie]

ok. сегодня поисщу, завтра с утра отправлю

[roof]

Сенк ю большое, докачал, буду разбираться.

[VLarin]

Я тоже писал программу расчета МКЭ на С++, пробовал разные библиотеки хранения разреженных матриц и решатели, пришел к тому, что:
1. Для хранения sparse-матриц хорошо подходит MTL (Matrix Template Library) (http://www.osl.iu.edu/research/mtl/).
2. У них же есть библиотека итерационных решателей ITL(http://www.osl.iu.edu/research/itl/) - там есть CG и набор preconditioner'ов, включая несколько модификаций предобуславливания неполным разложением Холесского, что вместе с CG дает используемый комерческими программами ICPCG с достаточно хорошим быстродействием.

Кстати, если не секрет, можете вкратце описать, что Вы пишите и какие элементы используете?

[roof]

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

[VLarin]

Вся ценность открытых библиотек - в просмотре кода. В MTL, на мой взгляд, есть достаточно грамотные реализации sparse-матриц, а в ITL - итерационные решатели. При этом с помощью Ctrl+C-Ctrl+V можно на первое время и свою библиотеку сотворить, а потом ее со временем усовершенствовать и переработать в пух и прах. А так, главное в MTL и ITL - это использование шаблонов, и то как их создатели все сделали - супер, просто и ясно.

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

[VLarin]

А, кстати, с генератором сетки к чему пришли?

[roof]

Ни к чему еще не пришли. Все еще только начинается. Читаю чужие исходники и умные книжки, поэтому за ссылки огромное СПАСИБО! Вы про сетки просто так спрашиваете или хотите что-то конкретное обсудить? У меня пока что вопросов гораздо больше чем ответов.

[Serguie]

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

[VLarin]

To roof:
Про сетки просто так спрашиваю. Интересно было, что Вы с ними решили - писать самим или брать что-то готовое. Я просто в свое время писал свою реализацию Advancing Front в 2D. Тут хорошо помогла книга Hand Book of Grig Generation (там глава 14 - структура данных, глава 17 - Advancing Front). Сам по себе Advancing Front медленный, но благодаря использованию деревьев для поиска, таких как Alternative Digital Tree (ADT) (в книге про них хорошо написано), получается очень даже ничего, достаточно шустро.
Делоне для МКЭ - это, по-моему, для МКЭ полная ерунда. Advancing Front позволяет генерить сразу качественные сетки, а самое главное позволяет вводить функции изменения шага сетки в расчетной области, то есть просто задать набор математических функций, которые будут определять изменение сетки во всей области, то есть максимально-допустимые размеры элементов в произвольных точках. И потом, если использовать adaptive h-refinement, то можно на грубой сетке сразу определить характер требуемой сетки - то есть где укрупнить, где уменьшить, а потом задать все это как функции для генерации и сгенерить нужную сетку.

А так генерация сетки Advancing Front'ом сводится к определению функций, с помощью которых сетка сгустится в местах сближения сегментов границ, плюс еще можно задать коэффифицент роста сетки - т.е. как будет расти размеры КЭ в процесе генерации для создания сеток с большим градиентом размеров КЭ. Тогда сетка сетка получится плотная в области сближения границ (т.е., например задать коэффициент плотности сетки, так чтобы в промежутке между двумя любыми объектами было не менее определенного числа элементов), а вдали от этих границ сетка вырастет и, там где не нужны мелкие элементы, они там сами получатся крупными благодаря их росту в соответсвии с коэффициентом роста.

Найти Hand Book of Grig Generation можно в любом файлообменнике, нужно только поискать. Я его в свое время в DC++ стянул.

А Делоне - в принципе его тоже можно до ума довести, но Advancing Front рулит в плане создания качественных сеток с первого раза в 2D.

[VLarin]

Кстати, если у Вас расчетное ядро еще в процессе разработки, советую посмотреть иерархию классов объектно-ориентированных FEM-программ: OOFEM, еще есть в инете диссер McKenna, посвященный объектно-оринетированныму программированию МКЭ. В диссере рассмотрены структура классов, заточенная под использование различных элементов (включая супер-элементы) + объектно-ориентированный Domain Decomposition.

И в продолжение темы по решателям и матрицам - на сайте OOFEM есть ссылки на IML++ (template library for numerical iterative methods) и SPOOLES (sparse direct solver). На сайте IML++, кроме кода, есть ссылка на книгу
Templates: Building Blocks for the Solution of Iterative Systems.

[roof]

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

[VLarin]

В общем, с расчетным ядром советую не торопиться, посмотреть другие реализации (тот же OOFEM + диссер McKenna), я просто тоже начал в лоб писать, в итоге плохо спроектировал иерархию классов и за год четыре раза расчетное ядро переколбашивал заново с учетом новых требований, которые в начале и не задумывал. Главное грамотно спроектировать ядро, а то потом все поползет. А самая тема - это сразу сделать заточку под супер-элементы и Domain Decomposition (или как его в ELCUT обозвали - метод геометрической декомпозиции), там просто и генератор сетки нужно подстраивать под это дело, и решатель, чтоб потом на больших задачах все летало.