MPIでいい感じに因数分解


# はじめに

これは[MPI Advent Calendar 2017](https://adventar.org/calendars/2548)の3日目の記事ということにします。

# スパコンとインターコネクトについて

昔はどうか知りませんが、現在のスパコンは、ほぼ普通のPCと同じ要素で構成されています。では普通のPCクラスタとスパコンの違いはなんでしょうか？普通のPCかそれに近いものをただ繋いだだけでスパコンになるのでしょうか？

「スーパーコンピュータ」という言葉がそもそも確固たる定義がない言葉であり、人によって定義も違うでしょうから、ただPCやノードを繋いだだけのものをスパコンと呼ぶ人もいるでしょうが、個人的に「スパコンとPCクラスタの違いはなにか？」と聞かれたら「インターコネクトと信頼性」と答えています。

信頼性についてはこのAdCのテーマから外れるので詳細は触れませんが、例えば通常のPCであれば1年から数年でどこか壊れるのが普通だと思います。ざっくりした話、だいたい一年に一度壊れるようなノードを1000ノード集めると、8〜9時間ごとにノードがどこか壊れることになります。従って、このシステムでは、1000ノードを24時間協調動作させることは絶望的です。HPLで性能出そうとすると数十時間とか計算する必要があったりするので、一年に一度壊れるようなノードを集めたシステムでは1000ノードでもまともなベンチマークすら取れないことになります。

もう一つの重要な要素はインターコネクトです。これも詳細は専門家がどこかに書くでしょうが、1000ノード以上がまともに協調動作するようなインターコネクトを組むのは全く非自明な技術です。ユーザから見て一番うれしいのは、「どのノードから見ても、別のノードが同じような距離にいる」ことです。こういうネットワーク(クロスバースイッチでファットツリーとか)を組むと、ユーザは自分のプロセスがどのノードにどういう形で割り当てられたか気にする必要がなくなります。また、ジョブスケジューラが楽になる、というのも見逃せないポイントです。ユーザのジョブのディスパッチにネットワーク形状を気にしなくてよくなるため、純粋に「空きノード数」だけを監視し、空いたところに放り込めば良いことになります[^1]。

[^1]: 逆に、6次元トーラスネットワークに、各ユーザが「Xノード取れればなんでもいい」「X*Yの形で欲しい」「X*Y*Zの形じゃないとイヤだ」とかてんでバラバラにジョブを投げてきた時にのジョブのディスパッチ(個人的に高次元テトリスと呼んでます)とか考えると、僕が担当者なら胃が痛くなります。

さて、良いことづくめのファットツリーですが、なぜ全部のサイトが採用せずに独自のネットワークを組んだりするかというと、それは高いからです。ノード数が増えると滅茶苦茶高額になります。そもそもスパコンを全系使うジョブというのは稀なので、そこにものすごくお金をかけるのはどうか、という気もしてきます。なので、例えばネットワーク形状をトーラスにしたり、もしくはドラゴンフライにしたりします。このあたりはベンダーの腕の見せ所です。

トーラスというのは、例えば高次元の周期境界でノードを繋いだものです。旧SGI(現HPE)は3次元トーラスを二枚重ねたハイパーキューブ構造を作りますし、富士通のTOFUは6次元トーラスを組むことになります。トーラスの頂点には、普通いくつかのノードをまとめたものがぶら下がります。この小さなまとまりがネットワークの最小構造になります。旧SGIのだと18ノードがひとまとまり、TOFUは12ノードがひとまとまりです。

ユーザから見た場合、プロセス数が2のベキだといろいろ便利です。典型的には高速フーリエ変換で、要素数が2のベキだとプログラムが簡単になります。それ以外でも、何か空間をツリー形状に再帰的に分割とかしようと思うと、2のベキじゃなくなった途端にいろいろ面倒になります。しかし、先程みたとおり、多くのサイトで「ネットワークの最小構成要素」のノード数が2のベキではなくなっています。従って、ユーザは素因数3を含むようなバタフライ演算的なものを組まなきゃいけなくなったりして面倒なわけですね。

話が逸れました。とにかく、現在のスパコンでは、ネットワークをいい感じに使おうと思うと、プロセス数が2のベキではなくなる、ということです。

# 因数分解

さて、面倒な話を全部忘れましょう。3次元の領域分割で計算しようと思った時、なるべくx, y, z方向の分割数をそろえたいですよね。例えば3360ノードの3次元分割なら16*15*14とか、2805なら17*15*11とか、そういう分割をプログラム側で自動でやりたくなりますが、ちょっと自分で組もうとするとわりと面倒です(興味のある人は組んでみると面白いと思います)。

で、そういうのはMPIで頻出するので、MPI側でそういう関数が用意されています、`MPI_Dims_create`という関数です。こんな感じに使います。

```test.cpp
#include <cstdio>
#include <cstdlib>
#include <mpi.h>
int
main(int argc, char**argv){
  MPI_Init(&argc, &argv);
  int n = 1;
  if(argc >1){
    n = atoi(argv[1]);
    int d2[2] = {}, d3[3] = {};
    MPI_Dims_create(n, 2, d2);
    MPI_Dims_create(n, 3, d3);
    printf("%d \t = %d * %d \t = %d * %d * %d\n",n, d2[0], d2[1], d3[0], d3[1], d3[2]);
  }else{
    printf("usage: ./a.out n\n");
  }
  MPI_Finalize();
}
```

`MPI_Dims_create`は、第一引数に因数分解する数、第二引数に次元、第三引数に配列を渡します。この時、配列に0以外の要素が入っていると、それを尊重して分割しようとします。

実行結果はこんな感じです。

```shell-session
$ ./a.out 3360
3360 	 = 60 * 56 	 = 20 * 14 * 12

$ ./a.out 2805 
2805 	 = 55 * 51 	 = 17 * 15 * 11
```

三次元分割の配列の最初に3を入れてみましょう。

```test2.cpp
#include <cstdio>
#include <cstdlib>
#include <mpi.h>
int
main(int argc, char**argv){
  MPI_Init(&argc, &argv);
  int n = 1;
  if(argc >1){
    n = atoi(argv[1]);
    int d2[2] = {}, d3[3] = {3,0,0}; // 最初の要素に3を入れた
    MPI_Dims_create(n, 2, d2);
    MPI_Dims_create(n, 3, d3);
    printf("%d \t = %d * %d \t = %d * %d * %d\n",n, d2[0], d2[1], d3[0], d3[1], d3[2]);
  }else{
    printf("usage: ./a.out n\n");
  }
  MPI_Finalize();
}
```

`MPI_Dims_create`はそれを尊重し、「3*Y*Z」という分割を試みます。

```shell-session
$ ./a.out 3360
3360 	 = 60 * 56 	 = 3 * 40 * 28

$ ./a.out 2805
2805 	 = 55 * 51 	 = 3 * 55 * 17
```

三次元分割が「3*Y*Z」という形になったのがわかると思います。ただし、ユーザの指定した分割が不可能だった場合、エラーとなります(処理系によってエラーメッセージは変わります)。

```shell-session
$ ./a.out 256
*** An error occurred in MPI_Dims_create
*** reported by process [2450391041,0]
*** on communicator MPI_COMM_WORLD
*** MPI_ERR_DIMS: invalid topology dimension
*** MPI_ERRORS_ARE_FATAL (processes in this communicator will now abort,
***    and potentially your MPI job)
```

# まとめ

プロセス数を食わすといい感じに因数分解してくれる関数、`MPI_Dims_create`を紹介しました。っていうか僕がこれ知ったの最近です。知ってる人は知ってるんでしょうが、有名なんですかね？ググっても2,790件しかひっかからないんですが。`MPI_Init`は6万件、`MPI_Sendrecv`が1万件ちょっとひっかかるので、相対的にマイナーなんですかね？似たような関数自作したりしてたので、もっと早く知りたかった・・・