diff options
| -rw-r--r-- | makefile | 2 | ||||
| -rw-r--r-- | src/cuda-fractals.cu | 55 | ||||
| -rw-r--r-- | src/fractals.c | 96 | ||||
| -rw-r--r-- | src/grids.c | 23 | ||||
| -rw-r--r-- | src/grids.h | 3 |
5 files changed, 122 insertions, 57 deletions
@@ -67,7 +67,7 @@ presentation/presentation.html: presentation/presentation.md analysis: analysis/analysis.html -analysis/analysis.html: analysis/analysis.Rmd # TODO: add compile command +analysis/analysis.html: analysis/analysis.Rmd analysis/data.csv Rscript -e "setwd('analysis'); rmarkdown::render('analysis.Rmd')" clean: diff --git a/src/cuda-fractals.cu b/src/cuda-fractals.cu index ef94a3f..52e7428 100644 --- a/src/cuda-fractals.cu +++ b/src/cuda-fractals.cu @@ -1,10 +1,10 @@ +/* + * Compute fractals using CUDA + */ #include <cuda_runtime.h> #include <thrust/complex.h> #include "fractals.h" - -//These won't work/ need to be modifiied for the cuda version -//#include "grids.h" -//#include "fractals.h" +#include "grids.h" /* * Macro for checking CUDA errors @@ -22,8 +22,13 @@ #define CBASE double #endif +#ifndef BLOCK_SIZE_X #define BLOCK_SIZE_X 16 +#endif + +#ifndef BLOCK_SIZE_Y #define BLOCK_SIZE_Y 16 +#endif // this ain't pretty, but this logic should never change for kernels #define SET_ROW_COL \ @@ -32,6 +37,9 @@ if(row >= rows || col >= cols) return +/* + * Device helper function to cconvert a grid point into its corresponding complex number + */ __device__ thrust::complex<CBASE> grid_to_complex(const thrust::complex<CBASE> lower_left, const thrust::complex<CBASE> upper_right, const size_t row, const size_t col, const size_t rows, const size_t cols){ const CBASE x_min = lower_left.real(); @@ -48,6 +56,9 @@ thrust::complex<CBASE> grid_to_complex(const thrust::complex<CBASE> lower_left, return z; } +/* + * Device function to compute if a point is or is not in the mandelbrot set + */ __device__ byte mandelbrot(const thrust::complex<CBASE> z0, const byte max_iterations){ thrust::complex<CBASE> z = z0; @@ -59,6 +70,9 @@ byte mandelbrot(const thrust::complex<CBASE> z0, const byte max_iterations){ return iteration; } +/* + * Kernel to run mandelbrot on an entire grid + */ __global__ void mandelbrot_kernel(byte* grid_data, const byte max_iterations, const thrust::complex<CBASE> lower_left, const thrust::complex<CBASE> upper_right, const size_t rows, const size_t cols){ SET_ROW_COL; @@ -68,6 +82,9 @@ void mandelbrot_kernel(byte* grid_data, const byte max_iterations, const thrust: grid_data[row*cols + col] = mandelbrot(z, max_iterations); } +/* + * Device function to compute if a point is or is not in the tricorn set + */ __device__ byte tricorn(const thrust::complex<CBASE> z0, const byte max_iterations){ thrust::complex<CBASE> z = z0; @@ -79,6 +96,9 @@ byte tricorn(const thrust::complex<CBASE> z0, const byte max_iterations){ return iteration; } +/* + * Kernel to run tricorn on an entire grid + */ __global__ void tricorn_kernel(byte* grid_data, const byte max_iterations, const thrust::complex<CBASE> lower_left, const thrust::complex<CBASE> upper_right, const size_t rows, const size_t cols){ SET_ROW_COL; @@ -88,6 +108,9 @@ void tricorn_kernel(byte* grid_data, const byte max_iterations, const thrust::co grid_data[row*cols + col] = tricorn(z, max_iterations); } +/* + * Device function to compute if a point is or is not in the burning ship set + */ __device__ byte burning_ship(const thrust::complex<CBASE> z0, const byte max_iterations){ thrust::complex<CBASE> z = z0; @@ -101,6 +124,9 @@ byte burning_ship(const thrust::complex<CBASE> z0, const byte max_iterations){ return iteration; } +/* + * Kernel to run burning_ship on an entire grid + */ __global__ void burning_ship_kernel(byte* grid_data, const byte max_iterations, const thrust::complex<CBASE> lower_left, const thrust::complex<CBASE> upper_right, const size_t rows, const size_t cols){ SET_ROW_COL; @@ -110,6 +136,9 @@ void burning_ship_kernel(byte* grid_data, const byte max_iterations, const thrus grid_data[row*cols + col] = burning_ship(z, max_iterations); } +/* + * Device function to compute if a point is or is not in the multibrot set + */ __device__ byte multibrot(const thrust::complex<CBASE> z0, const byte max_iterations, const double d){ thrust::complex<CBASE> z = z0; @@ -121,6 +150,9 @@ byte multibrot(const thrust::complex<CBASE> z0, const byte max_iterations, const return iteration; } +/* + * Kernel to run multibrot on an entire grid + */ __global__ void multibrot_kernel(byte* grid_data, const double degree, const byte max_iterations, const thrust::complex<CBASE> lower_left, const thrust::complex<CBASE> upper_right, const size_t rows, const size_t cols){ SET_ROW_COL; @@ -130,6 +162,9 @@ void multibrot_kernel(byte* grid_data, const double degree, const byte max_itera grid_data[row*cols + col] = multibrot(z, max_iterations, degree); } +/* + * Device function to compute if a point is or is not in the multicorn set + */ __device__ byte multicorn(const thrust::complex<CBASE> z0, const byte max_iterations, const double d){ thrust::complex<CBASE> z = z0; @@ -141,6 +176,9 @@ byte multicorn(const thrust::complex<CBASE> z0, const byte max_iterations, const return iteration; } +/* + * Kernel to run multicorn on an entire grid + */ __global__ void multicorn_kernel(byte* grid_data, const double degree, const byte max_iterations, const thrust::complex<CBASE> lower_left, const thrust::complex<CBASE> upper_right, const size_t rows, const size_t cols){ SET_ROW_COL; @@ -150,6 +188,9 @@ void multicorn_kernel(byte* grid_data, const double degree, const byte max_itera grid_data[row*cols + col] = multicorn(z, max_iterations, degree); } +/* + * Device function to compute if a point is or is not in the julia set + */ __device__ byte julia(const thrust::complex<CBASE> z0, const byte max_iterations, const thrust::complex<CBASE> c, const double R){ thrust::complex<CBASE> z = z0; @@ -161,6 +202,9 @@ byte julia(const thrust::complex<CBASE> z0, const byte max_iterations, const thr return iteration; } +/* + * Kernel to run julia on an entire grid + */ __global__ void julia_kernel(byte* grid_data, const thrust::complex<CBASE> constant, const double radius, const byte max_iterations, const thrust::complex<CBASE> lower_left, const thrust::complex<CBASE> upper_right, const size_t rows, const size_t cols){ SET_ROW_COL; @@ -170,6 +214,7 @@ void julia_kernel(byte* grid_data, const thrust::complex<CBASE> constant, const grid_data[row*cols + col] = julia(z, max_iterations, constant, radius); } +// prevent c++ name mangling so the grid functions can be linked against in standard c code extern "C" { void mandelbrot_grid(grid_t* grid, const grid_gen_params* params){ const size_t size = grid->size; @@ -181,9 +226,7 @@ void mandelbrot_grid(grid_t* grid, const grid_gen_params* params){ byte* d_grid_data; CHECK(cudaMalloc(&d_grid_data, size*sizeof(byte))); - //TODO: find good sizes dim3 block_size(BLOCK_SIZE_X, BLOCK_SIZE_Y); - //dim3 grid_size(0,0); dim3 grid_size((cols + block_size.x - 1) / block_size.x, (rows + block_size.y - 1) / block_size.y); mandelbrot_kernel<<<grid_size, block_size>>>(d_grid_data, max_iterations, lower_left, upper_right, rows, cols); CHECK(cudaDeviceSynchronize()); diff --git a/src/fractals.c b/src/fractals.c index 898f30d..1fcce3d 100644 --- a/src/fractals.c +++ b/src/fractals.c @@ -16,10 +16,16 @@ #define NUM_RUNS 5 #endif +/* + * Prints out usage information for the program + */ void print_usage(FILE* file, const char* program_name){ fprintf(file, "Usage: %s [-v] [-i iterations] [-x x_res] [-y y_res] [-z magnification] [-d degree] [-c constant] [-r radius] [-l lower_left] [-u upper_right] [-o output_grid] -f fractal\n", program_name); } +/* + * Print out the full help message for the program + */ void print_help(){ printf("Options:\n" " -i, --iterations <value> the number of iterations (default: 25, max 255)\n" @@ -41,6 +47,9 @@ void print_help(){ "\nExits with a status code of 1 if the program encounters an error, exits with 2 if an argument is incorrect\n"); } +/* + * Prints out additional information about the program, such as the precision it was compiled with + */ void print_info(const char* program_name){ #ifdef EXTENDED_PRECISION printf("Compiled with long double float precision\n"); @@ -50,6 +59,9 @@ void print_info(const char* program_name){ #endif } +/* + * Runs a fractal generator NUM_RUNS times and returns the average of those runs + */ double time_fractal(fractal_generator generator, grid_t* grid, grid_gen_params* params){ struct timespec start, end; clock_gettime(CLOCK_MONOTONIC, &start); @@ -67,6 +79,47 @@ static inline void parse_complex(const char* string, complex_t* z){ }; } +/* + * Parses a fractal generator form a string, exitting if it is supplied a parameter which it doesn't support + */ +fractal_generator parse_fractal_generator(const char* argument, const bool param_is_degree, const bool param_is_cr){ + if(strncmp(argument, "mandelbrot", strlen("mandelbrot")) == 0) { + return mandelbrot_grid; + } + else if(strncmp(argument, "tricorn", strlen("tricorn")) == 0) { + return tricorn_grid; + } + else if(strncmp(argument, "multibrot", strlen("multibrot")) == 0) { + if(param_is_cr){ + fprintf(stderr, "multibrot requires a degree, not constant and radius, exitting\n"); + exit(EXIT_BAD_ARGUMENT); + } + return multibrot_grid; + } + else if(strncmp(argument, "multicorn", strlen("multicorn")) == 0) { + if(param_is_cr){ + fprintf(stderr, "multicorn requires a degree, not constant and radius, exitting\n"); + exit(EXIT_BAD_ARGUMENT); + } + return multicorn_grid; + } + else if(strncmp(argument, "burning_ship", strlen("burning_ship")) == 0) { + return burning_ship_grid; + } + else if(strncmp(argument, "julia", strlen("julia")) == 0) { + if(param_is_degree){ + fprintf(stderr, "julia requires a constant and a radius, not a degree, exitting\n"); + exit(EXIT_BAD_ARGUMENT); + } + return julia_grid; + } + else { + fprintf(stderr, "Invalid fractal type: %s, see --help for a list of supported fractals\n", argument); + exit(EXIT_BAD_ARGUMENT); + } +} + + int main(const int argc, char *argv[]) { struct winsize w; ioctl(STDOUT_FILENO, TIOCGWINSZ, &w); @@ -82,6 +135,7 @@ int main(const int argc, char *argv[]) { bool performance = false; //degree is mutually exclusive with constant and radius + //this could be simplified if grid_gen_params was a struct instead of a union but ¯\_(ツ)_/¯ bool param_is_degree = false; bool param_is_cr = false; CBASE degree = 2; @@ -177,46 +231,8 @@ int main(const int argc, char *argv[]) { param_is_cr = true; break; case 'f': - if(strncmp(optarg, "mandelbrot", strlen("mandelbrot")) == 0) { - fractal_name = "mandelbrot"; - generator = mandelbrot_grid; - } - else if(strncmp(optarg, "tricorn", strlen("tricorn")) == 0) { - fractal_name = "tricorn"; - generator = tricorn_grid; - } - else if(strncmp(optarg, "multibrot", strlen("multibrot")) == 0) { - if(param_is_cr){ - fprintf(stderr, "multibrot requires a degree, not constant and radius, exitting\n"); - exit(EXIT_BAD_ARGUMENT); - } - fractal_name = "multibrot"; - generator = multibrot_grid; - } - else if(strncmp(optarg, "multicorn", strlen("multicorn")) == 0) { - if(param_is_cr){ - fprintf(stderr, "multicorn requires a degree, not constant and radius, exitting\n"); - exit(EXIT_BAD_ARGUMENT); - } - fractal_name = "multicorn"; - generator = multicorn_grid; - } - else if(strncmp(optarg, "burning_ship", strlen("burning_ship")) == 0) { - fractal_name = "burning ship"; - generator = burning_ship_grid; - } - else if(strncmp(optarg, "julia", strlen("julia")) == 0) { - if(param_is_degree){ - fprintf(stderr, "julia requires a constant and a radius, not a degree, exitting\n"); - exit(EXIT_BAD_ARGUMENT); - } - fractal_name = "julia"; - generator = julia_grid; - } - else { - fprintf(stderr, "Invalid fractal type: %s, see --help for a list of supported fractals\n", optarg); - exit(EXIT_BAD_ARGUMENT); - } + fractal_name = optarg; + generator = parse_fractal_generator(optarg, param_is_degree, param_is_cr); break; case 'z': if(sscanf(optarg, CFORMAT, &magnification) != 1){ diff --git a/src/grids.c b/src/grids.c index e65ba79..279fdb5 100644 --- a/src/grids.c +++ b/src/grids.c @@ -1,3 +1,6 @@ +/* + * Functions for the creation, manipulation, and evaluation of grids. + */ #include <complex.h> #include <stdio.h> #include <stdlib.h> @@ -92,6 +95,8 @@ bool grid_equal(const grid_t* grid1_p, const grid_t* grid2_p){ /* * Checks if two grids have a given maximum difference + * This function has no few if any uses, experimentall the grids are either exact or very different, + * indicating an error in their generation */ bool grid_allclose(const grid_t* restrict grid1, const grid_t* restrict grid2, const byte max_error){ if(grid1->x != grid2->x || grid1->y != grid2->y || @@ -101,8 +106,6 @@ bool grid_allclose(const grid_t* restrict grid1, const grid_t* restrict grid2, c } const size_t size = grid1->size; for(size_t i = 0; i < size; i++){ - //FIXME: figure out how to handle difference between two unsigned values - // possibly cast them to sized longs? if(abs(grid1->data[i] - grid2->data[i]) >= max_error) return false; } return true; @@ -172,12 +175,11 @@ void zoom_grid(grid_t* restrict grid, const CBASE magnification){ * * Returns 0 on success * - * FIXME: add dynamic bounding box resolution by storing sizeof points in the file * The .grid format is a binary file format - * The first 3 bytes of the file are a magic number defined in grids.h + * The first 3 bytes of the file are a magic number defined in grids.h (a goose !) * The next 16 bytes are the grid dimensions (x then y) - * The next 8 bytes is the max_iterations - * The next 8 bytes are the size of a grid point in bytes + * The next 1 byte is the max_iterations + * The next 8 bytes are the precision for the complex numbers * The next 2*precision bytes are the lower left and upper right corners * The rest of the file is the data for the grid, which should be exactly x*y*8 bytes */ @@ -239,14 +241,13 @@ void print_grid(FILE* file, const grid_t* grid){ const byte iterations = grid->max_iterations; const byte* data = grid->data; - //TODO: set values in output buffer rather than multiple printf calls - // the buffer needs to be larger to hold newlines char* output_buffer = malloc(size + grid->y-1); if(!output_buffer){ fprintf(stderr, "Failed to allocate output buffer for %zu points\n", size); return; } + // sets an output buffer to minimize number of syscalls setvbuf(file, output_buffer, _IOFBF, size + grid->y - 1); const char point_types[] = { ' ', '.', '*', '%', '#'}; @@ -323,8 +324,10 @@ grid_t* read_grid(FILE* restrict file){ if(fread(&lower_left, sizeof(complex_t), 1, file) != 1){ longjmp(file_read_error, 1); } if(fread(&upper_right, sizeof(complex_t), 1, file) != 1){ longjmp(file_read_error, 1); } - //TODO: look into mmaping the file to data, offseting by the bounding and resolution information - // this would likely require an alloc_grid function, similar to jeff's implementation in hw03 + // NOTE: for very large sizes it might make sense to memmap the file + // this would work on the CPU but break the cuda implmentation + // a potential fix would to just make multiple grids, aand stich them in + // after image rendering grid_t* grid = create_grid(x, y, max_iterations, lower_left, upper_right); if(!grid){ return NULL; diff --git a/src/grids.h b/src/grids.h index 6bf96f4..129cfd5 100644 --- a/src/grids.h +++ b/src/grids.h @@ -20,6 +20,9 @@ typedef struct { CBASE im; } complex_t; +// I've gone back and forth on the size for a grid point, eventually I settled on a byte +// ideally it should just be its own type which can be configured similarly to precision.h +// but it does not seem worth typedef unsigned char byte; typedef struct { |