オールドタイプ並列化

理系ネタ J

404 Blog Not Found:C - でOpenMP使ったら、+二行で範囲10億が2秒
ふむー。最近は便利なものが！
比較として馬鹿正直に分散メモリMPIで並列コード書いてみた。10行くらい追加。性能は同じくらい。BLOCKSIZEは1メガで。
1ブロック計算し終わったらMaster CPUに結果を転送開始し同時に次のブロックを計算とかして通信の隠蔽をしたらもっと速いかもしれない。

# gcc + MPICH (MPICHについては記事の最後に)
gcc -O6 pr-mpi.c -I/usr/local/include -L/usr/local/lib -lmpich -lpthread

time /usr/local/bin/mpirun -np 4 a.out
0.060u 0.017s 0:01.09 6.4%      0+0k 0+0io 0pf+0w

time /usr/local/bin/mpirun -np 4 a.out 0 10000000000
0.065u 0.016s 0:14.98 0.4%      0+0k 0+0io 0pf+0w
# Intel CC + INTEL MPI (計算は速いけどmpirunのオーバーヘッドが数秒ある)
icc -O -I/opt/intel/impi/3.2.1.009/include64 pr-mpi.c -L/opt/intel/impi/3.2.1.009/lib64 -lmpi

time /opt/intel/impi/3.2.1.009/bin64/mpirun -np 4 a.out
0.374u 0.129s 0:06.35 7.7%      0+0k 0+0io 0pf+0w

time /opt/intel/impi/3.2.1.009/bin64/mpirun -np 4 a.out 0 10000000000
0.364u 0.134s 0:12.59 3.8%      0+0k 0+0io 0pf+0w

以下並列版のdiff

--- prime.c	2010-08-06 16:18:30.000000000 +0900
+++ pr-mpi.c	2010-08-06 16:29:15.000000000 +0900
@@ -4,6 +4,7 @@
 #include <limits.h>
 #include <errno.h>
 #include <math.h>
+#include <mpi.h>
 
 #include <stdint.h>
 typedef uint32_t U32;
@@ -40,6 +41,12 @@
 }
 
 int main(int argc, char **argv) {
+    MPI_Init(&argc, &argv); //オマジナイ
+    int ncpu, cpu_id;
+    MPI_Comm_size(MPI_COMM_WORLD, &ncpu);//並列CPU数
+    MPI_Comm_rank(MPI_COMM_WORLD, &cpu_id); //自分の番号 = 0, 1, .. ,ncpu-1
+    int id_io = 0; // I/O 担当のCPU
+
     bitmap *sieve = bitmap_new(0, SIEVEFILE);
     if (!sieve){
         perror(SIEVEFILE);
@@ -71,11 +78,15 @@
 
     for(i = 0; i*block_size < size; i++)
     {
+	if ( (i % ncpu) != cpu_id) continue;
 	U64 size1 = block_size;
 	if (i*block_size + size1 > size) size1 = size - i*block_size;
         sieve_block(b, sieve, start, i*block_size, size1, pmax);
     }
 
+    bitmap *b_global = bitmap_new(size >> 1, NULL);
+    MPI_Reduce ( b->map, b_global->map, (size>>1)>>3, MPI_CHAR, MPI_BAND, id_io, MPI_COMM_WORLD); //take bitwise-AND of b->map among all CPU's, store result to b_global->map
+
 #ifdef VERBOSE
     // Check number of primes. See http://mathworld.wolfram.com/PrimeCountingFunction.html
     U64 np=0; 
@@ -96,6 +107,7 @@
 #ifdef VERBOSE
     printf("saving %s\n", filename);
 #endif
-    bitmap_save(b, filename);
+    if (cpu_id == id_io) bitmap_save(b_global, filename); // only one cpu does I/O
+     MPI_Finalize(); //オマジナイ
     return 0;
 }

MPICH