Fortran on Windowsで高精度なタイマーを使う

# 概要
Windowsでは，高精度なタイマーとして`QueryPerformanceCounter`が用意されています．Fortranから利用するには，インタフェースを作成します．
[Fortranにおける実行時間の測定](https://qiita.com/implicit_none/items/dd067e9bf5f7b49ce84b)の補足的な内容です．

# 環境
- Windows 10
- gfortran 10.3
- intel oneAPI 2021.1

# 背景
[Fortranにおける実行時間の測定](https://qiita.com/implicit_none/items/dd067e9bf5f7b49ce84b)でも述べましたが，Fortranには時間を計測するための手続がいくつかあります．

- `system_clock`
- `cpu_time`
- `date_and_time`

いずれも一長一短があり，並列計算を想定すると力不足でした．そこで，OpenMPの関数`omp_get_wtime`を用いることを推奨しました．

その後，色々と調べていると，Windowsでの時間計測の精度は仕様よりも悪いという情報がありました．Windowsのタイマーは1秒間に64回しか計測のタイミングがないので，例えば`cpu_time`がミリ秒単位で時間を返すとしても，測定精度は1 msよりも低いとのことです．

一方で，Windowsには[QueryPerformanceCounter](https://docs.microsoft.com/ja-jp/windows/win32/api/profileapi/nf-profileapi-queryperformancecounter)という高精度（正しくは高解像度）なタイマーが用意されています．Fortranの組込手続の中でこれが使われているかは判りませんが，`QueryPerformanceCounter`を呼ぶことができれば，計測精度が改善できると期待されます．

# QueryPerformanceCounterに対するインタフェース

## QueryPerformanceCounterの仕様
[`QueryPerformanceCounter`](https://docs.microsoft.com/ja-jp/windows/win32/api/profileapi/nf-profileapi-queryperformancecounter)のドキュメントによると，`QueryPerformanceCounter`の入出力は次のようになっています．

```C
BOOL QueryPerformanceCounter(
  [out] LARGE_INTEGER *lpPerformanceCount
);
```

`QueryPerformanceCounter`を呼び出すと，成否が`BOOL`として返ってきて，カウンタの値は`LARGE_INTEGER`型の変数に書き込まれます．`LARGE_INTEGER`はC++の共用体で，二つの構造体と一つの`LONGLONG`型変数でメモリを共用しています．`BOOL`と`LONGLONG`型変数が使えばいいので，インタフェースは容易に作成できます．

## QueryPerformanceFrequencyの仕様
カウンタの値を秒に変換するには，`QueryPerformanceFrequency`を呼び出すことによって得られるカウントの頻度[count/sec]値を利用します．[`QueryPerformanceFrequency`](https://docs.microsoft.com/en-us/windows/win32/api/profileapi/nf-profileapi-queryperformancefrequency)のドキュメントによると，次のように入出力が定義されています．

```C
BOOL QueryPerformanceFrequency(
  [out] LARGE_INTEGER *lpFrequency
);
```

利用しているハードウェアにおいて，高解像度のカウンタを利用できる場合，TRUEが返ります．カウント頻度はカウンタの値と同様に`LARGE_INTEGER`型の変数に書き込まれます．なお，このカウント頻度はシステム起動時に固定されるので，1度だけ呼び出してその値を再利用できます．

## インタフェース
仕様に沿ってインタフェースを作成します．C言語の論理値と`LONGLONG`を利用するには，`iso_c_binding`で定義されているkind定数`c_bool`と`c_long_long`を利用します．

```Fortran
interface
    !| パフォーマンスカウンタの現在のカウント値を取得する
    function QueryPerformanceCounter(lPerformanceCount_count) result(is_succeeded) bind(c, name="QueryPerformanceCounter") !&
        use, intrinsic :: iso_c_binding
        implicit none
        integer(c_long_long), intent(out) :: lPerformanceCount_count
            !! パフォーマンスカウンタの現在の値 [count]
        logical(c_bool) :: is_succeeded
            !! 関数が成功した場合は非ゼロ（真）を返す
    end function QueryPerformanceCounter

    !| パフォーマンスカウンタのカウント頻度を取得する
    function QueryPerformanceFrequency(lFrequency_countPerSec) result(is_supported) bind(c, name="QueryPerformanceFrequency") !&
        use, intrinsic :: iso_c_binding
        implicit none
        integer(c_long_long), intent(out) :: lFrequency_countPerSec
            !! パフォーマンスカウンタのカウント頻度 [count/sec]<br>
            !! プログラムを実行しているハードウェアで高解像度カウンタが利用可能な場合に非ゼロの値が書き込まれる
        logical(c_bool) :: is_supported
            !! ハードウェアで高解像度カウンタが利用可能な場合に非ゼロ（真）を返す
    end function QueryPerformanceFrequency
end interface
```

`QueryPerformanceCounter`および`QueryPerformanceFrequency`はkernel32.dllで定義されているので，ライブラリへのリンクは特段必要ありません．

インタフェースで定義したカウンタは，下記のように利用します．

```Fortran
integer(c_long_long) :: time_begin_qhc, time_end_qhc, freq
logical(c_bool) :: is_supported, is_succeeded

is_supported = QueryPerformanceFrequency(freq)
if(is_supported)then
    is_succeeded = QueryPerformanceCounter(time_begin_qhc)
    ! 測定したい処理
    is_succeeded = QueryPerformanceCounter(time_end_qhc)
    print *,real(time_end_qhc - time_begin_qhc)/real(freq), "sec"
end if
```

## 使用例

上記のインタフェースをモジュール`profile`に定義し，そのファイル名を`profile.f90`とします．

1秒だけ動作を停止するプログラムを実行してみます．

```Fortran
program main
    use,intrinsic :: iso_c_binding
    use,intrinsic :: iso_fortran_env
    use :: profile
    implicit none

    integer(c_long_long) :: time_begin_qhc, time_end_qhc, freq
    logical(c_bool) :: is_supported, is_succeeded

    is_supported = QueryPerformanceFrequency(freq)
    if(is_supported)then
        is_succeeded = QueryPerformanceCounter(time_begin_qhc)

        call sleep(1)

        is_succeeded = QueryPerformanceCounter(time_end_qhc)
        print *,real(time_end_qhc - time_begin_qhc)/real(freq), "sec"
    end if
end program main
```

コンパイルには，まず`profile.f90`をコンパイルしてオブジェクトファイルを作成し，それと共に`main.f90`をコンパイルします．

```console
gfortran -c profile.f90
gfortran main.f90 profile.o
```

```console
ifort /c profile.f90
ifort main.f90 profile.obj
```

Intel Fortranでコンパイルして作成された実行ファイルを5回ほど実行してみると，それぞれ1.000138, 1.000350, 1.000226, 1.000275という結果が得られました．
`system_clock`に置き換えて測定してみると，1.001000, 1.001000, 1.001000, 1.000000, 1.001000という結果でした．`system_clock`は1ミリ秒単位で計測でき，`QueryPerformanceCounter`はそれよりも高い解像度で計測できていると考えられます．

# 並列プログラムの実行時間の測定
[Fortranにおける実行時間の測定](https://qiita.com/implicit_none/items/dd067e9bf5f7b49ce84b)で時間を測定したプログラム（補足で言及している，ループの順序を修正した方）に`QueryPerformanceCounter`を追加して，実行時間を測定してみます．

<details><summary>main.f90</summary><div>

```Fortran
program main
    use,intrinsic :: iso_c_binding
    use,intrinsic :: iso_fortran_env
    use :: profile
    !$ use omp_lib
    implicit none

     enum, bind(c)
        enumerator :: Year = 1
        enumerator :: Month
        enumerator :: Day
        enumerator :: TimeDifference_min
        enumerator :: Hour
        enumerator :: Minute
        enumerator :: Second
        enumerator :: Millisecond
    end enum

    integer(int32),parameter :: N = 2**13
    real(real64),allocatable :: a(:,:),b(:,:),c(:,:)
    integer(int32) :: i,j

    ! system_clock
    integer(int64) :: time_begin_c,time_end_c,count_per_sec

    ! cpu_time
    real(real32)   :: time_begin_s,time_end_s

    ! omp_get_wtime
    !$ real(real64)   :: time_begin_ws,time_end_ws     ! w (means wall-clock time) is added to avoid name conflict

    ! date_and_time
    integer :: time_begin_values(8),time_end_values(8) ! integer with default kind
    real(real32),parameter :: to_second = 1e-3

    ! QueryPerformanceCounter
    integer(c_long_long) :: time_begin_qhc, time_end_qhc, freq
    logical(c_bool) :: is_supported, is_succeeded

    allocate(a(N,N))
    allocate(b(N,N))
    allocate(c(N,N))

    ! set number of threads
    !$ call omp_set_num_threads(4)

    !$omp parallel
    ! initialize arrays
    !$omp do
    do j = 1,N
    do i = 1,N
        a(i,j) = 1._real64
        b(i,j) = 2._real64
        c(i,j) = 0._real64
    end do
    end do
    !$omp end do

    ! system_clock
    !$omp master
    call system_clock(time_begin_c)
    !$omp end master
    !$omp do
    do j = 1,N
    do i = 1,N
        c(i,j) = a(i,j) + b(i,j)
    end do
    end do
    !$omp end do
    !$omp master
    call system_clock(time_end_c,count_per_sec)
    print *,real(time_end_c - time_begin_c)/count_per_sec,"sec",sum(c)/N**2
    !$omp end master

    ! cpu_time
    !$omp master
    call cpu_time(time_begin_s)
    !$omp end master
    !$omp do
    do j = 1,N
    do i = 1,N
        c(i,j) = a(i,j) + b(i,j)
    end do
    end do
    !$omp end do
    !$omp master
    call cpu_time(time_end_s)
    print *,time_end_s - time_begin_s,"sec",sum(c)/N**2
    !$omp end master

    ! date_and_time
    !$omp master
    call date_and_time(values = time_begin_values)
    !$omp end master
    !$omp do
    do j = 1,N
    do i = 1,N
        c(i,j) = a(i,j) + b(i,j)
    end do
    end do
    !$omp end do
    !$omp master
    call date_and_time(values = time_end_values)
    print *,time_end_values(     Second)-time_begin_values(     Second) &
          +(time_end_values(Millisecond)-time_begin_values(Millisecond))*to_second,"sec",sum(c)/N**2
    !$omp end master

    ! omp_get_wtime
    !$omp master
    !$ time_begin_ws = omp_get_wtime()
    !$omp end master
    !$omp do
    do j = 1,N
    do i = 1,N
        c(i,j) = a(i,j) + b(i,j)
    end do
    end do
    !$omp end do
    !$omp master
    !$ time_end_ws = omp_get_wtime()
    !$ print *,real(time_end_ws - time_begin_ws,real32),"sec",sum(c)/N**2
    !$omp end master

    ! query performance counter
    is_supported = QueryPerformanceFrequency(freq)
    if(is_supported)then
        !$omp master
        is_succeeded = QueryPerformanceCounter(time_begin_qhc)
        !$omp end master
        !$omp do
        do j = 1,N
        do i = 1,N
            c(i,j) = a(i,j) + b(i,j)
        end do
        end do
        !$omp end do
        !$omp master
        is_succeeded = QueryPerformanceCounter(time_end_qhc)
        print *,real(time_end_qhc - time_begin_qhc,real32)/real(freq,real32), "sec", sum(c)/N**2
        !$omp end master
    end if

    !$omp end parallel

    deallocate(a)
    deallocate(b)
    deallocate(c)
end program main
```
</div></details>

## 実行結果
5回測定してその結果を表にまとめました．実行時間は全て秒です．

### 1並列

|コンパイラ|時間測定の手続|run 1|run 2|run 3|run 4|run 5|average|
|:--|:--|:--|:--|:--|:--|:--|:--|
|Intel   |`system_clock`           |0.1260000|0.1280000|0.1280000|0.1260000|0.1220000    |1.270000E-01|
|        |`cpu_time`               |0.1250000|0.1093750|0.1250000|0.1250000|9.3750000E-02|1.156250E-01|
|        |`date_and_time`          |0.1250000|0.1240000|0.1260000|0.1250000|0.1289999    |1.250000E-01|
|        |`omp_get_wtime`          |0.1334495|0.1257561|0.1223831|0.1261075|0.1314339    |1.269241E-01|
|        |`QueryPerformanceCounter`|0.1245430|0.1262833|0.1275150|0.1266859|0.1280735    |1.262568E-01|
|gfortran|`system_clock`           |0.125433505|0.129628599|0.131128594|0.124462999|0.125900403    |1.276634E-01|
|        |`cpu_time`               |0.125000000|0.125000000|0.140625000|0.125000000| 9.37500000E-02|1.218750E-01|
|        |`date_and_time`          |0.133000001|0.130000010|0.146000013|0.131999999|0.130000010    |1.352500E-01|
|        |`omp_get_wtime`          |0.131000042|0.135999918|0.133999825|0.131999969|0.129999876    |1.332499E-01|
|        |`QueryPerformanceCounter`|0.131263107|0.135557204|0.134244293|0.127931297|0.130921006    |1.322490E-01|

### 2並列

|コンパイラ|時間測定の手続|run 1|run 2|run 3|run 4|run 5|average|
|:--|:--|:--|:--|:--|:--|:--|:--|
|Intel   |`system_clock`           |6.7000002E-02|6.8999998E-02|6.8999998E-02|7.1000002E-02|7.0000000E-02|6.920000E-02|
|        |`cpu_time`               |0.1562500    |0.1250000    |0.1250000    |0.1250000    |0.1250000    |1.312500E-01|
|        |`date_and_time`          |6.5000005E-02|6.6000000E-02|6.5000005E-02|6.8000004E-02|6.9000006E-02|6.660000E-02|
|        |`omp_get_wtime`          |6.3780501E-02|6.6513300E-02|6.6037901E-02|6.5664798E-02|6.7219697E-02|6.584324E-02|
|        |`QueryPerformanceCounter`|6.6515401E-02|6.5232597E-02|6.5742902E-02|6.9207802E-02|6.8397403E-02|6.701922E-02|
|gfortran|`system_clock`           |6.89190999E-02|6.87130988E-02|6.91405013E-02|7.36619979E-02|7.01569989E-02|7.011834E-02|
|        |`cpu_time`               |6.25000000E-02|7.81250000E-02|7.81250000E-02|6.25000000E-02|7.81250000E-02|7.187500E-02|
|        |`date_and_time`          |7.00000003E-02|6.99999332E-02|6.59999996E-02|7.10000023E-02|7.60000050E-02|7.059999E-02|
|        |`omp_get_wtime`          |6.80000782E-02|6.50000572E-02|6.90000057E-02|6.90000057E-02|6.69999123E-02|6.760001E-02|
|        |`QueryPerformanceCounter`|7.04241991E-02|6.68191984E-02|6.70834035E-02|6.53799027E-02|6.52005970E-02|6.698146E-02|

### 4並列

|コンパイラ|時間測定の手続|run 1|run 2|run 3|run 4|run 5|average|
|:--|:--|:--|:--|:--|:--|:--|:--|
|Intel   |`system_clock`           |5.0999999E-02|4.8000000E-02|4.8999999E-02|5.0000001E-02|5.0999999E-02|4.980000E-02|
|        |`cpu_time`               |9.3750000E-02|0.1093750    |0.1093750    |4.6875000E-02|9.3750000E-02|7.812500E-02|
|        |`date_and_time`          |3.2000002E-02|3.4000002E-02|3.1000001E-02|3.3000000E-02|3.3000000E-02|3.260000E-02|
|        |`omp_get_wtime`          |3.0819200E-02|2.9993400E-02|3.1744301E-02|3.2902401E-02|3.0556601E-02|3.120318E-02|
|        |`QueryPerformanceCounter`|3.0747401E-02|3.0455699E-02|3.0006699E-02|3.0873900E-02|3.1835299E-02|3.078380E-02|
|gfortran|`system_clock`           |5.00889011E-02|5.07357009E-02|5.16671017E-02|5.48393987E-02|5.71402013E-02|5.289426E-02|
|        |`cpu_time`               |1.56250000E-02|3.12500000E-02|3.12500000E-02|1.56250000E-02|3.12500000E-02|2.500000E-02|
|        |`date_and_time`          |3.40000018E-02|3.29999998E-02|3.50000001E-02|3.29999998E-02|3.29999998E-02|3.360000E-02|
|        |`omp_get_wtime`          |3.29999924E-02|3.39999199E-02|3.29999924E-02|3.29999924E-02|3.39999199E-02|3.339996E-02|
|        |`QueryPerformanceCounter`|3.27327996E-02|3.13258991E-02|3.51462997E-02|3.36192995E-02|3.35536003E-02|3.327558E-02|

### スレッド数による実行時間の変化
|コンパイラ|時間測定の手続|1並列|2並列|4並列|
|:--|:--|:--|:--|:--|
|Intel   |`system_clock`           |1.270000E-01|6.920000E-02|4.980000E-02|
|        |`cpu_time`               |1.156250E-01|1.312500E-01|7.812500E-02|
|        |`date_and_time`          |1.250000E-01|6.660000E-02|3.260000E-02|
|        |`omp_get_wtime`          |1.269241E-01|6.584324E-02|3.120318E-02|
|        |`QueryPerformanceCounter`|1.262568E-01|6.701922E-02|3.078380E-02|
|gfortran|`system_clock`           |1.276634E-01|7.011834E-02|5.289426E-02|
|        |`cpu_time`               |1.218750E-01|7.187500E-02|2.500000E-02|
|        |`date_and_time`          |1.352500E-01|7.059999E-02|3.360000E-02|
|        |`omp_get_wtime`          |1.332499E-01|6.760001E-02|3.339996E-02|
|        |`QueryPerformanceCounter`|1.322490E-01|6.698146E-02|3.327558E-02|

## 時間測定ルーチンの傾向

Intel Fortranの場合，`system_clock`および`date_and_time`は1ミリ秒単位で計測できるようですが，`cpu_time`の解像度はそこまで細かくないようです．これが背景で述べた1秒間に64回の解像度しかないということでしょうか．`omp_get_wtime`と`QueryPerformanceCounter`は同じような傾向を示しており，他3手続よりも高い解像度で測定できるようです．

gfortranの場合，`date_and_time`は1ミリ秒の解像度であること，`cpu_time`の解像度はそれよりも粗いことはIntel Fortranと同じです．`system_clock`は，Intel Fortranよりも解像度が高く，逆に`omp_get_wtime`は解像度が低いように見受けられます．

`QueryPerformanceCounter`は，同じ実装を呼び出しているので解像度は同じです．

# まとめ
Windowsで高精度なタイマーとして`QueryPerformanceCounter`を用いる方法を説明しました．インタフェースを作成するだけで，コンパイラに依存せずに利用できますし，使い勝手も`system_clock`とほぼ同じです．インタフェースの作成も，Fortran標準で用意されている`iso_c_binding`を用いれば，簡単に作成できます．マルチスレッド並列化した場合でも，正しく経過時間を測定できます．

Windowsで実行時間を測定する場合は，標準の方法として用いてもよいと思います．