CUDA-8

CUDA-8

这次把julia set 实现一下,再练习一下cuda相关的logic.

julia set 是复变函数中一个非常有意思的集合,所以应该先写一个关于复数的结构

struct cuComplex{
    float r;
    float i;
    
    cuComplex(float a, float b): r(a), i(b){}
    float magnitude2(void){return r*r+i*i;}   //返回其模长的平方
    cuComplex operator*(const cuComplex& a){
        return cuComplex(r*a.r-i*a.i, i*a.r+r*a.i);
    }
    cuComplex operator+(const cuComplex& a){
        return cuComplex(r+a.r, i+a.i);

    }

};

上面中又重写了加法和×法。

然后就是julia set

int julia(int x, int y){
    const float scale=1.5;
    float jx = scale*(float)(DIM/2-x)/(DIM/2);
    float jy = scale*(float)(DIM/2-y)/(DIM/2);
    
    cuComplex c(-0.8, 0.156);
    cuComplex a(jx,jy);

    int i=0;
    for(i=0; i<200; i++){
        a = a*a +c;
        if(a.magnitude2()>1000)   // divergence
            return 0;
    }
    return 1;   //convergent

}

上面是julia集的定义,这里iter了200次。

下面是主要的核函数(cpu)

void kernel(unsigned char *ptr){
    for(int y=0; y<DIM;y++){
        for(int x=0;x<DIM;x++){
            int offset = x+y*DIM;
            int juliaValue = julia(x,y);
            ptr[offset*4+0] = 255*juliaValue;
            ptr[offset*4+1] = 0;
            ptr[offset*4+2] = 0;
            ptr[offset*4+3] = 255;
        }

    }

}

然后是main函数

int main(void){
    CPUBitmap bitmap(DIM, DIM);
    unsigned char *ptr = bitmap.get_ptr();
    kernel(ptr);
    bitmap.display_and_exit();
}

结果如图

avator

下面是用gpu写的,可以体会一下差别,相当于以前是一个人来干活,现在有很多人来一起干的感觉。

#include "../common/book.h"
#include "../common/cpu_bitmap.h"
#define DIM 1000

struct cuComplex {
    float   r;
    float   i;
// 指定这些运算只能在device 上运行。
__device__ cuComplex( float a, float b ) : r(a), i(b)  {}
    __device__ float magnitude2( void ) {
        return r * r + i * i;
    }
    __device__ cuComplex operator*(const cuComplex& a) {
        return cuComplex(r*a.r - i*a.i, i*a.r + r*a.i);
    }
    __device__ cuComplex operator+(const cuComplex& a) {
        return cuComplex(r+a.r, i+a.i);
    }
};

__device__ int julia( int x, int y ) {
    const float scale = 1.5;
    float jx = scale * (float)(DIM/2 - x)/(DIM/2);
    float jy = scale * (float)(DIM/2 - y)/(DIM/2);

    cuComplex c(-0.8, 0.156);
    cuComplex a(jx, jy);

    int i = 0;
    for (i=0; i<200; i++) {
        a = a * a + c;
        if (a.magnitude2() > 1000)
            return 0;
    }   

    return 1;
}

__global__ void kernel( unsigned char *ptr ) { 
    // map from blockIdx to pixel position
    // 因为图像可以理解成两维的,所以这里也用2维的, 可以把y理解成第几行,但是其实因为这个图是方的,其
实理解成列也没事儿。
    int x = blockIdx.x;
    int y = blockIdx.y;
    int offset = x + y * gridDim.x;

    // now calculate the value at that position
int juliaValue = julia( x, y );
    ptr[offset*4 + 0] = 255 * juliaValue;
    ptr[offset*4 + 1] = 0;
    ptr[offset*4 + 2] = 0;
    ptr[offset*4 + 3] = 255;
}

// globals needed by the update routine
struct DataBlock {
    unsigned char   *dev_bitmap;
};

int main( void ) {
    DataBlock   data;
    CPUBitmap bitmap( DIM, DIM, &data );
    unsigned char    *dev_bitmap;

    HANDLE_ERROR( cudaMalloc( (void**)&dev_bitmap, bitmap.image_size() ) );
    data.dev_bitmap = dev_bitmap;

    dim3    grid(DIM,DIM);
    kernel<<<grid,1>>>( dev_bitmap );

    HANDLE_ERROR( cudaMemcpy( bitmap.get_ptr(), dev_bitmap,
                              bitmap.image_size(),
                              cudaMemcpyDeviceToHost ) );

    HANDLE_ERROR( cudaFree( dev_bitmap ) );

    bitmap.display_and_exit();
}


差别在于以前在host上运行的,现在要把这些运算放在device上去算,还有最后那个开blocks的地方。 作者直接开了一个dims grid(DIM, DIM), 我试过如果弄成

dims block(DIM, DIM), 然后把index 换成

int x = threadIdx.x;
int y = threadIdx.y;
int offset = x + y*blockDim.x;

<<<1, block>>>

结果却不对,不知道是什么原因,按书里面的话,说明是把代码copy了DIM*DIM份。而改了之后相当于

的话,相当于就没有copy,而只是在一个block中运行的,难道是有冲突了吗?

打赏,谢谢~~

取消

感谢您的支持,我会继续努力的!

扫码支持
扫码打赏,多谢支持~

打开微信扫一扫,即可进行扫码打赏哦