Aliasing-Optimierung?



  • @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×1n \times 1-Matrizen) x\mathbf{ x }, y\mathbf{ y } gilt nämlich:

    x,y=xTy=[x0x1x2xn][y0y1y3yn]=x0y0+x1y1+x2y2++xnyn\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 x\mathbf{ x } dann auch einfach und simpel:

    x2=xTx=x02+x12+x22++xn2\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×33\times3 und 4×44\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 😉



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



  • @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 😉



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



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



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



  • @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 xTx\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.



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



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



  • 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";
    }
    
    

Anmelden zum Antworten