loop optimieren

Nephelauxetic

Hallo

Ich habe etwas ähnliches wie diesen Code: In zwei verschachtelten for-loops (num ist 3-6, und wird aus einer Datei gelesen) führe ich Berechnungen durch.

void calc_pair(const unsigend int i, const unsigned int j, const unsigned int num) {
  for(unsigned int x = 0; x < num; ++x) {
    for(unsigned int y = 0; y < num; ++y) {
      calc(i,j,x,y);
    }
  }  
}

Die Berechnung lässt sich leider nicht wirklich mit SSE2 parallelisieren. Deshalb möchte ich paarweise über x loopen damit der Compiler dann den Code einfach mit SSE2 umsetzten kann.

void calc_pair(const unsigend int i, const unsigned int j, const unsigned int num) {
  for(unsigned int x = 0; x < num; x+=2) {
    for(unsigned int y = 0; y < num; ++y) {
      calc(i,j,x,y);
      calc(i,j,x+1,y);
    }
  }  
}

Nun meine Frage. num ist leider nicht immer gerade. Hat jemand eine Idee wie man das für ungerade num umsetzt?
die Funktion calc_pair wird Millionen mal aufgerufen und deshalb denke ich dass ein if nicht wirklich effizient ist. Was sind die Meinungen dazu? Ideen wie man das if vermeiden kann?

LG
Neph

Bashar

Duff's Device

Am Rande: Funktioniert das so mit SSE2, wie du dir das denkst? calc ist doch ein Funktionsaufruf und nicht irgendeine kleine, parallelisierbare Operation.

Nephelauxetic

Das ist nur verkürzt in pseudocode geschrieben. In Wirklichkeit ist calc eine Reihe von arithmetischen Operationen mit u.a. einer Division und einer Wurzel.

Shade Of Mine

mit openMP einfach mal die schleife parallelisieren?

Nephelauxetic

Das haben wir längst getan (OpenMP und MPI). Dies tun wir aber auf höhreren Ebenen (wie gesagt, der Loop wird Millionen mal aufgerufen). 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. Da komm ich doch nicht um SSE rum, oder?

LG
Neph

Walli

Hast Du mal gecheckt was so an Kommunikation zwischen den Nodes abläuft?

Jester

Walli schrieb:

Hast Du mal gecheckt was so an Kommunikation zwischen den Nodes abläuft?

Das halte ich auch für eine interessante Frage. Wenn es nicht skaliert, dann skaliert es wohl auch nicht, wenn Du es lokal schneller machst.

Nephelauxetic

Ja sicher. Das Problem ist in der jetzigen Form nicht sehr kommunikationslastig. Das Problem ist mehr dass Amdahl's Law zuschlägt und das Parallelisieren weiterer Komponenten des Programms hoch komplex ist.
Aber lasst uns doch bitte beim Thema bleiben und dass ist gute Performance auf einem Core Wenn man die man hat kann man immer noch auf mehrere ausweichen. Das hat auch eine ökologische Komponente: ist doch schwachsinnig auf 2 Nodes laufen zu lassen wenn man die selbe Performance auf 1 haben könnte.

Nephelauxetic

Jester schrieb:

Wenn es nicht skaliert, dann skaliert es wohl auch nicht, wenn Du es lokal schneller machst.

Nein logisch nicht aber wenn das Programm auf einem Core schnell läuft brauch ich gar nicht erst auf 16 Prozessoren...

Erst mal würde ich schauen ob man außerhalb was machen kann, damit es nicht mehr so oft aufgerufen wird.

Wie groß ist num eigentlich so im durchschnitt?

Nephelauxetic

num ist klein (Anzahl Atome in einem Lösungsmittelmolekül). In der Regel 3 kann aber auch bis ca. 6 werden.

Erst mal würde ich schauen ob man außerhalb was machen kann, damit es nicht mehr so oft aufgerufen wird.

Das ist bei physikalischen Berechnungen nicht so einfach. Wenn du 1'000 Teilchen hast zwischen denen du die Kräfte ausrechnen musst so das sind das halt schon 1'000'000 mögliche Paare.

void calc_pair(const unsigend int i, const unsigned int j, const unsigned int num) {
  unsigned int x;
  for(x = 0; x < num; x+=2) {
    for(unsigned int y = 0; y < num; ++y) {
      calc(i,j,x,y);
      calc(i,j,x+1,y);
    }
  }
  if (x != num) { // num ungerade
    for(unsigned int y = 0; y < num; ++y) {
      calc(i,j,num-1,y);
    }
  }
}

Das wäre eine mögliche Lösung bei der ich aber ein if verwenden muss. Zudem funktioniert sie für den Fall num == 1 nicht. Das ist aber egal: ich kann diesen sehr unwahrscheinlichen Fall einfach verbieten (da ist diese Art des Loops sowieso nicht die effizienteste).
Um das if mache ich mir aber Sorgen. Ich habe Angst dass es die SSE performance gewinne aufressen wird... Deshalb wäre irgendeine Lösung ohne if besser - vielleicht durch anderes loopen aber mir fällt irgendwie nichts gescheites ein.

LG
Neph

Shade Of Mine

loop unrolling sollte dein compiler koennen.
wenn nicht, dann wurde duffs device als loesung fuer haendisches unrolling schon genannt.

das ganze geht natuerlich auch mit templates um den loop unzurollen. aber ich denke nicht dass die schleife selber ein wirkliches optimierungspotential hat.

du kannst aber mal checken ob dein compiler die funktion inlined.
danach ueberpruefst du was genau der bottleneck ist. ist es die rohe cpu leistung? hast du zuviele cache misses? sind die branchpredictions zu oft falsch? oder kann einfach nicht schnell genug daten nachgeschoben werden...

der code so wie er ist, bietet er kein wirkliches potential. du kannst nur versuchen das zusammenspiel zwischen der schleife und den rest der anwendung zu optimieren.

uU kannst du aggressive caches einsetzen. kann man halt aktuell nur raten.

solche mikrooptimierungen wie du hier willst, bringen idR aber nichts. beachte die 80/20 Regel: 80% der zeit wird in 20% des codes verbracht. wenn die schleife wirklich zuviel zeit frisst, ist die richtige alternative sie seltener aufzurufen. da helfen uU gute caches oder vorkalkulierte tabellen.

maximAL

Nephelauxetic schrieb:

num ist klein (Anzahl Atome in einem Lösungsmittelmolekül). In der Regel 3 kann aber auch bis ca. 6 werden.

wenns garantiert nicht größer als 6 werden kann, könnte man doch ein switch von 1 - 6 nehmen und für jeden fall ausprogrammieren?

Nephelauxetic

maximAL schrieb:

wenns garantiert nicht größer als 6 werden kann, könnte man doch ein switch von 1 - 6 nehmen und für jeden fall ausprogrammieren?

dies hab ich für den wichtigsten Fall 3 bereits getan. Bei Wasser (ein 3er Fall) kann man dann noch weitere Vereinfachungen machen und spart weitere Performance.

Shade Of Mine schrieb:

loop unrolling sollte dein compiler koennen.
wenn nicht, dann wurde duffs device als loesung fuer haendisches unrolling schon genannt.

Der compiler macht unrolling. Ich optimiere auch mit dem profile. Das heisst aber noch nicht dass der Compiler kapiert dass er alles doppelt in der Schleife hat und deshalb alles mit SSE rechnen kann.

Shade Of Mine schrieb:

das ganze geht natuerlich auch mit templates um den loop unzurollen. aber ich denke nicht dass die schleife selber ein wirkliches optimierungspotential hat.

Nein die Schleife sicher nicht aber den inhalt der calc "Funktion" die ja ca. 30 Statements ist die man umordnen kann damit der compiler das mit SSE umsetzt. Ich denke im Moment daran dass ich Templates benutzen könnte je nachdem ob num gerade oder ungerade ist (funktion verdoppeln) und dann das if vor den loop über i und j schieben. So hab ich es wenigstens nicht im loop - wenn auch dann redundanten code.

Shade Of Mine schrieb:

du kannst aber mal checken ob dein compiler die funktion inlined.
danach ueberpruefst du was genau der bottleneck ist. ist es die rohe cpu leistung? hast du zuviele cache misses? sind die branchpredictions zu oft falsch? oder kann einfach nicht schnell genug daten nachgeschoben werden...

wie gesagt: es ist gar keine Funktion. Ich muss das mal durch cachegrind laufen lassen um zu sehen wie das mit dem Cache aussieht. Ich weiss es nicht genau. Dies zu optimieren wär aber nicht trivial (siehe unten wegen Paarliste).

Shade Of Mine schrieb:

der code so wie er ist, bietet er kein wirkliches potential. du kannst nur versuchen das zusammenspiel zwischen der schleife und den rest der anwendung zu optimieren.

Der Inhalt der calc "Funktion" sieht wie folgt aus.

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;
      assert(r2 != 0.0);
      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;

Das Problem ist das die ersten Zeilen bis ca double fx nicht wirklich SSE kompatibel sind. Würde ich aber Paarweise loopen könnte ich die Wurzel und die inverse doppelt rechnen etc.

Shade Of Mine schrieb:

solche mikrooptimierungen wie du hier willst, bringen idR aber nichts. beachte die 80/20 Regel: 80% der zeit wird in 20% des codes verbracht. wenn die schleife wirklich zuviel zeit frisst, ist die richtige alternative sie seltener aufzurufen. da helfen uU gute caches oder vorkalkulierte tabellen.

Die Schleife frisst 80% der Zeit (meint mein Profiler). Es gäbe eine möglichkeit den Loop anders zu machen (andere i, j) Paare aber diese cache optimierte Methode ist ein N^2 loop und das ist zum finden der Paare viel zu ineffizient. Dann hab ich zwar einen cache optimieren loop aber die Berechnung der Paarliste skaliert falsch.
Ich muss aber gestehen dass ich deinen Post nur ansatzweise verstehe. Ich werde mich da ein bisschen einlesen.

LG
Neph

rapso

ich glaube du optimierst am falschen ende, du solltest dafuer sorgen dass es nicht N^2 ist, statt den loop zu optimieren.

1.wuerdest du erstmal grob den raum einteilen, duerfte 99% der berechnungen wegfallen

(
2. wenn du bei "const double r2i = 1.0 / r2;" einfach pruefst ob r2 ueber einem gewissen wert ist und dann die berechnungen ueberspringen, wuerdest du mehr sparen als du mit SSE rausbekommen kannst. und da sich r2i durch die ganzen multiplikationen propagiert und alles gegen 0 laufen laesst, duerfte es kaum einen unterschied geben.
)

Nephelauxetic

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

rapso

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?

hustbaer

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.

rapso

hustbaer schrieb:

also SSE code mit hand schreiben.

hat er schon versucht, siehe assembler forum

Nephelauxetic

Das war aber ein anderer Loop