Visual Studio 2010 i C# 4.0 IV DIO: Parallelfx i rješavanje sistema linearnih jednačina


Najznačajnija novina u VS 2010 po mom mišljenju je ParallelFx biblioteka za paralelno programiranje koja iskorištava višejezgrene procesore. O ovoj novini pisao sam na nekoliko postova:
  1. ParallelFx nova generacija konkurentnog programiranja
  2. ParallelFx II dio-PLINQ
  3. Izračunavanje broja PI i implementacija prekoParallelFx
U ovom dijelu biće prezentirana primjena ParallelFx na rješavanje sistema linearnih jednačina. Više detalja oko same biblioteke možete pronaći na gornjim linkovima.
Uvod

Jako puno inženjerskih problema pri procesu rješavanja svodimo na rješavanje sistema linearnih jednačina, koje onda rješavamo jednom od standardnih matematičkih metoda. Rješavanje sistema linearnih jednačina zahtjeva znatno procesorsko vrijeme kao i znatnu količinu memorije za pohranjivanje. Kada rješavamo sisteme sa većim brojem nepoznatih 100, 1000, 10000 i više, vrijeme izvršavanja klasičnim metodama i algoritmima se smanjivalo se kako su se procesori razvijal. Obzirom da se takt procesora više ne može povećavati, a daljnje proširenje svodi se na višejezgrene procesore, naši algoritmi za rješavanje ostali su na istom nivou. Naime, svi ti algoritmi bazirani su na sekvencijalno načinu izvršavanja koda pa višejezgreni procesori nemaju efekta u rješavanju. Pojavom više jezgrenih procesora šire su se počeli razvijati algoritmi za paralelno procesuiranje, a tako i paralelni algoritmi za rješavanje sistema jednačina kao i drugih algoritama koji zahtjevaju značajno procesorsko vrijeme. Paralelni algoritmi za rješavanje određenih inženjerskih procesa pojavili su se početkom 90-tih godina i s njima su se koristile samo velike kompanije poput NASA-e i sličnih, a koji su imali znatne materijalne i financijske resurse. U tu svrhu povezivali su se hiljade procesora koji su onda paralelno izvršavali dijelove koda, a na kraju su se sublimirali i davali rješenja. Tada još nisu bili poznati višejezgreni procesori poput današnjih pa su se procesori povezivali klasično gdje su zauzimali znatno prostora.

servers

Slika 1. 1000-Pentium Beowulf-Style Cluster Computer (www.genetic-programming.org/)

Npr. John Koza tvorac genetskog programiranja, povezao je 1000 procesora gdje je vršio svoja ispitivanja i proračune sa paralelnim algoritmom genetskog programiranja (O genetskom programiranju već sam pisao ). Pojavom višejezgrenih procesora paralelno programiranje postalo je dostupno i običnim „smrtnicima“ a samim tim i šira upotreba paralelnog načina programiranja.

Gaussova metoda eliminacije

Jedna od najpopularnijih metoda rješavanja sistema linearnih jednačina je Gaussova metoda koja koristi princip eliminacije kojom se u konačnom broju iteracija na kraju dobija jedno rješenje, od kojeg se u povratnom smjeru dolazi do drugih rješenja. Linearni sistema od clip_image006 jednačina sa clip_image006[1] nepoznatih clip_image008 je skup jednačina clip_image010 slijedeće matrične forme:

clip_image012

clip_image014

clip_image016

Gdje su clip_image018i clip_image020 realni brojevi članovi matrice. Sistem zovemo homogenim ako su svi članovi matrice kolone slobodnih članova clip_image020[1] jednaki nuli; inače sistema zovemo nehomogeni. Korištenjem matrica i njihovih operacija, jednačine clip_image010[1] možemo posmatrati kao matrični jednačinu:

clip_image022

Gdje je matrica clip_image024 clip_image026 matrica

clip_image028,

dok su matrice kolona nepoznatih označavamo sa clip_image030 , a matricu kolona slobodnih članova označavamo sa clip_image032. Formalni zapisi ovih matrica dati su sljedećim izrazima respektivno:

clip_image034 ,

clip_image036.

Proširena matrica sistema clip_image038 je matrica reda clip_image040 koju definišemo kao:

clip_image042

Rješenje sistema linearnih jednačina je skup brojeva clip_image044 koji zadovoljavaju svih clip_image006[2] jednačina, a vektor rješenja sistema linearnih jednačina clip_image046 je vektor koji zadovoljava matričnu jednačinu. Način na koji se rješavaju sistemi linearnih jednačina Gaussovim sistemom eliminacije možete pogledati iz referenci koje su stavljene na kraju ovog članka. Da bi implementirali algoritma za rješavanje sistema linearnih jednačina potrebno je poznavati ovu metodu, te imati na umu procesorsko vrijeme i memoriju koju će algoritam zauzimati.

Algoritam Gausove metode eliminacije

Način na koji rješavamo sisteme jednačina „pješice“ najbolji je način kako da započnemo implementaciju algoritma. Za informaciju, na internetu postoje hiljade stranica na kojima se mogu naći algoritmi za rješavanje sistema linearnih jednačina u raznim programskim jezicima. S toga ova implementacija koja je data nije ništa revolucionarno.

Algoritam Gauss clip_image048 – proširena matrica sistema

ULAZ : Proširena matrica sistema koja definiše clip_image050 članova clip_image052, gdje je clip_image054, a clip_image056– predstavlja broj nepoznati. 1.

Nakon učitavanja matrice sistema potrebno je izvršiti triangulaciju matrice, gdje je po definiciji zadnjih član takve matrice rješenje clip_image058. 2. Obrnutim postupkom uz poznavanje rješenja clip_image058[1] dolazimo do rješenje clip_image060, i tako do rješenja clip_image062

IZLAZ : Vektor rješenja  clip_image064.gif. Izvorni kod ovog algoritma dat je u nastavku:

//Gaussova metoda eliminacije
for (int k = 0; k <= brKol - 2; k++)
 {
  // Prvo potrazimo najveci koeficijent u kolonama kao Pivot, zbog preciznosti rezultata
  double pivot = matSistema[k][k];
  for (int i = k + 1; i <= brVrsta - 1; i++)
   {
    if (pivot < Math.Abs(matSistema[i][k]))
    {
      double[] temp = matSistema[k];
      matSistema[k] = matSistema[i];
      matSistema[i] = temp;
    }
  }
 // Triangulacija
 for (int j = k + 1; j <= brVrsta - 1; j++)
  {
    m[j][k] = matSistema[j][k] / matSistema[k][k];
    for (int p = k; p <= brKol - 1; p++)
    {
      matSistema[j][p] = matSistema[j][p] - m[j][k] * matSistema[k][p];
    }
  }
 }
//
double[] rjesenje = new double[brKol - 1];
rjesenje[brKol - 2] = matSistema[brVrsta - 1][brKol - 1] / matSistema[brVrsta - 1][brKol - 2];
 for (int i = brVrsta - 2; i >= 0; i--)
 {
   double sum = 0.0;
   for (int j = brKol - 2; j >= i + 1; j--)
    {
      sum = sum + matSistema[i][j] * rjesenje[j];
    }
   rjesenje[i] = (1 / matSistema[i][i]) * (matSistema[i][brKol - 1] - sum);
 }
return rjesenje;

Kada smo implementirali algoritam za rješavanje sistema linearnih jednačina potrebno je sada analizirati kod i vidjeti gdje bi se mogla izvršiti paralelizacija da bi naš algoritam bio brzi od konvencionalnog odnosno da bi mogao iskorištavati današnja višejezgrene procesore. Prije nego se upustimo u analizu koda potrebno je pretpostaviti da će naš algoritam rješavati više od 500 jednačina odnosno tražiti više od 500 nepoznati. Ako bi algoritam rješavao manje jednačina paralelizacija ne bi imala smisla jer se vrijeme izvršavanja sekvencijalnog koda bilo brže od paralelnog. Ako pogledamo gornju implementaciju mjesta gdje bi paralelizirali kod je samo na dijelovima pri triangulaciji matrice i izračunavanju rješenja. Na ova dva mjesta imamo ugnježdene petlje, a one u algoritmu troše najviše procesorskog vremena.

Kako algoritam počinje…

On uzima prvu kolonu i provjerava da li u danoj koloni prvi član po apsolutnoj vrijednosti je najveći član od svih desnih kolona. Ako to slučaj nije tada se traži kolona u kojoj je prvi član najveći i te dvije kolone zamjenjuju mjesta. To možemo slobodno uraditi zbog osobina matrica. Kada se normalizira kolona tada se vrši triangulacija odnosno za svaku vrstu ispod posmatrane elementi kolone moraju biti nule (Gaussov postupak eliminacije). U ovom koraku se koriste dvije ugnježdene petlje (crveno obojen kod) kojim se ovo ostvaruje. Ovaj dio kod je najpogodniji za paralelnu implementaciju. Razlog zašto se ne ide na paralelnu implementacije prve petlje je taj što se mora ispoštovati poredak kolona, a u paralelnoj implementaciji to ne možemo postići jednostavno. Drugi dio koda pri kojem se na osnovu triangulacije matrice izračunavaju rješenja također koristi dvije ugnježdene petlje, a samim tim i dio koda koji bi mogao biti kandidat za paralelnu implementaciju. Ovaj dio koda također se mora izvršavati u striktno definisanom redu tj. prvo se izračunava zadnja vrsta, pa predzadnja i td. do prve vrste. Kada bi ovaj dio koda bio implementiran paralelno narušili bi poredak, a tako i rješenja ne bi bila izračunata tačno. Na osnovu analize paralelna implementacija kod izgleda na sljedeći način:

//Gaussova metoda eliminacije
for (int k = 0; k <= brKol - 2; k++)
{
  // Prvo potrazimo najveci koeficijent u kolonama kao Pivot, zbog preciznosti rezultata
  double pivot = matSistema[k][k];
  for (int i = k + 1; i <= brVrsta - 1; i++)
  {
    if (pivot < Math.Abs(matSistema[i][k]))
     {   double[] temp = matSistema[k];
       matSistema[k] = matSistema[i];
       matSistema[i] = temp;
     }
  }
 // Triangulacija
 Parallel.For(k + 1, brVrsta, j =>
 {
   m[j][k] = matSistema[j][k] / matSistema[k][k];
   for (int p = k; p <= brKol - 1; p++)
   {
     matSistema[j][p] = matSistema[j][p] - m[j][k] * matSistema[k][p];
   }
 }
             );
}
//
double[] rjesenje = new double[brKol - 1];
rjesenje[brKol - 2] = matSistema[brVrsta - 1][brKol - 1] / matSistema[brVrsta - 1][brKol - 2];
for (int i = brVrsta - 2; i >= 0; i--)
 {
   double sum = 0.0;
   for (int j = brKol - 2; j >= i + 1; j--)
    {
      sum = sum + matSistema[i][j] * rjesenje[j];
    }
   rjesenje[i] = (1 / matSistema[i][i]) * (matSistema[i][brKol - 1] - sum);
 }
 return rjesenje;

Vidimo da paralelna verzija algoritma se nije puno promijenila a što Parallelfx čini vrlo efikasnom bibliotekom.

Testiranje sekvencijalnog i paralelnog algoritma

Sada kada smo implementirali algoritam potrebno ga je i testirati na višejezgrenim PC mašinama da bi smo se uvjerili da naš paralelni algoritam brže radi od sekvencijalog. U tu svrhu generirano je 5 csv datoteka u kojim stoje proširene matrice sistema od 10, 50, 100, 500 i 1000 nepoznatih, kojim ćemo izvršiti testiranje algoritma. Sistemi su jednostavni čija su rješenja sve jedinice. Test aplikacija je koncipirana tako da učitava sisteme jednačina iz csv datoteka rješava sisteme koristeći sekvencijalni i paralelni algoritam te na konzolu izbacuje vremena koja su potrošena pri rješavanju. Test koji je prezentiran ovdje rađen je na Intel Core 2 Quad CPU Q6600 2,4 GZ:

CoreQuad

Rezultat testa je sljedeći:

restest

Iz testa možemo primijetiti da sistemi preko 100 nepoznati daju bolje rezultate kod paralelne implementacije i možemo smatrati okvirno da je to neka granica preko koje bi trebali koristiti paralelnu implementaciju algoritma. Kako broj nepoznatih raste tako raste razlika između paralelne i sekvencijalne implementacije algoritma. Kompletan izvorni kod algoritma i testa može se skinuti sa ovog linka. Projekt je rađen u Visual Studiu 2010 beta 1koji se može skinuti sa ovog linka.

Reference

  1. http://blogs.msdn.com/pfxteam/
  2. Ćamila Ljubović, Matematika, IPSvjetlost Sarajevo 1997
  3. http://geeks.netindonesia.net/blogs/norman

About Bahrudin Hrnjica

PhD in Mechanical Engineering, Microsoft MVP for .NET. Likes .NET, Math, Mechanical Engineering, Evolutionary Algorithms, Blogging.

Posted on 09/08/2009, in .NET, C# and tagged , , , . Bookmark the permalink. Leave a comment.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s