3. loop optimieren

loop optimieren

  • M

    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?

    0
  • N

    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

    0
  • R

    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.
    )

    0
  • N

    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

    0
  • R

    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?

    0

  • 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.

    0
  • R

    hustbaer schrieb:

    also SSE code mit hand schreiben.

    hat er schon versucht, siehe assembler forum 😉

    0
  • N

    Das war aber ein anderer Loop 🙂

    0
  • N

    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;
    }
  }
    0
  • ?

    Warum soll sich das so schlecht auf mehrere Cores aufteilen lassen?

    0
  • ?

    Bezahl lieber Leute die das können 😉

    0
  • N

    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. 🙄

    0
  • ?

    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. 😉

    0
  • ?

    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.

    0
  • N

    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

    0
  • ?

    Viele Cores verwendest du denn jetzt? Der Code da oben ist wohl kaum so ineffizient, das du da 50% rausholen kannst.

    0
  • ?

    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.

    0
  • N

    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

    0
  • S

    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?

    0
  • ?

    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

    0
Anmelden zum Antworten