Jak napisać kod agnostyczny wymiarowo?

Często piszę bardzo podobny kod dla jedno-, dwu- i trójwymiarowych wersji danej operacji / algorytmu. Utrzymanie wszystkich tych wersji może być nudne. Proste generowanie kodu działa dość dobrze, ale wydaje się, że istnieje lepszy sposób.

Czy istnieje względnie prosty sposób, aby napisać operację raz i generalizować ją do wyższych lub niższych wymiarów?

Jednym konkretnym przykładem jest: załóżmy, że muszę obliczyć gradient pola prędkości w przestrzeni spektralnej. W trzech wymiarach pętle Fortrana wyglądałyby mniej więcej tak:

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

gdzie ddxtablica jest odpowiednio zdefiniowana. (Można to również zrobić z mnożnikami macierzy.) Kod dwuwymiarowego przepływu jest prawie dokładnie taki sam, z wyjątkiem: trzeci wymiar jest usuwany z pętli, indeksów i liczby składników. Czy istnieje lepszy sposób na wyrażenie tego?

Innym przykładem jest: załóżmy, że mam prędkości płynu zdefiniowane punktowo na trójwymiarowej siatce. Aby interpolować prędkość do dowolnej lokalizacji (tj. Nie odpowiadającej punktom siatki), można zastosować jednowymiarowy algorytm Neville'a kolejno we wszystkich trzech wymiarach (tj. Redukcja wymiarów). Czy istnieje prosty sposób na zmniejszenie wymiarów, biorąc pod uwagę jednowymiarową implementację prostego algorytmu?

Patrzysz na to, jak robi to deal.II ( http://www.dealii.org/ ) - tam niezależność wymiarów leży w samym sercu biblioteki i jest modelowana jako szablonowy argument dla większości typów danych. Zobacz, na przykład, niezależny od wymiarów solver Laplace'a w programie samouczka krok 4:

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

Zobacz też

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

Pytanie podkreśla, że ​​większość „prostych” języków programowania (przynajmniej C, Fortran) nie pozwala ci robić tego w czysty sposób. Dodatkowym ograniczeniem jest to, że chcesz notatywnej wygody i dobrej wydajności.

Dlatego zamiast pisać kod specyficzny dla wymiaru, należy rozważyć napisanie kodu generującego kod specyficzny dla wymiaru. Ten generator jest niezależny od wymiarów, nawet jeśli kod obliczeniowy nie jest. Innymi słowy, dodajesz warstwę rozumowania między twoją notacją a kodem wyrażającym obliczenia. Szablony C ++ mają tę samą rzecz: do góry, są wbudowane bezpośrednio w język. Minusem, są nieco nieporęczne w pisaniu. Ogranicza to pytanie do praktycznej realizacji generatora kodu.

OpenCL pozwala generować kod w czasie wykonywania dość czysto. Zapewnia również bardzo czysty podział między „zewnętrznym programem sterującym” a „wewnętrznymi pętlami / jądrem”. Zewnętrzny, generujący program jest znacznie mniej ograniczony pod względem wydajności i dlatego równie dobrze może być napisany w wygodnym języku, takim jak Python. To moja nadzieja na wykorzystanie PyOpenCL - przepraszam za nową bezwstydną wtyczkę.

Można to osiągnąć w dowolnym języku za pomocą następującego szorstkiego prototypu mentalnego:

1. Utwórz listę zakresu każdego wymiaru (myślę, że coś takiego jak kształt () w MATLAB-ie)
2. Utwórz listę bieżącej lokalizacji w każdym wymiarze.
3. Napisz pętlę nad każdym wymiarem, zawierającą pętlę, której rozmiar zmienia się na podstawie zewnętrznej pętli.

Odtąd chodzi o walkę ze składnią określonego języka, aby zachować zgodność z kodem nd.

Po napisaniu n-wymiarowego solwera dynamiki płynów odkryłem, że pomocne jest posiadanie języka, który obsługuje rozpakowywanie listy podobnej do obiektu jako argumentów funkcji. Tj. A = (1,2,3) f (a *) -> f (1,2,3). Dodatkowo zaawansowane iteratory (takie jak ndenumerate in numpy) sprawiają, że kod jest o rząd wielkości czystszy.

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

Jasne odpowiedzi, jeśli chcesz utrzymać szybkość Fortran, to użyć języka, który ma prawidłowe generowanie kodu, takiego jak Julia lub C ++. Wzory C ++ zostały już wspomniane, więc wspomnę tutaj o narzędziach Julii. Funkcje generowane przez Julię umożliwiają wykorzystanie jej metaprogramowania do budowania funkcji na żądanie za pomocą informacji o typie. Zasadniczo możesz tutaj zrobić

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

a następnie używasz go Ndo programowego budowania kodu, który chcesz wykonać, biorąc pod uwagę jego Nwymiar. Następnie kartezjańską bibliotekę Julii lub pakiety, takie jak wyrażenia Einsum.jl, można łatwo zbudować dla Nfunkcji wymiarowej.

To, co jest miłe w Julii, to fakt, że ta funkcja jest kompilowana statycznie i zoptymalizowana dla każdej nowej macierzy wymiarów, której używasz, więc nie skompiluje więcej, niż potrzebujesz, ale zapewni Ci szybkość C / Fortran. W końcu jest to podobne do korzystania z szablonów C ++, ale jest to język wyższego poziomu z wieloma narzędziami, aby to ułatwić (na tyle łatwe, że byłby to miły problem dla studentów.

Kolejnym dobrym językiem jest Lisp, jak Common Lisp. Jest łatwy w użyciu, ponieważ podobnie jak Julia daje skompilowaną AST z wieloma wbudowanymi narzędziami introspekcji, ale w przeciwieństwie do Julii nie skompiluje jej automatycznie (w większości dystrybucji).

Jestem na tej samej łodzi (Fortran). Kiedy mam już elementy 1D, 2D, 3D i 4D (wykonuję geometrię rzutową), tworzę te same operatory dla każdego typu, a następnie piszę moją logikę za pomocą równań wysokiego poziomu, które wyjaśniają, co się dzieje. Nie jest tak wolny, jak mogłoby się wydawać, że ma osobne pętle dla każdej operacji i dużo kopii pamięci. Pozwoliłem kompilatorowi / procesorowi wykonać optymalizacje.

Na przykład

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

Do użycia w równaniach takich jak

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

gdzie ea ri gmoże mieć dowolną wymiarowości sprawia, że matematyczny sens.