Abstract

FriendlyElec's NanoPi-NEO4 SBC uses RK3399 SoC.

RK3399 has ARM Mali-T864 GPU, and this GPU can be used for GPGPU via OpenCL.

This article explains how to program ARM Mali-T864 GPU via OpenCL, and compares performance of GPGPU program with MPI-based CPU program.

はじめに

FriendlyElecのシングルボードコンピュータ、NanoPi-NEO4のSoCであるところのRK3399にはARM Mali-T864というGPUが搭載されています。このGPUはOpenCLに対応しています。

そこでOpenCLでのGPGPU計算に挑戦してみましたので報告いたします(実施例が少ないですし)。ARM Mali-T600シリーズ以降のGPUが載っているシングルボードコンピュータでも実施可能だと思います。

なお、NanoPi-NEO4のOSはFriendlyCoreを使用しました。ARM Mali用のOpenGL ESライブラリやOpenCLライブラリがあらかじめインストールされています。

OpenCLのプログラム

最初に今回OpenCLのプログラムの例題として作成した、Leibnizの級数によって円周率を求めるプログラムを掲載します。

leibniz_opencl.c


/*
    leibniz_opencl.c
    calculation of PI by Leibniz's formula
    (OpenCL version)
*/
#include <CL/cl.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>

#define MEM_SIZE (4096)
#define MAX_SOURCE_SIZE (0x100000)
#define N 2000000

float local_buf[MEM_SIZE];

double getSeconds(clock_t from, clock_t till)
{
    return (double)(till-from)/CLOCKS_PER_SEC;
}

int main(int argc, char *argv[])
{
    cl_device_id device_id = NULL;
    cl_context context = NULL;
    cl_command_queue command_queue = NULL;
    cl_mem memobj = NULL;
    cl_program program = NULL;
    cl_kernel kernel = NULL;
    cl_platform_id platform_id = NULL;
    cl_uint ret_num_devices;
    cl_uint ret_num_platforms;
    cl_int ret;

    clock_t start_calc_t,end_gpu_t,end_calc_t;

    long i,j;
    long n=N;
    float x,sum,pi;
    float tstart,tgpu,tend;

    printf("N              : %ld\n",n);

    for(i=0;i<MEM_SIZE;i++)
    {
        local_buf[i]=0.0f;
    }

    /* OpenCL setup */
    char *source_str=
    "__kernel void leibniz(__global float* X,const int n) \
    { \
        long i,j,k; \
        float x,y; \
        int id= get_global_id(0); \
        int isize = get_global_size(0); \
        for(i=n-id;i>0;i-=isize) \
        { \
            j = i - 1; \
            if((j%2)==0) y = 1.0; \
            else y = -1.0; \
            x = y*4.0/(float)(2*j+1); \
            X[id] = X[id] + x; \
        } \
    }";
    size_t source_size = strlen(source_str);

    ret = clGetPlatformIDs(1,&platform_id, &ret_num_platforms);
    if(ret != CL_SUCCESS)
    {
        printf("Error : clGetPlatformIDs()\n");
        return -1;
    }
    ret = clGetDeviceIDs(platform_id,CL_DEVICE_TYPE_DEFAULT,1,
        &device_id,&ret_num_devices);
    if(ret != CL_SUCCESS)
    {
        printf("Error : clGetDeviceIDs()\n");
        return -1;
    }
    context = clCreateContext(NULL,1,&device_id,NULL,NULL,&ret);
    command_queue = clCreateCommandQueue(context,device_id,0,&ret);
    memobj = clCreateBuffer(context,CL_MEM_READ_WRITE,
        MEM_SIZE*sizeof(float),NULL,&ret);
    program = clCreateProgramWithSource(context,1,
        (const char **)&source_str,(const size_t *)&source_size,&ret);
    ret = clBuildProgram(program, 1, &device_id, NULL, NULL, NULL);
    if(ret != CL_SUCCESS)
    {
        size_t len;
        char buffer[2048];
        printf("Error : clBuildProgram()\n");
        clGetProgramBuildInfo(program,device_id,CL_PROGRAM_BUILD_LOG,
            sizeof(buffer),buffer, &len);
        printf("\t%s\n",buffer);
        exit(1);
    }
    kernel = clCreateKernel(program, "leibniz",&ret);
    if(ret != CL_SUCCESS)
    {
        printf("Error : clCreateKernel()\n");
        exit(1);
    }
    ret = clSetKernelArg(kernel,0,sizeof(cl_mem),(void *)&memobj);
    ret |= clSetKernelArg(kernel,1,sizeof(int),(void *)&n);
    if(ret != CL_SUCCESS)
    {
        printf("Error : clSetKernelArg()\n");
        exit(1);
    }
    /* execute kernel */
    size_t globalWorkSize[] = {MEM_SIZE};
    size_t localWorkSize[] = {1};
    start_calc_t = clock();
    ret = clEnqueueNDRangeKernel(command_queue,kernel,1,NULL,
        globalWorkSize,localWorkSize,0,NULL,NULL);
    if(ret != CL_SUCCESS)
    {
        printf("Error : executing kernel\n");
        exit(1);
    }
    end_gpu_t = clock();
    /* reading result */
    ret = clEnqueueReadBuffer(command_queue,memobj,CL_TRUE,0,
        MEM_SIZE*sizeof(float),local_buf,0,NULL,NULL);
    pi = 0.0;
#pragma omp parallel for reduction(+ : pi)
    for(i=0;i<MEM_SIZE;i++)
    {
        pi = pi + local_buf[i];
    }
    end_calc_t = clock();   

    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(GPU TIME) : %10.7f\n",(float)getSeconds(start_calc_t,end_gpu_t));
    printf("TIME(CALC TIME): %10.7f\n",(float)getSeconds(start_calc_t,end_calc_t));

    /* Finalization */
    ret = clFlush(command_queue);
    ret = clFinish(command_queue);
    ret = clReleaseKernel(kernel);
    ret = clReleaseProgram(program);
    ret = clReleaseMemObject(memobj);
    ret = clReleaseCommandQueue(command_queue);
    ret = clReleaseContext(context);
}

GPU側プログラム(カーネル)構成の準備

まず、計算環境(コンテクスト)の設定を行います。

    ret = clGetPlatformIDs(1,&platform_id, &ret_num_platforms);
    if(ret != CL_SUCCESS)
    {
        printf("Error : clGetPlatformIDs()\n");
        return -1;
    }
    ret = clGetDeviceIDs(platform_id,CL_DEVICE_TYPE_DEFAULT,1,
        &device_id,&ret_num_devices);
    if(ret != CL_SUCCESS)
    {
        printf("Error : clGetDeviceIDs()\n");
        return -1;
    }
    context = clCreateContext(NULL,1,&device_id,NULL,NULL,&ret);
    command_queue = clCreateCommandQueue(context,device_id,0,&ret);
    memobj = clCreateBuffer(context,CL_MEM_READ_WRITE,
        MEM_SIZE*sizeof(float),NULL,&ret);

次に、GPU側のプログラムをコンパイルします。OpenGLES2のシェーダのプログラムと同様、OpenCLのGPU側のプログラムはソースをcharの配列として保持しており、実行時にコンパイルする形になっています。別ファイルにして実行時に読み出してコンパイルする方法もありますが、今回はcharの配列として埋め込んでいます。
GPU側のプログラムソースは次の通りです。ライプニッツの級数計算で円周率の値を求めるものです :

    char *source_str=
    "__kernel void leibniz(__global float* X,const int n) \
    { \
        long i,j,k; \
        float x,y; \
        int id= get_global_id(0); \
        int isize = get_global_size(0); \
        for(i=n-id;i>0;i-=isize) \
        { \
            j = i - 1; \
            if((j%2)==0) y = 1.0; \
            else y = -1.0; \
            x = y*4.0/(float)(2*j+1); \
            X[id] = X[id] + x; \
        } \
    }";
    size_t source_size = strlen(source_str);

GPU側プログラムはC言語とよく似ています(というよりそっくりです)ので、何をしているかはソースを見ればわかると思います。

本当は倍精度実数(double)で計算したいところですが、ARMの"Mali-T600 Series GPU OpenCL Developer Guide"によればdoubleは予約(使用できない)となっているので単精度実数(float)での計算としました。

get_global_id(0)はこのプログラムが実行されるときのスレッド番号(MPI並列におけるランク番号に相当)、get_global_size(0)はグループのサイズ(今回は1次元なのでスレッドの数に等しい。MPI並列におけるプロセス数に相当)をあらわします。

このプログラムをOpenCLの実行時にコンパイルします。

    program = clCreateProgramWithSource(context,1,
        (const char **)&source_str,(const size_t *)&source_size,&ret);
    ret = clBuildProgram(program, 1, &device_id, NULL, NULL, NULL);
    if(ret != CL_SUCCESS)
    {
        size_t len;
        char buffer[2048];
        printf("Error : clBuildProgram()\n");
        clGetProgramBuildInfo(program,device_id,CL_PROGRAM_BUILD_LOG,
            sizeof(buffer),buffer, &len);
        printf("\t%s\n",buffer);
        exit(1);
    }
    kernel = clCreateKernel(program, "leibniz",&ret);
    if(ret != CL_SUCCESS)
    {
        printf("Error : clCreateKernel()\n");
        exit(1);
    }

clCreateProgram()を先のGPU側プログラムソースの先頭の番地とサイズを指定して呼び出し、次にclBuileProgram()を呼び出してコンパイルします。

GPU側プログラムソースに構文上のエラーがある場合はclGetProgramBuildInfo()でエラーの内容を調べることもできます。
エラーなくコンパイルできた場合はclBuildProgramが値CL_SUCCESSを返しますので、clCreateKernel()でGPU側プログラムのカーネルを生成します。

以上で実行準備が整いました。

カーネルに引数を渡して実行する

生成したカーネルにclSetKernelArg()で実行時の引数を渡し、ワークサイズ(今回は並列スレッド数に等しい)を指定してclEnqueueNDRangeKernel()を呼び出すとGPU側でプログラムが走ります。

    ret = clSetKernelArg(kernel,0,sizeof(cl_mem),(void *)&memobj);
    ret |= clSetKernelArg(kernel,1,sizeof(int),(void *)&n);
    if(ret != CL_SUCCESS)
    {
        printf("Error : clSetKernelArg()\n");
        exit(1);
    }
    /* execute kernel */
    size_t globalWorkSize[] = {MEM_SIZE};
    size_t localWorkSize[] = {1};
    start_calc_t = clock();
    ret = clEnqueueNDRangeKernel(command_queue,kernel,1,NULL,
        globalWorkSize,localWorkSize,0,NULL,NULL);
    if(ret != CL_SUCCESS)
    {
        printf("Error : executing kernel\n");
        exit(1);
    }
    end_gpu_t = clock();
    /* reading result */
    ret = clEnqueueReadBuffer(command_queue,memobj,CL_TRUE,0,
        MEM_SIZE*sizeof(float),local_buf,0,NULL,NULL);

計算結果はclEnqueueReadBuffer()で読み出すことができます。今回は通常の方法(GPU側のバッファをCPU側にコピー)で読み出しましたが、ARM Maliの場合はコピーしなくても読み出すことができ、コピーに伴う性能低下を最小限にすることができるとARMの資料にあります。

ただし今回はコピーはGPUが計算した円周率の部分和の配列をCPU側にコピーする1回だけですので、通常のコピーで行いました。

なお、clEnqueueReadBuffer()の第3引数を上のようにCL_TRUEとすると読み出し完了まで待つブロッキング読み出しとなります。CL_FALSEとするとノンブロッキング読み出しとなり、バッファの内容はcl_Finish()を呼ぶまでは保証されません。

後始末

GPUでの計算が終わったら後始末します。

    /* Finalization */
    ret = clFlush(command_queue);
    ret = clFinish(command_queue);
    ret = clReleaseKernel(kernel);
    ret = clReleaseProgram(program);
    ret = clReleaseMemObject(memobj);
    ret = clReleaseCommandQueue(command_queue);
    ret = clReleaseContext(context);

性能比較

上記のプログラムをコンパイルし、CPUのみでの計算(NanoPi-NEO4は6コアですので、MPIで6並列実行とします)と所要時間を比較してみます。

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 2000000

int main(int argc, char *argv[])
{
    long long i,j;
    int irank,numprocs;
    char hostname[MPI_MAX_PROCESSOR_NAME];
    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              : %ld\n",n);
    }
    MPI_Bcast(&n,1,MPI_INT,0,MPI_COMM_WORLD);

    tstart = MPI_Wtime();
    pi = 0.0;
    sum = 0.0;
    for(i=n-irank;i>0;i-=numprocs)
    {
        j = i - 1;
        if((j%2)==0) y = 1.0;
        else y = -1.0;
        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           : %10.7f\n",tend-tstart);
    }

    MPI_Finalize();
}

コンパイルは次のコマンドで行います:

$ mpicc -o leibniz_c -Os leibniz_c.c -lm

OpenCLのプログラムは次のコマンドでコンパイルします:


$ cc -o leibniz_opencl -Os -fopenmp leibniz_opencl.c -lm -lOpenCL

CPUでの計算

CPUでの計算結果は次の通りです。

pi@NanoPi-NEO4:~$ mpiexec -np 1 ./leibniz_c
PROCESSORS     : 1
N              : 2000000
PI(Calculated) : 3.141592153589793490198189829243
PI(M_PI)       : 3.141592653589793115997963468544
Relative Error : 0.000000499999999625799773639301
TIME           :  0.0478882
pi@NanoPi-NEO4:~$ mpiexec -np 6 ./leibniz_c
PROCESSORS     : 6
N              : 2000000
PI(Calculated) : 3.141592153589775726629795826739
PI(M_PI)       : 3.141592653589793115997963468544
Relative Error : 0.000000500000017389368167641805
TIME           :  0.0085954

1コアのみのとき(逐次)の実行時間は約48ms程度かかります。200万回ループをまわしていることを考えるとNanoPi-NEO4の計算性能はなかなかのものです。

6コア総動員で実行したときは約8.5ms程度となり、大体コア数なりに速くなっていることがわかります。

GPUでの計算

GPUでの計算結果は次の通りです。現時点ではMaliのOpenCLが倍精度実数に対応していないので単精度実数での計算です。なお、部分和を集計して円周率の値を求める部分はCPU計算です。OpenMPで並列数を変えて実行時間を比較してみました。
GPUのスレッド並列数は4096としています。

pi@NanoPi-NEO4:~$ export OMP_NUM_THREADS=1;./leibniz_opencl
N              : 2000000
PI(Calculated) : 3.141593456268310546875000000000
PI(M_PI)       : 3.141592653589793115997963468544
Relative Error : 0.000000802678517430877036531456
TIME(GPU TIME) :  0.0000720
TIME(CALC TIME):  0.0020110
pi@NanoPi-NEO4:~$ export OMP_NUM_THREADS=6;./leibniz_opencl
N              : 2000000
PI(Calculated) : 3.141591310501098632812500000000
PI(M_PI)       : 3.141592653589793115997963468544
Relative Error : 0.000001343088694483185463468544
TIME(GPU TIME) :  0.0000750
TIME(CALC TIME):  0.0263910

所要時間の集計はGPUでの計算部分(部分和の配列を求めるところまで)と全体とを求めてあります。

全体での計算時間は集計部分を逐次としたとき(OMP_NUM_THREADS=1)のほうが圧倒的に短く、わずか2msで計算が終わっています。GPU計算部分の所要時間はわずか75usです。

集計部分を6スレッド並列としたとき(OMP_NUM_THREADS=6)は集計の所要時間が10倍以上かかっています。この程度の配列サイズであればスレッド並列化するより逐次計算のほうがスレッド並列化に伴うオーバーヘッド分得だということがわかります。

CPU計算と比べて約4倍の性能が出ていることがわかります。集計の計算部分もGPU化すればさらに高速化できる可能性があります。

まとめ

ARM Mali-T864でGPGPUしてみました。NanoPi-NEO4のCPU6コアを総動員してぶん回したときのさらに4倍程度の速度向上が見込めます。
C言語でMaliのOpenCLプログラムを書く例は少ないので、拙いプログラム例ですが参考になれば幸いです。

参照文献

Mali-T600 Series GPU OpenCL(version 1.0) Developer Guide(DUI0538E), ARM, 2012-2013
牛島省,OpenMPによる並列プログラミングと数値計算法, 丸善出版, 平成27年

GPGPU by NanoPi-NEO4 via OpenCL