Boost LambdaとCUDA/GPUについての雑感(2)

前の記事(Boost LambdaとCUDA/GPUについての雑感)が結構読まれているので、続編を書くことにした。

>ラムダ式は無名関数を提供する機能。
→ 個人的には、「いかにforループ式を書かないか」というコンセプトの１つ。

> 思えば、CUDA/GPU/Thrustも、いかに「C/C++的なforループを書かないか」という設計になっているとも言える。

であるが、ラムダ式（Boost.Lambda)もCUDA/GPUも、割と同じ感覚（イメージ）で書けるということでもある。

ラムダ式のコード

     1#include <boost/lambda/lambda.hpp>
     2#include <boost/lambda/if.hpp>
     3#include <vector>
     4#include <algorithm>
     5#include <iostream>
     6
     7int main()
     8{
     9  std::vector<int> v{1, 3, 2, 0, 4, 5, 2};
    10  std::for_each(v.begin(), v.end(),
    11    boost::lambda::if_then(boost::lambda::_1 > 2,
    12    std::cout << boost::lambda::_1 << "\n"));
    13}

https://github.com/RuoAndo/qiita/blob/master/lambda_if_then.cpp

これは、リストv{1, 3, 2, 0, 4, 5, 2}に対し、2よりも大きい値を表示する。
実行してみる。。

<pre>
# ./a.out 
3
4
5
</pre>

個人的には、リストの各要素にif_thenという関数を張り付ける（割り当てる）イメージである。

![qiita-ラムダ式.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/163445/31ed23e9-096e-8056-6ac3-4c992a95c432.png)

実は、CUDA/GPUも書いているときの感覚は同じである。

少し長いがコードを見てみる。。

<pre>
     1#include <cuda_runtime.h>
     2#include <stdio.h>
     3
     4void validateResult(float *hostRef, float *gpuRef, const int N)
     5{
     6    double epsilon = 1.0E-8;
     7    bool match = 1;
     8
     9    for (int i = 0; i < N; i++)
    10    {
    11        if (abs(hostRef[i] - gpuRef[i]) > epsilon)
    12        {
    13            match = 0;
    14            printf("Arrays do not match!\n");
    15            printf("host %5.2f gpu %5.2f at current %d\n", hostRef[i], gpuRef[i], i);
    16            break;
    17        }
    18    }
    19
    20    if (match) printf("Arrays match.\n\n");
    21
    22    return;
    23}
    24
    25void genData(float *ip, int size)
    26{
    27    time_t t;
    28    srand((unsigned) time(&t));
    29
    30    for (int i = 0; i < size; i++)
    31    {
    32        ip[i] = (float)( rand() & 0xFF ) / 10.0f;
    33    }
    34
    35    return;
    36}
    37
    38void sumArraysOnHost(float *A, float *B, float *C, const int N)
    39{
    40    for (int idx = 0; idx < N; idx++)
    41    {
    42        C[idx] = A[idx] + B[idx];
    43    }
    44}
    45__global__ void sumArraysOnGPU(float *A, float *B, float *C, const int N)
    46{
    47    int i = blockIdx.x * blockDim.x + threadIdx.x;
    48
    49    if (i < N) C[i] = A[i] + B[i];
    50}
    51
    52int main(int argc, char **argv)
    53{
    54    struct timespec startTime, endTime, sleepTime;
    55
    56    printf("%s Starting...\n", argv[0]);
    57
    58    int dev = 0;
    59    cudaDeviceProp deviceProp;
    60    cudaGetDeviceProperties(&deviceProp, dev);
    61    printf("Using Device %d: %s\n", dev, deviceProp.name);
    62    cudaSetDevice(dev);
    63
    64    int nElem = 1 << 28;
    65    printf("Vector size %d\n", nElem);
    66
    67    size_t nBytes = nElem * sizeof(float);
    68
    69    float *host_A, *host_B, *hostRef, *gpuRef;
    70    host_A     = (float *)malloc(nBytes);
    71    host_B     = (float *)malloc(nBytes);
    72    hostRef = (float *)malloc(nBytes);
    73    gpuRef  = (float *)malloc(nBytes);
    74
    75    genData(host_A, nElem);
    76    genData(host_B, nElem);
    77
    78    memset(hostRef, 0, nBytes);
    79    memset(gpuRef,  0, nBytes);
    80
    81    clock_gettime(CLOCK_REALTIME, &startTime);
    82    sleepTime.tv_sec = 0;
    83    sleepTime.tv_nsec = 123;    84
    85    sumArraysOnCPU(host_A, host_B, hostRef, nElem);
    86
    87    clock_gettime(CLOCK_REALTIME, &endTime);
    88    if (endTime.tv_nsec < startTime.tv_nsec) {
    89        printf("%ld.%09ld", endTime.tv_sec - startTime.tv_sec - 1
    90               ,endTime.tv_nsec + 1000000000 - startTime.tv_nsec);
    91    } else {
    92        printf("%ld.%09ld", endTime.tv_sec - startTime.tv_sec
    93               ,endTime.tv_nsec - startTime.tv_nsec);
    94    }
    95    printf(" sec\n");
    96
    97    float *device_A, *device_B, *device_C;
    98    cudaMalloc((float**)&device_A, nBytes);
    99    cudaMalloc((float**)&device_B, nBytes);
   100    cudaMalloc((float**)&device_C, nBytes);
   101
   102    cudaMemcpy(device_A, host_A, nBytes, cudaMemcpyHostToDevice);
   103    cudaMemcpy(device_B, host_B, nBytes, cudaMemcpyHostToDevice);
   104    cudaMemcpy(device_C, gpuRef, nBytes, cudaMemcpyHostToDevice);   105
   106    int iLen = 1024;
   107    dim3 block (iLen);
   108    dim3 grid  ((nElem + block.x - 1) / block.x);
   109
   110    clock_gettime(CLOCK_REALTIME, &startTime);
   111    sleepTime.tv_sec = 0;
   112    sleepTime.tv_nsec = 123;
   113
   114    sumArraysOnGPU<<<grid, block>>>(device_A, device_B, device_C, nElem);
   115    cudaDeviceSynchronize();
   116    cudaGetLastError() ;
   117
   118    clock_gettime(CLOCK_REALTIME, &endTime);
   119    if (endTime.tv_nsec < startTime.tv_nsec) {
   120        printf("%ld.%09ld", endTime.tv_sec - startTime.tv_sec - 1
   121               ,endTime.tv_nsec + 1000000000 - startTime.tv_nsec);
   122    } else {
   123        printf("%ld.%09ld", endTime.tv_sec - startTime.tv_sec
   124               ,endTime.tv_nsec - startTime.tv_nsec);
   125    }
   126    printf(" sec\n");
   127
   128    cudaMemcpy(gpuRef, device_C, nBytes, cudaMemcpyDeviceToHost);
   129
   130    validateResult(hostRef, gpuRef, nElem);
   131
   132    cudaFree(device_A);
   133    cudaFree(device_B);
   134    cudaFree(device_C);
   135
   136    free(host_A);
   137    free(host_B);
   138    free(hostRef);
   139    free(gpuRef);
   140
   141    return(0);
   142}
</pre>

https://github.com/RuoAndo/qiita/blob/master/simple.cu

ここで、この２つの関数であるが、なぜGPUの方の関数はforループがないかというと、上のラムダ式と似たような仕組みになっているからである。

<pre>
    38void sumArraysOnCPU(float *A, float *B, float *C, const int N)
    39{
    40    for (int idx = 0; idx < N; idx++)
    41    {
    42        C[idx] = A[idx] + B[idx];
    43    }
    44}
    45__global__ void sumArraysOnGPU(float *A, float *B, float *C, const int N)
    46{
    47    int i = blockIdx.x * blockDim.x + threadIdx.x;
    48
    49    if (i < N) C[i] = A[i] + B[i];
    50}
</pre>

図にするとこんな感じであらう。

![qiita-sumArrayOnGPU.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/163445/32548012-e0ea-f366-dabe-2a6fa13ff5cf.png)

つまり、GPUデバイス内の各スレッドに対して、sumArrayOnGPUを張り付ける（割り当てる）イメージである。

実行してみる。。

<pre>
# ./a.out 
./a.out Starting...
Using Device 0: Quadro GV100
Vector size 268435456
1.719520959 sec
0.005400906 sec
</pre>

（｀ー´）b

実は、Google発のMapRuduceも同じ感覚で書ける。。

https://research.google/pubs/pub62/

アレは、各データの要素にMap関数を張り付けていくイメージである。

このあたりのことは、LISPをやりまくると開眼するらしい。。 
 (ﾟДﾟ;)