FortranでC言語の関数を呼び出し例

# iso_c_binding
`iso_c_binding` を知っていますか?
なんと「FortranからC 言語の関数を呼び出すこと」ができてしまいます.
(私は 「**自作の**C 言語の関数を呼び出せる」と思っていましたが, 「[`printf`や`strcmp`のようなC言語の関数も呼び出せる](https://qiita.com/ijknabla/items/bd7fe074239741c7c7d9)」ということに気づきました...)

# やりたいこと
- **C 言語呼び出しの例を提示したい!**
### そのうちやりたいこと
- Fortran ではできないこと, 面倒なことを C 言語を通してやりたい!

# 実行環境
```bash
$ lsb_release -a
No LSB modules are available.
Distributor ID:	Ubuntu
Description:	Ubuntu 20.04.6 LTS
Release:	20.04
Codename:	focal
$ gfortran --version
GNU Fortran (GCC) 9.2.0
Copyright (C) 2019 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.
```

# 引数が可変長な関数を呼び出そう!
## Fortran から `printf` を呼び出してみよう!
`interface` を使ってFortranの関数の`cprintf_int`として宣言してみます.

```test_bind_printf.f90
module cprintf_m
  use, intrinsic :: iso_c_binding
  implicit none
  interface
    subroutine cprintf_int(s, n) bind(c, name = "printf") ! ②
      import c_char, c_int
      character(c_char), intent(in) :: s
      integer(c_int), value :: n
    end subroutine
  end interface
end module

program test_bind_printf
  use cprintf_m
  implicit none
  call cprintf_int("v = %d\n" // c_null_char, 42_c_int) ! ①
  write(*, *)
end program test_bind_printf

```
```sh
$ gfortran test_bind_printf.f90 && ./a.out
v = 42\n
```

① ちゃんと `v = 42\n` が出力されていますね... `\n`?
  + Fortranでエスケープシーケンスを使う場合は `-fbackslash` オプションでコンパイルしましょう.
  + もしくは, 改行する場合は `"v = %d" // achar(10) // c_null_char` で改行コード `achar(10)`を文字列結合 `//` しましょう.
  + `c_null_char` は `iso_c_binding` モジュールで宣言されている定数です.
    - C言語では, 文字列終端をヌル文字'\0'で表しますが, Fortran でC文字列を再現するなら, 文字列の後ろに `c_null_char` を `//` で結合させます. (ないと危険です.)

② c_biding関連
  + **iso_c_binding**
    - `c_char`, `c_int`, `c_null_char` が入っているモジュールです.
  + **subroutine cprintf_int(s, n) bind(c, name = "printf")**
    - `bind(c, name = "printf")` でC言語のprintf関数を参照するっぽいです.
    - C言語の `printf` 関数は可変長引数を取りますが, Fortranでは可変長引数はおそらく無理なので2変数にしています. 
    - 関数オーバーロードを使えば, 擬似的にn変数関数にできます(後述).
    - `interface` の `subroutine` の中では, 外で宣言された変数 `c_int`, `c_char` は `import` などをしないと見えません.
    - `c_int`, `c_char` はそれぞれC言語の `int` と `char *` に対応する `kind` です.
    - printf の第一引数は `const char *` なので, `intent(in)` でアドレス渡しをします.
      + `intent(inout)` を使っても(変更されないので)良いですが, `intent(in)`の方が文字列リテラルを渡せるので便利です.
    - 第二引数以降は値で渡すかアドレスで渡すかは, [第一引数のフォーマット](https://en.cppreference.com/w/c/io/fprintf)に依存します.
      + 今回は `int` 決め打ちなので `value` 属性で渡します.

`interface` でオーバーロード
```test_bind_printf_overload.f90
module cprintf_m
  use, intrinsic :: iso_c_binding
  implicit none
  interface cprintf
    subroutine cprintf_int(s, n) bind(c, name = "printf")
      import c_char, c_int
      character(c_char), intent(in) :: s
      integer(c_int), value :: n
    end subroutine
    subroutine cprintf_int_double(s, n, v) bind(c, name = "printf")
      import c_char, c_int, c_double
      character(c_char), intent(in) :: s
      integer(c_int), value :: n
      real(c_double), value :: v
    end subroutine
    subroutine cprintf_string(s, s2) bind(c, name = "printf")
      import c_char
      character(c_char), intent(in) :: s
      character(c_char), intent(in) :: s2
    end subroutine
  end interface cprintf
end module

program test_bind_printf
  use, intrinsic :: iso_c_binding
  use cprintf_m
  implicit none
  character, parameter :: return_char = achar(10)
  call cprintf("v = %d\n" // return_char // c_null_char, 42_c_int)
  call cprintf("s = %s\n" // return_char // c_null_char, "hello" // c_null_char)
  call cprintf("a = %05d, b = %.17f\n" // return_char // c_null_char, -10_c_int, 4.2d0)
end program test_bind_printf
```
```sh
$ gfortran test_bind_printf_overload.f90 && ./a.out
v = 42\n
s = hello\n
a = -0010, b = 4.20000000000000018\n
```

- 引数を3つにしても大丈夫です.
  + `interface` 内にひたすら色々書いていけば, 擬似的な printf が作れそうですね...
- printf 使うくらいなら, Fortran の write 文で良さそうですね.
  + むしろ, ヌル文字に警戒しないといけないため, やめた方が良いです.


# 引数がすごいやつを呼んでみたい!
## [qsort](https://en.cppreference.com/w/c/algorithm/qsort)を呼んでみたい!
qsort の引数には, 配列へのポインタと関数へのポインタが含まれるので, これを呼べたらなんでもできそうです.

https://en.cppreference.com/w/c/algorithm/qsort より, 引数:
> void qsort( void* ptr, size_t count, size_t size, int (*comp)(const void*, const void*) );

:::note warn
ポインタ関連の低レベルな場所に片足を突っ込んでいるので, 面倒なことが起きそうです.
素直に実装するか, `fortran-stdlib` のソート関数を使いましょう.
:::

<details open><summary> インターフェース部分 ソースコード(1/3) </summary>

```test_qsort.f90
module cqsort_m
  use, intrinsic :: iso_fortran_env
  use, intrinsic :: iso_c_binding
  implicit none
  interface
     subroutine qsort(ptr, cnt, siz, fn) bind(c, name = "qsort")
       import c_ptr, c_size_t, c_funptr
       type(c_ptr), value :: ptr
       integer(c_size_t), value :: cnt, siz
       type(c_funptr), value :: fn
     end subroutine qsort
  end interface
end module cqsort_m
```
- `type(c_ptr)` は `void *` です.
- `type(c_funptr)` は 関数ポインタです.
- `value` 属性で値渡しにします.
</details>
<details open><summary> ラッパールーチン ソースコード(2/3) </summary>

```test_qsort.f90
module qsort_int32_m
  use, intrinsic :: iso_fortran_env
  use, intrinsic :: iso_c_binding
  use cqsort_m
  interface fqsort
     module procedure :: fqsort_int32
  end interface fqsort
contains
  integer(c_int) function int_comp(a, b) result(res)
    integer(c_int32_t), intent(in) :: a, b
    res = int(a - b, c_int)
  end function int_comp

  subroutine fqsort_int32(arr)
    integer(int32), parameter :: bit_per_byte = 8
    integer(int32), intent(inout), target :: arr(:)
    integer(c_size_t) :: n, ss
    n = size(arr, kind = c_size_t)
    ss = storage_size(arr(1), kind = c_size_t) / bit_per_byte
    call qsort(c_loc(arr), n, ss, c_funloc(int_comp))
  end subroutine fqsort_int32
end module qsort_int32_m

module qsort_tuple2_m
  use, intrinsic :: iso_fortran_env
  use, intrinsic :: iso_c_binding
  use cqsort_m
  implicit none
  type :: tuple2
     integer(int32) :: x, y
  end type tuple2
  interface fqsort
     module procedure :: fqsort_tupel2
  end interface fqsort
contains
  integer(c_int) function tuple2_comp(a, b) result(res)
    type(tuple2), intent(in) :: a, b
    if (a%x < b%x) then
       res = -1_c_int
       return
    else if (a%x > b%x) then
       res = 1_c_int
       return
    end if
    !> a%x == b%x
    res = int(a%y - b%y, c_int)
  end function tuple2_comp

  subroutine fqsort_tupel2(arr)
    integer(int32), parameter :: bit_per_byte = 8
    type(tuple2), intent(inout), target :: arr(:)
    integer(c_size_t) :: n, ss
    n = size(arr, kind = c_size_t)
    ss = storage_size(arr(1), kind = c_size_t) / bit_per_byte
    call qsort(c_loc(arr), n, ss, c_funloc(tuple2_comp))
  end subroutine fqsort_tupel2
end module qsort_tuple2_m
```
- サブルーチン `fqsort_int32` 
  + 配列のサイズを `n`.
  + `storage_size` はビット数で返すので, `8` で割ってバイト数に変換しています.
  + `c_loc` は `type(c_ptr)` を返します. 詳しくは [nag](https://www.nag-j.co.jp/fortran/tips/tips_InteroperabilityWithC.html)
  + `c_funloc` は `type(c_funptr)` を返します.
- モジュール `qsort_tuple2_m`
  + ユーザ定義型 `tuple2` でも `int32` と同様に定義すればよいみたいです.
</details>
<details open><summary> main ソースコード(3/3) </summary>

```test_qsort.f90
program tets_fqsort
  use, intrinsic :: iso_fortran_env
  use cqsort_m
  use qsort_int32_m
  use qsort_tuple2_m
  implicit none
  test_int32: block
    integer(int32) :: arr(5) = [2, 5, 1, 4, 3]
    write(output_unit, '(a)') "before sort"
    write(output_unit, '(*(i0, 1x))') arr(:)
    call fqsort(arr)
    write(output_unit, '(a)') "after sort"
    write(output_unit, '(*(i0, 1x))') arr(:)
  end block test_int32
  write(error_unit, *)
  test_tuple2: block
    type(tuple2) :: arr(5) = &
         & [tuple2(100, 3), &
         &  tuple2(5, 1), &
         &  tuple2(4, 3), &
         &  tuple2(-1, -1), &
         & tuple2(5, 500)]
    integer(int32) :: i
    write(output_unit, '(a)') "before sort"
    do i = 1, size(arr)
       write(output_unit, '("(", i0, ", ", i0, ")")') arr(i)
    end do
    call fqsort(arr)
    write(output_unit, '(a)') "after sort"
    do i = 1, size(arr)
       write(output_unit, '("(", i0, ", ", i0, ")")') arr(i)
    end do
  end block test_tuple2
end program tets_fqsort
```
- `fqsort` を呼ぶだけです.
</details>

<details><summary> ソースコード全体 </summary>

```test_qsort.f90
module cqsort_m
  use, intrinsic :: iso_fortran_env
  use, intrinsic :: iso_c_binding
  implicit none
  interface
     subroutine qsort(ptr, cnt, siz, fn) bind(c, name = "qsort")
       import c_ptr, c_size_t, c_funptr
       type(c_ptr), value :: ptr
       integer(c_size_t), value :: cnt, siz
       type(c_funptr), value :: fn
     end subroutine qsort
  end interface
end module cqsort_m

module qsort_int32_m
  use, intrinsic :: iso_fortran_env
  use, intrinsic :: iso_c_binding
  use cqsort_m
  interface fqsort
     module procedure :: fqsort_int32
  end interface fqsort
contains
  integer(c_int) function int_comp(a, b) result(res)
    integer(c_int32_t), intent(in) :: a, b
    res = int(a - b, c_int)
  end function int_comp

  subroutine fqsort_int32(arr)
    integer(int32), parameter :: bit_per_byte = 8
    integer(int32), intent(inout), target :: arr(:)
    integer(c_size_t) :: n, ss
    n = size(arr, kind = c_size_t)
    ss = storage_size(arr(1), kind = c_size_t) / bit_per_byte
    call qsort(c_loc(arr), n, ss, c_funloc(int_comp))
  end subroutine fqsort_int32
end module qsort_int32_m

module qsort_tuple2_m
  use, intrinsic :: iso_fortran_env
  use, intrinsic :: iso_c_binding
  use cqsort_m
  implicit none
  type :: tuple2
     integer(int32) :: x, y
  end type tuple2
  interface fqsort
     module procedure :: fqsort_tupel2
  end interface fqsort
contains
  integer(c_int) function tuple2_comp(a, b) result(res)
    type(tuple2), intent(in) :: a, b
    if (a%x < b%x) then
       res = -1_c_int
       return
    else if (a%x > b%x) then
       res = 1_c_int
       return
    end if
    !> a%x == b%x
    res = int(a%y - b%y, c_int)
  end function tuple2_comp

  subroutine fqsort_tupel2(arr)
    integer(int32), parameter :: bit_per_byte = 8
    type(tuple2), intent(inout), target :: arr(:)
    integer(c_size_t) :: n, ss
    n = size(arr, kind = c_size_t)
    ss = storage_size(arr(1), kind = c_size_t) / bit_per_byte
    call qsort(c_loc(arr), n, ss, c_funloc(tuple2_comp))
  end subroutine fqsort_tupel2
end module qsort_tuple2_m

program tets_fqsort
  use, intrinsic :: iso_fortran_env
  use cqsort_m
  use qsort_int32_m
  use qsort_tuple2_m
  implicit none
  test_int32: block
    integer(int32) :: arr(5) = [2, 5, 1, 4, 3]
    write(output_unit, '(a)') "before sort"
    write(output_unit, '(*(i0, 1x))') arr(:)
    call fqsort(arr)
    write(output_unit, '(a)') "after sort"
    write(output_unit, '(*(i0, 1x))') arr(:)
  end block test_int32
  write(error_unit, *)
  test_tuple2: block
    type(tuple2) :: arr(5) = &
         & [tuple2(100, 3), &
         &  tuple2(5, 1), &
         &  tuple2(4, 3), &
         &  tuple2(-1, -1), &
         & tuple2(5, 500)]
    integer(int32) :: i
    write(output_unit, '(a)') "before sort"
    do i = 1, size(arr)
       write(output_unit, '("(", i0, ", ", i0, ")")') arr(i)
    end do
    call fqsort(arr)
    write(output_unit, '(a)') "after sort"
    do i = 1, size(arr)
       write(output_unit, '("(", i0, ", ", i0, ")")') arr(i)
    end do
  end block test_tuple2
end program tets_fqsort
```
</details>

実行結果
```bash
$ gfortran test_qsort.f90 && ./a.out
before sort
2 5 1 4 3
after sort
1 2 3 4 5

before sort
(100, 3)
(5, 1)
(4, 3)
(-1, -1)
(5, 500)
after sort
(-1, -1)
(4, 3)
(5, 1)
(5, 500)
(100, 3)
```
- ちゃんとソートされています.
- ユーザ定義型もちゃんとソートできています.
- 危険でも良いので雑にソートを定義したいときには便利なのかもしれません.
  + `c_loc` と `c_funloc` を経由しているため, 型や比較関数を間違えたら `segmentation falut` になるはずです.
  + つまり, 面倒くさいので, Fortran でソート関数を定義した方が楽な可能性(コンパイラが配列外参照のチェックや型チェックをできる)があると思います.


# マイクロ秒単位でプロセスを止めたい!
## \<unistd.h>
Fortran の sleep サブルーチンは[秒単位](https://gcc.gnu.org/onlinedocs/gfortran/SLEEP.html)でプロセスを止めます.
そして, GNU extension なので gfortran でしか使えません. (intel fortran にも同じ関数は[あります](https://www.intel.com/content/www/us/en/docs/fortran-compiler/developer-guide-reference/2023-0/sleep.html)が...)
なので, マイクロ秒単位で止められる `usleep` を C から借ります.

```test_usleep.f90
module usleep_m
  use, intrinsic :: iso_fortran_env
  use, intrinsic :: iso_c_binding
  implicit none
  interface
     subroutine usleep (x) bind(c, name = "usleep")
       import c_int
       integer(c_int), value :: x
     end subroutine usleep
  end interface
end module usleep_m

program test_usleep
  use, intrinsic :: iso_fortran_env
  use usleep_m
  implicit none
  call usleep(500000)
end program test_usleep
```
```bash
$ gfortran test_usleep.f90 && time ./a.out
real    0m0.503s
user    0m0.002s
sys     0m0.000s
```

- ちゃんと0.5秒くらい止まってます.
- マイクロ秒単位で止められると嬉しい場面は...
  + 0.5秒毎に結果を表示するとかでしょうか...
- \<unistd.h> は C の標準ではない?
  + Windows では使えないらしいです.

# まとめ
`iso_c_binding` で C 言語の関数を拝借すれば, Fortran 単体ではできないことや面倒なことから解放される可能性がありそうです.
今回は実用的な例を示せていない気がするので, そのうち, 有名なライブラリの関数を呼び出して有用性を示したいです.

# 参考
## c_binding の例 (nag)
https://www.nag-j.co.jp/fortran/fortran2003/Fortran2003_5_1.html

https://www.nag-j.co.jp/fortran/tips/tips_InteroperabilityWithC.html

## 他の例
https://qiita.com/i__am__kito/items/a61794bb3ad6e0ece67b

https://qiita.com/ijknabla/items/bd7fe074239741c7c7d9

https://qiita.com/sakamoti/items/0fdc989ee530a5639f55

## Fortran sleep
https://gcc.gnu.org/onlinedocs/gfortran/SLEEP.html

https://www.intel.com/content/www/us/en/docs/fortran-compiler/developer-guide-reference/2023-0/sleep.html

## C 言語
https://en.cppreference.com/w/c

https://en.cppreference.com/w/c/io/fprintf