Bagaimana cara menulis kode agnostik dimensi?

Saya sering menemukan diri saya menulis kode yang sangat mirip untuk versi satu, dua, dan tiga dimensi dari operasi / algoritma yang diberikan. Mempertahankan semua versi ini bisa menjadi membosankan. Pembuatan kode sederhana berfungsi cukup baik, tetapi tampaknya harus ada cara yang lebih baik.

Apakah ada cara yang relatif sederhana untuk menulis operasi sekali dan memilikinya generalisasi ke dimensi yang lebih tinggi atau lebih rendah?

Salah satu contoh konkret adalah: misalkan saya perlu menghitung gradien bidang kecepatan dalam ruang spektral. Dalam tiga dimensi, loop Fortran akan terlihat seperti:

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

di mana ddxarray didefinisikan dengan tepat. (Kita juga bisa melakukan ini dengan mengalikan matriks.) Kode untuk aliran dua dimensi hampir persis sama, kecuali: dimensi ketiga dijatuhkan dari loop, indeks, dan jumlah komponen. Apakah ada cara yang lebih baik untuk mengekspresikan ini?

Contoh lain adalah: misalkan saya memiliki kecepatan fluida yang ditentukan titik-bijaksana pada jaringan tiga dimensi. Untuk menginterpolasi kecepatan ke lokasi yang sewenang-wenang (yaitu, tidak sesuai dengan titik-titik grid), seseorang dapat menggunakan algoritma Neville satu dimensi secara berturut-turut di ketiga dimensi (yaitu, pengurangan dimensi). Apakah ada cara mudah untuk melakukan reduksi dimensi mengingat implementasi satu dimensi dari algoritma sederhana?

Anda melihat bagaimana deal.II ( http://www.dealii.org/ ) melakukannya - di sana, dimensi independensi terletak di jantung perpustakaan, dan dimodelkan sebagai argumen templat untuk sebagian besar tipe data. Lihat, misalnya, pemecah dimensi-agnostik Laplace dalam program tutorial langkah-4:

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

Lihat juga

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

Pertanyaannya menyoroti bahwa sebagian besar bahasa pemrograman "biasa" (C, Fortran, setidaknya) tidak memungkinkan Anda untuk melakukan ini dengan bersih. Kendala tambahan adalah bahwa Anda menginginkan kenyamanan notasi dan kinerja yang baik.

Oleh karena itu, alih-alih menulis kode dimensi spesifik, pertimbangkan untuk menulis kode yang menghasilkan kode dimensi spesifik. Generator ini adalah dimensi-independen, bahkan jika kode komputasi tidak. Dengan kata lain, Anda menambahkan lapisan alasan antara notasi Anda dan kode yang mengekspresikan perhitungan. C ++ templat sama dengan hal yang sama: Terbalik, mereka dibangun langsung ke dalam bahasa. Kelemahannya, mereka agak sulit untuk ditulis. Ini mengurangi pertanyaan untuk bagaimana mewujudkan pembuat kode secara praktis.

OpenCL memungkinkan Anda melakukan pembuatan kode pada waktu berjalan dengan cukup bersih. Itu juga membuat pemisahan yang sangat bersih antara 'program kontrol luar' dan 'loop internal / kernel'. Bagian luar, program penghasil jauh lebih sedikit kendala kinerja, dan oleh karena itu sebaiknya ditulis dalam bahasa yang nyaman, seperti Python. Itulah harapan saya untuk bagaimana PyOpenCL akan digunakan - maaf untuk plug tak tahu malu yang diperbarui.

Ini dapat dilakukan dalam bahasa apa pun dengan prototipe mental kasar berikut:

1. Buat daftar luasan dari setiap dimensi (sesuatu seperti bentuk () di MATLAB, saya pikir)
2. Buat daftar lokasi Anda saat ini di setiap dimensi.
3. Tulis loop di atas setiap dimensi, yang berisi loop di atas perubahan ukuran siapa yang didasarkan pada loop luar.

Dari sana, ini adalah pertanyaan untuk melawan sintaks dari bahasa tertentu Anda untuk menjaga kode Anda nd-compliant.

Setelah menulis pemecah dinamika dinamika fluida n-dimensi , saya menemukan bahwa sangat membantu untuk memiliki bahasa yang mendukung membongkar daftar seperti objek sebagai argumen dari suatu fungsi. Yaitu a = (1,2,3) f (a *) -> f (1,2,3). Selain itu iterator canggih (seperti ndenumerate in numpy) membuat kode urutan besarnya bersih.

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

Jawaban yang jelas jika Anda ingin menjaga kecepatan Fortran adalah menggunakan bahasa yang memiliki generasi kode yang tepat seperti Julia atau C ++. Templat C ++ telah disebutkan, jadi saya akan menyebutkan alat Julia di sini. Fungsi-fungsi Julia yang dihasilkan memungkinkan Anda menggunakan metaprogramming-nya untuk membangun fungsi sesuai permintaan melalui informasi jenis. Jadi pada dasarnya apa yang dapat Anda lakukan di sini adalah lakukan

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

dan kemudian Anda menggunakan Nuntuk secara program membangun kode yang ingin Anda jalankan mengingat bahwa itu Ndimensi. Kemudian pustaka Cartesian atau paket Julia seperti ekspresi Einsum.jl dapat dengan mudah dibangun untuk Nfungsi dimensi.

Apa yang baik tentang Julia di sini adalah bahwa fungsi ini dikompilasi dan dioptimalkan secara statis untuk setiap array dimensi baru yang Anda gunakan, sehingga tidak akan mengkompilasi lebih dari yang Anda butuhkan namun akan memberi Anda kecepatan C / Fortran. Pada akhirnya ini mirip dengan menggunakan template C ++, tetapi ini adalah bahasa tingkat yang lebih tinggi dengan banyak alat untuk membuatnya lebih mudah (cukup mudah bahwa ini akan menjadi masalah pekerjaan rumah yang bagus untuk seorang mahasiswa).

Bahasa lain yang bagus untuk ini adalah Lisp seperti Common Lisp. Mudah digunakan karena seperti Julia memberi Anda AST yang dikompilasi dengan banyak alat introspeksi yang dibangun, tetapi tidak seperti Julia itu tidak akan secara otomatis mengompilasinya (dalam sebagian besar distribusi).

Saya berada di kapal (Fortran) yang sama. Setelah saya memiliki elemen 1D, 2D, 3D dan 4D (saya melakukan geometri projektif), saya membuat operator yang sama untuk setiap jenis dan kemudian menulis logika saya dengan persamaan tingkat tinggi yang memperjelas apa yang sedang terjadi. Ini tidak lambat karena Anda mungkin berpikir untuk memiliki loop terpisah dari setiap operasi dan banyak salinan memori. Saya membiarkan kompiler / prosesor melakukan optimasi.

Sebagai contoh

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

Untuk digunakan dalam persamaan seperti

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

di mana edan rdan gdapat memiliki dimensi apa pun yang masuk akal secara matematis.