# include # include # include # include # include "cuda.h" int main ( ); void mandelbrot_image ( float xmin, float xmax, float ymin, float ymax, int xnum, int ynum, int count_max ); __global__ void mandelbrot_set ( float xmin, float xmax, float ymin, float ymax, int xnum, int ynum, int count_max, int count[] ); __device__ int mandelbrot_count ( float cx, float cy, int count_max ); void pgma_write ( char *file_name, char *comment, int xsize, int ysize, int maxval, int *gray ); /******************************************************************************/ int main ( ) /******************************************************************************/ /* Purpose: MANDELBROT_CUDA creates an image of the Mandelbrot set. Discussion: Over the rectangle [XMIN,XMAX] x [YMIN,YMAX], determine the Mandelbrot count for a grid of XNUMxYNUM points, using a particular value of COUNT_MAX. Licensing: This code is distributed under the GNU LGPL license. Modified: 15 May 2017 Author: John Burkardt Local Parameters: Local, float XMIN, XMAX, YMIN, YMAX, the physical limits of the rectangle. Local, int XNUM, YNUM, the number of points in X and Y directions. Local, int COUNT_MAX, the maximum number of iterations. */ { int count_max; float xmax; float xmin; int xnum; float ymax; float ymin; int ynum; printf ( "\n" ); printf ( "MANDELBROT_CUDA:\n" ); printf ( " CUDA version\n" ); printf ( " Compute the Mandelbrot set and save a graphic image of it.\n" ); xmin = -1.0; xmax = -0.6; ymin = 0.0; ymax = 0.4; xnum = 1000; ynum = 1000; count_max = 200; printf ( "\n" ); printf ( " X range: [ %g, %g ]\n", xmin, xmax ); printf ( " Y range: [ %g, %g ]\n", ymin, ymax ); printf ( " Xnum = %d x Ynum = %d = %d Pixels\n", xnum, ynum, xnum * ynum ); printf ( " Maximum number of iterations = %d\n", count_max ); mandelbrot_image ( xmin, xmax, ymin, ymax, xnum, ynum, count_max ); /* Terminate. */ printf ( "\n" ); printf ( "MANDELBROT_CUDA:\n" ); printf ( " Normal end of execution.\n" ); return 0; } /******************************************************************************/ void mandelbrot_image ( float xmin, float xmax, float ymin, float ymax, int xnum, int ynum, int count_max ) /******************************************************************************/ /* Purpose: MANDELBROT_IMAGE creates an image of the Mandelbrot set. Discussion: Over the rectangle [XMIN,XMAX] x [YMIN,YMAX], determine the Mandelbrot count for a grid of XNUMxYNUM points, using a particular value of COUNT_MAX. Licensing: This code is distributed under the GNU LGPL license. Modified: 15 May 2017 Author: John Burkardt Parameters: Input, float XMIN, XMAX, YMIN, YMAX, the physical limits of the rectangle. Input, int XNUM, YNUM, the number of points in X and Y directions. Input, int COUNT_MAX, the maximum number of iterations. */ { clock_t clock1; clock_t clock2; char comment[] = "PGM gray map created by mandelbrot_cuda.cu"; int *count_cpu; int *count_gpu; int count_size; char file_name[] = "mandelbrot_cuda.pgm"; struct timespec time1; struct timespec time2; double value; count_size = xnum * ynum * sizeof ( int ); /* Allocate space for CPU version of COUNT. */ count_cpu = ( int * ) malloc ( count_size ); /* Allocate space for GPU version of COUNT. */ cudaMalloc ( ( void ** ) &count_gpu, count_size ); /* Set timers and call kernel. */ clock_gettime ( CLOCK_REALTIME, &time1 ); clock1 = clock ( ); dim3 blocks ( xnum, ynum, 1 ); dim3 threads ( 1, 1, 1 ); mandelbrot_set <<< blocks, threads >>> ( xmin, xmax, ymin, ymax, xnum, ynum, count_max, count_gpu ); cudaMemcpy ( count_cpu, count_gpu, count_size, cudaMemcpyDeviceToHost ); clock2 = clock ( ); clock_gettime ( CLOCK_REALTIME, &time2 ); /* Report times. */ printf ( "\n" ); printf ( " Compute time:\n" ); value = ( double ) ( ( time2.tv_sec + time2.tv_nsec * 1.0E-09 ) - ( time1.tv_sec + time1.tv_nsec * 1.0E-09 ) ); printf ( " Wallclock: %.02f sec.\n", value ); value = ( double ) ( clock2 - clock1 ) / ( double ) CLOCKS_PER_SEC; printf ( " CPU: %.02f sec.\n", value ); /* Write data to PGM file. */ pgma_write ( file_name, comment, xnum, ynum, count_max, count_cpu ); printf ( "\n" ); printf ( " Plot saved in file '%s'\n", file_name ); /* Free memory. */ cudaFree ( count_gpu ); free ( count_cpu ); return; } /******************************************************************************/ __global__ void mandelbrot_set ( float xmin, float xmax, float ymin, float ymax, int xnum, int ynum, int count_max, int count[] ) /******************************************************************************/ /* Purpose: MANDELBROT_SET computes the Mandelbrot count for a grid of points. Discussion: Over the rectangle [XMIN,XMAX] x [YMIN,YMAX], determine the Mandelbrot count for a grid of XNUMxYNUM points, using a particular value of COUNT_MAX. Licensing: This code is distributed under the GNU LGPL license. Modified: 12 May 2017 Author: John Burkardt Parameters: Input, float XMIN, XMAX, YMIN, YMAX, the physical limits of the rectangle. Input, int XNUM, YNUM, the number of points in X and Y directions. Input, int COUNT_MAX, the maximum number of iterations. Output, int COUNT[XNUM*YNUM], the count for each point in the grid. */ { float cx; float cy; int i; int j; i = blockIdx.x; j = blockIdx.y; cx = ( ( float ) ( xnum - 1 - i ) * xmin + ( float ) ( i ) * xmax ) / ( float ) ( xnum - 1 ); cy = ( ( float ) ( ynum - 1 - j ) * ymax + ( float ) ( j ) * ymin ) / ( float ) ( ynum - 1 ); count[i+j*xnum] = mandelbrot_count ( cx, cy, count_max ); return; } /******************************************************************************/ __device__ int mandelbrot_count ( float cx, float cy, int count_max ) /******************************************************************************/ /* Purpose: MANDELBROT_COUNT returns the Mandelbrot count for a single point. Discussion: Starting with the value 0, repeatedly square and add complex value C. Return number of applications of this process at which the value exceeds 2 in norm. Repeat no more than COUNT_MAX times. CUDA's support of complex numbers and calling functions inside a kernel is bizarre enough that it's better to just use real and imaginary parts. Licensing: This code is distributed under the GNU LGPL license. Modified: 15 May 2017 Author: John Burkardt Parameters: Input, float CX, CY, the real and imaginary parts of the value added at each step. Input, int COUNT_MAX, the maximum number of iterations. Output, int MANDELBROT_COUNT, the number of operations at which the iterate first exceeded 2 in norm. If this never happens, return COUNT_MAX. */ { int k; int value; float zx; float zxnew; float zy; float zynew; value = count_max; zx = 0.0; zy = 0.0; for ( k = 0; k < count_max; k++ ) { if ( 4.0 <= ( zx * zx + zy * zy ) ) { value = k; break; } zxnew = zx * zx - zy * zy + cx; zynew = 2.0 * zx * zy + cy; zx = zxnew; zy = zynew; } return value; } /******************************************************************************/ void pgma_write ( char *file_name, char *comment, int xsize, int ysize, int maxval, int *gray ) /******************************************************************************/ /* Purpose: PGMA_WRITE writes the header and data for an ASCII PGM file. Example: P2 # feep.pgm 24 7 15 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 3 3 3 0 0 7 7 7 7 0 0 11 11 11 11 0 0 15 15 15 15 0 0 3 0 0 0 0 0 7 0 0 0 0 0 11 0 0 0 0 0 15 0 0 15 0 0 3 3 3 0 0 0 7 7 7 0 0 0 11 11 11 0 0 0 15 15 15 15 0 0 3 0 0 0 0 0 7 0 0 0 0 0 11 0 0 0 0 0 15 0 0 0 0 0 3 0 0 0 0 0 7 7 7 7 0 0 11 11 11 11 0 0 15 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Licensing: This code is distributed under the GNU LGPL license. Modified: 13 May 2017 Author: John Burkardt Parameters: Input, char *FILE_NAME, the name of the file. Input, char *COMMENT, a comment. Input, int XSIZE, YSIZE, the number of rows and columns of data. Input, int MAXVAL, the maximum value allowed. Input, int *GRAY, the array of XSIZE by YSIZE data values. */ { FILE *file_handle; int i; int *indexg; int j; /* Open the output file. */ file_handle = fopen ( file_name, "wt" ); if ( ! file_handle ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PGMA_WRITE - Fatal error!\n" ); fprintf ( stderr, " Cannot open the output file \"%s\".\n", file_name ); exit ( 1 ); } /* Write the header. */ fprintf ( file_handle, "P2\n" ); fprintf ( file_handle, "#%s\n", comment ); fprintf ( file_handle, "%d %d\n", xsize, ysize ); fprintf ( file_handle, "%d\n", maxval ); /* Write the data. */ indexg = gray; for ( j = 0; j < ysize; j++ ) { for ( i = 0; i < xsize; i++ ) { fprintf ( file_handle, " %d", *indexg ); indexg = indexg + 1; } fprintf ( file_handle, "\n" ); } /* Close the file. */ fclose ( file_handle ); return; }