Элементные операции с массивами

Существует три подхода к выполнению элементных операций с массивами при использовании процедур и функций:

• elemental (элементные) подпрограммы

• explicit-shape массивы (явной формы)

• реализация операции для векторов и написание простых процедур-обёрток (которые используют внутри операцию reshape) для каждой формы массива

В первом подходе используется ключевое слово elemental для создания функции, подобной этой:

real(dp) elemental function nroot(n, x) result(y)
  integer, intent(in) :: n
  real(dp), intent(in) :: x
  y = x**(1._dp / n)
end function

Все аргументы (in и out) должны быть скалярами. Вы можете использовать эту функцию с массивами любой (совместимой) формы, например:

print *, nroot(2, 9._dp)
print *, nroot(2, [1._dp, 4._dp, 9._dp, 10._dp])
print *, nroot(2, reshape([1._dp, 4._dp, 9._dp, 10._dp], [2, 2]))
print *, nroot([2, 3, 4, 5], [1._dp, 4._dp, 9._dp, 10._dp])
print *, nroot([2, 3, 4, 5], 4._dp)

Вывод будет следующим:

3.0000000000000000
1.0000000000000000        2.0000000000000000        3.0000000000000000        3.1622776601683795
1.0000000000000000        2.0000000000000000        3.0000000000000000        3.1622776601683795
1.0000000000000000        1.5874010519681994        1.7320508075688772        1.5848931924611136
2.0000000000000000        1.5874010519681994        1.4142135623730951        1.3195079107728942

В приведённом выше примере обычно n является параметром, а x – массивом произвольной формы, но, как видите, языку Fortran всё равно, пока конечная операция имеет (если один аргумент является массивом, то остальные аргументы должны быть либо массивами той же формы, либо скалярами). Если это не так, то вы получите сообщение компилятора об ошибке.

Ключевое слово elemental подразумевает ключевое слово pure, поэтому процедура должна быть чистой. Это приводит к тому, что elemental procedures (элементные подпрограммы) могут использовать только pure (чистые) подпрограммы и не имеют побочных эффектов.

Если алгоритм элементной подпрограммы можно сделать быстрее, используя внутри операции с массивами или если по каким-то причинам аргументы должны быть массивами несовместимых форм, то следует использовать два других подхода. Можно заставить nroot работать с вектором и написать простую обёртку для других форм массивов, например:

function nroot(n, x) result(y)
  integer, intent(in) :: n
  real(dp), intent(in) :: x(:)
  real(dp) :: y(size(x))
  y = x**(1._dp / n)
end function

function nroot_0d(n, x) result(y)
  integer, intent(in) :: n
  real(dp), intent(in) :: x
  real(dp) :: y
  real(dp) :: tmp(1)
  tmp = nroot(n, [x])
  y = tmp(1)
end function

function nroot_2d(n, x) result(y)
  integer, intent(in) :: n
  real(dp), intent(in) :: x(:, :)
  real(dp) :: y(size(x, 1), size(x, 2))
  y = reshape(nroot(n, reshape(x, [size(x)])), [size(x, 1), size(x, 2)])
end function

И использовать следующим образом:

print *, nroot_0d(2, 9._dp)
print *, nroot(2, [1._dp, 4._dp, 9._dp, 10._dp])
print *, nroot_2d(2, reshape([1._dp, 4._dp, 9._dp, 10._dp], [2, 2]))

Вывод будет следующим:

3.0000000000000000
1.0000000000000000        2.0000000000000000        3.0000000000000000        3.1622776601683795
1.0000000000000000        2.0000000000000000        3.0000000000000000        3.1622776601683795

Или можно использовать explicit-shape массивы (явной формы):

function nroot(n, k, x) result(y)
  integer, intent(in) :: n, k
  real(dp), intent(in) :: x(k)
  real(dp) :: y(k)
  y = x**(1._dp / n)
end function

Применяя следующим образом:

print *, nroot(2, 1, [9._dp])
print *, nroot(2, 4, [1._dp, 4._dp, 9._dp, 10._dp])
print *, nroot(2, 4, reshape([1._dp, 4._dp, 9._dp, 10._dp], [2, 2]))

Вывод будет таким же, как в прошлый раз:

3.0000000000000000
1.0000000000000000        2.0000000000000000        3.0000000000000000        3.1622776601683795
1.0000000000000000        2.0000000000000000        3.0000000000000000        3.1622776601683795