loop optimieren
-
pah wie fustrierend... das bringt alles gar nichts.
for(unsigned int atom_i = 0; atom_i < num_solvent_atoms; ++atom_i) { DEBUG(1, "atom i: " << atom_i); const double xi = (*(pos_i + atom_i))(0) + tx; __m128d reg_xi = _mm_set1_pd(xi); const double yi = (*(pos_i + atom_i))(1) + ty; __m128d reg_yi = _mm_set1_pd(yi); const double zi = (*(pos_i + atom_i))(2) + tz; __m128d reg_zi = _mm_set1_pd(zi); for(unsigned int atom_jn = 1; atom_jn < num_solvent_atoms; atom_jn += 2) { const unsigned int atom_j = atom_jn - 1; const unsigned int param = atom_i * 3 + atom_j; DEBUG(1, "param: " << param); DEBUG(1, "atom j: " << atom_j); __m128d tmp, reg; // for _mm_set_pd note that arguments are swapped!!! tmp = _mm_sub_pd(reg_xi, _mm_set_pd((*(pos_j+atom_jn))(0), (*(pos_j+atom_j))(0))); reg = _mm_mul_pd(tmp, tmp); _mm_store_pd(x, tmp); tmp = _mm_sub_pd(reg_yi, _mm_set_pd((*(pos_j+atom_jn))(1), (*(pos_j+atom_j))(1))); reg = _mm_add_pd(reg, _mm_mul_pd(tmp, tmp)); _mm_store_pd(y, tmp); tmp = _mm_sub_pd(reg_zi, _mm_set_pd((*(pos_j+atom_jn))(2), (*(pos_j+atom_j))(2))); reg = _mm_add_pd(reg, _mm_mul_pd(tmp, tmp)); _mm_store_pd(z, tmp); DEBUG(1, "0: x: " << x[0] << " y: " << y[0] << " z: " << z[0]); DEBUG(1, "1: x: " << x[1] << " y: " << y[1] << " z: " << z[1]); _mm_store_pd(r2, reg); tmp = _mm_div_pd(_mm_set1_pd(1.0), reg); _mm_store_pd(r2i, tmp); reg = _mm_sqrt_pd(tmp); _mm_store_pd(ri, reg); reg = _mm_mul_pd(tmp, tmp); reg = _mm_mul_pd(reg, tmp); _mm_store_pd(dist6i, reg); reg = _mm_mul_pd(reg, _mm_set_pd(pair_parameter[param+1].c12, pair_parameter[param].c12)); _mm_store_pd(dist6i_c12, reg); DEBUG(1, "0: r2: " << r2[0] << " r2i: " << r2i[0] << " ri: " << ri[0]); DEBUG(1, "0: dist6i: " << dist6i[0] << " dist6_c12: " << dist6i_c12[0]); DEBUG(1, "1: r2: " << r2[1] << " r2i: " << r2i[1] << " ri: " << ri[1]); DEBUG(1, "1: dist6i: " << dist6i[1] << " dist6_c12: " << dist6i_c12[1]); DEBUG(1, "0: r: " << sqrt(r2[0]) << " " << (sqrt(r2[0]) > 1.4 ? " cutoff! " : "")); DEBUG(1, "1: r: " << sqrt(r2[1]) << " " << (sqrt(r2[1]) > 1.4 ? " cutoff! " : "")); e_lj += (dist6i_c12[0] - pair_parameter[param].c6) * dist6i[0]; e_lj += (dist6i_c12[1] - pair_parameter[param+1].c6) * dist6i[1]; e_crf += pair_parameter[param].q * (ri[0] - m_crf_2cut3i * r2[0] - m_crf_cut); e_crf += pair_parameter[param+1].q * (ri[1] - m_crf_2cut3i * r2[1] - m_crf_cut); f[0] = (dist6i_c12[0] + dist6i_c12[0] - pair_parameter[param].c6) * 6.0 * dist6i[0] * r2i[0] + pair_parameter[param].q * (ri[0] * r2i[0] + m_crf_cut3i); f[1] = (dist6i_c12[1] + dist6i_c12[1] - pair_parameter[param+1].c6) * 6.0 * dist6i[1] * r2i[1] + pair_parameter[param+1].q * (ri[1] * r2i[1] + m_crf_cut3i); fx[0] = f[0] * x[0]; fx[1] = f[1] * x[1]; fy[0] = f[0] * y[0]; fy[1] = f[1] * y[1]; fz[0] = f[0] * z[0]; fz[1] = f[1] * z[1]; (*(force_i+atom_i))(0) += fx[0]; (*(force_i+atom_i))(0) += fx[1]; (*(force_j+atom_j))(0) -= fx[0]; (*(force_j+atom_jn))(0) -= fx[1]; (*(force_i+atom_i))(1) += fy[0]; (*(force_i+atom_i))(1) += fy[1]; (*(force_j+atom_j))(1) -= fy[0]; (*(force_j+atom_jn))(1) -= fy[1]; (*(force_i+atom_i))(2) += fz[0]; (*(force_i+atom_i))(2) += fz[1]; (*(force_j+atom_j))(2) -= fz[0]; (*(force_j+atom_jn))(2) -= fz[1]; storage.virial_tensor(0, 0) += x[0] * fx[0]; storage.virial_tensor(0, 1) += x[0] * fy[0]; storage.virial_tensor(0, 2) += x[0] * fz[0]; storage.virial_tensor(1, 0) += y[0] * fx[0]; storage.virial_tensor(1, 1) += y[0] * fy[0]; storage.virial_tensor(1, 2) += y[0] * fz[0]; storage.virial_tensor(2, 0) += z[0] * fx[0]; storage.virial_tensor(2, 1) += z[0] * fy[0]; storage.virial_tensor(2, 2) += z[0] * fz[0]; storage.virial_tensor(0, 0) += x[1] * fx[1]; storage.virial_tensor(0, 1) += x[1] * fy[1]; storage.virial_tensor(0, 2) += x[1] * fz[1]; storage.virial_tensor(1, 0) += y[1] * fx[1]; storage.virial_tensor(1, 1) += y[1] * fy[1]; storage.virial_tensor(1, 2) += y[1] * fz[1]; storage.virial_tensor(2, 0) += z[1] * fx[1]; storage.virial_tensor(2, 1) += z[1] * fy[1]; storage.virial_tensor(2, 2) += z[1] * fz[1]; } } // if num_solvent atoms odd only!!! if (num_solvent_atoms % 2) { const unsigned int atom_j = num_solvent_atoms - 1; for (unsigned int atom_i = 0; atom_i < num_solvent_atoms; ++atom_i) { const unsigned int param = atom_i * 3 + atom_j; const double x = (*(pos_i + atom_i))(0) - (*(pos_j + atom_j))(0) + tx; const double y = (*(pos_i + atom_i))(1) - (*(pos_j + atom_j))(1) + ty; const double z = (*(pos_i + atom_i))(2) - (*(pos_j + atom_j))(2) + tz; const double r2 = x * x + y * y + z * z; const double r2i = 1.0 / r2; const double ri = sqrt(r2i); const double dist6i = r2i * r2i * r2i; const double dist6i_c12 = pair_parameter[param].c12 * dist6i; e_lj += (dist6i_c12 - pair_parameter[param].c6) * dist6i; e_crf += pair_parameter[param].q * (ri - m_crf_2cut3i * r2 - m_crf_cut); const double f = (dist6i_c12 + dist6i_c12 - pair_parameter[param].c6) * 6.0 * dist6i * r2i + pair_parameter[param].q * (ri * r2i + m_crf_cut3i); const double fx = f * x; const double fy = f * y; const double fz = f * z; (*(force_i + atom_i))(0) += fx; (*(force_j + atom_j))(0) -= fx; (*(force_i + atom_i))(1) += fy; (*(force_j + atom_j))(1) -= fy; (*(force_i + atom_i))(2) += fz; (*(force_j + atom_j))(2) -= fz; storage.virial_tensor(0, 0) += x * fx; storage.virial_tensor(0, 1) += x * fy; storage.virial_tensor(0, 2) += x * fz; storage.virial_tensor(1, 0) += y * fx; storage.virial_tensor(1, 1) += y * fy; storage.virial_tensor(1, 2) += y * fz; storage.virial_tensor(2, 0) += z * fx; storage.virial_tensor(2, 1) += z * fy; storage.virial_tensor(2, 2) += z * fz; } }
-
Warum soll sich das so schlecht auf mehrere Cores aufteilen lassen?
-
Bezahl lieber Leute die das können
-
das lässt sich gut auf Cores aufteilen. Aber das ist nicht das Thema (siehe weiter oben)
Optimiser schrieb:
Bezahl lieber Leute die das können
Ich will das aber selbst lernen. Aber mach ruhig eine Offerte wenn du's so gut kannst.
-
Nephelauxetic schrieb:
Optimiser schrieb:
Bezahl lieber Leute die das können
Ich will das aber selbst lernen. Aber mach ruhig eine Offerte wenn du's so gut kannst.
Mit jeder Geschwindigkeitsverdoppelung doppelte Bezahlung, zur Basis 16 Euro.
-
Nephelauxetic schrieb:
Ich will nun aber die Performance auf einem einzigen Core verbessern, da der Code relativ schlecht hochskaliert (bis ca. 16 Prozessoren auf QsNET2) und wir die Performance brauchen.
Nephelauxetic schrieb:
das lässt sich gut auf Cores aufteilen.
Was den nun?
Wenn du es parallelisieren kannst, dann steck da arbeit rein. An der Stelle wo du jetzt rum machst wirst du nicht mehr viel rausholen.
-
freiefahrt schrieb:
Was den nun?
Wenn du es parallelisieren kannst, dann steck da arbeit rein. An der Stelle wo du jetzt rum machst wirst du nicht mehr viel rausholen.
Das eine schliesst das andere nicht aus. Nur weil du ineffizienten Code auf 10 Cores laufen lassen kannst heisst das noch nicht dass das sinnvoll ist. Insb. wenn man theoretisch die selbe Performance auf 5-8 Cores haben könnte. Stell dir mal vor die Benutzer eines grossen Computercluster würden alle so denken: so verschwendest du schnell mal gut 30-50% deines Clusters und das kostet Geld, Zeit und verschwendet Energie.
Zur Skalierung hab ich ja bereits was gesagt. Die könnte man schon verbessern: Domaindecomposition, Parallelisierung sämtlicher Algorithmen (Andahls Gesetz) und aktives Threadscheduling (Loadbalancing). Damit würde ich realistischerweise für den Standardalgorithmus auf 20-200 Cores kommen (je nach Algorithmus) - nach einem Entwicklungsaufwand von wahrscheinlich mehr als einem Jahr. Wir haben aber keinen Cluster bei dem man mehr als 64 Cores kriegen kann und bei 64 Cores gewinnt man da gar nichts weil die Zeit in der man nicht rechnet wartet man in der Warteschlange die für 64 Cores entsprechend lang ist und langsam rückt. Ich bin kein Freund von diesem Bluegene Ansatz...LG
Neph
-
Viele Cores verwendest du denn jetzt? Der Code da oben ist wohl kaum so ineffizient, das du da 50% rausholen kannst.
-
freiefahrt schrieb:
Viele Cores verwendest du denn jetzt? Der Code da oben ist wohl kaum so ineffizient, das du da 50% rausholen kannst.
Aber sicher, er sagt nichtmal an welcher Stelle die Geschwindigkeit verloren geht, nicht weshalb und wieviel vom theoretischen Maximum der CPU er momentan erreicht.
-
freiefahrt schrieb:
Viele Cores verwendest du denn jetzt? Der Code da oben ist wohl kaum so ineffizient, das du da 50% rausholen kannst.
Maximal 16 und das bringt dann einen Speedup von ca. 11. Bei anderen Algorithmen im Programm sogar nur 4-8.
Optimiser schrieb:
freiefahrt schrieb:
Viele Cores verwendest du denn jetzt? Der Code da oben ist wohl kaum so ineffizient, das du da 50% rausholen kannst.
Aber sicher, er sagt nichtmal an welcher Stelle die Geschwindigkeit verloren geht, nicht weshalb und wieviel vom theoretischen Maximum der CPU er momentan erreicht.
Ja das denke ich auch (auch wenn 50% vielleicht ein bisschen hoch gegriffen ist). Ich werd mal Montag nochmals Profilen und dann weiss ich dann wenigstens die Antwort auf "an welcher Stelle die Geschwindigkeit verloren geht". Dann schau ich mir mal noch hustbaers Vorschläge genauer an.
Das wird schon... irgendwann
LG
Neph
-
Nephelauxetic schrieb:
Das eine schliesst das andere nicht aus. Nur weil du ineffizienten Code auf 10 Cores laufen lassen kannst heisst das noch nicht dass das sinnvoll ist. Insb. wenn man theoretisch die selbe Performance auf 5-8 Cores haben könnte. Stell dir mal vor die Benutzer eines grossen Computercluster würden alle so denken: so verschwendest du schnell mal gut 30-50% deines Clusters und das kostet Geld, Zeit und verschwendet Energie.
Nur kannst du keine 30-50% aus dem Code herausholen. Du kannst maximal 2-3% rausholen und das waere schon gut.
Hast du die Idee mit den page faults/cache misses mal verfolgt? Du greifst auf viel speicher zu - ist das alles notwendig? wie sind die datenstrukturen aufgebaut?
Wir haben aber keinen Cluster bei dem man mehr als 64 Cores kriegen kann und bei 64 Cores gewinnt man da gar nichts weil die Zeit in der man nicht rechnet wartet man in der Warteschlange die für 64 Cores entsprechend lang ist und langsam rückt. Ich bin kein Freund von diesem Bluegene Ansatz...
Verstehe ich gerade nicht was du meinst...
Die 64 Cores lastest du aber schon zu 100% aus, oder?
-
Das einzige was bei dem Code noch groß was bringen könnte, wäre ein einfacheres und schneleres sqtr zu verwenden, aber dann wird die Sache wahrscheinlich zu ungenau
-
freiefahrt schrieb:
Das einzige was bei dem Code noch groß was bringen könnte, wäre ein einfacheres und schneleres sqtr zu verwenden, aber dann wird die Sache wahrscheinlich zu ungenau
Man könnte bei dem Code vielleicht noch mal versuchen, die ganzen inneren Schleifen loszuwerden. Wenn ich es richtig verstanden habe, gibt es irgendwo eine "große" Schleife, die irgendwo im 6 stelligen Bereich liegt und innerhalb dieser großen Schleife gibt es noch diese paar ineinander verschachtelten kleineren Schleifen. Wenn der Compiler vorhersehen könnte, wie viele Schritte die kleineren Schleifen bräuchten, würde er "loop-unrolling" machen. In der Regel bedeutet eine Schleife, dass eine bedingte Sprung-Anweisung generiert wird und die ist meistens schlecht für die Performanz (habe ich irgendwo gelesen), denn durch einen Sprung wird die CPU Pipeline "entleert". Bei einer großen Schleife sind diese verlorenen Takte nicht relevant. Bei kleinen Schleifen dagegen bemerkbar und man muss ja die ganzen Takte zum Auffüllen der Pipeline mit 6 stelliger Zahl multiplizieren - wenn ich es richtig sehe - sagen wir mal, Pi mal Daumen, 20 Takte zum Auffüllen der Pipeline mal 1E6 ergibt 20E6 verlorene Takte.
Ich weiss nicht, wie ich es am besten erklären könnte, hier ein künstliches Beispiel:for (i = 0; i < N; ++i) { for (j = 0; j < 8; ++j) { ... } }Die innere Schleife läuft zufällig bis 8, und das ist vielfaches von 2
Die innere Schleife lässt sich eliminieren:for (i = 0; i < (8 * N); ++i) { ...; ++j; j &= ~7; }Also ungefähr so, wenn ich da oben keinen Fehler gemacht habe... Ich muss ehrlich zugeben, ich habe keine Ahnung, wieviele Prozent es bringt und ob es sich wirklich lohnt, schließlich haben ja andere Leute 30 Jahre lang darüber nachgedacht, und diese Methode würde ja eine komplette Umdenke beim Algorithmus bedeuten. Und man müsste das Ganze ja noch debuggen können.
-
abc.w schrieb:
In der Regel bedeutet eine Schleife, dass eine bedingte Sprung-Anweisung generiert wird und die ist meistens schlecht für die Performanz (habe ich irgendwo gelesen), denn durch einen Sprung wird die CPU Pipeline "entleert".
Ne, die wird nur geleert, wenn die branch prediction falsch war, aber die soll bei Intel eigentlich schon ganz gut sein.
-
Shade Of Mine schrieb:
Nephelauxetic schrieb:
Das eine schliesst das andere nicht aus. Nur weil du ineffizienten Code auf 10 Cores laufen lassen kannst heisst das noch nicht dass das sinnvoll ist. Insb. wenn man theoretisch die selbe Performance auf 5-8 Cores haben könnte. Stell dir mal vor die Benutzer eines grossen Computercluster würden alle so denken: so verschwendest du schnell mal gut 30-50% deines Clusters und das kostet Geld, Zeit und verschwendet Energie.
Nur kannst du keine 30-50% aus dem Code herausholen. Du kannst maximal 2-3% rausholen und das waere schon gut.
Woher willst du das ohne profilen wissen? nur vom code her? sorry, nichts gegen dich persoenlich, aber euer "kann man doch"-"kann man nicht" geposte ist irgendwie fuer die katz!
-
rapso schrieb:
Woher willst du das ohne profilen wissen? nur vom code her? sorry, nichts gegen dich persoenlich, aber euer "kann man doch"-"kann man nicht" geposte ist irgendwie fuer die katz!
wissen kann ich es nicht, aber ich sehe einfach deutlich mehr potential es domain-maessig anzugehen.
du hast 1 sqrt und 1 division hier im code. der rest sind multiplikationen und additionen - da ist code maessig nicht viel zu machen (abgesehen von cache misses/page faults).
was dagegen sehr viel bringen wuerde, waere den code umzustrukturieren - sofern moeglich.
was dich vermutlich da verwirrt hat ist, dass ich zwischen code optimierungen und domain optimierungen trenne. ein
c=sqrt(a);
ist code maessig nicht optimierbar. domain maesig aber schon, wenn naemlich inv_sqrt reicht oder man gar mit a*a weiter reichnen kann ohne wurzel zu ziehen - das sind ebenfalls optimierungen - aber in der domain und nicht im code.
-
Ich muss gestehen, dass ich deinen calc Code nicht ganz blicke.
(*(pos_i+atom_i))(0)Ist das ein Funktionsaufruf der etwas macht oder nur die erste Koordinate liest?
Was vielleicht etwas bringt wäre:
void calc_pair(const unsigend int i, const unsigned int j, const unsigned int num) { for(unsigned int x = 0; x < num; ++x) { calc(i,j,x,x); for(unsigned int y = x+1; y < num; ++y) { calc(i,j,x,y); calc(i,j,y,x); } } }Vielleicht du kannst dir bei der Berechnung von calc(i,j,x,y) und calc(i,j,y,x) ein paar Zwischenwerte recyclen oder Symmetrien ausnutzen.
Du sagtest, dass num nur bis 6 gehen kann. Eventuell könntest etwas wie folgt versuchen:
void calc_pair_num6(const unsigend int i, const unsigned int j) { for(unsigned int x = 0; x < 6; ++x) { calc(i,j,x,0); calc(i,j,x,1); calc(i,j,x,2); calc(i,j,x,3); calc(i,j,x,4); calc(i,j,x,5); } }Vielleicht gibt es da wieder Sachen die doppel berechnet werden. Ist natürlich sehr aufwendig zu schreiben und nur schwer wartbar.
Vielleicht wartet auch einer von deinen Threads mehr als du erwartest und es würde was bringen den Code auf zwei Threads aufzuteilen damit der Scheduler was hat was er dem wartenden Core als Aufgabe geben kann.
Wenn du mehr Hilfe willst dann poste mal die gesamte Funktion (optimierte und nicht optimierte Version) damit man sieht was da eventuell mehrfach berechnet wird und am besten auch mit ein paar Kommentaren bzw mathematisch und/oder physikalischen Stichwörtern was du da machst.
Compiler und Prozessor anzugeben wäre auch sinnvoll. Wenn wirklich Micromircooptimierungen her müssen dann geht das nicht ohne das zu wissen.
-
Shade Of Mine schrieb:
rapso schrieb:
Woher willst du das ohne profilen wissen? nur vom code her? sorry, nichts gegen dich persoenlich, aber euer "kann man doch"-"kann man nicht" geposte ist irgendwie fuer die katz!
wissen kann ich es nicht, aber ich sehe einfach deutlich mehr potential es domain-maessig anzugehen.
du hast 1 sqrt und 1 division hier im code. der rest sind multiplikationen und additionen - da ist code maessig nicht viel zu machen (abgesehen von cache misses/page faults).
cache misses, Data hazards, pipeline switching. der code stallt vermutlich die meiste zeit ueber. ich wuerde nicht sagen dass da nicht viel zu machen ist.
was dagegen sehr viel bringen wuerde, waere den code umzustrukturieren - sofern moeglich.
was dich vermutlich da verwirrt hat ist, dass ich zwischen code optimierungen und domain optimierungen trenne. ein
c=sqrt(a);
ist code maessig nicht optimierbar. domain maesig aber schon, wenn naemlich inv_sqrt reicht oder man gar mit a*a weiter reichnen kann ohne wurzel zu ziehen - das sind ebenfalls optimierungen - aber in der domain und nicht im code.
ich denke schon dass das codemaessig optimierbar ist wenn du es im context siehst. die instruktion kann 5cycle oder ca 540cycles (auf nem AMD64) verbrauchen.
genau aus diesem grund sag ich die ganze zeit, das man sich mit nem profiler die kosten pro zeile anschauen muss. einfach nur aufgrund des sources kann man nur raten. klar findet mal ein blindes huhn das korn, aber es wird nicht wissen, ob und wieviele noch rumliegen und schon garnicht wo.er koennte den code auch posten und sagen 'er compiliert nicht', da wuerden sicher auch alle sagen, dass er den compiler output zeigen soll. ich wunder mich immer wieso das niemand fuers optimieren fragt, dass man ein wenig profiling daten sieht.
(und nein, das ist auch nichts gegen dich )
-
freiefahrt schrieb:
Das einzige was bei dem Code noch groß was bringen könnte, wäre ein einfacheres und schneleres sqtr zu verwenden, aber dann wird die Sache wahrscheinlich zu ungenau
Ja zu ungenau.
Shade Of Mine schrieb:
Hast du die Idee mit den page faults/cache misses mal verfolgt? Du greifst auf viel speicher zu - ist das alles notwendig? wie sind die datenstrukturen aufgebaut?
Notwendig ist es leider. Ich kann die Arrays auch nicht wirklich umsortieren damit sie cache optimaler waeren... Ich hab doch den valgrind output zu den cache misses gepostet (callgrind) aber ich finds sehr schwierig da was zu interpretieren...
Shade Of Mine schrieb:
Die 64 Cores lastest du aber schon zu 100% aus, oder?
Nein. Andahls Law und Load Balancing Probleme. Und die wird man so einfach auch nicht los.
abc.w schrieb:
Man könnte bei dem Code vielleicht noch mal versuchen, die ganzen inneren Schleifen loszuwerden.
Das ueberlasse ich dem Compiler. Der compiler optimiert mit den Profiling Informationen und kann also die Loops unrollen wenn er will. Bringen tut das aber nicht viel.
Ben04 schrieb:
Ich muss gestehen, dass ich deinen calc Code nicht ganz blicke.
(*(pos_i+atom_i))(0)Ist das ein Funktionsaufruf der etwas macht oder nur die erste Koordinate liest?
ist ein bisschen umstaendlich geschrieben aber es ist im Prinzip nur das lesen der ersten Koordinate. normal geht das so:
conf.current().pos(i)(0)}}}. Durch die Pointer sieht das dann halt ein bisschen haesslicher aus...
Ben04 schrieb:
Was vielleicht etwas bringt wäre:
void calc_pair(const unsigend int i, const unsigned int j, const unsigned int num) { for(unsigned int x = 0; x < num; ++x) { calc(i,j,x,x); for(unsigned int y = x+1; y < num; ++y) { calc(i,j,x,y); calc(i,j,y,x); } } }Vielleicht du kannst dir bei der Berechnung von calc(i,j,x,y) und calc(i,j,y,x) ein paar Zwischenwerte recyclen oder Symmetrien ausnutzen.
Nein leider kann man da keine Zwischenwerte recyclen. So hast du schlussendlich was aehnliches wie das von mir gepostete Beispiel... bringt nix.
Ben04 schrieb:
Du sagtest, dass num nur bis 6 gehen kann. Eventuell könntest etwas wie folgt versuchen:
void calc_pair_num6(const unsigend int i, const unsigned int j) { for(unsigned int x = 0; x < 6; ++x) { calc(i,j,x,0); calc(i,j,x,1); calc(i,j,x,2); calc(i,j,x,3); calc(i,j,x,4); calc(i,j,x,5); } }Vielleicht gibt es da wieder Sachen die doppel berechnet werden. Ist natürlich sehr aufwendig zu schreiben und nur schwer wartbar.
Ja so sehe ich das auch. Sowas will niemand warten muessen. Zudem kann das ganze bis 1000 gehen wenn jemand will... Also 6 war nur eine praktische Grenze weil mir grad kein gaengiges Loesungsmittel mit mehr als 6 Atomen einfiel. Aber im Prinzip waer das moeglich.
Ben04 schrieb:
Vielleicht wartet auch einer von deinen Threads mehr als du erwartest und es würde was bringen den Code auf zwei Threads aufzuteilen damit der Scheduler was hat was er dem wartenden Core als Aufgabe geben kann.
Ich weiss nicht genau wie die OpenMP Implementierung skaliert. Ich hab mich nur mit MPI beschaeftigt und da waer das nicht trivial zu implementieren. Aber daran gedacht hab ich auch schon
Ben04 schrieb:
Wenn du mehr Hilfe willst dann poste mal die gesamte Funktion (optimierte und nicht optimierte Version) damit man sieht was da eventuell mehrfach berechnet wird und am besten auch mit ein paar Kommentaren bzw mathematisch und/oder physikalischen Stichwörtern was du da machst.
Compiler und Prozessor anzugeben wäre auch sinnvoll. Wenn wirklich Micromircooptimierungen her müssen dann geht das nicht ohne das zu wissen.
g++ 4.2.3, Linux, -O3 -funroll-loops -fno-gcse
$ cat /proc/cpuinfo processor : 0 vendor_id : AuthenticAMD cpu family : 15 model : 67 model name : AMD Processor model unknown stepping : 3 cpu MHz : 1000.000 cache size : 1024 KB physical id : 0 siblings : 2 core id : 0 cpu cores : 2 fpu : yes fpu_exception : yes cpuid level : 1 wp : yes flags : fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush mmx fxsr sse sse2 ht syscall nx mmxext fxsr_opt rdtscp lm 3dnowext 3dnow rep_good pni cx16 lahf_lm cmp_legacy svm extapic cr8_legacy bogomips : 2011.93 TLB size : 1024 4K pages clflush size : 64 cache_alignment : 64 address sizes : 40 bits physical, 48 bits virtual power management: ts fid vid ttp tm stc processor : 1 vendor_id : AuthenticAMD cpu family : 15 model : 67 model name : AMD Processor model unknown stepping : 3 cpu MHz : 1000.000 cache size : 1024 KB physical id : 0 siblings : 2 core id : 1 cpu cores : 2 fpu : yes fpu_exception : yes cpuid level : 1 wp : yes flags : fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush mmx fxsr sse sse2 ht syscall nx mmxext fxsr_opt rdtscp lm 3dnowext 3dnow rep_good pni cx16 lahf_lm cmp_legacy svm extapic cr8_legacy bogomips : 2011.93 TLB size : 1024 4K pages clflush size : 64 cache_alignment : 64 address sizes : 40 bits physical, 48 bits virtual power management: ts fid vid ttp tm stcaufrufende Funktion (stub)
// rufe schnellen generischen lsm loop auf. if (sim.param().force.special_loop == simulation::special_loop_generic) { const unsigned int num_solvent_atoms = topo.solvent(0).num_atoms(); // Parameter vorberechnen (Lennard Jones Paare und Ladungsprodukt. // prepare parameters const unsigned int num_param = num_solvent_atoms * num_solvent_atoms; typename interaction::Nonbonded_Innerloop<t_interaction_spec>::solvent_pair_parameter pair_parameter[num_param]; unsigned int param = 0; // loop ueber atompaare for(unsigned int atom_i = 0; atom_i < num_solvent_atoms; ++atom_i) { for(unsigned int atom_j = 0; atom_j < num_solvent_atoms; ++atom_j, ++param) { assert(param < num_param); DEBUG(10, "\tsolvent pair parameter: " << param); // hole LJ Parameter (doubles) const lj_parameter_struct & lj = innerloop.param()->lj_parameter(topo.solvent(0).atom(atom_i).iac, topo.solvent(0).atom(atom_j).iac); pair_parameter[param].c12 = lj.c12; pair_parameter[param].c6 = lj.c6; // berechne Ladungsprodukt. pair_parameter[param].q = math::four_pi_eps_i * topo.solvent(0).atom(atom_i).charge * topo.solvent(0).atom(atom_j).charge; DEBUG(10, "\t\tc12: " << pair_parameter[param].c12 << " c6: " << pair_parameter[param].c6 << " q: " << pair_parameter[param].q); } } // use num_solvent_atoms-th atom (first of solvent molecule i) // loop ueber alle Molekuelpaare // liste enthaelt alle atompaare und deshalb werden immer num_solvent_atoms - 1 atome uebersprungen for (; i < size_i; i += num_solvent_atoms) { for (j_it = pairlist_solvent[i].begin(), j_to = pairlist_solvent[i].end(); j_it != j_to; j_it += num_solvent_atoms) { // use num_solvent_atoms-th atom (first of solvent molecule j) DEBUG(10, "\tsolvent_nonbonded_interaction: i " << i << " j " << *j_it); // Aufruf der zu optimierenden Funktion innerloop.solvent_innerloop(topo, pair_parameter, conf, num_solvent_atoms, i, *j_it, storage, periodicity); } } }Die zu optimierende Funktion:
template<typename t_nonbonded_spec> inline void interaction::Nonbonded_Innerloop<t_nonbonded_spec>::solvent_innerloop ( topology::Topology & topo, interaction::Nonbonded_Innerloop<t_nonbonded_spec>::solvent_pair_parameter * pair_parameter, configuration::Configuration & conf, const unsigned int num_solvent_atoms, const int i, const int j, Storage & storage, math::Periodicity<t_nonbonded_spec::boundary_type> const & periodicity ) { math::Vec r; // only one energy group // bestimme wo in array die Energie gespeichert werden soll const int egroup = topo.atom_energy_group(topo.num_solute_atoms()); DEBUG(8, "\tspc pair\t" << i << "\t" << j << " egroup " << egroup); // hole pointer auf position (const) und kraft von atom i math::Vec const * const pos_i = &conf.current().pos(i); math::Vec * const force_i = &storage.force(i); // dito fuer atom j math::Vec const * const pos_j = &conf.current().pos(j); math::Vec * const force_j = &storage.force(j); DEBUG(9, "i = " << i << " j = " << j); // first atom vs. first atom // berrechne naechstes periodisches Abild des Paars i, j (== Abstand zwischen atom i und j wird in vektor r gespeichert) periodicity.nearest_image(*pos_i, *pos_j, r); // berechne Periodizitaetsverschiebung fuer das Paar die dann auf // alle atom Paare der zwei Molekuele angewandt werden kann. const double tx = r(0) - (*pos_i)(0) + (*pos_j)(0); const double ty = r(1) - (*pos_i)(1) + (*pos_j)(1); const double tz = r(2) - (*pos_i)(2) + (*pos_j)(2); // speicher fuer die energie: Lennard-Jones und Coulombe double e_lj = 0.0, e_crf = 0.0; // loop over atom pairs for(unsigned int param = 0, atom_i = 0; atom_i < num_solvent_atoms; ++atom_i) { // berrechne verschiebung des atoms i const double xi = (*(pos_i+atom_i))(0) + tx; const double yi = (*(pos_i+atom_i))(1) + ty; const double zi = (*(pos_i+atom_i))(2) + tz; for(unsigned int atom_j = 0; atom_j < num_solvent_atoms; ++atom_j, ++param) { DEBUG(15, "\tatoms: i: " << atom_i << " j: " << atom_j); // berechne naechstes periodisches Abbild von i-j mit hilfe verschiebung von atom i const double x = xi - (*(pos_j+atom_j))(0); const double y = yi - (*(pos_j+atom_j))(1); const double z = zi - (*(pos_j+atom_j))(2); // rechne abstand im quadrat const double r2 = x * x + y * y + z * z; DEBUG(15, "\tr2: " << r2); assert(r2 != 0.0); // inverser abstand im quadrat const double r2i = 1.0 / r2; // inverser abstand const double ri = sqrt(r2i); // inverser abstand hoch 6 const double dist6i = r2i * r2i * r2i; // und multipliziert mit der konstante C12 const double dist6i_c12 = pair_parameter[param].c12 * dist6i; // berrechne (und addiere) LJ energie: (r^-6 * c12 - c6) * r^-6 e_lj += (dist6i_c12 - pair_parameter[param].c6) * dist6i; // berechne Coulombe und Fischer Reaction Field : // qi * qj / (4 Pi Eps0) * (r^-1 - b0 * r^2/cut^3 - cut^-1) e_crf += pair_parameter[param].q * (ri - m_crf_2cut3i * r2 - m_crf_cut); // Kraft: Ableitung nach r der beiden Energien const double f = (dist6i_c12 + dist6i_c12 - pair_parameter[param].c6) * 6.0 * dist6i * r2i + pair_parameter[param].q * (ri * r2i + m_crf_cut3i); // kraft als Vektor const double fx = f * x; const double fy = f * y; const double fz = f * z; // abspeichern der Kraft: actio = reactio (*(force_i+atom_i))(0) += fx; (*(force_j+atom_j))(0) -= fx; (*(force_i+atom_i))(1) += fy; (*(force_j+atom_j))(1) -= fy; (*(force_i+atom_i))(2) += fz; (*(force_j+atom_j))(2) -= fz; // berechnung des virial tensors: r (x) f // wobei (x) das inner product ist. storage.virial_tensor(0, 0) += x * fx; storage.virial_tensor(0, 1) += x * fy; storage.virial_tensor(0, 2) += x * fz; storage.virial_tensor(1, 0) += y * fx; storage.virial_tensor(1, 1) += y * fy; storage.virial_tensor(1, 2) += y * fz; storage.virial_tensor(2, 0) += z * fx; storage.virial_tensor(2, 1) += z * fy; storage.virial_tensor(2, 2) += z * fz; } } // abspeichern der Energien in den entsprechenden Energiegruppen. storage.energies.lj_energy[egroup][egroup] += e_lj; storage.energies.crf_energy[egroup][egroup] += e_crf; }Ich hoffe damit ist die Physik des Problems ein bisschen klarer...
LG
Neph
-
Hallo
Ich hab die cache misses ja gar nicht gepostet... wahrscheinlich hab ich nur von getraeumt
Events recorded abbreviations are: Ir : I cache reads (ie. instructions executed) I1mr: I1 cache read misses I2mr: L2 cache instruction read misses Dr : D cache reads (ie. memory reads) D1mr: D1 cache read misses D2mr: L2 cache data read misses Dw : D cache writes (ie. memory writes) D1mw: D1 cache write misses D2mw: L2 cache data write misses -------------------------------------------------------------------------------- Profile data file 'callgrind.out.7741' (creator: callgrind-3.3.0-Debian) -------------------------------------------------------------------------------- I1 cache: 65536 B, 64 B, 2-way associative D1 cache: 65536 B, 64 B, 2-way associative L2 cache: 1048576 B, 64 B, 8-way associative Timerange: Basic block 0 - 3810525588 Trigger: Program termination Profiled target: /home/nschmid/gromos/svn/gromosXX/x86_64/bin/md @f md.inp (PID 7741, part 1) Events recorded: Ir Dr Dw I1mr D1mr D1mw I2mr D2mr D2mw Events shown: Ir Dr Dw I1mr D1mr D1mw I2mr D2mr D2mw Event sort order: Ir Dr Dw I1mr D1mr D1mw I2mr D2mr D2mw Thresholds: 99 0 0 0 0 0 0 0 0 Include dirs: User annotated: Auto-annotation: on -------------------------------------------------------------------------------- Ir Dr Dw I1mr D1mr D1mw I2mr D2mr D2mw -------------------------------------------------------------------------------- 59,696,737,439 20,453,280,842 8,831,821,039 601,767 200,229,353 40,511,699 168,776 31,171,410 20,088,219 PROGRAM TOTALS -------------------------------------------------------------------------------- -- Auto-annotated source: /home/nschmid/gromos/svn/gromosXX/x86_64/src/interaction/../../../src/interaction/nonbonded/interaction/solvent_innerloop.cc -------------------------------------------------------------------------------- Ir Dr Dw I1mr D1mr D1mw I2mr D2mr D2mw -- line 27 ---------------------------------------- . . . . . . . . . math::Vec * const force_i = &storage.force(i); . . . . . . . . . . . . . . . . . . math::Vec const * const pos_j = &conf.current().pos(j); . . . . . . . . . math::Vec * const force_j = &storage.force(j); . . . . . . . . . DEBUG(9, "i = " << i << " j = " << j); . . . . . . . . . . . . . . . . . . // first atom vs. first atom . . . . . . . . . periodicity.nearest_image(*pos_i, *pos_j, r); 168,138,585 168,138,585 0 0 1,421 0 0 2 . const double tx = r(0) - (*pos_i)(0) + (*pos_j)(0); 100,883,151 100,883,151 0 0 184 . . . . const double ty = r(1) - (*pos_i)(1) + (*pos_j)(1); 100,883,151 100,883,151 0 100 90 0 100 . . const double tz = r(2) - (*pos_i)(2) + (*pos_j)(2); . . . . . . . . . . . . . . . . . . double e_lj = 0.0, e_crf = 0.0; . . . . . . . . . . . . . . . . . . // loop over atom pairs 638,926,623 100,883,151 168,138,585 100 0 272 100 0 10 for(unsigned int param = 0, atom_i = 0; atom_i < num_solvent_atoms; ++atom_i) { 403,532,604 201,766,302 100,883,151 100 134,553 0 100 1,613 . const double xi = (*(pos_i+atom_i))(0) + tx; 302,649,453 100,883,151 100,883,151 0 123,708 0 0 1,448 . const double yi = (*(pos_i+atom_i))(1) + ty; 706,182,057 201,766,302 201,766,302 0 125,634 0 0 1,000 . const double zi = (*(pos_i+atom_i))(2) + tz; 1,614,130,416 302,649,453 403,532,604 0 0 4,949 0 0 746 for(unsigned int atom_j = 0; atom_j < num_solvent_atoms; ++atom_j, ++param) { . . . . . . . . . DEBUG(15, "\tatoms: i: " << atom_i << " j: " << atom_j); 907,948,359 907,948,359 0 0 6,537,647 0 0 171,849 . const double x = xi - (*(pos_j+atom_j))(0); 605,298,906 605,298,906 0 100 6,587,194 0 100 173,798 . const double y = yi - (*(pos_j+atom_j))(1); 605,298,906 605,298,906 0 0 6,629,181 0 0 171,226 . const double z = zi - (*(pos_j+atom_j))(2); . . . . . . . . . 2,421,195,624 0 0 100 0 0 100 . . const double r2 = x * x + y * y + z * z; . . . . . . . . . DEBUG(15, "\tr2: " << r2); . . . . . . . . . assert(r2 != 0.0); 605,298,906 302,649,453 0 0 3,836 0 0 184 . const double r2i = 1.0 / r2; 1,210,597,812 . . . . . . . . const double ri = sqrt(r2i); 907,948,359 0 0 100 0 0 100 . . const double dist6i = r2i * r2i * r2i; 1,513,247,265 605,298,906 0 100 3,413 0 100 584 . const double dist6i_c12 = pair_parameter[param].c12 * dist6i; . . . . . . . . . 1,815,896,718 302,649,453 0 0 16,331 0 0 1,348 . e_lj += (dist6i_c12 - pair_parameter[param].c6) * dist6i; 2,118,546,171 907,948,359 0 100 8,117 0 100 1,083 . e_crf += pair_parameter[param].q * (ri - m_crf_2cut3i * r2 - m_crf_cut); . . . . . . . . . . . . . . . . . . const double f = (dist6i_c12 + dist6i_c12 - pair_parameter[param].c6) * 6.0 * dist6i * r2i + 2,723,845,077 605,298,906 0 2 5,440 0 2 684 . pair_parameter[param].q * (ri * r2i + m_crf_cut3i); . . . . . . . . . 605,298,906 . . . . . . . . const double fx = f * x; 605,298,906 . . . . . . . . const double fy = f * y; 302,649,453 . . . . . . . . const double fz = f * z; . . . . . . . . . 1,210,597,812 605,298,906 302,649,453 0 262,421 0 0 3,871 . (*(force_i+atom_i))(0) += fx; 907,948,359 302,649,453 302,649,453 100 33,260,417 0 100 971,453 . (*(force_j+atom_j))(0) -= fx; 907,948,359 302,649,453 302,649,453 0 193,492 0 0 4,409 . (*(force_i+atom_i))(1) += fy; 907,948,359 302,649,453 302,649,453 0 10,177,012 0 0 365,584 . (*(force_j+atom_j))(1) -= fy; 907,948,359 302,649,453 302,649,453 0 186,904 0 0 2,404 . (*(force_i+atom_i))(2) += fz; 907,948,359 302,649,453 302,649,453 0 9,921,693 0 0 255,507 . (*(force_j+atom_j))(2) -= fz; . . . . . . . . . 1,210,597,812 302,649,453 302,649,453 200 0 0 200 . . storage.virial_tensor(0, 0) += x * fx; 1,210,597,812 302,649,453 302,649,453 . . . . . . storage.virial_tensor(0, 1) += x * fy; 907,948,359 302,649,453 302,649,453 0 5,199 0 0 685 . storage.virial_tensor(0, 2) += x * fz; 1,210,597,812 302,649,453 302,649,453 100 0 0 100 . . storage.virial_tensor(1, 0) += y * fx; 1,210,597,812 302,649,453 302,649,453 . . . . . . storage.virial_tensor(1, 1) += y * fy; 907,948,359 302,649,453 302,649,453 100 0 0 100 . . storage.virial_tensor(1, 2) += y * fz; 907,948,359 302,649,453 302,649,453 0 2,577 0 0 259 . storage.virial_tensor(2, 0) += z * fx; 907,948,359 302,649,453 302,649,453 100 0 0 100 . . storage.virial_tensor(2, 1) += z * fy; 907,948,359 302,649,453 302,649,453 0 1,503 0 0 293 . storage.virial_tensor(2, 2) += z * fz; . . . . . . . . . } . . . . . . . . . } 100,883,151 67,255,434 33,627,717 0 283,981 0 0 261 . storage.energies.lj_energy[egroup][egroup] += e_lj; 168,138,585 100,883,151 67,255,434 0 4,627 0 0 1,017 . storage.energies.crf_energy[egroup][egroup] += e_crf; . . . . . . . . . } . . . . . . . . . . . . . . . . . . . . . . . . . . . template<typename t_nonbonded_spec> . . . . . . . . . inline void . . . . . . . . . interaction::Nonbonded_Innerloop<t_nonbonded_spec>::spc_innerloop . . . . . . . . . ( . . . . . . . . . topology::Topology & topo, -- line 97 ----------------------------------------Hier sind die Resultate: