2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




Начать новую тему Ответить на тему На страницу Пред.  1, 2
 
 
Сообщение03.04.2007, 18:10 


26/11/06
76
Цитата:
Да вроде разложилась


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

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


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

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

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


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

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

0.001

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


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

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

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

 Профиль  
                  
 
 
Сообщение04.04.2007, 16:37 


13/09/05
153
Москва
Цитата:
Кроме того, проверил по IC_Decomposition2 - нулевых подкоренных значений там нет

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

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

Хранение матриц у меня в программе -
Код:
typedef matrix<Type, symmetric<lower>, array< compressed<> >, column_major>::type Matrix;

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

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

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

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

 Профиль  
                  
 
 
Сообщение05.04.2007, 16:03 


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

Я тоже немного напутал.. :) под корнем отрицательные значения появляются это при разложении 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 - ом , так её называют). Но тоже поиски в инэте были без результатны

 Профиль  
                  
 
 
Сообщение05.04.2007, 18:45 


13/09/05
153
Москва
Цитата:
И нулей в конце нет.... странно у меня совсем другое разложение получается. Вы точно использовали мой 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 - ом , так её называют). Но тоже поиски в инэте были без результатны

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

 Профиль  
                  
 
 
Сообщение10.04.2007, 14:56 


26/11/06
76
Цитата:
А Вы размер массива на выходе не проверяли - исходный 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-х ситуациях получается немного разное решение….


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

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

 Профиль  
                  
 
 
Сообщение11.04.2007, 10:15 


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


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

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

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

 Профиль  
                  
 
 
Сообщение12.04.2007, 14:23 


26/11/06
76
Цитата:
Ну Вам в любом случае через некоторое время прийдется заморачиваться с генератором сетки

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

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

 Профиль  
                  
 
 Может быть полезным
Сообщение15.05.2007, 07:55 


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

 Профиль  
                  
 
 Re: Неполное разложение Холецкого разреженных матриц
Сообщение06.05.2011, 11:53 


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

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 24 ]  На страницу Пред.  1, 2

Модераторы: Karan, Toucan, PAV, maxal, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group