Паралельний алгоритм множення матриць

Размещено на http://www.allbest.ru/

Зміст

Вступ

1. Постановка задачі

2. Послідовний алгоритм

3. Множення матриць при стрічковій схемі розділення даних

3.1 Визначення підзадач

3.2 Виділення інформаційних залежностей

3.3 Масштабування і розприділення під задач по процесорам

3.4 Аналіз ефективності

3.5 Результати обчислювальних експериментів

4. Алгоритм Фокса множення матриць при блоковому розділенні даних

4.1 Визначення під задач

4.2 Виділення інформаційних залежностей

4.3 Масштабування і розподіл підзадач по процесорах

4.4 Аналіз ефективності

4.5 Програмна реалізація

4.6 Результати обчислювальних експериментів

5. Алгоритм Кэннона множення матриць при блоковому розділенні даних

5.1 Визначення під задач

5.2 Виділення інформаційних залежностей

5.3 Масштабування і розподіл підзадач по процесорах

5.4 Аналіз ефективності

5.5 Результати обчислювальних експериментів

6. Огляд літератури

Висновок

Список використаних джерел



Вступ

Операція множення матриць є одним із основних завдань матричних обчислень. У даному розділі розглядаються декілька різних паралельних алгоритмів для виконання цієї операції. Два з них засновані на стрічковій схемі розділення даних. Інші два методи використовують блокову схему розділення - це широко відомі алгоритми Фокса(Fox) і Кэннона(Cannon).



1. Постановка задачі

Множення матриці A розміру і матриці B розміру приводить до отримання матриці C розміру , кожен елемент якої визначається відповідно до виразу:

(1.1).

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

(1.2)

Цей алгоритм припускає виконання операцій множення і стільки ж операцій складання елементів початкових матриць. При множенні квадратних матриць розміру кількість виконаних операцій має порядок Відомі послідовні алгоритми множення матриць, що мають меншу обчислювальну складність(наприклад, алгоритм Страссена(Strassen's algorithm)), але ці алгоритми вимагають певних зусиль для їх освоєння і, тому, у даному розділі при розробці паралельних методів за основу було взято приведений вище послідовний алгоритм. Також мається на увазі, що всі матриці надалі будуть квадратними, і матимуть розмір .



2. Послідовний алгоритм

матриця інформаційний залежність алгоритм

Послідовний алгоритм множення матриць реалізується трьома вкладеними циклами:

// Алгоритм 1.2

// Послідовний алгоритм множення матриць

double MatrixA[Size][Size];

double MatrixB[Size][Size];

double MatrixC[Size][Size];

int i,j,k;

...

for (i=0; i<Size; i++){

for (j=0; j<Size; j++){

MatrixC[i][j] = 0;

for (k=0; k<Size; k++){

MatrixC[i][j] = MatrixC[i][j] + MatrixA[i][k]*MatrixB[k][j];

}

}

}

Алгоритм 1.1. Послідовний алгоритм множення двох квадратних матриць.

Цей алгоритм є ітеративним і орієнтований на послідовне обчислення рядків матриці С. Дійсно, при виконанні однієї ітерації зовнішнього циклу(по змінній ) обчислюється один рядок результуючої матриці(див. рис. 2.1)



Рис. 2.1. На першій ітерації циклу по змінній використовується перший рядок матриці A та усі стовпці матриці B для того, щоб вичислити елементи першого рядка результуючої матриці

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

(1.3)

де це час виконання однієї елементарної скалярної операції.



3. Множення матриць при стрічковій схемі розділення даних

Розглянемо два паралельні алгоритми множення матриць, в яких матриці A і B розбиваються на безперервні послідовності рядків або стовпців(смуги).

3.1 Визначення підзадач

З визначення операції матричного множення виходить, що обчислення усіх елементів матриці С може бути виконано незалежно один від одного. Як результат, імовірний підхід для організації паралельних обчислень полягає у використанні у якості базової підзадачі процедури визначення одного елементу результуючої матриці С. Для проведення усіх необхідних обчислень кожна підзадача повинна містити по одному рядку матриці А і одному стовпцю матриці В. Загальна кількість отримуваних при такому підході підзадач виявляється рівним n2(по числу елементів матриці С).

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

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

3.2 Виділення інформаційних залежностей

Для обчислення одного рядка матриці С необхідно, щоб в кожній підзадачі містився рядок матриці А і був забезпечений доступ до усіх стовпців матриці B. Можливі способи організації паралельних обчислень полягають в наступному.

1) Перший алгоритм. Алгоритм є ітераційною процедурою, кількість ітерацій якої співпадає з числом підзадач. На кожній ітерації алгоритму кожна підзадача містить по одному рядку матриці А і одному стовпцю матриці В. При виконанні ітерації проводиться скалярне множення рядків, що містяться в підзадачах, і стовпців, що призводить до отримання відповідних елементів результуючої матриці С. Після закінчення обчислень у кінці кожної ітерації стовпці матриці В мають бути передані між підзадачами так, щоб у кожній підзадачі опинилися нові стовпці матриці В і могли бути вичислені нові елементи матриці C. При цьому така передача стовпців між підзадачами має бути організована так, щоб після завершення ітерацій алгоритму в кожній підзадачі послідовно опинилися усі стовпці матриці В.

Схема організації необхідної послідовності передач стовпців матриці В між підзадачами полягає в уявленні топології інформаційних зв'язків підзадач у вигляді кільцевої структури. В цьому випадку, на кожній ітерації підзадача , , передаватиме свій стовпець матриці В підзадачі з номером (відповідно до кільцевої структури підзадача передає свої дані підзадачі з номером 0) - див. рис. 3.1. Після виконання усіх ітерацій алгоритму необхідна умова буде забезпечена - в кожній підзадачі по черзі опиняться усі стовпці матриці В.

Рис. 3.1. Загальна схема передачі даних для першого паралельного алгоритму матричного множення при ланцюговій схемі розділення даних.

На рис. 3.1 представлені ітерації алгоритму матричного множення для випадку, коли матриці складаються з чотирьох рядків і чотирьох стовпців(n=4). На початку обчислень в кожній підзадачі , розташовуються -ий рядок матриці A та -ий стовпець матриці B. В результаті їх перемножування підзадача отримує елемент результуючої матриці С. Далі підзадачі здійснюють обмін стовпцями, в ході якого кожна підзадача передає свій стовпець матриці B наступній підзадачі відповідно до кільцевої структури інформаційних взаємодій. Виконання описаних дій повторюється до завершення усіх ітерацій паралельного алгоритму.

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

При розглянутому способі розділення даних для виконання операції матричного множення треба забезпечити послідовне отримання в підзадачах усіх рядків матриці B, поелементне множення даних і підсумовування нових значень, із раніше вичисленими результатами. Організація необхідної послідовності передач рядків матриці B між підзадачами також може бути виконана з використанням кільцевої структури інформаційних зв'язків(див. рис. 3.2).

Рис. 3.2. Загальна схема передачі даних для другого паралельного алгоритму матричного множення при ланцюговій схемі розділення даних.

На рис. 3.2 представлені ітерації алгоритму матричного множення для випадку, коли матриці складаються з чотирьох рядків і чотирьох стовпців. На початку обчислень в кожній підзадачі, розташовуються i рядки матриці A і матриці B. В результаті їх перемножування підзадача визначає -ий рядок часткових результатів шуканої матриці C. Далі підзадачі здійснюють обмін рядками, в ході якого кожна підзадача передає свій рядок матриці B наступній підзадачі відповідно до кільцевої структури інформаційних взаємодій. Далі виконання описаних дій повторюється до завершення усіх ітерацій паралельного алгоритму.

3.3 Масштабування і роз приділення під задач по процесорам

Виділені базові підзадачі характеризуються однаковою обчислювальною складністю і рівним об'ємом передаваних даних. У разі, коли розмір матриць виявляється більшим, за число процесорів , базові підзадачі можна збільшити, об'єднавши у рамках однієї підзадачі декілька сусідніх рядків і стовпців перемножуваних матриць. В цьому випадку, початкова матриця A розбивається на ряд горизонтальних смуг, а матриця B представляється у вигляді набору вертикальних(для першого алгоритму) або горизонтальних(для другого алгоритму) смуг. Розмір смуг при цьому слід вибрати рівним k=n/p(при кратно), що дозволить як і раніше забезпечити рівномірність розподілу обчислювального навантаження на процесори, які складають багатопроцесорну обчислювальну систему.

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

3.4 Аналіз ефективності

Виконаємо аналіз ефективності першого паралельного алгоритму множення матриць.

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



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

та . (1.5)

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

З урахуванням числа і тривалості виконуваних операцій час виконання обчислень паралельного алгоритму може бути оцінений таким чином:

(1.6)

(тут, як і раніше, це час виконання однієї елементарної скалярної операції).

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



(1.7)

де б - латентність, в - пропускна здатність мережі, а w - розмір елементу матриці в байтах

Відповідно до отриманих співвідношень загальний час виконання паралельного алгоритму матричного множення визначається наступним вираженням:

3.5 Результати обчислювальних експериментів

Експерименти проводилися в обчислювальному кластері на базі процесорів Intel XEON 4 EM64T, 3000 Мгц і мереж Gigabit Ethernet під управлінням операційної системи Microsoft Windows Server 2003 Standart x64 Edition(див. п. 2.3).

Для оцінки тривалості базової скалярної операції відбувалося вирішення задачі множення матриць за допомогою послідовного алгоритму і отриманий таким чином час обчислень ділилвся на загальну кількість виконаних операцій, в результаті подібних експериментів для величини було отримано значення 6.4 нс. Експерименти, виконані для визначення параметрів мережі передачі даних, показали значення латентності і пропускної здатності відповідно 130 мкс і 53.29 Мбайт/с. Усі обчислення робилися над числовими значеннями типу double, т. е. величина дорівнює 8 байт.

Результати обчислювальних експериментів приведені в таблиці 3.1. Експерименти виконувалися з використанням двох, чотирьох і восьми процесорів.



Таблиця 3.1. Результати обчислювальних експериментів по дослідженню першого паралельного алгоритму матричного множення при стрічковій схемі розподілу даних

Розмір матриць	Послідовний алгоритм	2 процесори	4 ЦП	8 ЦП
		час	Прискорення	час	Прискорення	час	Прискорення
500	0,8752	0,3758	2,3287	0,1535	5,6982	0,0968	9,0371
1000	12,8787	5,4427	2,3662	2,2628	5,6912	0,6998	18,4014
1500	43,4731	20,9503	2,0750	11,0804	3,9234	5,1766	8,3978
2000	103,0561	45,7436	2,2529	21,6001	4,7710	9,4127	10,9485
2500	201,2915	99,5097	2,0228	56,9203	3,5363	18,3303	10,9813
3000	347,8434	171,9232	2,0232	111,9642	3,1067	45,5482	7,6368

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

Порівняння експериментального часу виконання експерименту і теоретичного часу з формули(1.8) представлені в таблиці 3.2 і на рис. 3.4.



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

Розмір матриць	2 процесори	4 ЦП	8 ЦП
	час	Прискорення	час	Прискорення	час	Прискорення
500	0,8243	0,3758	0,4313	0,1535	0,2353	0,0968
1000	6,51822	5,4427	3,3349	2,2628	1,7436	0,6998
1500	21,9137	20,9503	11,1270	11,0804	5,7340	5,1766
2000	51,8429	45,7436	26,2236	21,6001	13,4144	9,4127
2500	101,1377	99,5097	51,0408	56,9203	25,9928	18,3303
3000	174,6301	171,9232	87,9946	111,9642	44,6772	45,5482

Рис. 3.4. Графік залежності теоретичного і експериментального часу виконання паралельного алгоритму від об'єму початкових даних на двох процесорах(стрічкова схема розбиття даних)



4. Алгоритм Фокса множення матриць при блоковому розділенні даних

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

4.1 Визначення підзадач

Блокова схема розбиття матриць детально розглянута в попередньому підрозділі. При такому способі розділення даних початкові матриці , і результуюча матриця представляються у вигляді наборів блоків. Для простішого викладу матеріалу матимемо на увазі, що усі матриці є квадратними розміру кількість блоків по горизонталі і вертикалі є однаковою і рівною (тобто розмір усіх блоків рівний , ). При такому представленні даних операція множення матриць А і B у блоковому виді може бути представлена так:

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

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

Для виконання усіх необхідних обчислень базовим підзадачам мають бути доступні відповідні набори рядків матриці A і стовпців матриці B. Розміщення усіх необхідних даних в кожній підзадачі неминуче приведе до дублювання і до значного зростання об'єму використовуваної пам'яті. Як результат, обчислення мають бути організовані так, щоб в кожен теперішній момент часу підзадачі містили лише частину необхідних для проведення розрахунків даних, а доступ до іншої частини даних забезпечувався б за допомогою передачі повідомлень. Один з можливих підходів - це алгоритм Фокса(Fox). Другий спосіб - алгоритм Кэннона(Cannon) - описаний в підрозділі 5.

4.2 Виділення інформаційних залежностей

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

Можливий спосіб організації обчислень за таких умов полягає в застосуванні широко відомого алгоритму Фокса(Fox), - див, наприклад, Fox et al. (1987)Kumar et al. (1994).

Відповідно до алгоритму Фокса в ході обчислень на кожній базовій підзадачі() розташовується чотири матричні блоки:

· блок матриці , що обчислюється підзадачею;

· блок матриці , що розміщується в підзадачі перед початком обчислень;

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

Виконання паралельного методу включає:

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

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

- для кожного рядка , ,блок підзадачі() пересилається в усі підзадачі того ж рядка грат; індекс , що визначає положення підзадачі в рядку, обчислюється відповідно до виразу:

де - це операція отримання залишку від цілочисельного ділення;

- отримані в результати пересилок блоки , кожної під задачі () перемножуються і додаються до блоку

;

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

Для пояснення приведених правил паралельного методу на рис. 4.1 приведений стан блоків в кожній підзадачі в ході виконання ітерацій етапу обчислень(для грат підзадач ).



Рис. 4.1. Стан блоків в кожній підзадачі в ході виконання ітерацій алгоритму Фокса.

4.3 Масштабування і розподіл підзадач по процесорах

У розглянутій схемі паралельних обчислень кількість блоків може змінюватися залежно від вибору розміру блоків. Ці розміри можуть бути підібрані так, щоб загальна кількість базових підзадач співпадала з числом процесорів . Так, наприклад, в найбільш простому випадку, коли число процесорів можна спів ставити з (тобто є повним квадратом) можна вибрати кількість блоків в матрицях по вертикалі і горизонталі рівною (тобто ;). Такий спосіб визначення кількості блоків призводить до того, що об'єм обчислень в кожній підзадачі є однаковим і, тим самим, досягається повне балансування обчислювального навантаження між процесорами. У загальному випадку при довільних кількості процесорів і розмірах матриць балансування обчислень може бути не однаковим, але за належного вибору параметроів може бути оптимізованим.

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

4.4 Аналіз ефективності

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

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

(1.9)



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

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

(1.10)

(нагадаємо, що ; є час виконання однієї елементарної скалярної операції).

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

(1.11)

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

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

(1.12)

Підсумувавши усі отримані вирази, можна отримати, що загальний час виконання алгоритму Фокса може бути визначений за допомогою наступних співвідношень:

(1.13)

(нагадаємо, що параметр визначає розмір процесорних грат і .

4.5 Програмна реалізація

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

1) Головна функція програми.

Визначає основну логіку роботи алгоритму, послідовно викликає необхідні підпрограми.

// Програма 1.1

// Алгоритм Фокса множення матриць - блокове представлення даних

// Умови виконання програми : усі матриці квадратні, розмір блоків і їх

// кількість по горизонталі і вертикалі одинаковово, процесори утворюють

// квадратні грати

int ProcNum = 0; // Number of available processes

int ProcRank = 0; // Rank of current process

int GridSize; // Size of virtual processor grid

int GridCoords[2]; // Coordinates of current processor in grid

MPI_Comm GridComm; // Grid communicator

MPI_Comm ColComm; // Column communicator

MPI_Comm RowComm; // Row communicator

void main ( int argc, char * argv[] ) {

double* pAMatrix; // The first argument of matrix multiplication

double* pBMatrix; // The second argument of matrix multiplication

double* pCMatrix; // The result matrix

int Size; // Size of matricies

int BlockSize; // Sizes of matrix blocks on current process

double *pAblock; // Initial block of matrix A on current process

double *pBblock; // Initial block of matrix B on current process

double *pCblock; // Block of result matrix C on current process

double *pMatrixAblock;

double Start, Finish, Duration;

setvbuf(stdout, 0, _IONBF, 0);

MPI_Init(&argc, &argv);

MPI_Comm_size(MPI_COMM_WORLD, &ProcNum);

MPI_Comm_rank(MPI_COMM_WORLD, &ProcRank);

GridSize = sqrt((double)ProcNum);

if (ProcNum != GridSize*GridSize) {

if (ProcRank == 0) {

printf ("Number of processes must be a perfect square \n");

}

}

else {

if (ProcRank == 0)

printf("Parallel matrix multiplication program\n");

// Creating the cartesian grid, row and column communcators

CreateGridCommunicators();

// Memory allocation and initialization of matrix elements

ProcessInitialization ( pAMatrix, pBMatrix, pCMatrix, pAblock, pBblock,

pCblock, pMatrixAblock, Size, BlockSize );

DataDistribution(pAMatrix, pBMatrix, pMatrixAblock, pBblock, Size,

BlockSize);

// Execution of Fox method

ParallelResultCalculation(pAblock, pMatrixAblock, pBblock,

pCblock, BlockSize);

ResultCollection(pCMatrix, pCblock, Size, BlockSize);

TestResult(pAMatrix, pBMatrix, pCMatrix, Size);

// Process Termination

ProcessTermination (pAMatrix, pBMatrix, pCMatrix, pAblock, pBblock,

pCblock, pMatrixAblock);

}

MPI_Finalize();}

2) Функція CreateGridCommunicators.

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

Створення грат робиться за допомогою функції MPI _ Cart _ create(вектор Periodic визначає можливість передачі повідомлень між граничними процесами рядків і стовпців створюваних грат). Після створення грат кожен процес паралельної програми матиме координати свого положення в гратах; отримання цих координат забезпечується за допомогою функції MPI _ Cart _ coords.

Формування топологій завершується створенням безлічі комунікаторів для кожного рядка і кожного стовпця грат окремо(функція MPI _ Cart _ sub).

// Creation of two-dimensional grid communicator

// and communicators for each row and each column of the grid

void CreateGridCommunicators() {

int DimSize[2]; // Number of processes in each dimension of the grid

int Periodic[2]; // =1, if the grid dimension should be periodic

int Subdims[2]; // =1, if the grid dimension should be fixed

DimSize[0] = GridSize;

DimSize[1] = GridSize;

Periodic[0] = 0;

Periodic[1] = 0;

// Creation of the Cartesian communicator

MPI_Cart_create(MPI_COMM_WORLD, 2, DimSize, Periodic, 1, &GridComm);

// Determination of the cartesian coordinates for every process

MPI_Cart_coords(GridComm, ProcRank, 2, GridCoords);

// Creating communicators for rows

Subdims[0] = 0; // Dimensionality fixing

Subdims[1] = 1; // The presence of the given dimension in the subgrid

MPI_Cart_sub(GridComm, Subdims, &RowComm);

// Creating communicators for columns

Subdims[0] = 1;

Subdims[1] = 0;

MPI_Cart_sub(GridComm, Subdims, &ColComm);

3) Функція ProcessInitialization.

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

Ця функція визначає параметри вирішуваної задачі(розміри матриць і їх блоків), виділяє пам'ять для зберігання даних і здійснює введення початкових матриць(чи формує їх за допомогою якого-небудь датчика випадкових чисел). Всього в кожному процесі має бути виділена пам'ять для зберігання чотирьох блоків - для покажчиків використовуються змінні pAblock, pBblock, pCblock, pMatrixAblock. Перші три покажчики визначають блоки матриць , і відповідно. Слід зазначити, що вміст блоків pAblock і pBblock постійно змінюється відповідно до пересилки даних між процесами, тоді як блок pMatrixAblock матриці залишається незмінним і використовується при розсилках блоків по рядках грат процесів(див. функцію AblockCommunication).

Для визначення елементів початкових матриць використовуватимемо функцію RandomDataInitialization, реалізацію якої читачеві належить виконати самостійно.

// Function for memory allocation and data initialization

void ProcessInitialization (double* &pAMatrix, double* &pBMatrix,

double* &pCMatrix, double* &pAblock, double* &pBblock, double* &pCblock,

double* &pTemporaryAblock, int &Size, int &BlockSize ) {

if (ProcRank == 0) {

do {

printf("\nEnter size of the initial objects: ");

scanf("%d", &Size);

if (Size%GridSize != 0) {

printf ("Size of matricies must be divisible by the grid size! \n");

}

}

while (Size%GridSize != 0);

}

MPI_Bcast(&Size, 1, MPI_INT, 0, MPI_COMM_WORLD);

BlockSize = Size/GridSize;

pAblock = new double [BlockSize*BlockSize];

pBblock = new double [BlockSize*BlockSize];

pCblock = new double [BlockSize*BlockSize];

pTemporaryAblock = new double [BlockSize*BlockSize];

for (int i=0; i<BlockSize*BlockSize; i++) {

pCblock[i] = 0;

}

if (ProcRank == 0) {

pAMatrix = new double [Size*Size];

pBMatrix = new double [Size*Size];

pCMatrix = new double [Size*Size];

RandomDataInitialization(pAMatrix, pBMatrix, Size);

}

}

4) Функція DataDistribution для розподілу початкових даних і функція ResultCollection для збору результатів.

Після задання початкових матриць на нульовому процесі необхідно здійснити розподіл початкових даних. Для цього призначена функція DataDistribution. Можуть бути запропоновані два способи виконання блокового розділення матриць між процесорами, організованими в двовимірні квадратні грати. З одного боку, для організації передачі блоків у рамках однієї і тієї ж комунікаційної операції можна сформувати засобами MPI похідний тип даних. З іншого боку, можна організувати двохетапну процедуру. На першому етапі матриця розділяється на горизонтальні смуги. Ці смуги розподіляються на процеси, що складають нульовий стовпець процесорних грат. Далі кожна смуга розділяється на блоки між процесами, що становлять рядки процесорних грат.

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

Реалізація функцій DataDistribution і ResultCollection є завданням для самостійної роботи.

5) Функція AblockCommunication обміну блоками матриці A.

Функція виконує розсилку блоків матриці по рядках процесорних грат. Для цього в кожному рядку грат визначається провідний процес Pivot, що здійснює розсилку. Для розсилки використовується блок pMatrixAblock, переданий в процес у момент початкового розподілу даних. Виконання операції розсилки блоків здійснюється за допомогою функції MPI _ Bcast. Слід зазначити, що ця операція є колективною і її локалізація межами окремих рядків грат забезпечується за рахунок використання комунікаторів RowComm, визначених для набору процесів кожного рядка грат окремо.

// Broadcasting matrix A blocks to process grid rows

int BlockSize) {

// Defining the leading process of the process grid row

int Pivot = (GridCoords[0] + iter) % GridSize;

// Copying the transmitted block in a separate memory buffer

if (GridCoords[1] == Pivot) {

for (int i=0; i<BlockSize*BlockSize; i++)

pAblock[i] = pMatrixAblock[i];

} // Block broadcasting

MPI_Bcast(pAblock, BlockSize*BlockSize, MPI_DOUBLE, Pivot, RowComm);

}

6) Функція перемножування матричних блоків BlockMultiplication.

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

// Множення матричних блоків

void BlockMultiplication (double *pAblock, double *pBblock,

double *pCblock, int BlockSize) {

// вычисление произведения матричных блоков

for (int i=0; i<BlockSize; i++) {

for (int j=0; j<BlockSize; j++) {

double temp = 0;

for (int k=0; k<BlockSize; k++ )

temp += pAblock [i*BlockSize + k] * pBblock [k*BlockSize + j]

pCblock [i*BlockSize + j] += temp;

}

}

}

7) Функція BblockCommunication обміну блоками матриці B.

Функція виконує циклічне зрушення блоків матриці по стовпцях процесорних грат. Кожен процес передає свій блок наступному процесу NextProc у стовпці процесів і отримує блок, переданий з попереднього процесу PrevProc в стовпці грат. Виконання операцій передачі даних здійснюється за допомогою функції MPI _ SendRecv _ replace, яка забезпечує усі необхідні пересилки блоків, використовуючи при цьому один і той же буфер пам'яті pBblock. Крім того, ця функція гарантує відсутність можливої безвиході, коли операції передачі даних починають одночасно виконуватися декількома процесами при кільцевій топології мережі.

// Cyclic shift of matrix B blocks in the process grid columns

void BblockCommunication (double *pBblock, int BlockSize) {

MPI_Status Status;

int NextProc = GridCoords[0] + 1;

if ( GridCoords[0] == GridSize-1 ) NextProc = 0;

int PrevProc = GridCoords[0] - 1;

if ( GridCoords[0] == 0 ) PrevProc = GridSize-1;

MPI_Sendrecv_replace( pBblock, BlockSize*BlockSize, MPI_DOUBLE,

NextProc, 0, PrevProc, 0, ColComm, &Status);

}

8) Функція ParallelResultCalculation.

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

void ParallelResultCalculation(double* pAblock, double* pMatrixAblock,

double* pBblock, double* pCblock, int BlockSize) {

for (int iter = 0; iter < GridSize; iter ++) {

// Sending blocks of matrix A to the process grid rows

ABlockCommunication (iter, pAblock, pMatrixAblock, BlockSize);

// Block multiplication

BlockMultiplication(pAblock, pBblock, pCblock, BlockSize);

// Cyclic shift of blocks of matrix B in process grid columns

BblockCommunication(pBblock, BlockSize);

}

}

4.6 Результати обчислювальних експериментів

Обчислювальні експерименти для оцінки ефективності паралельного алгоритму проводилися за тих же умов, що і раніше (див. п. 3.5). Результати експериментів з використанням чотирьох і дев'яти процесорів приведені в таблиці 4.1.

Таблиця 4.1 Результатів обчислювальних експериментів по дослідженню паралельного алгоритму Фокса:

Розмір матриці	Послідовний алгоритм	Паралельний алгоритм
		4 ЦП	9 ЦП
		Час	Присокрення	Час	Присокрення
500	0,8527	0,2190	3,8925	0,1468	5,8079
1000	12,8787	3,0910	4,1664	2,1565	5,9719
1500	43,4731	10,8678	4,0001	7,2502	5,9960
2000	103,0561	24,1421	4,2687	21,4157	4,8121
2500	201,2915	51,4735	3,9105	41,2159	4,8838
3000	347,8434	87,0538	3,9957	58,2022	5,9764

Рис. 4.2. Залежність прискорення від розміру матриць при виконанні паралельного алгоритму Фокса.

Порівняння часу виконання експерименту і теоретичного часу , обчислено відповідно до виразу (1.13), представлено в таблиці 4.2 і на рис. 4.3.

Таблиця 4.2. Порівняння експериментального і теоретичного часу виконання паралельного алгоритму Фокса.

Розмір матриці	Паралельний алгоритм
	4 ЦП	9 ЦП
500	0,4217	0,2190	0,2200	0,1468
1000	3,2970	3,0910	1,5924	2,1565
1500	11,0419	10,8678	5,1920	7,2502
2000	26,0726	24,1421	12,0927	21,4157
2500	50,8049	51,4735	23,3682	41,2159
3000	87,6548	87,0538	40,0923	58,2022

Рис. 4.3. Графік залежності експериментального і теоретичного часу виконання алгоритму Фокса на чотирьох ЦП.



5. Алгоритм Кэннона множення матриць при блоковому розділенні даних

Розглянемо ще один паралельний алгоритм матричного множення, грунтований на блоковому розбитті матриць.

5.1 Визначення підзадач

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

5.2 Виділення інформаційних залежностей

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

Зважаючи на висловлені зауваження етап ініціалізації алгоритму Кэннона включає виконання наступних операцій обміну даними :

- у кожну підзадачу(i, j) передаються блоки Aij, Bij;

- для кожного рядка i грат підзадач блоки матриці A зрушуються на(i - 1) позицій вліво;

- для кожного стовпця j грат підзадач блоки матриці B зрушуються на(j - 1) позицій вгору.

Виконувані при перерозподілі матричних блоків процедури передачі даних є прикладом операції циклічного зрушення - див. розділ 3. На рис. 5.1 показаний приклад розташування блоків для грат підзадач .

Перерозподіл блоків матриць A и В

Рис. 5.1. Перерозподіл блоків початкових матриць між процесорами при виконанні алгоритму Кэннона

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



5.3 Масштабування і розподіл підзадач по процесорах

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

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

5.4 Аналіз ефективності

Перед проведенням аналізу ефективності слід зазначити, що алгоритм Кэннона відрізняється від методу Фокса тільки видом виконуваних в ході обчислень комунікаційних операцій. Як результат, використовуючи оцінки часу виконання обчислювальних операцій, приведених в п. 4.4, проведемо тільки аналіз комунікаційної складності алгоритму Кэннона.

Відповідно до правил алгоритму на етапі ініціалізації робиться перерозподіл блоків матриць і за допомогою циклічного зрушення матричних блоків по рядках і стовпцях процесорних грат. Трудомісткість виконання такої операції передачі даних істотним чином залежить від топології мережі. Для мережі із структурою повного графа усі необхідні пересилки блоків можуть бути виконані одночасно(тобто тривалість операції виявляється рівною часу передачі одного матричного блоку між сусідніми процесорами). Для мережі з топологією гіперкуба операція циклічного зрушення може потребувати виконання ітерацій. Для мережі з кільцевою структурою зв'язків необхідна кількість ітерацій виявляється рівною . Детальніше методи виконання операції циклічного зрушення розглянуті в розділі 3. Для побудови оцінки комунікаційної складності етапу ініціалізації використовується варіант топології повного графа. Час виконання початкового перерозподілу блоків може оцінюватися як:

(1.14)

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

Оцінимо тепер витрати на передачу даних між процесорами при виконанні основної частини алгоритму Кэннона. На кожній ітерації алгоритму після множення матричних блоків процесори передають свої блоки попереднім процесорам по рядках(для блоків матриці ) і стовпцях (для блоків матриці ) процесорних грат. Ці операції також можуть бути виконані процесорами паралельно і, тим самим, тривалість таких комунікаційних дій складає:

(1.15)

Оскільки кількість ітерацій алгоритму Кэннона являється рівним q, то з урахуванням оцінки (1.13) загальний час виконання паралельних обчислень може бути визначений за допомогою наступного співвідношення:

(1.16)

(у використовуваних виразах параметр визначає розмір процесорних грат).

5.5 Результати обчислювальних експериментів

Обчислювальні експерименти для оцінки ефективності паралельного алгоритму проводилися за тих же умов, що і раніше виконані експерименти(див. п. 3.5). Результати експериментів для випадків чотирьох і дев'яти процесорів приведені в таблиці 5.1.

Таблиця 5.1 Результатів обчислювальних експериментів по дослідженню паралельного алгоритму Кэннона.

Розмір об'єктів	Послідовний алгоритм	Паралельний алгоритм
		4 ЦП	9 ЦП
		Час	Присокрення	Час	Присокрення
1000	12,8787	3,0806	4,1805	1,1889	10,8324
1500	43,4731	11,1716	3,8913	4,6310	9,3872
2000	103,0561	24,0502	4,2850	14,4759	7,1191
2500	201,2915	53,1444	3,7876	23,5398	8,5511
3000	347,8434	88,2979	3,9394	36,3688	9,5643



Рис. 5.2. Залежність прискорення від розміру матриць при виконанні паралельного алгоритму Кэннона

Порівняння часу виконання експерименту і теоретичного часу, вичисленого відповідно до виразу (1.16), представлено в таблиці 5.2 і на рис. 5.3.

Таблиця 5.2. Порівняння експериментального і теоретичного часу виконання паралельного алгоритму Кэннона

Розмір матриці	Паралельний алгоритм
	4 ЦП	9 ЦП
1000	3,4485	3,0806	1,5669	1,1889
1500	11,3821	11,1716	5,1348	4,6310
2000	26,6769	24,0502	11,9912	14,4759
2500	51,7488	53,1444	23,2098	23,5398
3000	89,0138	88,2979	39,8643	36,3688



Рис. 5.3. Графік залежності експериментального і теоретичного часу виконання алгоритму Кэннона на чотирьох процесорах



6. Огляд літератури

Завдання множення матриць широко розглядається в літературі. В якості додаткового учбового матеріалу можуть бути рекомендовані роботи Воеводин В.В. і Воеводин Вл.В.(2002), Kumar(1994) і Quinn(2003). Широке обговорення питань паралельного виконання матричних обчислень виконане в роботі Dongarra, et al. (1999).

При розгляді питань програмної реалізації паралельних методів може бути рекомендована робота Blackford, et al. (1997). У цій роботі розглядається добре відома і широко використовувана в практиці паралельних обчислень програмна бібліотека чисельних методів ScaLAPACK.



Висновок

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

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

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

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



Рис. 6.1. Прискорення трьох паралельних алгоритмів при множенні матриць з використанням 4 процесорів.



Список використаних джерел

1. Гергель, В.П., Стронгин, Р.Г. (2001). Основы параллельных вычислений для многопроцессорных вычислительных систем. - Н.Новгород, ННГУ (2 изд., 2003).

2. Culler, D.E., et al. (1996). LogP: A practical model for parallel computation. - Comm. Of the ACM, 39, 11, pp. 75-85.

3. Hockney, R. (1994). The communication challenge for MPP: Intel Paragon and Meiko CS-2. - Parallel Computing, 20 (3), pp. 389-398.

4. Kumar V., Grama, A., Gupta, A., Karypis, G. (1994). Introduction to Parallel Computing. - The Benjamin/Cummings Publishing Company, Inc. (2nd edn., 2003)

5. Quinn, M. J. (2004). Parallel Programming in C with MPI and OpenMP. - New York, NY: McGraw-Hill.

6. Skillicorn, D.B., Talia, D. (1998). Models and languages for parallel computation. - ACM Computing surveys, 30, 2.

Размещено на Allbest.ru

...