Spis treści
Projekt dotyczy zrealizowania typowych operacji przeprowadzanych na macierzach. Przygotowalismy programy dzialajace na klastrze MPI i pod procesorem CELL/B.E. Naszym celem jest sprawdzenie jak duzy wplyw na czas wykonywanych operacji ma liczba uzytych procesorow.
a) Algorytm sekwencyjny - wyznacznik macierzy o zadanym wymiarze
Implementacja dla zadnanej macierzy 3x3
Wyznacznik dowolnego stopnia mozna obliczyc wprost z definicji, jednak jest to bardzo czasochlonne. Lepiej skorzystac z faktu, iz dodanie lub odjecie od dowolnego wiersza/kolumny innego wiersza/kolumny lub kombinacji liniowej tychze nie zmienia wartosci wyznacznika i tym sposobem sprowadzic jak najwiecej elementow pewnego wiersza/kolumny do zera. Zastosowanie twierdzenia Laplace'a powinnismy kontynuowac do uzyskania macierzy, ktorej wyznacznik mozna obliczyc wprost (drugiego, trzeciego stopnia)
Kod programu:
#include <iostream>
using namespace std;
int main() {
double det;
int n=3;
double a[3][3] = {{1,5,5},{2,2,3},{2,5,2}};
for(int i=0;i<n-1;i++)
for(int j=i+1;j<n;j++)
for(int k=i+1;k<n;k++) a[j][k]-= a[j][i] / a[i][i] * a[i][k]; det=1;
for(int i=0;i<n;i++) det *= a[i][i];
printf("%01.5f",det); cin.get();
return 0;
}b) Algorytm rownolegly
Mnozenie macierzy o zadanym rozmiarze
Mnozenie macierzy A i B, w wyniku, ktorego otrzymujemy iloczyn AB, jest wykonalne tylko wtedy, gdy liczba kolumn macierzy A jest rowna liczbie wierszy macierzy B. Jezeli macierz A z elementemi postaci aij ma wymiar mxp i macierz B z elementami postaci bij ma wymiar pxn, to ich iloczyn AB = C jest macierza wymiaru mxn o elementach cij bedacych suma iloczynow:
Mnozenie macierzy nie jest przemienne, tzn. na ogol AB jest rozne od BA.
Kod programu:
#include <stdio.h>;
#include <stdlib.h>;
#include <time.h>;
#include "mpi.h"
#define MASTER 0
double *initMatrix(int size) {
double *matrix = (double *)
malloc(sizeof(double) * size * size);
return matrix;
}
double *genMatrix(int n) {
int maxRand = 10;
int i;
double *matrix = initMatrix(n);
for (i = 0; i < n * n; i++) {
matrix[i] = rand() % maxRand;
}
return matrix;
}
int main(int argc, char** argv) {
int myRank, nrProcs;
int i, j, k;
int SIZE_OF_MATRIX;
double *matrixA,*matrixB, *matrixResult;
MPI_Datatype Column_VECTOR_TYPE, Column_TYPE;
double startTime, stopTime;
int lengths[2];
MPI_Aint disps[2];
MPI_Datatype types[2];
MPI_Init(&argc, &argv);
MPI_Comm_rank(MPI_COMM_WORLD, &myRank);
MPI_Comm_size(MPI_COMM_WORLD, &nrProcs);
if (argc != 2) {
printf("Usage: %s <Size of matrix>; \n", argv[0]);
return -1;
}
SIZE_OF_MATRIX = atoi(argv[1]);
srand(time(0) + myRank);
if (myRank == MASTER) {
matrixA = genMatrix(SIZE_OF_MATRIX);
matrixB = genMatrix(SIZE_OF_MATRIX);
matrixResult = initMatrix(SIZE_OF_MATRIX);
}
else {
matrixA = initMatrix(SIZE_OF_MATRIX);
matrixB = initMatrix(SIZE_OF_MATRIX);
matrixResult = initMatrix(SIZE_OF_MATRIX);
}
if (myRank == MASTER) {
startTime = MPI_Wtime();
}
MPI_Bcast(matrixA, SIZE_OF_MATRIX * SIZE_OF_MATRIX, MPI_DOUBLE, MASTER, MPI_COMM_WORLD);
MPI_Type_vector(SIZE_OF_MATRIX, 1, SIZE_OF_MATRIX, MPI_DOUBLE, &Column_VECTOR_TYPE);
MPI_Type_commit(&Column_VECTOR_TYPE);
types[0] = Column_VECTOR_TYPE;
types[1] = MPI_UB;
lengths[0] = 1;
lengths[1] = 1;
disps[0] = 0;
disps[1] = 2 * sizeof(void*);
MPI_Type_struct(2, lengths, disps, types, &Column_TYPE);
MPI_Type_commit(&Column_TYPE);
MPI_Scatter(matrixB, SIZE_OF_MATRIX / nrProcs, Column_TYPE, matrixB, SIZE_OF_MATRIX / nrProcs, Column_TYPE, MASTER, MPI_COMM_WORLD);
for (j = 0; j < SIZE_OF_MATRIX / nrProcs;j++) {
for (i = 0; i < SIZE_OF_MATRIX; i++) {
matrixResult[i * SIZE_OF_MATRIX + j] = 0.0;
for (k = 0; k < SIZE_OF_MATRIX; k++)
matrixResult[i * SIZE_OF_MATRIX + j] += matrixA[i * SIZE_OF_MATRIX+k] * matrixB[k * SIZE_OF_MATRIX + j];
}
}
MPI_Gather(matrixResult, SIZE_OF_MATRIX / nrProcs, Column_TYPE, matrixResult,
SIZE_OF_MATRIX/nrProcs, Column_TYPE, MASTER, MPI_COMM_WORLD);
if (myRank == MASTER) {
stopTime = MPI_Wtime();
printf("For %d processes, time:%.3f\n", nrProcs, stopTime - startTime);
}
MPI_Finalize();
return 0;
}c) Podsumowanie
Podczas operacji wykonywanych na duzych macierzach wykorzystanie najwiekszej liczby procesorow ma znaczacy wplyw na czas wykonania mnozenia macierzy im wieksza liczba procesorow tym czas byl krotszy.
Przy operacjach na malych macierzach uzycie wiekszej ilosci procesorow nie przekladalo sie znaczaco na czas wykonywania
Dla procesora CELL stworzony został program dodający do siebie 2 macierze o zadeklarowanych przez użytkownika rozmiarach.
Algorytm rownolegly
Program dodawania macierzy jest programem działającym równolegle. Daje możliwośc wyboru ilości wątków spe, które mają zostac utworzone. Macierze są dodawane wektorowo. Wektor składa się z 4 elementów, gdyż elementy macierzy są liczbami typu float, a typ wektora jest 16 bajtowy.
plik dma_example_multi_spe.cpp:
#include <stdlib.h>
#include <stdio.h>
#include <errno.h>
#include <string.h>
#include <libspe2.h>
#include <pthread.h>
#include <single_buffer.h>
#include <dma_example.h>
#include <sys/time.h>
typedef struct _ppu_thread_data
{
spe_context_ptr_t spe_context; //spe context
void* control_block_ptr; //pointer to the control block
} ppu_thread_data_t;
/* this is the pointer to the SPE code */
extern spe_program_handle_t dma_example_spu;
/* PPE pthread function that starts the SPE thread */
void *ppu_pthread_function(void *argp) {
ppu_thread_data_t *thread_data = (ppu_thread_data_t *)argp;
spe_context_ptr_t ctx;
unsigned int entry = SPE_DEFAULT_ENTRY;
ctx = thread_data->spe_context;
/* start the SPE thread. We want to pass the pointer to the control block to the SPE */
if (spe_context_run(ctx, &entry, 0, thread_data->control_block_ptr, NULL, NULL) < 0) {
perror ("Failed running context");
exit (1);
}
pthread_exit(NULL);
}
//powyzej jest standardowy program
int X = 1024*2 ,Y = 1024*2;
float *macierz1, *macierz2, *wynik;
void Pokaz_macierz(float **tabl);
void Zrob_macierz(float * mac, float **macierz);
void Przeksztalc_macierz(float **macierz, float *temp_macierz);
void Do_Everything(unsigned int number) ;
// *** MAIN *** //
int main()
{
printf("Podaj rozmiar macierzy : ");scanf("%d",&X);
Y=X;
int wybor;
//zmienne do obsługi pomiaru czasu
double tstart,tend;
timeval tim;
//rezerwacja pamieci dla operacji na macierzach, trzy ciagle obszary pamieci, na ktorych beda przeprowadzane‚ operacje
//to nie sa macierze, bo macierze sa tablicami, macierz jest kopiowana na taki obszar,
//potem sa operacje wykonywane i na koniec jest zamiana z powrotem na macierz
int rc = posix_memalign ((void**)(&macierz1), 128, X*Y * sizeof(float));
if (rc != 0)
{
fprintf (stderr, "Failed allocating space for control block\n");
exit (1);
}
rc = posix_memalign ((void**)(&macierz2), 128, X*Y * sizeof(float));
if (rc != 0)
{
fprintf (stderr, "Failed allocating space for control block\n");
exit (1);
}
rc = posix_memalign ((void**)(&wynik), 128, X*Y * sizeof(float));
if (rc != 0)
{
fprintf (stderr, "Failed allocating space for control block\n");
exit (1);
}
memset (wynik, 0, X*Y * sizeof (float));
//deklaracje macierzy tych w formie tablic, dla zachowania logiki programu,
//czyli dla zachowania motta operacje na macierzach(czyli musimy miec jakies macierze, a nie tylko ciagly obszar danych w pamieci)
float **mac1 = new float * [X];
float **mac2 = new float * [X];
float **mac_wynik = new float * [X];
for (int i=0; i<X; i++)
{
mac1[i] = new float[Y];
mac2[i] = new float[Y];
mac_wynik[i] = new float[Y];
}
for(int i = 0 ; i<X; i++)
for (int j=0; j<Y; j++)
{
mac1[i][j] = i+j;
mac2[i][j] = i+j;
mac_wynik[i][j] = 0;
}
//zamiana macierzy na ciagly obszar w pamieci , na ktorym bedziemy wykonywac wszystkie operacje
Przeksztalc_macierz(mac1,macierz1);
Przeksztalc_macierz(mac2,macierz2);
do {
printf("\nWitamy w programie Macierze !\n\n");
printf("1 - Pokaz macierz1\n");
printf("2 - Pokaz macierz2\n");
printf("3 - Pokaz wynik\n");
printf("4 - Dodaj \n");
printf("9 - Zmiana ilosci watkow SPE \n");
printf("0 - Zakonczenie programu\n");
scanf("%d",&wybor);
switch (wybor)
{
case 1: Pokaz_macierz(mac1);
break;
case 2: Pokaz_macierz(mac2);
break;
case 3: Pokaz_macierz(mac_wynik);
break;
case 4:
gettimeofday(&tim, NULL); //funkcja pobierania czasu i obliczania ile program pracuje
tstart=tim.tv_sec+(tim.tv_usec/1000000.0);
Do_Everything(2); //parametr nr2
Zrob_macierz(wynik,mac_wynik);
gettimeofday(&tim, NULL);
tend=tim.tv_sec+(tim.tv_usec/1000000.0);
printf("Czas wykonania: %.2lf s\n", tend-tstart);
break;
case 9:
printf("\n\nAktualna liczba watkow SPE to : %d",MAX_NUM_SPE_THREADS);
printf("\nPodaj nowa wartosc : ");scanf("%d",&MAX_NUM_SPE_THREADS);
break;
case 0: printf("Koniec Programu");
break;
}
} while (wybor!=0);
printf("\nThe End\n");
return (0);
}
// *** FUNKCJE POMOCNICZE *** //
//Funkcja, ktora wywolujemy z parametrem, ktory jest numerem operacji, ktora chcemy wykonac , tylko dodawanie (czyli nr4)
//tworzy ona watki spe i wywoluje program (czyli dma_example_spu.c) dla kazdego watku
void Do_Everything(unsigned int number)
{
int i, num_spe_threads, num_elements_per_spe, rc;
spe_context_ptr_t spe_contexts[MAX_NUM_SPE_THREADS];
pthread_t spe_threads[MAX_NUM_SPE_THREADS];
control_block_t* control_blocks[MAX_NUM_SPE_THREADS];
ppu_thread_data_t ppu_thread_data[MAX_NUM_SPE_THREADS];
union {
unsigned long long ull;
unsigned int ui[2];
} long_addr;
num_spe_threads = MAX_NUM_SPE_THREADS;
//liczba podzielna przez 4096, jesli nie jest podzielna to zostanie reszta, ktora jesto bliczana na koncu funkcji main
num_elements_per_spe = ((((int)(X*Y/num_spe_threads))/CHUNK_SIZE)*CHUNK_SIZE);
printf("\nLiczba elementow na 1 spe : %.2d\n", num_elements_per_spe);
int reszta_elem = X*Y - num_elements_per_spe*num_spe_threads;
for (i = 0; i < num_spe_threads; i++)
{
// create control block
rc = posix_memalign ((void**)(&control_blocks[i]), 128, sizeof (control_block_t));
if (rc != 0)
{
fprintf (stderr, "Failed allocating space for control block\n");
exit (1);
}
//przekazujemy adresy macierzy do watku
long_addr.ui[0] = 0;
long_addr.ui[1] = (unsigned int)macierz1 + i* num_elements_per_spe*sizeof(float);
control_blocks[i]->in_addr1 = long_addr.ull;
long_addr.ui[0] = 0;
long_addr.ui[1] = (unsigned int)macierz2 + i* num_elements_per_spe*sizeof(float);
control_blocks[i]->in_addr2 = long_addr.ull;
long_addr.ui[0] = 0;
long_addr.ui[1] = (unsigned int)wynik + i* num_elements_per_spe*sizeof(float);
control_blocks[i]->out_addr = long_addr.ull;
control_blocks[i]->num_elements_per_spe = num_elements_per_spe;
control_blocks[i]->id = i;
control_blocks[i]->id_task = number; //to jest nasz parametr funkcji
// Create SPE threads to execute 'dma_example_spu'.
spe_contexts[i]= spe_context_create (0, NULL);
if (spe_contexts[i] == NULL)
{
perror ("Failed creating SPE context");
exit (1);
}
// Load SPE program into context
if (spe_program_load (spe_contexts[i], &dma_example_spu))
{
perror ("Failed loading SPE program");
exit (1);
}
// load ppu_thread_data
ppu_thread_data[i].spe_context = spe_contexts[i];
ppu_thread_data[i].control_block_ptr = control_blocks[i];
// Create SPE thread
if (pthread_create (&spe_threads[i], NULL, &ppu_pthread_function, &ppu_thread_data[i]))
{
perror ("Failed creating SPE thread");
exit (1);
}
// SPE is executing
}
//zwalnienie pamieci dla poszczegolnych spe ta petla po zakonczeniu dzialania watkow kasuje je z pamieci
for (i = 0; i < num_spe_threads; i++)
{
// Wait for SPE thread to complete execution
if (pthread_join (spe_threads[i], NULL))
{
perror ("Failed pthread_join");
exit (1);
}
// Destroy context
if (spe_context_destroy (spe_contexts[i]) != 0)
{
perror ("Failed destroying SPE context");
exit (1);
}
free (control_blocks[i]);
}
//jak zostaly jakies nie obliczone dane (reszta z dzielenia przez 4096) to obliczenia sa dokanczane tutaj.
if (reszta_elem != 0 )
if(number == 2) //dodawanie
for (int i=X*Y - reszta_elem; i<X*Y; i++)
wynik[i] = macierz1[i]+macierz2[i];
printf("\nzakonczono\n");
}
//pokazuje wycinke macierzy 10x10 elementow, to tak dla nas zeby zobaczyc, czy cos sie dzieje
void Pokaz_macierz(float **tabl)
{
printf("\n\n");
for (int i=0 ; i<10 ;i++)
{
printf("\n");
for (int j=0; j<10; j++)
{
printf("%4.0f ",(tabl[i][j]));
}
}
printf("\nOd konca\n");
for (int i=X-10 ; i<X ;i++)
{
printf("\n");
for (int j=X-10; j<X; j++)
{
printf("%4.0f ",(tabl[i][j]));
}
}
}
//Tworzy macierz z ciaglego obszaru pamieci
void Zrob_macierz(float * mac, float **macierz)
{
for(int i = 0 ; i<X; i++)
for (int j=0; j<Y; j++)
macierz[i][j] = *(mac+X*i+j);
}
//twrzy ciagly obszar pamieci z macierzy
void Przeksztalc_macierz(float **macierz, float *temp_macierz)
{
for (int i=0 ; i<X ;i++)
for (int j=0; j<Y; j++)
*(temp_macierz+i*X+j) = macierz[i][j];
}plik dma_example_spu.c:
#include <stdlib.h>
#include <stdio.h>
#include <spu_mfcio.h>
#include <single_buffer.h>
//#define USE_TIMER 1
#ifdef USE_TIMER
#include <spu_timer.h>
#endif /* USE_TIMER */
/* here we allocate the data buffer in local memory to hold the input
* data as well as the resulting output data
*/
float local_buffer1[CHUNK_SIZE] __attribute__ ((aligned (128)));
float local_buffer2[CHUNK_SIZE] __attribute__ ((aligned (128)));
float local_buffer3[CHUNK_SIZE] __attribute__ ((aligned (128)));
/* allocating the control_block structure on the stack. We align it by */
control_block_t control_block __attribute__ ((aligned (128)));
void dodawanie_macierzy (float* buf1, float* buf2, float* buf3)
{
unsigned int i;
vector float* vbuf_in1; //wektor float ma 4 bajty
vector float* vbuf_in2;
vector float* vbuf_out;
//ustawienie wektorow na tablice z argumentow
vbuf_in1 = (vector float*)buf1;
vbuf_in2 = (vector float*)buf2;
vbuf_out = (vector float*)buf3;
vector float maxPixel = (vector float) {255,255,255,255};
for(i=0 ; i<CHUNK_SIZE/4; i++)
{
vbuf_out[i] = vbuf_in2[i] + vbuf_in2[i]; //dodawnie wektorow
}
/*
//wersja bez wektoryzacji
for(i=0 ; i<CHUNK_SIZE; i++)
{
*(buf3+i) = *(buf1+i) + *(buf2+i); //dodawnie wektorow
}
*/
}
int main(unsigned long long speid __attribute__ ((unused)),
unsigned long long argp,
unsigned long long envp __attribute__ ((unused)))
{
unsigned int tag;
unsigned long long in_addr1,in_addr2, out_addr;
int i, num_chunks;
#ifdef USE_TIMER
uint64_t start, time_working;
spu_slih_register (MFC_DECREMENTER_EVENT, spu_clock_slih);
spu_clock_start();
start = spu_clock_read();
#endif /* USE_TIMER */
/* First, we reserve a MFC tag for use */
tag = mfc_tag_reserve();
if (tag == MFC_TAG_INVALID)
{
printf ("SPU ERROR, unable to reserve tag\n");
return 1;
}
mfc_get (&control_block, argp, sizeof (control_block_t), tag, 0, 0);
mfc_write_tag_mask (1 << tag);
mfc_read_tag_status_all ();
num_chunks = control_block.num_elements_per_spe/CHUNK_SIZE;
printf("\nSPU %d Work",control_block.id);
if (control_block.id_task == 2)
{
for (i = 0; i < num_chunks; i++)
{
in_addr1 = control_block.in_addr1 + (i* CHUNK_SIZE * sizeof(float));
in_addr2 = control_block.in_addr2 + (i* CHUNK_SIZE * sizeof(float));
out_addr = control_block.out_addr + (i * CHUNK_SIZE * sizeof(float));
mfc_get (local_buffer1, in_addr1, CHUNK_SIZE * sizeof(float), tag, 0, 0);
mfc_write_tag_mask (1 << tag);
mfc_read_tag_status_all ();
mfc_get (local_buffer2, in_addr2, CHUNK_SIZE * sizeof(float), tag, 0, 0);
mfc_write_tag_mask (1 << tag);
mfc_read_tag_status_all ();
dodawanie_macierzy (local_buffer1,local_buffer2,local_buffer3);
mfc_put (local_buffer3, out_addr, CHUNK_SIZE * sizeof(float), tag, 0, 0);
mfc_write_tag_mask (1 << tag);
mfc_read_tag_status_all ();
}
}
#ifdef USE_TIMER
time_working = (spu_clock_read() - start);
spu_clock_stop();
printf ("SPE time_working = %lld\n", time_working);
#endif
return 0;
}Pomiary szybkości przetwarzania od wielkości macierzy wykazują, że wraz ze wzrostem rozmiaru macierzy czas przetwarzania wzrasta, a czas wykonania sekwencyjnego jest średnio o połowę większy niż czas wykonywania przez 6 wątków.
Przy zastosowaniu wektoryzacji zauważono, że program sekwencyjny wykonuje się 1/3 szybciej, niż w przypadku bez wektoryzacji, natomiast różnica czasu wykonania przy 6 wątkach jest nie zauażalna pomiędzy programem z wektoryzacją i bez wektoryzacji. Przy przetwarzaniu z wektoryzacją zauważono mniejsze różnice w czasie przetwarzania, niż w przypadku bez wektoryzacji.
Wykresy
Podczas badania zależności pomiędzy poszczególnymi czasami wykonania mierzono czas wykonania dla 1,2,3,4,5,6 wątków spe, dodatkowo zmianiano rozmiar macierzy od 2500*2500 elementów do 3000*3000.
Wykresy: