Merge pull request #10 from kaityo256/bugfix/7

fixes #7

day1/README.md

@@ -123,7 +123,7 @@ warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 
 ## はじめてのMPI
 
-環境構築ができたら、こんなコードを書いて、[hello.cpp](hello.cpp)という名前で保存してみよう。
+環境構築ができたら、こんなコードを書いて、`hello.cpp`という名前で保存してみよう。
 
 ```cpp
 #include <cstdio>
@@ -196,9 +196,7 @@ MPI_Comm_rank(MPI_COMM_WORLD, &rank);
 ```
 
 これを実行すると変数`rank`にランク番号が入る。N並列している場合、ランクは0からN-1までである。
-試してみよう。
-
-[rank.cpp](rank.cpp)
+試してみよう。以下を`rank.cpp`として保存し、コンパイル、実行してみよう。
 
 ```cpp
 #include <cstdio>
@@ -301,7 +299,7 @@ mpirun -np 4 ./a.out
 競合が起きないようにする。
 「ひとかたまりの処理」とは、例えば「`printf`で出力を始めてから終わるまで」である。
 
-例えば先程の[rank.cpp](rank.cpp)の例では、
+例えば先程の`rank.cpp`の例では、
 
 ```cpp
 printf("Hello! My rank is %d\n", rank);
@@ -330,11 +328,7 @@ puts("Hello! My rank is 3");
 
 とかになるだけで、さほど表示は乱れない。
 
-さて、同様なプログラムを`std::cout`で書いてみよう。
-
-こんな感じになると思う。
-
-[rank_stream.cpp](rank_stream.cpp)
+さて、同様なプログラムを`std::cout`で書いてみよう。以下を`rank_stream.cpp`という名前で保存、コンパイル、実行してみよう。
 
 ```cpp
 #include <iostream>
@@ -468,11 +462,7 @@ gdbで特定のプロセスにアタッチする。しかし、gdbでアタッ
 * gdbで変数をいじって無限ループを脱出させる
 * あとは好きなようにデバッグする
 
-という方針でいく。なお、なぜかMac OSではMPIプロセスへのgdbのアタッチがうまくいかなかったので、以下はCentOSで実行している。
-
-こんなコードを書く。
-
-[gdb_mpi.cpp](gdb_mpi.cpp)
+という方針でいく。なお、なぜかMac OSではMPIプロセスへのgdbのアタッチがうまくいかなかったので、以下はCentOSで実行している。以下を`gdb_mpi.cpp`という名前で保存しよう。
 
 ```cpp
 #include <cstdio>

day3/README.md

@@ -19,7 +19,7 @@
 まず、自明並列でよく出てくる例として、サンプリング数を並列化で稼ぐ方法を見てみよう。とりあえず定番の、
 モンテカルロ法で円周率を計算してみる。
 
-こんなコードを書いて、[calc_pi.cpp](calc_pi.cpp)という名前で保存してみよう。
+こんなコードを書いて、`calc_pi.cpp`という名前で保存してみよう。
 
 ```cpp
 #include <cstdio>
@@ -65,7 +65,7 @@ $ ./a.out
 4. ランク番号を乱数の種に使う
 5. そのまま`calc_pi`を呼ぶ。
 
-以上の修正をしたコードを[calc_pi_mpi.cpp](calc_pi_mpi.cpp)という名前で作成する。
+以上の修正をしたコードを`calc_pi_mpi.cpp`という名前で作成する。
 
 ```cpp
 #include <cstdio>
@@ -154,7 +154,7 @@ PID    COMMAND      %CPU TIME     #TH   #WQ  #PORT MEM    PURG   CMPRS  PGRP
 
 ## 自明並列テンプレート
 
-先程の並列プログラム[calc_pi_mpi.cpp](calc_pi_mpi.cpp)のmain関数はこうなっていた。
+先程の並列プログラム`calc_pi_mpi.cpp`のmain関数はこうなっていた。
 
 ```cpp
 int main(int argc, char **argv) {
@@ -242,9 +242,7 @@ MPI_Comm_size(MPI_COMM_WORLD, &procs)
 ```
 
 とすれば、`procs`にプロセス数が入る。これを使うと、先程のコードは
-こんな感じにかける。
-
-[processfiles.cpp](processfiles.cpp)
+こんな感じにかける(`processfiles.cpp`)。
 
 ```cpp
 #include <cstdio>
@@ -299,7 +297,7 @@ MPI_Allreduce(&pi, &pi_sum, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
 
 第一引数から、「和を取りたい変数」「和を受け取りたい変数」「変数の数」「変数の型」「やりたい演算」「コミュニケータ」の順番で指定する。ここでは一つの変数のみ総和演算を行っているが、配列を渡して一気に複数のデータについて総和を取ることもできる。また、総和だけでなく積や論理演算も実行できる。
 
-円周率の推定値`pi`と、その自乗`pi2 = pi*pi`について総和を取り、定義通りに期待値と標準偏差を求めるコードが[calc_pi_reduce.cpp](calc_pi_reduce.cpp)である。
+円周率の推定値`pi`と、その自乗`pi2 = pi*pi`について総和を取り、定義通りに期待値と標準偏差を求めるコードが以下の`calc_pi_reduce.cpp`である。
 
 ```cpp
 #include <cstdio>

day4/README.md

@@ -92,9 +92,74 @@ void dump(std::vector<double> &data) {
 }
 ```
 
-あとは適当な条件を与えれば時間発展させることができる。ここでは、「一様加熱」と「温度固定」の二通りを試してみよう。コードはこちら。
+あとは適当な条件を与えれば時間発展させることができる。ここでは、「一様加熱」と「温度固定」の二通りを試してみよう。以下のコードを`thermal.cpp`という名前で保存、実行してみよう。
 
-[thermal.cpp](thermal.cpp)
+```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](fig/thermal.png)
 
@@ -188,9 +253,7 @@ void fixed_temperature(std::vector<double> &lattice) {
 とりあえずメモリに問題なければ「3. 一度まとめてから吐く」が楽なので、今回はこれを採用しよう。メモリが厳しかったり、数万プロセスの計算とかする時にはなにか工夫してくださいまし。
 
 さて、「一度まとめてから吐く」ためには、「各プロセスにバラバラにあるデータを、どこかのプロセスに一括して持ってくる」必要があるのだが、MPIには
-そのものずばり`MPI_Gather`という関数がある。使い方は以下のサンプルを見たほうが早いと思う。
-
-[gather.cpp](gather.cpp)
+そのものずばり`MPI_Gather`という関数がある。使い方はサンプルを見たほうが早い。以下を`gather.cpp`という名前で保存、実行しよう。
 
 ```cpp
 #include <cstdio>
@@ -343,11 +406,104 @@ void fixed_temperature(std::vector<double> &lattice, int rank, int procs) {
 }
 ```
 
-これも一様加熱と同じで、「温度を固定している場所がどのプロセスが担当するどの場所か」を調べる必要があるが、それを考えるのはさほど難しくないだろう。
+これも一様加熱と同じで、「温度を固定している場所がどのプロセスが担当するどの場所か」を調べる必要があるが、それを考えるのはさほど難しくないだろう。そんなわけで完成した並列コードを`thermal_mpi.cpp`という名前で保存しよう。
 
-そんなわけで完成した並列コードがこちら。
+```cpp
+#include <cstdio>
+#include <fstream>
+#include <iostream>
+#include <mpi.h>
+#include <vector>
 
-[thermal_mpi.cpp](thermal_mpi.cpp)
+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();
+}
+```
 
 せっかく並列化したので、高速化したかどうか調べてみよう。一様加熱の計算をさせてみる。