ضرب دو ماتریس با cuda

ضرب ماتریس - ماتریس در GPU با  کمک Nvidia CUDA

 
ضرب ماتریس ماتریس
 

قبل از شروع، بهتر است به طور خلاصه چگونگی محاسبه ی ضرب ماتریس-ماترریس به صورت خلاصه توضیح داده شود. بگذارید بگوییم ما دو ماتریس A و B داریم. فرض کنید A یک ماتریس n × m است، به این معنی که آن دارای n ردیف و m ستون است. همچنین فرض کنید که B ماتریس m × w است. نتیجه ضرب ( A * B که متفاوت از B * A! است) یک ماتریس n × w است که ما آن را Mمی نامیم. یعنی تعداد ردیف ها در ماتریس حاصل، برابر با تعداد ردیف های ماتریس اول A و تعداد ستون ها برابر با ماتریس دوم B است.

  

چرا این اتفاق می افتد و چگونه کار می کند؟ پاسخ برای هر دو سوال در اینجا یکسان است. بیایید به سلول 1،1 (ردیف اول، ستون اول) از ماتریس M نگاه بندازیم. عدد داخل آن بعد از عملیات M = A * Bبرابر است با:  مجموع تمام عناصر عدد در A، ردیف 1، با عناصر عدد در B، ستون 1.  یعنی، در سلول  i,j، ازماتریس M ما مجموع سلول از تمام اعداد در ردیف iامین در ماتریس  A و ستون jامین در B داریم.

شکل زیر این ایده را توضیح می دهد:

اکنون باید بسیار واضح باشد چرا که ضرب ماتریس- ماتریس نمونه خوبی برای محاسبات موازی است. ما باید هر سلول را در c محاسبه کنیم، و هر کدام از آن ها مستقل از دیگری است، بنابراین ما می توانیم موازی سازی را به صورت موثری اجرا کنیم.

Grids And Blocks

هنگامی که ما یک kernel را که از دستورالعمل <<< >>> استفاده کرده است را فراخوانی می کنیم به طور خودکار یک متغیر از نوع سه بعدی(dim3) تعریف می کنیم  تعداد بلاک ها (blocks) در هر گرید (grid) و همچنین تعداد ترد ها (threads) در هر بلاک (block) تعریف می شود.

در حقیقت، grid ها و block ها آرایه ی سه بعدی به ترتیب از بلاک ها و ترد ها هستند. ما آن ها را قبل از فراخوانی کرنل(kernel) تعریف می کنیم:

dim3 blocksPerGrid(512, 1, 1)

dim3 threadsPerBlock(512, 1, 1)

kernel<<<blocksPerGrid, threadsPerBlock>>>()

 

 

در این روش، ما می توانیم به هر دو محور x و y به طریق مشابه مراجعه کنیم:

int row = blockIdx.y * blockDim.y + threadIdx.y;

int col = blockIdx.x * blockDim.x + threadIdx.x;

 

همان طور که می بینید، کد برای هر دو ی آن ها مشابه است. در کودا(cuda)، blockIdx، blockDim و threadIdx در تابع هایی با کمک مولفه های x،y  و z ساخته می شوند. آن ها مانند بردار های معمولی شماره گذاری می شوند، بنابراین بین 0 و بیشترین عدد منهای 1. برای مثال، اگر ما یک grid با مشخصات blocksPerGrid = (512, 1, 1) داشته باشیم، blockIdx.x بین 0 و 511 خواهد بود.

مجموع مقدار thread ها در یک block بیشتر از 1024 تجاوز نمی کند.

استفاده از یک block چند بعدی به این معنی است که شما باید در مورد توزیع این تعداد thread در میان تمام ابعاد مراقب باشید. در یک block یک بعدی، شما می توانید حداکثر 1024 thread را در محور x تنظیم کنید، اما در block دو بعدی، اگر شما 2 را به عنوان اندازه ی y تنظیم کنید، نمی توانید بیش از 512 برای  x در نظر بگیرید! برای مثال، dim3 threadsPerBlock(1024, 1, 1) مجاز است، و همین طور dim3 threadsPerBlock(512, 2, 1)، ولی نه dim3 threadsPerBlock(256, 3, 2).

آرایه های چند بعدی به صورت خطی

در این مقاله ما از آرایه های یک بعدی برای ماتریس های خود استفاده می کنیم. این ممکن است به نظر کمی گیج کننده باشد، اما خود  در زبان برنامه نویسی مسئله است. استانداردی که  CUDAبر اساس آن توسعه داده شده است احتیاج به دانشتن تعداد ستون ها قبل از کامپایل برنامه است. از این رو غیر ممکن است که آن را تغییر دهید یا آن را در وسط کد قرار دهید.

با این حال، کمی فکر نشان می دهد که این یک مسئله بزرگ نیست. ما نمیتوانیم از علامتگذاری  A [i] [j] به راحتی استفاده کنیم، اما ما زیاد تقلا نخواهیم کرد، همانطور که قبلا که به درستی شماره ی ردیفها و ستونها را می دانیم.

در حقیقت، ساده ترین راه برای خطی کردن یک آرایه ی دو بعدی این است که هر ردیف به صورت طولی، از اول تا آخرین، پشته شود. تصویر زیر این مفهوم را روشن تر می کند:

هسته(Kernel)

 

حالا که ما تمام اطلاعات لازم برای انجام این کار را داریم، بیایید نگاهی به کد kernel  هستیم. به خاطر ساده بودن ما در نمونه ما از ماتریس مربعی N × N استفاده می کنیم.

اولین چیزی که باید انجام داد، همانطور که قبلا مشاهده کردیم، تعیین شاخص محور x و y (به عنوان مثال شماره ردیف و ستون) است:

__global__ void multiplication(float *A, float* B, float *C, int N){

   int ROW = blockIdx.y*blockDim.y+threadIdx.y;

   int COL = blockIdx.x*blockDim.x+threadIdx.x;

 

سپس، ما بررسی می کنیم که مجموع ردیف ها و ستون ها از تعداد ردیف ها و ستون های واقعی در ماتریس تجاوز نمی کند. چون thread ها به حافظه به صورت تصادفی دسترسی پیدا می کنیم، ما باید این کار را برای جلوگیری از thread ها ی غیر ضروری از انجام عملیات بر روی ماتریس هایمان انجام دهیم. به این معنا، ما block ها و grid هایی با یک اندازه خاص را می سازیم، به طوری که اگر ما N را چند برابر آن اندازه قرار ندهیم، thread های بیشتری از نیازمان خواهیم داشت:

if (ROW < N && COL < N) {

 

در اینجا ما هیچ مشکلی با استفاده از ماتریسهای مربعی نداریم، اما همیشه این ایده خوبی است که در ذهن داشته باشیم. سپس ما باید متغیر موقت tmp_sum را برای جمع تمام سلول ها در ردیف و ستون انتخاب شده را مقدار دهی اولیه کنیم. همیشه تعیین نقاط اعشار و f، حتی اگر آن ها صفر باشند یک روش خوب است. بنابراین بنویسید:

 

  float tmp_sum = 0.0f;

 

حالا بخش ساده: به عنوان کد CPU، ما می توانیم از حلقه ی for برای محاسبه جمع استفاده کنیم و سپس آن را در سلول مربوطه C ذخیره کنیم.

 

     for (int i = 0; i < N; i++) {

            tmpSum += A[ROW * N + i] * B[i * N + COL];

        }

 

    C[ROW * N + COL] = tmpSum;

 

کد میزبان (CPU) با اعلان متغیرها به این صورت شروع می شود:

 

int main() {

    // declare arrays and variables

   //تعریف آرایه و متغیر ها

    int N = 16;

    int SIZE = N*N;

    vector<float> h_A(SIZE);

    vector<float> h_B(SIZE);

    vector<float> h_C(SIZE);

 

توجه داشته باشید که با توجه به شرط if در kernel، می توانیم مقادیر N را طوری تنظیم کنیم که لزوما چند برابر اندازه بلوک نباشد. همچنین، ما ازکتابخانه dev_array استفاده خواهیم کرد.

حالا اجازه دهید تا ماتریس ها را پر کنیم، چند راه برای انجام این کار وجود دارد، مانند ایجاد توابع برای ورودی دستی یا استفاده از اعداد تصادفی. در این مورد، ما به سادگی از یک حلقه برای پر کردن سلول ها با مقادیر مثلثاتی شاخص ها استفاده می کنیم:

for (int i=0; i<N; i++){

        for (int j=0; j<N; j++){

            h_A[i*N+j] = sin(i);

            h_B[i*N+j] = cos(j);

        }

    }

 

همان طور که در بخش "Grids and Blocks" اشاره شد، ما اکنون باید تنظیم متغیرهای dim3 را برای ابعاد بلوک ها و شبکه ها را بدانیم. grid فقط grid BLOCK_SIZE × BLOCK_SIZE است، بنابراین ما می توانیم بنویسیم:

dim3 threadsPerBlock (BLOCK_SIZE, BLOCK_SIZE)

 

همانطور که ما فقط با ماتریس هایی با اندازه های مختلف BLOCK_SIZE کار نمی کنیم، ما باید از دستورالعمل ceil استفاده کنیم تا شماره عدد صحیح بعدی را به عنوان اندازه ما بدست آوریم، همانطور که می بینید:

 

int n_blocks = ceil(N/BLOCK_SIZE);

dim3 blocksPerGrid (n_blocks, n_blocks)

 

توجه داشته باشید که به این ترتیب، از thread بیشتر از ضروریات استفاده خواهیم کرد، اما ما می توانیم با دستور شرطی  ifکه در kernel نوشتیم از وقوع آن در ماتریس جلوگیری کنیم.

این یک راه حل کلی که استدلال در پشت آن وجود دارد، اما در کد کامل شما یک نسخه کارآمد تر در kernel.cu پیدا می کنید.

اکنون ما فقط باید آرایه های دستگاه را ایجاد کنیم، حافظه را بر روی دستگاه اختصاص دهیم و kernel خود را فراخوانی می کنیم و در نتیجه ما یک برنامه ضرب ماتریس موازی داریم. با استفاده از dev_array، ما به سادگی می نویسیم:

    dev_array<float> d_A(SIZE);

    dev_array<float> d_B(SIZE);

    dev_array<float> d_C(SIZE);

 

برای اعلام آرایه ها و:

    d_A.set(&h_A[0], SIZE);

    d_B.set(&h_B[0], SIZE);

 

برای اختصاص حافظه. برای گرفتن نتیجه دستگاه و کپی کردن آن بر روی میزبان، ما از متد get استفاده می کنیم. یک بار دیگر، به این سادگی است:

    d_C.get(&h_C[0], SIZE);

 

در پایین این صفحه می توانید کد کامل را مشاهده کنید، از جمله مقایسه عملکرد و محاسبات خطا بین کد موازی و سریال.

کد کامل

dev_array.h

#ifndef _DEV_ARRAY_H_

#define _DEV_ARRAY_H_

 

#include <stdexcept>

#include <algorithm>

#include <cuda_runtime.h>

 

template <class T>

class dev_array

{

// public functions

public:

    explicit dev_array()

        : start_(0),

          end_(0)

    {}

 

    // constructor

    explicit dev_array(size_t size)

    {

        allocate(size);

    }

    // destructor

    ~dev_array()

    {

        free();

    }

 

    // resize the vector

    void resize(size_t size)

    {

        free();

        allocate(size);

    }

 

    // get the size of the array

    size_t getSize() const

    {

        return end_ - start_;

    }

 

    // get data

    const T* getData() const

    {

        return start_;

    }

 

    T* getData()

    {

        return start_;

    }

 

    // set

    void set(const T* src, size_t size)

    {

        size_t min = std::min(size, getSize());

        cudaError_t result = cudaMemcpy(start_, src, min * sizeof(T), cudaMemcpyHostToDevice);

        if (result != cudaSuccess)

        {

            throw std::runtime_error("failed to copy to device memory");

        }

    }

    // get

    void get(T* dest, size_t size)

    {

        size_t min = std::min(size, getSize());

        cudaError_t result = cudaMemcpy(dest, start_, min * sizeof(T), cudaMemcpyDeviceToHost);

        if (result != cudaSuccess)

        {

            throw std::runtime_error("failed to copy to host memory");

        }

    }

 

 

// private functions

private:

    // allocate memory on the device

    void allocate(size_t size)

    {

        cudaError_t result = cudaMalloc((void**)&start_, size * sizeof(T));

        if (result != cudaSuccess)

        {

            start_ = end_ = 0;

            throw std::runtime_error("failed to allocate device memory");

        }

        end_ = start_ + size;

    }

 

    // free memory on the device

    void free()

    {

        if (start_ != 0)

        {

            cudaFree(start_);

            start_ = end_ = 0;

        }

    }

 

    T* start_;

    T* end_;

};

 

#endif

 

 

matrixmul.cu

#include <iostream>

#include <vector>

#include <stdlib.h>

#include <time.h>

#include <cuda_runtime.h>

#include "kernel.h"

#include "kernel.cu"

#include "dev_array.h"

#include <math.h>

 

using namespace std;

 

int main()

{

    // Perform matrix multiplication C = A*B

    // where A, B and C are NxN matrices

    int N = 16;

    int SIZE = N*N;

 

    // Allocate memory on the host

    vector<float> h_A(SIZE);

    vector<float> h_B(SIZE);

    vector<float> h_C(SIZE);

 

    // Initialize matrices on the host

    for (int i=0; i<N; i++){

        for (int j=0; j<N; j++){

            h_A[i*N+j] = sin(i);

            h_B[i*N+j] = cos(j);

        }

    }

 

    // Allocate memory on the device

    dev_array<float> d_A(SIZE);

    dev_array<float> d_B(SIZE);

    dev_array<float> d_C(SIZE);

 

    d_A.set(&h_A[0], SIZE);

    d_B.set(&h_B[0], SIZE);

 

    matrixMultiplication(d_A.getData(), d_B.getData(), d_C.getData(), N);

    cudaDeviceSynchronize();

 

    d_C.get(&h_C[0], SIZE);

    cudaDeviceSynchronize();

 

    float *cpu_C;

    cpu_C=new float[SIZE];

 

    // Now do the matrix multiplication on the CPU

    float sum;

    for (int row=0; row<N; row++){

        for (int col=0; col<N; col++){

            sum = 0.f;

            for (int n=0; n<N; n++){

                sum += h_A[row*N+n]*h_B[n*N+col];

            }

            cpu_C[row*N+col] = sum;

        }

    }

 

    double err = 0;

    // Check the result and make sure it is correct

    for (int ROW=0; ROW < N; ROW++){

        for (int COL=0; COL < N; COL++){

            err += cpu_C[ROW * N + COL] - h_C[ROW * N + COL];

        }

    }

 

    cout << "Error: " << err << endl;

 

    return 0;

}

 

kernel.h

#ifndef KERNEL_CUH_

#define KERNEL_CUH_

 

void matrixMultiplication(float *A, float *B, float *C, int N);

 

#endif

kernel.cu

#include <math.h>

#include <iostream>

#include "cuda_runtime.h"

#include "kernel.h"

#include <stdlib.h>

 

using namespace std;

 

__global__ void matrixMultiplicationKernel(float* A, float* B, float* C, int N) {

 

    int ROW = blockIdx.y*blockDim.y+threadIdx.y;

    int COL = blockIdx.x*blockDim.x+threadIdx.x;

 

    float tmpSum = 0;

 

    if (ROW < N && COL < N) {

        // each thread computes one element of the block sub-matrix

        for (int i = 0; i < N; i++) {

            tmpSum += A[ROW * N + i] * B[i * N + COL];

        }

    }

    C[ROW * N + COL] = tmpSum;

}

 

 

void matrixMultiplication(float *A, float *B, float *C, int N){

 

    // declare the number of blocks per grid and the number of threads per block

    // use 1 to 512 threads per block

    dim3 threadsPerBlock(N, N);

    dim3 blocksPerGrid(1, 1);

        if (N*N > 512){

            threadsPerBlock.x = 512;

            threadsPerBlock.y = 512;

            blocksPerGrid.x = ceil(double(N)/double(threadsPerBlock.x));

            blocksPerGrid.y = ceil(double(N)/double(threadsPerBlock.y));

        }

 

    matrixMultiplicationKernel<<<blocksPerGrid,threadsPerBlock>>>(A, B, C, N);

}