微分方程式の数値解析 #7

RK4を使って一階連立微分方程式を解くためのモジュール。
サブルーチンset_funcでのを設定し、
関数nextで時刻tの値から時刻t+hの値を算出する。
FORTRAN 2003以降でないとコンパイルできないと思う。

module runge_kutta_4
  implicit none
  interface
    function func_interface(t, x)
      double precision :: func_interface
      double precision, intent(in) :: t
      double precision, dimension(:), intent(in) :: x
    end function func_interface
  end interface
  type func_type
    procedure (func_interface), pointer, nopass :: f
  end type func_type
  type rk4
    type (func_type), allocatable, dimension(:) :: f
  end type rk4
  type (rk4) f
contains
  subroutine set_func(fn)
    type (func_type), dimension(:) :: fn
    allocate(f%f(size(fn)))
    f%f = fn
  end subroutine set_func
  function next(t, x, h)
    double precision, intent(in) :: t, h
    double precision, dimension(:), intent(in) :: x
    double precision, dimension(size(x)) :: next
    double precision, dimension(size(x), 4) :: k
    integer :: i
    do i = 1, size(x)
      k(i, 1) = f%f(i)%f(t, x)
    end do
    do i = 1, size(x)
      k(i, 2) = f%f(i)%f(t + h / 2, x + h / 2 * k(:, 1))
    end do
    do i = 1, size(x)
      k(i, 3) = f%f(i)%f(t + h / 2, x + h / 2 * k(:, 2))
    end do
    do i = 1, size(x)
      k(i, 4) = f%f(i)%f(t + h, x + h * k(:, 3))
    end do
    next = x + h / 6 * (k(:, 1) + 2 * k(:, 2) + 2 * k(:, 3) + k(:, 4))
  end function next
end module runge_kutta_4

本当はparametrized derived typeでtype rk4を定義できればもう少しスマートになると思うが、
手元の環境にあるgfortranのバージョンではパラメータを認識できず未サポートっぽいのでとりあえず上のように。

ところで、x + h / 2 * k(:, 1)等のdoループ内の配列演算はループ外に掃き出す最適化は行われるのだろうか？
もっとコンパイラがヒントを得られるようにしないと最適化の対象にはならなさそうな気が。