Operações sobre Elementos em Matrizes

Há três abordagens para realizar operações sobre elementos em matrizes quando usando subrotinas e funções:

• Procedimentos do tipo elemental

• matrizes explicit-shape

• implementando a operação para vetores e escrevendo um empacotador simples de subrotinas (que use reshape internamente) para cada forma da matriz

Na primeira abordagem, usa-se a palavra-chave elemental para criar uma função como essa:

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

Todos os argumentos (dentro e fora) devem ser escalares. Você pode então usar essa função com matrizes de qualquer tipo (compatível), por exemplo:

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)

A saída será:

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

No código acima, tipicamente n é um parâmetro e x é uma matriz de formato arbitrário, mas como pode ver, Fortran não se importa desde que a operação final faça sentido (se um argumento é uma matriz, então os outros argumentos devem ser matrizes do mesmo formato ou escalares). Se não for assim, você receberá um erro do compilador.

A palavra-chave elemental implica na palavra-chave pure, então o procedimento deve ser puro. Isso resulta que elemental procedures só podem usar procedimentos pure e não há efeitos colaterais.

Se o algoritmo de procedimento elementar puder ser feito mais rápido usando operações de matrizes internamente, ou se por alguma razão os argumentos deverão ser matrizes de formatos incompatíveis, então deve-se usar uma das duas outras abordagens. Pode-se fazer o nroot operar em um vetor ou escrever um simples empacotador para os outros formatos das matrizes, p.ex.:

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

E use da seguinte forma:

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]))

O que será mostrado na tela como:

3.0000000000000000
1.0000000000000000        2.0000000000000000        3.0000000000000000        3.1622776601683795
1.0000000000000000        2.0000000000000000        3.0000000000000000        3.1622776601683795

Ou pode-se usar matrizes explicite-shape como apresentado a seguir:

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

Use como se segue:

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]))

A saída é a mesma que a anterior:

3.0000000000000000
1.0000000000000000        2.0000000000000000        3.0000000000000000        3.1622776601683795
1.0000000000000000        2.0000000000000000        3.0000000000000000        3.1622776601683795