Дипломная работа: Розробка алгоритму та програми чисельного розвязку систем лінійних алгебраїчних рівнянь з розрідженою

Реферат

Дипломна робота: 72 с., 8 рис., 5 табл., 1 додаток, 23 джерела.

Мета роботи – розробка алгоритму та програми чисельного розв'язку систем лінійних алгебраїчних рівнянь з розрідженою матрицею.

При виконанні роботи використані програми Microsoft Office Word 2007 ServicePack 1, Microsoft Office PowerPoint 2007 ServicePack 1 та Microsoft Visual Studio 2008 ServicePack 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 Контактна взаємодія оболонкової конструкції і ложемента

Висновки

Перелік посилань

Додаток А

Перелік умовних скорочень і термінів

1. Умовні скорочення

МСЕ – метод скінчених елементів.

СЛАР – система лінійних алгебраїчних рівнянь.

2. Терміни

Розмірність матриці – кількість невідомих в матриці.

Розріджена матриця – матриця, в якій більшість елементів дорівнює нулю.

Оперативна пам'ять – енергозалежний вид пам'яті комп'ютера, в якій зберігаються дані, необхідні процесору в даний час.

Крайова задача – система диференціальних рівнянь з визначеними лінійними співвідношеннями між значеннями шуканих функцій на кінцях інтервала інтегрування.

Вступ

Зараз у зв'язку з бурхливим розвитком обчислювальних засобів широке поширення одержали інформаційні технології, що мають різноманітну теоретичну і прикладну спрямованість. Серед них особливе місце займають системи автоматизованого проектування (САПР), невід'ємну частину яких становлять підсистеми математичного моделювання різних фізичних процесів. Одним із самих перспективних напрямків розвитку математичного моделювання є широке використання чисельних методів, таких як метод скінчених елементів (МСЕ). Це чисельний метод для диференціальних рівнянь, що зустрічаються у фізиці [1]. Виникнення цього методу пов'язане з розв’язком задач космічних досліджень. Уперше він був опублікований у роботі Тернера, Клужа, Мартіна і Топпа. Ця робота сприяла появі інших робіт; був опублікований ряд статей із застосуванням методу скінчених елементів до задач будівельної механіки і механіки суцільних середовищ. Важливий внесок у теоретичну розробку методу зробив в 1963 р. Мелош, який показав, що метод скінчених елементів можна розглядати як один з варіантів добре відомого методу Релєя-Рітца. У будівельній механіці метод скінчених елементів мінімізацією потенційної енергії дозволяє звести задачу до системи лінійних рівнянь рівноваги [2,3].

Однією з існуючих труднощів, що виникають при чисельній реалізації розв’язку контактних задач теорії пружності методом скінчених елементів, є розв’язоксистем лінійних алгебраїчних рівнянь (СЛАР) великої розмірності виду . За умови зберігання всіх елементів матриці A в оперативній пам'яті навіть комп'ютер з 24 ГБ оперативної пам'яті здатен вирішити СЛАР з кількістю невідомих не більше56755. Якщо ж невідомих більше, доведеться залучати суперкомп'ютер. Тому актуальною є проблема знаходження схеми компактнішого зберігання матриці СЛАР в оперативній пам'яті.

Зараз розроблена велика кількість програмних систем, що автоматизують різні аспекти чисельного моделювання проблем механіки (ANSYS, COSMOS, FORTU, МІРЕЛА, ПОЛЕ, ПРОЧНОСТЬ і ін.).

Одним з головних блоків усіх цих систем є модуль, відповідальний за розв’язок СЛАР. Слід зазначити, що звичайно розв’язок СЛАР займає більшість часу розрахунку задачі на комп'ютері. Крім того, від цього етапу прямо залежить якість і точність розв’язку задачі в цілому. Тому від вибору методу розв’язку СЛАР і алгоритму його реалізації багато в чому залежить подальший успіх чисельного розрахунку задачі МСЕ.

Більшість існуючих методів розв’язку СЛАР в МСЕ розроблені в припущенні того, що матриця A має стрічкову структуру, причому ширина стрічки m<<n, де n×n– розмірність A. Однак, при використанні МСЕ для чисельного розв’язку контактних задач можливі випадки, коли ширина стрічки m→n.

1. Огляд методів розв’язку СЛАР, що виникають у МСЕ

У процесі побудови дискретних аналогів крайових задач виникають великі системи в загальному випадку нелінійних алгебраїчних рівнянь, які вирішуються у два етапи: на першому етапі вони линеаризуються, а потім отримана система лінійних рівнянь вирішується за допомогою якого-небудь методу лінійної алгебри. Якщо збіжність не досягнута, то процес повторюється.

Основна ідея методу скінчених елементів полягає в тому, що будь-яку неперервну величину, таку як температура, тиск і переміщення, можна апроксимувати дискретною моделлю, яка будується на множині кусочно-неперервних функцій, визначених на кінцевім числі підобластей. Кусочно-неперервні функції визначаються за допомогою значень неперервної величини в кінцевім числі точок розглядаємої області [1,2,3].

У загальному випадку, неперервна величина заздалегідь невідома, і потрібно визначити значення цієї величини в деяких внутрішніх точках області. Дискретну модель, однак, дуже легко побудувати, якщо спочатку припустити, що числові значення цієї величини в кожній внутрішній точці області відомі. Після цього можна перейти до загального випадку. Отже, при побудові конкретної моделі неперервної величини роблять наступним чином:

а) у розглядаємій області фіксується кінцеве число точок. Ці точки називаються вузловими точками або просто вузлами;

б) значення неперервної величини в кожній вузловій точці вважається змінною, яка повинна бути визначена;

в) область визначення неперервної величини розбивається на кінцеве число підобластей, називаних елементами. Ці елементи мають загальні вузлові точки і у сукупності апроксимують форму області;

г) неперервна величина апроксимується на кожному елементі функцією, яка визначається за допомогою вузлових значень цієї величини. Для кожного елемента визначається своя функція, але функції підбираються таким чином, щоб зберігалася неперервність величини уздовж границь елемента.

Способи розв’язку СЛАР, застосовувані разом із МСЕ, можна умовно розбити на дві групи: прямі (точні) і ітераційні (наближені). Прийнято вважати, що прямі методи найбільш ефективні при рішенні СЛАР невисоких порядків (не більше 10⁵ невідомих), а ітераційні – для розв’язку СЛАР великих і надвеликих систем (понад 10⁵ невідомих).

Якщо всю матрицю коефіцієнтів СЛАР вдається розмістити в пам'яті, то розв’язок такої системи не представляє якої-небудь складності. Однак на практиці часто доводиться мати справу з матрицями, розміри яких суттєво перевищують об'єм оперативної пам'яті комп'ютера. У таких випадках і доводиться використовувати різні модифікації як прямих, так і ітераційних методів.

Остаточне рішення про застосування прямих або ітераційних методів розв’язку СЛАР необхідно ухвалювати на основі аналізу структури досліджуваної математичної задачі. Прямі методи розв’язку СЛАР вигідніше використовувати, якщо необхідно розв’язувати багато однакових систем з різними правими частинами, або якщо матриця А не є додатньо-визначеною. Крім того, існують задачі з такою структурою матриці, для якої прямі методи завжди кращі, ніж ітераційні.

1.1 Точні методи розв’язку СЛАР

Розглянемо ряд точних методів розв’язку СЛАР [4,5].

1.1.1 Метод Гауса

Нехай дана система , де A – матриця розмірності m×m (квадратна). У припущенні, що a₁₁ ≠0, перше рівняння системи

ділимо на коефіцієнт a₁₁ , у результаті одержуємо рівняння:

Потім з кожного з інших рівнянь віднімається перше рівняння, помножене на відповідний коефіцієнт a_i1 . У результаті ці рівняння перетворяться до виду:

Перше невідоме виявилося виключеним із усіх рівнянь, крім першого. Далі в припущенні, що , ділимо друге рівняння на коефіцієнт і виключаємо невідоме з усіх рівнянь, починаючи з третього і т.д. У результаті послідовного виключення невідомих, матриця системи рівнянь стане трикутною:

Сукупність проведених обчислень називається прямим ходом методу Гауса.

З m -го рівняння системи визначаємо x_m , з (m-1) -го рівняння визначаємо x_m-1 і т.д. до x₁ . Сукупність таких обчислень називають зворотнім ходом методу Гауса.

Реалізація прямого методу Гауса вимагає N~2m ² /3 арифметичних операцій, а зворотнього – N~m ² арифметичних операцій.

1.1.2 Метод Крамера

Нехай дана система лінійних рівнянь, в якій число рівнянь дорівнює числу невідомих:

Припустимо, що визначник системи d ≠0. Якщо тепер замінити послідовно у визначнику стовпці коефіцієнтів при невідомих x_j стовпцем вільних членів b _i , то вийдуть відповідно n визначників d₁ ,...,d_n .

Теорема Крамера. Система 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) перенесемо в праву частину всі члени з вільними невідомими x _r+1 ,...,x_n . Тоді одержимо систему, яка містить r рівнянь з r базисними невідомими. Оскільки визначник цієї системи є базисний мінор M_r , то система має єдинийрозв’язок щодо базисних невідомих, який можна знайти по формулах Крамера. Даючи вільним невідомим довільні числові значення, одержимо загальнийрозв’язок вихідної системи.

Однорідна система лінійних рівнянь. Нехай дана однорідна система лінійних рівнянь з n невідомими:

Оскільки додавання стовпця з нулів не змінює рангу матриці системи, то на підставі теореми Кронекера-Kaпeллі ця система завжди сумісна і має, принаймні, нульовийрозв’язок. Якщо визначник системи відмінний від нуля і число рівнянь системи дорівнює числу невідомих, то по теоремі Крамера нульовийрозв’язок є єдиним.

У тому випадку, коли ранг матриці системи менше числа невідомих, дана система крім нульового розв’язку буде мати і ненульові розв’язки. Для знаходження цих розв’язків у системі (1.9) виділяємо r< n лінійно незалежних рівнянь, інші відкидаємо. У виділених рівняннях у лівій частині залишаємо r базисних невідомих, а інші n-r вільних невідомих переносимо в праву частину. Тоді приходимо до системи, розв’язуючи яку по формулах Крамера, виразимо r базисних невідомих x₁ ,...,х_r через n-r вільних невідомих.

1.1. 3 Метод головних елементів

Нехай дана система n лінійних рівнянь з n невідомими:

де елементи a_ij ( i, j=1,…, n) утворюють розширену матрицю системи . Виберемо найбільший по модулю і не належачий стовпцю вільних членів елемент a_pq матриці, який називається головним елементом, і обчислимо множники m_i =- a_iq / a_pq для всіх рядків з номерами i≠ p (р -й рядок, що містить головний елемент, називається головним рядком).

Далі до кожного другорядного i -го рядку додамо головний рядок, помножений на відповідний множник m_i .

У результаті одержимо нову матрицю, усі елементи q -го стовпця якої, крім a_pq , складаються з нулів.

Відкинувши цей стовпець і головний p -й рядок, одержимо нову матрицю, число рядків і стовпців якої на одиницю менше. Повторюємо ті ж операції з отриманою матрицею, після чого одержуємо нову матрицю і т.д.

Таким чином, побудуємо послідовність матриць. Для визначення невідомих x_j поєднуємо в систему всі головні рядки, починаючи з останнього.

Викладений метод розв’язку систем лінійних рівнянь називається методом головних елементів. Необхідна умова його застосування полягає в тому, що визначник матриці не дорівнює нулю [6,7].

1.1. 4 Схема Халецького

Нехай система лінійних алгебраїчних рівнянь дана в матричному вигляді: , де А – квадратна матриця розмірності n ; , – вектори-стовпці.

Представимо матрицю А у вигляді добутку нижньої трикутної матриці С і верхньої трикутної матриці В з одиничною діагоналлю, тобто А=СВ , де

Причому елементи с_ij і b_ij визначаються по формулах:

Рівняння можна записати в наступному вигляді:

Добуток матриці B на вектор-стовпець є вектором-стовпцем, який позначимо через :

Тоді рівняння (1.13) перепишемо у вигляді:

Тут елементи с_ij відомі, тому що матриця А системи вважається вже розкладеною на добуток двох трикутних матриць С і В .

Перемноживши матриці в лівій частині рівності (1.15), одержуємо систему рівнянь, з якої одержимо наступні формули для визначення невідомих:

Невідомі y_i зручно обчислювати разом з елементами b_ij .

Після того, як усі y_i визначені по формулах (1.16), підставляємо їх у рівняння (1.14) [8].

Оскільки коефіцієнти b_ij визначені в (1.12), то значення невідомих, починаючи з останнього, обчислюємо по наступних формулах:

1.1.5 Метод квадратного кореня

Даний метод використовується для розв’язку систем лінійних алгебраїчних рівнянь виду

у яких матриця А симетрична, тобто А^Т =А , a_ij =a_ji (i=j=1,...,n) .

Розв’язок системи (1.18) здійснюється у два етапи.

Прямий хід. Перетворення матриці А і представлення її у вигляді добутку двох взаємно транспонованих трикутних матриць:

де

Перемножуючи S^Т і S , і дорівнюючи матриці А , одержимо наступні формули для визначення s_ij :

Після знаходження матриці S систему (1.18) заміняємо двома їй еквівалентними системами з трикутними матрицями (1.19):

Зворотній хід. Записуємо системи (1.21) у розгорнутому виді:

Використовуючи (1.22) і (1.23) послідовно знаходимо:

1.1.6 Метод прогону

Для розв’язку систем виду або, що те ж саме,

використовується метод прогону, заснований на припущенні, що шукані невідомі зв'язані рекурентним співвідношенням:

де i=1,… ,n-1 .

Використовуючи це співвідношення, виразимо x_i-1 і x_i через x_i+1 і підставимо в рівняння (1.26):

де F_i – права частина i -го рівняння. Це співвідношення буде виконуватися незалежно від розв’язку, якщо зажадати

Звідси випливає:

З першого рівняння одержимо:

Після знаходження прогоночних коефіцієнтів α і β , використовуючи рівняння (1.27), одержимо розв’язок системи. При цьому

1.1.7 Матричний метод

Матричний метод розв’язку систем лінійних алгебраїчних рівнянь з ненульовим визначником полягає в наступному.

Нехай дана система n лінійних рівнянь з n невідомими (над довільним полем):

Тоді її можна переписати в матричній формі:, де A – основна матриця системи, і – стовпці вільних членів і рішень системи відповідно:

Помножимо це матричне рівняння ліворуч на A^-1 – матрицю, зворотню до матриці A :

Оскільки A⁻¹ A=E (враховуючи асоциативність матричного добутку), одержуємо . Права частина цього рівняння дасть стовпець рішень вихідної системи. Умовою застосовності даного методу (як і взагалі існування розв’язку неоднорідної системи лінійних рівнянь із числом рівнянь, рівним числу невідомих) є невирідженість матриці A . Необхідною і достатньою умовою цього є нерівність нулю визначника матриці A : detA ≠0.

Для однорідної системи лінійних рівнянь, тобто коли , дійсно зворотнє правило: система має нетривіальний (тобто ненульовий) розв’язок тільки якщо det A =0. Такий зв'язок між рішеннями однорідних і неоднорідних систем лінійних рівнянь зветься альтернативою Фредгольма.

До прямих методів, що використовують властивість розрідженості матриці, можна віднести: алгоритм мінімального ступеня, алгоритм мінімального дефіциту, деревовидна блокова розбивка для асиметричного розкладання, методи вкладених або паралельних перетинів і ін.

1.2 Ітераційні методи розв’язку СЛАР

Наближені методи розв’язку систем лінійних рівнянь дозволяють одержувати значення коренів системи із заданою точністю у вигляді границі послідовності деяких векторів. Процес побудови такої послідовності називається ітераційним (повторюваним).

Ефективність застосування наближених методів залежить від вибору початкового вектора і швидкості збіжності процесу.

1.2.1 Метод простих ітерацій

Нехай дана система n лінійних рівнянь з n невідомими: . Припускаючи, що діагональні елементи a_ii ≠0 (i =1,...,n ) , виразимо x ₁ через перше рівняння системи, x ₂ – через друге рівняння і т.д. У результаті одержимо систему:

Позначимо b_i / a_ij =β _i - a_ij / a_ii =α _ij , де i,j =1,..., n . Тоді система запишеться, таким чином, у матричній формі: x =β +α x . Розв'яжемо цю систему методом простих ітерацій. За нульове наближення приймемо стовпець вільних членів. Кожне (k+1) -е наближення обчислюють по формулі:

Якщо послідовність наближень x⁽⁰⁾ ,..., x^{( k)} має границю , то ця границя є розв’язком системи, оскільки в силу властивості границі , тобто x=β+α x [4,6].

1.2.2 Метод Зейделя

Метод Зейделя являє собою модифікацію методу простих ітерацій.

Нехай дана СЛАР, приведена до нормального вигляду: x=β+αx . Вибираємо довільно початкові наближення невідомих і підставляємо в перше рівняння системи. Отримане перше наближення підставляємо в друге рівняння системи і так далі до останнього рівняння. Аналогічно будуємо другі, треті і т.д. ітерації.

Таким чином, припускаючи, що k -і наближення відомі, методом Зейделя будуємо (k+1) -і наближення по наступних формулах:

1.2.3 Метод релаксації

Метод релаксації – наближений метод розв’язку систем лінійних рівнянь.

Система лінійних рівнянь

приводиться до виду:

де

Знаходяться нев'язки R_j :

Вибирається початкове наближення x ⁽⁰⁾ =0. На кожному кроці необхідно перетворити на нуль максимальну нев'язку:

Умова зупинки:.

Відповідь знаходиться по формулі:

1.2.4 Багатосітковий метод

Багатосітковий метод – метод розв’язку систем лінійних алгебраїчних рівнянь, заснований на використанні послідовності зменшуваних сіток і операторів переходу від однієї сітки до іншої. Сітки будуються на основі більших значень у матриці системи, що дозволяє використовувати цей метод при рішенні еліптичних рівнянь навіть на нерегулярних сітках.

Припустимо, що необхідно розв'язати систему виду:

де A – n×n матриця з елементами a_ij . Для зручності зіставимо індекси з вузлами сітки, таким чином u_i являє собою значення u у вузлі i . Множину вузлів сітки позначимо як Ω={1,2,…,n }. Основна ідея багатосіткових методів полягає в тому, що помилка, яка не може бути усунута методами релаксації, повинна бути прибрана за допомогою корекції з розв’язку на грубій сітці.

Використовуючи верхній індекс як номер рівня введемо наступні позначення:

Послідовність сіток , причому

Сіткові оператори A=A¹ ,A² ,…,A^M .

Оператори інтерполяції P^k , k=1,2,…,M-1.

Оператори згладжування S^k , k=1,2,…,M-1.

Усі ці компоненти багатосіткового методу будуються на першому етапі, відомому як етап побудови.

Етап побудови.

Дорівняти k=1.

Розділити Ω^k на непересічні множини C^k і F^k .

Дорівняти Ω^k+1 =C^k .

Побудувати оператор інтерполяції P^k .

Побудувати A^{k +1} =(P^k )^T A^k P^k .

Побудувати якщо необхідно S^k .

Якщо сітка Ω^k досить мала, дорівняти M=k+1 і зупинитися. Інакше k=k+1 і перехід на крок 2.

Як тільки етап побудови завершений, може бути визначений рекурсивний цикл побудови розв’язку.

Алгоритм: MGV(A^k , P^k , S^k , u^k , f^k ).

Якщо k=M, розв'язати A^M u^M =f^M використовуючи прямий метод.

Інакше.

Застосувати метод релаксації S^k μ₁ раз до A^k u^k =f^k .

Зробити корекцію на грубій сітці.

Обчислити r^k =f^k −A^k u^k .

Обчислити r^{k +1} =(P^k )^T r^k .

Застосувати MGV(A^{k +1} , P^{k +1} , S^{k +1} , e^{k +1} , r^{k +1} ).

Обновити розв’язок u^k =u^k +P^k e^k +1.

Застосувати метод релаксації S^k μ₂ раз до A^k u^k =f^k .

Вищенаведений алгоритм описує V(µ₁ , µ₂ ) – цикл.

Вибір послідовності сіток і оператора інтерполяції є найбільш важливим елементом етапу побудови і багато в чому визначає якість багатосіткового методу. Критерієм якості є дві вимірювані величини: фактор збіжності – що показує, наскільки швидко сходиться метод, тобто яка кількість ітерацій потрібна для досягнення заданої точності; складність оператора – визначає кількість операцій і об'єм пам'яті, необхідної для кожної ітерації методу.

Складність оператора C_op розраховується як відношення кількості ненульових елементів у всіх матрицях A_k , k=1,2,…,n до кількості ненульових елементів у матриці першого рівня A¹ =A.

1.2. 5 Метод Ланцоша

Для розв’язку СЛАР високоїрозмірності, матриця коефіцієнтів якої зберігається в компактному нижчеописаному вигляді, найбільш зручним ітераційним методом є метод Ланцоша [4], схема якого має вигляд:

де a_i =(r_i-1 ,r_i-1 )/(As_i ,s_i ), b_i =(r_i ,r_i )/(r_i-1 ,r_i-1 ) .

Перевагою даного методу є його висока швидкість збіжності до точного розв’язку. Крім того, доведено, що він має властивість «квадратичного закінчення», тобто для додатньо визначеної матриці можна гарантовано одержати точнийрозв’язок при кількості ітерацій k≤n . Розмір необхідної пам'яті на кожній ітерації не змінюється, тому що не вимагає перетворення матриці A . За критерій зупинки даного ітераційного процесу звичайно використовують співвідношення:

де ε – задана точність. За інший критерій збіжності іноді зручніше використовувати середньоквадратичну різницю між рішеннями, отриманими на сусідніх ітераціях:

Середньоквадратичну різницю необхідно контролювати при виконанні кожних k наперед заданих ітерацій.

Недоліком методу Ланцоша є сильна залежність збіжності від точності обчислення напрямів спуску.

Практика показує, що при рішенні СЛАР із сильно розрідженими матрицями коефіцієнтів метод Ланцоша має досить високу збіжність, причому в якості початкового наближення можна брати будь-які числа (наприклад, нулі або праву частину).

2. Схеми компактного зберігання розріджених матриць

Одним з важливих достоїнств МСЕ є те, що в результаті його застосування, як правило, виникають позитивно визначені, симетричні і добре обумовлені матриці коефіцієнтів. Крім того, вони звичайно мають таку важливу властивість, як низька щільність, що при правильній нумерації вузлів кінцево-елементної розбивки тіла приводить до групування ненульових елементів стрічкою певної ширини уздовж головної діагоналі. Це і пояснює той факт, що більшість методів розв’язку СЛАР, застосовуваних разом із МСЕ, орієнтована на розв’язок систем зі стрічковою структурою матриці коефіцієнтів.

Симетричність матриці коефіцієнтів дозволяє не запам'ятовувати майже половину її елементів, а при рішенні СЛАР зменшити об'єм обчислень і, як наслідок, знизити кількість помилок округлення.

На жаль, існують цілі класи задач, для яких ширина стрічки прагне до її порядку. Наприклад, при дослідженні за допомогою МСЕ контактних задач механіки деформованого твердого тіла, у яких у взаємодії беруть участь більше двох тіл, ширина стрічки дорівнює розмірності глобальної матриці жорсткості всіх контактуючих тіл. Тому при більших порядках СЛАР ітераційні методи розв’язку часто виявляються кращими.

У деяких ситуаціях, наприклад, при моделюванні динамічних процесів, доводиться розв’язувати велику кількість систем рівнянь з однієї і тою же матрицею коефіцієнтів. У цьому випадку вартість прямої схеми може збігатися з вартістю розв’язку трикутної системи при відомому розкладанні, оскільки вартістю єдиного розкладання, віднесеного до всіх систем, можна знехтувати. У цих ситуаціях часто вірно і те, що число ітерацій, необхідних ітераційній схемі, дуже невелике, оскільки є хороші початкові вектори.

Усі ці факти, а також простота програмної реалізації більшості ітераційних методів розв’язку СЛАР і обумовлюють їхню популярність.

Застосування ітераційних методів дозволяє суттєво скоротити витрати оперативної пам'яті комп'ютера за рахунок того, що при формуванні глобальної матриці жорсткості досить фіксувати тільки її ненульові елементи. У зв'язку з цим при використанні кінцево-елементної технології виникає проблема розробки ефективних алгоритмів формування, зберігання і використання розріджених матриць.

Пам'ять, використовувана для зберігання розріджених матриць, складається з двох частин: основної пам'яті, у якій утримуються числові значення елементів матриць, і додаткової пам'яті, де зберігаються вказівники, індекси і інша інформація, необхідна для формування структури матриць і яка забезпечує доступ до числових значень їх елементів при виконанні процедур формування і розв’язку СЛАР, тобто так звані списки зв'язності. Способи зберігання і використання даних, що зберігаються в основній і додатковій пам'яті, досить різноманітні і визначаються, головним чином, обраним методом розв’язку СЛАР. Для структурної побудови інформації, що зберігається в додатковій пам'яті – списків зв'язності – широко використовується теорія графів.

Особливо складні алгоритми роботи з розрідженими матрицями виникають при використанні прямих методів розв’язку. У першу чергу це пов'язане з тим, що при розкладанні розрідженої матриці A=LL^T вона звичайно зазнає деякого заповнення, тобто нижній трикутний множник L має ненульові елементи в тих позиціях, де у вихідній матриці A стояли нулі (не виключено, що при цьому можуть з'явитися і нові нульові елементи). Крім того, прямі методи вимагають спеціальної нумерації вузлів кінцево-елементної сітки, яка забезпечує мінімальну ширину стрічки. Подібних проблем при використанні ітераційних методів не виникає [10,11].

У рамках даної роботи стосовно ітераційних методів розв’язку СЛАР викладені 2 алгоритми зберігання і використання розріджених матриць, які відрізняються простотою реалізації і високою обчислювальною ефективністю.

2.1 Перша схема

Матриця жорсткості, що виходить при застосуванні МСЕ, має симетричну структуру, що дозволяє в загальному випадку зберігати тільки верхню трикутну частину матриці. Однак для задач із великою кількістю невідомих це також приводить до проблеми нестачі пам'яті. Викладений у даній роботі метод дозволяє зберігати тільки ненульові члени матриці жорсткості. Суть його полягає в наступному.

Спочатку, з метою виявлення зв'язків кожного вузла з іншими, проводиться аналіз структури дискретизації області на КЕ. Наприклад, для КЕ-сітки, зображеної на рис. 2.1, відповідна структура зв'язків буде мати вигляд, зображений у табл. 2.1.

Таблиця 2.1 – Структура зв'язків КЕ-сітки

Тоді для зберігання матриці жорсткості необхідно построчно запам'ятовувати інформацію про коефіцієнти, відповідні до вузлів, з якими зв'язаний даний вузол. На рис. 2.2 наведені матриця жорсткості і її компактний вигляд для сітки, зображеної на рис. 2.1 [9].

Текст програми, що реалізує викладений алгоритм аналізу структури КЕ-розбивки тіла, наведений у додатку А [12,13].

Даний спосіб компактного зберігання матриці жорсткості дозволяє легко його використовувати разом з яким-небудь чисельним методом. Найбільш зручним для цієї мети здається використання вищевикладеного ітераційного методу Ланцоша, тому що на кожній ітерації потрібно тільки перемножувати матрицю коефіцієнтів СЛАР і заданий вектор. Отже, для використання запропонованого методу компактного зберігання матриці СЛАР необхідно побудувати пряме і зворотнє перетворення в первісну квадратну матрицю.

Нехай a_ij – елемент первісної квадратної матриці розмірністю 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 – Структура матриці СЛАР

Матриця A зберігається у вигляді двох векторів (рис. 2.4): – члени головної діагоналі, – ненульові елементи, розташовані над головною діагоналлю. Для формування і використання елементів вектора створюються додаткові чотири вектори, що містять інформацію про зв'язки вузлів кінцево-елементної моделі.

1	a 11	1	a 15
2	a 22	2	a 16
3	a 33	3	a 12
4	a 44	4	a 13
5	a 55	5	a 25
6	a 66	6	a 26
7	a 77	7	a 23
8	a 88	8	a 47
9	a 99	9	a 48
10	a 49
11	a 46
12	a 45
13	a 59
14	a 56
15	a 69
16	a 78
17	a 79
18	a 89

Рисунок 2.4 – Структура векторного зберігання матриці A

Розглянемо i-е рівняння СЛАР . У загальному випадку воно має вигляд:

де a_ij – елементи матриці A ; f_i – елемент вектора правої частини; x_j – елемент вектора невідомих; N – число невідомих.

У рівнянні (2.5) присутні суми добутків двох типів. Перша сума складається з добутків елементів a_ij – матриці 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	2	6
3	2	3	2	3	0	3	2
4	0	4	4	4	8	4	3
5	4	5	2	5	13	5	5
6	7	6	1	6	15	6	6
7	11	7	4	7	1	7	3
8	12	8	5	8	18	8	7
9	13	9	2	9	8
10	1	10	9
11	4	11	6
12	7	12	5
13	7	13	9
14	8	14	6
15	4	15	9
16	6	16	8
17	5	17	9
18	9

Рисунок 2.5 – Структура векторів-вказівників

У тих рідких випадках, коли той або інший розглядаємий вузол кінцево-елементної моделі не має зв'язку з вузлами, номера яких менше або більше номера даного вузла, тоді у відповідні гнізда векторів або заносяться нулі. Наприклад, вузол 4 (див. рис. 2.3) не має зв'язку з вузлами, номера яких менше 4, це враховується записом нуля в гніздо 4 вектора (рис. 2.5). Вузол 3 не має зв'язку з вузлами, номера яких більше 3 – цій обставині відповідає запис нуля в гніздо 3 вектора . Виключення становить запис в 1-е гніздо вектора , куди заноситься номер останнього гнізда вектора . Для всіх інших вузлів ці параметри визначаються тривіально. Для кожної кінцево-елементної моделі інформаційні вектори будуються один раз – відразу ж після того, як закінчена побудова масиву, що описує кінцеві елементи сітки.

Структура вектора (див. рис. 2.4) у значній мірі повторює структуру вектора . Ненульові наддіагональні елементи a_ij матриці A зберігаються у векторі порядково, а порядок їх розміщення в межах рядка визначається порядком розміщення індексів у векторі .

Тепер розглянемо побудову сум добутків, що входять у рівняння (2.5). Розгляд почнемо з другої суми, а саме:

При фіксованому індексі i визначаються номери першого і останнього гнізд вектора , у яких розміщаються номери вузлів, пов'язаних з j -м вузлом, і для яких j>i . Наприклад, для вузла 5 (див. рис. 2.3) маємо: номер 1-го гнізда 13, а номер останнього гнізда 14=15-1. Номера першого і останнього гнізд визначають у векторі діапазон розміщення ненульових елементів матриці A , які приймають участь у добутках суми (2.6). Крім того, ці ж номера гнізд визначають у векторі діапазон розміщення номерів відповідних елементів x_i вектора . Таким чином, наприклад, для вузла 5 сума (2.6) складається з двох додатків:

Побудова першої суми з рівняння (2.5)

алгоритмічно виконується складніше. При фіксованому номері i визначаються номери першого n і останнього m (n<m) гнізд вектора , у якому зберігаються номери j вузлів, пов'язаних з вузлом i , і для яких j<i . Цим самим задаються компоненти x_i вектора . Тепер необхідно визначити елементи a_ij = a_ji (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 у векторі можна знайти у векторі елемент a_ij =a_ji (i≠j, j<i) .

Розглянемо приклад. Побудуємо суму (2.8) для вузла 5 кінцево-елементної моделі, зображеної на рис. 2.3. Із цим вузлом, як ми вже відзначали вище, зв'язані вузли 4, 2 і 1. Номера цих вузлів зберігаються у векторі у гніздах з 4 по 6, причому їх номера менше 5. Таким чином, компоненти вектора відомі: х₄ , х₂ і x₁ . Визначимо елементи матриці A . Вузлу 4 у векторі відповідають гнізда з 8 по 12. У цих гніздах зберігаються номери вузлів, пов'язаних з номером 4, причому всі ці номери більше 4. Один із цих номерів – 5. Цей номер перебуває в гнізді вектора з номером 12. У гнізді з цим же номером 12 у векторі зберігається елемент a₄₅ =a₅₄ (j=4<i=5) матриці А . Аналогічно відшукуються елементи a₂₅ (j=2<i=5) і a₁₅ (j=1<i=5) . Таким чином, одержуємо вираз для :

Остаточно з урахуванням (2.7) уся сума S₅ має такий вигляд:

Із представленого вище алгоритму видно, що обчислювальні витрати на побудову сум добутків і різні. Друга сума будується значно швидше. Це пов'язано з тим, що при побудові першої суми необхідно проводити додатковий порівняльний аналіз номерів вузлів, розміщених у векторі . При розгляді процедур алгоритму передбачалося, що всі необхідні дані розташовуються в оперативній пам'яті комп'ютера. Довжина одномірних масивів, у яких розміщаються вектори , і , залежить від структури розглядаємої кінцево-елементної моделі. У табл. 2.3 представлені дані по розмірності цих векторів для трьох типів кінцево-елементних моделей, побудованих у результаті розбивки одиничного куба на 20-вузлові квадратичні елементи. Куб розбивався в наступних пропорціях: 5:5:5 (КЕМ № 1), 10:10:10 (КЕМ № 2) і 20:20:20 (КЕМ № 3).

Таблиця 2.3 – Розмірність векторів КЕМ

Номер КЕМ	Число вузлів	Число елементів
1	756	125	16079	16071	16079
2	4961	1000	121709	121691	121709
3	35721	8000	946619	946619	946619

Порівняння даного алгоритму з алгоритмами, розглянутими в роботах [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м, L₁ =1,98м, L₂ =0,01м. Пружні властивості матеріалу рами: E=10⁶ Па, ν=0,3. Зовнішнє навантаження P=10² Н.

Кінцево-елементна модель усієї конструкції містила 1449604 вузла і 4374320 скінчених елементів у формі тетраедра. СЛАР розв’язувалася з використанням першої схеми. Час розв’язку склав близько двох годин.

4.2 Контактна взаємодія оболонкової конструкції і ложемента

Була вирішена задача про контактну взаємодію оболонкової конструкції і ложемента (рис. 4.2).

Дана задача часто виникає на практиці. При транспортуванні або зберіганні з горизонтальним розташуванням осі оболонкові конструкції встановлюються на кругові опори – ложементи. Взаємодія підкріплених оболонкових конструкцій і ложементів здійснюється через опорні шпангоути, довжина яких уздовж вісі оболонки порівняна з шириною ложементів і багато менше радіуса оболонки і величини зони контакту.

Дискретна модель ложемента (у тривимірній постановці) представлена на рис. 4.3.

При побудові даної КЕ-моделі було використано 6689вузлів і 15323 КЕ у формі тетраедра. Розмір нестиснутої матриці жорсткості для такої задачі становить (6689*3)² *8=3221475912 байт, що приблизно дорівнює 3ГБ оперативної пам'яті.

Задача розв’язуваласядвома способами – методом Гауса і методом Ланцоша. Порівняння результатів розв’язку наведено в табл. 4.2.

Таблиця 4.2 – Результати розв’язку

З табл. 4.2 легко бачити, що час розв’язку методом Ланцоша значно менше, ніж у випадку використання методу Гауса.

Висновки

У даній роботі розглянуті способи компактного зберігання матриці коефіцієнтів системи лінійних алгебраїчних рівнянь і методи розв’язку СЛАР. Розглянуті алгоритми компактного зберігання матриці жорсткості дозволяють значно скоротити об'єм необхідної пам'яті для зберігання матриці жорсткості.

Класичні методи зберігання, що враховують симетричну і стрічкову структуру матриць жорсткості, що виникають при застосуванні методу скінчених елементів, як правило, не застосовні при рішенні контактних задач, тому що при їхньому рішенні матриці жорсткості декількох тіл поєднуються в одну загальну матрицю, ширина стрічки якої може прагнути до розмірності системи.

Викладені у роботі методи компактного зберігання матриці коефіцієнтів СЛАР дозволяють на прикладі розв’язку контактних задач досягти значної економії оперативної пам'яті. Це дозволяє розв’язувати системи з мільйонами невідомих на сучасних персональних комп'ютерах.

Вважаю, що для оптимальної продуктивності при розв’язку СЛАР високих ступенів треба використовувати 24 ГБ оперативної пам'яті стандарту DDRIII 1333, мікропроцесор Intel Core i7-965 ExtremeEdition, системну плату ASUS RampageIIExtreme, дві відеокарти AMD Radeon HD 4870X2 та твердотільний накопичувачOCZSummit. Потужність такої конфігурації може скласти близько 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 %s\n\n", sw);

return;

}

if (argc==6)

{

op=argv[5];

if (strcmp(op, "/8") && strcmp(op, "/6"))

{

printf("Unknown options %s\n\n", op);

return;

}

}

if (CreateFile(input1, input2, input3, op)) printf("OK\n");

}

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: %ld\r", 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: %ld\r", 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: %ld\r", Current);

}

NumFE=Current;

for (DWORD i=0; i<NumFE; i++)

for (DWORD j=0; j<10; j++) FE[i][j]--;

return true;

}