Неполное разложение Холецкого разреженных матриц

[vitaly333]

Цитата:

Да вроде разложилась

Да. А можете выложить файл с разложенным вектором AN. И ещё выложите пожалуйста, файл с разложенной той 1-ой матрицей 7x7.

Цитата:

Кроме того, проверил по IC_Decomposition2 - нулевых подкоренных значений там нет

Не может быть Несколько раз раскладывал IC_ Decomposition2 в самом конце получаются нули. А при решении происходит на них деление и ……

Как я понял вы используете тот же самый алгоритм из ITL. Он же ни чем не отличается от моего аналога IC_Decomposition2 . Тогда почему у вас она раскладывается а у меня нет?
Вы ничего в нем не меняли, просто взяли ctrl+c и ctrl+v? И ещё как у вас в программе организовано хранение разреженных матриц. А именно:
Какую треугольную часть, включая диагональ, храните над диагональную или под диагональную? И как храните по столбцам или по строкам?

Цитата:

73659 - это уже перебор. Так обычно максимум за 1000 итераций все сходится

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

Цитата:

Потом с какой точностью сошлось - относительной/абсолютной.

0.001

Цитата:

И, кстати, почему Холецкий, когда он вроде Холесский?

Незнаю…. В одних книгах пишут Холецкий в других Холесский. Мне больше понравилось Холецкий

Цитата:

Еще хороший вариант проверить решение Вашей матрицы - это MATLAB?

Попробовал. Не получилось создать там sparse – матрицу из моих 3-х массивов. Я тогда сделал немного по другому. Я в своей программе из этих 3-х массивов воссоздал полную матрицу , загрузил её в матлаб и там уже из неё создал sparse. Далее пробовал решить, воспользовавшись командой pcg . Матлаб долго думал и выдал флаг = 1 , означающий то что количество итераций превысило максимальное количество, которое я задал равным 10000.

[VLarin]

Цитата:

Кроме того, проверил по IC_Decomposition2 - нулевых подкоренных значений там нет

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

Матрицу AN2 после IC_Decomposition2 выкладываю тут.

Хранение матриц у меня в программе -

typedef matrix<Type, symmetric<lower>, array< compressed<> >, column_major>::type Matrix;

Такой же тип матриц приводится в примерах по ITL для решателя CG в связке с cholesky.

IC_Decomposition2 скопировал Ваш.

МАТЛАБ - попробуйте стандартный решатель.
Плюс сделать код, который по сжатой матрице создает в матлабе полную:)

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

[vitaly333]

Цитата:

У Вас ранее было написано, что под корнем отрицательные значения появляются. Я проверил, у меня их нет.

Я тоже немного напутал.. под корнем отрицательные значения появляются это при разложении IC_Decomposition (мой) . А при разложении IC_Decomposition2 в самом конце появляются нули.
Посмотрел ваш файл AN2_mod3000 действительно разложил!!!
И нулей в конце нет.... странно у меня совсем другое разложение получается. Вы точно использовали мой IC_Decomposition2?А ваш ITL -овский дает тоже самое?

Цитата:

typedef matrix<Type, symmetric<lower>, array< compressed<> >, column_major>::type Matrix;

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



Попробовал использовать ваш разложенный AN_mod3000 в качестве предобуславливателя, но cg выдал вектор неизвестных весь NAN после 1073 итерации. Решил разобраться из-за чего, оказывается при нахождении скалярного произведения в CG производиться умножение примерно таких огромных чисел:

-2.574249076548091E299 * 1.9501677570630539E304

после котого результат такого произведения естественно не влезает в формат double и записывает infinity. А потом ещё в CG при нахождении alfa происходит деление double/infinity = NAN. Вот, откуда и появляется вектор, где все элементы его NAN. И поэтому и у вас наверно и не сходится. Только тут вопрос откуда взялись такие огромные числа?

Цитата:

Может попробывать перенумерацию узлов сделать.

Где то яслышал что пед тем как применять алгоритм Неполного Холецкого нужно предварительно сделать перестановку строк , чтобы максимально сократить появление новых ненулевых элементов в процессе факторизации в тех местах где раньше стояли нули в исходной матрице до факоризации. Один из самых лучших таких алгоритмов это Minimal Degree (Алгоритм минимальной степени). Но если это делать то это уже получается не итерационное решение а прямое, потому что при применении такого алгоритма Неполное разложение будет максимально приближено к полному, а в некоторых случаях будет полным. И как вы ранее говорили

Цитата:

это значит, что всю работу делает разложение, а не решатель.
На то он и итерационный, чтоб решать за несколько итераций

Тем более в этом примере мне кажется что и даже перенумерация не поможет так как расчетная схема - это арка. Узел 1 связан с узлом 2 - конечным элементом 1 ; Узел 2 связан с узлом 3 -конечным элементом 2 ........ Узел 999 связан с узлом 1000 - конечным элементом 999. Тут как узлы не перенумеровывай всё равно получается 6 ненулевых в строке. Хотя я могу ошибаться.

Кроме того я читал в одном форуме что в переобуславливателе заполнение нужно контролировать порогами.Т.е ненулевые элементы динамически отбрасываются по порогу и при факторизации должна постоянно поддерживаться положительно определенность(алгоритм Азиза Дженнингса). Но ничего по этому в инэте не нашел. Может вы что - то знаете по этому поводу?

Кроме того есть ещё какая-то модернизация неполного разложения Холецкого со сдвигом (со Shift - ом , так её называют). Но тоже поиски в инэте были без результатны

[VLarin]

Цитата:

И нулей в конце нет.... странно у меня совсем другое разложение получается. Вы точно использовали мой IC_Decomposition2?А ваш ITL -овский дает тоже самое?

Точно, и ITL такое же дает:)

Цитата:

А при разложении IC_Decomposition2 в самом конце появляются нули.

А Вы размер массива на выходе не проверяли - исходный 14978, а Ваш - 17978. Последние 3000 элементов - это нули:))))
Причем судя по последним (ненулевым) элементам - вроде то, что должно быть.

Цитата:

Попробовал использовать ваш разложенный AN_mod3000 в качестве предобуславливателя, но cg выдал вектор неизвестных весь NAN после 1073 итерации. Решил разобраться из-за чего, оказывается при нахождении скалярного произведения в CG производиться умножение примерно таких огромных чисел:

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

Цитата:

Тем более в этом примере мне кажется что и даже перенумерация не поможет так как расчетная схема - это арка. Узел 1 связан с узлом 2 - конечным элементом 1 ; Узел 2 связан с узлом 3 -конечным элементом 2 ........ Узел 999 связан с узлом 1000 - конечным элементом 999. Тут как узлы не перенумеровывай всё равно получается 6 ненулевых в строке. Хотя я могу ошибаться.

Просто я не очень понял, что за задача и что Вы понимаете под аркой (что такое арка я знаю:)) - я механикой не занимался. Это 1D? Тогда почему 6 ненулевых элементов матрицы в строке, когда их должно быть 3*(число степеней свободы)?
Просто странно все это, что задача не решается. Может быть взять другой тестовый пример, более приближенный к реальности?:)
Во всем мире считают cholesky+CG и вроде не жалуются:)

Цитата:

Кроме того я читал в одном форуме что в переобуславливателе заполнение нужно контролировать порогами.Т.е ненулевые элементы динамически отбрасываются по порогу и при факторизации должна постоянно поддерживаться положительно определенность(алгоритм Азиза Дженнингса). Но ничего по этому в инэте не нашел. Может вы что - то знаете по этому поводу?

Кроме того есть ещё какая-то модернизация неполного разложения Холецкого со сдвигом (со Shift - ом , так её называют). Но тоже поиски в инэте были без результатны

Таких тонкостей я честно сказать не знаю:(

[vitaly333]

Цитата:

А Вы размер массива на выходе не проверяли - исходный 14978, а Ваш - 17978. Последние 3000 элементов - это нули )))
Причем судя по последним (ненулевым) элементам - вроде то, что должно быть.

Да точно 17978. Вот где прокол был.

Цитата:

Тогда почему 6 ненулевых элементов матрицы в строке, когда их должно быть 3*(число степеней свободы)?

Ну как бы вам объяснить…. давайте на примере. Пусть дана арка состоящая из 3 – х узлов. Тогда глобальная матрица жесткости будет размером 3*3 = 9x9.
Т.к узел 1 соединен с узлом 2 конечным элементом 1 то в первых – 3 –х строчках матрицы будет по 6 ненулевых элементов = 3 от узла 1 (т.к 3 степени свободы в узле 1) и 3 от узла 2 , т.к (во 2 –ом узле также 3 степени свободы.) Т.к 1 –ый узел не соединен с 3-им то последние 3-ри элемента нулевые,и.т.д



1 2 3 4 5 6 7 8 9
1 * * * * * * 0 0 0
2 * * * * * * 0 0 0
3 * * * * * * 0 0 0
4 * * * * * * * * *
5 * * * * * * * * *
6 * * * * * * * * *
7 0 0 0 * * * * * *
8 0 0 0 * * * * * *
9 0 0 0 * * * * * *


Да кстати, взял тут одинпример, так не реальной задачи а сам придумал на обум. Матрицу взял 10000x10000, все диагональные элементы равны 2. в каждой строке по 16 ненулевых, включая диагональ. Попробовал её решить cg без предобуславливания , у меня cg сошелся за 7 итераций. С предобуславливателем IC_Decomposition cg сошелся за 1 итерацию – опять всю работу сделало разложение , продолжительность которого составило аж 70 секунд!!! С предобуславливанием IC_Decomposition2 cg сошелся на 3997 итерации. И тут я совсем запутался какой лучше предобуславливатель использовать? А может лучше без него? Конечно такая задача которую я придумал вред ли встретиться , но всё таки почему ITL – овский предобуславливатель даёт такую плохую сходимость. Попробуйте у себя решить, может быть я опять что – то напутал….Да кстати во всех 3-х ситуациях получается немного разное решение….

Цитата:

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

Да я бы с радостью только надо заморачиваться с генератором сетки. Сгенерируйте мне, если вам не трудно, хотя бы один примерчик. И мы сверим заодно результаты!!!

[VLarin]

Цитата:

Ну как бы вам объяснить…. давайте на примере. Пусть дана арка состоящая из 3 – х узлов. Тогда глобальная матрица жесткости будет размером 3*3 = 9x9.
Т.к узел 1 соединен с узлом 2 конечным элементом 1 то в первых – 3 –х строчках матрицы будет по 6 ненулевых элементов = 3 от узла 1 (т.к 3 степени свободы в узле 1) и 3 от узла 2 , т.к (во 2 –ом узле также 3 степени свободы.) Т.к 1 –ый узел не соединен с 3-им то последние 3-ри элемента нулевые,и.т.д

Только 6 ненулевых будет для первых и последних трех строк, а для всей остальной части матрицы - 3*(число степеней свободы).
Как я уже писал ранее, число ненулевых элементов в строке матрицы определяется просто: это 1 + число узлов, с которыми соединяется некоторый узел посредством конечных элементов, и все это дело унижаем на число степеней свободы в узле.

Цитата:

Да я бы с радостью только надо заморачиваться с генератором сетки. Сгенерируйте мне, если вам не трудно, хотя бы один примерчик. И мы сверим заодно результаты!!!

Ну Вам в любом случае через некоторое время прийдется заморачиваться с генератором сетки:)
Попробуте GeomPack++ - там достаточно просто формат входых и выходных данных.

[vitaly333]

Цитата:

Ну Вам в любом случае через некоторое время прийдется заморачиваться с генератором сетки

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

Ну вы попробовали примерчик с матрицей 10000x1000?

[Comga]

2-я МЕЖДУНАРОДНАЯ КОНФЕРЕНЦИЯ
"МАТРИЧНЫЕ МЕТОДЫ И ОПЕРАТОРНЫЕ УРАВНЕНИЯ"
23-28 июля 2007 года, Москва
http://bach.inm.ras.ru/album_vi/conf2007/zero_r.htm

[metallica4ever]

Здравствуйте! Скажите, пожалуйста, а где можно взять текст библиотеки ITL или IML? Конкретно интересует текст алгоритма CG.