Comment écrire du code dimensionnellement agnostique?

Je me retrouve souvent à écrire du code très similaire pour les versions à une, deux et trois dimensions d'une opération / algorithme donné. La maintenance de toutes ces versions peut devenir fastidieuse. La génération de code simple fonctionne assez bien, mais il semble que l'on pense qu'il doit y avoir un meilleur moyen.

Existe-t-il un moyen relativement simple d'écrire une opération une fois et de la généraliser à des dimensions supérieures ou inférieures?

Un exemple concret est: supposons que j'ai besoin de calculer le gradient d'un champ de vitesse dans l'espace spectral. En trois dimensions, les boucles Fortran ressembleraient à quelque chose comme:

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

où le ddxtableau est correctement défini. (On pourrait également le faire avec des multiplications matricielles.) Le code d'un flux bidimensionnel est presque exactement le même, sauf: la troisième dimension est supprimée des boucles, des index et du nombre de composants. Y a-t-il une meilleure façon d'exprimer cela?

Un autre exemple est: supposons que j'ai des vitesses de fluide définies point par point sur une grille tridimensionnelle. Pour interpoler la vitesse à un emplacement arbitraire (c'est-à-dire ne correspondant pas aux points de la grille), on peut utiliser successivement l'algorithme Neville unidimensionnel sur les trois dimensions (c'est-à-dire la réduction dimensionnelle). Existe-t-il un moyen facile de faire une réduction dimensionnelle étant donné l'implémentation unidimensionnelle d'un algorithme simple?

Vous regardez comment deal.II ( http://www.dealii.org/ ) le fait - là, l'indépendance de dimension est au cœur même de la bibliothèque et est modélisée comme un argument de modèle pour la plupart des types de données. Voir, par exemple, le solveur Laplace indépendant des dimensions dans le didacticiel de l'étape 4:

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

Voir également

https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#why-use-templates-for-the-space-dimension

La question souligne que la plupart des langages de programmation "simples" (C, Fortran, au moins) ne vous permettent pas de le faire proprement. Une contrainte supplémentaire est que vous souhaitez une commodité de notation et de bonnes performances.

Par conséquent, au lieu d' écrire un code spécifique à une dimension, envisagez d'écrire un code qui génère un code spécifique à une dimension. Ce générateur est indépendant des dimensions, même si le code de calcul ne l'est pas. En d'autres termes, vous ajoutez une couche de raisonnement entre votre notation et le code exprimant le calcul. Les modèles C ++ reviennent à la même chose: à l'envers, ils sont intégrés directement dans le langage. Inconvénient, ils sont un peu lourds à écrire. Cela réduit la question de savoir comment réaliser pratiquement le générateur de code.

OpenCL vous permet de générer du code au moment de l'exécution assez proprement. Il permet également une séparation très nette entre «programme de contrôle externe» et «boucles / noyaux internes». Le programme externe de génération est beaucoup moins contraint aux performances, et pourrait donc aussi bien être écrit dans un langage confortable, comme Python. C'est mon espoir pour la façon dont PyOpenCL sera utilisé - désolé pour la fiche sans vergogne renouvelée.

Cela peut être accompli dans n'importe quelle langue avec le prototype mental approximatif suivant:

1. Créez une liste des étendues de chaque dimension (quelque chose comme shape () dans MATLAB je pense)
2. Créez une liste de votre emplacement actuel dans chaque dimension.
3. Écrivez une boucle sur chaque dimension, contenant une boucle sur laquelle la taille change en fonction de la boucle externe.

De là, il s'agit de lutter contre la syntaxe de votre certain langage pour garder votre code conforme nd.

Après avoir écrit un solveur de dynamique des fluides à n dimensions , j'ai trouvé qu'il était utile d'avoir un langage qui supporte le décompactage d'un objet de type liste comme arguments d'une fonction. Soit a = (1,2,3) f (a *) -> f (1,2,3). De plus, des itérateurs avancés (tels que ndenumerate dans numpy) font du code un ordre de grandeur plus propre.

$n_1 \times n_2 \times n_3$$n_j=1$

Les réponses claires si vous voulez garder la vitesse de Fortran sont d'utiliser un langage qui a une bonne génération de code comme Julia ou C ++. Les modèles C ++ ont déjà été mentionnés, je vais donc mentionner les outils de Julia ici. Les fonctions générées par Julia vous permettent d'utiliser sa métaprogrammation pour créer des fonctions à la demande via des informations de type. Donc, essentiellement, ce que vous pouvez faire ici, c'est faire

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

puis vous utilisez le Npour créer par programme le code que vous souhaitez exécuter étant donné qu'il est Ndimensionnel. Ensuite, la bibliothèque cartésienne de Julia ou des packages comme les expressions Einsum.jl peuvent être facilement construits pour la Nfonction dimensionnelle.

Ce qui est bien avec Julia ici, c'est que cette fonction est compilée statiquement et optimisée pour chaque nouveau tableau dimensionnel que vous utilisez, donc elle ne compilera pas plus que ce dont vous avez besoin, mais elle vous donnera la vitesse C / Fortran. En fin de compte, cela est similaire à l'utilisation de modèles C ++, mais c'est un langage de niveau supérieur avec beaucoup d'outils pour le rendre plus facile (assez facile que ce serait un joli problème de devoirs pour un étudiant de premier cycle).

Un autre langage qui est bon pour cela est un Lisp comme Common Lisp. Il est facile à utiliser car comme Julia, il vous donne l'AST compilé avec beaucoup d'outils d'introspection intégrés, mais contrairement à Julia, il ne le compilera pas automatiquement (dans la plupart des distributions).

Je suis dans le même bateau (Fortran). Une fois que j'ai mes éléments 1D, 2D, 3D et 4D (je fais de la géométrie projective), je crée les mêmes opérateurs pour chaque type, puis j'écris ma logique avec des équations de haut niveau qui rendent clair ce qui se passe. Ce n'est pas aussi lent qu'on pourrait le penser d'avoir des boucles distinctes pour chaque opération et beaucoup de copie de mémoire. Je laisse le compilateur / processeur faire les optimisations.

Par exemple

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

À utiliser dans des équations comme

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

où eet ret gpeut avoir toute dimensionnalité qui a un sens mathématique.