数组上的逐元素操作#

在使用子程序和函数时,有三种方法可以对数组执行逐元素操作

  • elemental 过程

  • 显式形状数组

  • 为向量实现操作并编写简单的包装子程序(在内部使用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

所有参数(输入和输出)必须是标量。然后可以使用此函数与任何(兼容)形状的数组一起使用,例如

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

或者可以使用显式形状数组,如下所示

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