Klassentemplate fuer Matrizen, Spezialisierung fuer quadratische Matrizen
-
Ich habe ein Matrixtemplate, dass durch Typ und Dimension parametrisiert wird:
template<typename T, int r, int c> class Matrix { ..Sie implementiert beispielsweise die Matrixmultiplikation typsicher.
- Manche Operationen wie
operator*=(..)oderdiag()(Matrix mit Diagonalelementen zurueckgeben) machen nur Sinn im Zusammenhang fuer quadratische Matrizen. Dabei gibt es verschiedene Moeglichkeiten:
a)
Matrix<T,r,r>& operator*=(const Matrix<T,r,r>& m);b)
template<typename T, int r, int c> Matrix<T,r,c>& Matrix<T,r,c>::operator*=(const Matrix<T,r,c>& m) { static_assert(r == c, "Matrix::operator*=(..): matrix has to be quadratic");c)
template<T,n> class QuadMatrix : public Matrix<T,n,n> { ..Vererbung wollte ich nicht benutzen, da durch das Template sowieso immer ein neuer Typ angelegt wird. Zwischen a) und b) besteht prinzipiell nur ein Unterschied in der Fehlermeldung. Variante b) kann natuerlich weiter ausgebaut werden, indem die Parameter der zeiten Matrix ebenfalls frei gewaehlt werden um mittels
static_asserteine bessere Fehlermeldung auszugeben. Ich habe mich fuer Variante b) entschieden, da die Fehlermeldung problembezogener ist. Das ist aber nicht die eigentliche vorgehensweise bei Templatespezialisierung. Welche Variante ist im Allgemeinen zu bevorzugen? Gibt es noch weitere?- Mit der Klasse sollen natuerlich auch lineare Gleichungssysteme der Form Ax = b geloest werden und die (eindeutige) Loesung in Form eines Spaltenvektors (Matrix<T,r,1>) zurueckgegeben werden. Auch hier soll es nur fuer quadratische Matrizen moeglich sein. Es gibt die Moeglichkeit a) das als Methode
Matrix<T,d,1> solve(const Matrix<T,d,1>& v) constoder b) als freie Funktion
Matrix<T,d,1> solve(const Matrix<T,r,c>& m, Matrix<T,d,1>)zu implementieren. Auch hier kann mit den Moeglichkeiten aus 1.a) bzw. 1.b) die korrekten Matrixdimensionen sichergestellt werden.
Konkret: Soll
solveals Methode oder freie Funktion implementiert werden? Welche Moeglich ist zu bevorzugen, wenn die korrekten Matrixdimensionen schon beim Uebersetzen sicherzustellen sind?edit:
Matrix.h: http://codepad.org/IZICoNH5
Zero.h: http://codepad.org/WoYzH0oQ
main.cpp: http://codepad.org/98JVUefUUebersetzt mit g++ -std=c++0x main.cpp
- Manche Operationen wie
-
bei 1) würde ich Methode 2 wählen, das gibt aussagekräftigere Fehlermeldungen.
bei 2) würde ich die freien FUnktionen wählen. Dies ist zum einen Konvention, zum anderen weil ich alle 3 teile von Solve als gleichberechtigt erachte. typischerweise wird solve aber als Funktion mit 3 Argumenten implementiert - hauptsächlich um die Kopie zu sparen. Könnte man aber auch mit einem expression-template lösen, falls dein System das kann.
Was sind denn die Anwendungsgebiete der Bibliothek? Spontan würde ich behaupten, dass feste Größen eher kleine Matrizen bedeutet?
-
knivil schrieb:
template<T,n> class QuadMatrix : public Matrix<T,n,n>Im nächsten C++-Standard kannst du dafür
usingbenutzen (Template-Aliases), allerdings weiss ich nicht, wie gut der Compilersupport dafür bereits ist.knivil schrieb:
Welche Variante ist im Allgemeinen zu bevorzugen? Gibt es noch weitere?
Du könntest z.B. das Klassentemplate spezialisieren. Ob das besonders sinnvoll ist, hm. Jedenfalls hättest du Codeduplizierung (ein bisschen weniger, wenn du eine gemeinsame Basis verwendest).
template <typename T, int R, int C> class Matrix { ... }; template <typename T, int N> class Matrix<T, N, N> { ... };otze schrieb:
typischerweise wird solve aber als Funktion mit 3 Argumenten implementiert - hauptsächlich um die Kopie zu sparen. Könnte man aber auch mit einem expression-template lösen, falls dein System das kann.
Am naheliegendsten wäre dafür wohl Move-Semantik, wenn knivil schon C++0x benutzt.
-
Also: Ich spezialiesiere mit
static_assertwegen der besseren Fehlermeldung undsolvehabe ich als freie Funktion implementiert (ok, LU-Zerlegung war dann doch besser). Template-Aliases werden leider vom g++ 4.7 noch nicht unterstuetzt. Danke fuer eure Antworten.
-
Ich denke, man sollte in dem Zusammenhang unbedingt boost.uBLAS erwähnen

-
dot schrieb:
Ich denke, man sollte in dem Zusammenhang unbedingt boost.uBLAS erwähnen

naja, ich habs in meinem jetzigen Projekt verwendet und aktuell würde ich lieber Eigen verwenden. Ublas optimiert einfach das Falsche. In ublas ist nichtmal
A=prod(B,C)fix. Faktor 2 langsamer als die handgeschriebene Schleife. Ja, ich hab das ausgemessen. Man muss das ekelhafte axpy_prod verwenden um schnell zu sein. Yeah, Syntax wie zu C-Zeiten. Und komplexere Ausdrücke wie:
A = prod(B,C+D+E)fliegen einem komplett um die Ohren, insbesondere wenn man mal teureres als nur ne Addition hat.
uBLAS ist definitiv kein Positivbeispiel für ein gutes LinAlg-Design.@Nexus nein, ein selbst evaluierender Kernel für die Berechnungen, der das Ergebnis direkt in den Ergebnisvektor schreibt ist hier wohl angesagt. Bei Matrizen dynamischer Größe hast du sonst den Allokationsoverhead und statische Arrays (und danach siehts mir hier aus) kannst du nicht gut verschieben.
-
Ja, ich habe auch ueberlegt, Eigen oder MTL zu benutzen. Aber fuer meine einfachen Operationen, wollte ich keine ganze Library haben.