Effiziente Kalkulation von Produkten von Kronecker in C

8

Ich bin ziemlich neu in C, da ich für die meiste Zeit meiner Forschung nicht viel schneller als Python brauche. Es stellt sich jedoch heraus, dass neuere Arbeiten, die ich gemacht habe, die Berechnung von ziemlich großen Vektoren / Matrizen erforderten, und daher könnte eine C + MPI-Lösung in Ordnung sein.

Mathematisch gesehen ist die Aufgabe sehr einfach. Ich habe viele Vektoren der Dimensionalität ~ 40k und möchte das Kronecker Produkt ausgewählter Paare dieser Vektoren berechnen und dann summe diese kronecker-produkte.

Die Frage ist, wie man das effizient macht? Ist etwas falsch mit der folgenden Code-Struktur, die For-Schleifen verwenden oder den Effekt erhalten?

Die unten beschriebene Funktion kron übergibt die Vektoren A und B der Längen vector_size und berechnet ihr kronecker-Produkt, das sie in C , a vector_size*vector_size matrix speichert.

%Vor%

Das scheint mir gut zu sein, und sicherlich (wenn ich keinen dummen Syntaxfehler gemacht habe), das richtige Ergebnis zu erzielen, aber ich habe den schleichenden Verdacht, dass eingebettete for-Schleifen nicht optimal sind. Wenn es noch einen anderen Weg geben sollte, lass es mich wissen. Vorschläge willkommen.

Ich danke Ihnen für Ihre Geduld und Ihre Ratschläge. Noch einmal, ich bin sehr unerfahren mit C, aber Googeln herum hat mir wenig Freude für diese Frage gebracht.

    
Edward Grefenstette 08.02.2011, 21:50
quelle

6 Antworten

3

Für Vektoren mit doppelter Genauigkeit (Einzelpräzision und komplex sind ähnlich) können Sie die BLAS-Routine DGER (Rang-1-Aktualisierung) oder ähnliches verwenden, um die Produkte nacheinander auszuführen, da dies der Fall ist alles auf Vektoren. Wie viele Vektoren multiplizieren Sie? Denken Sie daran, dass das Hinzufügen einer Menge Vektor-Außenprodukte (mit denen Sie die Kronecker-Produkte behandeln können) zu einer Matrix-Matrix-Multiplikation führt, mit der BLAS DGEMM effizient umgehen kann. Möglicherweise müssen Sie jedoch Ihre eigenen Routinen schreiben, wenn Sie tatsächlich Integer-Operationen benötigen.

    
Jeremiah Willcock 08.02.2011, 22:05
quelle
6

Da Ihre Schleifenkörper alle völlig unabhängig sind, gibt es sicherlich eine Möglichkeit, dies zu beschleunigen. Am einfachsten wäre es schon, mehrere Kerne zu nutzen, bevor man an MPI denkt. OpenMP sollte das ganz gut machen.

%Vor%

Dies wird heutzutage von vielen Compilern unterstützt.

Sie könnten auch versuchen, einige gebräuchliche Ausdrücke aus der inneren Schleife zu ziehen, aber ordentliche Compiler, wie gcc, icc oder clang, sollten dies ganz alleine machen:

%Vor%

BTW, Indizierung mit int ist normalerweise nicht das Richtige. size_t ist die korrekte typedef für alles, was mit Indizierung und Größe von Objekten zu tun hat.

    
Jens Gustedt 08.02.2011 22:07
quelle
2

Wenn Ihr Compiler C99 unterstützt (und Sie niemals den gleichen Vektor wie A und B übergeben), sollten Sie in einem C99-unterstützenden Modus kompilieren und Ihre Funktionssignatur ändern in:

%Vor%

Das Schlüsselwort restrict verspricht dem Compiler, dass die Arrays, auf die von A , B und C gezeigt wird, keinen Alias ​​(Overlap) haben. Mit Ihrem Code, wie er geschrieben wurde, muss der Compiler A[i] bei jeder Ausführung der inneren Schleife neu laden, da er konservativ sein muss und davon ausgehen muss, dass Ihre Stores in C[] Werte in A[] ändern können. Unter restrict kann der Compiler davon ausgehen, dass dies nicht passieren wird.

    
caf 09.02.2011 00:18
quelle
2

Lösung gefunden (Danke an @Jeremiah Willcock): GSLs BLAS-Bindungen scheinen so zu sein mach den Trick schön. Wenn wir fortlaufend Paare von Vektoren A und B auswählen und sie zu einer 'laufenden Summe' vector / matrix C hinzufügen, die folgende modifizierte Version der obigen kron-Funktion

%Vor%

entspricht genau der BLAS DGER-Funktion (zugänglich als gsl_blas_dger ), funktional gesehen. Die anfängliche Funktion kron ist DGER mit alpha = 0 und C ist eine nicht initialisierte (zeroed) Matrix / Vektor der korrekten Dimensionalität.

Es stellt sich heraus, dass es am Ende einfacher ist, einfach Python-Bindungen für diese Bibliotheken zu verwenden. Ich denke jedoch, dass ich viel gelernt habe, als ich versuchte, dieses Zeug herauszufinden. Es gibt einige weitere hilfreiche Vorschläge in den anderen Antworten, überprüfen Sie sie, wenn Sie das gleiche Problem haben, mit dem Sie sich befassen müssen. Danke an alle!

    
Edward Grefenstette 09.02.2011 00:34
quelle
1

Dies ist ein häufig genug auftretendes Problem in numerischen Berechnungskreisen, das wirklich die beste Sache wäre, ein gut debuggtes Paket wie Matlab (oder eines seiner Klone für Freie Software ).

Sie könnten wahrscheinlich sogar eine Python-Bindung finden, damit Sie loswerden können von C.

Alles oben genannte wird (wahrscheinlich) schneller sein als Code, der ausschließlich in Python geschrieben wurde. Wenn Sie mehr Geschwindigkeit brauchen, würde ich ein paar Dinge vorschlagen:

    Sehen Sie sich Fortran anstelle von C an. Fortran-Compiler tendieren dazu, numerische Berechnungen besser zu optimieren (eine Ausnahme wäre, wenn Sie gcc verwenden, da sowohl C- als auch Fortran-Compiler dasselbe Backend verwenden).
  1. Erwägen Sie, Ihren Algorithmus zu parallelisieren. Es gibt Varianten von Fortran, von denen ich weiß, dass sie parallele Schleifenanweisungen haben. Ich denke, es gibt einige C-Addons, die das gleiche machen. Wenn Sie einen PC (und Single-Precision) verwenden, können Sie auch die Grafikkarte Ihrer Grafikkarte in Betracht ziehen, die im Grunde genommen ein sehr günstiger Array-Prozessor ist.
T.E.D. 08.02.2011 22:17
quelle
1

Eine weitere leicht zu implementierende Optimierung besteht darin, dass Sie, wenn Sie wissen, dass die innere Dimension Ihrer Arrays durch n teilbar ist, dem Schleifenkörper n Zuweisungsanweisungen hinzufügen und die Anzahl der notwendigen Iterationen mit entsprechenden Änderungen reduzieren zur Schleifenzählung.

Diese Strategie kann verallgemeinert werden, indem eine switch-Anweisung um die äußere Schleife mit Fällen für Array-Größen verwendet wird, die durch zwei, drei, vier und fünf oder was auch immer am üblichsten ist. Dies kann zu einem recht großen Leistungsgewinn führen und ist kompatibel mit den Vorschlägen 1 und 3 zur weiteren Optimierung / Parallelisierung. Ein guter Compiler kann sogar so etwas für Sie erledigen (alias Loop entrolling).

Eine weitere Optimierung wäre, die Zeigerarithmetik zu verwenden, um die Array-Indizierung zu vermeiden. So etwas sollte den Trick machen:

%Vor%

Dadurch wird auch vermieden, mehrfach auf den Wert von A [i] zuzugreifen, indem Sie ihn in einer lokalen Variablen zwischenspeichern, was zu einer geringen Geschwindigkeitssteigerung führen kann. (Beachten Sie, dass diese Version nicht parallelisierbar ist, da sie den Wert der Zeiger ändert, aber beim Schleifenabwickeln immer noch funktioniert.)

    
Keith 10.02.2011 02:17
quelle

Tags und Links