[菜鳥每天來段CUDA_C]基于GPU的Julia集

Julia集是一個著名的分形集，它是複數經過疊代得到的，是滿足某個複數計算函數的所有點構成的邊界。 算法思想： 通過一個簡單的疊代等式對複平面中的點求值，如果在計算某個點時疊代的結果是發散的，那麼這個點就不屬于Julia集。相反，如果疊代計算得到的一系列值都位于某個邊界範圍之内，那麼這個點就屬于Julia集。疊代計算的公式為： Z(n+1) = Z(n)^2 + C  (C為常數) 基于GPU的Julia集通過DIM*DIM個thread計算平面中各點的值，最終用OpenGL顯示計算的結果如下：

主要程式代碼：

#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>
#include <cutil_inline.h>
#include "CPUBitmap.h"
//#include <device_launch_parameters.h>

#define DIM 512

struct myComplex{
	float r;
	float i;
	__device__ myComplex(float a, float b): r(a), i(b){}

	__device__ float magnitude2(){
		return r*r+i*i;
	}

	__device__ myComplex operator*(const myComplex& a){
		return myComplex(r*a.r - i*a.i, i*a.r + r*a.i);
	}

	__device__ myComplex operator+(const myComplex& a){
		return myComplex(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);

	myComplex c(-0.8, 0.156);
	myComplex 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)
{
	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;
	ptr[offset*4 + 1] = 255 *juliaValue;
	ptr[offset*4 + 2] = 255 *juliaValue;
	ptr[offset*4 + 3] = 255;
}

int main(int argc, char* argv[])
{
	if(!InitCUDA()) {
		return 0;
	}

	CPUBitmap bitmap(DIM, DIM);

	unsigned char *dev_bitmap;
	cutilSafeCall(cudaMalloc((void**)&dev_bitmap, bitmap.image_size()));
	dim3 grid(DIM, DIM);
	kernel<<<grid, 1>>>(dev_bitmap);

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

	bitmap.display_and_exit();

	cutilSafeCall(cudaFree(dev_bitmap));

	return 0;
}