CUDA 行列の内積での疑問点

Question

疑問点のあるコード

c言語歴5ヶ月,cuda歴5日の初心者です。現在c++で多層パーセプトロンのライブラリを
作成中なのですがcpuで計算させるためのコードは完成し、GPUバージョンを作ろうと
cudaに手を出したのですが内積の処理を行うサンプルコードがわけわからなくて
大変苦しんでおります。
ある程度qiitaなどを参考に理解できましたがどうしてもわからない点が二点あるので
大変初歩的な質問では御座いますがどうかご教授よろしくお願いします。

/////////////////////////////////////////////////////////////////////////
//
//      このコードでは入力配列Adと入力配列Bdの内積を行い結果を配列Cdへ格納する
//
//                       Ad  dot  Bd  =  Cd
//                        ↑        ↑      ↑
//                   (200,300)(300,500)(200,500)
//
//////////////////////////////////////////////////////////////////////////

#include <stdio.h>
#include <stdlib.h>
#include <cuda.h>
#include <cuda_runtime.h>
#include <sys/time.h>

//////////////////////////////////////////////////////////
//     ↓GPU側で内積を行う関数
//////////////////////////////////////////////////////////

__global__ 
void matmul(float *A, float *B, float *C, int l, int m, int n)
{
    int i, j, k;
    float sum;

    i = blockIdx.y * blockDim.y + threadIdx.y;
    j = blockIdx.x * blockDim.x + threadIdx.x;

    sum = 0.0;
    for (k = 0; k < m; k++) {
        sum += A[i * m + k] * B[k * n + j];
    }
    C[i*n+j] = sum;
}

//////////////////////////////////////////////////////////
//     ↓行列に乱数代入
//////////////////////////////////////////////////////////

void init_matrix(float *m, int h, int w)
{
    int i, j;
    for (i = 0; i < h; i++)
        for (j = 0; j < w; j++)
            m[i * w + j] = (float)random();
}

//////////////////////////////////////////////////////////
//     ↓処理時間計算
//////////////////////////////////////////////////////////
double get_elapsed_time(struct timeval *begin, struct timeval *end)
{
    return (end->tv_sec - begin->tv_sec) * 1000
            + (end->tv_usec - begin->tv_usec) / 1000.0;
}

//////////////////////////////////////////////////////////
//     ↓メイン関数
//////////////////////////////////////////////////////////
int main(int argc, char *argv[])
{
    float *Ad, *Bd, *Cd;//GPU側のメモリ(配列)
    float *Ah, *Bh, *Ch;//CPU側のメモリ(配列)
    struct timeval t1, t2;//処理時間測定用変数

    //配列Ad,Ahの形状定義
    int a_col = 200;//行
    int a_row = 300;//列

    //配列Bd,Bhの形状定義
    int b_col = 300;//行
    int b_row = 500;//列

    //配列Cd,Chの形状定義
    int c_col = 200;//行
    int c_row = 500;//列

    //配列Ah,Adを動的確保し乱数を流し込む(CPU側GPU側共に)
    Ah = (float *)malloc(sizeof(float) * a_col * a_row);
    cudaMalloc((void **)&Ad, sizeof(float) * a_col * a_row);
    init_matrix(Ah, a_col, a_row);

    //乱数を流しこんだAhの内容をデバイス側のAdに代入?
    cudaMemcpy(Ad, Ah, sizeof(float) * a_col * a_row,cudaMemcpyHostToDevice);

    //配列Bh,Bdを動的確保し乱数を流し込む(CPU側GPU側共に)
    Bh = (float *)malloc(sizeof(float) * b_col * b_row);
    cudaMalloc((void **)&Bd, sizeof(float) * b_col * b_row);
    init_matrix(Bh, b_col, b_row);

    //乱数を流しこんだBhの内容をデバイス側のBdに代入?
    cudaMemcpy(Bd, Bh, sizeof(float) * b_col * b_row,cudaMemcpyHostToDevice);

    //配列Ch,Cdを動的確保(CPU側GPU側共に)
    Ch = (float *)malloc(sizeof(float) * c_col * c_row);
    cudaMalloc((void **)&Cd, sizeof(float) * c_col * c_row);

    //スレッドの同期?
    cudaDeviceSynchronize();

    //内積を開始した時間代入
    gettimeofday(&t1, NULL);

    // カーネル関数呼び出し
    matmul<<<dim3(a_row / 16, b_row / 16),dim3(16, 16)>>>(Ad, Bd, Cd, b_row, a_col, a_row);

    //スレッドの同期?
    cudaDeviceSynchronize();

    //内積を終了した時間代入
    gettimeofday(&t2, NULL);

    //処理にかかった時間表示
    printf("Elapsed time: %f msec\n", get_elapsed_time(&t1, &t2));
    
    // 計算結果をGPUからCPUのメモリに代入
    cudaMemcpy(Ch, Cd, sizeof(float) * b_row * a_row, cudaMemcpyDeviceToHost);

    //GPU側に確保したメモリ解放
    cudaFree(Ad);
    cudaFree(Bd);
    cudaFree(Cd);

    return 0;
}

疑問点その1:GPU側で内積を行う関数について。

・普通内積を計算するならこうなるはずです。

for(int i=0;i<l;i++){       ←ない
    for(int j=0;j<n;j++){   ←ない
        sum = 0.0;
        for (int k = 0; k < m; k++) {
            sum += A[i * m + k] * B[k * n + j];
        }
    C[i*n+j] = sum;
    }
}

しかし上２つのfor文が無く代わりに
i = blockIdx.y * blockDim.y + threadIdx.y;
j = blockIdx.x * blockDim.x + threadIdx.x;
が記述されています。
なぜこのような書き方ができるのでしょうか？

疑問点その2:カーネル関数呼び出し時のグリッド数とブロック数

今回のコードでは下記のように記述されています。

matmul<<<dim3(a_row / 16, b_row / 16),dim3(16, 16)>>>(Ad, Bd, Cd, b_row, a_col, a_row);

カーネル関数を呼び出すときは下記のように呼びだせば良いと知っています。
しかし適切なグリッド数とブロック数はどのように求めればよいのでしょうか？

func<<<grid(dim3型), block(dim3型)>>>();

Comments

[takabosoft]

ちょっと世代が古いですが、CUDAの書籍で日本語のものがありますので、そちらで勉強された方がトータルではいいかもしれません。
https://www.amazon.co.jp/CUDA-%E3%83%97%E3%83%AD%E3%83%95%E3%82%A7%E3%83%83%E3%82%B7%E3%83%A7%E3%83%8A%E3%83%AB-%E3%83%97%E3%83%AD%E3%82%B0%E3%83%A9%E3%83%9F%E3%83%B3%E3%82%B0-impress-gear/dp/4844338919/ref=pd_lpo_sbs_14_img_0?_encoding=UTF8&psc=1&refRID=SCAHGQ3C91S24W175RPV

Answer 1

しかし上２つのfor文が無く代わりに
i = blockIdx.y * blockDim.y + threadIdx.y;
j = blockIdx.x * blockDim.x + threadIdx.x;
が記述されています。
なぜこのような書き方ができるのでしょうか？

大量のCUDA-coreが総数l*n個動くから。てかそうすることで高速に演算してる。

しかし適切なグリッド数とブロック数はどのように求めればよいのでしょうか？
func<<<grid(dim3型), block(dim3型)>>>();

両者のx,yごとの積が省略されたloopの繰り返し回数に相当する。
適切な数はcase-by-case。CUDAのドキュメントを読むがよい。

[追記]

グリッド数とブロック数は

正しくは「グリッドあたりのブロック数 と ブロックあたりのスレッド数」な。

Comments

[退会済みユーザー]

ご回答ありがとうございました。