10
Реферат
Дипломна робота: 72 с., 8 рис., 5 табл., 1 додаток, 23 джерела.
Мета роботи - розробка алгоритму та програми чисельного розвязку систем лінійних алгебраїчних рівнянь з розрідженою матрицею.
При виконанні роботи використані програми Microsoft Office Word 2007 Service Pack 1, Microsoft Office PowerPoint 2007 Service Pack 1 та Microsoft Visual Studio 2008 Service Pack 1, а також основні поняття лінійної алгебри та математичного моделювання.
В результаті роботи розроблена програма чисельного розвязку систем лінійних алгебраїчних рівнянь з розрідженою матрицею, яка економно витрачає оперативну память, що дозволяє розвязувати багато систем високих ступенів за допомогою персональних компютерів. Для розвязання таких систем при класичній схемі зберігання всіх елементів матриці в оперативній памяті довелось би залучати суперкомпютери.
СИСТЕМА ЛІНІЙНИХ АЛГЕБРАЇЧНИХ РІВНЯНЬ, ОПЕРАТИВНА ПАМЯТЬ, МАТРИЦЯ, АЛГОРИТМ, ПРОГРАМА.
Зміст
Перелік умовних скорочень і термінів
Вступ
1 Огляд методів розвязку СЛАР, що виникають у МСЕ
1.1 Точні методи розвязку СЛАР
1.1.1 Метод Гауса
1.1.2 Метод Крамера
1.1.3 Метод головних елементів
1.1.4 Схема Халецького
1.1.5 Метод квадратного кореня
1.1.6 Метод прогону
1.1.7 Матричний метод
1.2 Ітераційні методи розвязку СЛАР
2.1 Метод простих ітерацій
1.2.2 Метод Зейделя
1.2.3 Метод релаксації
1.2.4 Багатосітковий метод
1.2.5 Метод Ланцоша
2 Схеми компактного зберігання розріджених матриць
2.1 Перша схема
2.2 Друга схема
3 Оптимізація обчислень
4 Чисельні експерименти
4.1 Пружне деформування тонкостінної просторової рами
4.2 Контактна взаємодія оболонкової конструкції і ложемента
Висновки
Перелік посилань
Додаток А
49
Перелік умовних скорочень і термінів
1. Умовні скорочення
МСЕ - метод скінчених елементів.
СЛАР - система лінійних алгебраїчних рівнянь.
2. Терміни
Розмірність матриці - кількість невідомих в матриці.
Розріджена матриця - матриця, в якій більшість елементів дорівнює нулю.
Оперативна память - енергозалежний вид памяті компютера, в якій зберігаються дані, необхідні процесору в даний час.
Крайова задача - система диференціальних рівнянь з визначеними лінійними співвідношеннями між значеннями шуканих функцій на кінцях інтервала інтегрування.
Вступ
Зараз у звязку з бурхливим розвитком обчислювальних засобів широке поширення одержали інформаційні технології, що мають різноманітну теоретичну і прикладну спрямованість. Серед них особливе місце займають системи автоматизованого проектування (САПР), невідємну частину яких становлять підсистеми математичного моделювання різних фізичних процесів. Одним із самих перспективних напрямків розвитку математичного моделювання є широке використання чисельних методів, таких як метод скінчених елементів (МСЕ). Це чисельний метод для диференціальних рівнянь, що зустрічаються у фізиці [1]. Виникнення цього методу повязане з розвязком задач космічних досліджень. Уперше він був опублікований у роботі Тернера, Клужа, Мартіна і Топпа. Ця робота сприяла появі інших робіт; був опублікований ряд статей із застосуванням методу скінчених елементів до задач будівельної механіки і механіки суцільних середовищ. Важливий внесок у теоретичну розробку методу зробив в 1963 р. Мелош, який показав, що метод скінчених елементів можна розглядати як один з варіантів добре відомого методу Релєя-Рітца. У будівельній механіці метод скінчених елементів мінімізацією потенційної енергії дозволяє звести задачу до системи лінійних рівнянь рівноваги [2,3].
Однією з існуючих труднощів, що виникають при чисельній реалізації розвязку контактних задач теорії пружності методом скінчених елементів, є розвязок систем лінійних алгебраїчних рівнянь (СЛАР) великої розмірності виду . За умови зберігання всіх елементів матриці A в оперативній памяті навіть компютер з 24 ГБ оперативної памяті здатен вирішити СЛАР з кількістю невідомих не більше 56755. Якщо ж невідомих більше, доведеться залучати суперкомпютер. Тому актуальною є проблема знаходження схеми компактнішого зберігання матриці СЛАР в оперативній памяті.
Зараз розроблена велика кількість програмних систем, що автоматизують різні аспекти чисельного моделювання проблем механіки (ANSYS, COSMOS, FORTU, МІРЕЛА, ПОЛЕ, ПРОЧНОСТЬ і ін.).
Одним з головних блоків усіх цих систем є модуль, відповідальний за розвязок СЛАР. Слід зазначити, що звичайно розвязок СЛАР займає більшість часу розрахунку задачі на компютері. Крім того, від цього етапу прямо залежить якість і точність розвязку задачі в цілому. Тому від вибору методу розвязку СЛАР і алгоритму його реалізації багато в чому залежить подальший успіх чисельного розрахунку задачі МСЕ.
Більшість існуючих методів розвязку СЛАР в МСЕ розроблені в припущенні того, що матриця A має стрічкову структуру, причому ширина стрічки m<<n, де n?n - розмірність A. Однак, при використанні МСЕ для чисельного розвязку контактних задач можливі випадки, коли ширина стрічки m>n.
1. Огляд методів розвязку СЛАР, що виникають у МСЕ
У процесі побудови дискретних аналогів крайових задач виникають великі системи в загальному випадку нелінійних алгебраїчних рівнянь, які вирішуються у два етапи: на першому етапі вони линеаризуються, а потім отримана система лінійних рівнянь вирішується за допомогою якого-небудь методу лінійної алгебри. Якщо збіжність не досягнута, то процес повторюється.
Основна ідея методу скінчених елементів полягає в тому, що будь-яку неперервну величину, таку як температура, тиск і переміщення, можна апроксимувати дискретною моделлю, яка будується на множині кусочно-неперервних функцій, визначених на кінцевім числі підобластей. Кусочно-неперервні функції визначаються за допомогою значень неперервної величини в кінцевім числі точок розглядаємої області [1,2,3].
У загальному випадку, неперервна величина заздалегідь невідома, і потрібно визначити значення цієї величини в деяких внутрішніх точках області. Дискретну модель, однак, дуже легко побудувати, якщо спочатку припустити, що числові значення цієї величини в кожній внутрішній точці області відомі. Після цього можна перейти до загального випадку. Отже, при побудові конкретної моделі неперервної величини роблять наступним чином:
а) у розглядаємій області фіксується кінцеве число точок. Ці точки називаються вузловими точками або просто вузлами;
б) значення неперервної величини в кожній вузловій точці вважається змінною, яка повинна бути визначена;
в) область визначення неперервної величини розбивається на кінцеве число підобластей, називаних елементами. Ці елементи мають загальні вузлові точки і у сукупності апроксимують форму області;
г) неперервна величина апроксимується на кожному елементі функцією, яка визначається за допомогою вузлових значень цієї величини. Для кожного елемента визначається своя функція, але функції підбираються таким чином, щоб зберігалася неперервність величини уздовж границь елемента.
Способи розвязку СЛАР, застосовувані разом із МСЕ, можна умовно розбити на дві групи: прямі (точні) і ітераційні (наближені). Прийнято вважати, що прямі методи найбільш ефективні при рішенні СЛАР невисоких порядків (не більше 105 невідомих), а ітераційні - для розвязку СЛАР великих і надвеликих систем (понад 105 невідомих).
Якщо всю матрицю коефіцієнтів СЛАР вдається розмістити в памяті, то розвязок такої системи не представляє якої-небудь складності. Однак на практиці часто доводиться мати справу з матрицями, розміри яких суттєво перевищують обєм оперативної памяті компютера. У таких випадках і доводиться використовувати різні модифікації як прямих, так і ітераційних методів.
Остаточне рішення про застосування прямих або ітераційних методів розвязку СЛАР необхідно ухвалювати на основі аналізу структури досліджуваної математичної задачі. Прямі методи розвязку СЛАР вигідніше використовувати, якщо необхідно розвязувати багато однакових систем з різними правими частинами, або якщо матриця А не є додатньо-визначеною. Крім того, існують задачі з такою структурою матриці, для якої прямі методи завжди кращі, ніж ітераційні.
1.1 Точні методи розвязку СЛАР
Розглянемо ряд точних методів розвязку СЛАР [4,5].
1.1.1 Метод Гауса
Нехай дана система , де A - матриця розмірності m?m (квадратна). У припущенні, що a11?0, перше рівняння системи
ділимо на коефіцієнт a11, у результаті одержуємо рівняння:
Потім з кожного з інших рівнянь віднімається перше рівняння, помножене на відповідний коефіцієнт ai1. У результаті ці рівняння перетворяться до виду:
Перше невідоме виявилося виключеним із усіх рівнянь, крім першого. Далі в припущенні, що , ділимо друге рівняння на коефіцієнт і виключаємо невідоме з усіх рівнянь, починаючи з третього і т.д. У результаті послідовного виключення невідомих, матриця системи рівнянь стане трикутною:
Сукупність проведених обчислень називається прямим ходом методу Гауса.
З m-го рівняння системи визначаємо xm, з (m-1)-го рівняння визначаємо xm-1 і т.д. до x1. Сукупність таких обчислень називають зворотнім ходом методу Гауса.
Реалізація прямого методу Гауса вимагає N~2m2/3 арифметичних операцій, а зворотнього - N~m2 арифметичних операцій.
1.1.2 Метод Крамера
Нехай дана система лінійних рівнянь, в якій число рівнянь дорівнює числу невідомих:
Припустимо, що визначник системи d?0. Якщо тепер замінити послідовно у визначнику стовпці коефіцієнтів при невідомих xj стовпцем вільних членів bi, то вийдуть відповідно n визначників d1,...,dn.
Теорема Крамера. Система n лінійних рівнянь з n невідомими, визначник якої відмінний від нуля, завжди сумісна і має єдиний розвязок, який обчислюється по формулах:
Розвязок довільних систем лінійних рівнянь. Нехай
довільна система лінійних рівнянь, де число m рівнянь системи менше числа n невідомих. Тоді в матриці СЛАР знайдуться r?m лінійно незалежних рядків, а інші m-r рядків виявляться їхніми лінійними комбінаціями. Перестановкою рівнянь можна добитися того, що ці r лінійно незалежних рядків займуть перші r місць.
Звідси випливає, що кожне з останніх m-r рівнянь системи (1.7) можна представити як суму перших r рівнянь (які називаються лінійно незалежними або базисними), узятих з деякими коефіцієнтами. Тоді система еквівалентна наступній системі r рівнянь з n невідомими:
Припустимо, що мінор r-ї розмірності, складений з коефіцієнтів при перших r невідомих, відмінний від нуля: Мr?0, тобто є базисним мінором. У цьому випадку невідомі, коефіцієнти при яких складають базисний мінор, називаються базисними невідомими, а інші n-r - вільними невідомими.
У кожному з рівнянь системи (1.8) перенесемо в праву частину всі члени з вільними невідомими xr+1,...,xn. Тоді одержимо систему, яка містить r рівнянь з r базисними невідомими. Оскільки визначник цієї системи є базисний мінор Mr, то система має єдиний розвязок щодо базисних невідомих, який можна знайти по формулах Крамера. Даючи вільним невідомим довільні числові значення, одержимо загальний розвязок вихідної системи.
Однорідна система лінійних рівнянь. Нехай дана однорідна система лінійних рівнянь з n невідомими:
Оскільки додавання стовпця з нулів не змінює рангу матриці системи, то на підставі теореми Кронекера-Kaпeллі ця система завжди сумісна і має, принаймні, нульовий розвязок. Якщо визначник системи відмінний від нуля і число рівнянь системи дорівнює числу невідомих, то по теоремі Крамера нульовий розвязок є єдиним.
У тому випадку, коли ранг матриці системи менше числа невідомих, дана система крім нульового розвязку буде мати і ненульові розвязки. Для знаходження цих розвязків у системі (1.9) виділяємо r<n лінійно незалежних рівнянь, інші відкидаємо. У виділених рівняннях у лівій частині залишаємо r базисних невідомих, а інші n-r вільних невідомих переносимо в праву частину. Тоді приходимо до системи, розвязуючи яку по формулах Крамера, виразимо r базисних невідомих x1,...,хr через n-r вільних невідомих.
1.1.3 Метод головних елементів
Нехай дана система n лінійних рівнянь з n невідомими:
де елементи aij (i,j=1,…,n) утворюють розширену матрицю системи . Виберемо найбільший по модулю і не належачий стовпцю вільних членів елемент apq матриці , який називається головним елементом, і обчислимо множники mi=-aiq /apq для всіх рядків з номерами i?p (р-й рядок, що містить головний елемент, називається головним рядком).
Далі до кожного другорядного i-го рядку додамо головний рядок, помножений на відповідний множник mi.
У результаті одержимо нову матрицю, усі елементи q-го стовпця якої, крім apq, складаються з нулів.
Відкинувши цей стовпець і головний p-й рядок, одержимо нову матрицю, число рядків і стовпців якої на одиницю менше. Повторюємо ті ж операції з отриманою матрицею, після чого одержуємо нову матрицю і т.д.
Таким чином, побудуємо послідовність матриць. Для визначення невідомих xj поєднуємо в систему всі головні рядки, починаючи з останнього.
Викладений метод розвязку систем лінійних рівнянь називається методом головних елементів. Необхідна умова його застосування полягає в тому, що визначник матриці не дорівнює нулю [6,7].
1.1.4 Схема Халецького
Нехай система лінійних алгебраїчних рівнянь дана в матричному вигляді: , де А - квадратна матриця розмірності n; , - вектори-стовпці.
Представимо матрицю А у вигляді добутку нижньої трикутної матриці С і верхньої трикутної матриці В з одиничною діагоналлю, тобто А=СВ, де
Причому елементи сij і bij визначаються по формулах:
Рівняння можна записати в наступному вигляді:
Добуток матриці B на вектор-стовпець є вектором-стовпцем, який позначимо через :
Тоді рівняння (1.13) перепишемо у вигляді:
Тут елементи сij відомі, тому що матриця А системи вважається вже розкладеною на добуток двох трикутних матриць С і В.
Перемноживши матриці в лівій частині рівності (1.15), одержуємо систему рівнянь, з якої одержимо наступні формули для визначення невідомих:
Невідомі yi зручно обчислювати разом з елементами bij.
Після того, як усі yi визначені по формулах (1.16), підставляємо їх у рівняння (1.14) [8].
Оскільки коефіцієнти bij визначені в (1.12), то значення невідомих, починаючи з останнього, обчислюємо по наступних формулах:
1.1.5 Метод квадратного кореня
Даний метод використовується для розвязку систем лінійних алгебраїчних рівнянь виду
у яких матриця А симетрична, тобто АТ=А, aij=aji (i=j=1,...,n).
Розвязок системи (1.18) здійснюється у два етапи.
Прямий хід. Перетворення матриці А і представлення її у вигляді добутку двох взаємно транспонованих трикутних матриць:
де
Перемножуючи SТ і S, і дорівнюючи матриці А, одержимо наступні формули для визначення sij:
Після знаходження матриці S систему (1.18) заміняємо двома їй еквівалентними системами з трикутними матрицями (1.19):
Зворотній хід. Записуємо системи (1.21) у розгорнутому виді:
Використовуючи (1.22) і (1.23) послідовно знаходимо:
1.1.6 Метод прогону
Для розвязку систем виду або, що те ж саме,
використовується метод прогону, заснований на припущенні, що шукані невідомі звязані рекурентним співвідношенням:
де i=1,…,n-1.
Використовуючи це співвідношення, виразимо xi-1 і xi через xi+1 і підставимо в рівняння (1.26):
де Fi - права частина i-го рівняння. Це співвідношення буде виконуватися незалежно від розвязку, якщо зажадати
Звідси випливає:
З першого рівняння одержимо:
Після знаходження прогоночних коефіцієнтів ? і ?, використовуючи рівняння (1.27), одержимо розвязок системи. При цьому
1.1.7 Матричний метод
Матричний метод розвязку систем лінійних алгебраїчних рівнянь з ненульовим визначником полягає в наступному.
Нехай дана система n лінійних рівнянь з n невідомими (над довільним полем):
Тоді її можна переписати в матричній формі: , де A - основна матриця системи, і - стовпці вільних членів і рішень системи відповідно:
Помножимо це матричне рівняння ліворуч на A-1 - матрицю, зворотню до матриці A:
Оскільки A?1A=E (враховуючи асоциативність матричного добутку), одержуємо . Права частина цього рівняння дасть стовпець рішень вихідної системи. Умовою застосовності даного методу (як і взагалі існування розвязку неоднорідної системи лінійних рівнянь із числом рівнянь, рівним числу невідомих) є невирідженість матриці A. Необхідною і достатньою умовою цього є нерівність нулю визначника матриці A: det A?0.
Для однорідної системи лінійних рівнянь, тобто коли , дійсно зворотнє правило: система має нетривіальний (тобто ненульовий) розвязок тільки якщо det A=0. Такий звязок між рішеннями однорідних і неоднорідних систем лінійних рівнянь зветься альтернативою Фредгольма.
До прямих методів, що використовують властивість розрідженості матриці, можна віднести: алгоритм мінімального ступеня, алгоритм мінімального дефіциту, деревовидна блокова розбивка для асиметричного розкладання, методи вкладених або паралельних перетинів і ін.
1.2 Ітераційні методи розвязку СЛАР
Наближені методи розвязку систем лінійних рівнянь дозволяють одержувати значення коренів системи із заданою точністю у вигляді границі послідовності деяких векторів. Процес побудови такої послідовності називається ітераційним (повторюваним).
Ефективність застосування наближених методів залежить від вибору початкового вектора і швидкості збіжності процесу.
1.2.1 Метод простих ітерацій
Нехай дана система n лінійних рівнянь з n невідомими: . Припускаючи, що діагональні елементи aii?0 (i=1,...,n), виразимо x1 через перше рівняння системи, x2 - через друге рівняння і т.д. У результаті одержимо систему:
Позначимо bi /aij=?i-aij /aii=?ij, де i,j=1,...,n. Тоді система запишеться, таким чином, у матричній формі: x=?+?x. Розвяжемо цю систему методом простих ітерацій. За нульове наближення приймемо стовпець вільних членів. Кожне (k+1)-е наближення обчислюють по формулі:
Якщо послідовність наближень x(0),...,x(k) має границю , то ця границя є розвязком системи, оскільки в силу властивості границі , тобто x=?+?x [4,6].
1.2.2 Метод Зейделя
Метод Зейделя являє собою модифікацію методу простих ітерацій.
Нехай дана СЛАР, приведена до нормального вигляду: x=?+?x. Вибираємо довільно початкові наближення невідомих і підставляємо в перше рівняння системи. Отримане перше наближення підставляємо в друге рівняння системи і так далі до останнього рівняння. Аналогічно будуємо другі, треті і т.д. ітерації.
Таким чином, припускаючи, що k-і наближення відомі, методом Зейделя будуємо (k+1)-і наближення по наступних формулах:
1.2.3 Метод релаксації
Метод релаксації - наближений метод розвязку систем лінійних рівнянь.
Система лінійних рівнянь
приводиться до виду:
де
Знаходяться невязки Rj:
Вибирається початкове наближення x(0)=0. На кожному кроці необхідно перетворити на нуль максимальну невязку:
Умова зупинки: .
Відповідь знаходиться по формулі:
1.2.4 Багатосітковий метод
Багатосітковий метод - метод розвязку систем лінійних алгебраїчних рівнянь, заснований на використанні послідовності зменшуваних сіток і операторів переходу від однієї сітки до іншої. Сітки будуються на основі більших значень у матриці системи, що дозволяє використовувати цей метод при рішенні еліптичних рівнянь навіть на нерегулярних сітках.
Припустимо, що необхідно розвязати систему виду:
де A - n?n матриця з елементами aij. Для зручності зіставимо індекси з вузлами сітки, таким чином ui являє собою значення u у вузлі i. Множину вузлів сітки позначимо як ?={1,2,…,n}. Основна ідея багатосіткових методів полягає в тому, що помилка, яка не може бути усунута методами релаксації, повинна бути прибрана за допомогою корекції з розвязку на грубій сітці.
Використовуючи верхній індекс як номер рівня введемо наступні позначення:
Послідовність сіток , причому
Сіткові оператори A=A1,A2,…,AM.
Оператори інтерполяції Pk, k=1,2,…,M-1.
Оператори згладжування Sk, k=1,2,…,M-1.
Усі ці компоненти багатосіткового методу будуються на першому етапі, відомому як етап побудови.
Етап побудови.
Дорівняти k=1.
Розділити ?k на непересічні множини Ck і Fk.
Дорівняти ?k+1=Ck.
Побудувати оператор інтерполяції Pk.
Побудувати Ak+1=(Pk)TAkPk.
Побудувати якщо необхідно Sk.
Якщо сітка ?k досить мала, дорівняти M=k+1 і зупинитися. Інакше k=k+1 і перехід на крок 2.
Як тільки етап побудови завершений, може бути визначений рекурсивний цикл побудови розвязку.
Алгоритм: MGV(Ak, Pk, Sk, uk, fk).
Якщо k=M, розвязати AMuM=fM використовуючи прямий метод.
Інакше.
Застосувати метод релаксації Sk ?1 раз до Akuk=fk.
Зробити корекцію на грубій сітці.
Обчислити rk=fk?Akuk.
Обчислити rk+1=(Pk)Trk.
Застосувати MGV(Ak+1, Pk+1, Sk+1, ek+1, rk+1).
Обновити розвязок uk=uk+Pkek+1.
Застосувати метод релаксації Sk ?2 раз до Akuk=fk.
Вищенаведений алгоритм описує V(µ1, µ2) - цикл.
Вибір послідовності сіток і оператора інтерполяції є найбільш важливим елементом етапу побудови і багато в чому визначає якість багатосіткового методу. Критерієм якості є дві вимірювані величини: фактор збіжності - що показує, наскільки швидко сходиться метод, тобто яка кількість ітерацій потрібна для досягнення заданої точності; складність оператора - визначає кількість операцій і обєм памяті, необхідної для кожної ітерації методу.
Складність оператора Cop розраховується як відношення кількості ненульових елементів у всіх матрицях Ak, k=1,2,…,n до кількості ненульових елементів у матриці першого рівня A1=A.
1.2.5 Метод Ланцоша
Для розвязку СЛАР високої розмірності, матриця коефіцієнтів якої зберігається в компактному нижчеописаному вигляді, найбільш зручним ітераційним методом є метод Ланцоша [4], схема якого має вигляд:
де ai=(ri-1,ri-1)/(Asi,si), bi=(ri,ri)/(ri-1,ri-1).
Перевагою даного методу є його висока швидкість збіжності до точного розвязку. Крім того, доведено, що він має властивість «квадратичного закінчення», тобто для додатньо визначеної матриці можна гарантовано одержати точний розвязок при кількості ітерацій k?n. Розмір необхідної памяті на кожній ітерації не змінюється, тому що не вимагає перетворення матриці A. За критерій зупинки даного ітераційного процесу звичайно використовують співвідношення:
де ? - задана точність. За інший критерій збіжності іноді зручніше використовувати середньоквадратичну різницю між рішеннями, отриманими на сусідніх ітераціях:
Середньоквадратичну різницю необхідно контролювати при виконанні кожних k наперед заданих ітерацій.
Недоліком методу Ланцоша є сильна залежність збіжності від точності обчислення напрямів спуску.
Практика показує, що при рішенні СЛАР із сильно розрідженими матрицями коефіцієнтів метод Ланцоша має досить високу збіжність, причому в якості початкового наближення можна брати будь-які числа (наприклад, нулі або праву частину).
2. Схеми компактного зберігання розріджених матриць
Одним з важливих достоїнств МСЕ є те, що в результаті його застосування, як правило, виникають позитивно визначені, симетричні і добре обумовлені матриці коефіцієнтів. Крім того, вони звичайно мають таку важливу властивість, як низька щільність, що при правильній нумерації вузлів кінцево-елементної розбивки тіла приводить до групування ненульових елементів стрічкою певної ширини уздовж головної діагоналі. Це і пояснює той факт, що більшість методів розвязку СЛАР, застосовуваних разом із МСЕ, орієнтована на розвязок систем зі стрічковою структурою матриці коефіцієнтів.
Симетричність матриці коефіцієнтів дозволяє не запамятовувати майже половину її елементів, а при рішенні СЛАР зменшити обєм обчислень і, як наслідок, знизити кількість помилок округлення.
На жаль, існують цілі класи задач, для яких ширина стрічки прагне до її порядку. Наприклад, при дослідженні за допомогою МСЕ контактних задач механіки деформованого твердого тіла, у яких у взаємодії беруть участь більше двох тіл, ширина стрічки дорівнює розмірності глобальної матриці жорсткості всіх контактуючих тіл. Тому при більших порядках СЛАР ітераційні методи розвязку часто виявляються кращими.
У деяких ситуаціях, наприклад, при моделюванні динамічних процесів, доводиться розвязувати велику кількість систем рівнянь з однієї і тою же матрицею коефіцієнтів. У цьому випадку вартість прямої схеми може збігатися з вартістю розвязку трикутної системи при відомому розкладанні, оскільки вартістю єдиного розкладання, віднесеного до всіх систем, можна знехтувати. У цих ситуаціях часто вірно і те, що число ітерацій, необхідних ітераційній схемі, дуже невелике, оскільки є хороші початкові вектори.
Усі ці факти, а також простота програмної реалізації більшості ітераційних методів розвязку СЛАР і обумовлюють їхню популярність.
Застосування ітераційних методів дозволяє суттєво скоротити витрати оперативної памяті компютера за рахунок того, що при формуванні глобальної матриці жорсткості досить фіксувати тільки її ненульові елементи. У звязку з цим при використанні кінцево-елементної технології виникає проблема розробки ефективних алгоритмів формування, зберігання і використання розріджених матриць.
Память, використовувана для зберігання розріджених матриць, складається з двох частин: основної памяті, у якій утримуються числові значення елементів матриць, і додаткової памяті, де зберігаються вказівники, індекси і інша інформація, необхідна для формування структури матриць і яка забезпечує доступ до числових значень їх елементів при виконанні процедур формування і розвязку СЛАР, тобто так звані списки звязності. Способи зберігання і використання даних, що зберігаються в основній і додатковій памяті, досить різноманітні і визначаються, головним чином, обраним методом розвязку СЛАР. Для структурної побудови інформації, що зберігається в додатковій памяті - списків звязності - широко використовується теорія графів.
Особливо складні алгоритми роботи з розрідженими матрицями виникають при використанні прямих методів розвязку. У першу чергу це повязане з тим, що при розкладанні розрідженої матриці A=LLT вона звичайно зазнає деякого заповнення, тобто нижній трикутний множник L має ненульові елементи в тих позиціях, де у вихідній матриці A стояли нулі (не виключено, що при цьому можуть зявитися і нові нульові елементи). Крім того, прямі методи вимагають спеціальної нумерації вузлів кінцево-елементної сітки, яка забезпечує мінімальну ширину стрічки. Подібних проблем при використанні ітераційних методів не виникає [10,11].
У рамках даної роботи стосовно ітераційних методів розвязку СЛАР викладені 2 алгоритми зберігання і використання розріджених матриць, які відрізняються простотою реалізації і високою обчислювальною ефективністю.
2.1 Перша схема
Матриця жорсткості, що виходить при застосуванні МСЕ, має симетричну структуру, що дозволяє в загальному випадку зберігати тільки верхню трикутну частину матриці. Однак для задач із великою кількістю невідомих це також приводить до проблеми нестачі памяті. Викладений у даній роботі метод дозволяє зберігати тільки ненульові члени матриці жорсткості. Суть його полягає в наступному.
49
Спочатку, з метою виявлення звязків кожного вузла з іншими, проводиться аналіз структури дискретизації області на КЕ. Наприклад, для КЕ-сітки, зображеної на рис. 2.1, відповідна структура звязків буде мати вигляд, зображений у табл. 2.1.
Таблиця 2.1 - Структура звязків КЕ-сітки
|
№ вузла |
1 |
2 |
3 |
4 |
5 |
6 |
7 |
|
|
Звязки |
1,2,5,6,7 |
1,2,3,6 |
2,3,4,6 |
3,4,5,6,7 |
1,4,5,7 |
1,2,3,4,6,7 |
1,4,5,6,7 |
|
Тоді для зберігання матриці жорсткості необхідно построчно запамятовувати інформацію про коефіцієнти, відповідні до вузлів, з якими звязаний даний вузол. На рис. 2.2 наведені матриця жорсткості і її компактний вигляд для сітки, зображеної на рис. 2.1 [9].
Текст програми, що реалізує викладений алгоритм аналізу структури КЕ-розбивки тіла, наведений у додатку А [12,13].
Даний спосіб компактного зберігання матриці жорсткості дозволяє легко його використовувати разом з яким-небудь чисельним методом. Найбільш зручним для цієї мети здається використання вищевикладеного ітераційного методу Ланцоша, тому що на кожній ітерації потрібно тільки перемножувати матрицю коефіцієнтів СЛАР і заданий вектор. Отже, для використання запропонованого методу компактного зберігання матриці СЛАР необхідно побудувати пряме і зворотнє перетворення в первісну квадратну матрицю.
Нехай aij - елемент первісної квадратної матриці розмірністю n?n, а - її компактний вигляд. Тоді для зворотнього перетворення будуть слушні наступні співвідношення:
де m - кількість ступенів волі (m=1,2,3).
Для прямого перетворення будуть слушні співвідношення, зворотні до співвідношень (2.1) - (2.4).
2.2 Друга схема
При реалізації ітераційного розвязку системи рівнянь досить частою є ситуація, коли необхідно виконати множення матриці системи A на який-небудь вектор, наприклад, на вектор вузлових невідомих або на вектор невязки , де [20]. Побудова добутків типу або є одним з вузьких місць ефективної реалізації всього ітераційного процесу розвязку системи рівнянь, оскільки вимагає найбільших витрат процесорного часу. Нижче представлений алгоритм побудови добутку і безпосередньо повязаний з ним спосіб компактного зберігання матриці A [14-16].
Рисунок 2.3 - Кінцево-елементна модель
Розглянемо як приклад кінцево-елементну модель, зображену на рис. 2.3. З кожним вузлом сітки звязане деяке число вузлів, які можуть бути розбиті на дві множини: першу множину становлять вузли, номери яких менше i - номера розглядаємого вузла, а друга множина - вузли, номери яких більше i. У табл. 2.2 зображена структура матриці СЛАР, відповідна до вузлових звязків кінцево-елементної моделі, представленої на рис. 2.3.
Таблиця 2.2 - Структура матриці СЛАР
|
1 |
2 |
3 |
4 |
5 |
6 |
7 |
8 |
9 |
||
|
1 |
a11 |
a12 |
a13 |
0 |
a15 |
a16 |
0 |
0 |
0 |
|
|
2 |
a21 |
a22 |
a23 |
0 |
a25 |
a26 |
0 |
0 |
0 |
|
|
3 |
a31 |
a32 |
a33 |
0 |
0 |
0 |
0 |
0 |
0 |
|
|
4 |
0 |
0 |
0 |
a44 |
a45 |
a46 |
a47 |
a48 |
a49 |
|
|
5 |
a51 |
a52 |
0 |
a54 |
a55 |
a56 |
0 |
0 |
a59 |
|
|
6 |
a61 |
a62 |
0 |
a64 |
a65 |
a66 |
0 |
0 |
a69 |
|
|
7 |
0 |
0 |
0 |
a74 |
0 |
0 |
a77 |
a78 |
a79 |
|
|
8 |
0 |
0 |
0 |
a84 |
0 |
0 |
a87 |
a88 |
a89 |
|
|
9 |
0 |
0 |
0 |
a94 |
a95 |
a96 |
a97 |
a98 |
a99 |
|
Матриця A зберігається у вигляді двох векторів (рис. 2.4): - члени головної діагоналі, - ненульові елементи, розташовані над головною діагоналлю. Для формування і використання елементів вектора створюються додаткові чотири вектори, що містять інформацію про звязки вузлів кінцево-елементної моделі.
|
|
|
||||
|
1 |
a11 |
1 |
a15 |
||
|
2 |
a22 |
2 |
a16 |
||
|
3 |
a33 |
3 |
a12 |
||
|
4 |
a44 |
4 |
a13 |
||
|
5 |
a55 |
5 |
a25 |
||
|
6 |
a66 |
6 |
a26 |
||
|
7 |
a77 |
7 |
a23 |
||
|
8 |
a88 |
8 |
a47 |
||
|
9 |
a99 |
9 |
a48 |
||
|
10 |
a49 |
||||
|
11 |
a46 |
||||
|
12 |
a45 |
||||
|
13 |
a59 |
||||
|
14 |
a56 |
||||
|
15 |
a69 |
||||
|
16 |
a78 |
||||
|
17 |
a79 |
||||
|
18 |
a89 |
||||
Рисунок 2.4 - Структура векторного зберігання матриці A
Розглянемо i-е рівняння СЛАР . У загальному випадку воно має вигляд:
де aij - елементи матриці A; fi - елемент вектора правої частини ; xj - елемент вектора невідомих ; N - число невідомих.
У рівнянні (2.5) присутні суми добутків двох типів. Перша сума складається з добутків елементів aij - матриці A і елементів вектора невідомих , у яких j<i, а для добутків другої суми маємо j>i. Для побудови добутків першої суми використовуються вектори і , приклад структури яких показаний на рис. 2.5. У векторі зберігаються номери тих вузлів, які повязані з даним розглядаємим вузлом, і їх номера менше номера даного вузла. Наприклад, розглянемо вузол з номером 5 (див. рис. 2.3). Цей вузол входить в елементи з номерами 2 і 3. У силу цього він має звязок з вузлами 4, 2 і 1, номера яких менше 5. Номера цих вузлів записуються у вектор у послідовності, визначеній номерами елементів і правилом обходу їх вузлів. На рис. 2.3 стрілками показаний початок і напрямок обходу вузлів елементів. Цим пояснюється порядок появи номерів 4, 2 і 1 у векторі . У векторі зберігаються номери індексів, що координують у векторі розташування номера першого вузла всього масиву номерів вузлів, повязаних з розглядаємим вузлом. Наприклад, як видно на рис. 2.3, масив вузлів, повязаних з вузлом 5, починається з індексу 6. Останній індекс масиву визначається початковим індексом вузла, наступного за розглянутим мінус одиницю. Наприклад, для вузла 5 останній індекс масиву у векторі рівний 8=9-1, де номер 9 визначає початок масиву вузлів для вузла 6 (див. рис. 2.3).
Зовсім аналогічно будуються ще два вектори і , що містять інформацію про номери вузлів, повязаних з даним розглядаємим вузлом, але номера яких більше номера даного вузла. На рис. 2.5 зображено приклад конструкції зазначених векторів стосовно до кінцево-елементної моделі, зображеної на рис. 2.3. У векторі вказується початковий індекс, що фіксує початок списку вузлів, розміщених у векторі і повязаних з розглядаємим вузлом.
|
|
|
|
|
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
1 |
17 |
1 |
1 |
1 |
1 |
1 |
5 |
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
2 |
1 |
2 |
1 |
2 |
5 |
</ td>
Рисунок 2.5 - Структура векторів-вказівників У тих рідких випадках, коли той або інший розглядаємий вузол кінцево-елементної моделі не має звязку з вузлами, номера яких менше або більше номера даного вузла, тоді у відповідні гнізда векторів або заносяться нулі. Наприклад, вузол 4 (див. рис. 2.3) не має звязку з вузлами, номера яких менше 4, це враховується записом нуля в гніздо 4 вектора (рис. 2.5). Вузол 3 не має звязку з вузлами, номера яких більше 3 - цій обставині відповідає запис нуля в гніздо 3 вектора . Виключення становить запис в 1-е гніздо вектора , куди заноситься номер останнього гнізда вектора . Для всіх інших вузлів ці параметри визначаються тривіально. Для кожної кінцево-елементної моделі інформаційні вектори будуються один раз - відразу ж після того, як закінчена побудова масиву, що описує кінцеві елементи сітки. Структура вектора (див. рис. 2.4) у значній мірі повторює структуру вектора . Ненульові наддіагональні елементи aij матриці A зберігаються у векторі порядково, а порядок їх розміщення в межах рядка визначається порядком розміщення індексів у векторі . Тепер розглянемо побудову сум добутків, що входять у рівняння (2.5). Розгляд почнемо з другої суми, а саме: При фіксованому індексі i визначаються номери першого і останнього гнізд вектора , у яких розміщаються номери вузлів, повязаних з j-м вузлом, і для яких j>i. Наприклад, для вузла 5 (див. рис. 2.3) маємо: номер 1-го гнізда 13, а номер останнього гнізда 14=15-1. Номера першого і останнього гнізд визначають у векторі діапазон розміщення ненульових елементів матриці A, які приймають участь у добутках суми (2.6). Крім того, ці ж номера гнізд визначають у векторі діапазон розміщення номерів відповідних елементів xi вектора . Таким чином, наприклад, для вузла 5 сума (2.6) складається з двох додатків: Побудова першої суми з рівняння (2.5) алгоритмічно виконується складніше. При фіксованому номері i визначаються номери першого n і останнього m (n<m) гнізд вектора , у якому зберігаються номери j вузлів, повязаних з вузлом i, і для яких j<i. Цим самим задаються компоненти xi вектора . Тепер необхідно визначити елементи aij=aji (i?j) матриці A, що зберігаються у векторі . Ця процедура виконується в такий спосіб. Послідовно з гнізд із n по m у векторі беруться номери вузлів j (j<i). По цих номерах j з вектора визначаються номери першого k і останнього l (k<l) гнізд вектора . Для кожного вузла j у діапазоні гнізд (l-k) вектора перебуває номер i, оскільки, по-перше, вузол i повязаний з кожним вузлом j, а, по-друге, i>j. Переглядаючи гнізда вектора і порівнюючи номера вузлів, що перебувають у цих гніздах, з номером i, можна визначити номер гнізда, у якому для даного вузла j і його діапазону номерів гнізд (l-k) перебуває номер вузла i. По номеру гнізда вузла i у векторі можна знайти у векторі елемент aij=aji (i?j, j<i). Розглянемо приклад. Побудуємо суму (2.8) для вузла 5 кінцево-елементної моделі, зображеної на рис. 2.3. Із цим вузлом, як ми вже відзначали вище, звязані вузли 4, 2 і 1. Номера цих вузлів зберігаються у векторі у гніздах з 4 по 6, причому їх номера менше 5. Таким чином, компоненти вектора відомі: х4, х2 і x1. Визначимо елементи матриці A. Вузлу 4 у векторі відповідають гнізда з 8 по 12. У цих гніздах зберігаються номери вузлів, повязаних з номером 4, причому всі ці номери більше 4. Один із цих номерів - 5. Цей номер перебуває в гнізді вектора з номером 12. У гнізді з цим же номером 12 у векторі зберігається елемент a45=a54 (j=4<i=5) матриці А. Аналогічно відшукуються елементи a25 (j=2<i=5) і a15 (j=1<i=5). Таким чином, одержуємо вираз для : Остаточно з урахуванням (2.7) уся сума S5 має такий вигляд: Із представленого вище алгоритму видно, що обчислювальні витрати на побудову сум добутків і різні. Друга сума будується значно швидше. Це повязано з тим, що при побудові першої суми необхідно проводити додатковий порівняльний аналіз номерів вузлів, розміщених у векторі . При розгляді процедур алгоритму передбачалося, що всі необхідні дані розташовуються в оперативній памяті компютера. Довжина одномірних масивів, у яких розміщаються вектори , і , залежить від структури розглядаємої кінцево-елементної моделі. У табл. 2.3 представлені дані по розмірності цих векторів для трьох типів кінцево-елементних моделей, побудованих у результаті розбивки одиничного куба на 20-вузлові квадратичні елементи. Куб розбивався в наступних пропорціях: 5:5:5 (КЕМ № 1), 10:10:10 (КЕМ № 2) і 20:20:20 (КЕМ № 3). Таблиця 2.3 - Розмірність векторів КЕМ
Порівняння даного алгоритму з алгоритмами, розглянутими в роботах [17-19], показує наступне. Для реалізації запропонованого алгоритму потрібно більше оперативної памяті, оскільки для координування ненульових елементів використовуються не два масиви (вектори номерів рядків і стовпців), а чотири. Додаткові вектори і , розмірності яких не перевищують числа діагональних елементів, дозволяють локалізувати зони пошуку елементів-співмножників при побудові сум (2.6) і (2.8). За рахунок цього досягається економія часових витрат на розвязок СЛАР. На закінчення відзначу, що викладений вище алгоритм формування сум у рівнянні (2.5) може бути повністю застосований і для побудови інших аналогічних матричних добутків, необхідних для реалізації ітераційних процесів. 3 Оптимізація обчислень
Задача вирішення систем рівнянь виду досить поширена. Вона виникає при описі декількох процесів, кожний з яких описується лінійним рівнянням. Сюди відноситься вивчення процесів, що протікають у лінійних електричних схемах, апроксимація систем диференціальних рівнянь, зокрема, які описують процеси теплопровідності. Процедура розвязку великої системи рівнянь на одному процесорі приводить до значних витрат часу. Побудова багатопроцесорної обчислювальної системи не завжди є доступним і оптимальним рішенням. З підвищенням надійності і швидкодії сучасних компютерних мереж зявилася можливість розвязку подібних завдань у розподіленім обчислювальнім середовищі [21,22]. Оскільки компоненти середовища фізично рознесені в просторі, то кожний процес має локальну память і використання глобальних структур даних неможливо. Така схема взаємодії припускає обмін повідомленнями між процесами. Пропонується будувати взаємодію між процесами по наступному принципу. Є деякий процес, який одержує дані від користувача. Він проводить попередню обробку цих даних. Потім процес формує повідомлення і розсилає їх необхідному числу допоміжних процесів, які виконують безпосередні обчислення. Відповідь повертається першому процесу у вигляді повідомлень від допоміжних процесів. Тут відбувається аналіз результатів обчислень і, при необхідності, процедура повторюється. Такий принцип взаємодії процесів зветься «керуючий-робітники». Перший процес є керуючим, допоміжні процеси - це робітники [23]. Існує кілька технологій побудови розподілених обчислювальних систем. Одна з перших, технологія RPC (Remote Procedure Call), заснована на віддалених викликах процедур. Процесу дозволяється викликати процедури, реалізовані на іншому, віддаленому компютері. При цьому параметри процедури передаються на віддалений компютер, там здійснюється виконання процедури і результати вертаються в місце виклику. Технологія COM (Component Object Model) - це обєктна модель компонентів. Вона застосовується для звязку між компонентами і програмами, а також при описі API і взаємодії різних мов і середовищ програмування. Одним з головних аспектів технології COM є реалізація клієнт-серверних додатків за допомогою інтерфейсів. В COM-програму входять COM-сервер, COM-клієнт і COM-інтерфейс. COM-сервер - програма або бібліотека, яка надає деякі служби програмі-клієнтові. COM-сервер містить один або кілька COM-обєктів. Кожний COM-обєкт реалізує одну або кілька служб. Кожна служба описується інтерфейсом. COM-клієнт - програма, яка використовує служби COM-сервера, запитуючи при цьому описані інтерфейси. COM-інтерфейс - абстрактний клас, що описує властивості і методи COM-обєкта. Клієнт через інтерфейс обєкта може використовувати його методи і властивості як свої власні. Існує також технологія CORBA (Common Object Request Broker Architecture) - узагальнена архітектура брокера обєктних запитів. Дана технологія також заснована на ключовому понятті інтерфейсу. Обєкт надає службу. Клієнт звертається до неї, використовуючи інтерфейс. У технологіях COM і CORBA реалізована ідея розподілених обєктів (distributed objects). При цьому інтерфейс служить універсальним засобом взаємодії. Розподілені обєкти розміщаються на різних компютерах. Коли процес викликає метод, реалізація інтерфейсу на машині з процесом перетворює виклик методу в повідомлення, передане обєкту. Обєкт виконує запит і повертає результати. Реалізація інтерфейсу перетворює відповідне повідомлення в повертаєме значення, яке передається викликаючому процесу. Реалізація розподіленої системи, що розвязує систему рівнянь виду , можлива засобами технології COM. Перенос моделі «керуючий-робітники» на COM означає наступне. Керуючий процес повинен бути реалізований як COM-клієнт, робочі процеси, що виконують обчислення, реалізуються як COM-сервера послуги, що надають свої послуги через інтерфейси. COM-клієнт приймає вихідну систему від користувача, виконує перетворення до виду, зручного для ітерації, ділить систему на блоки, розсилає частини вихідних даних COM-серверам, використовуючи їх інтерфейси. Коли відповідь отримана, клієнт виконує збирання результатів і аналіз точності отриманого розвязку. При досягненні заданої точності процес обчислень припиняється. Важливо, що в цьому випадку COM-клієнт реалізується як одиночна програма, що використовує послуги декількох віддалених COM-серверів. У свою чергу, COM-сервера містять COM-обєкти, які і реалізують функціональність, що забезпечує розвязок систем виду ітераційними методами. Слід зазначити, що при такій реалізації розвязок окремих блоків рівнянь вихідної системи здійснюється на різних хостах розподіленої системи. Це дозволяє робити обчислення паралельно і незалежно, що значно скорочує час розвязку в порівнянні з обчисленнями на однопроцесорній системі. 4. Чисельні експерименти
Була розвязана задача з кількістю невідомих, що перевищує чотири мільйона та проведене порівняння ефективності прямих і ітераційних методів. Конфігурація тестового компютера наведена в табл. 4.1. Таблиця 4.1 - Конфігурація тестового компютера
4.1 Пружне деформування тонкостінної просторової рами
Була розглянута задача про пружне деформування тонкостінної просторової рами (рис. 4.1) при наступних геометричних параметрах: L=H=2м, L1=1,98м, L2=0,01м. Пружні властивості матеріалу рами: E=106Па, ?=0,3. Зовнішнє навантаження P=102Н. Кінцево-елементна модель усієї конструкції містила 1449604 вузла і 4374320 скінчених елементів у формі тетраедра. СЛАР розвязувалася з використанням першої схеми. Час розвязку склав близько двох годин. 4.2 Контактна взаємодія оболонкової конструкції і ложемента
Була вирішена задача про контактну взаємодію оболонкової конструкції і ложемента (рис. 4.2). Дана задача часто виникає на практиці. При транспортуванні або зберіганні з горизонтальним розташуванням осі оболонкові конструкції встановлюються на кругові опори - ложементи. Взаємодія підкріплених оболонкових конструкцій і ложементів здійснюється через опорні шпангоути, довжина яких уздовж вісі оболонки порівняна з шириною ложементів і багато менше радіуса оболонки і величини зони контакту. Дискретна модель ложемента (у тривимірній постановці) представлена на рис. 4.3. При побудові даної КЕ-моделі було використано 6689 вузлів і 15323 КЕ у формі тетраедра. Розмір нестиснутої матриці жорсткості для такої задачі становить (6689*3)2*8=3221475912 байт, що приблизно дорівнює 3 ГБ оперативної памяті. Задача розвязувалася двома способами - методом Гауса і методом Ланцоша. Порівняння результатів розвязку наведено в табл. 4.2. Таблиця 4.2 - Результати розвязку
З табл. 4.2 легко бачити, що час розвязку методом Ланцоша значно менше, ніж у випадку використання методу Гауса. Висновки
У даній роботі розглянуті способи компактного зберігання матриці коефіцієнтів системи лінійних алгебраїчних рівнянь і методи розвязку СЛАР. Розглянуті алгоритми компактного зберігання матриці жорсткості дозволяють значно скоротити обєм необхідної памяті для зберігання матриці жорсткості. Класичні методи зберігання, що враховують симетричну і стрічкову структуру матриць жорсткості, що виникають при застосуванні методу скінчених елементів, як правило, не застосовні при рішенні контактних задач, тому що при їхньому рішенні матриці жорсткості декількох тіл поєднуються в одну загальну матрицю, ширина стрічки якої може прагнути до розмірності системи. Викладені у роботі методи компактного зберігання матриці коефіцієнтів СЛАР дозволяють на прикладі розвязку контактних задач досягти значної економії оперативної памяті. Це дозволяє розвязувати системи з мільйонами невідомих на сучасних персональних компютерах. Вважаю, що для оптимальної продуктивності при розвязку СЛАР високих ступенів треба використовувати 24 ГБ оперативної памяті стандарту DDRIII 1333, мікропроцесор Intel Core i7-965 Extreme Edition, системну плату ASUS Rampage II Extreme, дві відеокарти AMD Radeon HD 4870X2 та твердотільний накопичувач OCZ Summit. Потужність такої конфігурації може скласти близько 5 терафлопс, що досить непогано. Для пришвидшення розрахунків вважаю доцільним розгін мікропроцесора, використання у програмі спеціалізованих наборів інструкцій процесорів по роботі з числами з плаваючою крапкою та використання у програмі ресурсів відеопроцесора відеокарти, за умови оптимізації програми під такі нововведення. Також значно пришвидшити розрахунки допоможе використання більшої кількості процесорів (або процесорних ядер). Наведені алгоритми розвязку СЛАР високих ступенів у задачах механіки можуть ефективно застосовуватися для аналізу тривимірних задач високої розмірності, що особливо актуально в задачах сучасного машинобудування. Перелік посилань
1. Зенкевич О., Морган К. Конечные методы и аппроксимация. - М.: Мир, 1980. - 479 с. 2. Зенкевич О. Метод конечных элементов. - М.: Мир, 1975. - 217 с. 3. Стрэнг Г., Фикс Дж. Теория метода конечных элементов. - М.: Мир, 1977. - 346 с. 4. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. - М.: Наука, 1987. - 632 с. 5. Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления. - М.: Наука, 1984. - 176 с. 6. Бахвалов Н.С. Численные методы. - М.: Наука, 1975. - 153 с. 7. Годунов С.К. Решение систем линейных уравнений. - Новосибирск: Наука, 1980. - 81 с. 8. Gustavson F.G. Some basic techniques for solving sparse matrix algorithms. - New York: Plenum Press, 1972. - 215 p. 9. Джордж А., Лиу Дж. Численное решение больших разрежённых систем уравнений. - М.: Мир, 1984. - 333 с. 10. Rose D.J. A graph theoretic study of the numerical solution of sparse positive definite system of linear equations. - New York: Academic Press, 1972. - 473 p. 11. Мосаковский В.И., Гудрамович В.С., Макеев Е.М. Контактные задачи теории оболочек и стержней. - М.: Машиностроение, 1978. - 718 с. 12. Рвачев В.Л., Шевченко А.Н. Проблемно-ориентированные языки и системы для инженерных расчётов. - К.: Техніка, 1988. - 198 с. 13. Голуб Дж., Ван Лоун Ч. Матричные вычисления. - М.: Мир, 1999. - 548 с. 14. Ортега Дж., Пул У. Введение в численные методы решения дифференциальных уравнений. - М.: Наука, 1986. - 288 с. 15. Тьюарсон Р. Разрежённые матрицы. - М.: Мир, 1977. - 191 с. 16. Брамеллер А., Аллан Р., Хэмэм Я. Слабозаполненные матрицы. - М.: Энергия, 1979. - 192 с. 17. Эстербю О., Златев 3. Прямые методы для разрежённых матриц. - М.: Мир, 1987. - 120 с. 18. Писсанецки С. Технология разреженных матриц. - М.: Мир, 1988. - 410 с. 19. Сабоннадьер Ж.-К., Кулон Ж.-Л. Метод конечных элементов и САПР. - М.: Мир, 1989. - 192 с. 20. Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. - М.: Наука, 1978. - 592 с. 21. Ортега Дж. Введение в параллельные и векторные методы решения линейных систем. - М.: Мир, 1991. - 347 с. 22. Валях Е. Последовательно-параллельные вычисления. - М.: Мир, 1985. - 216 с. 23. Таненбаум Э., Ван Стеен М. Распределённые системы. Принципы и парадигмы. - СПб.: Питер, 2003. - 481 с. Додаток А
Вихідний текст програми розвязку СЛАР з розрідженою матрицею
#include <math.h> #include <stdio.h> #include <string.h> #include <stdlib.h> #include <fstream.h> #include "matrix.h" #include "smatrix.h" #define BASE3D_4 4 #define BASE3D_8 8 #define BASE3D_10 10 const double Eps=1.0E-10; DWORD CurrentType=BASE3D_10; class RVector { private: Vector<double> Buffer; public: RVector() {} ~RVector() {} RVector(DWORD Size) { Buffer.ReSize(Size); } RVector(RVector &right) { Buffer=right.Buffer; } RVector(Vector<double> &right) { Buffer=right; } DWORD Size(void) { return Buffer.Size(); } void ReSize(DWORD Size) { Buffer.ReSize(Size); } double& operator[] (DWORD i) { return Buffer[i]; } RVector& operator=(RVector &right) { Buffer=right.Buffer; return * this; } RVector& operator=(Vector<double> &right) { Buffer=right; return * this; } void Sub(RVector&); void Sub(RVector&, double); void Mul(double); void Add(RVector&); friend double Norm(RVector&, RVector&); }; class TSMatrix { private: Vector<double> Right; Vector<double> * Array; Vector<DWORD> * Links; int Dim; DWORD Size; public: TSMatrix(void) { Size=0; Dim=0; Array=NULL; Links=NULL; } TSMatrix(Vector<DWORD> *, DWORD, int); ~TSMatrix(void) { if (Array) delete [] Array; } Vector<double> &GetRight(void) { return Right; } DWORD GetSize(void) { return Size; } int GetDim(void) { return Dim; } Vector<double> &GetVector(DWORD i) { return Array[i]; } Vector<DWORD> * GetLinks(void) { return Links; } void SetLinks(Vector<DWORD> * l) { Links=l; } void Add(Matrix<double>&,Vector<DWORD>&); void Add(DWORD I, DWORD L, DWORD J, DWORD K, double v) { DWORD Row=I, Col=L*Links[I].Size()*Dim+Find(I, J)*Dim+K; Array[Row][Col]+ v; } void Add(DWORD I, double v) { Right[I]+=v; } DWORD Find(DWORD, DWORD); void Restore(Matrix<double>&); void Set(DWORD, DWORD, double, bool); void Set(DWORD Index1, DWORD Index2, double value) { DWORD I=Index1/Dim, L=Index1%Dim, J=Index2/Dim, K=Index2%Dim, Pos=Find(I, J), Row=I, Col; if (Pos==DWORD(-1)) return; Col=L*Links[I].Size()*Dim+Find(I, J)*Dim+K; Array[Row][Col]=value; } bool Get(DWORD Index1, DWORD Index2, double &value) { DWORD I=Index1/Dim, L=Index1%Dim, J=Index2/Dim, K=Index2%Dim, Pos=Find(I, J), Row=I, Col; value=0; if (Pos==DWORD(-1)) return false; Col=L*Links[I].Size()*Dim+Find(I, J)*Dim+K; value=Array[Row][Col]; return true; } void Mul(RVector&, RVector&); double Mul(DWORD, RVector&); void write(ofstream&); void read(ifstream&); }; class RMatrix { private: Vector<double> Buffer; DWORD size; public: RMatrix(DWORD sz) { size=sz; Buffer.ReSize(size*(size+1)*0.5); } ~RMatrix() {} DWORD Size(void) { return size; } double& Get(DWORD i, DWORD j) { return Buffer[(2*size+1-i)*0.5*i+j-i]; } }; void RVector::Sub(RVector &Right) { for (DWORD i=0; i<Size(); i++) (*this)[i]-=Right[i]; } void RVector::Add(RVector &Right) { for (DWORD i=0; i<Size(); i++) (*this)[i]+=Right[i]; } void RVector::Mul(double koff) { for (DWORD i=0; i<Size(); i++) (*this)[i]*=koff; } void RVector::Sub(RVector &Right, double koff) { for (DWORD i=0; i<Size(); i++) (*this)[i]-=Right[i]*koff; } TSMatrix::TSMatrix(Vector<DWORD> * links, DWORD size, int dim) { Dim=dim; Links=links; Size=size; Right.ReSize(Dim*Size); Array=new Vector<double>[Size]; for (DWORD i=0; i<Size; i++) Array[i].ReSize(Links[i].Size()*Dim*Dim); } void TSMatrix::Add(Matrix<double> &FEMatr, Vector<DWORD> &FE) { double Res; DWORD RRow; for (DWORD i=0L; i<FE.Size(); i++) for (DWORD l=0L; l<Dim; l++) for (DWORD j=0L; j<FE.Size(); j++) for (DWORD k=0L; k<Dim; k++) { Res=FEMatr[i*Dim+l][j*Dim+k]; if (Res) Add(FE[i], l, FE[j], k, Res); } for (DWORD i=0L; i<FE.Size(); i++) for (DWORD l=0L; l<Dim; l++) { RRow=FE[INT(i%(FE.Size()))]*Dim+l; Res=FEMatr[i*Dim+l][FEMatr.Size1()]; if (Res) Add(RRow,Res); } } DWORD TSMatrix::Find(DWORD I, DWORD J) { DWORD i; for (i=0; i<Links[I].Size(); i++) if (Links[I][i]==J) return i; return DWORD(-1); } void TSMatrix::Restore(Matrix<double> &Matr) { DWORD i, j, NRow, NPoint, NLink, Pos; Matr.ReSize(Size*Dim, Size*Dim+1); for (i=0; i<Size; i++) for (j=0; j<Array[i].Size(); j++) { NRow=j/(Array[i].Size()/Dim); NPoint=(j-NRow*(Array[i].Size()/Dim))/Dim; NLink=j%Dim; Pos=Links[i][NPoint]; Matr[i*Dim+NRow][Pos*Dim+NLink]=Array[i][j]; } for (i=0; i<Right.Size(); i++) Matr[i][Matr.Size1()]=Right[i]; } void TSMatrix::Set(DWORD Index, DWORD Position, double Value, bool Case) { DWORD Row=Index, Col=Position*Links[Index].Size()*Dim+Find(Index, Index)*Dim+Position, i; double koff=Array[Row][Col], val; if (!Case) Right[Dim*Index+Position]=Value; else { Right[Index*Dim+Position]=Value*koff; for (i=0L; i<Size*Dim; i++) if (i!=Index*Dim+Position) { Set(Index*Dim+Position, i, 0); Set(i, Index*Dim+Position, 0); if (Get(i, Index*Dim+Position, val) Right[i]-=val*Value; } } } void TSMatrix::Mul(RVector &Arr, RVector &Res) { DWORD i, j, NRow, NPoint, NLink, Pos; Res.ReSize(Arr.Size()); for (i=0; i<Size; i++) for (j=0; j<Array[i].Size(); j++) { NRow=j/(Array[i].Size()/Dim); NPoint=(j-NRow*(Array[i].Size()/Dim))/Dim; NLink=j%Dim; Pos=Links[i][NPoint]; Res[i*Dim+NRow]+=Arr[Pos*Dim+NLink]*Array[i][j]; } } double TSMatrix::Mul(DWORD Index, RVector &Arr) { DWORD j, I=Index/Dim, L=Index%Dim, Start=L*(Array[I].Size()/Dim), Stop=Start+(Array[I].Size()/Dim), NRow, NPoint, NLink, Pos; double Res=0; for (j=Start; j<Stop; j++) { NRow=j/(Array[I].Size()/Dim); NPoint=(j-NRow*(Array[I].Size()/Dim))/Dim; NLink=j%Dim; Pos=Links[I][NPoint]; Res+Arr[Pos*Dim+NLink]*Array[I][j]; } return Res; } void TSMatrix::write(ofstream& Out) { DWORD ColSize; Out.write((char *)&(Dim), sizeof(DWORD)); Out.write((char *)&(Size), sizeof(DWORD)); for (DWORD i=0; i<Size; i++) { ColSize=Array[i].Size(); Out.write((char *)&(ColSize), sizeof(DWORD)); for (DWORD j=0; j<ColSize; j++) Out.write((char *)&(Array[i][j]), sizeof(double)); } for (DWORD i=0; i<Size*Dim; i++) Out.write((char *)&(Right[i]), sizeof(double)); } void TSMatrix::read(ifstream& In) { DWORD ColSize; In.read((char *)&(Dim), sizeof(DWORD)); In.read((char *)&(Size), sizeof(DWORD)); if (Array) delete [] Array; Array=new Vector<double>[Size]; Right.ReSize(Size*Dim); for (DWORD i=0; i<Size; i++) { In.read((char *)&(ColSize), sizeof(DWORD)); Array[i].ReSize(ColSize); for (DWORD j=0; j<ColSize; j++) In.read((char *)&(Array[i][j]), sizeof(double)); } for (DWORD i=0; i<Size*Dim; i++) In.read((char *)&(Right[i]), sizeof(double)); } bool Output(char * fname, Vector<double> &x, Vector<double> &y, Vector<double> &z, Matrix<DWORD> &tr, DWORD n, DWORD NumNewPoints, DWORD ntr, Matrix<DWORD> &Bounds, DWORD CountBn) { char * Label="NTRout"; int type=CurrentType, type1=(type==BASE3D_4 || type==BASE3D_10) ? 3 : 4; DWORD NewSize, i, j; ofstream Out; if (type==BASE3D_4) type1=3; else if (type==BASE3D_8) type1=4; else if (type==BASE3D_10) type1=6; Out.open(fname, ios::out | ios::binary); if (Out.bad()) return true; Out.write((const char *) Label, 6*sizeof(char)); if (Out.fail()) return true; Out.write((const char *) &type, sizeof(int)); if (Out.fail()) return true; Out.write((const char *) &CountBn, sizeof(DWORD)); if (Out.fail()) { Out.close(); return true; } Out.write((const char *) &(NewSize=n+NumNewPoints), sizeof(DWORD)); if (Out.fail()) return true; Out.write((const char *) &(NumNewPoints), sizeof(DWORD)); if (Out.fail()) { Out.close(); return true; } for (DWORD i=0; i<n; i++) { Out.write((const char *) &x[i], sizeof(double)); Out.write((const char *) &y[i], sizeof(double)); Out.write((const char *) &z[i], sizeof(double)); if (Out.fail()) { Out.close(); return true; } } for (i=0; i<NumNewPoints; i++) { Out.write((const char *) &x[n+i], sizeof(double)); Out.write((const char *) &y[n+i], sizeof(double)); if (Out.fail()) { Out.close(); return true; } } Out.write((const char *) &(ntr), sizeof(DWORD)); if (Out.fail()) { Out.close(); return true; } for (i=0; i<ntr; i++) for (j=0; j<(DWORD) type; j++) { DWORD out=tr[i][j]; Out.write((const char *) &out, sizeof(DWORD)); if (Out.fail()) { Out.close(); return true; } } for (i=0; i<CountBn; i++) for (j=0; j<(DWORD) type1; j++) { DWORD out=Bounds[i][j]; Out.write((const char *) &out, sizeof(DWORD)); if (Out.fail()) { Out.close(); return true; } } } bool Test(DWORD * a, DWORD * b) { bool result; int NumPoints=3; if (CurrentType==BASE3D_8) NumPoints=4; else if (CurrentType==BASE3D_10) NumPoints=6; for (int i=0; i<NumPoints; i++) { result=false; for (int j=0; j<NumPoints; j++) if (b[j]==a[i]) { result=true; break; } if (result==false) return false; } return true; } void Convert(Vector<double> &X, Vector<double> &Y, Vector<double> &Z, Matrix<DWORD> &FE, DWORD NumTr, Matrix<DWORD> &Bounds, DWORD &BnCount) { int cData8[6][5]={{0, 4, 5, 1, 7}, {6, 2, 3, 7, 0}, {4, 6, 7, 5, 0}, {2, 0, 1, 3, 5}, {1, 5, 7, 3, 4}, {6, 4, 0, 2, 1}}, cData4[4][4]={{0, 1, 2, 3}, {1, 3, 2, 0}, {3, 0, 2, 1}, {0, 3, 1, 2}}, cData10[4][7]={{0, 1, 2, 4, 5, 6, 3}, {0, 1, 3, 4, 8, 7, 2}, {1, 3, 2, 8, 9, 5, 0}, {0, 2, 3, 6, 9, 7, 1}}, cData[6][7], Data[6], l, Num1, Num2, m; DWORD i, j, p[6], pp[6], Index; Matrix<DWORD> BoundList(4*NumTr, 6); double cx, cy, cz, x1, y1, z1, x2, y2, z2, x3, y3, z3; Bounds.ReSize(4*NumTr, 6); switch (CurrentType) { case BASE3D_4: Num1=4; Num2=3; for (l=0; l<Num1; l++) for (m=0; m<=Num2; m++) cData[l][m]=cData4[l][m]; break; case BASE3D_8: Num1=6; Num2=4; for (l=0; l<Num1; l++) for (m=0; m<=Num2; m++) cData[l][m]=cData8[l][m]; break; case BASE3D_10: Num1=4; Num2=6; for (l=0; l<Num1; l++) for (m=0; m<=Num2; m++) cData[l][m]=cData10[l][m]; } printf("Create bounds...r"); for (i=0; i<NumTr-1; i++) for (int j=0; j<Num1; j++) if (!BoundList[i][j]) { for (l=0; l<Num2; l++) p[l]=FE[i][cData[j][l]]; for (DWORD k=i+1; k<NumTr; k++) for (int m=0; m<Num1; m++) if (!BoundList[k][m]) { for (int l=0; l<Num2; l++) pp[l]=FE[k][cData[m][l]]; if (Test(p, pp)) BoundList[i][j]=BoundList[k][m]=1; } } x1=X[FE[i][cData[j][0]]]; y1=Y[FE[i][cData[j][0]]]; z1=Z[FE[i][cData[j][0]]]; x2=X[FE[i][cData[j][1]]]; y2=Y[FE[i][cData[j][1]]]; z2=Z[FE[i][cData[j][1]]]; x3=X[FE[i][cData[j][2]]]; y3=Y[FE[i][cData[j][2]]]; z3=Z[FE[i][cData[j][2]]]; for (l=0; l<Num2; l++) Data[l]=cData[j][l]; if (((cx-x1)*(y2-y1)*(z3-z1)+(cy-y1)*(z2-z1)*(x3-x1)+(y3-y1)*(cz-z1)*(x2-x1)-(x3-x1)*(y2-y1)*(cz-z1)-(y3-y1)*(z2-z1)*(cx-x1)-(cy-y1)*(z3-z1)*(x2-x1))>0) { if (CurrentType==BASE3D_4) { Data[0]=cData[j][0]; Data[1]=cData[j][2]; Data[2]=cData[j][1]; } else if (CurrentType==BASE3D_10) { Data[0]=cData[j][0]; Data[1]=cData[j][2]; Data[2]=cData[j][1]; Data[3]=cData[j][5]; Data[5]=cData[j][3]; } else { Data[0]=cData[j][0]; Data[1]=cData[j][3]; Data[2]=cData[j][2]; Data[3]=cData[j][1]; } } for (l=0; l<Num2; l++) Bounds[BnCount][l]=FE[i][Data[l]]; BnCount++; } void main(int argc,char * argv) { char * input1, * input2, * input3, * op="", * sw; if (argc<5 || argc>6) return; sw=argv[1]; input1=argv[2]; input2=argv[3]; input3=argv[4]; if (!strcmp(sw, "-t10")) CurrentType=BASE3D_10; else if (!strcmp(sw, "-c8")) CurrentType=BASE3D_8; else { printf("Unknown switch %snn", sw); return; } if (argc==6) { op=argv[5]; if (strcmp(op, "/8") && strcmp(op, "/6")) { printf("Unknown options %snn", op); return; } } if (CreateFile(input1, input2, input3, op)) printf("OKn"); } bool CreateFile(char * fn1, char * fn2, char * fn3, char * Op) { FILE * in1, * in2, * in3; Vector<double> X(1000), Y(1000), Z(1000); DWORD NumPoints, NumFE, NumBounds=0, tmp; Matrix<DWORD> FE(1000, 10), Bounds; if ((in1=fopen(fn1, "r"))==NULL) { printf("Unable open file %s", fn1); return false; } if ((in2=fopen(fn2, "r"))==NULL) { printf("Unable open file %s", fn2); return false; } if (CurrentType==BASE3D_10) { if (!ReadTetraedrData(fn1, fn2, in1, in2, X, Y, Z, FE, NumPoints, NumFE)) return false; if (!strcmp(Op, "/8")) Convert1024(FE, NumFE); Convert(X, Y, Z, FE, NumFE, Bounds, NumBounds); return !Output(fn3, X, Y, Z, FE, NumPoints, 0, NumFE, Bounds, NumBounds); } if (CurrentType==BASE3D_8) { if (!ReadCubeData(fn1, fn2, in1, in2, X, Y, Z, FE, NumPoints, NumFE)) return false; if (!strcmp(Op, "/6")) Convert824(FE, NumFE); Convert(X, Y, Z, FE, NumFE, Bounds, NumBounds); return !Output(fn3, X, Y, Z, FE, NumPoints, 0, NumFE, Bounds, NumBounds); } return false; } void Convert824(Matrix<DWORD> &FE, DWORD &NumFE) { Matrix<DWORD> nFE(6*NumFE, 4); DWORD data[][4]={{0, 2, 3, 6}, {4, 5, 1, 7}, {0, 4, 1, 3}, {6, 7, 3, 4}, {1, 3, 7, 4}, {0, 4, 6, 3}}, Current=0; for (DWORD i=0; i<NumFE; i++) for (DWORD j=0; j<6; j++) { for (DWORD k=0; k<4; k++) nFE[Current][k]=FE[i][data[j][k]]; Current++; } CurrentType=BASE3D_4; NumFE=Current; FE=nFE; } void Convert1024(Matrix<DWORD> &FE, DWORD &NumFE) { Matrix<DWORD> nFE(8*NumFE, 4); DWORD data[][4]={{3, 7, 8, 9}, {0, 4, 7, 6}, {2, 5, 9, 6}, {7, 9, 8, 6}, {4, 8, 5, 1}, {4, 5, 8, 6}, {7, 8, 4, 6}, {8, 9, 5, 6}}, Current=0; for (DWORD i=0; i<NumFE; i++) for (DWORD j=0; j<8; j++) { for (DWORD k=0; k<4; k++) nFE[Current][k]=FE[i][data[j][k]]; Current++; } CurrentType=BASE3D_4; NumFE=Current; FE=nFE; } bool ReadTetraedrData(char * fn1, char * fn2, FILE * in1, FILE * in2, Vector<double> &X, Vector<double> &Y, Vector<double> &Z, Matrix<DWORD> &FE, DWORD &NumPoints, DWORD &NumFE) { double tx, ty, tz; char TextBuffer[1001]; DWORD Current=0L, tmp; while (!feof(in1)) { if (fgets(TextBuffer, 1000, in1)==NULL) { if (feof(in1)) break; printf("Unable read file %s", fn1); fclose(in1); fclose(in2); return false; } if (sscanf(TextBuffer, "%ld %lf %lf %lf", &NumPoints, &tx, &ty, &tz)!=4) continue; X[Current]=tx; Y[Current]=ty; Z[Current]=tz; if (++Current==999) { Vector<double> t1=X, t2=Y, t3=Z; X.ReSize(2*X.Size()); Y.ReSize(2*X.Size()); Z.ReSize(2*X.Size()); for (DWORD i=0; i<Current; i++) { X[i]=t1[i]; Y[i]=t2[i]; Z[i]=t3[i]; } } if (Current%100==0) printf("Line: %ldr", Current); } fclose(in1); NumPoints=Current; Current=0L; while (!feof(in2)) { if (fgets(TextBuffer, 1000, in2)==NULL) { if (feof(in2)) break; printf("Unable read file %s", fn2); fclose(in2); return false; } if (sscanf(TextBuffer, "%d %d %d %d %d %ld %ld %ld %ld %ld %ld %ld %ld", &tmp, &tmp, &tmp, &tmp, &tmp, &FE[Current][0], &FE[Current][1], &FE[Current][2], &FE[Current][3], &FE[Current][4], &FE[Current][5], &FE[Current][6], &FE[Current][7])!=13) continue; if (fgets(TextBuffer, 1000, in2)==NULL) { printf("Unable read file %s", fn2); fclose(in2); return false; } if (sscanf(TextBuffer, "%ld %ld", &FE[Current][8], &FE[Current][9])!=2) { printf("Unable read file %s", fn2); fclose(in2); return false; } } } bool ReadCubeData(char * fn1, char * fn2, FILE * in1, FILE * in2, Vector<double> &X, Vector<double> &Y, Vector<double> &Z, Matrix<DWORD> &FE, DWORD &NumPoints, DWORD &NumFE) { double tx, ty, tz; char TextBuffer[1001]; DWORD Current=0L, tmp; while (!feof(in1)) { if (fgets(TextBuffer, 1000, in1)==NULL) { if (feof(in1)) break; printf("Unable read file %s", fn1); fclose(in1); fclose(in2); return false; } if (sscanf(TextBuffer, "%ld %lf %lf %lf", &NumPoints, &tx, &ty, &tz)!=4) continue; X[Current]=tx; Y[Current]=ty; Z[Current]=tz; if (++Current==999) { Vector<double> t1=X, t2=Y, t3=Z; X.ReSize(2*X.Size()); Y.ReSize(2*X.Size()); Z.ReSize(2*X.Size()); for (DWORD i=0; i<Current; i++) { X[i]=t1[i]; Y[i]=t2[i]; Z[i]=t3[i]; } } if (Current%100==0) printf("Line: %ldr", Current); } fclose(in1); NumPoints=Current; Current=0L; while (!feof(in2)) { if (fgets(TextBuffer, 1000, in2)==NULL) { if (feof(in2)) break; printf("Unable read file %s", fn2); fclose(in2); return false; } if (sscanf(TextBuffer, "%d %d %d %d %d %ld %ld %ld %ld %ld %ld %ld %ld", &tmp, &tmp, &tmp, &tmp, &tmp, &FE[Current][0], &FE[Current][1], &FE[Current][3], &FE[Current][2], &FE[Current][4], &FE[Current][5], &FE[Current][7], &FE[Current][6])!=13) continue; if (++Current==999) { Matrix<DWORD> t=FE; FE.ReSize(2*FE.Size1(), 10); for (DWORD i=0; i<Current; i++) for (DWORD j=0; j<10; j++) FE[i][j]=t[i][j]; } if (Current%100==0) printf("Line: %ldr", Current); } NumFE=Current; for (DWORD i=0; i<NumFE; i++) for (DWORD j=0; j<10; j++) FE[i][j]--; return true; } |
| ! | Как писать дипломную работу Инструкция и советы по написанию качественной дипломной работы. |
| ! | Структура дипломной работы Сколько глав должно быть в работе, что должен содержать каждый из разделов. |
| ! | Оформление дипломных работ Требования к оформлению дипломных работ по ГОСТ. Основные методические указания. |
| ! | Источники для написания Что можно использовать в качестве источника для дипломной работы, а от чего лучше отказаться. |
| ! | Скачивание бесплатных работ Подводные камни и проблемы возникающие при сдаче бесплатно скачанной и не переработанной работы. |
| ! | Особенности дипломных проектов Чем отличается дипломный проект от дипломной работы. Описание особенностей. |
| → | по экономике Для студентов экономических специальностей. |
| → | по праву Для студентов юридических специальностей. |
| → | по педагогике Для студентов педагогических специальностей. |
| → | по психологии Для студентов специальностей связанных с психологией. |
| → | технических дипломов Для студентов технических специальностей. |
| → | выпускная работа бакалавра Требование к выпускной работе бакалавра. Как правило сдается на 4 курсе института. |
| → | магистерская диссертация Требования к магистерским диссертациям. Как правило сдается на 5,6 курсе обучения. |
| Дипломная работа | Формирование устных вычислительных навыков пятиклассников при изучении темы "Десятичные дроби" |
| Дипломная работа | Технологии работы социального педагога с многодетной семьей |
| Дипломная работа | Человеко-машинный интерфейс, разработка эргономичного интерфейса |
| Дипломная работа | Организация туристско-экскурсионной деятельности на т/к "Русский стиль" Солонешенского района Алтайского края |
| Дипломная работа | Разработка мероприятий по повышению эффективности коммерческой деятельности предприятия |
| Дипломная работа | Совершенствование системы аттестации персонала предприятия на примере офиса продаж ОАО "МТС" |
| Дипломная работа | Разработка системы менеджмента качества на предприятии |
| Дипломная работа | Организация учета и контроля на предприятиях жилищно-коммунального хозяйства |
| Дипломная работа | ЭКСПРЕСС-АНАЛИЗ ФИНАНСОВОГО СОСТОЯНИЯ ООО «АКТ «ФАРТОВ» |
| Дипломная работа | Психическая коммуникация |