Matrixklasse - kleiner Fehler
-
Hallo.
Ich habe folgende Matrixklasse im Internet entdeckt. Sie läuft ohne Probleme durch, nur wie geschrieben, sind die Ergebnisse falsch. Beim zweiten Beitrag weiß ich nicht, wo ich das einbetten soll. Hier erstmal der Link:
http://www.cplusplus.com/forum/general/17908/
Ich habe es z.B. für die Summe mit folgendem versucht:
//To compute the summation matrix matrix::operator+(matrix& a){ matrix sum(nrows, ncols); for (int i=0; i<nrows; i++) for(int j=0; j<ncols; j++) { this->operator()(i,j); // Das habe ich hinzugefügt sum(i,j) = (i,j) + a(i,j); } return sum; }Das gibt aber keine Veränderung. Ich hoffe, ihr könnt mir helfen. Die Klasse an sich find ich nämlich ganz gut und einfach aufgebaut.
-
wenn ich mir das so auf den ersten blick angucke, würde ich dir empfehlen die finger davon zu lassen
schreib dir selbst ne kleine klasse für matrizen rechnung, ist für dich ne gute übung. sollte mindestens genausoviel können wie die, die du uns gezeigt hast. und danach weist du, dass sie klappt.
wenn mehr willst oder besseres, dann google nach ner richtigen mathe bibliothek, die solche klassen anbietet...
-
Die Klasse ist schlecht. Die Hälfte der Member hier kann wahrscheinlich in 20 Minuten eine (viel, viel) bessere schreiben.
Wozu brauchst du diese Klasse btw? Spiele? Sehr große Matrizen? Ist Performance wichtig? Sicher, dass dir die naiven Algorithmen reichen?
-
Aber warum ist die denn schlecht? Sehe doch da eig. nur das nötigste.
Wichtig sind für mich große Matrizen und Performance.Und weitere Algorithmen würde ich dann selbst dazuschreiben. Nur komm ich von Matlab, und da brauche ich ja das ganze nicht. Daher ist für mich das ganze Neuland und hätte gerne eine Vorlage gehabt.
Könnt ihr mir aber trotzdem nicht kurz sagen, wo man den Befehl einbetten müsste?
-
Wenn du wirklich Performance brauchst, dann nimm eine der großen Mathebibliotheken. (Eigen, MTL, http://scicomp.stackexchange.com/questions/351/recommendations-for-a-usable-fast-c-matrix-library)
Dass die Klasse schlecht ist sieht man schon daran, dass sie kein Template ist, int statt size_t und rohe Pointer statt std::vector benutzt werden. Das ist nur der erste kleine Blick.
Ansonsten hier eine Vorlage für eine ordentliche Klasse:
template <typename T> class matrix { std::vector<T> data_; std::size_t y_size_; public: matrix(std::size_t x_size, std::size_t y_size) : data_(x_size * y_size) , y_size_(y_size) {} const T& operator () (std::size_t i, std::size_t j) const { return data_[i * y_size_ + j]; } T& operator () (std::size_t i, std::size_t j) { return data_[i * y_size_ + j]; } };Tadaa.. Operatoren kannste selbst basteln.
-
Danke schonmal. Welches von den ganzen würdest Du denn nehmen? Du kennst Dich da ja besser aus als ich.
-
Nä, ich hab keine Ahnung von Matrixlibs im wissenschaftlichen Umfeld, ich weiß nur wo ich suchen muss. Aber ich würde MTL wohl als erstes testen.
-
Die gibs wohl nur als Testversion. Aber werde dann mal die anderen versuchen. Am liebsten würd ich das ja selbst programmieren. Nur, das wäre kein optimaler Code. Und ich leide leider unter etwas Zeitdruck. Ich habe zwar ein Grundlagenbuch durchgemacht und weiß meistens auch, von was die Rede ist. Aber alles selbst programmieren, ne, das kann ich noch nicht.
Ich habe dazu noch eine ganz andere Frage: Ist es eine gute Idee, Matrizen als einfache Vektoren zu übergeben? Ich meine damit, dass ich die einzelnen Zeilen der Matrix nacheinanderhäng, und z.B. mit
#define A(Reihe, Spalte) = A[Reihe + irgendwas*Spalte]dass dann definiere?
-
Keine Ahnung was du meinst, aber #define ist (fast) nie eine gute Idee.
- Du stehst unter Zeitdruck? Was genau musst du denn eigentlichen machen? (Wie groß sind deine Matrizen?)
-
Muss viel mit Matrizen hantieren, die wiederrum Ein- und Ausgabe von Funktionen sind. Habe mir auch schon mal das mex-Interface angeschaut. Möchte Matlab und C++ koppeln.
Die Matrizen sind, mh... in Matlab benutze ich schon Sparematrizen, weil diese schon 240.000x240.000 groß sein können. Das ist aber eher schon das Maximale, was ich benötige.
Und wegen den Templates: wenn ich weiß, dass ich nur double-Zahlen habe, kann ich mir die ja sparen. Für mich sind Templates da nur logisch, wenn ich ab und zu auch Integer-Zahlen haben, und damit weniger Speicherplatz benötige.
-
Bei der Matrixgröße kannst du jede konventionelle Matrixklasse vergessen. Überleg doch mal, wie viel Speicher das ist. 460800000000 Byte = 429 Gigabyte. Hast du so viel RAM? Wohl nicht. Also das was du da aus dem Forum kopiert hast (und leider auch meine Variante), sind schon mal völlig aus dem Spiel. Die großen Bibliotheken sollten allerdings uU damit zurecht kommen, die müssten das komprimieren. Was genau du machen willst, und wozu du C++ brauchst wenn du Matlab hast, hast du uns allerdings noch nicht verraten.
-
RobertB schrieb:
Und wegen den Templates: wenn ich weiß, dass ich nur double-Zahlen habe, kann ich mir die ja sparen. Für mich sind Templates da nur logisch, wenn ich ab und zu auch Integer-Zahlen haben, und damit weniger Speicherplatz benötige.
Vllt, möchtes du ja irgendwann ma was genaueres benutzen als double ... Dann fängst du wieder von vorne an?
-
cooky451 schrieb:
Wenn du wirklich Performance brauchst, dann nimm eine der großen Mathebibliotheken. (Eigen, MTL, http://scicomp.stackexchange.com/questions/351/recommendations-for-a-usable-fast-c-matrix-library)
Dass die Klasse schlecht ist sieht man schon daran, dass sie kein Template ist, int statt size_t und rohe Pointer statt std::vector benutzt werden. Das ist nur der erste kleine Blick.
Ansonsten hier eine Vorlage für eine ordentliche Klasse:
template <typename T> class matrix { std::vector<T> data_; std::size_t y_size_; public: matrix(std::size_t x_size, std::size_t y_size) : data_(x_size * y_size) , y_size(y_size) {} const T& operator () (std::size_t i, std::size_t j) const { return data_[i * y_size_ + j]; } T& operator () (std::size_t i, std::size_t j) { return data_[i * y_size_ + j]; } };Tadaa.. Operatoren kannste selbst basteln.
Ich habe mal diese Klasse ausprobiert, aber es tritt ein Fehler auf. Folgendes habe ich gemacht:
#include <iostream> using namespace std; template <typename T> class matrix { vector<T> data_; size_t y_size_; public: matrix(size_t x_size, size_t y_size) : data_(x_size * y_size) , y_size(y_size) {} const T& operator () (size_t i, size_t j) const { return data_[i * y_size_ + j]; } T& operator () (size_t i, size_t j) { return data_[i * y_size_ + j]; } }; int main() { return 0; };Könnt ihr mir bitte sagen, was ich falsch gemacht habe?
-
<vector> nicht #included. https://ideone.com/p6RKY0
Und meinen Beitrag oben nicht gelesen. Bei der Matrixgröße kannst du diesen Ansatz sicher vergessen.
-
cooky451 schrieb:
<vector> nicht #included. https://ideone.com/p6RKY0
Und meinen Beitrag oben nicht gelesen. Bei der Matrixgröße kannst du diesen Ansatz sicher vergessen.Deinen Beitrag habe ich gelesen. Der Ansatz war nur aus Interesse. Aber dieser ist dann allen Anschein nicht geeignet.
-
Hallo RobertB,
bei der angestrebten Größe der Matrix kommst Du nur mit einer Sparse Matrix weiter - vorausgesetzt natürlich, Dein Problem, welches Du lösen willst, erzeugt auch nur dünn besetzte Matrizen.
boost.ublas bietet sparse matrix an. Ich selber habe aber keine Erfahrungen damit.Wenn Du an einer einfach gestrickten Lösung interessiert bist, so kannst Du Dir auch selber was bauen. Im Prinzip lässt sich so eine Matrix mit einem vector von maps lösen. So wäre zumindest die Implementierung relativ straight forward.
#include <iostream> #include <map> #include <vector> #include <cassert> #include <iomanip> // setw #include <numeric> // accumulate template< typename T > class SparseMatrix { typedef std::map< std::size_t, T > Row; // key -> column typedef std::vector< Row > Field; // index -> row struct Ref { Ref( Row& row, std::size_t columnIdx ) : row_( &row ) , column_( columnIdx ) {} Ref& operator=( const T& element ) { if( element == T() ) { auto iElement = row_->find( column_ ); if( iElement != end(*row_) ) row_->erase( iElement ); } else (*row_)[column_] = element; return *this; } operator T() const { auto iElement = row_->find( column_ ); return iElement == end(*row_)? T(): iElement->second; } private: Row* row_; std::size_t column_; }; public: SparseMatrix( std::size_t rows, std::size_t columns ) : columns_( columns ) , field_( rows ) {} std::size_t Columns() const { return columns_; } std::size_t Rows() const { return field_.size(); } // -- Element Zugriff T operator()( std::size_t rowIdx, std::size_t columnIdx ) const { assert( rowIdx < Rows() ); assert( columnIdx < Columns() ); auto iElement = field_[rowIdx].find( columnIdx ); return iElement == end(field_[rowIdx])? T(): iElement->second; } Ref operator()( std::size_t rowIdx, std::size_t columnIdx ) { assert( rowIdx < Rows() ); assert( columnIdx < Columns() ); return Ref( field_[rowIdx], columnIdx ); } // -- Arithmetik SparseMatrix& operator+=( const SparseMatrix& b ); SparseMatrix& operator-=( const SparseMatrix& b ); // -- Ausgabe template< typename E, typename Traits > friend std::basic_ostream< E, Traits >& operator<<( std::basic_ostream< E, Traits >& out, const SparseMatrix& mx ) { std::streamsize width = out.width(); if( width == 0 ) width = 4; for( auto iRow = begin(mx.field_); iRow != end(mx.field_); ++iRow ) { auto iElement = begin(*iRow); for( std::size_t columnIdx = 0; columnIdx < mx.columns_; ++columnIdx ) { out << " " << std::setw( width ); if( iElement != end(*iRow) && iElement->first == columnIdx ) { out << iElement->second; ++iElement; } else out << T(); } out << "\n"; } return out << std::flush; } std::size_t NoElements() const // Anzahl belegter Elemente; für Debugzwecke { return std::accumulate( begin(field_), end(field_), std::size_t(0), []( std::size_t sum, const Row& row )-> std::size_t { return sum + row.size(); } ); } private: std::size_t columns_; Field field_; };Was bei dieser Matrixgröße ein No-GO ist, wenn Du jetzt zwei Schleifen baust, die über alle Zeilen und Spalte laufen und dann irgendwelche Operationen ausführen. Bei einer Kantenlänge von 240.000 Elementen einer quadratischen Matrix schlägt jede NanoSekunde der Operation mit einer knappen Minute Laufzeit zu Buche.
D.h. man muss sich auch bei der Operation die dünne Besetzung der Matrix zu Nutze machen.Als Demonstration habe ich die Addition zweier Matrizen implementiert:
template< typename T > SparseMatrix< T >& SparseMatrix< T >::operator+=( const SparseMatrix< T >& b ) { assert( Rows() == b.Rows() ); assert( Columns() == b.Columns() ); auto iRowDst = begin(field_); for( auto iRow = begin(b.field_); iRow != end(b.field_); ++iRow, ++iRowDst ) { auto iElementDst = begin(*iRowDst); for( auto iElement = begin(*iRow); iElement != end(*iRow); ++iElement ) { for( ; iElementDst != end(*iRowDst) && iElementDst->first < iElement->first; ++iElementDst ) ; if( iElementDst != end(*iRowDst) && iElementDst->first == iElement->first ) { if( (iElementDst->second += iElement->second) == T() ) // Auslöschung beachten iElementDst = iRowDst->erase( iElementDst ); else ++iElementDst; } else iRowDst->insert( iElementDst, *iElement ); } } return *this; }Die Subtraktion ist identisch - nur in Zeile 15 das += durch ein -= ersetzen.
Bei mir (VS10) benötigt eine leere map<size_t,double> 20 Bytes Speicher. Damit hat eine leere SparseMatrix mit 240.000 Zeilen einem Speicherbedarf von knapp 5MByte. Für jedes Element kommen dann in der Größenordnung 20Byte hinzu (ein Knoten in der Map).
Wahrscheinlich ist es auch schlauer die Map durch eine sortierte Liste von pairstd::size_t,T (evt. forward_list) zu ersetzen. Aber mit der Map war es einfacher zu codieren.Abschließend eine Demo, die die Funktion zeigt.
int main() { using namespace std; SparseMatrix< double > mx( 8, 5 ); mx( 0, 1 ) = 3.4; mx( 7, 4 ) = 99.; cout << fixed << setprecision( 1 ); cout << "\n" << mx << "Anzahl Elemente !=0: " << mx.NoElements() << endl; for( std::size_t c = 0; c < mx.Columns(); ++c ) mx( 6, c ) = (c+1)*1.1; cout << "\n" << mx << "Anzahl Elemente !=0: " << mx.NoElements() << endl; mx( 6, 1 ) = 0.; cout << "\n" << mx << "Anzahl Elemente !=0: " << mx.NoElements() << endl; cout << "Der Wert von mx(6,2) ist: " << mx(6,2) << endl; SparseMatrix< double > m2( 8, 5 ); m2(2,0) = 2.4; m2(6,2) = 3.3; m2(6,4) = -5.5; mx += m2; cout << "\n" << mx << "Anzahl Elemente !=0: " << mx.NoElements() << endl; return 0; }Wie viele Elemente!=0 enthält Deine Matrix so in etwa? Und welche Operationen benötigst Du für die Matrix?
Gruß
Werner
-
Wow, ich bin beeindruckt vom Code! Vielen Dank!
Dünn besetze Matrizen habe ich. Muss dann noch ein Gleichungssystem lösen. Aber dafür werde ich wahrscheinlich Matlab nutzen und das rüberziehen. Oder ich benutze mit dem Mex-Interface eine Build-in Funktion von Matlab in C++ für das Lösen.
Weitere Operationen brauche ich eig. nicht. Es geht dann doch hauptsächlich um das Füllen und Bearbeiten einzelner Einträge dieser großen Matrix. Das geht ja alles mit dem Code, vielen Dank. Ich muss diesen nur noch ganz verstehen :).
Und leider brauche ich eine ganz große Schleife. Die erste geht nämlich über die Dimension der Länge bzw. Breite (habe eine quadratische Matrix) der Matrix, die zweite dann von 1 bis 3. Also im schlimmsten Fall habe ich 240.000*3 Durchläufe.
-
Wenn du nen matlab interface hast bietet sich Werners Implementierung aus Performance Gründen nicht an. Soweit ich weiss brauchst du für Matlab alles in einem zusammenhängenden Speicher. D.h. ein haufen Kopieroperation werden fällig.
(nanosekunden und so ... :D)Google einfach mal nach "Compressed row format" bzw. "Compressed column format".
-
Du hast uns immer noch nicht verraten warum du C++ benutzen willst wenn du doch matlab hast?
-
cooky451 schrieb:
Du hast uns immer noch nicht verraten warum du C++ benutzen willst wenn du doch matlab hast?
Aus Performancegründen soll ich C++ benutzen. Man sagte mir, dass wegen den Schleifen C++ schneller ist. Aber wegen den Kopiervorgängen habe ich mir auch schon "Sorgen" gemacht, wie das beim MEX-Interface abläuft.
Ich soll das Programm so schnell wie möglich machen, das ist das Ziel. Ich würde lieber bei Matlab bleiben :).