Aliasing-Optimierung?

Quiche Lorraine

@Steffo
Ich habe mir mal den Code in dem Link angeschaut.

Die erste Optimierung scheint wie schon gesagt der Austausch des Skalarproduktes mit der quadratischen Länge zu sein, sofern beide Vektoren identisch sind. So richtig überzeugt mich aber diese Optimierung nicht.

Die zweite Optimierung steckt in der Umformung der MItternachtsformel, bekannt auch als ABC Formel. Da b hier immer 2.0 * dot(oc, r.direction()) ist, lässt sich die Formel in eine vereinfachte Form bringen. Dadurch spart man 3 Multiplikationen. Das könnte man merken.

Steffo

@Quiche-Lorraine Danke, aber die Optimierung der Mitternachtsformel war mir schon klar. Die length_squared() Optimierung lässt sich zumindest in nicht optimierten Assembler-Code nachvollziehen. Beim optimierten Assembler-Code ist mir nichts aufgefallen, aber ich bin auch kein Assembler-Experte.

Schlangenmensch

Bei der ersten Optimierung hat man einen Aufruf von direction() weniger und damit eine vec3 Kopie weniger. Ob das im optimierten Code mit RVO und so einen Unterschied macht, weiß ich nicht, wahrscheinlich eher nicht. Hab ich mir aber nicht angeguckt

Steffo

@Schlangenmensch Danke, dann seht ihr beide das wohl wie ich. Kann das Tutorial übrigens nur empfehlen! Cool ist auch, dass überhaupt keine Libs benötigt werden!

Finnegan

@Steffo Vielleicht ist die optimierte Variante auch noch interessant, da sieht man nämlich, dass der Compiler genug Code-Verständnis hat, um für beide Varianten exakt den selben Code zu generieren:

https://godbolt.org/z/v1sY4TGT7

Meiner Meinung nach macht es aber vor allem für für simpleren Code Sinn, die Längenberechnung über das Skalarprodukt zu definieren. Ich würde eventuell sogar noch einen Schritt weiter gehen und bereits das Skalarprodukt über die Matrixmultiplikation definieren. Für die Spaltenvektoren (das sind de facto $n \times 1$ -Matrizen) $\mathbf{ x }$ , $\mathbf{ y }$ gilt nämlich:

$\langle \mathbf{ x }, \mathbf{ y } \rangle = \mathbf{ x }^T\mathbf{ y }=\begin{bmatrix} x_0 & x_1 & x_2 & \ldots & x_n\end{bmatrix} \cdot \begin{bmatrix} y_0 \cr y_1 \cr y_3 \cr \vdots \cr y_n \end{bmatrix}=x_0 y_0 + x_1 y_1 + x_2 y_2 + \ldots + x_n y_n$

Somit ist die Länge des Spaltenvektors $\mathbf{ x }$ dann auch einfach und simpel:

$\lvert \lvert \mathbf{ x } \rvert \rvert_2=\sqrt{ \mathbf{ x }^T\mathbf{ x } } = \sqrt{ x_0^2+x_1^2+x_2^2 + \ldots + x_n^2}$

Das ist aber sehr stark davon abhängig, wie flexibel deine Matrixmultiplikation ist und ob das Transponieren in diesem Fall etwas kostet (ähnlicher Code von mir würde da einfach die Matrix anders abschreiten, was letztendlich zu dem selben effizienten Maschinencode führt). Falls du diese nur für quadratische $3\times3$ und $4\times4$ -Matrizen implementieren willst, dann ist der zusätzliche Skalarprodukt-Code so oder so notwendig. Ansonsten wirds IMHO eleganter und übersichtlicher, wenn mans so wie ich beschrieben habe macht

Steffo

@Finnegan Matrixmultiplikation ist natürlich extrem generisch, aber ich habe ein wenig meine Zweifel, dass der Compiler dann genau so gut optimieren kann wie spezialisierter Code. Wenn er das kann, dann hat er wahrscheinlich mehr Aufwand investiert.

Finnegan

@Steffo sagte in Aliasing-Optimierung?:

@Finnegan Matrixmultiplikation ist natürlich extrem generisch, aber ich habe ein wenig meine Zweifel, dass der Compiler dann genau so gut optimieren kann wie spezialisierter Code. Wenn er das kann, dann hat er wahrscheinlich mehr Aufwand investiert.

Ja, ich habe meine Aussage mit Edits weiter eingeschränkt. Ich hab da eher eine Implementierung von mir im Kopf, die für die Matrixmultiplikation in diesem Fall den selben Code generieren würde, wie das handgeschriebene Skalarprodukt. Da bin ich etwas über die Stränge geschlagen das so direkt zu empfehlen - wenn man eine flexible und gut implementierte Matrixmultiplikation hat, dann geht das natürlich und vereinfacht die Codebasis.

Also: Bitte nur als Anregung unter Vorbehalt verstehen... mir gings auch vornehmlich um den Link zum compiler-optimierten Code von dir

Steffo

@Finnegan Danke, aber den Godbolt-Link findest du in meinem zweiten Beitrag.

Finnegan

@Steffo sagte in Aliasing-Optimierung?:

@Finnegan Danke, aber den Godbolt-Link findest du in meinem zweiten Beitrag.

Es ging mir um meinen Link oben in meinem (längeren) Beitrag. Da habe ich mir das mal mit Optimierungen angesehen und GCC generiert da für beide Varianten exakt den selben Code. Just FYI. Ich hatte den Eindruck dir war das nicht ganz bewusst:

https://godbolt.org/z/v1sY4TGT7

*john 0

@Finnegan sagte in Aliasing-Optimierung?:

Ansonsten wirds IMHO eleganter und übersichtlicher, wenn mans so wie ich beschrieben habe macht

Wenn man so weit in das Thema einsteigen will, ist es sinnvoller gleich eine optimierte BLAS Library zu nutzen.

Finnegan

@john-0 sagte in Aliasing-Optimierung?:

@Finnegan sagte in Aliasing-Optimierung?:

Ansonsten wirds IMHO eleganter und übersichtlicher, wenn mans so wie ich beschrieben habe macht

Wenn man so weit in das Thema einsteigen will, ist es sinnvoller gleich eine optimierte BLAS Library zu nutzen.

Ich denke das sind sehr unterscheidliche Anwendungsfälle. BLAS-Bibliotheken sind eher für sehr große, dynamische Matrizen ausgelegt. Für Computergrafik will man eher kleine "flache" (automatische/Stack-Objekte) Matrizen haben, die man auch mal einfach in ein zusammenhängendes Array für einen (Vertex-) Buffer oder ähnliches packen kann. Das muss nicht so ein schweres Geschütz werden wie die wissenschaftlichen Lineare-Algebra-Bibliotheken, aber ein bisschen Feintuning ist da nicht verkehrt.

Ich arbeite da z.B. mit Code, der mit Row-Major und Colum-Major-Matrizen gleichzeitig umgehen kann und bei der Multiplikation die Matrizen entsprechend anders abschreitet, je nachdem welche Speicherordnung diese haben. Transponieren ist da de facto eine No-Op, da die Matrix-Daten dann einfach als in der jeweils anderen Ordnung vorliegend interpretiert werden. Daher wird bei mir für $\mathbf x^T\mathbf x$ auch der selbe Code generiert, als hätte ich das Skalarprodukt manuell implementiert.

Der eigentliche Grund für die Unterstützung beider Element-Ordnungen war eigentlich Kompatibilität mit möglichst vielen Bibliotheken und APIs, die das nicht alle einheitlich handhaben. Mir ist aber schon klar, dass das man das nicht in jedem Hobby-Projekt so machen möchte - ich habe den Vorschlag etwas vorschnell und unter dem Eindruck meines eigenen Codes gemacht und beim Schreiben auch ehrlich gesagt nicht daran gedacht, dass ein naives Transponieren höhere CPU- und Speicherzugriffs-Kosten mit sich bringt.

Finnegan

@Steffo

https://godbolt.org/z/v1sY4TGT7

Ich muss allerdings noch anmerken, dass hier mit -O3 immer noch 3 Multiplikationen durchgeführt werden. Eigentlich können moderne CPUs gerade diese Multiplikationen in nur eine einzigen AVX-Instruktion durchführen. Die 256-Bit AVX-Register können schliesslich mit bis zu 4 doubles gleichzeitig rechnen. Ein zusätzlicher -mavx/-mavx2-Parameter hat allerdings keinen Unterschied gemacht, ~~erst ein sehr spezifisches Tuning für eine entsprechende CPU hat die Multiplikationen dann zu einer Instruktion zusammengefasst~~:

https://godbolt.org/z/Ma4xG5b97

Korrektur: fmadd ist natürlich ein "Fused Multiply-Add", also Multiplikation und Addition in einer Instruktion. Gut, so geht es auch, auch wenn ich es "per Hand" wahrscheinlich so gemacht hätte die Multiplikationen in einer Instruktion zu machen und dann die Elemente aufzuaddieren (gibts dafür nicht auch eine Instruktion? Summe über AVX-Vektor-Elemente?). Ich vertrau hier aber erstmal dem Compiler, dass er weiss, warum er das so macht - hier sind quasi die 2 extra Additionen in die 3 Multiplikationen "hineingewoben". Auch ne Strategie

znver1 steht für AMD Ryzen der ersten Generation, ich denke aber solcher Code dürfte sich auch noch mit einer Reihe anderer Tuning-Parameter erzeugen lassen. Das zeigt einmal mehr, dass man sich immer den Maschinencode anschauen sollte. Im Zweifelsfall erzeugt man sogar besser statt generischem Code, der überall läuft, solchen, der eine minimale CPU-Generation/Features voraussetzt (wie eben AVX/AVX2). Zumindest wenn man es möglichst effizient haben will.

Steffo

@Finnegan Interessant zu sehen, dass das so einen Unterschied macht. Danke.

*john 0

Die AVX Befehle sind zum Teil etwas sonderbar. Die folgende Abfolge kommt mit AVX und SSE aus, und liefert das Ergebnis.

#include <iostream>
#include <immintrin.h>

class vec3 {
public:
    __m256d e;

    vec3() : e (_mm256_setzero_pd ()) {}
    vec3(const double e0, const double e1, const double e2) : e (_mm256_set_pd (0.0, e2, e1, e0)) {}
    // ...
    double length_squared() const;
};

double vec3::length_squared () const {
    double ret[2] = {0.0, 0.0};

    __m256d t    = _mm256_mul_pd (e, e);
    __m256d r    = _mm256_hadd_pd (t, t);
    __m128d lo = _mm256_extractf128_pd (r, 0);
    __m128d hi = _mm256_extractf128_pd (t, 1);
    __m128d tt = _mm_blend_pd (lo, hi, 1);
    __m128d rr = _mm_hadd_pd(tt, tt);

    _mm_store_pd (ret, rr);

    return ret[0];
}

double dot(const vec3 &u, const vec3 &v) {
    return u.e[0] * v.e[0]
      + u.e[1] * v.e[1]
      + u.e[2] * v.e[2];
}

int main () {
    vec3 a(1.0, 1.0, 1.0), b(2.0, 2.0, 2.0);

    std::cout << dot (a, a) << "\n";
    std::cout << dot (b, b) << "\n";
    std::cout << b.length_squared() << "\n";
}