Wie schreibe ich dimensionsunabhängigen Code?

19

Oft schreibe ich sehr ähnlichen Code für ein-, zwei- und dreidimensionale Versionen einer bestimmten Operation / eines Algorithmus. Das Verwalten all dieser Versionen kann mühsam werden. Einfache Code-Generierung funktioniert ziemlich gut, aber es scheint, als gäbe es einen besseren Weg.

Gibt es eine relativ einfache Möglichkeit, eine Operation einmal zu schreiben und auf höhere oder niedrigere Dimensionen zu verallgemeinern?

Ein konkretes Beispiel ist: Angenommen, ich muss den Gradienten eines Geschwindigkeitsfeldes im Spektralraum berechnen. In drei Dimensionen würden die Fortran-Schleifen ungefähr so ​​aussehen:

do k = 1, n
  do j = 1, n
    do i = 1, n
      phi(i,j,k) = ddx(i)*u(i,j,k) + ddx(j)*v(i,j,k) + ddx(k)*w(i,j,k)
    end do
  end do
end do

wo das ddxArray entsprechend definiert ist. (Man könnte dies auch mit Matrixmultiplikatoren tun.) Der Code für einen zweidimensionalen Fluss ist nahezu identisch, mit der Ausnahme, dass die dritte Dimension aus den Schleifen, Indizes und der Anzahl der Komponenten entfernt wird. Gibt es eine bessere Möglichkeit, dies auszudrücken?

Ein anderes Beispiel ist: Angenommen, ich habe Fluidgeschwindigkeiten, die punktweise auf einem dreidimensionalen Gitter definiert sind. Um die Geschwindigkeit auf einen beliebigen Ort zu interpolieren (dh nicht den Gitterpunkten zu entsprechen), kann der eindimensionale Neville- Algorithmus nacheinander über alle drei Dimensionen hinweg verwendet werden (dh Dimensionsreduktion). Gibt es eine einfache Möglichkeit zur Dimensionsreduktion bei eindimensionaler Implementierung eines einfachen Algorithmus?

Matthew Emmett
quelle

Antworten:

13

Sie sehen, wie deal.II ( http://www.dealii.org/ ) dies tut - dort liegt die Dimensionsunabhängigkeit im Herzen der Bibliothek und wird als Vorlagenargument für die meisten Datentypen modelliert. Siehe beispielsweise den dimensionsunabhängigen Laplace-Solver im Lernprogramm für Schritt 4:

http://www.dealii.org/developer/doxygen/deal.II/step_4.html

Siehe auch

https://github.com/dealii/dealii/wiki/Häufig gestellte Fragen#why-use-templates-for-the-space-dimension

Wolfgang Bangerth
quelle
Ich bin sehr einverstanden. Ich habe keinen besseren Ansatz gefunden als Deal.II. Sie verwenden Vorlagen auf sehr interessante Weise, um dieses Problem zu umgehen.
Eldila
1
Eine gute Ressource, aber ziemlich einschüchternd, wenn Sie keine C ++ - Vorlagen verwenden.
Meawoppl
@Wolfgang Bangerth: Definiert deal.ii Iteratoren auch über Templates?
Matthew Emmett
@MatthewEmmett: Ja.
Wolfgang Bangerth
@meawoppl: Eigentlich nein. Ich unterrichte regelmäßig Unterrichtsstunden in deal.II und sage den Schülern zu Beginn einfach, dass alles, was ClassWhatever <2> sagt, in 2d steht, ClassWhatever <3> in 3d und ClassWhatever <dim> in dim-d. Ich bringe die Lektion über Vorlagen irgendwann in Woche 3 und während es wahrscheinlich ist, dass Schüler nicht verstehen, wie es vorher funktioniert , sind sie trotzdem voll funktionsfähig , wenn sie es verwenden.
Wolfgang Bangerth
12

Die Frage hebt hervor, dass die meisten "einfachen" Programmiersprachen (zumindest C, Fortran) es Ihnen nicht erlauben, dies sauber zu machen. Eine zusätzliche Einschränkung besteht darin, dass Sie eine einfache Schreibweise und eine gute Leistung wünschen .

Statt also das Schreiben eine Dimension spezifischen Code, sollten Sie schreiben einen Code, erzeugt eine Dimension spezifischen Code. Dieser Generator ist dimensionsunabhängig, auch wenn der Berechnungscode dies nicht ist. Mit anderen Worten, Sie fügen eine Argumentationsebene zwischen Ihrer Notation und dem Code ein, der die Berechnung ausdrückt. C ++ - Vorlagen haben dieselbe Bedeutung: Sie sind direkt in die Sprache integriert. Nachteil, sie sind etwas umständlich zu schreiben. Dies reduziert die Frage, wie der Codegenerator praktisch zu realisieren ist.

Mit OpenCL können Sie Code zur Laufzeit ziemlich sauber generieren. Es sorgt auch für eine sehr saubere Trennung zwischen 'äußerem Steuerungsprogramm' und 'inneren Schleifen / Kerneln'. Das äußere, generierende Programm ist weit weniger leistungsbeschränkt und kann daher genauso gut in einer bequemen Sprache wie Python geschrieben werden. Das ist meine Hoffnung, wie PyOpenCL verwendet wird - Entschuldigung für den erneuerten schamlosen Stecker.

Andreas Klöckner
quelle
Andreas! Willkommen bei Scicomp! Ich bin froh, Sie auf der Seite zu haben. Ich denke, Sie wissen, wie Sie sich mit mir in Verbindung setzen können, wenn Sie Fragen haben.
Aron Ahmadia
2
+10000 für die automatische Codegenerierung als Lösung für dieses Problem anstelle von C ++ Magic.
Jeff
9

Dies kann in jeder Sprache mit dem folgenden groben mentalen Prototyp erreicht werden:

  1. Erstellen Sie eine Liste der Ausmaße jeder Dimension (so etwas wie shape () in MATLAB, denke ich)
  2. Erstellen Sie eine Liste Ihres aktuellen Standorts in jeder Dimension.
  3. Schreiben Sie eine Schleife über jede Dimension, die eine Schleife enthält, deren Größe sich basierend auf der äußeren Schleife ändert.

Von da an geht es darum, gegen die Syntax Ihrer bestimmten Sprache zu kämpfen, um Ihren Code nd-konform zu halten.

Nachdem ich einen n-dimensionalen Fluiddynamik-Löser geschrieben habe , habe ich festgestellt, dass es hilfreich ist, eine Sprache zu haben, die das Entpacken eines listenähnlichen Objekts als Argument einer Funktion unterstützt. Dh a = (1,2,3) f (a *) -> f (1,2,3). Durch erweiterte Iteratoren (z. B. ndenumerate in numpy) wird der Code um eine Größenordnung übersichtlicher.

meawoppl
quelle
Die Python-Syntax dafür sieht nett und prägnant aus. Ich frage mich, ob es einen schönen Weg gibt, dies mit Fortran zu tun ...
Matthew Emmett
1
Es ist ein bisschen schmerzhaft, sich in Fortran mit dem dynamischen Gedächtnis auseinanderzusetzen. Wahrscheinlich meine Hauptbeschwerde mit der Sprache.
Meawoppl
5

n1×n2×n3nj=1

Arnold Neumaier
quelle
Um dimensionsunabhängig zu sein, muss Ihr Code für maxdim + 1-Dimensionen geschrieben werden, wobei maxdim die maximal mögliche Dimension ist, auf die der Benutzer jemals stoßen kann. Angenommen, maxdim = 100. Wie nützlich ist der resultierende Code?
Jeff
4

Die klaren Antworten, wenn Sie die Geschwindigkeit von Fortran beibehalten möchten, sind die Verwendung einer Sprache, die über die richtige Codegenerierung wie Julia oder C ++ verfügt. C ++ - Vorlagen wurden bereits erwähnt, daher erwähne ich hier Julias Tools. Mit den von Julia generierten Funktionen können Sie mithilfe der Metaprogrammierung Funktionen nach Bedarf über Typinformationen erstellen. Im Wesentlichen können Sie hier also Folgendes tun

@generated function f(x)
   N = ndims(x)
   quote
     # build the code for the function
   end
end

Anschließend Nerstellen Sie mit dem Befehl programmgesteuert den Code, den Sie ausführen möchten, da er Ndimensioniert ist. Dann können Julias kartesische Bibliothek oder Pakete wie Einsum.jl- Ausdrücke einfach für die NDimensionsfunktion erstellt werden.

Das Schöne an Julia ist, dass diese Funktion statisch kompiliert und für jedes neue dimensionale Array optimiert wird, sodass nicht mehr kompiliert wird, als Sie benötigen, und Sie trotzdem die C / Fortran-Geschwindigkeit erreichen. Am Ende ähnelt dies der Verwendung von C ++ - Vorlagen, es ist jedoch eine höhere Sprache mit vielen Werkzeugen, die es einfacher machen (leicht genug, dass dies ein nettes Hausaufgabenproblem für einen Studenten wäre).

Eine andere Sprache, die dafür gut ist, ist ein Lisp wie Common Lisp. Es ist einfach zu benutzen, da es Ihnen wie Julia den kompilierten AST mit einer Menge eingebauter Introspection-Tools bietet, aber im Gegensatz zu Julia wird es nicht automatisch kompiliert (in den meisten Distributionen).

Chris Rackauckas
quelle
1

Ich bin im selben (Fortran) Boot. Sobald ich meine 1D-, 2D-, 3D- und 4D-Elemente (ich mache projektive Geometrie) habe, erstelle ich für jeden Typ die gleichen Operatoren und schreibe dann meine Logik mit Gleichungen auf hoher Ebene, die klar machen, was los ist. Es ist nicht so langsam, wie Sie vielleicht glauben, separate Schleifen für jede Operation und viel Speicherkopie zu haben. Ich überlasse es dem Compiler / Prozessor, die Optimierungen vorzunehmen.

Beispielsweise

interface operator (.x.)
    module procedure cross_product_1x2
    module procedure cross_product_2x1
    module procedure cross_product_2x2
    module procedure cross_product_3x3
end interface 

subroutine cross_product_1x2(a,b,c)
    real(dp), intent(in) :: a(1), b(2)
    real(dp), intent(out) :: c(2)

    c = [ -a(1)*b(2), a(1)*b(1) ]
end subroutine

subroutine cross_product_2x1(a,b,c)
    real(dp), intent(in) :: a(2), b(1)
    real(dp), intent(out) :: c(2)

    c = [ a(2)*b(1), -a(1)*b(1) ]
end subroutine

subroutine cross_product_2x2(a,b,c)
    real(dp), intent(in) :: a(2), b(2)
    real(dp), intent(out) :: c(1)

    c = [ a(1)*b(2)-a(2)*b(1) ]
end subroutine

subroutine cross_product_3x3(a,b,c)
    real(dp), intent(in) :: a(3), b(3)
    real(dp), intent(out) :: c(3)

    c = [a(2)*b(3)-a(3)*b(2), a(3)*b(1)-a(1)*b(3), a(1)*b(2)-a(2)*b(1)]
end subroutine

Zur Verwendung in Gleichungen wie

m = e .x. (r .x. g)  ! m = e×(r×g)

wo eund rund gkann jede Dimension haben, die mathematisch Sinn macht.

ja72
quelle