并行与分布式计算（2）

#include <stdio.h>
#include <pthread.h>
// Use this code to time your threads
// 我加入的代码在126~134行和164~174行
#include "CycleTimer.h"

/*

  15418 Spring 2012 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.
*/

// Core computation of Mandelbrot set membershop
// Iterate complex number c to determine whether it diverges
static inline 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;
}

//
// 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, endRow describe how much of the image to compute
void mandelbrotSerial(
	float x0, float y0, float x1, float y1,
	int width, int height,
	int startRow, int endRow,
	int maxIterations,
	int output[])
{
	float dx = (x1 - x0) / width;
	float dy = (y1 - y0) / height;

	for (int j = startRow; j < endRow; j++)
	{
		for (int i = 0; i < width; ++i)
		{
			float x = x0 + i * dx;
			float y = y0 + j * dy;

			int index = (j * width + i);
			output[index] = mandel(x, y, maxIterations);
		}
	}
}

// Struct for passing arguments to thread routine
typedef struct
{
	float x0, x1;
	float y0, y1;
	unsigned int width;
	unsigned int height;
	int maxIterations;
	int *output;
	int threadId;
	int numThreads;
} WorkerArgs;

//
// workerThreadStart --
//
// Thread entrypoint.
void *workerThreadStart(void *threadArgs)
{

	WorkerArgs *args = static_cast<WorkerArgs *>(threadArgs);

	// TODO: Implement worker thread here.

	printf("Hello world from thread %d\n", args->threadId);
	//以下是我加的第一部分内容
	mandelbrotSerial(
		args->x0, args->y0, args->x1, args->y1,
		args->width, args->height,
		args->height / args->numThreads * args->threadId,																//startRow
		args->threadId == args->numThreads - 1 ? args->height : args->height / args->numThreads * (args->threadId + 1), //endRow
		args->maxIterations,
		args->output);
	//增加完毕，只需要调用串行部分的代码即可
	return NULL;
}

//
// MandelbrotThread --
//
// Multi-threaded implementation of mandelbrot set image generation.
// Multi-threading performed via pthreads.
void mandelbrotThread(
	int numThreads,
	float x0, float y0, float x1, float y1,
	int width, int height,
	int maxIterations, int output[])
{
	const static int MAX_THREADS = 32;

	if (numThreads > MAX_THREADS)
	{
		fprintf(stderr, "Error: Max allowed threads is %d\n", MAX_THREADS);
		exit(1);
	}

	pthread_t workers[MAX_THREADS];
	WorkerArgs args[MAX_THREADS];

	for (int i = 0; i < numThreads; i++)
	{
		// TODO: Set thread arguments here.
		args[i].threadId = i;
		//以下是我加的第二段内容
		args[i].x0 = x0;
		args[i].x1 = x1;
		args[i].y0 = y0;
		args[i].y1 = y1;
		args[i].width = width;
		args[i].height = height;
		args[i].maxIterations = maxIterations;
		args[i].output = output;
		args[i].numThreads = numThreads;
		//增加完毕，把所有信息传给worker线程即可
	}

	// Fire up the worker threads.  Note that numThreads-1 pthreads
	// are created and the main app thread is used as a worker as
	// well.

	for (int i = 1; i < numThreads; i++)
		pthread_create(&workers[i], NULL, workerThreadStart, &args[i]);

	workerThreadStart(&args[0]);

	// wait for worker threads to complete
	for (int i = 1; i < numThreads; i++)
		pthread_join(workers[i], NULL);
}

#include <stdio.h>
#include <algorithm>
#include <getopt.h>

#include "CycleTimer.h"

extern void mandelbrotSerial(
    float x0, float y0, float x1, float y1,
    int width, int height,
    int startRow, int endRow,
    int maxIterations,
    int output[]);

extern void mandelbrotThread(
    int numThreads,
    float x0, float y0, float x1, float y1,
    int width, int height,
    int maxIterations,
    int output[]);

extern void writePPMImage(
    int *data,
    int width, int height,
    const char *filename,
    int maxIterations);

void scaleAndShift(float &x0, float &x1, float &y0, float &y1,
                   float scale,
                   float shiftX, float shiftY)
{

    x0 *= scale;
    x1 *= scale;
    y0 *= scale;
    y1 *= scale;
    x0 += shiftX;
    x1 += shiftX;
    y0 += shiftY;
    y1 += shiftY;
}

void usage(const char *progname)
{
    printf("Usage: %s [options]\n", progname);
    printf("Program Options:\n");
    printf("  -t  --threads <N>  Use N threads\n");
    printf("  -v  --view <INT>   Use specified view settings (1-6)\n");
    printf("  -?  --help         This message\n");
}

bool verifyResult(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 0;
            }
        }
    }

    return 1;
}

#define VIEWCNT 6

int main(int argc, char **argv)
{

    const int width = 1440;
    const int height = 900;
    const int maxIterations = 256;
    int numThreads = 2;

    float x0 = -2.167;
    float x1 = 1.167;
    float y0 = -1;
    float y1 = 1;

    // Support VIEWCNT views
    float scaleValues[VIEWCNT + 1] = {1.0f, 1.0f, 0.015f, 0.02f, 0.02f, 0.02f, 0.002f};
    float shiftXs[VIEWCNT + 1] = {0.0f, 0.0f, -0.98f, 0.35f, 0.0f, -1.5f, -1.4f};
    float shiftYs[VIEWCNT + 1] = {0.0f, 0.0f, 0.30f, 0.05f, 0.73f, 0.0f, 0.0f};

    // parse commandline options ////////////////////////////////////////////
    int opt;
    static struct option long_options[] = {
        {"threads", 1, 0, 't'},
        {"view", 1, 0, 'v'},
        {"help", 0, 0, '?'},
        {0, 0, 0, 0}};

    while ((opt = getopt_long(argc, argv, "t:v:?", long_options, NULL)) != EOF)
    {

        switch (opt)
        {
        case 't':
        {
            numThreads = atoi(optarg);
            break;
        }
        case 'v':
        {
            int viewIndex = atoi(optarg);
            // change view settings
            if (viewIndex >= 1 && viewIndex <= VIEWCNT)
            {
                float scaleValue = scaleValues[viewIndex];
                float shiftX = shiftXs[viewIndex];
                float shiftY = shiftYs[viewIndex];
                scaleAndShift(x0, x1, y0, y1, scaleValue, shiftX, shiftY);
            }
            else
            {
                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.
    //
    memset(output_serial, 0, width * height * sizeof(int));
    double minSerial = 1e30;
    for (int i = 0; i < 3; ++i)
    {
        double startTime = CycleTimer::currentSeconds();
        mandelbrotSerial(x0, y0, x1, y1, width, height, 0, height, maxIterations, output_serial);
        double endTime = CycleTimer::currentSeconds();
        minSerial = std::min(minSerial, endTime - startTime);
    }

    printf("[mandelbrot serial]:\t\t[%.3f] ms\n", minSerial * 1000);
    writePPMImage(output_serial, width, height, "mandelbrot-serial.ppm", maxIterations);

    //
    // Run the threaded version
    //
    memset(output_thread, 0, width * height * sizeof(int));
    double minThread = 1e30;
    for (int i = 0; i < 5; ++i)
    {
        double startTime = CycleTimer::currentSeconds();
        mandelbrotThread(numThreads, x0, y0, x1, y1, width, height, maxIterations, output_thread);
        double endTime = CycleTimer::currentSeconds();
        minThread = std::min(minThread, endTime - startTime);
    }

    printf("[mandelbrot thread]:\t\t[%.3f] ms\n", minThread * 1000);
    writePPMImage(output_thread, width, height, "mandelbrot-thread.ppm", maxIterations);

    if (!verifyResult(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", minSerial / minThread, numThreads);

    delete[] output_serial;
    delete[] output_thread;

    return 0;
}