Julia曲线的生成

/*
生成Julia曲线，与GPU高性能编程CUDA实战这本书相对应，
其中用到了作者提供的与书配合的库：book.h,cpu_bitmap.h

本程序是在cuda3.0上调通的。
 */

#include "book.h"
#include "cpu_bitmap.h"

#define DIM 1000
struct cuComplex{
 float r;
 float i;
 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);      //Zc为复数
 //cuComplex c(-0.8,0);      //Zc为实数
 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)
{
 //将threadIdx/BlockIdx映射到像素位置
 int x=blockIdx.x;
 int y=blockIdx.y;
 int offset=x+y*gridDim.x;
 //顼在计算这个位置上的值
 int juliaValue=julia(x,y);
 //每个颜色用四个字节表示，给每个字节赋值
 ptr[offset*4+0]=255*juliaValue;  //确定图像的颜色的Red分量
 ptr[offset*4+1]=0;   //确定图像的颜色的Green分量
 ptr[offset*4+2]=0;  //确定图像的颜色的Blue分量
 ptr[offset*4+3]=255;
}
int main(void)
{
 CPUBitmap bitmap(DIM,DIM);
 unsigned char *dev_bitmap;

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

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

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