自力でFortranからFDPSを使ってみる

# 概要

粒子シミュレーション開発用フレームワークである[FDPS](https://github.com/FDPS/FDPS)はFortranからも利用できますが，本記事では開発元が提供するインタフェースを使わず，自力でインタフェースを書いてFDPSをつかってみました．その方法は開発元のインタフェースと本質的には同じで，本記事はただの車輪の再発明ですが，Fortran 2003で導入されたC-Fortran相互利用の学習にはよいと思います．

# FDPS

FDPSとは，Framework for Developing Particle Simulatorの略で，粒子系シミュレーションに分類される方法（N体やSPH，渦法など）の高性能化を支援するフレームワークです．理化学研究所計算科学研究機構粒子系シミュレータ研究チームが開発しています．C++のヘッダーライブラリとして開発されていますが，FDPSバージョン3からはFortranからも利用できるようになりました．

粒子系のシミュレーション方法の多くは，基本的なアルゴリズムは単純であるものの，粒子数の増加に伴って計算量が急激に増加します．遠方の粒子の影響をまとめるなど様々な高速化の工夫があるのですが，実装するのはかなり手間がかかります．並列化（プロセス並列）も，差分法など格子系の方法と比較すると結構大変です．FDPSは，それらを支援してくれるかなり素晴らしいフレームワークです．

使い方については[チュートリアル](https://github.com/FDPS/FDPS/blob/master/doc/doc_tutorial_cpp_ja.pdf?raw=true)や[他の方の記事](https://qiita.com/kaityo256/items/c4955adaad718fd399bc)を参考にしてください．

## FDPSをFortranから利用するには

詳しい説明は，これも[チュートリアル](https://github.com/FDPS/FDPS/blob/master/doc/doc_tutorial_ftn_ja.pdf?raw=true)に譲りますが，簡単にいうと，

- 粒子の型およびそれらの相互作用を定義したuser_defined.f90
- 全体の処理の流れを記述したf_main.f90

の二つを作成し，Pythonスクリプトを利用してC-Fortranのインタフェースを作成します．

C言語の`main`関数からf_main.f90に記述したサブルーチンを呼出し，当該サブルーチンではC-Fortranのインタフェースを介してFDPSの関数を呼び出すという流れで実行されます．

# 自作のインタフェース

FDPSの動作を理解するために，[この記事](https://qiita.com/kaityo256/items/c4955adaad718fd399bc)に沿って`Initialize`と`Finalize`をしようと試みたのですが，Pythonスクリプトでインタフェースを生成するには粒子情報を書く必要があったため，開発元が提供するインタフェースを利用せず，自力で書くことにしました．

## C側

Fortranは名前空間に相当する機能がないので，FDPSの関数がFortranから見えるよう，まずC側で対応する必要があります．そこで，`extern "C"`でラップしただけの関数を定義しました．

```C:fdps_c_if.cpp
#include<particle_simulator.hpp>

extern "C" {
    void PS_Initialize(int argc, char **argv){
        PS::Initialize(argc, argv);
    }
    void PS_Finalize(void){
        PS::Finalize();
    }
}
```

## Fortran側

FortranからCの関数を利用するには，interfaceを記述します．Fortran 2003では，`bind(c, name="C側の関数名")`でC言語の関数名を指定します．昔は関数名にアンダースコアを付けるなどの対処が必要だったようですが，Modern Fortranではその必要はありません．

```Fortran:fdps_f_if.f90
module FDPS_Fortran_Interface
    use, intrinsic :: iso_c_binding
    implicit none
    public

    interface
        subroutine PS_initialize(argc, argv) bind(c, name="PS_Initialize")
            use, intrinsic :: iso_c_binding
            implicit none
            integer(c_int),value :: argc
            type(c_ptr),intent(in) :: argv 
        end subroutine PS_initialize

        subroutine PS_finalize() bind(c, name="PS_Finalize")
            use,intrinsic :: iso_c_binding
            implicit none
        end subroutine PS_finalize
    end interface

end module FDPS_Fortran_Interface
```

次に，Cの`main`関数から呼ぶためのサブルーチン`f_main`を記述します．`PS_Initialize`に`argc`と`argv`を渡すために，`f_main`にも引数として`argc`と`argv`を設けます．

`argc`は値渡しである必要があるので，`value`属性が必要です．`argv`の型は`char **`ですが，Fortranが標準でポインタ渡しであることを考慮すると，ダブルポインタではなくポインタを渡せばよいことになります．FortranでC言語のポインタを利用するには，派生型`type(c_ptr)`を利用します．なお，`type(c_ptr)`は`integer(c_intptr_t)`をプライベートな成分として持つ派生型でしかなく，アドレスの先にある数値が何型かを調べることはできません．

サブルーチン名の後ろにある`bind(c, name = "f_main")`は，上のモジュールとは逆に，C側に見える名前を設定します．

```Fortran:f_main.f90
subroutine f_main(argc, argv) bind(c, name = "f_main")
    use,intrinsic :: iso_c_binding
    use :: FDPS_Fortran_Interface
    implicit none
    integer(c_int),value :: argc
    type(c_ptr),intent(in) :: argv

    call PS_Initialize(argc, argv)
    call PS_Finalize()
end subroutine f_main
```

## Cのmain関数

最後に，Cの`main`関数を書いて，`f_main`を呼びます．

```C:main.cpp
extern "C" void f_main(int argc, char **argv);

int main(int argc, char **argv){
    f_main(argc, argv);

    return 0;
}
```

# コンパイルと実行

作成したプログラムのコンパイルと実行は，以下の環境で行いました．

- Ubuntu 16.04 64bit
- g++ 5.4.0
- gfortran 5.4.0
- FDPS 4.1b （~/FDPS以下に展開）

コンパイルは，とりあえず愚直に以下のように行いました．

```console
g++ -c fdps_c_if.cpp -I ~/FDPS/src
gfortran -c fdps_f_if.f90
gfortran -c f_main.f90
g++ -c main.cpp
g++ main.o f_main.o fdps_f_if.o fdps_c_if.o -lgfortran
```

作成された実行ファイルを実行すると，無事にFPDSが実行されていることが確認できました．

```
     //==================================\\
     ||                                  ||
     || ::::::: ::::::. ::::::. .::::::. ||
     || ::      ::    : ::    : ::       ||
     || ::::::  ::    : ::::::'  `:::::. ||
     || ::      ::::::' ::      `......' ||
     ||     Framework for Developing     ||
     ||        Particle Simulator        ||
     ||     Version 4.1b (2018/08)       ||
     \\==================================//

       Home   : https://github.com/fdps/fdps 
       E-mail : fdps-support@mail.jmlab.jp
       Licence: MIT (see, https://github.com/FDPS/FDPS/blob/master/LICENSE)
       Note   : Please cite Iwasawa et al. (2016, Publications of the Astronomical Society of Japan, 68, 54)

       Copyright (C) 2015 
         Masaki Iwasawa, Ataru Tanikawa, Natsuki Hosono,
         Keigo Nitadori, Takayuki Muranushi, Daisuke Namekata,
         Kentaro Nomura, Junichiro Makino and many others
******** FDPS has successfully begun. ********
******** FDPS has successfully finished. ********
```

# まとめ

FDPSの`Initialize`と`Finalize`を実行したいがためだけに，C-Fortranのインタフェースを作成しました．これまで，FortranからCを呼ぶことはあっても，CからFortranを呼ぶことがなかったのでよい経験になりました．開発元のインタフェースを全て置き換える気はありません（できません）が，車輪の再発明といっても，車輪を発明する過程は勉強になります．これからMPI実行なども試していきたいと思います．

書くだけ書いてから，Cの`main`を書かなくても，Fortranから`Initialize`と`Finalize`を実行するだけでよかったのではないかと気付きました．

## メモ
Fortranの`bind(c, name="")`は，interfaceで用いる場合にはC言語の関数名を指定し，サブルーチンの実装で用いる場合には，C言語側に見せるサブルーチン名前を指定する．