C++ - paralelizace pomocí OpenMP

C++ - paralelizace pomocí OpenMP

Příspěvekod Ariczek » 21 říjen 2011 12:25:25

Studijni materialy: http://openmp.org/wp/ https://computing.llnl.gov/tutorials/openMP/

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: Obrázek

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(&lt);
   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(&lt);
   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(&lt);
   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
Ariczek
 
Příspěvky: 51
Registrován: 21 říjen 2011 10:54:52

Re: C++ - paralelizace pomocí OpenMP

Příspěvekod Wlezley » 27 květen 2013 07:35:24

Dovolím si nesouhlasit s tím že to Visual Studio Express Edition nepodporuje, OpenMP tam jde narvat násilím. Stejně jako se tvrdilo, že v EE edicích nelze kompilovat x64 aplikace a já to normálně dělám. :)
Uživatelský avatar
Wlezley
 
Příspěvky: 316
Registrován: 24 září 2011 22:54:46
Bydliště: Plzeň
Projekt: Wlezley EU


Zpět na C/C++

Kdo je online

Uživatelé procházející toto fórum: Žádní registrovaní uživatelé a 1 návštěvník


cron