*
@Finnegan sagte in Matrix Template Klasse als Beispiel für die Verwendung von Allokatoren:
Das mit dem "nicht kopieren" ist auch eine Motivation, eine transponierte Matrix einfach als eine der anderen Elementorndung zu betrachten. Transponieren ist dann nämlich ebenfalls erstmal eine No-Op, die Element-Daten werden einfach in eine Matrix der jeweils anderen Ordnung gemoved, ohne sie (obendrein auch noch cache-ungünstig) kopieren zu müssen.
Ehrlich gesagt hatte ich das Thema noch nicht vertiefend betrachtet, da man in der Realität das nicht so häufig braucht, und ganz andere Optimierungsstrategien beim HPC eine Rolle spielen. Die trivial Operationen wie +, -, und Hadamard-Multiplikation kann man zwar für diese Fälle optimieren. Allerdings sind das Operationen, die im großen ganzen nur vom Speicherdurchsatz abhängen, und man nutzt in der Regel die gleichen Matrizen. Das betrifft sowohl die Speicherorganisation wie auch die eigentlichen Variablen in Programmcode.
Dazu braucht man allerdings auch Matrix-Rechenoperationen, die mit gemischten Odnungen umgehen können. Das eigentliche "Arbeit" des Transponierens würde dann z.B. in einer folgenden Matrix-Addition stattfinden, die bei gemischten Ordnungen eine der Matrizen anders durchlaufen müsste. Bei der Matrix-Multiplikation könnten entgegengesetzte Ordnungen sogar förderlich sein, da man beide Matrizen in der Speicher-Reihenfolge durchlaufen könnte, anstatt eine Zeilenweise und die andere Spaltenweise, was - zumindest bei naiver Implementierung - nicht so günstig für den CPU-Cache ist.
Der triviale Ansatz ist in der Realität eher unwichtig, da der Geschwindigkeitsunterschied so groß ist, dass man mit Sicherheit in real existierendem Code eine der linearen Algebra Bibliotheken nutzen wird. Als Standard hat sich dabei die BLAS und LAPACK Bibliothek etabliert. Es gibt eine Implementation GotoBLAS, die auf den Autoren Kazushige Goto zurückgeht. Er hat einen Artikel in der ACM TOMS veröffentlicht, in dem er darstellt wie man eine simple Matrizenmultiplikation auf einem modernen Prozessor optimieren kann. TLDR neben der notwendigen drei Schleifen für den eigentlichen Algorithmus muss man je Cacheebene drei weitere Schleifen in den Code einbauen. Die Erfahrung zeigt allerdings, dass die für den L1 Cache der Compiler das selbst ganz gut hinbekommt, bzw. man besser von Hand SIMD-Code dafür schreibt. Aber die anderen Schleifen muss man umsetzen, und so die Anzahl der Cache Misses minimieren.
Daher war mein Gedankengang, dass ich explizit nicht versuche das so allgemein wie möglich zu halten, sondern die Eigenschaften Speicherordnung und Transposition als Typen kodierte, um dann die Möglichkeit zu haben über eine Template-Spezialisierung direkt eine BLAS- bzw. LAPACK-Implementation zu nutzen. Da kommt man dann bei einem Aufruf ohne großen Aufwand aus, und übergibt direkt den Zeiger auf den Speicher an die entsprechende Routine.
Ich habe das mit den unterschiedlichen Ordnungen mal sehr generisch gelöst, indem ich den Matrizen row_stride und col_stride Integer-Attribute gegeben habe. Das sind die Anzahl der Elemente, die man überspringen muss, um in die jeweils nächste Zeile oder Spalte zu kommen. Für eine m×nm \times nm×n Matrix ist bei Row-Major-Ordnung row_stride =nnn und col_stride = 1 und bei Colum-Major-Anordung row_stride = 1 und col_stride =mmm. Damit wird die Index-Berechnung zu col_stride * j + row_stride * i, egal in welcher Ordung die Matrix vorliegt - so kann man an einigen Stellen auch auf einen Branch zur Fallunterscheidung zwischen Row/Column Major verzichten. Auch eventuelles Padding, z.B. so dass Zeilen/Spalten im Speicher immer ein vielfaches der Breite von SIMD-Datentypen sind, lässt sich damit gut handhaben.
Der Ansatz ist interessant, wobei für das SIMD-Speicherlayout bereits die Leading Dimension ausreichen ist. Ich werde das im Hinterkopf behalten und dann schauen was man daraus machen kann.
P.S. Danke fürs Feedback
