je to knihovna, pro paralelizace - primarne urcena pro matematicke vypocty... navrzena pro C/C++ a Fortran !
Ja se budu zabyvat tou C++ casti...
zakladni myslenka:
Takze na cem si to ukazeme ... co treba rovnice ?
mejme soustavu mnoha rovnic o mnoha neznamych. Idealne pro pocitacove zpracovani zapsano v maticich !
Takže co budeme potřebovat - pár podpůrných funkcí:
- Kód: Vybrat vše
#include <iostream>
#include <fstream>
#include <windows.h>
//#include <ctime>
#include <omp.h>
using namespace std;
//funkce pro nacteni matice ze souboru filename
void nacti(double** pole, int& radku, int& sloupcu, char* filename) {
ifstream ifs(filename);
for (int i = 0; i < radku; i++) {
for (int j = 0; j < sloupcu; j++) {
ifs >> pole[i][j];
}
}
}
// vypis matice
//na std vystup
// matice pole, pocet radku rad, pocet sloupcu slo
void vypis(double** pole,int rad,int slo) {
for (int i = 0; i < rad; i++) {
for (int j = 0; j < slo; j++) {
cout << pole[i][j] << " ";
}
cout << endl;
}
cout << endl;
}
// funkce pro prohozeni radku a za b
void prohodRadky(double** pole,int rad, int slo, int a, int b) {
//cout << "menim radky " << a << " a " << b << endl;
double pom = 0.0;
for (int i = 0; i < slo; i++) {
pom = pole[a][i];
pole[a][i] = pole[b][i];
pole[b][i] = pom;
}
}
//funkce pro odecteni radku ktery od radku od od pozice poz
void odecti(double** pole,int rad, int slo, int od, int ktery, int poz) {
//cout << "odecitam od radku " << od << " radek " << ktery << endl;
double pomoc = (pole[od][poz] / pole[ktery][poz]);
for (int i = poz; i < slo; i++) {
pole[od][i] -= pole[ktery][i]*pomoc;
}
}
//funkce pro ziskani kopie matice
double** copyMat(double** zdroj, int radku, int sloupcu) {
double** pole = new double*[radku];
for (int i = 0; i < radku; i++) {
pole[i] = new double[sloupcu];
}
for (int i = 0; i < radku; i++) {
for (int j = 0; j < sloupcu; j++) {
pole[i][j] = zdroj[i][j];
}
}
return pole;
}
Takže teď neparalelní verze algoritmu:
- Kód: Vybrat vše
//neparalelni verze
long gaussSer(double** k, int radku, int sloupcu) {
//vypis(k,radku,sloupcu);
SYSTEMTIME st, lt;
GetSystemTime(&st);
int pozice = 0;
for (int i = 0; i < radku; i++) {
if (pozice > sloupcu) break;
if (k[i][pozice] == 0) {
for (int j = i+1; j < radku; j++) {
if (k[j][pozice] != 0) {
prohodRadky(k,radku, sloupcu,i,j);
//vypis(k, radku, sloupcu);
break;
}
}
if (k[i][pozice] == 0) {
pozice++;
i--;
continue;
}
}
double pom = k[i][pozice];
for (int j = pozice; j < sloupcu; j++) {
k[i][j] /= pom;
}
//vypis(k,radku,sloupcu);
for (int j = 0; j < radku; j++) {
if (j!=i) {
odecti(k,radku,sloupcu,j,i,pozice);
//vypis(k, radku, sloupcu);
}
}
pozice++;
}
GetSystemTime(<);
int a = st.wMinute*60000 + st.wSecond * 1000 + st.wMilliseconds;
int b = lt.wMinute*60000 + lt.wSecond * 1000 + lt.wMilliseconds;
//cout << d << endl;
//vypis(k,radku,sloupcu);
return b-a;
}
jenom doplnim ze tu pomoci WINAPI merim cas, a nevracim samotnou matici potreboval jsem jen rychlost, ne primo kolik to vyjde
tak ted moznost jedna, zparalelizujeme for cykly, a pripadne pouzijeme direktivu pro kritickou sekci!
- Kód: Vybrat vše
//paralelizace pomoci parallel for a critical
long gaussPar1(double** pole, int radku, int sloupcu) {
omp_set_dynamic(0);
//vypis(pole,radku,sloupcu);
SYSTEMTIME st, lt;
GetSystemTime(&st);
//cout << omp_get_wtime() << " " << omp_get_wtick() << endl;
int radekKVymene = 0;
int pozice = 0;
for (int i = 0; i < radku; i++) {
//vypis(pole,radku,sloupcu);
if (pozice > sloupcu) break;
if (pole[i][pozice] == 0) {
radekKVymene = i;
omp_set_num_threads(4);
#pragma omp parallel for schedule(dynamic,1)
for (int j = i+1; j < radku; j++) {
if (pole[j][pozice] != 0) {
#pragma omp critical (vymena)
{
radekKVymene = j;
}
}
}
if (radekKVymene == i) {
pozice++;
i--;
continue;
}
prohodRadky(pole,radku,sloupcu,i,radekKVymene);
}
double pom = pole[i][pozice];
omp_set_num_threads(4);
#pragma omp parallel for schedule(dynamic,1)
for (int j = pozice; j < sloupcu; j++) {
pole[i][j] /= pom;
}
//vypis(k,radku,sloupcu);
#pragma omp parallel for
for (int j = 0; j < radku; j++) {
if (j!=i ) {
//cout << omp_get_thread_num() << endl;
odecti(pole,radku,sloupcu,j,i,pozice);
//vypis(k, radku, sloupcu);
}
}
pozice++;
}
GetSystemTime(<);
long a = st.wMinute*60000 + st.wSecond * 1000 + st.wMilliseconds;
long b = lt.wMinute*60000 + lt.wSecond * 1000 + lt.wMilliseconds;
//cout << d << endl;
//vypis(pole,radku,sloupcu);
return b-a;
}
a ted moznost 2 - pomoci direktiv pro sekce, pripadne funkci pro zamek:
- Kód: Vybrat vše
//paralelizace pomoci parallel sections a funkci pro zamek
long gaussPar2(double** pole, int radku, int sloupcu) {
omp_set_dynamic(0);
//vypis(pole,radku,sloupcu);
SYSTEMTIME st, lt;
GetSystemTime(&st);
omp_lock_t * Lock = new omp_lock_t;
omp_init_lock(Lock);
int radekKVymene = 0;
int pozice = 0;
int j = 0;
for (int i = 0; i < radku; i++) {
//vypis(pole,radku,sloupcu);
if (pozice > sloupcu) break;
if (pole[i][pozice] == 0) {
radekKVymene = i;
omp_set_num_threads(4);
#pragma omp parallel sections private(j) shared(Lock)
{
#pragma omp section
{
for (j = i+1; j < radku; j = j + 4) {
if (pole[j][pozice] != 0) {
omp_set_lock(Lock);
radekKVymene = j;
omp_unset_lock(Lock);
break;
}
}
}
#pragma omp section
{
for (j = i+2; j < radku; j = j + 4) {
if (pole[j][pozice] != 0) {
omp_set_lock(Lock);
radekKVymene = j;
omp_unset_lock(Lock);
break;
}
}
}
#pragma omp section
{
for (j = i+3; j < radku; j = j + 4) {
if (pole[j][pozice] != 0) {
omp_set_lock(Lock);
radekKVymene = j;
omp_unset_lock(Lock);
break;
}
}
}
#pragma omp section
{
for (j = i+4; j < radku; j = j + 4) {
if (pole[j][pozice] != 0) {
omp_set_lock(Lock);
radekKVymene = j;
omp_unset_lock(Lock);
break;
}
}
}
}
if (radekKVymene == i) {
pozice++;
i--;
continue;
}
prohodRadky(pole,radku,sloupcu,i,radekKVymene);
}
double pom = pole[i][pozice];
omp_set_num_threads(4);
#pragma omp parallel sections private(j)
{
#pragma omp section
{
for (j = pozice; j < sloupcu; j=j+4) {
pole[i][j] /= pom;
}
}
#pragma omp section
{
for (j = pozice+1; j < sloupcu; j=j+4) {
pole[i][j] /= pom;
}
}
#pragma omp section
{
for (j = pozice+2; j < sloupcu; j=j+4) {
pole[i][j] /= pom;
}
}
#pragma omp section
{
for (j = pozice+3; j < sloupcu; j=j+4) {
pole[i][j] /= pom;
}
}
}
//vypis(k,radku,sloupcu);
omp_set_num_threads(4);
#pragma omp parallel sections private(j)
{
#pragma omp section
{
for (j = 0; j < sloupcu; j=j+4) {
if (j!=i) {
odecti(pole,radku,sloupcu,j,i,pozice);
}
}
}
#pragma omp section
{
for (j = 1; j < sloupcu; j=j+4) {
if (j!=i) {
odecti(pole,radku,sloupcu,j,i,pozice);
}
}
}
#pragma omp section
{
for (j = 2; j < sloupcu; j=j+4) {
if (j!=i) {
odecti(pole,radku,sloupcu,j,i,pozice);
}
}
}
#pragma omp section
{
for (j = 3; j < sloupcu; j=j+4) {
if (j!=i) {
odecti(pole,radku,sloupcu,j,i,pozice);
}
}
}
}
pozice++;
}
GetSystemTime(<);
long a = st.wMinute*60000 + st.wSecond * 1000 + st.wMilliseconds;
long b = lt.wMinute*60000 +lt.wSecond * 1000 + lt.wMilliseconds;
//cout << d << endl;
//vypis(pole,radku,sloupcu);
omp_destroy_lock(Lock);
delete Lock;
return b-a;
}
Tak tolik zatim k openmp
P.S. je potreba pouzit kompilator co openmp umi, napriklad g++ od nejake verze tusim 4.4 ? ... ma parametr -fopenmp
http://openmp.org/wp/openmp-compilers/ - visual studio v EE verzi toto nepodporuje !
Ariczek