Matrixklasse - kleiner Fehler

RobertB

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.

ScottZhang

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".

cooky451

Du hast uns immer noch nicht verraten warum du C++ benutzen willst wenn du doch matlab hast?

RobertB

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 :).

cooky451

Mein Tipp: Bleib einfach bei Matlab. Nicht, dass ich groß Ahnung davon hätte, aber ich weiß ein paar Dinge:
- Die Leute die Matlab implementiert haben, wissen bestimmt was sie tun.
- Du hast keine Ahnung von C++.
-> Es dauert wahrscheinlich mindestens einige Wochen bis du gut genug bist, um Matlab performancemäßig zu schlagen, auch mit Bibliotheken die für dich das Rechnen übernehmen.

Die zeit ist recht grob geschätzt, da ich weder dich noch Matlab gut genug kenne, aber das sollte man trotzdem mal erwähnen. Abgesehen davon: Wenn du wirklich, wirklich Performance brauchst, warum lädst du dann wieder Sachen nach Matlab? Das bringt doch nichts. Wenn, dann rechne gleiche alles in C++ aus. Nimm Eigen oder MTL und lass die möglichst alles rechnen, was die können. Denn deren Algorithmen und Implementierungen sind mit an Sicherheit grenzender Wahrscheinlichkeit deutlich schneller als alles was du oder wir hier mal eben zusammenbasteln könnten.

RobertB

Das mit Matlab entscheide leider nicht ich allein, sondern meine "Chefs". Aber Danke für Deinen Rat.

Denn, ich habe mir etwas überlegt: ich schmeiße ganz am Anfang in Matlab erstellte Daten in C++ (mit build-in Funktionen von Matlab), und am Ende hole ich mir das Ergebnis und lasse es in Matlab plotten. Dann habe ich nur die beiden Schnittstellen.

Und, was hälts Du davon, dass ich mein ganzes Matlabprogramm am Ende einfach mit einem Coder in C++-Code übersetze? Matlab bietet dafür eigens eine Lösung. Ob das dann schneller ist?

Auf jeden Fall ziehe ich den Hut vor euch C++-Programmieren! Das ist doch die wahre Kunst.

ScottZhang

Das wollt ich gerade vorschlagen. Das moderne Matlab kann C++ Code exportieren (vorausgesetzt ihr habt die entsprechenden Lizenzen).

Ich habe zwar nur seeehhr wenig ahnung vom mex-Interface. Aber könntest du nicht rauskriegen wie Matlab die Sparse-Matrizen implementiert (ich tippe auf eine von den von mir vorgeschlagenden Formaten) und einfach für diese Implementierung die entsprechend Funktionen schreiben die du brauchst (Addition , Multiplikation, Assemblierung und was sonst noch anfällt. Bei den ersten beiden würd ich ma vermuten das sind schon build-ins).

Werner Salomon

ScottZhang schrieb:

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".

Dass es dafür auch Lösungen gibt, die nur mit Arrays und Indizes funktionieren war ja klar; die gute alte Fortran-Welt lässt grüßen. Aber dass es eine Datenstruktur als Quasi-Standard bis in die Wikipedia geschafft; damit hätte ich jetzt nicht gerechnet.
Ich habe nicht herausbekommen, ob MatLab die Sparse matrix im CRS abspeichert, aber ich habe mir mal die Mühe gemacht, meine class SparseMatrix auf dieses Format umzustellen.
Es ist nicht wirklich mehr Code geworden, aber es war schwieriger zu implementieren. Ich habe bestimmt mehr als die doppelte Zeit gebraucht wie für die std::map-Lösung.

Keine Ahnung, ob das anderen auch so geht, aber mit den STL-Container kann man auf einem höheren Abstraktionsniveau codieren, und das geht bei mir wesentlich flüssiger.

Die class SparseMatrixc im CRS Datenformat:

#include <algorithm> // for_each, find_if,..
#include <cassert>
#include <iterator> // distance
#include <ostream>
#include <vector>
#include <utility> // std::pair<>

template< typename T >
class SparseMatrix
{
    struct Ref
    {
        Ref( SparseMatrix& mx, std::size_t rowIdx, std::size_t columnIdx )
            : mx_( &mx )
            , rowIdx_( rowIdx )
            , columnIdx_( columnIdx )
        {}
        Ref& operator=( const T& val )
        {
            auto element = mx_->find_element( rowIdx_, columnIdx_ );
            if( element.second )
            {
                if( (mx_->val_[element.first] = val) == T() )
                    mx_->erase_( rowIdx_, element.first );
            }
            else if( val != T() )
                mx_->insert_( rowIdx_, columnIdx_, element.first, val );
            return *this;
        }
        operator T() const
        {
            auto element = mx_->find_element( rowIdx_, columnIdx_ );
            return element.second? mx_->val_[element.first]: T();
        }
    private:
        SparseMatrix* mx_;
        std::size_t rowIdx_;
        std::size_t columnIdx_;
    };
    friend struct Ref;

public:
    SparseMatrix( std::size_t rows, std::size_t columns )
        : columns_( columns )
        , val_()
        , colInd_()
        , rowPtr_( rows+1 )
    {}
    std::size_t Columns() const { return columns_; }
    std::size_t Rows() const { return rowPtr_.size()-1; }

    // --  Element Zugriff
    T operator()( std::size_t rowIdx, std::size_t columnIdx ) const
    {
        assert( rowIdx < Rows() );
        assert( columnIdx < Columns() );
        auto element = find_element( rowIdx, columnIdx );
        return element->second? val_[element->first]: T();
    }
    Ref operator()( std::size_t rowIdx, std::size_t columnIdx )
    {
        assert( rowIdx < Rows() );
        assert( columnIdx < Columns() );
        return Ref( *this, rowIdx, columnIdx );
    }

    // --  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;
        std::size_t i = 0; // Index in mx.colInd_[] & mx.val_[]
        for( auto iRow = begin(mx.rowPtr_); iRow != end(mx.rowPtr_)-1; ++iRow )
        {
            for( std::size_t columnIdx = 0; columnIdx < mx.columns_; ++columnIdx )
            {
                (out << " ").width( width );
                if( i < *(iRow+1) && mx.colInd_[i] == columnIdx )
                    out << mx.val_[i++];
                else
                    out << T();
            }
            out << "\n";
        }
        return out << std::flush;
    }

    std::size_t NoElements() const // Anzahl belegter Elemente; für Debugzwecke
    {
        return val_.size();
    }

private:
    void insert_( std::size_t rowIdx, std::size_t columnIdx, std::size_t i, const T& val )
    {
        assert( rowPtr_[rowIdx] <= i );
        assert( i <= rowPtr_[rowIdx+1] );
        assert( std::none_of( begin(colInd_) + rowPtr_[rowIdx], begin(colInd_) + rowPtr_[rowIdx+1],
            [columnIdx]( std::size_t c )->bool { return c == columnIdx; } ) );
        if( val != T() )
        {
            val_.insert( begin(val_) + i, val );
            colInd_.insert( begin(colInd_) + i, columnIdx );
            std::for_each( begin(rowPtr_) + rowIdx+1, end(rowPtr_), []( std::size_t& idx ) { ++idx; } );
        }
    }
    void erase_( std::size_t rowIdx, std::size_t i )
    {
        assert( i < val_.size() );
        assert( rowPtr_[rowIdx] <= i );
        assert( i < rowPtr_[rowIdx+1] );
        val_.erase( begin(val_) + i );
        colInd_.erase( begin(colInd_) + i );
        std::for_each( begin(rowPtr_) + rowIdx+1, end(rowPtr_), []( std::size_t& idx ) { --idx; } );
    }
    std::pair< std::size_t, bool > find_element( std::size_t rowIdx, std::size_t columnIdx ) const
    {
        auto iColumn = std::find_if( begin(colInd_) + rowPtr_[rowIdx], begin(colInd_) + rowPtr_[rowIdx+1],
            [columnIdx]( std::size_t c )->bool { return columnIdx <= c; } );
        return std::make_pair( std::size_t( std::distance( begin(colInd_), iColumn ) )
            , iColumn != begin(colInd_) + rowPtr_[rowIdx+1] && *iColumn == columnIdx );
    }

    // --   Members
    std::size_t columns_;
    // --   Datenstruktur der Sparse Matrix im CRS -> <http://de.wikipedia.org/wiki/Compressed_Row_Storage>
    std::vector< T > val_;
    std::vector< std::size_t > colInd_;
    std::vector< std::size_t > rowPtr_;
};

Der Speicherbedarf dürfte niedriger sein (ca. bei 20%-30% der vorherigen Lösung) dafür ist das Einfügen neuer Elemente nicht mehr performant. Zumindest nicht wenn man sie einzeln zuweist.

.. es fehlt noch die Addition. Wenn ihr das Demo-main aus meinem vorherigen Posting benutzt, so lasst die Addition einfach weg, oder schreibt Euch eine - viel Spaß damit.

Gruß
Werner

ScottZhang

Sieht wohl so aus als würde Matlab die columns compressen: Matlab: mxCreateSparse

otze

CSR ist auch in ublas implmentiert. das Problem an CSR ist wie gesagt das Einfügen bei matrix-vektor und matrix-matrix produkten.

Ein netter Trick dafür ist etwas mehr Speicher als benötigt zu reservieren und die Elemente bei der Berechnung des Produkts am Ende unsortiert einzufügen. Am Schluss wird dann gemerged.