[C++] Optimalizace kódu (násobení matic)

C++, C#, Visual Basic, Delphi, Perl a ostatní

Moderátor: Moderátoři Živě.cz

Odeslat příspěvekod Prochycz 14. 4. 2016 22:11

Zdravím,

řeším problém s násobením matic. Mám následující vzorec
-2*M'*(setPoint-y-Mp*dup)
Rozměry matic:
M(x,y)
y,setPoint - skalár
Mp(y,z)
dup(z,1)
Což počítám následujícím způsobem:
Kód: Vybrat vše
   for (int j = 0; j < N-N1+1; j++)
   {
      temp = 0;
      for (int m = 0; m < horizontPrediction-1; m++)
      {
         temp = temp + Mp[j][m] * dup[m];
      }
      tempMatrix[j] = setPoint - y - temp;
   }
   for (int j = 0; j < horizontControl; j++)
   {
      temp = 0;
      for (int m = 0; m < N-N1+1; m++)
      {
         temp = temp + Mtranspose[j][m] * tempMatrix[m];
      }
      g[j] = temp*-2;
   }


Je možné to nějak líp napsat, aby to procesor líp přechroustal (rychleji) ? Např. při maticích:
Mt=15x615;
Mp=615x599;
dup=599x1;
tak výpočet trvá přibližně 10ms. Sám nevím, jestli to je či není pomalé, proto bych poprosil o Váš názor. Rád bych to stáhnul ještě o polovinu, jestli to je možné.
Vím, že pro násobení matic je výhodne použít Strassenův algoritmus, ale ten se vyplatí při větší maticích, navíc bych je musel doplnit na čtvercový matice, což si myslím, že taky v mém případě není ideální.

Děkuji za případné rady

Edit: Zapomněl jsem dodat sestavu:
i7-4702MQ 2,20GHz
8GB RAM
Win 8.1
Program překládán ve Visual Studiu.

Edit2:
Ještě dodám, že ve vzorci -2*M'*(setPoint-y-Mp*dup) zůstává matice M a Mp konstantní.
Prochycz
Junior

Odeslat příspěvekod Nargon 15. 4. 2016 00:00

Můžeš použít ukazatele. Nevím jak dnes pokročily optimalizace v překladači ale z toho co si pamatuji tak v tomto stavu tam probíhá spousta zbytečných matematických operací jen pro určení které dvě proměnné se budou násobit. Takže určitě nahradit to [m] ukazatelem, tím by sis mohl docela pomoci. A i ostatní hodnoty, které jsou stále stejné nemá cenu počítat znovu a znovu, takže si je uložit do proměnné a tu pak jen použít.
Ten kód ber s rezervou, V c++ nedělám a ukazatele jsem nepoužil několik let tak si nejsem jistý jestli to píšu správně. A netuším jaký datový typ mají tvoje matice, takže tipuju int.
Kód: Vybrat vše
   int tmp = setPoint - y;
   int tmpn = N-N1+1;
   for (int j = 0; j < tmpn; j++)
   {
      temp = 0;
      int* pMpj = &Mp[j][0];
      int* pdup = &dup[0];
      for (int m = 0; m < horizontPrediction-1; m++)
      {
         temp = temp + (*pMpj++) * (*pdup++);
      }
      tempMatrix[j] = tmp - temp;
   }
   for (int j = 0; j < horizontControl; j++)
   {
      temp = 0;
      int* pMtransposej = &Mtranspose[j][0];
      int* ptempMatrix = &tempMatrix[0];
      for (int m = 0; m < tmpn; m++)
      {
         temp = temp + (*pMtransposej++) * (ptempMatrix++);
      }
      g[j] = temp*-2;
   }

Pokud jsem se s tím datovým typem int trefil tak možná ještě ten poslední výpočet násobení by se dal nahradit bitovým posunem, který bude snad rychlejší než to násobení.
Kód: Vybrat vše
      g[j] = -(temp<<1);


Ovšem to píšu z hlavy a v c++ nedělám takže jestli tam jsou nějaké chyby tak si je musíš opravit sám.
Desktop: Ryzen 7 1800X (3.95GHz, 1.35V), Asus Crosshair VI Hero, 16GB DDR4 Ram (3200MHz), 128GB SSD + 3TB HDD, Nvidia GTX 1080
Notebook: Asus UL50VT 15.6" (SU7300@1.7GHz, 4GB ram, 500GB HDD, Intel GMA 4500MHD + nVidia G210M, dlouha vydrz cca 7+ hod)
Nargon
Moderátor

Odeslat příspěvekod filipxsikora 15. 4. 2016 06:51

Optimalizace a automatický vektorizér kódu máš zaplý ? to by mohlo dost pomoct. Ale ani vektorizér kódu neumí zázraky, pak už jen inline assembly a využít AVX instrukce.

EDIT: A úplně jsem zapomněl na základní otázku, jak moc (jestli vůbec) máš kód paralelizovaný ?
filipxsikora
Junior

Odeslat příspěvekod borekz 15. 4. 2016 09:41

filipxsikora píše:Optimalizace a automatický vektorizér kódu máš zaplý ? to by mohlo dost pomoct. Ale ani vektorizér kódu neumí zázraky, pak už jen inline assembly a využít AVX instrukce.


Pochybuju, že nějaký kompilátor umí vektorizovat cyklus s jedním násobením. Tipuju, že by se to muselo přepsat nějak takto:
Kód: Vybrat vše
for (int j = 0; j < N-N1+1; j++)
   {
      temp = 0;
      auto *p1 = Mp[j]; // je to double ?
      auto *p2 = dup;
      int m = horizontPrediction - 1;
      for (; m >= 4; m -= 4)
      {
            // tento radek by kompilator mohl vektorizovat
            temp = temp + p1[0] * p2[0] + p1[1] * p2[1] + p1[2] * p2[2] + p1[3] * p2[3];
            p1 += 4; p2 += 4;
      }
      // zbytek, co nejde vektorizovat
      while (--m >= 0)
      {
            temp = temp + *p1++ * *p2++;
      }
      tempMatrix[j] = setPoint - y - temp;
   }
borekz
Junior

Odeslat příspěvekod Prochycz 15. 4. 2016 10:02

Děkuji všem za odpovědi, večer až se k tomu zase dostanu, tak na to opět mrknu.

Koukám, že jsem zapomněl dodat, že to je všechno v double. Při použití integeru by to bylo asi daleko rychlejší, si myslím.

filipxsikora píše:EDIT: A úplně jsem zapomněl na základní otázku, jak moc (jestli vůbec) máš kód paralelizovaný ?


Nad tím jsem přemýšlel, pak jsem na to zapomněl. Osobně s paralelizací nemám žádné zkušenosti. Je to možné udělat nějak automaticky? Pokud ne, rozumím tomu správně, že to násobení matic bych musel rozdělit např. do 4 vláken, kde 1 vlákno by počítalo 1/4 násobení matic, 2 vlákno 2/4 násobení matic atd. A pro vytvoření vlákna použít funkci např. CreateThread()?
Prochycz
Junior

Odeslat příspěvekod Ertefol 15. 4. 2016 10:07

Tak spočítá se počet dostupných vláken, tím se vydělí pole, provede se funkce a pak se to (už v jednom vlákně) zase spojí dohromady.
PC, Wii U, 3DS
"Slabé pohlaví je silnějším pohlavím z důvodu slabosti silnějšího pohlaví k slabšímu pohlaví." - Přísloví
Jakub Kovář, externí redaktor Games, LEVELu, Pevnosti
Ertefol
Junior
Uživatelský avatar

Odeslat příspěvekod Prochycz 15. 4. 2016 10:08

Jasný děkuji, večer otestuji.
Prochycz
Junior

Odeslat příspěvekod subdivider 15. 4. 2016 10:24

Este je tu moznost, ze vyuzijes nejaku uz existujucu kniznicu pre pracu s maticami, napr. MTL (Matrix Template Library)
http://www.osl.iu.edu/research/mtl/
Open-Source Edition je zdarma...
http://www.simunova.com/node/135
Supercomputing Edition s multi-threadingom je uz platena...
http://www.simunova.com/pmtl4
Lepsie povedane je draha ako svina, permanentna licencia stoji 5.000,00 €... To je zamerane na univerzitnych vyskumnikov...
subdivider
Junior

Odeslat příspěvekod borekz 15. 4. 2016 10:33

Prochycz píše:Koukám, že jsem zapomněl dodat, že to je všechno v double. Při použití integeru by to bylo asi daleko rychlejší, si myslím.

S tím moc nepočítej. To platilo možná na 80486 nebo na AMD. Intel Pentium (originál z roku 1993) násobí float, double i extended double za 3 strojové cykly s propustností 1 cyklu pokud následuje sčítání a s propustností 1 cyklů, pokud následuje násobení. V tvojem případě tedy je propustnost na Pentiu 1 cyklus, protože lze střídát násobení a sčítání. ALU stejného procesoru násobila integer cca 10 cyklů. V následujících procesorech Intel ten rozdíl postupně mizel, ale existuje dodnes.
Naposledy upravil borekz dne 15. 4. 2016 11:31, celkově upraveno 1
borekz
Junior

Odeslat příspěvekod Prochycz 15. 4. 2016 10:44

subdivider
Já prvně používal knihovnu eigen. Ale právě to bylo o hodně pomalejší než u toho mýho programu. Ale u placených knihoven to bude určitě jiný.
Prochycz
Junior

Odeslat příspěvekod filipxsikora 15. 4. 2016 15:12

Prochycz píše:Nad tím jsem přemýšlel, pak jsem na to zapomněl. Osobně s paralelizací nemám žádné zkušenosti. Je to možné udělat nějak automaticky? Pokud ne, rozumím tomu správně, že to násobení matic bych musel rozdělit např. do 4 vláken, kde 1 vlákno by počítalo 1/4 násobení matic, 2 vlákno 2/4 násobení matic atd. A pro vytvoření vlákna použít funkci např. CreateThread()?


Asi nejjednodušší způsob paralelizace je využití knihovny OpenMP, která je součástí Visual C++.
Více zde: https://msdn.microsoft.com/cs-cz/library/fw509c3b.aspx a https://msdn.microsoft.com/cs-cz/library/tt15eb9t.aspx
filipxsikora
Junior

Odeslat příspěvekod logout 15. 4. 2016 17:15

Pokud potřebuješ řešit matice, tak nepoužívej vlastní funkce ale optimalizovanej BLAS/LAPACK - implementace je daleko rychlejší než "naivní", paralelizace je zadarmo. Napsat dobrou knihovnu pro lingebru je práce na několik let.... pro tým lidí nebo pro génia.

Doporučuju OpenBLAS, je rychlej, zdarmo, ve většíně linuxovejch distribucí je přidanej, popř. jde snadno zkompilovat. Další z variant je Atlas. Pokud bys chtěl fakt rychlost, tak si sežeň MKL od Intelu - ta bejvá ještě o něco rychlejší než OpenBlas, ale nevím, jakou má teď intel licenční politiku a jestli je k dispozici zdarma.

Ty knihovny maj stanardizovaný rozhraní, nápověda k němu je tady:
http://www.netlib.org/blas/
Tebe bude zajímat funkce
DGEMM
http://www.netlib.org/lapack/explore-ht ... 12878226e1

-- 15. 4. 2016 16:18 --

PS: Daleko víc, než na tom, jak trvá která instrukce je u takovýchto úloh důležité, jak se využívá cache procesoru, aby se minimalizovalo čekání na paměť. A to je magie, kterou normální člověk fakt nedá.

Atlas se to snaží řešit tím, že si provede testy a podle toho si to zkompiluje - a je dost pomalej. OpenBlas ručně tunil jakýsi japonec Goto několik let a udělal daleko rychlejší kus kódu.
logout
Junior

Odeslat příspěvekod PiranhaGreg 15. 4. 2016 22:51

lol, koukám teď na zdrojáky OpenBLAS s takhle šílený optimalizace jsem ještě neviděl :shock: .
PiranhaGreg
Mírně pokročilý
Uživatelský avatar

Odeslat příspěvekod Prochycz 18. 4. 2016 19:28

Tak jsem testoval jak OpenBLAS, tak MKL (je zdrama pro univerzitní účely a nekomerční využití) a u testu násobení matice A(500x500)*B(500*500), tak mkl mi to zpočítalo za 20ms u OpenBLAS za cca 120ms. Ale možná nemám OpenBLAS optimálně nastavený, jelikož tady je uváděno, že OpenBLAS by měl být rychlejší (z 2013).

Ale nastal problém, jelikož používám trochu upravený Visual Studio, abych mohl vytvářet real-time aplikace, tak se mi nedaří použít tyto knihovny. Normálně mám nalinkovaný soubor libopenblas.lib, ale i tak mi poté překladač vyhazuje následující error:error LNK2019: unresolved external symbol cblas_dgemm referenced in function wmain
Může být ještě někde jinde něco špatně nastaveného, že mi to nechce načíst tu danou funkci?

Podotýkám, že pokud vytvořím klasický konzolový projekt, tak po nalinkování knihovny, funkce fungují správně. Používám pak ještě knihovnu qpOASES.lib a tam mi to funguje správně. Říkám si, jestli není problém, že k tomu (OpenBLAS) patří ještě libopenblas.dll soubor a možná s tím neumí upravený překladač pracovat, nebo mám prostě něco chybně nastaveného. :-(

Pokud se to nepovede zprovonit, je to možné ještě nějak obejít bez nalinkování toho souboru? Přemýšlel jsem, že bych tu funkci cblas_dgemm normálně při debugování vykrokoval, jestli by se mi to povedlo a bylo to možné. Pokud ne, tak budu muset zkusit tu paralelizaci, jestli to vylepší rychlost.

Děkuji
Prochycz
Junior

Odeslat příspěvekod logout 19. 4. 2016 14:11

Takhle velkej rozdíl je divnej. Kompiloval jsi openblas sám? Máš u něj zapnutej multithread (poznáš i podle vytížení CPU). U mkl imo bude defaultně zapnutá.

Ohledně nenajití knihovny, to bude nějakej problém linkování ve Visual Studiu, v tom Ti asi moc nepomůžu. Kouknul bych se, jestli ten výslednej soubor má v závislostech tu openblas.dll.

Taky jestlis kompiloval openblas sám, tak se kouknout, jestli je tam zapnutý cblas interface, nebo ještě líp se kouknout do tý dll, jestli je tam ta fce, co chybí.
logout
Junior

Další stránka

Kdo je online

Uživatelé procházející toto fórum: Žádní registrovaní uživatelé a 0 návštevníků