Day 1: 環境構築
-スパコン上で実行されるプログラムは並列プログラムである。したがって「スパコンを使う」ということは、 狭義には「並列化されたプログラムを実行する」ということを意味する。したがって、誰かが作った並列プログラムをスパコン上で実行すれば、スパコンは使えることになる。 それはそれでOKなのだが、本稿のタイトルは「一週間でなれる！スパコンプログラマ」であるから、スパコン上で動くコードを開発できるようになることを目的とする。 それはすなわち、「並列プログラミングをする」ということである。「並列プログラミング」という字面を見ると「難しそう」という印象を持つかもしれない。 しかし、(世の中の多くの「一見難しそうなもの」がそうであるように)並列プログラミングはさほど難しくない。 「一週間でなれる！スパコンプログラマ」の初日は、まずローカルに並列プログラミング環境を構築し、並列プログラミングに触れてみるところからはじめてみよう。
+ +スパコン上で実行されるプログラムは並列プログラムである。したがって「スパコンを使う」ということは、 狭義には「並列化されたプログラムを実行する」ということを意味する。したがって、誰かが作った並列プログラムをスパコン上で実行すれば、スパコンは使えることになる。 それはそれでOKなのだが、本稿のタイトルは「一週間でなれる！スパコンプログラマ」であるから、スパコン上で動くコードを開発できるようになることを目的とする。 それはすなわち、「並列プログラミングをする」ということである。「並列プログラミング」という字面を見ると「難しそう」という印象を持つかもしれない。 しかし、(世の中の多くの「一見難しそうなもの」がそうであるように)並列プログラミングはさほど難しくない。 「一週間でなれる！スパコンプログラマ」の初日は、まずローカルに並列プログラミング環境を構築し、並列プログラミングに触れてみるところからはじめてみよう。
MPIとは
一口に「並列化」といっても、様々な種類がありえる。一般に使われている並列プログラミングモデルは、「データ並列」「共有メモリ並列」「分散メモリ並列」の三種類であろう。 以後、プロセスやスレッドといった単語についてかなりいい加減な言葉遣いをするため、ちゃんと学びたい人はちゃんとした書籍を参考にされたい。特にWindowsとLinuxのプロセスの違いとか言い出すと話が長くなるので、ここでは説明しない。また、データ並列についてはとりあえずおいておく。
「共有メモリ並列」とは、並列単位がメモリを共有する並列化方法である。 通常は並列単位としてスレッドを用いるので、ここでは「スレッド並列」と呼ぶ。 逆に「分散メモリ並列」とは、並列単位がメモリを共有しない並列化方法である。 通常は並列単位としてプロセスを用いるので、ここでは「プロセス並列」と呼ぶ。 また、「プロセス並列」と「スレッド並列」を両方行う「ハイブリッド並列」という並列化もある。
@@ -134,7 +135,7 @@MPIのインストール
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.したがって、インクルードパスやリンクの設定を明示的にするならば、mpic++を呼び出す必要はない。 スパコンサイトによっては、環境変数でMPIのインクルードパスが設定されている場合もあるだろう。その場合は単にg++でもicpcでも、MPIを用いたコードがそのままコンパイルできる。ただし、リンクのために-lmpiの指定が(場合によっては-lmpi_cxxも)必要なので注意。
はじめてのMPI
-環境構築ができたら、こんなコードを書いて、hello.cppという名前で保存してみよう。
+環境構築ができたら、こんなコードを書いて、hello.cppという名前で保存してみよう。
#include <cstdio>
#include <mpi.h>
@@ -174,8 +175,7 @@ ランク
MPIでは、起動したプロセスに通し番号が振られる。その通し番号のことを ランク(rank) と呼ぶ。 ランクの取得にはMPI_Comm_rank関数を使う。
int rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
-これを実行すると変数rankにランク番号が入る。N並列している場合、ランクは0からN-1までである。 試してみよう。
-
+これを実行すると変数rankにランク番号が入る。N並列している場合、ランクは0からN-1までである。 試してみよう。以下をrank.cppとして保存し、コンパイル、実行してみよう。
#include <cstdio>
#include <mpi.h>
@@ -235,7 +235,7 @@ 標準出力について
mpirun -np 4 ./a.out
などとしてMPIプログラムを実行したとする。この場合は4プロセス立ち上がり、それぞれにPIDが与えられ、固有のメモリ空間を持つ。しかし、これらのプロセスは標準出力は共有している。 したがって、「せーの」で標準出力に出力しようとしたら競合することになる。 この時、例えば他のプロセスが出力している時に他のプロセスが書き込んだり、出力が混ざったりしないように、 後ろで交通整理が行われる。そもそも画面になにかを表示する、というのはわりと奥が深いのだが、 そのあたりの話は tanakamura さんの 実践的低レベルプログラミング とかを読んでほしい。
さて、とにかく今は標準出力というリソースは一つしかないのに、4つのプロセスがそこを使う。 この時、「あるひとかたまりの処理については、一つのプロセスが独占して使う」ようにすることで 競合が起きないようにする。 「ひとかたまりの処理」とは、例えば「printfで出力を始めてから終わるまで」である。
-例えば先程のrank.cppの例では、
+例えば先程のrank.cppの例では、
printf("Hello! My rank is %d\n", rank);
という命令があった。ここでは、まず出力すべき文字列、例えば Hello! My rank is 0を作る。そして、せーので書き出す。 イメージとしては
puts("Hello! My rank is 0");
@@ -248,9 +248,7 @@ 標準出力について
puts("Hello! My rank is 1");
puts("Hello! My rank is 3");
とかになるだけで、さほど表示は乱れない。
-さて、同様なプログラムをstd::coutで書いてみよう。
-こんな感じになると思う。
-
+さて、同様なプログラムをstd::coutで書いてみよう。以下をrank_stream.cppという名前で保存、コンパイル、実行してみよう。
#include <iostream>
#include <mpi.h>
@@ -340,9 +338,7 @@ GDBによるMPIプログ
gdbで変数をいじって無限ループを脱出させる
あとは好きなようにデバッグする
-という方針でいく。なお、なぜかMac OSではMPIプロセスへのgdbのアタッチがうまくいかなかったので、以下はCentOSで実行している。
-こんなコードを書く。
-
+という方針でいく。なお、なぜかMac OSではMPIプロセスへのgdbのアタッチがうまくいかなかったので、以下はCentOSで実行している。以下をgdb_mpi.cppという名前で保存しよう。
#include <cstdio>
#include <sys/types.h>
#include <unistd.h>
diff --git a/day2/README.md b/day2/README.md
index 064de85..ff833c8 100644
--- a/day2/README.md
+++ b/day2/README.md
@@ -1,8 +1,8 @@
# Day 2 : スパコンの使い方
-## はじめに
-
+
スパコンを使うのに、必ずしもスパコンがどのように構成されているかを知る必要はない。しかし、せっかくスパコンを使うのだから、スパコンとは何かについて簡単に知っておいても良いであろう。ただし、こういう単語にありがちだが「何がスパコンか」は人によって大きく異なる。ここで紹介するのはあくまで「執筆者が思うスパコンの定義」の説明であり、他の人は他の定義があることを承知されたい。ここは、「読むとなにかができるようになる」というよりは、「スパコンを使ったことがない人が、将来スパコンを使うにあたって知っておくと良さそうなこと」を書いておく。特に手を動かすところはない。読み物として流して読んでいただければ良い。
+
## スパコンとは
diff --git a/day2/index.html b/day2/index.html
index 2378a93..f818014 100644
--- a/day2/index.html
+++ b/day2/index.html
@@ -72,8 +72,8 @@
Day 2 : スパコンの使い方
-はじめに
-スパコンを使うのに、必ずしもスパコンがどのように構成されているかを知る必要はない。しかし、せっかくスパコンを使うのだから、スパコンとは何かについて簡単に知っておいても良いであろう。ただし、こういう単語にありがちだが「何がスパコンか」は人によって大きく異なる。ここで紹介するのはあくまで「執筆者が思うスパコンの定義」の説明であり、他の人は他の定義があることを承知されたい。ここは、「読むとなにかができるようになる」というよりは、「スパコンを使ったことがない人が、将来スパコンを使うにあたって知っておくと良さそうなこと」を書いておく。特に手を動かすところはない。読み物として流して読んでいただければ良い。
+
+スパコンを使うのに、必ずしもスパコンがどのように構成されているかを知る必要はない。しかし、せっかくスパコンを使うのだから、スパコンとは何かについて簡単に知っておいても良いであろう。ただし、こういう単語にありがちだが「何がスパコンか」は人によって大きく異なる。ここで紹介するのはあくまで「執筆者が思うスパコンの定義」の説明であり、他の人は他の定義があることを承知されたい。ここは、「読むとなにかができるようになる」というよりは、「スパコンを使ったことがない人が、将来スパコンを使うにあたって知っておくと良さそうなこと」を書いておく。特に手を動かすところはない。読み物として流して読んでいただければ良い。
スパコンとは
普通のPCは、CPU、メモリ、ネットワーク、ディスクなどから構成されている。スパコンも全く同様に、CPU、メモリ、ネットワーク、ディスクがある。 それぞれちょっと高級品を使っているだけで、基本的には普通のPCと同じと思って良い。ただし、PCとはつなぎ方がちょっと異なる。 スパコンは、CPUとメモリをまとめたものを「ノード」と呼ぶ。このノードをたくさん集めて高速なネットワークでつないだものがスパコン本体である。普通のPCではCPUの近くにディスクがあるが、最近のスパコンのノードはディスクレスの構成にすることが多い。そのかわり、大きなファイルシステムとネットワークでつなぐ。
diff --git a/day3/README.md b/day3/README.md
index 7301b0c..2edb064 100644
--- a/day3/README.md
+++ b/day3/README.md
@@ -1,7 +1,6 @@
# Day 3 : 自明並列
-## 自明並列、またの名を馬鹿パラとは
-
+
例えば、100個の画像データがあるが、それらを全部リサイズしたい、といったタスクを考える。
それぞれのタスクには依存関係が全くないので、全部同時に実行してもなんの問題もない。
したがって、100並列で実行すれば100倍早くなる。
@@ -12,7 +11,8 @@
その意義は大きい。
なにはなくとも、まず馬鹿パラができないことには非自明並列もできないわけだし、馬鹿パラができるだけでも、できない人に比べて
圧倒的な攻撃力を持つことになる。
-ここでは、まず馬鹿パラのやり方を見てみよう。
+本章では、まず馬鹿パラのやり方を見てみよう。
+
## 自明並列の例1: 円周率
diff --git a/day3/index.html b/day3/index.html
index f0ad5ae..cc9c3ff 100644
--- a/day3/index.html
+++ b/day3/index.html
@@ -73,11 +73,11 @@
Day 3 : 自明並列
-自明並列、またの名を馬鹿パラとは
-例えば、100個の画像データがあるが、それらを全部リサイズしたい、といったタスクを考える。 それぞれのタスクには依存関係が全くないので、全部同時に実行してもなんの問題もない。 したがって、100並列で実行すれば100倍早くなる。 このように、並列タスク間で依存関係や情報のやりとりが発生しない並列化のことを自明並列と呼ぶ。 英語では、Trivial Parallelization(自明並列)とか、Embarrassingly parallel(馬鹿パラ)などと表現される。 「馬鹿パラ」とは「馬鹿でもできる並列化」の略で(諸説あり)、その名の通り簡単に並列化できるため、 文字通り馬鹿にされることも多いのだが、並列化効率が100%であり、最も効率的に計算資源を利用していることになるため、 その意義は大きい。 なにはなくとも、まず馬鹿パラができないことには非自明並列もできないわけだし、馬鹿パラができるだけでも、できない人に比べて 圧倒的な攻撃力を持つことになる。 ここでは、まず馬鹿パラのやり方を見てみよう。
+
+例えば、100個の画像データがあるが、それらを全部リサイズしたい、といったタスクを考える。 それぞれのタスクには依存関係が全くないので、全部同時に実行してもなんの問題もない。 したがって、100並列で実行すれば100倍早くなる。 このように、並列タスク間で依存関係や情報のやりとりが発生しない並列化のことを自明並列と呼ぶ。 英語では、Trivial Parallelization(自明並列)とか、Embarrassingly parallel(馬鹿パラ)などと表現される。 「馬鹿パラ」とは「馬鹿でもできる並列化」の略で(諸説あり)、その名の通り簡単に並列化できるため、 文字通り馬鹿にされることも多いのだが、並列化効率が100%であり、最も効率的に計算資源を利用していることになるため、 その意義は大きい。 なにはなくとも、まず馬鹿パラができないことには非自明並列もできないわけだし、馬鹿パラができるだけでも、できない人に比べて 圧倒的な攻撃力を持つことになる。 本章では、まず馬鹿パラのやり方を見てみよう。
自明並列の例1: 円周率
まず、自明並列でよく出てくる例として、サンプリング数を並列化で稼ぐ方法を見てみよう。とりあえず定番の、 モンテカルロ法で円周率を計算してみる。
-こんなコードを書いて、calc_pi.cppという名前で保存してみよう。
+こんなコードを書いて、calc_pi.cppという名前で保存してみよう。
#include <cstdio>
#include <random>
#include <algorithm>
@@ -113,7 +113,7 @@ 自明並列の例1: 円周率
ランク番号を乱数の種に使う
そのままcalc_piを呼ぶ。
-以上の修正をしたコードをcalc_pi_mpi.cppという名前で作成する。
+以上の修正をしたコードをcalc_pi_mpi.cppという名前で作成する。
#include <cstdio>
#include <random>
#include <algorithm>
@@ -178,7 +178,7 @@ 自明並列の例1: 円周率
45162 a.out 89.1 00:12.47 3/1 0 15 2620K 0B 0B 45162
4並列実行したので、45162から45165まで4つのプロセスが起動され、実行していることがわかる。 このように、なにか統計平均を取りたい時、並列化によってサンプル数を稼ぐ並列化をサンプル並列と呼ぶ。
自明並列テンプレート
-先程の並列プログラムcalc_pi_mpi.cppのmain関数はこうなっていた。
+先程の並列プログラムcalc_pi_mpi.cppのmain関数はこうなっていた。
int main(int argc, char **argv) {
MPI_Init(&argc, &argv);
int rank;
@@ -227,8 +227,7 @@ 自明並列の例2:
さて、ファイル数はともかく、プロセス数がハードコーディングされているのが気になる。 MPIのプログラムは、実行時にプロセス数を自由に指定することができる。 実行するプロセス数を変えるたびにコンパイルし直すのは面倒だ。 というわけで、実行時に総プロセス数を取得する関数MPI_Comm_sizeが用意されている。 使い方はMPI_Comm_rankと同じで、
int procs;
MPI_Comm_size(MPI_COMM_WORLD, &procs)
-とすれば、procsにプロセス数が入る。これを使うと、先程のコードは こんな感じにかける。
-
+とすれば、procsにプロセス数が入る。これを使うと、先程のコードは こんな感じにかける(processfiles.cpp)。
#include <cstdio>
#include <mpi.h>
@@ -264,7 +263,7 @@ 自明並列の例3: 統計処理
double pi_sum = 0.0;
MPI_Allreduce(&pi, &pi_sum, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
第一引数から、「和を取りたい変数」「和を受け取りたい変数」「変数の数」「変数の型」「やりたい演算」「コミュニケータ」の順番で指定する。ここでは一つの変数のみ総和演算を行っているが、配列を渡して一気に複数のデータについて総和を取ることもできる。また、総和だけでなく積や論理演算も実行できる。
-円周率の推定値piと、その自乗pi2 = pi*piについて総和を取り、定義通りに期待値と標準偏差を求めるコードがcalc_pi_reduce.cppである。
+円周率の推定値piと、その自乗pi2 = pi*piについて総和を取り、定義通りに期待値と標準偏差を求めるコードが以下のcalc_pi_reduce.cppである。
#include <cstdio>
#include <random>
#include <algorithm>
diff --git a/day4/README.md b/day4/README.md
index 74e3cad..2dc2f54 100644
--- a/day4/README.md
+++ b/day4/README.md
@@ -1,14 +1,17 @@
# Day 4 : 領域分割による非自明並列
+
+これまでは並列化として、自明並列を扱ってきた。自明並列はほとんど通信を必要とせず、並列化効率が高いため計算資源を最も有効に使える計算方法である。しかし、せっかくスパコンを使うのであるから、本格的に通信を伴う非自明並列に挑戦してみよう。
+
+
## 非自明並列
-Day 3では自明並列を扱ってきた。自明並列は別名「馬鹿パラ」と呼ばれ、馬鹿にされる傾向にあるのだが、並列化効率が高いため、「計算資源は」最も有効に使える計算方法である。さて、「スパコンはノードを束ねたもの」であり、「ノードとは本質的にはPCと同じもの」であることは既に述べた。しかし「普通のPCを多数束ねたらスパコンになるか」というとそうではなく、スパコンとして動作をするためには「ネットワーク」と「信頼性」が重要なのであった。実は、馬鹿パラは「ネットワーク」と「信頼性」のどちらも必要としない。
+「スパコンはノードを束ねたもの」であり、「ノードとは本質的にはPCと同じもの」であることは既に述べた。しかし「普通のPCを多数束ねたらスパコンになるか」というとそうではなく、スパコンとして動作をするためには「ネットワーク」と「信頼性」が重要である。これまで扱ってきた自明並列(通称「馬鹿パラ」)は、「ネットワーク」と「信頼性」のどちらも必要としない。
-パラメタ並列の場合、一番最初に「どのパラメタをどのプロセスが担当すべきか」をばらまくのに通信したあとは通信不要である(計算が終わったら結果をファイルに吐いてしまえばよい)。したがって、各ノードが高速なネットワークで接続されている必要はなく、たとえばイーサネットなどでつないでしまって全く問題ない。
-また、大規模な非自明並列計算を実行するには高い信頼性が求められるが、馬鹿パラは信頼性も要求しない。計算途中でノードが壊れてしまっても、そのノードでしていた計算だけやり直せばよいだけのことである。
-つまり馬鹿パラとは最も計算資源は有効に使えるものの、「ネットワーク」と「信頼性」という、スパコンの重要な特性を全く使わない計算方法なのであった。なので、主に馬鹿パラで計算する場合には、「普通のPCを多数束ねたPCクラスタ」で全く構わない。
+パラメタ並列の場合、一番最初に「どのパラメタをどのプロセスが担当すべきか」をばらまくのに通信したあとは通信不要である(計算が終わったら結果をファイルに吐いてしまえばよい)。したがって、各ノードが高速なネットワークで接続されている必要はなく、たとえばイーサネットなどでつないでしまって全く問題ない。また、大規模な非自明並列計算を実行するには高い信頼性が求められるが、馬鹿パラは信頼性も要求しない。計算途中でノードが壊れてしまっても、そのノードでしていた計算だけやり直せばよいだけのことである。
+つまり馬鹿パラとは最も計算資源は有効に使えるものの、「ネットワーク」と「信頼性」という、スパコンの重要な二大特性を全く使わない計算方法なのであった。なので、主に馬鹿パラで計算する場合には「普通のPCを多数束ねたPCクラスタ」で全く構わない。
さて、馬鹿パラであろうとなんであろうと、スパコンを活用していることにはかわりないし、それで良い科学的成果が出るのならそれで良いのだが、せっかくスパコンを使うのなら、もう少し「スパコンらしさ」を活用してみたい。というわけで、「ネットワーク」と「信頼性」をどちらも要求する **非自明並列 (non-trivial parallel)** に挑戦してみよう。
diff --git a/day4/index.html b/day4/index.html
index c86c480..927b5ea 100644
--- a/day4/index.html
+++ b/day4/index.html
@@ -73,13 +73,15 @@
Day 4 : 領域分割による非自明並列
+
+これまでは並列化として、自明並列を扱ってきた。自明並列はほとんど通信を必要とせず、並列化効率が高いため計算資源を最も有効に使える計算方法である。しかし、せっかくスパコンを使うのであるから、本格的に通信を伴う非自明並列に挑戦してみよう。
非自明並列
-Day 3では自明並列を扱ってきた。自明並列は別名「馬鹿パラ」と呼ばれ、馬鹿にされる傾向にあるのだが、並列化効率が高いため、「計算資源は」最も有効に使える計算方法である。さて、「スパコンはノードを束ねたもの」であり、「ノードとは本質的にはPCと同じもの」であることは既に述べた。しかし「普通のPCを多数束ねたらスパコンになるか」というとそうではなく、スパコンとして動作をするためには「ネットワーク」と「信頼性」が重要なのであった。実は、馬鹿パラは「ネットワーク」と「信頼性」のどちらも必要としない。
+「スパコンはノードを束ねたもの」であり、「ノードとは本質的にはPCと同じもの」であることは既に述べた。しかし「普通のPCを多数束ねたらスパコンになるか」というとそうではなく、スパコンとして動作をするためには「ネットワーク」と「信頼性」が重要である。これまで扱ってきた自明並列(通称「馬鹿パラ」)は、「ネットワーク」と「信頼性」のどちらも必要としない。
fig/bakapara.png
-パラメタ並列の場合、一番最初に「どのパラメタをどのプロセスが担当すべきか」をばらまくのに通信したあとは通信不要である(計算が終わったら結果をファイルに吐いてしまえばよい)。したがって、各ノードが高速なネットワークで接続されている必要はなく、たとえばイーサネットなどでつないでしまって全く問題ない。 また、大規模な非自明並列計算を実行するには高い信頼性が求められるが、馬鹿パラは信頼性も要求しない。計算途中でノードが壊れてしまっても、そのノードでしていた計算だけやり直せばよいだけのことである。 つまり馬鹿パラとは最も計算資源は有効に使えるものの、「ネットワーク」と「信頼性」という、スパコンの重要な特性を全く使わない計算方法なのであった。なので、主に馬鹿パラで計算する場合には、「普通のPCを多数束ねたPCクラスタ」で全く構わない。
+パラメタ並列の場合、一番最初に「どのパラメタをどのプロセスが担当すべきか」をばらまくのに通信したあとは通信不要である(計算が終わったら結果をファイルに吐いてしまえばよい)。したがって、各ノードが高速なネットワークで接続されている必要はなく、たとえばイーサネットなどでつないでしまって全く問題ない。また、大規模な非自明並列計算を実行するには高い信頼性が求められるが、馬鹿パラは信頼性も要求しない。計算途中でノードが壊れてしまっても、そのノードでしていた計算だけやり直せばよいだけのことである。 つまり馬鹿パラとは最も計算資源は有効に使えるものの、「ネットワーク」と「信頼性」という、スパコンの重要な二大特性を全く使わない計算方法なのであった。なので、主に馬鹿パラで計算する場合には「普通のPCを多数束ねたPCクラスタ」で全く構わない。
さて、馬鹿パラであろうとなんであろうと、スパコンを活用していることにはかわりないし、それで良い科学的成果が出るのならそれで良いのだが、せっかくスパコンを使うのなら、もう少し「スパコンらしさ」を活用してみたい。というわけで、「ネットワーク」と「信頼性」をどちらも要求する 非自明並列 (non-trivial parallel) に挑戦してみよう。
馬鹿パラではほとんど通信が発生しなかったのに対して、非自明並列は頻繁に通信が必要とする。 科学計算はなんらかの繰り返し計算(例えば時間発展)をすることが多いが、意味のある並列計算を行う場合、毎ステップ通信が必要となる。この時、「計算に関わる全ノードと毎回通信が発生する」タイプと、「論理的に距離が近いノードのみと通信が必要となる」タイプにわかれる。
@@ -132,8 +134,70 @@ 一次元拡散方程式 (シ
}
index++;
}
-あとは適当な条件を与えれば時間発展させることができる。ここでは、「一様加熱」と「温度固定」の二通りを試してみよう。コードはこちら。
-
+あとは適当な条件を与えれば時間発展させることができる。ここでは、「一様加熱」と「温度固定」の二通りを試してみよう。以下のコードをthermal.cppという名前で保存、実行してみよう。
+#include <cstdio>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+const int L = 128;
+const int STEP = 100000;
+const int DUMP = 1000;
+
+void onestep(std::vector<double> &lattice, const double h) {
+ static std::vector<double> orig(L);
+ std::copy(lattice.begin(), lattice.end(), orig.begin());
+ for (int i = 1; i < L - 1; i++) {
+ lattice[i] += h * (orig[i - 1] - 2.0 * orig[i] + orig[i + 1]);
+ }
+ // For Periodic Boundary
+ lattice[0] += h * (orig[L - 1] - 2.0 * lattice[0] + orig[1]);
+ lattice[L - 1] += h * (orig[L - 2] - 2.0 * lattice[L - 1] + orig[0]);
+}
+
+void dump(std::vector<double> &data) {
+ static int index = 0;
+ char filename[256];
+ sprintf(filename, "data%03d.dat", index);
+ std::cout << filename << std::endl;
+ std::ofstream ofs(filename);
+ for (int i = 0; i < data.size(); i++) {
+ ofs << i << " " << data[i] << std::endl;
+ }
+ index++;
+}
+
+void fixed_temperature(std::vector<double> &lattice) {
+ const double h = 0.01;
+ const double Q = 1.0;
+ for (int i = 0; i < STEP; i++) {
+ onestep(lattice, h);
+ lattice[L / 4] = Q;
+ lattice[3 * L / 4] = -Q;
+ if ((i % DUMP) == 0) dump(lattice);
+ }
+}
+
+void uniform_heating(std::vector<double> &lattice) {
+ const double h = 0.2;
+ const double Q = 1.0;
+ for (int i = 0; i < STEP; i++) {
+ onestep(lattice, h);
+ for (auto &s : lattice) {
+ s += Q * h;
+ }
+ lattice[0] = 0.0;
+ lattice[L - 1] = 0.0;
+ if ((i % DUMP) == 0) dump(lattice);
+ }
+}
+
+int main() {
+ std::vector<double> lattice(L, 0.0);
+ //uniform_heating(lattice);
+ fixed_temperature(lattice);
+}
+実行結果は以下のようになる。
fig/thermal.png
@@ -201,8 +265,7 @@ 一次元拡散方程式 (並列版)
「一度まとめてから吐く」いちど、ルートプロセス(ランク0番)に通信でデータを集めてしまってから、ルートプロセスが責任を持って一気に吐く。数千プロセスでも速度面で問題なくファイル出力できたが、全プロセスが保持する状態を一度一つのノードに集めるため、数万プロセス実行時にメモリ不足で落ちた。
とりあえずメモリに問題なければ「3. 一度まとめてから吐く」が楽なので、今回はこれを採用しよう。メモリが厳しかったり、数万プロセスの計算とかする時にはなにか工夫してくださいまし。
-さて、「一度まとめてから吐く」ためには、「各プロセスにバラバラにあるデータを、どこかのプロセスに一括して持ってくる」必要があるのだが、MPIには そのものずばりMPI_Gatherという関数がある。使い方は以下のサンプルを見たほうが早いと思う。
-
+さて、「一度まとめてから吐く」ためには、「各プロセスにバラバラにあるデータを、どこかのプロセスに一括して持ってくる」必要があるのだが、MPIには そのものずばりMPI_Gatherという関数がある。使い方はサンプルを見たほうが早い。以下をgather.cppという名前で保存、実行しよう。
#include <cstdio>
#include <mpi.h>
#include <vector>
@@ -314,9 +377,101 @@ 一次元拡散方程式 (並列版)
if ((i % DUMP) == 0) dump_mpi(lattice, rank, procs);
}
}
-これも一様加熱と同じで、「温度を固定している場所がどのプロセスが担当するどの場所か」を調べる必要があるが、それを考えるのはさほど難しくないだろう。
-そんなわけで完成した並列コードがこちら。
-
+これも一様加熱と同じで、「温度を固定している場所がどのプロセスが担当するどの場所か」を調べる必要があるが、それを考えるのはさほど難しくないだろう。そんなわけで完成した並列コードをthermal_mpi.cppという名前で保存しよう。
+#include <cstdio>
+#include <fstream>
+#include <iostream>
+#include <mpi.h>
+#include <vector>
+
+const int L = 128;
+const int STEP = 100000;
+const int DUMP = 1000;
+
+void dump(std::vector<double> &data) {
+ static int index = 0;
+ char filename[256];
+ sprintf(filename, "data%03d.dat", index);
+ std::cout << filename << std::endl;
+ std::ofstream ofs(filename);
+ for (unsigned int i = 0; i < data.size(); i++) {
+ ofs << i << " " << data[i] << std::endl;
+ }
+ index++;
+}
+
+void dump_mpi(std::vector<double> &local, int rank, int procs) {
+ static std::vector<double> global(L);
+ MPI_Gather(&(local[1]), L / procs, MPI_DOUBLE, global.data(), L / procs, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ if (rank == 0) {
+ dump(global);
+ }
+}
+
+void onestep(std::vector<double> &lattice, double h, int rank, int procs) {
+ const int size = lattice.size();
+ static std::vector<double> orig(size);
+ std::copy(lattice.begin(), lattice.end(), orig.begin());
+ // ここから通信のためのコード
+ const int left = (rank - 1 + procs) % procs; // 左のランク番号
+ const int right = (rank + 1) % procs; // 右のランク番号
+ MPI_Status st;
+ // 右端を右に送って、左端を左から受け取る
+ MPI_Sendrecv(&(lattice[size - 2]), 1, MPI_DOUBLE, right, 0, &(orig[0]), 1, MPI_DOUBLE, left, 0, MPI_COMM_WORLD, &st);
+ // 左端を左に送って、右端を右から受け取る
+ MPI_Sendrecv(&(lattice[1]), 1, MPI_DOUBLE, left, 0, &(orig[size - 1]), 1, MPI_DOUBLE, right, 0, MPI_COMM_WORLD, &st);
+
+ //あとはシリアル版と同じ
+ for (int i = 1; i < size - 1; i++) {
+ lattice[i] += h * (orig[i - 1] - 2.0 * orig[i] + orig[i + 1]);
+ }
+}
+
+void uniform_heating(std::vector<double> &lattice, int rank, int procs) {
+ const double h = 0.2;
+ const double Q = 1.0;
+ for (int i = 0; i < STEP; i++) {
+ onestep(lattice, h, rank, procs);
+ for (auto &s : lattice) {
+ s += Q * h;
+ }
+ if (rank == 0) {
+ lattice[1] = 0.0;
+ }
+ if (rank == procs - 1) {
+ lattice[lattice.size() - 2] = 0.0;
+ }
+ if ((i % DUMP) == 0) dump_mpi(lattice, rank, procs);
+ }
+}
+
+void fixed_temperature(std::vector<double> &lattice, int rank, int procs) {
+ const double h = 0.01;
+ const double Q = 1.0;
+ const int s = L / procs;
+ for (int i = 0; i < STEP; i++) {
+ onestep(lattice, h, rank, procs);
+ if (rank == (L / 4 / s)) {
+ lattice[L / 4 - rank * s + 1] = Q;
+ }
+ if (rank == (3 * L / 4 / s)) {
+ lattice[3 * L / 4 - rank * s + 1] = -Q;
+ }
+ if ((i % DUMP) == 0) dump_mpi(lattice, rank, procs);
+ }
+}
+
+int main(int argc, char **argv) {
+ MPI_Init(&argc, &argv);
+ int rank, procs;
+ MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+ MPI_Comm_size(MPI_COMM_WORLD, &procs);
+ const int mysize = L / procs + 2;
+ std::vector<double> local(mysize);
+ uniform_heating(local, rank, procs);
+ //fixed_temperature(local, rank, procs);
+ MPI_Finalize();
+}
せっかく並列化したので、高速化したかどうか調べてみよう。一様加熱の計算をさせてみる。
まずはシリアル版の速度。
$ clang++ -O3 -std=c++11 thermal.cpp
diff --git a/day5/README.md b/day5/README.md
index d90d49e..46d21f2 100644
--- a/day5/README.md
+++ b/day5/README.md
@@ -1,7 +1,13 @@
# Day 5 :二次元反応拡散方程式
-Day 4で一次元拡散方程式を領域分割により並列化した。後はこの応用で相互作用距離が短いモデルはなんでも領域分割できるのだが、二次元、三次元だと、一次元よりちょっと面倒くさい。後、熱伝導方程式は、「最終的になにかに落ち着く」方程式なので、シミュレーションしててあまりおもしろいものではない。そこで、二次元で、差分法で簡単に解けて、かつ結果がそこそこ面白い題材として反応拡散方程式(reaction-diffusion system)を取り上げる。反応拡散方程式とは、拡散方程式に力学系がくっついたような系で、様々なパターンを作る。例えば「reaction-diffusion system」でイメージ検索してみて欲しい。生物の模様なんかがこの方程式系で説明されたりする。
+
+Day 4では一次元拡散方程式を領域分割により並列化した。後はこの応用で相互作用距離が短いモデルはなんでも領域分割できるのだが、二次元、三次元だと、一次元よりちょっと面倒くさい。後、熱伝導方程式は、「最終的になにかに落ち着く」方程式なので、シミュレーションしててあまりおもしろいものではない。そこで、二次元で、差分法で簡単に解けて、かつ結果がそこそこ面白い題材として反応拡散方程式を取り上げる。
+
+
+## 反応拡散方程式
+
+反応拡散方程式(reaction-diffusion system)とは、拡散方程式に力学系がくっついたような系で、様々なパターンを作る。例えば「reaction-diffusion system」でイメージ検索してみて欲しい。生物の模様なんかがこの方程式系で説明されたりする。
世の中には様々な反応拡散方程式があるのだが、ここでは[Gray-Scottモデル](https://groups.csail.mit.edu/mac/projects/amorphous/GrayScott/)と呼ばれる、以下の方程式系を考えよう。
diff --git a/day5/index.html b/day5/index.html
index 2b319dc..7a982fa 100644
--- a/day5/index.html
+++ b/day5/index.html
@@ -73,7 +73,10 @@
Day 5 :二次元反応拡散方程式
-Day 4で一次元拡散方程式を領域分割により並列化した。後はこの応用で相互作用距離が短いモデルはなんでも領域分割できるのだが、二次元、三次元だと、一次元よりちょっと面倒くさい。後、熱伝導方程式は、「最終的になにかに落ち着く」方程式なので、シミュレーションしててあまりおもしろいものではない。そこで、二次元で、差分法で簡単に解けて、かつ結果がそこそこ面白い題材として反応拡散方程式(reaction-diffusion system)を取り上げる。反応拡散方程式とは、拡散方程式に力学系がくっついたような系で、様々なパターンを作る。例えば「reaction-diffusion system」でイメージ検索してみて欲しい。生物の模様なんかがこの方程式系で説明されたりする。
+
+Day 4では一次元拡散方程式を領域分割により並列化した。後はこの応用で相互作用距離が短いモデルはなんでも領域分割できるのだが、二次元、三次元だと、一次元よりちょっと面倒くさい。後、熱伝導方程式は、「最終的になにかに落ち着く」方程式なので、シミュレーションしててあまりおもしろいものではない。そこで、二次元で、差分法で簡単に解けて、かつ結果がそこそこ面白い題材として反応拡散方程式を取り上げる。
+反応拡散方程式
+反応拡散方程式(reaction-diffusion system)とは、拡散方程式に力学系がくっついたような系で、様々なパターンを作る。例えば「reaction-diffusion system」でイメージ検索してみて欲しい。生物の模様なんかがこの方程式系で説明されたりする。
世の中には様々な反応拡散方程式があるのだが、ここではGray-Scottモデルと呼ばれる、以下の方程式系を考えよう。
\[
\frac{\partial u}{\partial t} = D_u \Delta u + u^2 v - (F+k)u
@@ -152,8 +155,96 @@ シリアル版
}
先程述べたように、偶数時刻と奇数時刻で二本の配列を使い分けているのに注意。
save_as_datは、呼ばれるたびに配列を連番のファイル名で保存する関数である。
-全体のコードはこんな感じになる。
-
+全体のコードはこんな感じになる(gs.cpp)。
+#include <cstdio>
+#include <iostream>
+#include <vector>
+#include <fstream>
+
+const int L = 128;
+const int TOTAL_STEP = 20000;
+const int INTERVAL = 200;
+const double F = 0.04;
+const double k = 0.06075;
+const double dt = 0.2;
+const double Du = 0.05;
+const double Dv = 0.1;
+
+typedef std::vector<double> vd;
+
+void init(vd &u, vd &v) {
+ int d = 3;
+ for (int i = L / 2 - d; i < L / 2 + d; i++) {
+ for (int j = L / 2 - d; j < L / 2 + d; j++) {
+ u[j + i * L] = 0.7;
+ }
+ }
+ d = 6;
+ for (int i = L / 2 - d; i < L / 2 + d; i++) {
+ for (int j = L / 2 - d; j < L / 2 + d; j++) {
+ v[j + i * L] = 0.9;
+ }
+ }
+}
+
+double calcU(double tu, double tv) {
+ return tu * tu * tv - (F + k) * tu;
+}
+
+double calcV(double tu, double tv) {
+ return -tu * tu * tv + F * (1.0 - tv);
+}
+
+double laplacian(int ix, int iy, vd &s) {
+ double ts = 0.0;
+ ts += s[ix - 1 + iy * L];
+ ts += s[ix + 1 + iy * L];
+ ts += s[ix + (iy - 1) * L];
+ ts += s[ix + (iy + 1) * L];
+ ts -= 4.0 * s[ix + iy * L];
+ return ts;
+}
+
+void calc(vd &u, vd &v, vd &u2, vd &v2) {
+ for (int iy = 1; iy < L - 1; iy++) {
+ for (int ix = 1; ix < L - 1; ix++) {
+ double du = 0;
+ double dv = 0;
+ const int i = ix + iy * L;
+ du = Du * laplacian(ix, iy, u);
+ dv = Dv * laplacian(ix, iy, v);
+ du += calcU(u[i], v[i]);
+ dv += calcV(u[i], v[i]);
+ u2[i] = u[i] + du * dt;
+ v2[i] = v[i] + dv * dt;
+ }
+ }
+}
+
+void save_as_dat(vd &u) {
+ static int index = 0;
+ char filename[256];
+ sprintf(filename, "conf%03d.dat", index);
+ std::cout << filename << std::endl;
+ std::ofstream ofs(filename, std::ios::binary);
+ ofs.write((char *)(u.data()), sizeof(double)*L * L);
+ index++;
+}
+
+int main() {
+ const int V = L * L;
+ vd u(V, 0.0), v(V, 0.0);
+ vd u2(V, 0.0), v2(V, 0.0);
+ init(u, v);
+ for (int i = 0; i < TOTAL_STEP; i++) {
+ if (i & 1) {
+ calc(u2, v2, u, v);
+ } else {
+ calc(u, v, u2, v2);
+ }
+ if (i % INTERVAL == 0) save_as_dat(u);
+ }
+}
コンパイル、実行してみよう。
$ g++ -O3 gs.cpp
$ time ./a.out
@@ -457,8 +548,8 @@ 並列化ステップ2: デ
}
}
送信前や送信後にデータの処理が必要となるので、やってることが単純なわりにコード量がそこそこの長さになる。 このあたりが「MPIは面倒くさい」と言われる所以かもしれない。筆者も「MPIは面倒くさい」ことは否定しない。 しかし、ここまで読んでくださった方なら「MPIは難しくはない」ということも同意してもらえると思う。 MPIは書いた通りに動く。なので、通信アルゴリズムが決まっていれば、その手順どおりに書くだけである。 実際面倒なのは通信そのものよりも、通信の前処理と後処理だったりする(そもそも今回も通信は一行だけだ)。
-以上をすべてまとめたコードは以下の通り。
-
+以上をすべてまとめたコードをgather2d.cppとしよう。やや大きいので、ウェブへのリンクを貼っておく。
+https://github.com/kaityo256/sevendayshpc/blob/master/day5/gather2d.cpp
main関数だけ書いておくとこんな感じ。
int main(int argc, char **argv) {
MPI_Init(&argc, &argv);
@@ -536,7 +627,7 @@ 並列化ステップ2:
}
全く同様にy方向の通信も書けるが、先に述べたように「左右からもらったデータも転送」するため、その分がちょっとだけ異なる。
このアルゴリズムを実装するとこんな感じになる。
-
+https://github.com/kaityo256/sevendayshpc/blob/master/day5/sendrecv.cpp
実行結果はこんな感じ。
$ mpic++ sendrecv.cpp
$ mpirun -np 4 ./a.out
diff --git a/day6/README.md b/day6/README.md
index e313b51..d33371e 100644
--- a/day6/README.md
+++ b/day6/README.md
@@ -1,7 +1,6 @@
# Day 6 : ハイブリッド並列
-## ハイブリッド並列とは
-
+
これまで、並列化の手段としてMPIを使った「プロセス並列」を行ってきた。
最初に述べたように、並列化には他にも「スレッド並列」という手段がある。
プロセス並列が分散メモリ型、スレッド並列が共有メモリ型であり、
@@ -13,10 +12,11 @@
面倒になるので、できることならやりたくない。しかし、アプリケーションや
サイズによっては、ハイブリッド並列を選択せざるを得ない場合もあるだろう。
ここでは、スレッド並列を行うときの注意点や、ハイブリッド並列の実例について見てみよう。
+
## 仮想メモリとTLB
-さて、プロセス並列ではあまり気にしなくてよかったが、スレッド並列を行う時には
+プロセス並列ではあまり気にしなくてよかったが、スレッド並列を行う時には
気にしなければいけないものとして「NUMA」というものがある。
「NUMA」を気にするためには、仮想メモリについて知らないといけない。
というわけで、仮想メモリについて見てみよう。
diff --git a/day6/index.html b/day6/index.html
index 2ecc231..7d7d1cb 100644
--- a/day6/index.html
+++ b/day6/index.html
@@ -72,10 +72,10 @@
Day 6 : ハイブリッド並列
-ハイブリッド並列とは
-これまで、並列化の手段としてMPIを使った「プロセス並列」を行ってきた。 最初に述べたように、並列化には他にも「スレッド並列」という手段がある。 プロセス並列が分散メモリ型、スレッド並列が共有メモリ型であり、 スレッド並列だけではノードをまたぐことができないので、普通「スパコンを使う」 というとプロセス並列が必須になる。 さて、MPIを使ったプロセス並列「だけ」による並列化を「flat-MPI」と呼ぶ。 一方、プロセス並列とスレッド並列を併用する並列化を「ハイブリッド並列」と呼ぶ。 当然のことながら、ハイブリッド並列は、プロセス並列単体、スレッド並列単体よりも 面倒になるので、できることならやりたくない。しかし、アプリケーションや サイズによっては、ハイブリッド並列を選択せざるを得ない場合もあるだろう。 ここでは、スレッド並列を行うときの注意点や、ハイブリッド並列の実例について見てみよう。
+
+これまで、並列化の手段としてMPIを使った「プロセス並列」を行ってきた。 最初に述べたように、並列化には他にも「スレッド並列」という手段がある。 プロセス並列が分散メモリ型、スレッド並列が共有メモリ型であり、 スレッド並列だけではノードをまたぐことができないので、普通「スパコンを使う」 というとプロセス並列が必須になる。 さて、MPIを使ったプロセス並列「だけ」による並列化を「flat-MPI」と呼ぶ。 一方、プロセス並列とスレッド並列を併用する並列化を「ハイブリッド並列」と呼ぶ。 当然のことながら、ハイブリッド並列は、プロセス並列単体、スレッド並列単体よりも 面倒になるので、できることならやりたくない。しかし、アプリケーションや サイズによっては、ハイブリッド並列を選択せざるを得ない場合もあるだろう。 ここでは、スレッド並列を行うときの注意点や、ハイブリッド並列の実例について見てみよう。
仮想メモリとTLB
-さて、プロセス並列ではあまり気にしなくてよかったが、スレッド並列を行う時には 気にしなければいけないものとして「NUMA」というものがある。 「NUMA」を気にするためには、仮想メモリについて知らないといけない。 というわけで、仮想メモリについて見てみよう。
+プロセス並列ではあまり気にしなくてよかったが、スレッド並列を行う時には 気にしなければいけないものとして「NUMA」というものがある。 「NUMA」を気にするためには、仮想メモリについて知らないといけない。 というわけで、仮想メモリについて見てみよう。
OSは実に様々なことをやっているが、特に重要な仕事に「メモリ管理」がある。 物理的には「メモリ」はマザーボードに刺さったDRAMを指すが、 OSの管理下で動くプロセスから見える「メモリ」は、それを仮想化したものである。 プロセスにとっては連続に見えるメモリも、実はDRAM上にバラバラに割り付けられて いるかもしれない。OSは、「プロセスから見えるアドレス」と「物理的にDRAMに割り当てられたアドレス」を うまいこと変換して、プロセスが物理メモリを意識しないで済むようにしている。 このような仕組みを「仮想メモリ (virtual memory)」と呼ぶ。 仮想メモリを扱う利点としては、
- OSがメモリを管理してくれるので複数のプロセスがお互いのメモリを気にしなくて良くなる(セキュリティ上も好ましい)
@@ -83,8 +83,7 @@ 仮想メモリとTLB
- メモリが足りない時にハードディスクなどにスワップすることで、物理メモリより大きな論理メモリ空間がとれる
などが挙げられる。なお、Windowsでは「ハードディスクにスワップする領域の上限」のことを「仮想メモリ」と呼んでいるようなので注意。
-実際に、プロセスごとに固有の仮想メモリが与えられているのを見てみよう。こんなコードを書いてみる。
-
+実際に、プロセスごとに固有の仮想メモリが与えられているのを見てみよう。こんなコード(vmem.cpp)を書いてみる。
#include <cstdio>
#include <mpi.h>
@@ -182,7 +181,9 @@ OpenMPの例
- Intel(R) Xeon(R) CPU E5-2680 v3 @ 2.50GHz 12コア x 2ソケット
-まず、シリアルコードとしてDay 4で使ったGray Scottモデルの計算を使おう。純粋に計算のみをカウントするため、途中のファイル出力を削除し、また実行時間を測定するようにしたのがgs.cppである。ただし、デバッグのために最終結果だけファイルに出力している。コンパイルしてperfでプロファイルをとってみよう。まず、perf recordで記録を取る。
+まず、シリアルコードとしてDay 4で使ったGray Scottモデルの計算を使おう。純粋に計算のみをカウントするため、途中のファイル出力を削除し、また実行時間を測定するようにしたものがgs.cppである。
+https://github.com/kaityo256/sevendayshpc/blob/master/day6/gs.cpp
+ただし、デバッグのために最終結果だけファイルに出力している。コンパイルしてperfでプロファイルをとってみよう。まず、perf recordで記録を取る。
$ g++ -O3 -mavx2 -std=c++11 -fopenmp gs.cpp -o gs.out
$ perf record ./gs.out
2527 [ms]
@@ -217,7 +218,7 @@ OpenMPの例
}
二重ループになっている。OpenMPは、並列実行したいループの直前にディレクティブを入れて、「このループを並列化してください」と指示することで並列化する。スレッド並列する時には、ループインデックス間に依存性がないか確認しなければならないのだが、今回はたまたまループインデックス間に全く依存関係がないので、好きなように並列化してよい(たまたまというか、そうなるように題材を選んだわけだが)。
まずは内側のループにディレクティブを入れてみよう。#pragma omp parallel forというディレクティブを対象ループの直前に入れるだけでよい。
-
+https://github.com/kaityo256/sevendayshpc/blob/master/day6/gs_omp1.cpp
void calc(vd &u, vd &v, vd &u2, vd &v2) {
for (int iy = 1; iy < L - 1; iy++) {
#pragma omp parallel for
@@ -243,7 +244,7 @@ OpenMPの例
diff conf000.org conf000.dat
問題なさそうですね。
次に、外側を並列化してみよう。
-
+https://github.com/kaityo256/sevendayshpc/blob/master/day6/gs_omp2.cpp
void calc(vd &u, vd &v, vd &u2, vd &v2) {
#pragma omp parallel for
for (int iy = 1; iy < L - 1; iy++) {
@@ -412,7 +413,9 @@ ハイブリッド並列の実例
// なにか時間を計測したい処理
const auto e = std::chrono::system_clock::now();
const auto elapsed = std::chrono::duration_cast<std::chrono::milliseconds>(e - s).count();
-これでelapsedにミリ秒単位の値が入る。このようにして作ったハイブリッド版の反応拡散方程式ソルバがgs_hybrid.cppである。 筆者の環境ではMPIにパスが通してあるので、以下のようなオプションでコンパイルできる。
+これでelapsedにミリ秒単位の値が入る。このようにして作ったハイブリッド版の反応拡散方程式ソルバがgs_hybrid.cppである。
+https://github.com/kaityo256/sevendayshpc/blob/master/day6/gs_hybrid.cpp
+筆者の環境ではMPIにパスが通してあるので、以下のようなオプションでコンパイルできる。
g++ -fopenmp -O3 -mavx2 gs_hybrid.cpp -lmpi -lmpi_cxx
手元のMacで2プロセス x 2スレッドで実行すると以下のような結果を吐く。
$ OMP_NUM_THREADS=2 mpiexec -np 2 ./a.out
diff --git a/day7/README.md b/day7/README.md
index 33a28ef..96c0d17 100644
--- a/day7/README.md
+++ b/day7/README.md
@@ -1,9 +1,8 @@
# Day 7 : SIMD化
-## はじめに
-
-ここまで読んだ人、お疲れ様です。ここから読んでいる人、それでも問題ありません。
-これまで、主に並列化についてだらだら書いてきたが、最後はシングルコアでの最適化技術であるSIMD化について説明してみたいと思う。
+
+ここまで読んだ人、お疲れ様です。ここから読んでいる人、問題ありません。スパコンにおける並列化には、プロセス並列、スレッド並列、データ並列の三種類がある。これまでプロセス並列、スレッド並列について述べたので、最後はデータ並列であるSIMD化について紹介しよう。
+
## SIMDとは
diff --git a/day7/index.html b/day7/index.html
index d290c89..214feb9 100644
--- a/day7/index.html
+++ b/day7/index.html
@@ -73,8 +73,8 @@
Day 7 : SIMD化
-はじめに
-ここまで読んだ人、お疲れ様です。ここから読んでいる人、それでも問題ありません。 これまで、主に並列化についてだらだら書いてきたが、最後はシングルコアでの最適化技術であるSIMD化について説明してみたいと思う。
+
+ここまで読んだ人、お疲れ様です。ここから読んでいる人、問題ありません。スパコンにおける並列化には、プロセス並列、スレッド並列、データ並列の三種類がある。これまでプロセス並列、スレッド並列について述べたので、最後はデータ並列であるSIMD化について紹介しよう。
SIMDとは
スパコンプログラミングに興味があるような人なら、「SIMD」という言葉を聞いたことがあるだろう。SIMDとは、「single instruction multiple data」の略で、「一回の命令で複数のデータを同時に扱う」という意味である。先に、並列化は大きく分けて「データ並列」「共有メモリ並列」「分散メモリ並列」の三種類になると書いたが、SIMDはデータ並列(Data parallelism)に属す。現在、一般的に数値計算に使われるCPUにはほとんどSIMD命令が実装されている。後述するが、SIMDとは1サイクルに複数の演算を同時に行う技術であり、CPUの「理論ピーク性能」は、SIMDの能力を使い切った場合の性能を指す。したがって、まったくSIMD化できなければ、ピーク性能が数分の1になることと等価である。ここでは、なぜSIMDが必要になるか、そしてSIMDとは何かについて見てみよう。
計算機というのは、要するにメモリからデータと命令を取ってきて、演算器に投げ、結果をメモリに書き戻す機械である。CPUの動作単位は「サイクル」で表される。演算器に計算を投げてから、結果が返ってくるまでに数サイクルかかるが、現代のCPUではパイプライン処理という手法によって事実上1サイクルに1個演算ができる。1サイクル1演算できるので、あとは「1秒あたりのサイクル数=動作周波数」を増やせば増やすほど性能が向上することになる。
@@ -97,8 +97,7 @@ SIMDレジスタを触ってみるvoid print256d(__m256d x) {
printf("%f %f %f %f\n", x[3], x[2], x[1], x[0]);
}
-_m256d xが、そのままdouble x[4]として使えているのがわかると思う。この時、x[0]が一番下位となる。 先程の代入と合わせるとこんな感じになる。
-
+_m256d xが、そのままdouble x[4]として使えているのがわかると思う。この時、x[0]が一番下位となる。 先程の代入と合わせるとこんな感じになる(print.cpp)。
#include <cstdio>
#include <x86intrin.h>
@@ -130,8 +129,7 @@ SIMDレジスタを触ってみるret
vaddpdはSIMDの足し算を行う命令であり、ちゃんとYMMレジスタの足し算が呼ばれていることがわかる。
-実際に4要素同時に足し算できることを確認しよう。
-
+実際に4要素同時に足し算できることを確認しよう(add.cpp)。
#include <cstdio>
#include <x86intrin.h>
@@ -148,8 +146,7 @@ SIMDレジスタを触ってみる$ g++ -mavx2 add.cpp
$ ./a.out
10.000000 8.000000 6.000000 4.000000
-(0,1,2,3)というベクトルと、(4,5,6,7)というベクトルの和をとり、(4,6,8,10)というベクトルが得られた。 このように、ベクトル同士の演算に見えるので、SIMD化のことをベクトル化と呼んだりする。ただし、線形代数で出てくる ベクトルの積とは違い、SIMDの積は単に要素ごとの積になることに注意。実際、さっきの和を積にするとこうなる。
-
+(0,1,2,3)というベクトルと、(4,5,6,7)というベクトルの和をとり、(4,6,8,10)というベクトルが得られた。 このように、ベクトル同士の演算に見えるので、SIMD化のことをベクトル化と呼んだりする。ただし、線形代数で出てくる ベクトルの積とは違い、SIMDの積は単に要素ごとの積になることに注意。実際、さっきの和を積にするとこうなる(mul.cpp)。
#include <cstdio>
#include <x86intrin.h>
@@ -167,8 +164,7 @@ SIMDレジスタを触ってみる./a.out
21.000000 12.000000 5.000000 0.000000
それぞれ、0*0、1*5、2*6、3*7が計算されていることがわかる。
-あとSIMD化で大事なのは、SIMDレジスタへのデータの読み書きである。先程はデバッグのために_mm256_set_pdを使ったが、これは極めて遅い。 どんな動作をするか見てみよう。
-
+あとSIMD化で大事なのは、SIMDレジスタへのデータの読み書きである。先程はデバッグのために_mm256_set_pdを使ったが、これは極めて遅い。どんな動作をするか見てみよう(setpd.cpp)。
#include <x86intrin.h>
__m256d setpd(double a, double b, double c, double d) {
@@ -189,8 +185,7 @@ SIMDレジスタを触ってみる
ということをしている。ここで、xmm0レジスタとymm0レジスタの下位128ビットは共有していることに注意。 つまり、xmm0レジスタに読み書きすると、ymm0レジスタの下位128ビットも影響を受ける。 上記の例はそれを利用して、最終的に欲しい情報、つまり4要素をパックしたレジスタを作っている。
とりあえず4つの要素をYMMレジスタに載せることができれば、あとは4要素同時に計算ができるようになるのだが、 4要素をパックする際に_mm256_set_pdを使うとメモリアクセスが多くなって性能が出ない。 そのため、メモリから連続するデータをごそっとレジスタにとってきたり、書き戻したりする命令がある。 例えば、_mm256_load_pdは、指定されたポインタから連続する4つの倍精度実数をとってきて YMMレジスタに入れてくれる。ただし、そのポインタの指すアドレスは32バイトアラインされていなければならない。
-利用例はこんな感じになる。
-
+利用例はこんな感じになる(load.cpp)。
#include <cstdio>
#include <x86intrin.h>
@@ -211,8 +206,7 @@ SIMDレジスタを触ってみる$ g++ -mavx2 load.cpp
$ ./a.out
10.000000 8.000000 6.000000 4.000000
-_mm256_load_pdが何をやっているか(どんなアセンブリに対応するか)も見てみよう。こんなコードのアセンブリを見てみる。
-
+_mm256_load_pdが何をやっているか(どんなアセンブリに対応するか)も見てみよう。こんなコードのアセンブリを見てみる(loadasm.cpp)。
#include <x86intrin.h>
__m256d load(double *a, int index) {
return _mm256_load_pd(a + index);
@@ -232,8 +226,7 @@ 余談：アセンブ
繰り返しになるが、現在は「アセンブリ言語 (assembly language)」の方が一般的な用語であると思われるので、「アセンブリ言語をアセンブラがアセンブルして機械語にする」と表現することになんの問題もない。 しかし、誰かが「アセンブラを書く」もしくは「アセンブラ言語」と言ったときに、脊髄反射で「アセンブリ言語が正しい」とマウントを取る前に、上記のような事情を思い出していただけたらと思う。
余談の余談となるが、アセンブリで書かれたものを手で機械語に翻訳する作業を「ハンドアセンブル」と呼ぶ。昔のアセンブリはほぼ機械語と一対一対応しており、「便利なマクロ付き機械語」といった趣であったため、ハンドアセンブルはさほど難しい作業ではなかった。しかし、現在の機械語、特にx86の機械語はかなり複雑になっており、アセンブリから機械語に翻訳するのはかなり大変になっている。そのあたりは例えばx86_64機械語入門なんかを参照してほしい。
簡単なSIMD化の例
-では、実際にSIMD化をやってみよう。こんなコードを考える。 一次元の配列の単純な和のループである。
-
+では、実際にSIMD化をやってみよう。こんなコードを考える。 一次元の配列の単純な和のループである(func.cpp)。
const int N = 10000;
double a[N], b[N], c[N];
@@ -263,8 +256,7 @@ 簡単なSIMD化の例
二つのレジスタを足す
結果のレジスタを配列cのしかるべき場所に保存する
-ということをすればSIMD化完了である。コードを見たほうが早いと思う。
-
+ということをすればSIMD化完了である。コードを見たほうが早いと思う(func_simd.cpp)。
#include <x86intrin.h>
void func_simd() {
for (int i = 0; i < N; i += 4) {
@@ -320,8 +312,61 @@ 簡単なSIMD化の例
printf("%s is NG\n", type);
}
}
-全部まとめたコードはこちら。
-
+全部まとめたコードをsimdcheck.cppとして保存、実行してみよう。
+#include <cstdio>
+#include <random>
+#include <algorithm>
+#include <x86intrin.h>
+
+const int N = 10000;
+
+double a[N], b[N], c[N];
+double ans[N];
+
+void check(void(*pfunc)(), const char *type) {
+ pfunc();
+ unsigned char *x = (unsigned char *)c;
+ unsigned char *y = (unsigned char *)ans;
+ bool valid = true;
+ for (int i = 0; i < 8 * N; i++) {
+ if (x[i] != y[i]) {
+ valid = false;
+ break;
+ }
+ }
+ if (valid) {
+ printf("%s is OK\n", type);
+ } else {
+ printf("%s is NG\n", type);
+ }
+}
+
+void func() {
+ for (int i = 0; i < N; i++) {
+ c[i] = a[i] + b[i];
+ }
+}
+
+void func_simd() {
+ for (int i = 0; i < N; i += 4) {
+ __m256d va = _mm256_load_pd(&(a[i]));
+ __m256d vb = _mm256_load_pd(&(b[i]));
+ __m256d vc = va + vb;
+ _mm256_store_pd(&(c[i]), vc);
+ }
+}
+
+int main() {
+ std::mt19937 mt;
+ std::uniform_real_distribution<double> ud(0.0, 1.0);
+ for (int i = 0; i < N; i++) {
+ a[i] = ud(mt);
+ b[i] = ud(mt);
+ ans[i] = a[i] + b[i];
+ }
+ check(func, "scalar");
+ check(func_simd, "vector");
+}
実際に実行してテストしてみよう。
$ g++ -mavx2 -O3 simdcheck.cpp
$ ./a.out
@@ -565,7 +610,7 @@ もう少し実戦的なSIMD化
std::cout << r[i].z << std::endl;
}
}
-時間発展後に座標をダンプして、その結果を比較しよう。シリアル版をmag.cpp、SIMD版をmag_simd.cppとしておき、以下のようにコンパイル、実行、結果の比較をする。
+時間発展後に座標をダンプして、その結果を比較しよう。シリアル版をmag.cpp、SIMD版をmag_simd.cppとしておき、以下のようにコンパイル、実行、結果の比較をする。
$ g++ -std=c++11 -O3 -mavx2 -mfma mag.cpp -o a.out
$ g++ -std=c++11 -O3 -mavx2 -mfma mag_simd.cpp -o b.out
$ time ./a.out > a.txt
@@ -668,7 +713,7 @@ もう少し実戦的なSIMD化
cmpq %rcx, %rax
jne L13
vpermpdがシャッフル命令である。ループボディがかなり小さいが、このループは100000回まわるため、25000回しかまわらないコンパイラによる自動SIMD化ルーチンには勝てない。大雑把な話、ループボディの計算コストが半分だが、回転数が4倍なので2倍負けた、という感じである。
-上記の例のように、「いま手元にあるコード」をがんばって「そのままSIMD化」して高速化しても、データ構造を変えるとコンパイラがあっさり自動SIMD化できて負けることがある。多くの場合「SIMD化」はデータ構造のグローバルな変更を伴う。先のコードのAoS版であるmd.cppと、SoA版であるmd_soa.cppは、全く同じことをしているが 全書き換え になっている。今回はコードが短いから良いが、10万行とかあるコードだと「やっぱりSoAの方が早いから全書き換えで！」と気軽には言えないだろう。また、デバイスによってデータ構造の向き不向きもある。例えば「CPUではAoSの方が早いが、GPGPUではSoAの方が早い」なんてこともざらにある。こういう場合には、ホットスポットルーチンに入る前にAoS←→SoAの相互変換をしたりすることも検討するが、もちろんその分オーバーヘッドもあるので面倒くさいところである。
+上記の例のように、「いま手元にあるコード」をがんばって「そのままSIMD化」して高速化しても、データ構造を変えるとコンパイラがあっさり自動SIMD化できて負けることがある。多くの場合「SIMD化」はデータ構造のグローバルな変更を伴う。先のコードのAoS版であるmag.cppと、SoA版であるmag_soa.cppは、全く同じことをしているが 全書き換え になっている。今回はコードが短いから良いが、10万行とかあるコードだと「やっぱりSoAの方が早いから全書き換えで！」と気軽には言えないだろう。また、デバイスによってデータ構造の向き不向きもある。例えば「CPUではAoSの方が早いが、GPGPUではSoAの方が早い」なんてこともざらにある。こういう場合には、ホットスポットルーチンに入る前にAoS←→SoAの相互変換をしたりすることも検討するが、もちろんその分オーバーヘッドもあるので面倒くさいところである。
まぁ、以上のようにいろいろ面倒なことを書いたが、ちゃんと手を動かして上記を試してみた方には「SIMD化は(原理的には)簡単だ」ということには同意してもらえると思う。MPIもSIMD化も同じである。いろいろ考えることがあって面倒だが、やること自体は単純なので難しくはない。今回はシャッフル命令を取り上げたが、他にもマスク処理やgather/scatter、pack/unpackなど、SIMDには実に様々な命令がある。しかし、「そういう命令欲しいな」と思って調べたらたいがいある。あとは対応する組み込み関数を呼べばよい。要するに「やるだけ」である。ただし、MPI化は「やれば並列計算ができ、かつプロセスあたりの計算量を増やせばいくらでも並列化効率を上げられる」ことが期待されるのに対して、SIMD化は「やっても性能が向上するかはわからず、下手に手を出すよりコンパイラに任せた方が早い」なんてこともある。全くSIMD化されていないコードに対してSIMD化で得られるゲインは、256bitなら4倍、512ビットでも8倍程度しかなく、現実にはその半分も出れば御の字であろう。SIMD化はやってて楽しい作業であるが、手間とコストが釣り合うかどうかは微妙だな、というのが筆者の実感である。
diff --git a/makefile b/pdf.mk
similarity index 100%
rename from makefile
rename to pdf.mk
diff --git a/review/config.yml b/review/config.yml
index 0d33197..c258504 100644
--- a/review/config.yml
+++ b/review/config.yml
@@ -79,7 +79,7 @@ date: 2020-02-29
# 複数指定する場合は次のように記述する
# [["初版第1刷の日付", "初版第2刷の日付"], ["第2版第1刷の日付"]]
# 日付の後ろを空白文字で区切り、任意の文字列を置くことも可能。
-history: [["2020-02-25 ver 1.1"]]
+history: [["2020-03-03 ver 1.2"]]
# [experimental] 新刊を頒布したイベント名（例：「技術書典6（2019年春）新刊」）
pubevent_name:
# 権利表記(配列で複数指定可)
diff --git a/review/pre.rb b/review/pre.rb
index 6be6556..d314b52 100644
--- a/review/pre.rb
+++ b/review/pre.rb
@@ -1,6 +1,5 @@
-
def escape_underscore(str)
- str.gsub('_','@')
+ str.gsub("_", "@")
end
def escape_inline_math(str)
@@ -11,8 +10,18 @@ def escape_inline_math(str)
str
end
+def replace_review_command(line)
+ return line if line !~ /^$/
+
+ key = $1.strip
+ return "//}" if key == "end"
+
+ line = "//#{key}\{"
+ line
+end
+
def in_math
- while line=gets
+ while (line = gets)
if line=~/\$\$/
puts "//}"
return
@@ -22,11 +31,12 @@ def in_math
end
end
-while line=gets
+while (line = gets)
if line=~/\$\$/
puts "//texequation{"
in_math
else
+ line = replace_review_command(line)
puts escape_inline_math(line)
end
end
diff --git a/sevendayshpc.pdf b/sevendayshpc.pdf
index 5c4e189..0f82afc 100644
Binary files a/sevendayshpc.pdf and b/sevendayshpc.pdf differ