init: HW2

part2/Makefile

@@ -0,0 +1,34 @@
+APP_NAME = mandelbrot
+OBJDIR = objs
+COMMONDIR = ./common
+
+CXX = g++ -m64
+CXXFLAGS = -I$(COMMONDIR) -I$(OBJDIR) -O3 -std=c++17 -Wall
+LDLIBS = -lm -lpthread
+
+PPM_CXX = $(COMMONDIR)/ppm.cpp
+PPM_OBJ = $(addprefix $(OBJDIR)/, $(subst $(COMMONDIR)/,, $(PPM_CXX:.cpp=.o)))
+
+
+default: $(APP_NAME)
+
+dirs:
+	mkdir -p $(OBJDIR)/
+.PHONY: dirs
+
+clean:
+	$(RM) -r $(OBJDIR) *.ppm *~ $(APP_NAME)
+.PHONY: clean
+
+OBJS = $(OBJDIR)/main.o $(OBJDIR)/mandelbrot_serial.o $(OBJDIR)/mandelbrot_thread.o $(PPM_OBJ)
+
+$(APP_NAME): dirs $(OBJS)
+	$(CXX) $(CXXFLAGS) -o $@ $(OBJS) $(LDLIBS)
+
+$(OBJDIR)/%.o: %.cpp
+	$(CXX) $< $(CXXFLAGS) -c -o $@
+
+$(OBJDIR)/%.o: $(COMMONDIR)/%.cpp
+	$(CXX) $< $(CXXFLAGS) -c -o $@
+
+$(OBJDIR)/main.o: $(COMMONDIR)/cycle_timer.h

part2/common/cycle_timer.h

@@ -0,0 +1,195 @@
+#ifndef SYRAH_CYCLE_TIMER_H
+#define SYRAH_CYCLE_TIMER_H
+
+#if defined(__APPLE__)
+#if defined(__x86_64__)
+#include <sys/sysctl.h>
+#else
+#include <mach/mach.h>
+#include <mach/mach_time.h>
+#endif // __x86_64__ or not
+
+#include <cstdio>  // fprintf
+#include <cstdlib> // exit
+
+#elif _WIN32
+#include <time.h>
+#include <windows.h>
+#else
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <sys/time.h>
+#endif
+
+// This uses the cycle counter of the processor.  Different
+// processors in the system will have different values for this.  If
+// you process moves across processors, then the delta time you
+// measure will likely be incorrect.  This is mostly for fine
+// grained measurements where the process is likely to be on the
+// same processor.  For more global things you should use the
+// Time interface.
+
+// Also note that if you processors' speeds change (i.e. processors
+// scaling) or if you are in a heterogenous environment, you will
+// likely get spurious results.
+class CycleTimer
+{
+  public:
+    using SysClock = unsigned long long;
+
+    //////////
+    // Return the current CPU time, in terms of clock ticks.
+    // Time zero is at some arbitrary point in the past.
+    static SysClock current_ticks()
+    {
+#if defined(__APPLE__) && !defined(__x86_64__)
+        return mach_absolute_time();
+#elif defined(_WIN32)
+        LARGE_INTEGER qwTime;
+        QueryPerformanceCounter(&qwTime);
+        return qwTime.QuadPart;
+#elif defined(__x86_64__)
+        unsigned int a, d;
+        asm volatile("rdtsc" : "=a"(a), "=d"(d));
+        return static_cast<unsigned long long>(a) | (static_cast<unsigned long long>(d) << 32);
+#elif defined(__ARM_NEON__) && 0 // mrc requires superuser.
+        unsigned int val;
+        asm volatile("mrc p15, 0, %0, c9, c13, 0" : "=r"(val));
+        return val;
+#else
+        timespec spec;
+        clock_gettime(CLOCK_THREAD_CPUTIME_ID, &spec);
+        return CycleTimer::SysClock(static_cast<float>(spec.tv_sec) * 1e9
+                                    + static_cast<float>(spec.tv_nsec));
+#endif
+    }
+
+    //////////
+    // Return the current CPU time, in terms of seconds.
+    // This is slower than current_ticks().  Time zero is at
+    // some arbitrary point in the past.
+    static double current_seconds()
+    {
+        return static_cast<double>(current_ticks()) * seconds_per_tick();
+    }
+
+    //////////
+    // Return the conversion from seconds to ticks.
+    static double ticks_per_second()
+    {
+        return 1.0 / seconds_per_tick();
+    }
+
+    static const char *tick_units()
+    {
+#if defined(__APPLE__) && !defined(__x86_64__)
+        return "ns";
+#elif defined(__WIN32__) || defined(__x86_64__)
+        return "cycles";
+#else
+        return "ns"; // clock_gettime
+#endif
+    }
+
+    //////////
+    // Return the conversion from ticks to seconds.
+    static double seconds_per_tick()
+    {
+        static bool initialized = false;
+        static double seconds_per_tick_val;
+        if (initialized)
+            return seconds_per_tick_val;
+#if defined(__APPLE__)
+#ifdef __x86_64__
+        int args[] = {CTL_HW, HW_CPU_FREQ};
+        unsigned int Hz;
+        size_t len = sizeof(Hz);
+        if (sysctl(args, 2, &Hz, &len, NULL, 0) != 0)
+        {
+            fprintf(stderr, "Failed to initialize seconds_per_tick_val!\n");
+            exit(-1);
+        }
+        seconds_per_tick_val = 1.0 / (double)Hz;
+#else
+        mach_timebase_info_data_t time_info;
+        mach_timebase_info(&time_info);
+
+        // Scales to nanoseconds without 1e-9f
+        seconds_per_tick_val
+            = (1e-9 * static_cast<double>(time_info.numer)) / static_cast<double>(time_info.denom);
+#endif // x86_64 or not
+#elif defined(_WIN32)
+        LARGE_INTEGER qwTicksPerSec;
+        QueryPerformanceFrequency(&qwTicksPerSec);
+        seconds_per_tick_val = 1.0 / static_cast<double>(qwTicksPerSec.QuadPart);
+#else
+        FILE *fp = fopen("/proc/cpuinfo", "r");
+        char input[1024]; // NOLINT(modernize-avoid-c-arrays): Extensively working with C libraries.
+        if (!fp)
+        {
+            fprintf(stderr, "%s failed: couldn't find /proc/cpuinfo.", __PRETTY_FUNCTION__);
+            exit(-1);
+        }
+        // In case we don't find it, e.g. on the N900
+        seconds_per_tick_val = 1e-9;
+        while (!feof(fp) && fgets(input, 1024, fp))
+        {
+            // NOTE(boulos): Because reading cpuinfo depends on dynamic
+            // frequency scaling it's better to read the @ sign first
+            float ghz, mhz;
+            if (strstr(input, "model name"))
+            {
+                char *at_sign = strstr(input, "@");
+                if (at_sign)
+                {
+                    char *after_at = at_sign + 1;
+                    char *ghz_str = strstr(after_at, "ghz");
+                    char *mhz_str = strstr(after_at, "mhz");
+                    if (ghz_str)
+                    {
+                        *ghz_str = '\0';
+                        if (1 == sscanf(after_at, "%f", &ghz))
+                        {
+                            // printf("GHz = %f\n", ghz);
+                            seconds_per_tick_val = 1e-9f / ghz;
+                            break;
+                        }
+                    }
+                    else if (mhz_str)
+                    {
+                        *mhz_str = '\0';
+                        if (1 == sscanf(after_at, "%f", &mhz))
+                        {
+                            // printf("MHz = %f\n", mhz);
+                            seconds_per_tick_val = 1e-6f / mhz;
+                            break;
+                        }
+                    }
+                }
+            }
+            else if (1 == sscanf(input, "cpu MHz : %f", &mhz))
+            {
+                // printf("MHz = %f\n", mhz);
+                seconds_per_tick_val = 1e-6f / mhz;
+                break;
+            }
+        }
+        fclose(fp);
+#endif
+
+        initialized = true;
+        return seconds_per_tick_val;
+    }
+
+    //////////
+    // Return the conversion from ticks to milliseconds.
+    static double ms_per_tick()
+    {
+        return seconds_per_tick() * 1000.0;
+    }
+
+    CycleTimer() = delete;
+};
+
+#endif // SYRAH_CYCLE_TIMER_H

part2/common/ppm.cpp

@@ -0,0 +1,38 @@
+#include <algorithm>
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+
+void write_ppm_image(int *data, int width, int height, const char *filename, int max_iterations)
+{
+    FILE *fp = fopen(filename, "wb");
+    if (!fp)
+    {
+        perror("fopen() failed");
+        exit(EXIT_FAILURE);
+    }
+
+    // write ppm header
+    fprintf(fp, "P6\n");
+    fprintf(fp, "%d %d\n", width, height);
+    fprintf(fp, "255\n");
+
+    for (int i = 0; i < width * height; ++i)
+    {
+
+        // Clamp iteration count for this pixel, then scale the value
+        // to 0-1 range.  Raise resulting value to a power (<1) to
+        // increase brightness of low iteration count
+        // pixels. a.k.a. Make things look cooler.
+
+        float mapped = std::pow(
+            std::min(static_cast<float>(max_iterations), static_cast<float>(data[i])) / 256.f, .5f);
+
+        // convert back into 0-255 range, 8-bit channels
+        unsigned char result = static_cast<unsigned char>(255.f * mapped);
+        for (int j = 0; j < 3; ++j)
+            fputc(result, fp);
+    }
+    fclose(fp);
+    printf("Wrote image file %s\n", filename);
+}

part2/main.cpp

@@ -0,0 +1,190 @@
+#include <algorithm>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <getopt.h>
+
+#include "cycle_timer.h"
+
+extern void mandelbrot_serial(float x0,
+                              float y0,
+                              float x1,
+                              float y1,
+                              int width,
+                              int height,
+                              int start_row,
+                              int num_rows,
+                              int max_iterations,
+                              int *output);
+
+extern void mandelbrot_thread(int num_threads,
+                              float x0,
+                              float y0,
+                              float x1,
+                              float y1,
+                              int width,
+                              int height,
+                              int max_iterations,
+                              int *output);
+
+extern void
+write_ppm_image(int *data, int width, int height, const char *filename, int max_iterations);
+
+void scale_and_shift(
+    float &x0, float &x1, float &y0, float &y1, float scale, float shift_x, float shift_y)
+{
+
+    x0 *= scale;
+    x1 *= scale;
+    y0 *= scale;
+    y1 *= scale;
+    x0 += shift_x;
+    x1 += shift_x;
+    y0 += shift_y;
+    y1 += shift_y;
+}
+
+void usage(const char *progname)
+{
+    printf("Usage: %s [options]\n", progname);
+    printf("Program Options:\n");
+    printf("  -t  --threads <N>  Use N threads (Default = 2)\n");
+    printf("  -v  --view <INT>   Use specified view settings (Default = 1)\n");
+    printf("  -?  --help         This message\n");
+}
+
+bool verify_result(int *gold, int *result, int width, int height)
+{
+
+    int i, j;
+
+    for (i = 0; i < height; i++)
+    {
+        for (j = 0; j < width; j++)
+        {
+            if (gold[(i * width) + j] != result[(i * width) + j])
+            {
+                printf("Mismatch : [%d][%d], Expected : %d, Actual : %d\n", i, j,
+                       gold[(i * width) + j], result[(i * width) + j]);
+                return false;
+            }
+        }
+    }
+
+    return true;
+}
+
+int main(int argc, char **argv)
+{
+
+    const unsigned int width = 1600;
+    const unsigned int height = 1200;
+    const int max_iterations = 256;
+    int num_threads = 2;
+
+    float x0 = -2;
+    float x1 = 1;
+    float y0 = -1;
+    float y1 = 1;
+
+    // parse commandline options ////////////////////////////////////////////
+    int opt;
+    // NOLINTNEXTLINE(modernize-avoid-c-arrays): Required by C function.
+    static struct option long_options[] = {{"threads", 1, nullptr, 't'},
+                                           {"view", 1, nullptr, 'v'},
+                                           {"help", 0, nullptr, '?'},
+                                           {nullptr, 0, nullptr, 0}};
+
+    while ((opt = getopt_long(argc, argv, "t:v:?", long_options, nullptr)) != EOF)
+    {
+
+        switch (opt)
+        {
+            case 't':
+            {
+                num_threads = atoi(optarg);
+                break;
+            }
+            case 'v':
+            {
+                int view_index = atoi(optarg);
+                // change view settings
+                if (view_index == 2)
+                {
+                    float scale_value = .015f;
+                    float shift_x = -.986f;
+                    float shift_y = .30f;
+                    scale_and_shift(x0, x1, y0, y1, scale_value, shift_x, shift_y);
+                }
+                else if (view_index > 1)
+                {
+                    fprintf(stderr, "Invalid view index\n");
+                    return 1;
+                }
+                break;
+            }
+            case '?':
+            default:
+                usage(argv[0]);
+                return 1;
+        }
+    }
+    // end parsing of commandline options
+
+    int *output_serial = new int[width * height];
+    int *output_thread = new int[width * height];
+
+    //
+    // Run the serial implementation.  Run the code three times and
+    // take the minimum to get a good estimate.
+    //
+
+    double min_serial = 1e30;
+    for (int i = 0; i < 5; ++i)
+    {
+        memset(output_serial, 0, width * height * sizeof(int));
+        double start_time = CycleTimer::current_seconds();
+        mandelbrot_serial(x0, y0, x1, y1, width, height, 0, height, max_iterations, output_serial);
+        double end_time = CycleTimer::current_seconds();
+        min_serial = std::min(min_serial, end_time - start_time);
+    }
+
+    printf("[mandelbrot serial]:\t\t[%.3f] ms\n", min_serial * 1000);
+    write_ppm_image(output_serial, width, height, "mandelbrot-serial.ppm", max_iterations);
+
+    //
+    // Run the threaded version
+    //
+
+    double min_thread = 1e30;
+    for (int i = 0; i < 5; ++i)
+    {
+        memset(output_thread, 0, width * height * sizeof(int));
+        double start_time = CycleTimer::current_seconds();
+        mandelbrot_thread(num_threads, x0, y0, x1, y1, width, height, max_iterations,
+                          output_thread);
+        double end_time = CycleTimer::current_seconds();
+        min_thread = std::min(min_thread, end_time - start_time);
+    }
+
+    printf("[mandelbrot thread]:\t\t[%.3f] ms\n", min_thread * 1000);
+    write_ppm_image(output_thread, width, height, "mandelbrot-thread.ppm", max_iterations);
+
+    if (!verify_result(output_serial, output_thread, width, height))
+    {
+        printf("Error : Output from threads does not match serial output\n");
+
+        delete[] output_serial;
+        delete[] output_thread;
+
+        return 1;
+    }
+
+    // compute speedup
+    printf("\t\t\t\t(%.2fx speedup from %d threads)\n", min_serial / min_thread, num_threads);
+
+    delete[] output_serial;
+    delete[] output_thread;
+
+    return 0;
+}

part2/mandelbrot_serial.cpp

@@ -0,0 +1,102 @@
+/*
+
+  Note: This code was modified from example code
+  originally provided by Intel.  To comply with Intel's open source
+  licensing agreement, their copyright is retained below.
+
+  -----------------------------------------------------------------
+
+  Copyright (c) 2010-2011, Intel Corporation
+  All rights reserved.
+
+  Redistribution and use in source and binary forms, with or without
+  modification, are permitted provided that the following conditions are
+  met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+
+    * Neither the name of Intel Corporation nor the names of its
+      contributors may be used to endorse or promote products derived from
+      this software without specific prior written permission.
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+   IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+   TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+   PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+   OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+namespace
+{
+
+int mandel(float c_re, float c_im, int count)
+{
+    float z_re = c_re, z_im = c_im;
+    int i;
+    for (i = 0; i < count; ++i)
+    {
+
+        if (z_re * z_re + z_im * z_im > 4.f)
+            break;
+
+        float new_re = (z_re * z_re) - (z_im * z_im);
+        float new_im = 2.f * z_re * z_im;
+        z_re = c_re + new_re;
+        z_im = c_im + new_im;
+    }
+
+    return i;
+}
+
+} // namespace
+
+//
+// MandelbrotSerial --
+//
+// Compute an image visualizing the mandelbrot set.  The resulting
+// array contains the number of iterations required before the complex
+// number corresponding to a pixel could be rejected from the set.
+//
+// * x0, y0, x1, y1 describe the complex coordinates mapping
+//   into the image viewport.
+// * width, height describe the size of the output image
+// * startRow, totalRows describe how much of the image to compute
+void mandelbrot_serial(float x0,
+                       float y0,
+                       float x1,
+                       float y1,
+                       int width,
+                       int height,
+                       int start_row,
+                       int total_rows,
+                       int max_iterations,
+                       int *output)
+{
+    float dx = (x1 - x0) / (float)width;
+    float dy = (y1 - y0) / (float)height;
+
+    int end_row = start_row + total_rows;
+
+    for (int j = start_row; j < end_row; j++)
+    {
+        for (int i = 0; i < width; ++i)
+        {
+            float x = x0 + ((float)i * dx);
+            float y = y0 + ((float)j * dy);
+
+            int index = ((j * width) + i);
+            output[index] = mandel(x, y, max_iterations);
+        }
+    }
+}

part2/mandelbrot_thread.cpp

@@ -0,0 +1,107 @@
+#include <array>
+#include <cstdio>
+#include <cstdlib>
+#include <thread>
+
+struct WorkerArgs
+{
+    float x0, x1;
+    float y0, y1;
+    unsigned int width;
+    unsigned int height;
+    int maxIterations;
+    int *output;
+    int threadId;
+    int numThreads;
+};
+
+extern void mandelbrot_serial(float x0,
+                              float y0,
+                              float x1,
+                              float y1,
+                              int width,
+                              int height,
+                              int start_row,
+                              int num_rows,
+                              int max_iterations,
+                              int *output);
+
+//
+// worker_thread_start --
+//
+// Thread entrypoint.
+void worker_thread_start(WorkerArgs *const args)
+{
+
+    // TODO FOR PP STUDENTS: Implement the body of the worker
+    // thread here. Each thread could make a call to mandelbrot_serial()
+    // to compute a part of the output image. For example, in a
+    // program that uses two threads, thread 0 could compute the top
+    // half of the image and thread 1 could compute the bottom half.
+    // Of course, you can copy mandelbrot_serial() to this file and
+    // modify it to pursue a better performance.
+
+    printf("Hello world from thread %d\n", args->threadId);
+}
+
+//
+// mandelbrot_thread --
+//
+// Multi-threaded implementation of mandelbrot set image generation.
+// Threads of execution are created by spawning std::threads.
+void mandelbrot_thread(int num_threads,
+                       float x0,
+                       float y0,
+                       float x1,
+                       float y1,
+                       int width,
+                       int height,
+                       int max_iterations,
+                       int *output)
+{
+    static constexpr int max_threads = 32;
+
+    if (num_threads > max_threads)
+    {
+        fprintf(stderr, "Error: Max allowed threads is %d\n", max_threads);
+        exit(1);
+    }
+
+    // Creates thread objects that do not yet represent a thread.
+    std::array<std::thread, max_threads> workers;
+    std::array<WorkerArgs, max_threads> args = {};
+
+    for (int i = 0; i < num_threads; i++)
+    {
+        // TODO FOR PP STUDENTS: You may or may not wish to modify
+        // the per-thread arguments here.  The code below copies the
+        // same arguments for each thread
+        args[i].x0 = x0;
+        args[i].y0 = y0;
+        args[i].x1 = x1;
+        args[i].y1 = y1;
+        args[i].width = width;
+        args[i].height = height;
+        args[i].maxIterations = max_iterations;
+        args[i].numThreads = num_threads;
+        args[i].output = output;
+
+        args[i].threadId = i;
+    }
+
+    // Spawn the worker threads.  Note that only numThreads-1 std::threads
+    // are created and the main application thread is used as a worker
+    // as well.
+    for (int i = 1; i < num_threads; i++)
+    {
+        workers[i] = std::thread(worker_thread_start, &args[i]);
+    }
+
+    worker_thread_start(&args[0]);
+
+    // join worker threads
+    for (int i = 1; i < num_threads; i++)
+    {
+        workers[i].join();
+    }
+}