オブジェクト指向Fortranで関数ポインタを使う

概要

Fortranの関数ポインタ（手続ポインタ）を用いてコールバック処理を実装します。これを用いると数値計算の初期値代入などを簡単に書けるようになります。

導入

古典的なFORTRANでは手続き型プログラミングのパラダイムで書かれる一方で、Fortran 2003より新しい規格では、オブジェクト指向プログラミング（object oriented programming, OOP）と関数ポインタ（より正確には手続ポインタ（procedure pointer）と呼ばれる）の機能が導入されています。

OOPで数値計算プログラムを記述することは、『Fortranによる実践オブジェクト指向プログラミング』（暗黙の型宣言 編、2019年第3版）で指摘されている通り、その変更容易性と人為的ミスの排除ができる点のメリットが大きいです。

FortranでのOOPは、派生型の宣言のcontains以下にサブルーチンや関数を宣言することで、それらの手続をその派生型に所属させることによって実現されます（これを型束縛手続（type-bound procedure）と言います ）。型束縛手続には、通常の手続と同様、手続ポインタを引数に取ることができます。

Fortranでの手続ポインタは、関数についての明示的なinterfaceブロックとprocedure文によって定義されます。手続ポインタを使用する方法は通常の関数呼び出しと同様で、別モジュールの手続に、メインプログラムで定義した関数などを渡すことができます。

これを利用すると、いわゆるコールバック処理を記述することができ、初期値代入や境界条件の処理を簡潔に記述することができます。そのため手続ポインタは、数値計算プログラムをOOPで書く際の変更容易性を補完することができるでしょう。

本稿では、手続ポインタを初期値代入に使う方法を述べた上で、OOPで書かれたモジュールに手続ポインタを渡す／受け取る方法について説明します。

初期値代入に手続ポインタを渡す

例えば、単純なsin関数を初期値に与えるとしましょう。

この初期値を与えるためには、次節で説明するOOPで実装されたモジュールを使用すると、以下の文を実行することで実現できます。

fptr => dsin            ! ポインタ変数fptrに倍精度型の正弦関数をアサインする
call u%allocate()       ! 配列を割り付ける型束縛手続を呼び出す
call u%initialize(fptr) ! 初期化のための型束縛手続に手続ポインタを渡して呼び出す

call data_output()      ! データ出力するサブルーチンを呼び出す
call u%deallocate()     ! 配列の割り付け解除をする型束縛手続を呼び出す

これをcos関数にしたい場合は、次のようにポインタが参照する関数を切り替えるだけです。

fptr => dcos            ! dcosは倍精度型の余弦関数
call u%initialize(fptr)

ユーザー定義の関数も渡すことができます。

   ...
   fptr => cliff          ! ユーザ定義の内部関数cliffを参照する
   call u%allocate()
   call u%initialize(fptr)
   call data_output()
   
contains

   ! 内部関数cliff
   function cliff(x)
      implicit none
      real(real64) :: cliff
      real(real64), intent(in) :: x
      real(real64), parameter :: pi_half = acos(-1d0)/2d0 
      
      ! 階段関数を返す
      if (x <= pi_half ) then
         cliff = 1d0   ! π/2以下ならば1.0
      else
         cliff = 0d0   ! それ以外は0.0
      end if
   
   end function cliff
   ...

これをプロットすると、次のようになります。

このような例から、手続ポインタを渡して配列の初期化を行う方法は、シンプルにコードを書くことができ、なおかつ変更が容易であることが分かります。次のセクションではこれらを実現するオブジェクト指向で書かれたモジュールについてを解説します。

モジュールをOOPで書く

本稿で使う例では、簡単のため、以下の3つのコードから構成される場合を考えます。

• 配列に対する初期化操作を行うkernelモジュール
• 派生型Variableを定義するclass_variableモジュール
• mainプログラム

kernelモジュール：kernel.f90

kernelモジュールは、必要な定数の定義と、コールバック処理を含みます。モジュール手続としてサブルーチンinitializeを実装しています。この手続きは、引数に配列と手続ポインタを受け取り、手続ポインタの関数をコールして結果をその配列に代入します。

module kernel
   use :: iso_fortran_env, only: real64
   implicit none
   private
   public :: initialize
   
   ! 定数を定義する
   integer, parameter, public      :: nx = 2**8+1       ! 離散点の数
   real(real64), parameter, public :: PI = acos(-1d0)   ! 円周率
   real(real64), parameter, public :: L = 2*PI          ! 計算領域の長さ、範囲[0,2π]を仮定   
   real(real64), parameter, public :: dx = L/dble(nx-1) ! 離散点の間隔

contains
   
   ! 関数ポインタを受け取ったときの初期化関数のコールバック処理を実装する。
   subroutine initialize(array, fp)
      use, intrinsic :: iso_fortran_env, only: real64
      implicit none
      real(real64), intent(out) :: array(:)
   
      ! ポインタを宣言する。
      pointer :: fp
      !手続ポインタfpのインターフェースブロックを記述する。
      interface
         function fp(x) result(res)
            use, intrinsic :: iso_fortran_env, only: real64
            real(real64), intent(in) :: x
            real(real64) :: res
         end function fp
      end interface

      integer :: i
      real(real64) :: x
      !--------ここまで宣言部-----------!

      do i = 1, nx
         ! x座標を計算する。
         x = dble(i-1)*dx
         
         ! 手続きポインタfpに引数xを渡して、結果をarray(i)に代入する。
         array(i) = fp(x)
      end do
   end subroutine initialize

end module kernel

モジュールの宣言部にinterfaceブロックを記述した場合は、内部手続きで以下のようなprocedure文で

procedure(hoge), pointer :: fp

と手続ポインタを宣言することができます。この場合は、interfaceブロック内で定義される名前をprocedure文の型パラメーター（この例ではhoge）に渡します。この名前は、変数名などと衝突するものを使用することはできません。

派生型の定義：class_Variable.f90

派生型を定義するモジュールclass_Variableです。オブジェクト指向プログラミングの概念を取り入れて実装された派生型Variableは、割り付け可能な成分value配列と、allocate、initialize、deallocateの型束縛手続を含みます。

3つの型束縛手続はvalueに対する操作を行います。

• allocateサブルーチンは、成分valueが割り付けられていない場合、これを割り付けます。
• initializeサブルーチンは、引数に手続ポインタを受け取り、自身のオブジェクトのvalue成分と手続ポインタをkernelモジュールのinitializeサブルーチンに渡します。
• deallocateサブルーチンは、成分valueの割り付けを解除します。

module class_Variable
   use, intrinsic :: iso_fortran_env
   private
   public :: Variable
   
   type Variable
      real(real64), allocatable, public :: value(:)
   contains
      procedure, public, pass :: allocate
      procedure, public, pass :: initialize
      procedure, public, pass :: deallocate
   end type Variable
 
contains
   
   ! 派生型selfの成分value配列を割り付ける。
   subroutine allocate(self)
      use :: kernel, only:Nx
      implicit none
      class(Variable), intent(inout) :: self
      
      if (.not. allocated(self%value)) then
         allocate(self%value(Nx), source=0d0)
      end if
   end subroutine allocate


   ! kernelの手続initializeに、手続ポインタfpを渡して初期化を実行するラッパー手続
   subroutine initialize (self, fp)
      use, intrinsic :: iso_fortran_env, only: real64
      use :: kernel, only: kernel_init => initialize
      implicit none
      class(Variable), intent(inout) :: self

      pointer :: fp
      interface
         function fp(x) result(res)
            use, intrinsic :: iso_fortran_env, only: real64
            real(real64), intent(in) :: x
            real(real64) :: res
         end function fp
      end interface

      call kernel_init(self%value, fp)
   end subroutine initialize


   ! 派生型selfの成分valueの割り付けを解除する。
   subroutine deallocate(self)
      implicit none
      class(Variable), intent(inout) :: self

      if (allocated(self%value)) then
          deallocate(self%value)
      end if
   end subroutine deallocate
   
end module class_Variable

手続ポインタを受け取るinitializeサブルーチンの宣言部には、その手続ポインタについてのinterfaceブロックを記述する必要があります。前のセクションと同様に、モジュールの宣言部にinterfaceブロックを置いて、procedure文で手続ポインタを宣言することができます。

メインプログラム：main.f90

以下はメインプログラムのコードです。use文でkernelとclass_Variableのモジュールを読み込んでいます。

実行部では、最初に出力用のx座標軸を表すVariable型変数xを作成した後、手続ポインタfptrに内部関数cliffを参照させています。次に型束縛手続u%allocate()とu%initialize(fptr)を呼び出し、手続ポインタを使用してvalue成分に初期値を代入しています。

最後にサブルーチンdate_outputを呼び出してデータ出力したあとに、型束縛手続deallcoateでu, x各変数の成分の割り付けを解除しています。

program main
   use :: iso_fortran_env, only:real64
   use :: kernel, only: nx, dx
   use :: class_Variable, only: Variable
   implicit none

   ! 派生型の宣言
   type(Variable) :: u, x
   
   ! 手続ポインタのインタフェースを定義する
   interface
      function fp(x)
         use, intrinsic :: iso_fortran_env, only: real64
         real(real64), intent(in) :: x
         real(real64) :: fp
      end function fp
   end interface
   ! 手続ポインタ変数fptrを宣言する
   procedure(fp), pointer ::  fptr => null()

   integer :: i
   
   ! x軸の座標値を作る。
   call x%allocate()
   do i = 1, nx
      x%value(i) = 0d0 + dx*dble(i-1)
   end do

   call u%allocate()

   ! 手続ポインタfptrに内部関数cliffを関連付ける。
   fptr => cliff
!  fptr => dsin
!  fptr => dcos
   
   ! Variable派生型変数uの型束縛手続initializeを、手続ポインタfptrを引数に渡して呼び出す。
   call u%initialize(fptr)
  
   call data_output(u)
   
   call u%deallocate()
   call x%deallocate()

contains

   ! 階段関数
   pure function cliff(x)
      implicit none
      real(real64) :: cliff
      real(real64), intent(in) :: x
      
      if (x <= 1.57d0) then
         cliff = 1d0
      else
         cliff = 0d0
      end if
   
   end function cliff

	! データ出力用のサブルーチン
   subroutine data_output(v)
      use :: kernel, only: nx
      implicit none
      type(Variable), intent(in) :: v
      
      integer, save :: uni = 10
      integer, save :: count = 1
      character(7), save :: filename = 'out.dat'
      logical :: isOpened

      inquire(file=filename, opened=isOpened)

      if (.not. isOpened) then
      	open(uni, file=filename, form='formatted', action='write', status='replace')
      end if
      
		write(uni, '(a, i0)') '> -Z', count  ! header for GMT6

      do i = 1, nx
         write(uni, '(e10.3, 1x, e10.3)') x%value(i), v%value(i)
      end do
      
   end subroutine data_output

end program main

このプログラムを実行すると、データファイルout.datが出力され、上のような画像をプロットすることができます。

以上のプログラムは次のようなシーケンス図で表現されるでしょう。

まとめ

オブジェクト指向プログラミングの方法を用いてFortranプログラムを構築し、配列の初期化に手続ポインタを渡して実行する方法について述べました。この方法を使うと、メインプログラムとモジュールの役割分担を明確にすることができ、コードの再利用性を高めることができます。

OOPのFortranの情報はネットでも少ないですが、手続ポインタの情報は極少ないです（Modern Fortran Explained 2018でも6ページしか分量がない）。本稿では、実際にどのように手続ポインタを使うことができるのかを述べ、その使い方を学ぶための情報を増やすことができたと思います。

補遺

• データのプロットにはGeneric Mapping Tools v6.4.0を使用しました。

• ifortでのポインタを使用したプログラムの実行は結構遅いので、問題がなければ最適化オプション-O2をつけるのをおすすめします。

テスト環境

OSはGentoo Linux、コンパイラのバージョンは以下の通りです。

• % gfortran --version
  GNU Fortran (Gentoo 12.3.1_p20230526 p2) 12.3.1 20230526
  Copyright (C) 2022 Free Software Foundation, Inc.
  This is free software; see the source for copying conditions.  There is NO
  warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  

• % ifort --version
  ifort (IFORT) 2021.9.0 20230302
  Copyright (C) 1985-2023 Intel Corporation.  All rights reserved.