Лабораторная работа 3 (Умножение двух матриц методом статического разделения на полосы (Захаров))

Посмотреть архив целиком

Национальный Исследовательский Университет

МОСКОВСКИЙ ЭНЕРГЕТИЧЕСКИЙ ИНСТИТУТ

Институт автоматики и вычислительной техники

Кафедра прикладной математики

Лабораторная работа № 3

Умножение двух матриц

методом статического разделения на полосы

Курс «Параллельные системы и параллельные вычисления»




Выполнил

студент 5 курса группы А-13-08

Захаров Антон

Преподаватель

Панков Николай Александрович


Постановка задачи

Пусть даны две прямоугольные матрицы и размерности и
соответственно:

Требуется найти матрицу (произведением) размерности :

Для нахождения произведения матриц методом статического разделения на полосы необходимо составить последовательно-параллельную программу на языке C или C++, использующую принципы нитевого распараллеливания, а также исследовать характеристики разработанной программы в зависимости от числа потоков.


Тестирование проводились на компьютере со следующей конфигурацией:

ПРОЦЕССОР Intel Core i5 2500MHz Ivy Bridge

ОПЕРАТИВНАЯ ПАМЯТЬ 16Gb DDR3 1600MHz

ФИЗИЧЕСКИЙ НАКОПИТЕЛЬ OCZ-VERTEX3 (120Gb, SATA600, SSD)

ГРАФИЧЕСКИЙ ПРОЦЕССОР AMD Radeon HD 7700 (1Gb DDR5 4.6GHz)

ОПЕРАЦИОННАЯ СИСТЕМА Windows 7 Ultimate x64 (SP1)



Последовательный алгоритм решения

Умножение матриц по определению

В соответствии с определением, произведение матриц состоит из всех возможных комбинаций скалярных произведений строк матрицы и столбцов матрицы . Элемент матрицы с индексами (i, j) есть скалярное произведение i-ой строки матрицы и j-го столбца матрицы .


  1. for (i = 0; i < m; i++) {

  2. for (j = 0; j < q; j++) {

  3. C[i][j] = 0;

  4. for (k = 0; k < n; k++)

  5. C[i][j] += A[i][k] * B[k][j];

  6. }

  7. }


На первый взгляд это минимальный объем работы, необходимый для перемножения двух матриц. Однако исследователям не удалось доказать минимальность, и в результате они обнаружили другие алгоритмы, умножающие матрицы более эффективно.


Алгоритм Штрассена

Первый алгоритм быстрого умножения матриц был разработан В. Штрассеном в 1969. В основе алгоритма лежит рекурсивное разбиение матриц на блоки. Недостатком данного метода является большая сложность программирования по сравнению со стандартным алгоритмом, численная неустойчивость и большой объём используемой памяти.

Разработано большое количество алгоритмов на основе метода Штрассена, которые улучшают его численную устойчивость и уменьшают объём используемой памяти.


Алгоритм Копперсмита-Винограда

В 1990 Копперсмит и Виноград опубликовали алгоритм, умножающий матрицы со сложностью . Этот алгоритм использует идеи, схожие с алгоритмом Штрассена. На сегодняшний день алгоритм Копперсмита-Винограда является наиболее асимптотически быстрым, но он эффективен только на очень больших матрицах и поэтому не применяется.

В 2003 Кох и др. рассмотрели в своих работах алгоритмы Штрассена и Копперсмита-Винограда в контексте теории групп. Они показали возможность существования алгоритмов умножения матриц со сложностью .


Параллельный алгоритм решения

В предлагаемой реализации метода статического разделения на полосы исходная матрица A разбивается на горизонтальные полосы. Все полосы матрицы распределяются между потоками, а обращение к элементам матрицы , расположенной в общей памяти, осуществляется по мере необходимости. При этом каждый из имеющихся потоков использует только одну полосу матрицы и всю матрицу . Перемножение полосы на матрицу (а данная операция может быть выполнена параллельно) приводит к получению частей (полос) результирующей матрицы , которые затем в совокупности и дадут искомую матрицу.

– полоса матрицы ;

– число потоков.



Результаты вычислительного эксперимента


Матрицы 1 000 × 1 000


Число
потоков

Время
решения1
(сек)

Ускорение

1

5,7562

1,0000

2

3,5368

1,6275

3

2,7934

2,0606

4

2,4386

2,3605

5

2,4630

2,3371

6

2,4670

2,3333

7

2,4676

2,3327

8

2,4740

2,3267

9

2,4746

2,3261

10

2,4860

2,3154

11

2,4990

2,3034

12

2,6638

2,1609








Время

(сек)

Зависимость времени решения задачи от числа потоков

Число

потоков





Ускорение

Зависимость ускорения параллельного решения от числа потоков

Число

потоков


Матрицы 10 000 × 10 000


Число
потоков

Время
решения2
(сек)

Ускорение

1

1447,8922

1,0000

2

784,3633

1,8459

3

561,0955

2,5805

4

428,4089

3,3797

5

407,7476

3,5510

6

379,4659

3,8156

7

385,1177

3,7596

8

390,7400

3,7055

9

408,0124

3,5486

10

412,3470

3,5113

11

446,2310

3,2447

12

491,2599

2,9473





Время

(сек)

Зависимость времени решения задачи от числа потоков

Число

потоков





Ускорение

Зависимость ускорения параллельного решения от числа потоков

Число

потоков




Матрицы 20 000 × 20 000


Число
потоков

Время
решения
(сек)

Ускорение

1

6005,4665

1,0000

2

3401,2681

1,7657

3

2375,1655

2,5284

4

1820,5416

3,2987

5

1873,0987

3,2062

6

1891,9394

3,1742

7

1873,6492

3,2052

8

1831,9423

3,2782

9

1951,0563

3,0781

10

2089,2749

2,8744

11

2116,7367

2,8371

12

2085,5324

2,8796





Время

(сек)

Зависимость времени решения задачи от числа потоков

Число

потоков





Ускорение

Зависимость ускорения параллельного решения от числа потоков

Число

потоков




Приложение. Код программы.


// Отключение предупреждений безопасности

#define _CRT_SECURE_NO_WARNINGS

#define _CRT_SECURE_NO_DEPRECATE

#define _CRT_NONSTDC_NO_DEPRECATE


// Коды ошибок

#define E_THREAD_ALLOC 1

#define E_THREAD_CREATE 2


#include

#include

#include

#include

#include

#include

#include

#include

#include

#include



// --- Настройки программы ---


// Включение режима отладки

// #define DEBUG


// Число запусков теста

#define TEST_COUNT 5



// --- Настройки потоков ---


// Ограничения на число потоков

#define MIN_THREADS 1

#define MAX_THREADS 12



int gen_random(float, float);

int matrix_create(const char*, int, int);

void print_log(const char*, ...);

void print_err(const char*);


DWORD WINAPI MyThreadFunction(LPVOID lpParam);


// Структура данных потока

typedef struct ThreadData {

int i; // Номер потока

} TDATA, *PTDATA;


int **A;

int **B;

int **C;


// Размерности матриц А{mxn}, B{nxq}, C{mxq}

int n = 20000, m = 20000, q = 20000;


// Пути к файлам исходных данных и результата

const char *A_path = "A.txt"; FILE *Astream;

const char *B_path = "B.txt"; FILE *Bstream;

const char *C_path = "C.txt"; FILE *Cstream;

const char *L_path = "lab.log"; FILE *Lstream;


// Числа строк матрицы A, обрабатываемые потоками

int Arows[MAX_THREADS];


// Номера строк матрицы A, обрабатываемые потоками

int Nrows[MAX_THREADS];


int _tmain()

{

srand((UINT32) time(NULL));

int i, j, h, d;


// Время выполнения каждого из запусков программы

double time[TEST_COUNT];


PTDATA pDataArray[MAX_THREADS];

DWORD dwThreadIdArray[MAX_THREADS];

HANDLE hThreadArray[MAX_THREADS];


// Запускаем программу для 1, 2, ..., MAX_THREADS потоков

for(int threads = MIN_THREADS; threads <= MAX_THREADS; threads++)

{

// Запускаем тест TEST_COUNT раз

for(h = 0; h < TEST_COUNT; h++)

{

// Фиксируем время начала работы программы

time[h] = (double) clock();


// Пытаемся открыть файлы с матрицами A и B

Astream = fopen(A_path, "r");

Bstream = fopen(B_path, "r");

// Обработка ошибок при открытии файлов

if(Astream == NULL || Bstream == NULL)

{

if((Astream != NULL && fclose(Astream)) || (Bstream != NULL && fclose(Bstream)))

{

print_log("[ERROR] Не удалось закрыть один из исходных файлов");

return 0;

}


// Генереция новых файлов с матрицами A и B

if(!matrix_create(A_path, m, n) || !matrix_create(B_path, n, q))

{

print_log("[ERROR] Не удалось сгенерировать матрицы A и B");

return 0;

}

// Не учитываем время, затраченное на генерацию матриц

time[h] = (double) clock();


Astream = fopen(A_path, "r");

Bstream = fopen(B_path, "r");


if(Astream == NULL || Bstream == NULL)

{

print_log("[ERROR] Не удалось сгенерировать матрицы A и B");

return 0;

}

}


// Считываем размерности матриц

if(fscanf(Astream, "%d %d\n", &m, &n) <= 0 || fscanf(Bstream, "%d %d\n", &d, &q) <= 0)

{

print_log("[ERROR] Не удалось считать размерности матриц");

fclose(Astream);

fclose(Bstream);

return 0;

}


if(n != d || n <= 0 || m <= 0 || q <= 0)

{

if(n != d)

print_log("[ERROR] Не верно заданы размерности матриц:

Число столбцов в матрице A не совпадает с числом строк в матрице B");


if(n <= 0 || m <= 0 || q <= 0)

print_log("[ERROR] Не верно заданы размерности матриц:

Размерность не может быть меньше либо равна нулю");


fclose(Astream);

fclose(Bstream);

return 0;

}


// Делим строки матрицы A между потоками

d = (int) floor(double(m / threads));

for(i = 0; i < threads; i++)

{

Arows[i] = d;

Nrows[i] = i * d;

}

d = m % threads;

for(i = 0; d != 0; d--)

{

Arows[i]++;

for(j = i + 1; j < threads; j++)

Nrows[j]++;

i = (i < threads) ? i + 1 : 0;

}


// Запись информации о группе тестов в лог-файл

if(!h)

{

time[h] = (double) clock() - time[h];

if(threads == MIN_THREADS)

print_log("----------------------------------------------\n

Размерности матриц m=%d n=%d q=%d\n\n", m, n, q);

print_log("Число потоков %d\n\n", threads);

print_log("[DEBUG] Распределение строк матрицы A между потоками:\n{");

for(i = 0; i < threads - 1; i++)

print_log("%d (%d), ", Arows[i], Nrows[i]);

print_log("%d (%d)}\n\n", Arows[threads - 1], Nrows[threads - 1]);

time[h] = (double) clock() - time[h];

}


// Выделение памяти

A = (int **) malloc(m * sizeof(int *));

B = (int **) malloc(n * sizeof(int *));

C = (int **) malloc(m * sizeof(int *));


if(A == NULL || B == NULL || C == NULL)

{

print_err("Не удалось выделить память");

fclose(Astream);

fclose(Bstream);

if(A != NULL) free(A);

if(B != NULL) free(B);

if(C != NULL) free(C);

return 0;

}


for(i = 0; i < m; i++)

{

A[i] = (int *) malloc(n * sizeof(int));

C[i] = (int *) malloc(q * sizeof(int));


if(A[i] == NULL || C[i] == NULL)

{

print_err("Не удалось выделить память");

fclose(Astream);

fclose(Bstream);

free(A); free(B); free(C);

return 0;

}

}


for(i = 0; i < n; i++)

{

B[i] = (int *) malloc(q * sizeof(int));

if(B[i] == NULL)

{

print_err("Не удалось выделить память");

fclose(Astream);

fclose(Bstream);

free(A); free(B); free(C);

return 0;

}

}


// Чтение матрицы A

for(i = 0; i < m; i++)

{

for(j = 0; j < n; j++)

{

if(fscanf(Astream, "%d", &A[i][j]) <= 0)

{

print_err("Не удалось считать элемент матрицы A");

fclose(Astream);

fclose(Bstream);

free(A); free(B); free(C);

return 0;

}

}

}

fclose(Astream);

// Чтение матрицы B

for(i = 0; i < n; i++)

{

for(j = 0; j < q; j++)

{

if(fscanf(Bstream, "%d", &B[i][j]) <= 0)

{

print_err("Не удалось считать элемент матрицы B");

fclose(Bstream);

free(A); free(B); free(C);

return 0;

}

}

}

fclose(Bstream);


// Обнуление результирующей матрицы C

for(i = 0; i < m; i++)

for(j = 0; j < q; j++)

C[i][j] = 0;


for(i = 0; i < threads; i++)

{

// Выделяем память для потока

pDataArray[i] = (PTDATA) HeapAlloc(GetProcessHeap(), HEAP_ZERO_MEMORY, sizeof(TDATA));


// Если не удалось выделить память, то прерываем выполнение программы

if(pDataArray[i] == NULL)

{

print_log("[ERROR] Не удалось выделить память для потока\n");

free(A); free(B); free(C);

ExitProcess(E_THREAD_ALLOC);

}

pDataArray[i]->i = i;


// Создаём поток

hThreadArray[i] = CreateThread(

NULL, // атрибуты безопасности по умолчанию

0, // используем размер стека по умолчанию

MyThreadFunction, // имя функции потока

pDataArray[i], // аргумент функции потока

0, // используем стандартные флаги

&dwThreadIdArray[i] // получаем идентификатор потока

);


// Если не удалось запустить поток, то прерываем выполнение программы

if(hThreadArray[i] == NULL)

{

ExitProcess(E_THREAD_CREATE);

}

}


// Ждём, пока не выполнятся все потоки

WaitForMultipleObjects(threads, hThreadArray, TRUE, INFINITE);


// Закрываем все дескрипторы потоков и освобождаем выделенную память

for(i = 0; i < threads; i++)

{

CloseHandle(hThreadArray[i]);

if(pDataArray[i] != NULL)

{

HeapFree(GetProcessHeap(), 0, pDataArray[i]);

pDataArray[i] = NULL;

}

}

free(A);

free(B);


// Записываем результат в файл

Cstream = fopen(C_path, "w");


if (Cstream == NULL)

{

print_err("Не удалось создать файл для произведения матриц");

return 0;

}


for(i = 0; i < m; i++)

{

for(j = 0; j < q - 1; j++)

{

if(fprintf(Cstream, "%d ", C[i][j]) < 0)

{

fclose(Cstream);

return 0;

}


}

if(fprintf(Cstream, "%d\n", C[i][q - 1]) < 0)

{

fclose(Cstream);

return 0;

}

}


free(C);

if(fclose(Cstream))

{

print_err("Не удалось закрыть поток");

return 0;

}


// Записываем время в лог-файл

time[h] = ((double) clock() - time[h]) / CLOCKS_PER_SEC;

print_log("Время (%d): %f сек\n", h + 1, time[h]);

}


// Записываем среднее время в лог-файл

time[0] = time[0] / TEST_COUNT;

for(i = 1; i < TEST_COUNT; i++)

time[0] += time[i] / TEST_COUNT;

print_log("Среднее : %f сек\n\n", time[0]);

}

return 0;

}


DWORD WINAPI MyThreadFunction(LPVOID lpData)

{

PTDATA data = (PTDATA) lpData;

int t = data->i;


int first = Nrows[t];






Чтобы не видеть здесь видео-рекламу достаточно стать зарегистрированным пользователем.
Чтобы не видеть никакую рекламу на сайте, нужно стать VIP-пользователем.
Это можно сделать совершенно бесплатно. Читайте подробности тут.