Подскажите с методом оптимизации матрицы жесткости

[vitaly333]

Пишу математическое ядро для программы расчета линейных стержневых систем(фермы,балки,рамы) методом конечных элементов.
Написал несколько прямых методов для решение СЛАУ. Но они без оптимизации. Хотелось бы еще написать оптимизацию матрицы жесткости т.к она может быть просто огромна.
Почитал тут на вашем форуме про различные методы и решил что буду исходную матрицу превращать в ленточную или профильную а потом с помощью обратного алгоритма Катхилла-Макки уменьшать ширину ленты+буду работать только с одной треугольной частью матрицы жесткости.
Почитал книжку А.Джордж , Дж.Лю - Численное решение больших разряженных систем - но там ни хрена не понял как это делать.Можете кто-нибудь привести конкретный пример как преобразовать обычную матрицу в ленточную (или профильную) а потом уменьшить ширину ленты с пояснениями пожалуста а то я только начинаю работать в этом направлении еще много чего не знаю.
И ещё меня тревожат такие вопросы:
1.На каком этапе выполнять эту оптимизацию матрицы?
На этапе сборки глобальной матрицы?На этапе сборки локальных матриц? После сборки глобальной матрицы?
Если это делать после сборки глобальной матрицы то под глобальную ведь нужно выделять память - а это мо-моему уже не рационально, т.к как если например в системе 10тыс. узлов *3 степеней свободы =Глобальная матрица размером 30000X30000 , 900 000 000*4байт = 3600 000 000байт = 3,4Гб (Если использовать симметрию , то 1,7Гб) - это ведь памяти нехватит!!!
2.Слышал что есть метод Слоана , который получше чем Катхилла-Маки, ну чё то про него нигде ничего не нашел.

[VLarin]

Reverse Cuthill-McKee (вроде так пишется:)) предназначен для перенумерации узлов сетки так, чтобы потом ширина лента была меньше исходной и работает он с графом, получаемом как из сетки, так уже и из глобальной матрицы. Таким образом, им можно пройтись по сетке КЭ сразу после ее создания и перенумеровать узлы один-единственный раз, а потом считать, считать и считать...

Но это что касается ленты и применения заточенных под это решателей.
А почему бы просто не воспользоваться готовыми (или написать по их подобию своё методом Ctrl+C + Ctrl+V ) ) библиотеками хранения и решения спарс-матриц (итерационными методами, сопряженным градиентом в частности). Тогда все Ваши вопросы просто отпадают, для них перенумерация и прочее не требуется вовсе.

И потом 4 байта - это float, а по-хорошему нужен все-таки double, для резконеоднородных полей с большими градиентами решения.

Кроме того, ленточная форма хранения не такая уж и компактная, так как в строке всего 10-20 ненулевых элементов, а ширина лента может быть и более. Да и симметрия сильно помогает в случае спарс-матриц.

В этом топике уже был разговор и числе ненулевых элементов в матрице жесткости и ее компактном хранении.

[vitaly333]

Цитата:

А почему бы просто не воспользоваться готовыми (или написать по их подобию своё методом Ctrl+C + Ctrl+V ) ) библиотеками хранения и решения спарс-матриц (итерационными методами, сопряженным градиентом в частности). Тогда все Ваши вопросы просто отпадают, для них перенумерация и прочее не требуется вовсе.

Всё таки хотелось самому написать.

[VLarin]

Тогда Вам нужно написать:
1. Библиотеку работы со спарс-матрицами (хранение + базовые операции с матрицами)
2. Библиотеку для решения матриц

Таких в свободном доступе куча, а создавать свое - это, ИМХО, разработка велосипедов:)
В МКЭ есть много других задач, на которые можно убить много своего времени.

Второй момент, так к сведению, в коммерческих программах зачастую используют обычные спарс-матрицы для хранения и метод сопряженных градиетов с предобуславливанием неполным разложением Холеского (ICPCG, Incomplete Cholesky Preconditioned Conjaguate Gradient).
Это связано с тем, что матрица жесткости сильно разряжена, очень плохо обусловлена, и применение итерационных решателей дает сильный выигрыш.
В этом случае ширина ленты не отражается ни на хранении, ни на решении, и задачи уменьшения ширины ленты просто отпадают.

[vitaly333]

Цитата:

спарс-матрицы

Что такое спарс-матрицы.

Цитата:

Таких в свободном доступе куча

Нашел только одну на офф. сайте МГУ.

Цитата:

Второй момент, так к сведению, в коммерческих программах зачастую используют обычные спарс-матрицы для хранения и метод сопряженных градиетов с предобуславливанием неполным разложением Холеского (ICPCG, Incomplete Cholesky Preconditioned Conjaguate Gradient).
Это связано с тем, что матрица жесткости сильно разряжена, очень плохо обусловлена, и применение итерационных решателей дает сильный выигрыш.
В этом случае ширина ленты не отражается ни на хранении, ни на решении, и задачи уменьшения ширины ленты просто отпадают.

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

Цитата:

Reverse Cuthill-McKee (вроде так пишется) предназначен для перенумерации узлов сетки так, чтобы потом ширина лента была меньше исходной и работает он с графом, получаемом как из сетки, так уже и из глобальной матрицы. Таким образом, им можно пройтись по сетке КЭ сразу после ее создания и перенумеровать узлы один-единственный раз, а потом считать, считать и считать...

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

[незваный гость]

vitaly333 писал(а):

Что такое спарс-матрицы.

разреженные матрицы.

[VLarin]

Спарс-матрица - от англ. Sparse-matrix:)

Цитата:

Но всё таки всё равно надо делать оптимизацию матрицы (приведение к ленточному виду без уменьшения ширины ленты). Не хранить же в памяти матрицу размером 3,4 Гб(для случая с 10-ю тысячами узлов)

Нет, ни какой оптимизации делать не нужно. Спрас-матрицы хранят только ненулевые элементы.
В каждой строке матрицы будет всего N-ненулевых элемента, где N, как я уже писал, число узлов, с которыми связан i-узел посредством конечных элементов + 1, умноженное на число степеней свободы в узле.
Если у Вас конечные элементы 1 порядка четырехугольные - то это (8+1)*3, треугольные (6+1)*3. На основании этого можно прикинуть сколько у Вас будет ненулевых элементов в матрице и сколько для этого нужно памяти, чтобы хранить их в виде спарс-матрицы. А еще и симметрия матрицы жесткости.

Цитата:

Дайте ссылочку на этот метод сопряженных градиетов с предобуславливанием неполным разложением Холеского.

Библиотека на шаблонах для хранения и работы с матрицами
Matrix Template Library (ITL)
Итерационные решатели с различными предобуславливателями
Iterative Template Library (ITL)

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

Если интересует сам метод, то можно посмотреть
Методы решения СЛАУ большой размерности

Цитата:

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

В принципе да. Если все-таки интересует уменьшение ширины ленты - то все как я уже говорил - сначала делаем сетку, потом перенумеровываем узлы и сохраняем. Специальный граф сущестует только на этапе перенумерации, потом он просто не нужен.
В итоге до стадии сборки матрицы будем иметь сетку конечных элементов такую, что она априори дает хорошую ширины ленты. Но опять таки - это нужно только для хранения в ленточном виде и применения соответствующих решателей.
Для связки спарс-матрица+CG (Conjaguate Gradient) этого не требуется:)

Главное хорошо продумать иерархию классов. По этому вопросу в инете есть куча работ типа "Объектно-ориентированный МКЭ" - диссеры, статьи, сайты с иерархией классов. Мой Вам совет - не нужно изобретать велосипед, нужно хорошенько ознакомиться с тем что уже есть и на основании этого додумав и доделав что-то сотворить своё. Неплохо было бы реализовать супер-элементы и domain decomposition, это дает ускорение решения.

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

Добавлено спустя 4 минуты 47 секунд:

Цитата:

Почитал книжку А.Джордж , Дж.Лю - Численное решение больших разряженных систем

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

[vitaly333]

Цитата:

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

Например. . ... можете дать ссылки на электронные(бесплатные(pdf,djvu)) версии.
Есть ещё две книжки по Sparse-матриам:
1. С.Писсанецки, Технология разреженных матриц, М, Мир, 1988
2. Тьюарссон «Разреженные матрицы», М,Мир,1977

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

Цитата:

Спрас-матрицы хранят только ненулевые элементы.
В каждой строке матрицы будет всего N-ненулевых элемента, где N, как я уже писал, число узлов, с которыми связан i-узел посредством конечных элементов + 1, умноженное на число степеней свободы в узле.
Если у Вас конечные элементы 1 порядка четырехугольные - то это (8+1)*3, треугольные (6+1)*3. На основании этого можно прикинуть сколько у Вас будет ненулевых элементов в матрице и сколько для этого нужно памяти, чтобы хранить их в виде спарс-матрицы. А еще и симметрия матрицы жесткости.

Для хранения Sparse -матриц как я понял также есть множество способов:
(старые)
1.РСФ(SCR)(Разряженный строчный формат. Матрица представляется в виде 3-х одномерных массивов).
2.РCтФ(Разряженный столбцовый вариант).
(более новые)
1.Хэш-таблицы.
2.Линейный список из ненулевых элементов.
3. Массив указателей на строки. В строках линейные списки.
Так какой метод хранения лучше использовать для итерационного решателя ICPCG?

Цитата:

[Comga]

По итерационным методам разреженных систем см. книгу Saad. Iterative Methods for Sparse Linear Systems, 2nd edition.
Первое издание книги Saad. Iterative Methods for Sparse Linear Systems
Там же имеется пакет для разреженных матриц SPARSKIT. Весьма эффективен.

[vitaly333]

Цитата:

По итерационным методам разреженных систем см. книгу Saad. Iterative Methods for Sparse Linear Systems, 2nd edition.
Первое издание книги Saad. Iterative Methods for Sparse Linear Systems

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

Добавлено спустя 16 минут 11 секунд:

Saad. Iterative Methods for Sparse Linear Systems, 2nd edition. - посмотрел, очень хорошая книжка, жаль что на англисском, может где-то в переводе есть?

[VLarin]

Цитата:

Например. . ... можете дать ссылки на электронные(бесплатные(pdf,djvu)) версии.

А для чего по Вашему я давал ссулку на Методы решения СЛАУ большой размерности

Там есть ссыла на книгу в формате pdf:)
В книге рассмотрены вопросы хранения матриц, простых операций, предобуславливателей и итерационных решателей. Ничего лучше я на русском в инете не видел. Что самое главное - книга написана применительно к МКЭ, акцент делается на том, что матрицы симметричные, сильно разрежены.

Для хранения спарс-матриц обычно используют три массива. Хранить можно по строкам, можно по столбцам. Если храним по строкам, то имеем - массив ненулевых элементов, массив строк (с какого элемента начинается новая строка) и массив индексов столбцов (индексы ненулевых элементов). Это главная идея и сделать ее можно по-разному.
Так матрицу мы создаем динамически и можем и не знать заранее число ненулевых элементов и их положение, то тут нужно иметь в виду списки.

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

Все-таки еще раз советую посмотреть код MTL и других библиотек.

[Comga]

Подкину еще ссылочку на английском

H. Elman, D. Silvester, A. Wathen. Finite Elements and Fast Iterative Solvers: with Applications in Incompressible Fluid Dynamics. Oxford University Press, 2005.

http://www.manchester.ac.uk/ifiss
http://www.cs.umd.edu/~elman/ifiss.html

[vitaly333]

Цитата:

В каждой строке матрицы будет всего N-ненулевых элемента, где N, как я уже писал, число узлов, с которыми связан i-узел посредством конечных элементов + 1, умноженное на число степеней свободы в узле.
Если у Вас конечные элементы 1 порядка четырехугольные - то это (8+1)*3, треугольные (6+1)*3. На основании этого можно прикинуть сколько у Вас будет ненулевых элементов в матрице и сколько для этого нужно памяти, чтобы хранить их в виде спарс-матрицы.

У меня в програме пока не предусмотрен генеранотор сетки конечных элементов и расчетную схему пользователь создает вручную введя сначала координаты узлов а потом связывает узлы конечными элементами. Система считает стержневые системы(балки, рамы , фермы)
Никахих трехугольгных и четырехугольных элементов нет.
Как определить заранее, перед сборкой матрицы жесткости количество ненулевых элементов в такой системе. (N-1)*кол.степеней свободы - эта формула не прокатывает в этом случае. Пользователь ведь может создать любую расчетную схему. Например пользователь введет три узла , свяжет 1-ый узел со 2-ым посредством 1 -ого конечного элемента, свяжет 2-ой узел с 3-им посредством конечного 2-ого конечного элемета и свяжет 1-ый узел с 3-им посредством 3-его конечного элемента. Тогда матрица жесткости получается плотная без ненулевых элементов.

[VLarin]

vitaly333 писал(а):

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

На момент сборки матрицы сетка КЭ уже создана - в том, или ином виде. Иначе о какой сборки может идти речь.

Проходим по сетке. Определяем, с какими элементами связан i-узел, какие узлы содержат эти элементы, что в итоге позволяет оценить число уникальных узлов, с которыми связан i-узел. Соответственно, как я уже не раз отмечал ранее, в i-строке число ненулевых элементов равно числу узлов, с которыми связан текущий i-узел. Проход по всем узлам сетки КЭ позволит определить количество ненулевых элементов и их положение в матрице жесткости, то есть портрет матрицы жесткости.

[Cow-Brat]

VLarin в сообщении #42550 писал(а):

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

Здравствуйте!
Я получил матрицу жесткости для стержневого элемента, несколько отличающуюся от привычной. В ней я пытаюсь учесть деформацию депланации, характерную для тонкостенных стержней. Размерность матрицы в общем случае 14х14. Часть отвечающая только за кручение и депланацию имеет размерность 4х4.

Собственно вопрос. Хотелось бы опробовать полученную матрицу "в деле". Как я понял, почитав эту тему, и ряд смежных, написание своей программы очень трудоемкий процесс. Но я никак не могу понять, как в моем случае можно использовать уже готовые решения, подменив в них "чужую" матрицу жесткости на свою?