はじめに

こないだ買ってきたNVIDIAのjetson nanoでπの値を計算し、CPU演算とGPU演算とでどれくらい速度が違うのか比較してみました。
なお、jetson nanoにはopenMPI(バージョン4.0)をインストールしてあります。これは将来複数台のjetson nanoをクラスタ(ミニスパコン)化する計画があるためです。

CPUでの計算

jetson nanoにはARM Cortex-A57相当のコアが4コア載っています。まずはこのCPUコアのみでの計算速度を試してみます。

使用した公式は収束が遅いことで有名なライプニッツの公式です。

leibniz_c.c

/*
	leibniz_c.c
	calculation of PI by Leibniz's formula
	(MPI version)
*/
#include <mpi.h>
#include <stdio.h>
#include <math.h>

#define N 1000000000

int main(int argc, char *argv[])
{
	long long i,j;
	int irank,numprocs;
	char hostname[MPI_MAX_PROCESSOR_NAME];
	long long n=N;
	double x,y,sum,pi;
	double tstart,tend;

	MPI_Init(&argc,&argv);
	MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
	MPI_Comm_rank(MPI_COMM_WORLD,&irank);
	if(irank==0)
	{
		printf("PROCESSORS      : %d\n",numprocs);
		printf("N               : %lld\n",n);
	}
	MPI_Bcast(&n,1,MPI_INT,0,MPI_COMM_WORLD);

	tstart = MPI_Wtime();
	pi = 0.0;
	sum = 0.0;
	for(i=irank+1;i<=n;i+=numprocs)
	{
		j = i - 1;
		if((i%2)==0)
		{
			y = -1.0;
		}
		else
		{
			y = 1.0;
		}
		//x = (double)pow(-1,j)*4.0/(double)(2*j+1);
		x = y*4.0/(double)(2*j+1);
		sum = sum + x;
	}
	MPI_Reduce(&sum,&pi,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
	tend = MPI_Wtime();

	if(irank==0)
	{
		printf("PI(Calculated)  : %.30f\n",pi);
		printf("PI(M_PI)        : %.30f\n",M_PI);
		printf("Relative Error  : %.30f\n",fabs(pi-M_PI));
		printf("TIME            : %9.3f\n",tend-tstart);
	}

	MPI_Finalize();
}

MPIを使用して4コア並列実行ができるようにしました。
プログラムのコンパイルは次のコマンドで行います。

$ /usr/local/bin/mpicc -o leibniz_c -Os leibniz_c.c -lm

コンパイルできたら実行してみます。mpiexecコマンドの-npオプションで実行時の並列数を指定します。

1コア(逐次計算)のとき

shingo@jetson-1:~/cuda$ /usr/local/bin/mpiexec -np 1 leinbniz_c
PROCESSORS      : 1
N               : 1000000000
PI(Calculated)  : 3.141592652588050427198140823748
PI(M_PI)        : 3.141592653589793115997963468544
Relative Error  : 0.000000001001742688799822644796
TIME            :    21.063

4コア並列実行のとき

shingo@jetson-1:~/cuda$ mpiexec -np 4 leibniz_c
PROCESSORS      : 4
N               : 1000000000
PI(Calculated)  : 3.141592652591443268761395302135
PI(M_PI)        : 3.141592653589793115997963468544
Relative Error  : 0.000000000998349847236568166409
TIME            :     5.282

比較

1コアの逐次計算と比べて、ほぼコア数なりに速くなっていることがわかります。

GPUでの計算(CUDA+MPI)

次に上記コードをCUDA用に書き換えてみました。GPUで級数展開の部分を行って配列にストアし、計算が終わったらMPI_Reduceで部分和の集計をとってπの値を求めます。

leibniz_cuda.cu

/*
  	leibniz_cuda.cu
  	calculation of PI by Leibniz's formula
	(CUDA/MPI version)
*/
#include <mpi.h>
#include <stdio.h>
#include <math.h>
#include <cuda_runtime.h>

#define NP 512
#define NB 8
#define N  1000000000

#define CHECK(call) \
{ \
	const cudaError_t error = call; \
	if(error != cudaSuccess) \
	{ \
		printf("Error : %s:%d, ",__FILE__,__LINE__); \
		printf("code : %d, reason : %s\n",error, cudaGetErrorString(error)); \
		exit(1); \
	} \
}

/*
	CUDA kernel routine
	calculates partial sum of Leibniz's formula
*/
__global__ void leibniz(double *X,long long n)
{
	int k;
	long long i,j;
	double x,y;
	for(k=0;k<NB;k++)
	{
		i=n+threadIdx.x+1+k*NP;
		j = i - 1;
		if((i%2)==0)
		{
			y=-1.0;
		}
		else
		{
			y=1.0;
		}
		x = y*4.0/(double)(2*j+1);
		X[threadIdx.x] = X[threadIdx.x] + x;
	}
}

int main(int argc, char *argv[])
{
	long long i,j;
	int irank,numprocs;
	char hostname[MPI_MAX_PROCESSOR_NAME];
	long long n=N;
	double pi,sum0;
	double tstart,tend;
	/* setup device */
	int dev = 0;
	cudaDeviceProp deviceProp;
	CHECK(cudaGetDeviceProperties(&deviceProp,dev));
	CHECK(cudaSetDevice(dev));

	MPI_Init(&argc,&argv);
	MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
	MPI_Comm_rank(MPI_COMM_WORLD,&irank);
	if(irank==0)
	{
		printf("USING DEVICE %d : %s\n",dev,deviceProp.name);
		printf("PROCESSORS(CPU) : %d\n",numprocs);
		printf("GPU THREADS     : %d\n",NP);
		printf("N               : %lld\n",n);
	}
	MPI_Bcast(&n,1,MPI_INT,0,MPI_COMM_WORLD);

	tstart = MPI_Wtime();

	/* setup array */
	size_t nBytes = NP * sizeof(double);
	double *h_X,*d_X,*sum;
	h_X = (double *)malloc(nBytes);
	sum = (double *)malloc(nBytes);
	CHECK(cudaMallocHost((double**)&d_X,nBytes));
	memset(h_X,0,nBytes);
	CHECK(cudaMemcpy(d_X,h_X,nBytes,cudaMemcpyHostToDevice));

	/* start kernel */
	cudaStream_t stream;
	cudaStreamCreateWithFlags(&stream,cudaStreamNonBlocking);

	for(i=irank*NP*NB;i<n;i+=NP*numprocs*NB)
	{
		leibniz<<<1,NP,nBytes,stream>>>(d_X,i);
	}

	CHECK(cudaStreamSynchronize(stream));
	CHECK(cudaMemcpy(h_X,d_X,nBytes,cudaMemcpyDeviceToHost));
	CHECK(cudaStreamDestroy(stream));
	CHECK(cudaFreeHost(d_X));
	CHECK(cudaDeviceReset());

	MPI_Allreduce(h_X,sum,NP,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
	sum0 = 0.0;
	for(i=irank;i<NP;i+=numprocs)
	{
		sum0 = sum0 + sum[i];
	}
	pi = 0.0;
	MPI_Reduce(&sum0,&pi,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
	tend = MPI_Wtime();

	if(irank==0)
	{
		printf("PI(Calculated)  : %.30f\n",pi);
		printf("PI(M_PI)        : %.30f\n",M_PI);
		printf("Relative Error  : %.30f\n",fabs(pi-M_PI));
		printf("TIME            : %9.3f\n",tend-tstart);
	}

	MPI_Finalize();
	free(sum);
	free(h_X);
}

MPIを使用してCPU4コアでの並列実行(1台のGPUを4つのプロセスで分け分けして使用)ができるようにしました。次のコマンドでコンパイルします。

$ nvcc -o leibniz_cuda_1 -arch=sm_53 -Xcompiler -Os -I/usr//local/include -L/usr/local/lib leibniz_cuda_1.cu -lm -lmpi

1プロセス(GPU512スレッド並列)実行のとき

shingo@jetson-1:~/cuda$ /usr/local/bin/mpiexec -np 1 leibniz_cuda_1
USING DEVICE 0 : NVIDIA Tegra X1
PROCESSORS(CPU) : 1
GPU THREADS     : 512
N               : 1000000000
PI(Calculated)  : 3.141592652589856538014601028408
PI(M_PI)        : 3.141592653589793115997963468544
Relative Error  : 0.000000000999936577983362440136
TIME            :     3.965

4プロセス(GPU512スレッドx4プロセス並列)実行のとき

shingo@jetson-1:~/cuda$ /usr/local/bin/mpiexec -np 4 leibniz_cuda_1
USING DEVICE 0 : NVIDIA Tegra X1
PROCESSORS(CPU) : 4
GPU THREADS     : 512
N               : 1000000000
PI(Calculated)  : 3.141592652590006196078320499510
PI(M_PI)        : 3.141592653589793115997963468544
Relative Error  : 0.000000000999786919919642969035
TIME            :     2.987

比較

4並列のときに4倍速くなっていませんが、これはGPUのコア数が足りないために実行待ちが発生しているものと思われます。

まとめ

CUDAへの移植としてはごくごく素朴なコードですが、それでも倍くらいは速くなっています。