Блог

Я хочу вычислить ‘out = alpha * px beta * py’, ‘px’ и ‘py’ — это массив.*

У меня простое ядро:

 __global__ 
void saxpyGPU2( float *out, const float *px, const float *py, size_t N, float alpha,float beta )
{
    size_t i = blockDim.x*blockIdx.x   threadIdx.x;
    while (i < N)
    {
        out[i] = alpha * px[i]   beta * py[i];
        i  = blockDim.x*gridDim.x;
    } 
}
 

Это работает, поэтому я хочу развернуть цикл.

Код в cuda-handbook:

 template<const int n> 
__device__ 
void saxpy_unrolled(float *out, const float *px, const float *py, size_t N, float alpha,float beta)
{
    float x[n], y[n];
    size_t i;
    for ( i = n*blockIdx.x*blockDim.x threadIdx.x; i < N-n*blockDim.x*gridDim.x; i  = n*blockDim.x*gridDim.x ) {
        for ( int j = 0; j < n; j   ) {
            size_t index = i j*blockDim.x;
            x[j] = px[index];
            y[j] = py[index];
        }
        for ( int j = 0; j < n; j   ) {
            size_t index = i j*blockDim.x;
            out[index] = alpha*x[j] beta* y[j];
        }
    }
    // to avoid the (index<N) conditional in the inner loop, 
    // we left off some work at the end
    for ( int j = 0; j < n; j   ) {
        for ( int j = 0; j < n; j   ) {
            size_t index = i j*blockDim.x;
            if ( index<N ) {
                x[j] = px[index];
                y[j] = py[index];
            }
        }
        for ( int j = 0; j < n; j   ) {
            size_t index = i j*blockDim.x;
            if ( index<N ) out[index] = alpha*x[j] beta* y[j];
        }
    }
}

__global__ 
void saxpyGPU( float *out, const float *px, const float *py, size_t N, float alpha,float beta )
{
    saxpy_unrolled<4>( out, px, py, N, alpha ,beta);
}
 

Я не понимаю во второй ветке, когда i> N-n * blockDim.x * gridDim.x. зачем использовать внешний цикл

for ( int j = 0; j < n; j ) {
for ( int j = 0; j < n; j )....}

И я тестирую эти два ядра, первое в порядке, но второе, которое я копирую из книги, неверно.

Я создаю два массива while(i<1024) a[i] = i; b[i] = 10*i;i , и я хочу вычислить c = alpha * a beta * b, используя два ядра выше, но результат в развернутом ядре цикла равен 4.3e8 для всех элементов в c.

Это мой тестовый код:

 int main(){
int arraySize = 1024;
float* a =new float[arraySize];
float* b =new float[arraySize];
float* c =new float[arraySize];
for (int i =0;i<arraySize;i  )
{
    a[i] = 1.0* i;
    b[i] = 10.0*i;
    c[i] = 0.0;
}
float* d_a;
float* d_b;
float* d_c;
cudaMalloc((void**)amp;d_a,sizeof(float)*arraySize);
cudaMemcpy(d_a,a,sizeof(float)*arraySize,cudaMemcpyHostToDevice);
cudaMalloc((void**)amp;d_b,sizeof(float)*arraySize);
cudaMemcpy(d_b,b,sizeof(float)*arraySize,cudaMemcpyHostToDevice);
cudaMalloc((void**)amp;d_c,sizeof(float)*arraySize);
for (int i=0;i<arraySize;i  )
{
    c[i] = a[i]   b[i];
}
dim3 block_size(256,1,1);
dim3 grid_size((arraySize -1)/block_size.x 1,1,1);
float alpha = 1.0;
float beta = 1.0;
bool flag = true;
if(flag)
{
    saxpyGPU<<<grid_size,block_size>>>(d_c,d_a,d_b,arraySize,alpha,beta);
    float* temp = new float[arraySize];
    cudaMemcpy(temp,d_c,arraySize*sizeof(float),cudaMemcpyDeviceToHost);
    for (int i = 0;i<arraySize;i  )
    {
        cout<<(temp[i] - c[i])<<",";
    }
}
else
{
    saxpyGPU2<<<grid_size,block_size>>>(d_c,d_a,d_b,arraySize,alpha,beta);
    cudaMemcpy(temp,d_c,arraySize*sizeof(float),cudaMemcpyDeviceToHost);
    for (int i = 0;i<arraySize;i  )
    {
        cout<<(temp[i] - c[i])<<",";
    }
 

Эти два ядра показывают разный результат

Опубликованный вами код ядра абсолютно корректен и дает ожидаемые результаты. Это можно продемонстрировать с помощью следующего кода:

 #include <thrust/random.h>
#include <thrust/device_vector.h>
#include <thrust/transform.h>
#include <thrust/copy.h>
#include <thrust/iterator/counting_iterator.h>
#include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>

template<const int n> 
__device__ 
void saxpy_unrolled(float *out, const float *px, const float *py, 
        size_t N, float alpha,float beta) {
    float x[n], y[n];
    size_t i;
    for ( i = n*blockIdx.x*blockDim.x threadIdx.x; 
        i < N-n*blockDim.x*gridDim.x; 
        i  = n*blockDim.x*gridDim.x ) {
        for ( int j = 0; j < n; j   ) {
            size_t index = i j*blockDim.x;
            x[j] = px[index];
            y[j] = py[index];
        }
        for ( int j = 0; j < n; j   ) {
            size_t index = i j*blockDim.x;
            out[index] = alpha*x[j] beta* y[j];
        }
    }
    for ( int j = 0; j < n; j   ) {
        for ( int j = 0; j < n; j   ) {
            size_t index = i j*blockDim.x;
            if ( index<N ) {
                x[j] = px[index];
                y[j] = py[index];
            }
        }
        for ( int j = 0; j < n; j   ) {
            size_t index = i j*blockDim.x;
            if ( index<N ) {
                out[index] = alpha*x[j]   beta*y[j];
            }
        }
    }
}

__global__ 
void saxpyGPU( float *out, const float *px, const float *py, 
        size_t N, float alpha,float beta ) {
    saxpy_unrolled<4>( out, px, py, N, alpha ,beta);
}

struct prg {
    float a, b;
    __host__ __device__
    prg(float _a=0.f, float _b=1.f) : a(_a), b(_b) {};

    __host__ __device__
    float operator()(const unsigned int n) const {
        thrust::default_random_engine rng;
        thrust::uniform_real_distribution<float> dist(a, b);
        rng.discard(n);
        return dist(rng);
    }
};

int main(void) {
    const int N = 100000;
    const float alpha = 0.12345f, beta = 0.9876f;

    prg gen(1.f, 2.f);
    thrust::device_vector<float> x(N), y(N), z(N);
    thrust::counting_iterator<unsigned int> iseqx(0);
    thrust::counting_iterator<unsigned int> iseqy(N);
    thrust::transform(iseqx, iseqx   N, x.begin(), gen);
    thrust::transform(iseqy, iseqy   N, y.begin(), gen);

    float *xp = thrust::raw_pointer_cast(amp;x[0]);
    float *yp = thrust::raw_pointer_cast(amp;y[0]);
    float *zp = thrust::raw_pointer_cast(amp;z[0]);
    dim3 blockdim(128);
    dim3 griddim(16);
    saxpyGPU<<<griddim, blockdim>>>(zp, xp, yp, N, alpha, beta);
    cudaDeviceSynchronize();

    std::vector<float> xh(N), yh(N), zh(N);
    thrust::copy(x.begin(), x.end(), xh.begin());
    thrust::copy(y.begin(), y.end(), yh.begin());
    thrust::copy(z.begin(), z.end(), zh.begin());

    float maxabserr = -1.f, maxrelerr = -1.f;
    for(int i=0; i<N; i  ) {
        float saxpyval = alpha * xh[i]   beta * yh[i];
        float abserr = fabs(zh[i]-saxpyval);
        float relerr = abserr / fmaxf(fabs(zh[i]), fabs(saxpyval));
        maxabserr = fmaxf(abserr, maxabserr);
        maxrelerr = fmaxf(relerr, maxrelerr);
    }
    std::cout.precision(10);
    std::cout << "Maximum absolute error = " << maxabserr << std::endl;
    std::cout << "Maximum relative error = " << maxrelerr << std::endl;

    return 0;
}
 

что дает мне следующее:

 $ nvcc -arch=sm_30 -o unrolled_saxpy unrolled_saxpy.cu
$ ./unrolled_saxpy 
Maximum absolute error = 2.384185791e-07
Maximum relative error = 1.1920676e-07
 

Если вы (все еще) не понимаете, почему ядро написано так, как оно есть, следуйте тому, что я показал вам в вашем предыдущем вопросе, и вручную разверните функцию saxpy. Начните с n = 1 и убедитесь, что он функционально совпадает с развернутым эквивалентом, а затем попробуйте n = 2, n = 4 и т. Д. чтобы увидеть, каково действие развертывания цикла.