Как написать код, не зависящий от размеров?

Я часто нахожу себя пишущим очень похожий код для одной, двух и трехмерных версий данной операции / алгоритма. Поддержание всех этих версий может стать утомительным. Простая генерация кода работает довольно хорошо, но кажется, что должен быть лучший способ.

Есть ли относительно простой способ написать операцию один раз и сделать ее обобщенной для более высоких или более низких измерений?

Один конкретный пример: предположим, мне нужно вычислить градиент поля скорости в спектральном пространстве. В трех измерениях петли Фортрана будут выглядеть примерно так:

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

где ddxмассив правильно определен. (Можно также сделать это с умножением матрицы.) Код для двумерного потока почти точно такой же, за исключением: третье измерение отбрасывается из циклов, индексов и количества компонентов. Есть ли лучший способ выразить это?

Другой пример: предположим, у меня есть скорости жидкости, определенные точечно на трехмерной сетке. Чтобы интерполировать скорость в произвольное место (т. Е. Не соответствующее точкам сетки), можно использовать одномерный алгоритм Невилла последовательно по всем трем измерениям (т. Е. Уменьшению размеров). Есть ли простой способ сделать размерную редукцию при одномерной реализации простого алгоритма?

Вы посмотрите как дела.II ( http://www.dealii.org/ ) это делает - там независимость измерений лежит в самом сердце библиотеки и моделируется как аргумент шаблона для большинства типов данных. См., Например, не зависящий от измерения решатель Лапласа в учебной программе шага 4:

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

Смотрите также

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

Вопрос подчеркивает, что большинство «простых» языков программирования (по крайней мере, C, Fortran) не позволяют вам делать это чисто. Дополнительным ограничением является то, что вы хотите удобство записи и хорошую производительность.

Поэтому вместо написания кода, специфичного для измерения, подумайте о написании кода, который генерирует специфичный для измерения. Этот генератор не зависит от размеров, даже если вычислительный код - нет. Другими словами, вы добавляете слой рассуждений между вашей нотацией и кодом, выражающим вычисления. Шаблоны C ++ означают одно и то же: с другой стороны, они встроены прямо в язык. Оборотная сторона, они несколько громоздки для написания. Это сводит вопрос к тому, как практически реализовать генератор кода.

OpenCL позволяет вам делать генерацию кода во время выполнения довольно чисто. Это также обеспечивает очень четкое разделение между «внешней управляющей программой» и «внутренними циклами / ядрами». Внешняя генерирующая программа гораздо менее ограничена в производительности и поэтому может быть написана на удобном языке, таком как Python. Это моя надежда на то, как PyOpenCL будет использоваться - извините за обновленный бесстыдный плагин.

Это может быть выполнено на любом языке с помощью следующего грубого умственного прототипа:

1. Создайте список экстентов каждого измерения (что-то вроде shape () в MATLAB, я думаю)
2. Создайте список вашего текущего местоположения в каждом измерении.
3. Напишите цикл для каждого измерения, содержащий цикл для изменения размера чьего-либо в зависимости от внешнего цикла.

Оттуда, это вопрос борьбы с синтаксисом вашего определенного языка, чтобы сохранить ваш код nd-совместимым.

Написав n-мерный решатель гидродинамики , я обнаружил, что полезно иметь язык, который поддерживает распаковку списка, подобного объекту, в качестве аргументов функции. Т.е. a = (1,2,3) f (a *) -> f (1,2,3). Дополнительно продвинутые итераторы (такие как ndenumerate в numpy) делают код на порядок чище.

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

Четкие ответы, если вы хотите сохранить скорость Fortran, - это использовать язык с надлежащей генерацией кода, такой как Julia или C ++. Шаблоны C ++ уже упоминались, поэтому я упомяну инструменты Джулии здесь. Сгенерированные функции Джулии позволяют вам использовать ее метапрограммирование для создания функций по запросу через информацию о типе. Так что, по сути, вы можете сделать здесь

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

а затем вы используете Nпрограммно для сборки кода, который вы хотите выполнить, учитывая его Nразмерность. Затем можно легко построить картезианскую библиотеку Джулии или пакеты, такие как выражения Einsum.jl, дляN размерной функции.

Что хорошо в Джулии, так это то, что эта функция статически компилируется и оптимизируется для каждого используемого вами нового размерного массива, поэтому она не будет компилировать больше, чем вам нужно, но при этом вы получите скорость C / Fortran. В конце концов, это похоже на использование шаблонов C ++, но это язык более высокого уровня с множеством инструментов, облегчающих его (достаточно просто, чтобы это было хорошей домашней работой для старшекурсника).

Другой язык, который хорош для этого, - Лисп, как Common Lisp. Он прост в использовании, поскольку, как и Julia, он предоставляет вам скомпилированный AST с большим количеством встроенных инструментов для самоанализа, но в отличие от Julia он не будет автоматически компилировать его (в большинстве дистрибутивов).

Я в той же (Фортран) лодке. Когда у меня есть элементы 1D, 2D, 3D и 4D (я занимаюсь проективной геометрией), я создаю одинаковые операторы для каждого типа, а затем пишу свою логику с помощью уравнений высокого уровня, которые проясняют, что происходит. Это не так медленно, как вы могли бы подумать, чтобы иметь отдельные циклы каждой операции и много копий памяти. Я позволил компилятору / процессору выполнить оптимизацию.

Например

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

Быть использованным в уравнениях как

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

где eи rи gможет иметь любую размерность, которая имеет математический смысл.