Como escrevo código agnóstico dimensionalmente?

Costumo encontrar-me escrevendo código muito semelhante para uma, duas e três versões dimensionais de uma dada operação / algoritmo. Manter todas essas versões pode se tornar tediosa. A geração simples de código funciona razoavelmente bem, mas parece que deve haver uma maneira melhor.

Existe uma maneira relativamente simples de escrever uma operação uma vez e generalizá-la para dimensões superiores ou inferiores?

Um exemplo concreto é: suponha que eu precise calcular o gradiente de um campo de velocidade no espaço espectral. Em três dimensões, os loops Fortran se pareceriam com:

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

onde a ddxmatriz é definida adequadamente. (Também é possível fazer isso com multiplicações de matrizes.) O código para um fluxo bidimensional é quase exatamente o mesmo, exceto: a terceira dimensão é eliminada dos loops, índices e número de componentes. Existe uma maneira melhor de expressar isso?

Outro exemplo é: suponha que eu tenha velocidades de fluido definidas ponto a ponto em uma grade tridimensional. Para interpolar a velocidade para um local arbitrário (isto é, não corresponde aos pontos da grade), pode-se usar o algoritmo unidimensional de Neville sucessivamente nas três dimensões (isto é, redução dimensional). Existe uma maneira fácil de fazer redução dimensional, dada a implementação unidimensional de um algoritmo simples?

Você vê como o deal.II ( http://www.dealii.org/ ) faz isso - aí, a independência de dimensão está no coração da biblioteca e é modelada como um argumento de modelo para a maioria dos tipos de dados. Consulte, por exemplo, o solucionador de Laplace independente de dimensão no programa do tutorial da etapa 4:

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

Veja também

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

A questão destaca que a maioria das linguagens de programação "simples" (C, Fortran, pelo menos) não permite que você faça isso corretamente. Uma restrição adicional é que você deseja conveniência notável e bom desempenho.

Portanto, em vez de escrever um código específico da dimensão, considere escrever um código que gere um código específico da dimensão. Este gerador é independente da dimensão, mesmo que o código de computação não seja. Em outras palavras, você adiciona uma camada de raciocínio entre sua notação e o código que expressa a computação. Os modelos C ++ equivalem à mesma coisa: de cabeça para baixo, eles estão embutidos na linguagem. Desvantagem, eles são um tanto complicados de escrever. Isso reduz a questão de como realizar praticamente o gerador de código.

O OpenCL permite gerar o código em tempo de execução de maneira bastante limpa. Ele também cria uma divisão muito clara entre o 'programa de controle externo' e os 'loops / kernels internos'. O programa externo de geração é muito menos restrito ao desempenho e, portanto, também pode ser escrito em uma linguagem confortável, como Python. Essa é a minha esperança de como o PyOpenCL será usado - desculpe pelo plugue descarado e renovado.

Isso pode ser realizado em qualquer idioma com o seguinte protótipo mental bruto:

1. Crie uma lista das extensões de cada dimensão (algo como shape () no MATLAB, eu acho)
2. Crie uma lista da sua localização atual em cada dimensão.
3. Escreva um loop sobre cada dimensão, contendo um loop sobre quem muda de tamanho com base no loop externo.

A partir daí, é uma questão de combater a sintaxe do seu idioma específico para manter seu código e compatível.

Depois de escrever um solucionador de dinâmica de fluidos n-dimensional , descobri que é útil ter uma linguagem que suporte a descompactação de uma lista como objeto como argumento de uma função. Ou seja, a = (1,2,3) f (a *) -> f (1,2,3). Além disso, iteradores avançados (como ndenumerate in numpy) tornam o código uma ordem de magnitude mais limpa.

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

As respostas claras, se você deseja manter a velocidade do Fortran, devem usar uma linguagem com geração de código adequada, como Julia ou C ++. Os modelos C ++ já foram mencionados, portanto, mencionarei as ferramentas de Julia aqui. As funções geradas por Julia permitem que você use sua metaprogramação para criar funções sob demanda por meio de informações de tipo. Então, essencialmente, o que você pode fazer aqui é fazer

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

e então você usa o Npara criar programaticamente o código que você deseja executar, pois é Ndimensional. A biblioteca ou pacotes cartesianos de Julia, como as expressões Einsum.jl, podem ser facilmente construídas para a Nfunção dimensional.

O que há de bom em Julia aqui é que essa função é estaticamente compilada e otimizada para cada nova matriz dimensional que você usa, portanto não compilará mais do que o necessário, mas obterá a velocidade do C / Fortran. No final, isso é semelhante ao uso de modelos C ++, mas é uma linguagem de nível superior com muitas ferramentas para facilitar (fácil o suficiente para que este seja um bom problema de lição de casa para um estudante de graduação).

Outra linguagem que é boa para isso é um Lisp como o Common Lisp. É fácil de usar, pois, como Julia, fornece o AST compilado com muitas ferramentas de introspecção integradas, mas, ao contrário de Julia, não o compila automaticamente (na maioria das distribuições).

Estou no mesmo barco (Fortran). Depois de ter meus elementos 1D, 2D, 3D e 4D (eu faço geometria projetiva), crio os mesmos operadores para cada tipo e depois escrevo minha lógica com equações de alto nível que deixam claro o que está acontecendo. Não é tão lento quanto você imagina ter loops separados para cada operação e muita cópia da memória. Deixei o compilador / processador fazer as otimizações.

Por exemplo

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

Para ser usado em equações como

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

onde ee re gpode ter qualquer dimensionalidade que faz sentido matemático.