resetreboot/mandelbrot_study
1
0
Fork 0
Code Issues Pull requests Releases Wiki Activity

First commit

Browse source
This commit is contained in:
José Carlos Cuevas 2013-08-29 19:31:20 +02:00
commit a7db7bced3
18 changed files with 1106 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

BIN
.Makefile.swp Normal file
View file

Binary file not shown.

BIN
.main.c.swp Normal file
View file

Binary file not shown.

BIN
.mandel_classic.c.swp Normal file
View file

Binary file not shown.

BIN
.mandelbrot_kernel.cl.swp Normal file
View file

Binary file not shown.

BIN
.sum_kernel.cl.swp Normal file
View file

Binary file not shown.

BIN
.vector_add_kernel.cl.swp Normal file
View file

Binary file not shown.

1
CL Symbolic link
View file

@ -0,0 +1 @@
/opt/intel/opencl-1.2-3.0.67279/include/CL

38
Makefile Normal file
View file

@ -0,0 +1,38 @@
# MandelCL Makefile
CC=gcc
# Flags!
SDLFLAGS=$(shell sdl-config --cflags)
# Comment this line and uncomment the next to get Julia fractals
CFLAGS=-c -Wall -O2 $(SDLFLAGS)
# CFLAGS=-c -Wall -O2 -DJULIA $(SDLFLAGS)
# CFLAGS=-c -Wall -ggdb $(SDLFLAGS)
# Libs!
SDLLIBS=$(shell sdl-config --libs)
LIBS=-lm -lpthread $(SDLLIBS)
# Includes!
INCLUDE=-I/usr/include/SDL
all: mandelclassic test
mandelclassic: mandel_classic.o
$(CC) $(INCLUDE) mandel_classic.o $(LIBS) -o mandelclassic
mandelclassic.o: mandel_classic.c
$(CC) $(CFLAGS) $(INCLUDE) $(LIBS) mandel_classic.c -o mandel_classic.o
test: test.o
$(CC) $(INCLUDE) test.o $(LIBS) -o test
test.o: test.c
$(CC) $(CFLAGS) $(INCLUDE) $(LIBS) test.c -o test.o
.PHONY: clean
clean:
@rm *.o mandelclassic test

38
frag_shader.glsl Normal file
View file

@ -0,0 +1,38 @@
void main(void)
{
float zoom = float(iMouse.x) / float(iResolution.x);
float pos_x = ((float(gl_FragCoord.x) / float(iResolution.x)) * 3.5 * zoom) - (2.5 - (1.0 - zoom));
float pos_y = ((float(gl_FragCoord.y) / float(iResolution.y)) * 2.0 * zoom) - (1.0 - (1.0 - zoom));
float x = 0.0;
float y = 0.0;
int iteration = 0;
float normal_iter = 0.0;
int max_iteration = 255;
float xtemp;
while (iteration < max_iteration)
{
xtemp = x * x - y * y + pos_x;
y = 2.0 * x * y + pos_y;
if ((x * x) + (y * y) >= (4.0)) break;
x = xtemp;
iteration++;
}
if (iteration >= max_iteration)
{
gl_FragColor = vec4(0.0,0.0,0.0,1.0);
}
else if (iteration < 128)
{
normal_iter = float(iteration) / 255.0;
gl_FragColor = vec4(0,0.1 + normal_iter,normal_iter,1.0);
}
else
{
normal_iter = float(iteration) / 255.0;
gl_FragColor = vec4(1.0 - normal_iter, 1.0, 1.0 - normal_iter, 1.0);
}
}

209
main.c Normal file
View file

@ -0,0 +1,209 @@
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#ifdef __APPLE__
#include <OpenCL/opencl.h>
#else
#include <CL/cl.h>
#endif
#include <SDL.h>
#define MAX_SOURCE_SIZE (0x100000)
int main(int argn, char **argv) {
// Init SDL
if(SDL_Init(SDL_INIT_VIDEO) != 0)
fprintf(stderr, "Could not initialize SDL: %s\n", SDL_GetError());
printf("SDL Initialized\n");
// Create screen surface
SDL_Surface *screen;
int res_x = 800;
int res_y = 600;
int current_line = 0;
int total_res = res_x * res_y;
screen = SDL_SetVideoMode(res_x, res_y, 0, SDL_HWSURFACE|SDL_DOUBLEBUF);
if(!screen)
fprintf(stderr,"Could not set video mode: %s\n",SDL_GetError());
// Prepare the resolution and sizes and colors...
int i;
int temp;
const int ITERATIONS = 256;
int *red_scale = (int*)malloc(sizeof(int)*ITERATIONS);
int *blue_scale = (int*)malloc(sizeof(int)*ITERATIONS);
for(i = 0; i < ITERATIONS; i++) {
red_scale[i] = i;
blue_scale[i] = 255 - i;
}
// Load the kernel source code into the array source_str
FILE *fp;
char *source_str;
size_t source_size;
fp = fopen("mandelbrot_kernel.cl", "r");
if (!fp) {
fprintf(stderr, "Failed to load kernel.\n");
exit(1);
}
source_str = (char*)malloc(MAX_SOURCE_SIZE);
source_size = fread( source_str, 1, MAX_SOURCE_SIZE, fp);
fclose( fp );
// Get platform and device information
cl_platform_id platform_id = NULL;
cl_device_id device_id = NULL;
cl_uint ret_num_devices;
cl_uint ret_num_platforms;
cl_int ret = clGetPlatformIDs(1, &platform_id, &ret_num_platforms);
ret = clGetDeviceIDs( platform_id, CL_DEVICE_TYPE_GPU, 1,
&device_id, &ret_num_devices);
// Create an OpenCL context
cl_context context = clCreateContext( NULL, 1, &device_id, NULL, NULL, &ret);
// Create a command queue
cl_command_queue command_queue = clCreateCommandQueue(context, device_id, 0, &ret);
// Create memory buffers on the device for returning iterations
cl_mem kernel_res_x = clCreateBuffer(context, CL_MEM_READ_ONLY,
sizeof(cl_int), NULL, &ret);
cl_mem kernel_res_y = clCreateBuffer(context, CL_MEM_READ_ONLY,
sizeof(cl_int), NULL, &ret);
cl_mem kernel_current_line = clCreateBuffer(context, CL_MEM_READ_ONLY,
sizeof(cl_int), NULL, &ret);
cl_mem graph_mem_obj = clCreateBuffer(context, CL_MEM_WRITE_ONLY,
res_x * sizeof(cl_int), NULL, &ret);
// Copy resolution x and y for the kernel
ret = clEnqueueWriteBuffer(command_queue, kernel_res_x, CL_TRUE, 0,
sizeof(cl_int), &res_x, 0, NULL, NULL);
ret = clEnqueueWriteBuffer(command_queue, kernel_res_y, CL_TRUE, 0,
sizeof(cl_int), &res_y, 0, NULL, NULL);
ret = clEnqueueWriteBuffer(command_queue, kernel_current_line, CL_TRUE, 0,
sizeof(cl_int), &current_line, 0, NULL, NULL);
// Create a program from the kernel source
cl_program program = clCreateProgramWithSource(context, 1,
(const char **)&source_str, (const size_t *)&source_size, &ret);
// Build the program
ret = clBuildProgram(program, 1, &device_id, NULL, NULL, NULL);
// Create the OpenCL kernel
cl_kernel kernel = clCreateKernel(program, "mandelbrot_point", &ret);
// Our screen in a linear array
int *graph_dots = (int*)malloc(total_res * sizeof(int));
cl_int *graph_line = (cl_int*)malloc(res_x * sizeof(cl_int));
for (current_line = 0; current_line < 600; current_line++)
{
// Set the arguments of the kernel
ret = clEnqueueWriteBuffer(command_queue, kernel_current_line, CL_TRUE, 0,
sizeof(cl_int), &current_line, 0, NULL, NULL);
ret = clSetKernelArg(kernel, 0, sizeof(cl_mem), &kernel_res_x);
ret = clSetKernelArg(kernel, 1, sizeof(cl_mem), &kernel_res_y);
ret = clSetKernelArg(kernel, 2, sizeof(cl_mem), &kernel_current_line);
ret = clSetKernelArg(kernel, 3, sizeof(cl_mem), graph_mem_obj);
// Execute the OpenCL kernel on the list
size_t global_item_size = res_x; // Process the entire screen
size_t local_item_size = 64; // Process in groups of 64
ret = clEnqueueNDRangeKernel(command_queue, kernel, 1, NULL,
&global_item_size, &local_item_size, 0, NULL, NULL);
// Read the memory buffer graph_mem_obj on the device to the local variable graph_dots
ret = clEnqueueReadBuffer(command_queue, graph_mem_obj, CL_TRUE, 0,
res_x * sizeof(cl_int), graph_line, 0, NULL, NULL);
for (i = 0; i < 800; i++)
{
graph_dots[(current_line * 800) + i] = graph_line[i];
}
}
// Display the result to the screen
/* for(i = 0; i < 3078; i++)
printf("Linear: %d -> %d\n", i, graph_dots[i]); */
printf("Rendering...\n");
int iteration;
Uint32 *pixel;
// Lock surface
SDL_LockSurface(screen);
// rank = screen->pitch/sizeof(Uint32);
pixel = (Uint32*)screen->pixels;
/* Draw all dots */
for(i = 0;i < total_res;i++)
{
// Get the iterations for the point
// printf("Point %d\n", i);
iteration = graph_dots[i];
if ((iteration < 1000) && (iteration >= 0)) {
pixel[i] = SDL_MapRGBA(screen->format,
red_scale[iteration],
0,
blue_scale[iteration],
255);
}
else
{
pixel[i] = SDL_MapRGBA(screen->format,
0,
0,
0,
255);
}
}
// Unlock surface
SDL_UnlockSurface(screen);
// Draw to the scree
SDL_Flip(screen);
// Clean up
ret = clFlush(command_queue);
ret = clFinish(command_queue);
ret = clReleaseKernel(kernel);
ret = clReleaseProgram(program);
// ret = clReleaseMemObject(a_mem_obj);
// ret = clReleaseMemObject(b_mem_obj);
ret = clReleaseMemObject(graph_mem_obj);
ret = clReleaseCommandQueue(command_queue);
ret = clReleaseContext(context);
// free(A);
// free(B);
free(graph_dots);
SDL_Event ev;
int active;
active = 1;
while(active)
{
/* Handle events */
while(SDL_PollEvent(&ev))
{
if(ev.type == SDL_QUIT)
active = 0; /* End */
}
}
SDL_Quit();
return 0;
}

462
mandel_classic.c Normal file
View file

@ -0,0 +1,462 @@
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <pthread.h>
#include <math.h>
#include <SDL.h>
#define MAX_SOURCE_SIZE (0x100000)
#ifdef CACHE
int** cached_points;
int** cached_x;
int** cached_y;
#endif
int *iteration_pixels;
typedef struct point_args point_args;
struct point_args
{
int res_x;
int res_y;
int image_x;
int image_y;
float zoom;
int max_iteration;
int thread_number;
};
typedef struct piece_args piece_args;
struct piece_args
{
int res_x;
int res_y;
float zoom;
int max_iteration;
int total_threads;
int thread_number;
};
int get_x (int linear_point, int width)
{
return linear_point % width;
}
int get_y (int linear_point, int height)
{
return floor(linear_point / height);
}
float map_x(int x, int width, float zoom)
{
#ifndef JULIA
return (((float)x / (float)width) * (3.5 * zoom)) - (2.5 - (1.0 - zoom));
#else
return (((float)x / (float)width) * (3.5 * zoom)) - (1.75 - (1.0 - zoom));
#endif
}
float map_y(int y, int height, float zoom)
{
return (((float)y / (float)height) * (2.0 * zoom)) - (1.00001 - (1.0 - zoom));
}
#ifdef CACHE
int cached_iteration(float pos_x, float pos_y)
{
float centered_x = pos_x + 2.5;
float centered_y = pos_y + 1.0;
float temp_x = floor(centered_x * 1000.0);
float temp_y = floor(centered_y * 1000.0);
int trs_pos_x = (int)temp_x;
int trs_pos_y = (int)temp_y;
return cached_points[trs_pos_x][trs_pos_y];
}
float get_cached_x(float pos_x, float pos_y)
{
float centered_x = pos_x + 2.5;
float centered_y = pos_y + 1.0;
float temp_x = floor(centered_x * 1000.0);
float temp_y = floor(centered_y * 1000.0);
int trs_pos_x = (int)temp_x;
int trs_pos_y = (int)temp_y;
return cached_x[trs_pos_x][trs_pos_y];
}
float get_cached_y(float pos_x, float pos_y)
{
float centered_x = pos_x + 2.5;
float centered_y = pos_y + 1.0;
float temp_x = floor(centered_x * 1000.0);
float temp_y = floor(centered_y * 1000.0);
int trs_pos_x = (int)temp_x;
int trs_pos_y = (int)temp_y;
return cached_y[trs_pos_x][trs_pos_y];
}
void store_iteration(float pos_x, float pos_y, int iteration, float x, float y)
{
float centered_x = pos_x + 2.5;
float centered_y = pos_y + 1.0;
float temp_x = floor(centered_x * 1000.0);
float temp_y = floor(centered_y * 1000.0);
int trs_pos_x = (int)temp_x;
int trs_pos_y = (int)temp_y;
cached_points[trs_pos_x][trs_pos_y] = iteration;
cached_x[trs_pos_x][trs_pos_y] = x;
cached_y[trs_pos_x][trs_pos_y] = y;
}
#endif
int mandelbrot_point(int res_x, int res_y, int image_x, int image_y, float zoom, int max_iteration)
{
// Get the index of the current element
float pos_x = map_x(image_x, res_x, zoom);
float pos_y = map_y(image_y, res_y, zoom);
float x = 0.0;
float y = 0.0;
float q, x_term;
float xtemp, xx, yy;
#ifdef CACHE
int storeable = 1;
#endif
int iteration = 0;
yy = y * y;
// Period-2 bulb check
if (((x + 1) * (x + 1) + yy) < 0.0625) return 0;
// Cardioid check
x_term = x - 0.25;
q = x_term * x_term + yy;
q = q * (q + x_term);
if (q > (0.25 * yy)) return 0;
#ifdef CACHE
// Look up our cache
iteration = cached_iteration(pos_x, pos_y);
if (iteration > 0)
{
x = get_cached_x(pos_x, pos_y);
y = get_cached_y(pos_x, pos_y);
yy = y * y;
}
if (iteration < 0) storeable = 0;
#endif
while (iteration < max_iteration)
{
xx = x * x;
if ((xx) + (yy) > (4.0)) break;
y = (x + y) * (x + y) - xx - yy;
y = y + pos_y;
xtemp = xx - yy + pos_x;
x = xtemp;
yy = y * y;
iteration++;
}
if (iteration >= max_iteration)
{
return 0;
}
else
{
#ifdef CACHE
if (storeable == 1)
{
store_iteration(pos_x, pos_y, iteration, x, y);
}
#endif
return iteration;
}
}
int julia_point(int res_x, int res_y, int image_x, int image_y, float zoom, int max_iteration)
{
// Get the index of the current element
float pos_x = map_x(image_x, res_x, 1.0);
float pos_y = map_y(image_y, res_y, 1.0);
float x = pos_x;
float y = pos_y;
float xtemp, xx, yy;
#ifdef CACHE
int storeable = 1;
#endif
int iteration = 0;
#ifdef CACHE
// Look up our cache
iteration = cached_iteration(pos_x, pos_y);
if (iteration > 0)
{
x = get_cached_x(pos_x, pos_y);
y = get_cached_y(pos_x, pos_y);
yy = y * y;
}
if (iteration < 0) storeable = 0;
#endif
while (iteration < max_iteration)
{
xx = x * x;
yy = y * y;
if ((xx) + (yy) > (4.0)) break;
y = pow((x + y), 2) - xx - yy;
y = y + 0.288;
xtemp = xx - yy + 0.353 + zoom;
x = xtemp;
iteration++;
}
if (iteration >= max_iteration)
{
return 0;
}
else
{
#ifdef CACHE
if (storeable == 1)
{
store_iteration(pos_x, pos_y, iteration, x, y);
}
#endif
return iteration;
}
}
void *thread_launcher(void *arguments)
{
piece_args *args;
args = (piece_args *) arguments;
int x,y, small_res_x, small_res_y, init_x, init_y, limit_x, limit_y;
int iteration, split, piece_x, piece_y;
if(args->total_threads != 1)
{
split = args->total_threads / 2;
}
else
{
split = 1;
}
if (args->thread_number > 0)
{
piece_x = args->thread_number % split;
piece_y = floor((float)args->thread_number / (float)split);
}
else
{
piece_x = 0;
piece_y = 0;
}
small_res_x = floor((float)args->res_x / (float)split);
small_res_y = floor((float)args->res_y / (float)split);
init_x = small_res_x * piece_x;
init_y = small_res_y * piece_y;
limit_x = init_x + small_res_x;
limit_y = init_y + small_res_y;
for (y = init_y; y < limit_y; y++)
{
for (x = init_x; x < limit_x; x++)
{
#ifndef JULIA
iteration_pixels[x + (y * args->res_x)] = mandelbrot_point(args->res_x, args->res_y, x, y, args->zoom, args->max_iteration);
#else
iteration_pixels[x + (y * args->res_x)] = julia_point(args->res_x, args->res_y, x, y, args->zoom, args->max_iteration);
#endif
}
}
}
int get_cpus()
{
int number_of_cores = 0;
number_of_cores = sysconf(_SC_NPROCESSORS_ONLN);
return number_of_cores;
}
int main(int argn, char **argv)
{
// Init SDL
if(SDL_Init(SDL_INIT_VIDEO) != 0)
fprintf(stderr, "Could not initialize SDL: %s\n", SDL_GetError());
printf("SDL Initialized\n");
// Create screen surface
SDL_Surface *screen;
int res_x = 800;
int res_y = 600;
int number_threads = get_cpus();
printf("Number of threads autodetect: %d\n", number_threads);
#ifdef CACHE
// Init our cached points
cached_points = malloc(res_y * 1000 * sizeof(int *));
cached_x = malloc(res_y * 1000 * sizeof(float *));
cached_y = malloc(res_y * 1000 * sizeof(float *));
if (cached_points == NULL)
{
fprintf(stderr, "Bad luck, out of memory\n");
return 2;
}
int count;
for (count = 0; count < res_y * 1000; count++)
{
cached_points[count] = malloc(res_x * 1000 * sizeof(int));
if(cached_points[count] == NULL)
{
fprintf(stderr, "Bad luck, out of memory\n");
return 2;
}
cached_x[count] = malloc(res_x * 1000 * sizeof(float));
cached_y[count] = malloc(res_x * 1000 * sizeof(float));
/*for (count2 = 0; count2 < res_x * 100; count2++)
{
cached_points[count][count2] = -1;
}*/
}
printf("Cache ready\n");
#endif
// screen = SDL_SetVideoMode(res_x, res_y, 0, SDL_HWSURFACE|SDL_DOUBLEBUF);
screen = SDL_SetVideoMode(res_x, res_y, 0, SDL_DOUBLEBUF);
if(!screen)
fprintf(stderr,"Could not set video mode: %s\n",SDL_GetError());
// Prepare the resolution and sizes and colors, threads...
int i;
iteration_pixels = malloc(res_x * res_y * sizeof(int));
pthread_t threads[number_threads];
piece_args arguments[number_threads];
printf("Rendering...\n");
float zoom;
#ifndef JULIA
for (zoom = 1.0; zoom > 0.0001 ; zoom = zoom * 0.98)
#else
for (zoom = 1.0; zoom > -2.5 ; zoom -= 0.01)
#endif
{
i = 0;
int iteration, max_iteration, x, y, res;
if((zoom < -0.02) && (zoom > -1.0))
{
max_iteration = 100;
}
else
{
max_iteration = 170;
}
int thread_count;
for(thread_count = 0; thread_count < number_threads; thread_count++)
{
arguments[thread_count].res_x = res_x;
arguments[thread_count].res_y = res_y;
arguments[thread_count].zoom = zoom;
arguments[thread_count].max_iteration = max_iteration;
arguments[thread_count].total_threads = number_threads;
arguments[thread_count].thread_number = thread_count;
pthread_create( &threads[thread_count], NULL, thread_launcher, (void*) &arguments[thread_count]);
}
for(thread_count = 0; thread_count < number_threads; thread_count++)
{
res = pthread_join(threads[thread_count], NULL);
if (res != 0)
{
printf("Error in %d thread\n", thread_count);
}
}
int rank;
Uint32 *pixel;
rank = screen->pitch/sizeof(Uint32);
pixel = (Uint32*)screen->pixels;
for(y = 0; y < res_y ; y++)
{
for(x = 0; x < res_x; x++)
{
iteration = iteration_pixels[x + y * res_x];
if ((iteration < 128) && (iteration > 0)) {
pixel[x + y * rank] = SDL_MapRGBA(screen->format,
0,
20 + iteration,
0,
255);
}
else if ((iteration >= 128) && (iteration < max_iteration))
{
pixel[x + y * rank] = SDL_MapRGBA(screen->format,
iteration,
148,
iteration,
255);
}
else
{
pixel[x + y * rank] = SDL_MapRGBA(screen->format,
0,
0,
0,
255);
}
}
}
SDL_Flip(screen);
}
// printf("Max Iteration value: %d\n", max_iter);
SDL_Event ev;
int active;
active = 1;
while(active)
{
/* Handle events */
while(SDL_PollEvent(&ev))
{
if(ev.type == SDL_QUIT)
active = 0; /* End */
}
}
SDL_Quit();
return 0;
}

BIN
mandel_classic.o Normal file
View file

Binary file not shown.

45
mandelbrot_kernel.cl Normal file
View file

@ -0,0 +1,45 @@
float map_x(int x, int width)
{
return (((float)x / (float)width) * 3.5) - 2.5;
}
float map_y(int y, int height)
{
return (((float)y / (float)height) * 2.0) - 1.0;
}
__kernel void mandelbrot_point(__global const int res_x, __global const int res_y, __global const int line, __global int *graph_line)
{
// Get the index of the current element
int image_x = get_global_id(0);
int image_y = line;
float pos_x = map_x(image_x, res_x);
float pos_y = map_y(image_y, res_y);
float x = 0.0;
float y = 0.0;
int iteration = 0;
int max_iteration = 255;
float xtemp;
while (iteration < max_iteration)
{
xtemp = x * x - y * y + pos_x;
y = 2.0 * x * y + pos_y;
x = xtemp;
iteration++;
if ((x * x) + (y * y) >= (4.0)) break;
}
if (iteration >= max_iteration)
{
graph_line[image_x] = 0;
}
else
{
graph_line[image_x] = iteration;
}
}

BIN
mandelclassic Executable file
View file

Binary file not shown.

BIN
test Executable file
View file

Binary file not shown.

301
test.c Normal file
View file

@ -0,0 +1,301 @@
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <SDL.h>
#define MAX_SOURCE_SIZE (0x100000)
#ifdef CACHE
int** cached_points;
int** cached_x;
int** cached_y;
#endif
int get_x (int linear_point, int width)
{
return linear_point % width;
}
int get_y (int linear_point, int height)
{
return floor(linear_point / height);
}
double map_x(int x, int width, double zoom)
{
#ifndef JULIA
return (((double)x / (double)width) * (3.5 * zoom)) - (2.5 - (1.0 - zoom));
#else
return (((double)x / (double)width) * (3.5 * zoom)) - (1.7 - (1.0 - zoom));
#endif
}
double map_y(int y, int height, double zoom)
{
return (((double)y / (double)height) * (2.0 * zoom)) - (1.0 - (1.0 - zoom));
}
#ifdef CACHE
int cached_iteration(double pos_x, double pos_y)
{
double centered_x = pos_x + 2.5;
double centered_y = pos_y + 1.0;
double temp_x = floor(centered_x * 1000.0);
double temp_y = floor(centered_y * 1000.0);
int trs_pos_x = (int)temp_x;
int trs_pos_y = (int)temp_y;
return cached_points[trs_pos_x][trs_pos_y];
}
double get_cached_x(double pos_x, double pos_y)
{
double centered_x = pos_x + 2.5;
double centered_y = pos_y + 1.0;
double temp_x = floor(centered_x * 1000.0);
double temp_y = floor(centered_y * 1000.0);
int trs_pos_x = (int)temp_x;
int trs_pos_y = (int)temp_y;
return cached_x[trs_pos_x][trs_pos_y];
}
double get_cached_y(double pos_x, double pos_y)
{
double centered_x = pos_x + 2.5;
double centered_y = pos_y + 1.0;
double temp_x = floor(centered_x * 1000.0);
double temp_y = floor(centered_y * 1000.0);
int trs_pos_x = (int)temp_x;
int trs_pos_y = (int)temp_y;
return cached_y[trs_pos_x][trs_pos_y];
}
void store_iteration(double pos_x, double pos_y, int iteration, double x, double y)
{
double centered_x = pos_x + 2.5;
double centered_y = pos_y + 1.0;
double temp_x = floor(centered_x * 1000.0);
double temp_y = floor(centered_y * 1000.0);
int trs_pos_x = (int)temp_x;
int trs_pos_y = (int)temp_y;
cached_points[trs_pos_x][trs_pos_y] = iteration;
cached_x[trs_pos_x][trs_pos_y] = x;
cached_y[trs_pos_x][trs_pos_y] = y;
}
#endif
int mandelbrot_point(int res_x, int res_y, int image_x, int image_y, double zoom, int max_iteration)
{
return abs(floor(sin(((double)image_x + zoom) * 0.1) * 127) + floor(cos(((double)image_y + zoom) * 0.1) * 127));
}
int julia_point(int res_x, int res_y, int image_x, int image_y, double zoom, int max_iteration)
{
// Get the index of the current element
double pos_x = map_x(image_x, res_x, 1.0);
double pos_y = map_y(image_y, res_y, 1.0);
double x = pos_x;
double y = pos_y;
double q, x_term;
double xtemp, xx, yy;
#ifdef CACHE
int storeable = 1;
#endif
int iteration = 0;
#ifdef CACHE
// Look up our cache
iteration = cached_iteration(pos_x, pos_y);
if (iteration > 0)
{
x = get_cached_x(pos_x, pos_y);
y = get_cached_y(pos_x, pos_y);
yy = y * y;
}
if (iteration < 0) storeable = 0;
#endif
while (iteration < max_iteration)
{
xx = x * x;
yy = y * y;
if ((xx) + (yy) > (4.0)) break;
y = pow((x + y), 2) - xx - yy;
y = y + 0.288;
xtemp = xx - yy + 0.353 + zoom;
x = xtemp;
iteration++;
}
if (iteration >= max_iteration)
{
return 0;
}
else
{
#ifdef CACHE
if (storeable == 1)
{
store_iteration(pos_x, pos_y, iteration, x, y);
}
#endif
return iteration;
}
}
int main(int argn, char **argv) {
// Init SDL
if(SDL_Init(SDL_INIT_VIDEO) != 0)
fprintf(stderr, "Could not initialize SDL: %s\n", SDL_GetError());
printf("SDL Initialized\n");
// Create screen surface
SDL_Surface *screen;
int res_x = 800;
int res_y = 600;
int total_res = res_x * res_y;
#ifdef CACHE
// Init our cached points
cached_points = malloc(res_y * 1000 * sizeof(int *));
cached_x = malloc(res_y * 1000 * sizeof(double *));
cached_y = malloc(res_y * 1000 * sizeof(double *));
if (cached_points == NULL)
{
fprintf(stderr, "Bad luck, out of memory\n");
return 2;
}
int count;
for (count = 0; count < res_y * 1000; count++)
{
cached_points[count] = malloc(res_x * 1000 * sizeof(int));
if(cached_points[count] == NULL)
{
fprintf(stderr, "Bad luck, out of memory\n");
return 2;
}
cached_x[count] = malloc(res_x * 1000 * sizeof(double));
cached_y[count] = malloc(res_x * 1000 * sizeof(double));
/*for (count2 = 0; count2 < res_x * 100; count2++)
{
cached_points[count][count2] = -1;
}*/
}
printf("Cache ready\n");
#endif
// screen = SDL_SetVideoMode(res_x, res_y, 0, SDL_HWSURFACE|SDL_DOUBLEBUF);
screen = SDL_SetVideoMode(res_x, res_y, 0, SDL_DOUBLEBUF);
if(!screen)
fprintf(stderr,"Could not set video mode: %s\n",SDL_GetError());
// Prepare the resolution and sizes and colors...
int i;
int temp;
const int ITERATIONS = 256;
printf("Rendering...\n");
double zoom;
#ifndef JULIA
for (zoom = 1.0; zoom < 200.0 ; zoom += 1.0)
#else
for (zoom = 1.0; zoom > -2.5 ; zoom -= 0.01)
#endif
{
i = 0;
int iteration, max_iteration, x, y;
if((zoom < -0.02) && (zoom > -1.0))
{
max_iteration = 100;
}
else
{
max_iteration = 255;
}
int col_value;
Uint32 *pixel;
int rank;
// Lock surface
// SDL_LockSurface(screen);
rank = screen->pitch/sizeof(Uint32);
pixel = (Uint32*)screen->pixels;
/* Draw all dots */
for(y = 0;y < res_y;y++)
{
for(x = 0;x < res_x;x++)
{
#ifndef JULIA
iteration = mandelbrot_point(res_x, res_y, x, y, zoom, max_iteration);
#else
iteration = julia_point(res_x, res_y, x, y, zoom, max_iteration);
#endif
if ((iteration < 128) && (iteration > 0)) {
pixel[x + y * rank] = SDL_MapRGBA(screen->format,
0,
20 + iteration,
0,
255);
}
else if ((iteration >= 128) && (iteration < max_iteration))
{
pixel[x + y * rank] = SDL_MapRGBA(screen->format,
iteration,
200,
iteration,
255);
}
else
{
pixel[x + y * rank] = SDL_MapRGBA(screen->format,
0,
0,
0,
255);
}
i++;
}
}
// Unlock surface
// SDL_UnlockSurface(screen);
// Draw to the screen
SDL_Flip(screen);
}
// printf("Max Iteration value: %d\n", max_iter);
SDL_Event ev;
int active;
active = 1;
while(active)
{
/* Handle events */
while(SDL_PollEvent(&ev))
{
if(ev.type == SDL_QUIT)
active = 0; /* End */
}
}
SDL_Quit();
return 0;
}

BIN
test.o Normal file
View file

Binary file not shown.

12
vector_add_kernel.cl Normal file
View file

@ -0,0 +1,12 @@
int vec_add (int a, int b) {
return (a + b) * 3;
}
__kernel void vector_add(__global int *A, __global int *B, __global int *C) {
// Get the index of the current element
int i = get_global_id(0);
// Do the operation
C[i] = vec_add (A[i], B[i]);
}