cleanup of .o and emacs temporary copies uploaded by mistake

# OpenMP Mandelbrot renderer

#

#   make            build mandelbrot_omp (optimised, -fopenmp)

#   make profile    build mandelbrot_omp.profile with -DPROFILE: it also writes

#                   <output>.dat, a per-pixel thread map for plot_thread_map.py

#   make debug      build with -O0 -g

#   make clean      remove objects and executables

#

# No external library is required: PNG output goes through the

# single-header stb_image_write.h, and libm is not needed.

CC      ?= gcc

CFLAGS  ?= -O3 -march=native -mtune=native -Wall -Wextra

OMP      = -fopenmp

LDFLAGS ?=

TARGET  = mandelbrot_omp

SRCS    = mandelbrot.omp.c mandelbrot_utilities.omp.c

OBJS    = $(SRCS:.c=.o)

HEADERS = mandelbrot_omp.h mandelbrot_loop.omp.h stb_image_write.h

all: $(TARGET)

$(TARGET): $(OBJS)

	$(CC) $(CFLAGS) $(OMP) -o $@ $(OBJS) $(LDFLAGS)

%.o: %.c $(HEADERS)

	$(CC) $(CFLAGS) $(OMP) -c -o $@ $<

profile: $(SRCS) $(HEADERS)

	$(CC) $(CFLAGS) $(OMP) -DPROFILE -o $(TARGET).profile $(SRCS) $(LDFLAGS)

debug: CFLAGS = -O0 -g -Wall -Wextra

debug: clean $(TARGET)

clean:

	rm -f $(TARGET) $(TARGET).profile $(OBJS)

.PHONY: all profile debug clean

/* ────────────────────────────────────────────────────────────────────────── *

│                                                                            │

│ This file is part of the exercises for the Lectures on                     │

│   "Foundations of High Performance Computing"                              │

│ given at                                                                   │

│   Master in HPC and                                                        │

│   Master in Data Science and Scientific Computing                          │

│ @ SISSA, ICTP and University of Trieste                                    │

│                                                                            │

│ contact: luca.tornatore@inaf.it                                            │

│                                                                            │

│     This is free software; you can redistribute it and/or modify           │

│     it under the terms of the GNU General Public License as published by   │

│     the Free Software Foundation; either version 3 of the License, or      │

│     (at your option) any later version.                                    │

│     This code is distributed in the hope that it will be useful,           │

│     but WITHOUT ANY WARRANTY; without even the implied warranty of         │

│     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          │

│     GNU General Public License for more details.                           │

│                                                                            │

│     You should have received a copy of the GNU General Public License      │

│     along with this program.  If not, see <http://www.gnu.org/licenses/>   │

│                                                                            │

* ────────────────────────────────────────────────────────────────────────── */

/*

  OpenMP Mandelbrot renderer for the HPC final-exam exercise.

  This is the OpenMP solution of the serial baseline in ../ : the pixel loop

  (mandelbrot_loop.h) is distributed among threads with a scheduling policy

  chosen at run time through --mode, and per-thread timings are collected to

  discuss load balancing.  Compile with -DPROFILE to also dump, next to the

  image, a map of which thread computed each pixel (<output>.dat), to be

  plotted with plot_thread_map.py.

  Command-line parsing and image writers (PPM or PNG, through the stb header)

  are in mandelbrot_utilities.c; shared declarations are in mandelbrot_omp.h.

  Have fun.

*/

// hint for timing

// you can use either the TIME macro that you find in the examples/timing.h header

// or

//   omp_get_wtime()

// that returns the wall-clock time per thread

//

#include "mandelbrot.omp.h"

/* ========================================================================================

   : ------------------------------------------------------ :

   :  CALCULATION FUNCTIONS                                 :

   : ------------------------------------------------------ :

 */ 

/*

  Compute the Mandelbrot escape iteration for one point c = cr + i ci.  This

  is the direct scalar kernel, clean, keeps the operation count easy to read.

  This can be heavility optimized.

*/

static unsigned int

mandelbrot_escape (double       cr,      // real part of the pixel centre

		   double       ci,      // imaginary part of the pixel centre

		   unsigned int kmax     // maximum number of iterations

		   )

{

           double zr;

           double zi;

           double zr_next;

           double zr2;

           double zi2;

  unsigned int    k;

  zr = 0.0;

  zi = 0.0;

  k = 0;

  while (k < kmax)

    {

      zr2 = zr * zr;

      zi2 = zi * zi;

      if (zr2 + zi2 > 4.0)

        break;

      zr_next = zr2 - zi2 + cr;

      zi = 2.0 * zr * zi + ci;

      zr = zr_next;

      k += 1;

    }

  return k;

}

/*

  Map one iteration count to an RGB colour.  Interior points are black.  The

  outside palette is deliberately simple.

  It might be fun experimenting by advanced colouring or smooth-potential formulas.

  Don't be distracted by this, though

*/

static void

colour_from_iteration (unsigned int    iteration,    // escape iteration count

		       unsigned int    kmax,         // maximum iteration count

		       unsigned char  *rgb           // output triplet: red, green, blue

		       )

{

  double t;

  double one_minus_t;

  if (iteration >= kmax)

    {

      rgb[0] = 0u;

      rgb[1] = 0u;

      rgb[2] = 0u;

      return;

    }

  t = (double) iteration / (double) kmax;

  one_minus_t = 1.0 - t;

  rgb[0] = (unsigned char) (255.0 * 9.0 * one_minus_t * t * t * t);

  rgb[1] = (unsigned char) (255.0 * 15.0 * one_minus_t * one_minus_t * t * t);

  rgb[2] = (unsigned char) (255.0 * 8.5 * one_minus_t * one_minus_t * one_minus_t * t);

}

static void

render_image (options_t        *options,     // view, factors, and kmax

	      image_geometry_t *geometry,    // image dimensions and pixel spacing

	      unsigned char    *image        // complete row-major RGB image buffer

	     #if defined(PROFILE)

	      , unsigned char    *timage

	     #endif

	      )

{

  uint64_t Nrows = geometry->height;

  uint64_t Ncols = geometry->width;

  char schedule_type[20] = {0};

  int nthreads;

 #pragma omp parallel

 #pragma omp single

  nthreads = omp_get_num_threads();

  double threads_atimers[nthreads];

 #if defined(PROFILE)

  double threads_ctimers[nthreads];

  memset( threads_ctimers, 0, nthreads*sizeof(double));

 #endif

  double time_all = CPU_TIME_W;

 #pragma omp parallel

  {

    int    myid         = omp_get_thread_num();

    double mytime       = CPU_TIME_T;

    double my_calc_time = 0;

    if ( options->mode == 0)

      {

       #pragma omp single

	sprintf( schedule_type, "static, max chunk" );

       #pragma omp for nowait

       #include "mandelbrot_loop.omp.h"

      }

    else if (options->mode > 0)

      {

       #pragma omp single

	sprintf( schedule_type, "static, %d-line chunk", options->mode);

       #pragma omp for schedule(static, options->mode) nowait

       #include "mandelbrot_loop.omp.h"

      }

    else

      {

       #pragma omp single

	{

	  options->mode * -1;

	  sprintf( schedule_type, "dynamic, %d points chunk", options->mode );

	}

       #pragma omp for collapse(2) schedule(dynamic, options->mode) nowait

       #include "mandelbrot_loop.omp.h"

      }

    mytime = CPU_TIME_T - mytime;

    threads_atimers[ myid ] = mytime;

   #if defined(PROFILE)

    threads_ctimers[myid] = my_calc_time;

   #endif

  }

  time_all = CPU_TIME_W - time_all;

  double atmin = threads_atimers[0];

  double atmax = threads_atimers[0];

 #if defined(PROFILE)

  double ctmin = threads_ctimers[0];

  double ctmax = threads_ctimers[0];

 #endif

  for ( int t = 1; t < nthreads; t++ )

    {

      atmin = ( atmin < threads_atimers[t] ? atmin : threads_atimers[t] );      

      atmax = ( atmax > threads_atimers[t] ? atmax : threads_atimers[t] );

     #if defined(PROFILE)

      ctmin = ( ctmin < threads_ctimers[t] ? ctmin : threads_ctimers[t] );

      ctmax = ( ctmax > threads_ctimers[t] ? ctmax : threads_ctimers[t] );

     #endif

    }

  printf("\n\nMandelbrot set with %d threads and scheduling %s\n"

	 "took %g sec overall\n"

	 " ·································\n"

	 "thread min time: %g\n"

	 "thread max time: %g\n"

	#if defined(PROFILE)

	 " ·································\n"

	 "thread calc min time: %g\n"

	 "thread calc max time: %g\n"

	#endif

	 " ·································\n",

	 nthreads, schedule_type,

	 time_all,

	 atmin, atmax

	#if defined(PROFILE)

	 , ctmin, ctmax

	#endif

	 );

  return;

}

/* ======================================================================================== */

/*

  Program entry point.

*/

int

main (int    argc,    // number of command-line tokens

      char **argv     // command-line token vector

      )

{

  options_t         options;

  image_geometry_t  geometry;

  render_stats_t    stats;

  unsigned char    *image;

 #if defined(PROFILE)

  unsigned char    *timage;

 #endif

  size_t            image_bytes;

  int               parse_status;

  double            inside_fraction;

  double            average_iterations;

  uint64_t          pixel_count;

  set_default_options (&options);

  parse_status = parse_command_line (argc, argv, &options);

  if (parse_status > 0)

    return EXIT_SUCCESS;

  if (parse_status < 0 || validate_options (&options) != 0

      || make_geometry (&options, &geometry) != 0)

    return EXIT_FAILURE;

  if ((size_t) geometry.width > SIZE_MAX / (size_t) geometry.height

      || (size_t) geometry.width * (size_t) geometry.height > SIZE_MAX / 3u)

    {

      fprintf (stderr, "Image is too large for this platform\n");

      return EXIT_FAILURE;

    }

  image_bytes = (size_t) geometry.width * (size_t) geometry.height * 3u;

  image = (unsigned char *) malloc (image_bytes);

  if (image == NULL)

    {

      fprintf (stderr, "Unable to allocate image buffer (%zu bytes)\n", image_bytes);

      return EXIT_FAILURE;

    }

  fprintf (stderr,

           "Rendering %u x %u pixels, kmax = %u, view = [%.17g, %.17g] x [%.17g, %.17g]\n",

           geometry.width, geometry.height, options.kmax,

           options.xmin, options.xmax, options.ymin, options.ymax);

  if (options.dx_factor != DEFAULT_DX_FACTOR)

    {

      fprintf (stderr,

               "Note: --dx-factor is accepted for assignment-interface compatibility; "

               "this OpenMP version does not use it.\n");

    }

 #if defined(PROFILE)

  timage = (unsigned char *) malloc (image_bytes);

  render_image (&options, &geometry, image, timage);

 #else

  render_image (&options, &geometry, image);

 #endif

  int errio;

  if ( options.format == _MANDEL_PPM_ )

    errio = write_ppm_image (image, geometry.width, geometry.height, options.output_path);

  else if ( options.format == _MANDEL_PNG_ )

    errio = write_png_rgb8 (options.output_path, geometry.width, geometry.height, image);

  else

    {

      fprintf (stderr, "Unknown output format %u\n", options.format);

      errio = -1;

    }

  free (image);

  if ( errio != 0 )

    return EXIT_FAILURE;

 #if defined(PROFILE)

  char buffer[strlen(options.output_path)+5];

  sprintf( buffer, "%s.dat", options.output_path);

  if (write_ppm_image (timage, geometry.width, geometry.height, buffer) != 0)

    {

      free (timage);

      return EXIT_FAILURE;

    }  

  free (timage);

 #endif

  if ( options.mode == 0)

    {

      pixel_count        = (uint64_t) geometry.width * (uint64_t) geometry.height;

      inside_fraction    = (double) stats.inside_pixels / (double) pixel_count;

      average_iterations = (double) stats.total_iterations / (double) pixel_count;

      printf ("output_file              %s\n", options.output_path);

      printf ("width                    %u\n", geometry.width);

      printf ("height                   %u\n", geometry.height);

      printf ("kmax                     %u\n", options.kmax);

      printf ("iteration_checksum       %016llx\n", (unsigned long long) stats.checksum);

      printf ("inside_pixels            %llu\n", (unsigned long long) stats.inside_pixels);

      printf ("inside_fraction          %.17g\n", inside_fraction);

      printf ("average_iterations       %.17g\n", average_iterations);

      printf ("total_iterations         %llu\n", (unsigned long long) stats.total_iterations);

    }

  return EXIT_SUCCESS;

}

/* ────────────────────────────────────────────────────────────────────────── *

│                                                                            │

│ This file is part of the exercises for the Lectures on                     │

│   "Foundations of High Performance Computing"                              │

│ given at                                                                   │

│   Master in HPC and                                                        │

│   Master in Data Science and Scientific Computing                          │

│ @ SISSA, ICTP and University of Trieste                                    │

│                                                                            │

│ contact: luca.tornatore@inaf.it                                            │

│                                                                            │

│     This is free software; you can redistribute it and/or modify           │

│     it under the terms of the GNU General Public License as published by   │

│     the Free Software Foundation; either version 3 of the License, or      │

│     (at your option) any later version.                                    │

│     This code is distributed in the hope that it will be useful,           │

│     but WITHOUT ANY WARRANTY; without even the implied warranty of         │

│     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          │

│     GNU General Public License for more details.                           │

│                                                                            │

│     You should have received a copy of the GNU General Public License      │

│     along with this program.  If not, see <http://www.gnu.org/licenses/>   │

│                                                                            │

* ────────────────────────────────────────────────────────────────────────── */

/*

  OpenMP Mandelbrot renderer for the HPC final-exam exercise.

  This is the OpenMP solution of the serial baseline in ../ : the pixel loop

  (mandelbrot_loop.h) is distributed among threads with a scheduling policy

  chosen at run time through --mode, and per-thread timings are collected to

  discuss load balancing.  Compile with -DPROFILE to also dump, next to the

  image, a map of which thread computed each pixel (<output>.dat), to be

  plotted with plot_thread_map.py.

  Command-line parsing and image writers (PPM or PNG, through the stb header)

  are in mandelbrot_utilities.c; shared declarations are in mandelbrot_omp.h.

  Have fun.

*/

// hint for timing

// you can use either the TIME macro that you find in the examples/timing.h header

// or

//   omp_get_wtime()

// that returns the wall-clock time per thread

//

#include "mandelbrot_omp.h"

/* ========================================================================================

   : ------------------------------------------------------ :

   :  CALCULATION FUNCTIONS                                 :

   : ------------------------------------------------------ :

 */ 

/*

  Compute the Mandelbrot escape iteration for one point c = cr + i ci.  This

  is the direct scalar kernel, clean, keeps the operation count easy to read.

  This can be heavility optimized.

*/

static unsigned int

mandelbrot_escape (double       cr,      // real part of the pixel centre

		   double       ci,      // imaginary part of the pixel centre

		   unsigned int kmax     // maximum number of iterations

		   )

{

           double zr;

           double zi;

           double zr_next;

           double zr2;

           double zi2;

  unsigned int    k;

  zr = 0.0;

  zi = 0.0;

  k = 0;

  while (k < kmax)

    {

      zr2 = zr * zr;

      zi2 = zi * zi;

      if (zr2 + zi2 > 4.0)

        break;

      zr_next = zr2 - zi2 + cr;

      zi = 2.0 * zr * zi + ci;

      zr = zr_next;

      k += 1;

    }

  return k;

}

/*

  Map one iteration count to an RGB colour.  Interior points are black.  The

  outside palette is deliberately simple.

  It might be fun experimenting by advanced colouring or smooth-potential formulas.

  Don't be distracted by this, though

*/

static void

colour_from_iteration (unsigned int    iteration,    // escape iteration count

		       unsigned int    kmax,         // maximum iteration count

		       unsigned char  *rgb           // output triplet: red, green, blue

		       )

{

  double t;

  double one_minus_t;

  if (iteration >= kmax)

    {

      rgb[0] = 0u;

      rgb[1] = 0u;

      rgb[2] = 0u;

      return;

    }

  t = (double) iteration / (double) kmax;

  one_minus_t = 1.0 - t;

  rgb[0] = (unsigned char) (255.0 * 9.0 * one_minus_t * t * t * t);

  rgb[1] = (unsigned char) (255.0 * 15.0 * one_minus_t * one_minus_t * t * t);

  rgb[2] = (unsigned char) (255.0 * 8.5 * one_minus_t * one_minus_t * one_minus_t * t);

}

static void

render_image (options_t        *options,     // view, factors, and kmax

	      image_geometry_t *geometry,    // image dimensions and pixel spacing

	      unsigned char    *image        // complete row-major RGB image buffer

	     #if defined(PROFILE)

	      , unsigned char    *timage

	     #endif

	      )

{

  uint64_t Nrows = geometry->height;

  uint64_t Ncols = geometry->width;

  char schedule_type[20] = {0};

  int nthreads;

 #pragma omp parallel

 #pragma omp single

  nthreads = omp_get_num_threads();

  double threads_atimers[nthreads];

 #if defined(PROFILE)

  double threads_ctimers[nthreads];

  memset( threads_ctimers, nthreads*sizeof(double), 0);

 #endif

  double time_all = CPU_TIME_W;

 #pragma omp parallel

  {

    int    myid         = omp_get_thread_num();

    double mytime       = CPU_TIME_T;

    double my_calc_time = 0;

    if ( options->mode == 0)

      {

	sprintf( schedule_type, "static, max chunk" );

       #pragma omp for nowait

       #include "mandelbrot_loop.h"

      }

    else if (options->mode > 0)

      {

	sprintf( schedule_type, "static, %d-line chunk", options->mode);

       #pragma omp for schedule(static, options->mode) nowait

       #include "mandelbrot_loop.h"

      }

    else

      {

	options->mode *= -1;

	sprintf( schedule_type, "dynamic, %d points chunk", options->mode );

       #pragma omp for collapse(2) schedule(dynamic, options->mode) nowait

       #include "mandelbrot_loop.h"

      }

    mytime = CPU_TIME_T - mytime;

    threads_atimers[ myid ] = mytime;

  }

  time_all = CPU_TIME_W - time_all;

  double atmin = threads_atimers[0];

  double atmax = threads_atimers[0];

 #if defined(PROFILE)

  double ctmin = threads_ctimers[0];

  double ctmax = threads_ctimers[0];

 #endif

  for ( int t = 1; t < nthreads; t++ )

    {

      atmin = ( atmin < threads_atimers[t] ? atmin : threads_atimers[t] );      

      atmax = ( atmax > threads_atimers[t] ? atmax : threads_atimers[t] );

     #if defined(PROFILE)

      ctmin = ( ctmin < threads_ctimers[t] ? ctmin : threads_ctimers[t] );

      ctmax = ( ctmax > threads_ctimers[t] ? ctmax : threads_ctimers[t] );

     #endif

    }

  printf("\n\nMandelbrot set with %d threads and scheduling %s\n"

	 "took %g sec overall\n"

	 " ·································\n"

	 "thread min time: %g\n"

	 "thread max time: %g\n"

	#if defined(PROFILE)

	 " ·································\n"

	 "thread calc min time: %g\n"

	 "thread calc max time: %g\n"

	#endif

	 " ·································\n",

	 nthreads, schedule_type,

	 time_all,

	 atmin, atmax

	#if defined(PROFILE)

	 , ctmin, ctmax

	#endif

	 );

  return;

}

/* ======================================================================================== */

/*

  Program entry point.

*/

int

main (int    argc,    // number of command-line tokens

      char **argv     // command-line token vector

      )

{

  options_t         options;

  image_geometry_t  geometry;

  render_stats_t    stats;

  unsigned char    *image;

 #if defined(PROFILE)

  unsigned char    *timage;

 #endif

  size_t            image_bytes;

  int               parse_status;

  double            inside_fraction;

  double            average_iterations;

  uint64_t          pixel_count;

  set_default_options (&options);

  parse_status = parse_command_line (argc, argv, &options);

  if (parse_status > 0)

    return EXIT_SUCCESS;

  if (parse_status < 0 || validate_options (&options) != 0

      || make_geometry (&options, &geometry) != 0)

    return EXIT_FAILURE;

  if ((size_t) geometry.width > SIZE_MAX / (size_t) geometry.height

      || (size_t) geometry.width * (size_t) geometry.height > SIZE_MAX / 3u)

    {

      fprintf (stderr, "Image is too large for this platform\n");

      return EXIT_FAILURE;

    }

  image_bytes = (size_t) geometry.width * (size_t) geometry.height * 3u;

  image = (unsigned char *) malloc (image_bytes);

  if (image == NULL)

    {

      fprintf (stderr, "Unable to allocate image buffer (%zu bytes)\n", image_bytes);

      return EXIT_FAILURE;

    }

  fprintf (stderr,

           "Rendering %u x %u pixels, kmax = %u, view = [%.17g, %.17g] x [%.17g, %.17g]\n",

           geometry.width, geometry.height, options.kmax,

           options.xmin, options.xmax, options.ymin, options.ymax);

  if (options.dx_factor != DEFAULT_DX_FACTOR)

    {

      fprintf (stderr,

               "Note: --dx-factor is accepted for assignment-interface compatibility; "

               "this OpenMP version does not use it.\n");

    }

 #if defined(PROFILE)

  timage = (unsigned char *) malloc (image_bytes);

  render_image (&options, &geometry, image, timage);

 #else

  render_image (&options, &geometry, image);

 #endif

  int errio;

  if ( options.format == _MANDEL_PPM_ )

    errio = write_ppm_image (image, geometry.width, geometry.height, options.output_path);

  else if ( options.format == _MANDEL_PNG_ )

    errio = write_png_rgb8 (options.output_path, geometry.width, geometry.height, image);

  else

    {

      fprintf (stderr, "Unknown output format %u\n", options.format);

      errio = -1;

    }

  free (image);

  if ( errio != 0 )

    return EXIT_FAILURE;

 #if defined(PROFILE)

  char buffer[strlen(options.output_path)+5];

  sprintf( buffer, "%s.dat", options.output_path);

  if (write_ppm_image (timage, geometry.width, geometry.height, buffer) != 0)

    {

      free (timage);

      return EXIT_FAILURE;

    }  

  free (timage);

 #endif

  if ( options.mode == 0)

    {

      pixel_count        = (uint64_t) geometry.width * (uint64_t) geometry.height;

      inside_fraction    = (double) stats.inside_pixels / (double) pixel_count;

      average_iterations = (double) stats.total_iterations / (double) pixel_count;

      printf ("output_file              %s\n", options.output_path);

      printf ("width                    %u\n", geometry.width);

      printf ("height                   %u\n", geometry.height);

      printf ("kmax                     %u\n", options.kmax);

      printf ("iteration_checksum       %016llx\n", (unsigned long long) stats.checksum);

      printf ("inside_pixels            %llu\n", (unsigned long long) stats.inside_pixels);

      printf ("inside_fraction          %.17g\n", inside_fraction);

      printf ("average_iterations       %.17g\n", average_iterations);

      printf ("total_iterations         %llu\n", (unsigned long long) stats.total_iterations);

    }

  return EXIT_SUCCESS;

}