loop optimieren
-
Ja dieses grid-based pairlisting (Einteilung des Raums in Zellen) mache ich natuerlich (das ist der N*cutoff^3 Algorithmus). Alle Paare die ausserhalb eines gegebenen Abstand (cutoff) liegen werden gar nicht erst berechnet (gehen gar nicht in diesen Loop). Aber selbst mit dieser Paarliste muss ich trotzdem noch Millionen von Paaren beruecksichtigen damit die Physik des Problems auch nur einigermassen stimmt.
Ihr koennt mir wirklich glauben: ich muss diesen einen Loop Millionen mal aufrufen. Es geht einfach nicht anders (zumindest haben 30 Jahren Forschung auf dem Gebiet keinen alternativen Algorithmus fuer die Berechnung von Shortrange-Kraeften hervorgebracht).LG
Neph
-
kannst du eventuell posten welche zeile wieviel zeit zieht?
bist du auf double angewiesen? reicht nicht eventuell float?wieviel schneller muss es am ende sein?
-
also die einzelnen arrays die du da verwendest sehen für mich nicht wirklich cache-freundlich aus.
pack mal alle daten zu einem element (atom?) in eine struktur, und verwende dann ein array dieses struktur-typs.dann könntest du noch versuchen SSE zu verwenden. einfach nur beim compiler SSE zu aktivieren wird aber nix bringen, da die SSE vectorizer in C++ compilern nicht gerade mächtig sind. also SSE code mit hand schreiben. MSVC bietet einen haufen intrinsic funktionen an die man dazu verwenden kann.
falls "param" konstant ist kannst du den lookup pair_parameter[param].cXX noch aus der schleife rausziehen. kann durchaus was bringen, vor allem wenn der compiler nicht checkt dass es kein aliasing zwischen "pair_parameter" und anderen zeigerd (über die daten geändert werden) gibt.
was gleich der nächste punkt wäre: es könnte was bringen dem compiler mitzuteilen dass bestimmte zeiger "aliasing-frei" sind (wenn sie das sind, weiss ich natürlich nicht). dazu gibts bei MSVC und GCC extensions (__declspec und sowas).
was deine "if-sorge" angeht: bei dem haufen code in "calc" tut dir das if nicht weh, das merkst du nichtmal.
-
hustbaer schrieb:
also SSE code mit hand schreiben.
hat er schon versucht, siehe assembler forum
-
Das war aber ein anderer Loop
-
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!
