[Programmiererstammtisch] "Zum ächzenden Compiler"

[alpha civ]

Der error-Vektor könnte mit .clear oder .fill ohne Schleife gesetzt werden, aber die Stelle bringt nicht so viel, da es ja eher an der doppelten for-Schleife hängt.

Die doppelte Schleife ist die Hauptarbeit: O(n^2) Laufzeit. Das bekomme ich aber nicht weg.
Ich habe mir schon überlegt, die innere Schleife zu parallelisieren, aber der zusätzliche Aufwand bezüglich der Thread-Syncronisation macht es vermutlich zu einem Null-Summen-Spiel.

Ich bin mir unsicher, wie C++ mehrdimensionale Arrays auflöst. Aber die [i][j] Zugriffe erfordern doch immer eine Multiplikation, die du doch auflösen könntest?!

Ich verwende doch Vektoren von Vektoren? Ansonsten wüsste ich hier nicht, wie ich was auflösen sollte.

Falls die Eingabe von round positiv ist, könnte man das auch durch (int)(val+0.5) ersetzen.

Davon kann ich nicht ausgehen. Die Mu-Matrix kann sowohl positiv als auch negative Einträge haben. Ich weis nur, dass die Einträge betragsmäßig kleiner oder gleich 0.5 ist, außer auf der Hauptdiagonalen, die besteht nur aus Nuller. Mu ist eine untere Dreiecksmatrix (und stammt aus der Gram-Schmidt-Orthogonalisierung).

[Harleen]

Ich habe mir schon überlegt, die innere Schleife zu parallelisieren, aber der zusätzliche Aufwand bezüglich der Thread-Syncronisation macht es vermutlich zu einem Null-Summen-Spiel.

Das kann man ja schnell ausprobieren. Einfach mal "#pragma omp parallel for" vor die Schleife schreiben und OpenMP in den Compilereinstellungen aktivieren- Fertig!

[alpha civ]

Matrix muss aber irgendwie anders definiert sein oder du hast noch eigene Operatoren definiert. Jedenfalls sehen einige Anweisungen wie mu[i] -= d; merkwürdig aus.

Mu ist die Matrix, mu ist ein Vektor.

[Harleen]

Für die Variablenbenamsung:


Falls die Iteratoren zu unübersichtlich sind, kannst du zumindest den zweifachen []-Operator aus der inneren Schleife killen:
Edit: Oder wird zuerst die rechte Seite [j] ausgewertet?

            for (int j=0; j < i; ++j) 
            {
                const Vector &MuRow = Mu[i];
                y = - d * MuRow[j] + err[j];
                t = mu[j] + y;
                err[j] = (t - mu[j]) - y;
                mu[j] = t;
            } 

Außerdem ist es einen Versuch wert, y und t innerhalb dieser Schleife als const zu deklarieren statt außerhalb. Mit eingeschalteter Compileroptimierung fallen damit hoffentlich die Kopien weg (d.h. y und t werden wegoptimiert). Vielleicht macht der Compiler das aber auch schon. Aber mit der jetzigen Struktur ist es für ihn schwerer zu erkennen, dass die Variablen wegoptimiert werden können.

[[VK]]

Der einzige Tipp den ich bei Multiplen Aufrufen geben kann ist es Cacheing auszunutzen oder den 2D Array in ein 1D Array zu schreiben und dann beim Aufruf die Position durch Multiplikation zu bestimmen...

[alpha civ]

y und t sollen nicht wegotimiert werden! Das dient ja gerade dazu, den durch Gleitkommaarithmetik bedingten Rundungsfehler zu berechnen, um diesen dann auszugleichen.
(Das ist der Kahan-Algortihmus.)

Wie kann MuRow als const deklariert werden, wenn es pro Schleifendurchlauf einen anderen Wert hat (oder haben Referenzen hier eine andere Bedeutung?).

[alpha civ]

Der einzige Tipp den ich bei Multiplen Aufrufen geben kann ist es Cacheing auszunutzen oder den 2D Array in ein 1D Array zu schreiben und dann beim Aufruf die Position durch Multiplikation zu bestimmen...

Ach das meint ihr mit Multiplikation. Die Linearisierung der Matrix. Aber macht das tatsächlich so einen Unterschied?

Cacheing ist hier nicht möglich. Der Algorthmus erzeugt sehr viele verschiedene Vektoren (mehrere Millionen), und die Länge eines Vektors ist von der Größenordung 300.

[Peregrin_Tooc]

Vielleicht kann man von der mathematischen Seite noch ran. Algorithmisch sehe ich nichts.

[alpha civ]

Werde das mal ausprobieren. Kannte ich bisher gar nicht.

[Harleen]

Wie kann MuRow als const deklariert werden, wenn es pro Schleifendurchlauf einen anderen Wert hat (oder haben Referenzen hier eine andere Bedeutung?).

MuRow muss natürlich vor der Schleife deklariert werden, um die Anzahl der Zugriffe zu verringern. So wie ich das hingeschrieben habe, bringt das gar nichts.
Als const kann es deklariert werden, weil diese Referenz bei jedem Schleifendurchlauf neu angelegt wird. Es ist ja eine ganz normale Deklaration. Durch die Referenz erspart man sich die Kopie.

[Zulan]

Folgende Funktion ist die Bremse meines Programms, bzw. hier wird die Hauptarbeit erledigt. Die Frage ist, was kann ich noch optimieren:

(Vector und Matrix ist einfach die std::vector<double> bzw. std::vector<std::vector<double>>.)

void nearestPlane(Vector& mu, const Matrix& Mu, int k)
{
    double d, y, t;

    int n = mu.size();

    Vector temp(k);
    Vector err(n);  // Error vector
    for (int i=0; i<n; ++i) {err[i] = 0;}
    for (int i=0; i<k; ++i) {temp[i] = mu[n-k+i];}
    
    for (int i = n-1; i >= 0; --i)
    {
        d = round(mu[i]);        
        if (d != 0)
        {
            mu[i] -= d;                        // |mu[i]| <= 0.5
            for (int j=0; j < i; ++j) 
            {
                y = - d * Mu[i][j] + err[j];
                t = mu[j] + y;
                err[j] = (t - mu[j]) - y;
                mu[j] = t;
            }
        }
        if (i >= n-k) {mu[i] = temp[i-n+k] - mu[i];}
        else {mu[i] *= -1;}
    }
    
} 

Interessante Sache, erstmal die Kleinigkeiten:

double y/t als const nach innen ziehen ist gut, weil besser lesbar, sollte aber fuer den Compiler keinen Unterschied machen. Fuer die Parallelisierung ist es aber wichtig.

err und temp solltest du mit dem fill / range constructor initialisieren - hilft vorallem der Lesbarkeit.

Mit der Variablenbenennung Mu/mu fang ich mal nicht an

Womit verbringt der Code hier wohl seine Zeit?
Da du mit n*n durch ne n*n Matrix gehst haengst du hauptsaechlich am Speicher.

Wie gross ist n? Wie gross deine Caches? --> Passen mu/tmp/err in den L1/L2/L3? Sollte bei n=300 gegeben sein -> ansonsten Blocking.

Hilft es die Matrix zu linearisieren? Vielleicht, aber vermutlich nicht viel da es das eigentliche Speicherproblem nicht angeht. Kompakte Speicherung ist wichtig. Es ist ja eine Dreiecksmatrix - probier mal boost:triangular_matrix.

Du hast unschoene Datenabhaenigkeiten zwischen den Iterationen: In jeder i-Iteration wird das err und mu der vorhergehenden Iteration gelesen und fuer die naechste geschrieben. In der inneren Schleife hast du die Abhaenigkeit ueber die Summe.

Das macht die Parallelisierung schwierig. Es gibt in OpenMP reduction Schleifen, glaube aber auch nicht, dass die das bei durchschnittlich ~150 Iterationen auszahlen.

Aber ich muss auch gestehen, ich seh nicht so richtig durch was der code letztlich macht. Sieht aus wie eine Matrix-Vektor Multiplikation, aber nicht wirklich. Idealer Weise findest du was in einer guten BLAS Bibliothek dafuer - die ist dann schon bis ins letzte Handoptimiert. MKL / ACML fuer C sind etabliert. Fuer C++ kannst du dir mal http://www.simunova.com/mtl4 anschauen.

Es lohnt aber auch ein Schitt zurueck - wie performt dein Code so wie er ist? In welchem Verhaeltnis steht das zu der theoretischen Maximalperformance (unter der Annahme, das Mu nicht aus einem Cache kommen kann - limitiert durch die Speicherbandbreite). Gibt es da noch viel rauszuholen?

[Ramkhamhaeng]

Ich verwende doch Vektoren von Vektoren? Ansonsten wüsste ich hier nicht, wie ich was auflösen sollte.

Selbst wenn die Vektoren im Speicher verstreut werden, müsste doch der Abstand zwischen aufeinander folgende Matrixelemente in beide Richtungen konstant sein!?
D.h. in Pointerschreibweise kann man sich so ohne Multiplikationen und ohne Indexvariablen durchhangeln (Pseudocode)

double *pMu_ij = &Mu[i,0];
const double d = &Mu[i,1]-&Mu[i,0]; //bin mir gerade unsicher, ob das hier 1 ist....
const double pMu_ijEnd = &Mu[i,i];
while( pMu_ij < pMu_ijEnd ){
[...]
pMu_ij += d
}

Davon kann ich nicht ausgehen. Die Mu-Matrix kann sowohl positiv als auch negative Einträge haben. Ich weis nur, dass die Einträge betragsmäßig kleiner oder gleich 0.5 ist, außer auf der Hauptdiagonalen, die besteht nur aus Nuller. Mu ist eine untere Dreiecksmatrix (und stammt aus der Gram-Schmidt-Orthogonalisierung).

Dann sind ja eigentlich nur die Werte -0.5 und 0.5 (plus Rundungsfehler weg von Null) relevant. D.h. der Aufruf der Rundungsfunktion kann durch eine Bitmaske ersetzt werden, welche den Exponenten überprüft. Der Exponent muss ja 2^-1 sein, damit zu +1 oder -1 gerundet wird. Leider nervt dann noch das Vorzeichen.

Achtung Spoiler:

#include "stdio.h"
int main(){

  double testvariablen[6] = { 
    0.500001, 0.5, 0.49999999,
    -0.500001, -0.5, -0.49999999
  };  

  long long mask1 = 0x3FE0000000000000;//0.5
  long long mask2 = 0xBFE0000000000000;//-0.5

  size_t i;
  for( i=0; i<6U; ++i ){
    const unsigned long long dBits = *((unsigned long long*)&testvariablen[i] );
    switch( (dBits&mask2)>>52 ){
      case 0x3FE: printf("+1 ");
                  break;
      case 0xBFE: printf("-1 ");
                  break;
      default: printf(" 0 ");
    }   
  }

  return 0;
}

Ist aber mehr akademischer Natur, denn eigentlich müssten heutige CPUs fürs Runden einen eigenen Maschinenbefehl haben, der schneller ist. Da kostet das Switch dann mehr als man spart (ungetestet)

Bei der Verarbeitung von d fällt noch auf, dass du die zwei Fälle d=1 und d=-1 in zwei If-Zweigen abhandeln könntest. So sind sogar alle Multiplikationen komplett eliminierbar. (Was kein Optimierer hinbekommt, da die Werte für d zur Kompilezeit nicht bekannt sind.)


Edit: Der Verweis auf BLAS ist nat. besser und sicher nützlicher als mein Zeug Dreiecksmatrix ist auch gut.
Irgendwo kommt mir der Codeschnipsel in seiner Struktur auch bekannt vor. Kann sein dass es da wirklich schon einen handoptimierten Code gibt. Dann besteht aber immer noch die Schwierigkeit die Eingabedaten in die passende Form für diese Blas-Methode zu bekommen.

[alpha civ]

Interessante Sache, erstmal die Kleinigkeiten:

double y/t als const nach innen ziehen ist gut, weil besser lesbar, sollte aber fuer den Compiler keinen Unterschied machen. Fuer die Parallelisierung ist es aber wichtig.

Warum const? Hier gibt es doch pro Schleifendurchlauf eine neue Zuordnung. Wichtig ist, dass in der inneren Schleife y und t nicht wegoptimiert werden, weil sonst der Rundungsfehler nicht berechnet wird.

err und temp solltest du mit dem fill / range constructor initialisieren - hilft vorallem der Lesbarkeit.

Stimmt.

Mit der Variablenbenennung Mu/mu fang ich mal nicht an

Mathenotation. Eigentlich ja [math]\mu_{i,j}[\math] statt Mu, aber da ich schon mu als Orthogonalvektor verwende...

Womit verbringt der Code hier wohl seine Zeit?
Da du mit n*n durch ne n*n Matrix gehst haengst du hauptsaechlich am Speicher.

Wie gross ist n? Wie gross deine Caches? --> Passen mu/tmp/err in den L1/L2/L3? Sollte bei n=300 gegeben sein -> ansonsten Blocking.

n ist zwischen 150 und 400, das hängt von meiner Eingabe ab.
Was ist L1/L2/L3 ?

Hilft es die Matrix zu linearisieren? Vielleicht, aber vermutlich nicht viel da es das eigentliche Speicherproblem nicht angeht. Kompakte Speicherung ist wichtig. Es ist ja eine Dreiecksmatrix - probier mal boost:triangular_matrix.

Auch das verwerde ich mal tun. Aber nicht Boost verwenden, deren Matrizen sind langsam.

Du hast unschoene Datenabhaenigkeiten zwischen den Iterationen: In jeder i-Iteration wird das err und mu der vorhergehenden Iteration gelesen und fuer die naechste geschrieben. In der inneren Schleife hast du die Abhaenigkeit ueber die Summe.

Das macht die Parallelisierung schwierig. Es gibt in OpenMP reduction Schleifen, glaube aber auch nicht, dass die das bei durchschnittlich ~150 Iterationen auszahlen.

Ich befürchte auch, das da wenig bei rauskommt. Bei meinen ersten Test gestern abern noch hatte ich eher das Gefühl, das es langsamer geworden ist. Kann mich aber auch ihren. Muss mich noch mehr mit OpenMP beschäftigen.

Aber ich muss auch gestehen, ich seh nicht so richtig durch was der code letztlich macht. Sieht aus wie eine Matrix-Vektor Multiplikation, aber nicht wirklich. Idealer Weise findest du was in einer guten BLAS Bibliothek dafuer - die ist dann schon bis ins letzte Handoptimiert. MKL / ACML fuer C sind etabliert. Fuer C++ kannst du dir mal http://www.simunova.com/mtl4 anschauen.

Das ist die Nearest Plane Methode. Gehört zur Gitter-Theorie. Damit kann man zu einem gegebenen reellwertigen Vektor einen nahen Gittervektor approximieren.
Ich verwende das, um sehr viele verschiedene kurze Gittervektoren zu erzeugen.

Es lohnt aber auch ein Schitt zurueck - wie performt dein Code so wie er ist? In welchem Verhaeltnis steht das zu der theoretischen Maximalperformance (unter der Annahme, das Mu nicht aus einem Cache kommen kann - limitiert durch die Speicherbandbreite). Gibt es da noch viel rauszuholen?

Das muss ich endlich in der Tat mal testen, wo da überhaupt noch Optimierungsmöglichkeiten bestehen.

[alpha civ]

Dann sind ja eigentlich nur die Werte -0.5 und 0.5 (plus Rundungsfehler weg von Null) relevant. D.h. der Aufruf der Rundungsfunktion kann durch eine Bitmaske ersetzt werden, welche den Exponenten überprüft. Der Exponent muss ja 2^-1 sein, damit zu +1 oder -1 gerundet wird. Leider nervt dann noch das Vorzeichen.

d kann auch andere ganzzahlige Werte annehmen, je nachdem, welchen Wert mu[i] hat.

Edit: Der Verweis auf BLAS ist nat. besser und sicher nützlicher als mein Zeug Dreiecksmatrix ist auch gut.
Irgendwo kommt mir der Codeschnipsel in seiner Struktur auch bekannt vor. Kann sein dass es da wirklich schon einen handoptimierten Code gibt. Dann besteht aber immer noch die Schwierigkeit die Eingabedaten in die passende Form für diese Blas-Methode zu bekommen.

Ich denke nicht, das ich BLAS brauchen werde. Ich brauche ja eigentlich keine richtige Matrix (im mathematischen Sinne), sondern nur einen Container mit Matrix-Struktur. Mathe-Operationen brauche ich nicht. Ein simples Array will ich aus den üblichen Gründen aber auch nicht nehmen.

[Ramkhamhaeng]

L1, L2 und L3 sind die Cachestufen der CPU. Ganz grob gesagt hat L1 die geringsten Zugriffszeiten, hat aber auch den wenigsten Platz.