МІНІСТЕРСТВО ОСВІТИ І НАУКИ, МОЛОДІ ТА СПОРТУ УКРАЇНИ НУ «ЛЬВІВСЬКА ПОЛІТЕХНІКА» Кафедра ЕОМ ЛАБОРАТОРНА РОБОТА №6 з предмету “ Паралельні та розподілені обчислення ” на тему: “ ПАРАЛЕЛЬНІ АЛГОРИТМИ МНОЖЕННЯ МАТРИЦІ НА ВЕКТОР ” Варіант №9 Мета. Ознайомитись з методами організації паралельного множення матриці на вектор та розробити паралельну програму з використанням технології MPI. ТЕОРЕТИЧНІ ВІДОМОСТІ Матриці та операції над ними широко використовуються при математичному моделюванні найрізноманітніших процесів, явищ і систем. Матричні обчислення складають основу багатьох наукових і інженерних розрахунків. Враховуючи важливість ефективного виконання матричних обчислень багато стандартних бібліотек програм містять процедури для різних матричних операцій. Об'єм програмного забезпечення для операцій над матрицями постійно збільшується - розробляються нові оптимальні структури даних для зберігання матриць спеціального типу (трикутних, стрічкових, розріджених тощо), створюються різні високоефективні машинно-залежні реалізації алгоритмів, проводяться теоретичні дослідження для пошуку швидших методів матричних обчислень. Будучи обчислювально трудомісткими, матричні обчислення є класичною областю застосування паралельних обчислень. З одного боку, використання високопродуктивних багатопроцесорних систем дозволяє істотно підвищити складність завдань, які розв’язуються. З іншого боку, через своє достатньо просте формулювання матричні операції надають прекрасну можливість для демонстрації багатьох прийомів і методів паралельного програмування. У даній лабораторній роботі аналізуються методи паралельних обчислень для операції векторно-матричного множення. 1. Принципи розпаралелювання Для багатьох методів матричних обчислень характерним є повторення одних і тих же операцій для різних даних. Дана властивість свідчить про наявність паралелізму за даними при виконанні матричних обчислень, і, як результат, розпаралелювання матричних операцій зводиться, в більшості випадків, до розбиття оброблюваних матриць між процесорами використовуваної обчислювальної системи. Вибір способу поділу матриць приводить до визначення конкретного методу паралельних обчислень; існування різних схем розподілу даних породжує ряд паралельних алгоритмів матричних обчислень. Найбільш загальні і широко використовувані способи поділу матриць полягають в розбитті даних на смуги (по вертикалі або горизонталі) або на прямокутні фрагменти (блоки). 1.1. Стрічкове розбиття матриці. При стрічковому (block-striped) розбитті кожному процесору виділяється певна підмножина рядків (rowwise або горизонтальне розбиття) або стовпців (columnwise або вертикальне розбиття) матриці. При такому підході для горизонтального розбиття по рядках, наприклад, матриця A представляється у вигляді (1)
, де (1) ai = (ai1,ai2,...,ain), 0 <= i < m i-й рядок матриці A (передбачається, що кількість рядків m кратна кількості процесорів p, тобто m = k·p). Інший можливий підхід до формування смуг полягає в застосуванні тієї чи іншої схеми чергування (циклічності) рядків або стовпців. Як правило, для чергування використовується кількість процесорів p - в цьому випадку при горизонтальному розбитті матриця A приймає вигляд
(2) 1.2. Блокове розбиття матриці. При блоковому (chessboard block) розбитті матриця ділиться на прямокутні набори. Хай кількість процесорів складає p = s·q, кількість рядків матриці є кратним s, а кількість стовпців - кратним q, тобто m = k·s і n = l·q. Представимо початкову матрицю A у вигляді набору прямокутних блоків таким чином (3):
(3) де Aij - блок матриці, що складається з елементів:
(4) При такому підході доцільно, щоб обчислювальна система мала фізичну або, принаймні, логічну топологію процесорних граток з s рядків і q стовпців. При такому розбитті даних, сусідні в структурі граток процесори, обробляють суміжні блоки початкової матриці. Треба зазначити, що і для блокової схеми може бути застосоване циклічне чергування рядків і стовпців. У лабораторній роботі розглядаються три паралельні алгоритми для множення квадратної матриці на вектор. Кожен підхід заснований на різному типі розбиття початкових даних (елементів матриці і вектора) між процесорами. Розбиття даних міняє схему взаємодії процесорів, тому кожен з представлених методів істотним чином відрізняється від решти.
а) б) в) Рис. 1. Способи розбиття елементів матриці між процесорами обчислювальної системи На рис. 1 схематично неведені способи розбиття матриць між процесорами: а) горизонтальне розбиття, б) вертикальне розбиття та в) блокове розбиття матриці. 2. Постановка задачі В результаті перемноження матриці А розмірності m х n і вектора b, що складається з n елементів, отримується вектор розміру m, кожен i-й елемент якого є результат скалярного множення i-того рядка матриці А (позначимо цей рядок aі) і вектора b:

(4) Тим самим отримання результуючого вектора С припускає повторення m однотипних операцій по множенню рядків матриці A і вектора b. Кожна така операція включає множення елементів рядка матриці і вектора b (n операцій) і подальше підсумовування отриманих множень (n - 1 операцій). Загальна кількість необхідних скалярних операцій є величина T1 = m·(2n-1). (5) 2.1. Послідовний алгоритм. Послідовний алгоритм перемноження матриці на вектор може бути представлений таким чином. // Послідовний алгоритм множення матриці на вектор for (i = 0; i < m; i++){ с[i]= 0; for (j = 0; j < n; j++){ с[i]+= A[i][j]*b[j] } } Векторно-матричне множення - це послідовність обчислення скалярних добутків. Оскільки кожне обчислення скалярного добутку векторів довжини n вимагає виконання n операцій множення і (n-1) операцій додавання, його трудомісткість становить O(n). Для виконання векторно-матричного множення необхідно здійснити m операцій обчислення скалярного добутку, таким чином, алгоритм має трудомісткість порядку O(mn). 2.2. Розділення даних. При виконанні паралельних алгоритмів перемноження матриці на вектор, окрім матриці А, необхідно розбити вектор b і вектор результату с. Елементи векторів можна продублювати, тобто скопіювати всі елементи вектора на всі процесори, складові багатопроцесорної обчислювальної системи, або розділити між процесорами. При блоковому розбитті вектора з n елементів кожен процесор обробляє безперервну послідовність із k елементів вектора (припускається, що розмірність вектора n без остачі ділиться на число процесорів, тобто n= k·p). Пояснимо, чому дублювання векторів b і с між процесорами є допустимим рішенням (далі для простоти викладу вважатимемо, що m=n). Вектори b і с складаються з n елементів, тобто містять стільки ж даних, скільки і один рядок або один стовпець матриці. Якщо процесор зберігає рядок або стовпець матриці і одиночні елементи векторів b і с, то загальне число елементів, що зберігаються, має трудомісткість порядку O(n). Якщо процесор зберігає рядок (стовпець) матриці і всі елементи векторів b і с, то загальна кількість елементів, що зберігаються, також порядку O(n). Таким чином, при дублюванні і при розбитті векторів вимоги до об'єму пам'яті з одного класу складності. ЗАВДАННЯ № варіанту Розмір матриці Тип розбиття Кількість процесорів  9 250 750 блокове 25   ВИКОНАННЯ ЗАВДАННЯ Виходячи з вище вказаних тверджень можливі розміри блоків 24 blocks: 125x60 1 block: 250x30 Розібємо дані блоки між процесорами
З даної схеми видно, що матриця розбита на 25 блоків однакового розміру, тільки остання підматриця має інші розміри, але кількість елементів, які входять в усі блоки однакова.
Отже на кожний процесор буде переслано (250*750)/25=7500 елементів. Кількість операцій на кожному процесорі буде рівна: (750/25*250)(2*250-1)= 3742500 ТЕКСТ ПРОГРАМИ #include<stdio.h> #include<stdlib.h> #include<mpi.h> #include<time.h> #include <time.h> void ProcessInitialization (double* &pMatrix, double* &pVector, double* &pResult, double* &pProcRows, double* &pProcResult, int &RowNum); void DataDistribution(double* pMatrix, double* pProcRows, double* pVector, int RowNum); void ParallelResultCalculation(double* pProcRows, double* pVector, double* pProcResult, int RowNum); void ResultReplication(double* pProcResult, double* pResult, int RowNum); void RandomDataInitialization(double * &pMatrix, double * &pVector); void ProcessTermination(double * pMatrix, double * pVector, double * pResult, double * pProcRows, double * pProcResult); void SeperateMatrix(double* pNewMatrix, double* pMatrix, int * pSendInd, int ColNum, int RowNum); void AddSubMatrix(double *pProcResult,int ColNum, int RowNum); int ProcNum, ProcRank; int Row = 250,Col = 750; // Розміри початкової матриці і вектора double* pProcRows; int ARow, ACol, RowNum1, ColNum1; int bFlag = 0; // Множення матриці на вектор - стрічкове горизонтальне розбиття // (початковий і результуючий вектори дублюються між процесами) void main(int argc, char* argv[]) { double* pMatrix; // Перший аргумент - початкова матриця double* pVector; // Другий аргумент - початковий вектор double* pResult; // Результат множення матриці на вектор double* pProcRows; double* pProcResult; int RowNum; double Start, Finish, Duration; MPI_Init(&argc, &argv); MPI_Comm_size(MPI_COMM_WORLD, &ProcNum); MPI_Comm_rank(MPI_COMM_WORLD, &ProcRank); // Виділення пам'яті і ініціалізація початкових даних ProcessInitialization(pMatrix, pVector, pResult, pProcRows, pProcResult, RowNum); // Розподіл початкових даних між процесами DataDistribution(pMatrix, pProcRows, pVector, RowNum); // Паралельне виконання множення матриці на вектор ParallelResultCalculation(pProcRows, pVector, pProcResult, RowNum); // Збір результуючого вектора на всіх процесах ResultReplication(pProcResult, pResult, RowNum); MPI_Finalize(); // Завершення процесу обчислень ProcessTermination(pMatrix, pVector, pResult, pProcRows, pProcResult); } // Функція для виділення пам'яті і ініціалізації початкових даних void ProcessInitialization (double* &pMatrix, double* &pVector, double* &pResult, double* &pProcRows, double* &pProcResult, int &RowNum) { int Matrix; int SubMatrix; int i; Matrix = Row*Col; SubMatrix = Matrix/ProcNum; if (SubMatrix%Row==0) { ColNum1 = SubMatrix/(Row/2); RowNum1 = Row/2; } else if (SubMatrix%Col==0) { RowNum1 = SubMatrix/(Col/2); ColNum1 = Col/2; } if(Row%RowNum1==0) { if(Col%ColNum1==0) { printf("%d blocks: %dx%d\n",ProcNum,RowNum1,ColNum1); } else { printf("%d blocks: %dx%d\n",ProcNum-1,RowNum1,ColNum1); printf("1 block: %dx%d\n",Row,ColNum1/2); ARow = Row; ACol = ColNum1/2; bFlag = 1; } } else { if(Col%ColNum1==0) { printf("%d blocks: %dx%d\n",ProcNum-1,RowNum1,ColNum1); printf("1 block: %dx%d\n",RowNum/2,Col); ARow = RowNum/2; ACol = Col; bFlag = 1; } } int RestRows; // Кількість рядків матриці, які ще // не розподілені MPI_Bcast(&Row, 1, MPI_INT, 0, MPI_COMM_WORLD); RestRows = Row; for (int i=0; i<ProcRank; i++) RestRows = RestRows-RestRows/(ProcNum-i); RowNum = RestRows/(ProcNum-ProcRank); pVector = new double [Col]; pResult = new double [Row]; pProcRows = new double [RowNum*Col]; pProcResult = new double [RowNum]; if (ProcRank == 0) { pMatrix = new double [Row*Col]; RandomDataInitialization(pMatrix, pVector); } } // Функція для розбиття початкових даних між процесами void DataDistribution(double* pMatrix, double* pProcRows, double* pVector, int RowNum) { int *pSendNum; // Кількість елементів, що посилаються процесу int *pSendInd; // Індекс першого елементу даних // посиланого процесу int RestRows=Row; // Кількість рядків матриці, які ще // не розподілені MPI_Bcast(pVector, Col, MPI_DOUBLE, 0, MPI_COMM_WORLD); // Виділення пам'яті для зберігання тимчасових об'єктів pSendInd = new int [ProcNum]; pSendNum = new int [ProcNum]; // Визначення положення рядків матриці, призначених // кожному процесу RowNum = (Row/ProcNum); pSendNum[0] = RowNum*Col; pSendInd[0] = 0; for (int i=1; i<ProcNum; i++) { RestRows -= RowNum; RowNum = RestRows/(ProcNum-i); pSendNum[i] = RowNum*Col; pSendInd[i] = pSendInd[i-1]+pSendNum[i-1]; } // Розсилка рядків матриці MPI_Scatterv(pMatrix, pSendNum, pSendInd, MPI_DOUBLE, pProcRows, pSendNum[ProcRank], MPI_DOUBLE, 0, MPI_COMM_WORLD); // Звільнення пам'яті delete [] pSendNum; delete [] pSendInd; } // Функція для обчислення частини результуючого вектора void ParallelResultCalculation(double* pProcRows, double* pVector, double* pProcResult, int RowNum) { int i, j; for (i=0; i<RowNum; i++){ pProcResult[i] = 0; for (j=0; j<Col; j++) pProcResult[i] += pProcRows[i*Col+j]*pVector[j]; } AddSubMatrix(pProcResult, RowNum, ColNum1); } // Функція для збору результуючого вектора на всіх процесах void ResultReplication(double* pProcResult, double* pResult, int RowNum) { int *pReceiveNum; // Кількість елементів, що посилаються процесом int *pReceiveInd; // Індекс елементу даних в результуючому // векторі int RestRows=Row; // Кількість рядків матриці, які ще не // розподілені int i; // Виділення пам'яті для тимчасових об'єктів pReceiveNum = new int [ProcNum]; pReceiveInd = new int [ProcNum]; // Визначення положення блоків результуючого вектора pReceiveInd[0] = 0; pReceiveNum[0] = Row/ProcNum; for (i=1; i<ProcNum; i++) { RestRows -= pReceiveNum[i-1]; pReceiveNum[i] = RestRows/(ProcNum-i); pReceiveInd[i] = pReceiveInd[i-1]+pReceiveNum[i-1]; } // Збір всього результуючого вектора на всіх процесах MPI_Allgatherv(pProcResult, pReceiveNum[ProcRank], MPI_DOUBLE, pResult, pReceiveNum, pReceiveInd, MPI_DOUBLE, MPI_COMM_WORLD); // Звільнення пам'яті delete [] pReceiveNum; delete [] pReceiveInd; } void RandomDataInitialization(double * &pMatrix, double * &pVector) { for(int i=0;i<Row*Col;i++) { if(i<Col) pVector[i] = (double)(rand()%10+1); pMatrix[i] = (double)(rand()%10+1); } } void ProcessTermination(double * pMatrix, double * pVector, double * pResult, double * pProcRows, double * pProcResult) { for(int i=0;i<Row;i++) printf("C[%d] = %6.2f\n",i,pResult[i]); delete [] pMatrix; delete [] pVector; delete [] pResult; delete [] pProcRows; delete [] pProcResult; } void SeperateBlocks(double* &pNewMatrix, double* pMatrix, int * &pSendInd, int ColNum, int RowNum) { int l=0; if(bFlag && ARow==Row) { for(int m=0;m<2;m++) for(int k=0;k<ProcNum/2;k++){ for(int i=m*RowNum;i<(m+1)*RowNum;i++) for(int j=k*ColNum;j<(k+1)*ColNum;j++,l++) pNewMatrix[l] = pMatrix[i*Col+j]; pSendInd[m+k] = (m+k)*ColNum*RowNum; } for(int i=0;i<Row;i++) for(int j=Col-ColNum/2;j<Col;j++,l++) pNewMatrix[l]=pMatrix[i*Col+j]; pSendInd[ProcNum-1] = (ProcNum-1)*ColNum*RowNum; } if(bFlag && ACol==Col) { for(int m=0;m<ProcNum/2;m++) for(int k=0;k<2;k++){ for(int i=m*RowNum;i<(m+1)*RowNum;i++) for(int j=k*ColNum;j<(k+1)*ColNum;j++,l++) pNewMatrix[l] = pMatrix[i*Col+j]; pSendInd[m+k] = (m+k)*ColNum*RowNum; } for(int i=Row-RowNum/2;i<Row;i++) for(int j=0;j<Col;j++,l++) pNewMatrix[l]=pMatrix[i*Col+j]; pSendInd[ProcNum-1] = (ProcNum-1)*ColNum*RowNum; } if(!bFlag) { for(int m=0;m<2;m++) for(int k=0;k<ProcNum/2;k++) for(int i=m*RowNum;i<(m+1)*RowNum;i++){ for(int j=k*ColNum;j<(k+1)*ColNum;j++,l++) pNewMatrix[l] = pMatrix[i*Col+j]; pSendInd[m+k] = (m+k)*ColNum*RowNum; } } } void AddSubMatrix(double *pProcResult, int ColNum, int RowNum) { double *pNewResult = new double[Col]; for(int m=Col;m<2;m++) for(int k=0;k<ProcNum/2;k++) for(int i=m*RowNum;i<(m+1)*RowNum;i++) for(int j=k*ColNum;j<(k+1)*ColNum;j++) pNewResult[j] += pProcResult[i*Col+j]; } РЕЗУЛЬТАТИ ВИКОНАННЯ ПРОГРАМИ
ВИСНОВКИ