WindowsでQuantum ESPRESSOがMPI並列で動かない原因

# はじめに
Quantum ESPRESSO(以下QEと表記)はオープンソースの第一原理計算ソフトである。シェアはVASPの次くらいあるのではないだろうか。Linuxでは簡単にビルドでき、実行にも問題ない。
Windows MSYS2ではソースに３か所ほど手を加える必要があった。cygwinではソースそのままでコンパイルできる。
しかしビルドできてもMPI並列で動かなかった。どこに原因があるか調べたので報告する。

# 結論
バージョンによるのかもしれないが、特定の環境ではMPI並列で動かない。原因はQEではなく、特定の環境でMPI_ALLREDUCEの挙動がおかしいことにあるようだ。MPIなしでビルドすればこのようなことは起こらない。
QEに限らず、fortranで書かれたMPI並列プログラムは、同様の問題が起きる可能性がある。

|System|MPI|fortran|MPI実行|
|:-:|:-:|:-:|:-:|
|MSYS2 UCRT64|msmpi 10.1.1-10|gfortran 13.2.0|NG|
|cygwin64   |OpenMPI 4.1.5-1|gfortran 11.4.0|NG|
|WSL2 Ubuntu 22.04 LTS|OpenMPI 4.1.2-2ubuntu1|gfortran 11.4.0   |OK|

# ビルド
MSYS2 UCRT64で行った。
QEのドキュメントにも書いてあるが、
./bash_profileに
export LC_ALL=C
を入れておく。
予めWindowsでMicrosoft MPIをインストールしてMSYS2のパスに"/c/Program Files/Microsoft MPI/Bin/"を追加しておく。QEの実行にはWindowsのmpiexecを使うためだ。MSYS2にはmpirunやmpiexecは入っていない。

コンパイラ等のツールのインストール
`pacman -S mingw-w64-ucrt-x86_64-toolchain`
msmpiのツールをインストール
`pacman -S mingw-w64-ucrt-x86_64-msmpi`

お好みで、
`pacman -S mingw-w64-ucrt-x86_64-openblas`
`pacman -S mingw-w64-ucrt-x86_64-fftw`
`pacman -S mingw-w64-ucrt-x86_64-scalapack`

QEのソースを展開し、
qe-7.2\external\devxlib\src\timer.cの
`#include <sys/times.h>`
を削除
qe-7.2\XClib\beefun.cに
`#define random rand`
`#define srandom srand`
を追加。

`./configure`
の後、make.incのFFLAGSに-fallow-invalid-bozを追加。
これは、gfortranでmpif.hのインクルード時にエラーが出るためである。

`make all`
`make install`

ビルド、インストールはできるのだが、実行するとエラーが出る。
cygwinでも試してみたが同じだった。cygwinのときはtimer.cとbeefun.cの修正は必要なかった。

# デバッグ
gdbで追いかけたところ、qe-7.2/FFTXlib/src/fft_types.f90 1177行目のCALL MPI_ALLREDUCEの結果がおかしいことがわかった。以下のようにMPI_ALLREDUCE前後で変数nbの変化を見てみる。

```diff_fortran
@@ -1174,7 +1174,9 @@
       END DO

 #if defined(__MPI)
+      write(6,*) '***nb=', nb(1), nb(2), nb(3)
       CALL MPI_ALLREDUCE( MPI_IN_PLACE, nb, 3, MPI_INTEGER, MPI_MAX, dfft%comm, i )
+      write(6,*) '***max=', nb(1), nb(2), nb(3), 'err=', i
 #endif

 ! ... the size of the 3d FFT matrix depends upon the maximum indices. With
```

mpiexec -n 4で実行すると、
```
 ***nb=          40          63          88
 ***nb=          40          63          89
 ***nb=          40          63          89
 ***nb=          40          63          90
 ***max=           0           0           0 err=           0
 ***max=           0           0           0 err=           0
 ***max=           0           0           0 err=           0
 ***max=           0           0           0 err=           0
```

ここでは、`***max= 40 63 90`が4つ出るはずなのだ。

# MPI＿ALLREDUCEの検証
```Fortran
MPI_ALLREDUCE(SENDBUF, RECVBUF, COUNT, DATATYPE, OP, COMM, IERROR)
        <type> SENDBUF(*), RECVBUF(*)
        INTEGER COUNT, DATATYPE, OP, COMM, IERROR
```
MPI_ALLREDUCEは、全てのプロセスのSENDBUFに演算OPをして、結果を全てのプロセスのRECVBUFに入れるものである。SENDBUF,RECVBUFはDATATYPE型で個数COUNTの配列である。SENDBUFがMPI_IN_PLACEの時には、RECVBUFに演算OPをして結果で置き換える。

![mpi_allreduce.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/487217/a89a0901-6d8a-c85b-cadd-fd4ee4b2b696.png)


この通りの動作をするかをテストするためのプログラムを作って挙動を見てみる。

## C
```allreduce.c
#include <stdio.h>
#include "mpi.h"

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

  int my_rank, num_proc, i, ierr;
  int nb[3];

  MPI_Init(&argc, &argv);
  MPI_Comm_size(MPI_COMM_WORLD, &num_proc);
  MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);

  for(i=0;i<3;i++) nb[i] = i*10 + my_rank;
  printf("my_rank=%d, nb= %d %d %d\n",my_rank, nb[0], nb[1],nb[2]);
  ierr = MPI_Allreduce(MPI_IN_PLACE, nb, 3, MPI_INT, MPI_MAX, MPI_COMM_WORLD);
  printf("my_rank=%d, max= %d %d %d ierr=%d\n",my_rank, nb[0], nb[1],nb[2], ierr);

  MPI_Finalize();

  return 0;
}
```
MSYS2でコンパイルして、実行
```
$ mpicc allreduce.c
```

```
$ mpiexec -n 4 ./a.exe
my_rank=2, nb= 2 12 22
my_rank=2, max= 3 13 23 ierr=0
my_rank=3, nb= 3 13 23
my_rank=3, max= 3 13 23 ierr=0
my_rank=0, nb= 0 10 20
my_rank=0, max= 3 13 23 ierr=0
my_rank=1, nb= 1 11 21
my_rank=1, max= 3 13 23 ierr=0
```
仕様通りの挙動である。
## Fortran
```allreduce.f90
program main

implicit none
include "mpif.h"
integer :: my_rank, num_proc, i, ierr
integer :: nb(3)

call MPI_Init(ierr)
call MPI_Comm_size(MPI_COMM_WORLD, num_proc, ierr)
call MPI_Comm_rank(MPI_COMM_WORLD, my_rank, ierr)

do i=1,3
    nb(i) = i*10 + my_rank
end do
write(6,*) "my_rank=", my_rank, "nb=", (nb(i), i=1,3)
call MPI_Allreduce(MPI_IN_PLACE, nb, 3, MPI_INTEGER, MPI_MAX, MPI_COMM_WORLD, ierr)
write(6,*) "my_rank=", my_rank, "max=", (nb(i), i=1,3), "ierr=", ierr

call MPI_Finalize(ierr)

stop
end
```
コンパイルして、実行
```
$ mpif90 -fallow-invalid-boz allreduce.f90
```

```
$ mpiexec -n 4 ./a.exe
 my_rank=           1 nb=          11          21          31
 my_rank=           2 nb=          12          22          32
 my_rank=           3 nb=          13          23          33
 my_rank=           0 nb=          10          20          30
 my_rank=           0 max=           0           0           0 ierr=           0
 my_rank=           2 max=           0           0           0 ierr=           0
 my_rank=           1 max=           0           0           0 ierr=           0
 my_rank=           3 max=           0           0           0 ierr=           0
```
全てのrankでmax= 13 23 33になるべきなのであるが、全て0になっている。
つまりfortranではMPI_ALLREDUCEが正常に実行されないということだ。
cygwin64でも試してみたが、同じくダメだった。
問題はQEではないので、ここで断念。