/* -*- Mode: C; c-basic-offset:4 ; indent-tabs-mode:nil ; -*- */

/*

 *  (C) 2001 by Argonne National Laboratory.

 *      See COPYRIGHT in top-level directory.

 */

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include <string.h>

#ifdef HAVE_WINDOWS_H

#include <process.h> /* getpid() */

#include <winsock2.h>

#else

#include <unistd.h>

#include <sys/types.h>

#include <sys/socket.h>

#include <sys/select.h>

#include <netinet/in.h>

#include <errno.h>

#endif

#include "mpi.h"

/* definitions */

#define PROMPT 1

#define USE_PPM /* define to output a color image, otherwise output is a grayscale pgm image */

#ifndef MIN

#define MIN(x, y)	(((x) < (y))?(x):(y))

#endif

#ifndef MAX

#define MAX(x, y)	(((x) > (y))?(x):(y))

#endif

#define DEFAULT_NUM_SLAVES 3 /* default number of children to spawn */

#define DEFAULT_PORT 7470   /* default port to listen on */

#define NOVALUE 99999       /* indicates a parameter is as of yet unspecified */

#define MAX_ITERATIONS 10000 /* maximum 'depth' to look for mandelbrot value */

#define INFINITE_LIMIT 4.0  /* evalue when fractal is considered diverging */

#define MAX_X_VALUE 2.0     /* the highest x can go for the mandelbrot, usually 1.0 */

#define MIN_X_VALUE -2.0    /* the lowest x can go for the mandelbrot, usually -1.0 */

#define MAX_Y_VALUE 2.0     /* the highest y can go for the mandelbrot, usually 1.0 */

#define MIN_Y_VALUE -2.0    /* the lowest y can go for the mandelbrot, usually -1.0 */

#define IDEAL_ITERATIONS 100 /* my favorite 'depth' to default to */

/* the ideal default x and y are currently the same as the max area */

#define IDEAL_MAX_X_VALUE 1.0

#define IDEAL_MIN_X_VALUE -1.0

#define IDEAL_MAX_Y_VALUE 1.0

#define IDEAL_MIN_Y_VALUE -1.0

#define MAX_WIDTH 2048       /* maximum size of image, across, in pixels */

#define MAX_HEIGHT 2048      /* maximum size of image, down, in pixels */

#define IDEAL_WIDTH 768       /* my favorite size of image, across, in pixels */

#define IDEAL_HEIGHT 768      /* my favorite size of image, down, in pixels */

#define RGBtocolor_t(r,g,b) ((color_t)(((unsigned char)(r)|((unsigned short)((unsigned char)(g))<<8))|(((unsigned long)(unsigned char)(b))<<16)))

#define getR(r) ((int)((r) & 0xFF))

#define getG(g) ((int)((g>>8) & 0xFF))

#define getB(b) ((int)((b>>16) & 0xFF))

typedef int color_t;

typedef struct complex_t

{

  /* real + imaginary * sqrt(-1) i.e.   x + iy */

  double real, imaginary;

} complex_t;

/* prototypes */

void read_mand_args(int argc, char *argv[], int *o_max_iterations,

		    int *o_pixels_across, int *o_pixels_down,

		    double *o_x_min, double *o_x_max,

		    double *o_y_min, double *o_y_max, int *o_julia,

		    double *o_julia_real_x, double *o_julia_imaginary_y,

		    double *o_divergent_limit, int *o_alternate,

		    char *filename, int *num_colors, int *use_stdin,

		    int *save_image, int *num_children);

void check_mand_params(int *m_max_iterations,

		       int *m_pixels_across, int *m_pixels_down,

		       double *m_x_min, double *m_x_max,

		       double *m_y_min, double *m_y_max,

		       double *m_divergent_limit,

		       int *m_num_children);

void check_julia_params(double *julia_x_constant, double *julia_y_constant);

int additive_mandelbrot_point(complex_t coord_point, 

			    complex_t c_constant,

			    int Max_iterations, double divergent_limit);

int additive_mandelbrot_point(complex_t coord_point, 

			    complex_t c_constant,

			    int Max_iterations, double divergent_limit);

int subtractive_mandelbrot_point(complex_t coord_point, 

			    complex_t c_constant,

			    int Max_iterations, double divergent_limit);

complex_t exponential_complex(complex_t in_complex);

complex_t multiply_complex(complex_t in_one_complex, complex_t in_two_complex);

complex_t divide_complex(complex_t in_one_complex, complex_t in_two_complex);

complex_t inverse_complex(complex_t in_complex);

complex_t add_complex(complex_t in_one_complex, complex_t in_two_complex);

complex_t subtract_complex(complex_t in_one_complex, complex_t in_two_complex);

double absolute_complex(complex_t in_complex);

void dumpimage(char *filename, int in_grid_array[], int in_pixels_across, int in_pixels_down,

	  int in_max_pixel_value, char input_string[], int num_colors, color_t colors[]);

int single_mandelbrot_point(complex_t coord_point, 

			    complex_t c_constant,

			    int Max_iterations, double divergent_limit);

color_t getColor(double fraction, double intensity);

int Make_color_array(int num_colors, color_t colors[]);

void output_data(int *in_grid_array, int coord[4], int *out_grid_array, int width, int height);

void PrintUsage(void);

static int sock_write(int sock, void *buffer, int length);

static int sock_read(int sock, void *buffer, int length);

#ifdef USE_PPM

const char *default_filename = "pmandel.ppm";

#else

const char *default_filename = "pmandel.pgm";

#endif

/* Solving the Mandelbrot set

   Set a maximum number of iterations to calculate before forcing a terminate

   in the Mandelbrot loop.

*/

/* Command-line parameters are all optional, and include:

   -xmin # -xmax #      Limits for the x-axis for calculating and plotting

   -ymin # -ymax #      Limits for the y-axis for calculating and plotting

   -xscale # -yscale #  output pixel width and height

   -depth #             Number of iterations before we give up on a given

                        point and move on to calculate the next

   -diverge #           Criteria for assuming the calculation has been

                        "solved"-- for each point z, when

                             abs(z=z^2+c) > diverge_value

   -julia        Plot a Julia set instead of a Mandelbrot set

   -jx # -jy #   The Julia point that defines the set

   -alternate #    Plot an alternate Mandelbrot equation (varient 1 or 2 so far)

*/

int myid;

int use_stdin = 0;

int sock;

void swap(int *i, int *j)

{

    int x;

    x = *i;

    *i = *j;

    *j = x;

}

typedef struct header_t

{

    int data[5];

} header_t;

int main(int argc, char *argv[])

{

    complex_t coord_point, julia_constant;

    double x_max, x_min, y_max, y_min, x_resolution, y_resolution;

    double divergent_limit;

    char file_message[160];

    char filename[100];

    int icount, imax_iterations;

    int ipixels_across, ipixels_down;

    int i, j, k, julia, alternate_equation;

    int imin, imax, jmin, jmax;

    int work[5];

    header_t *result = NULL;

    /* make an integer array of size [N x M] to hold answers. */

    int *in_grid_array, *out_grid_array = NULL;

    int numprocs;

    int  namelen;

    char processor_name[MPI_MAX_PROCESSOR_NAME];

    int num_colors;

    color_t *colors = NULL;

    MPI_Status status;

    int listener;

    int save_image = 0;

    int optval;

    int num_children = DEFAULT_NUM_SLAVES;

    int master = 1;

    MPI_Comm parent, *child_comm = NULL;

    MPI_Request *child_request = NULL;

    int error_code, error;

    char error_str[MPI_MAX_ERROR_STRING];

    int length;

    int index;

    int pid;

    MPI_Init(&argc, &argv);

    MPI_Comm_size(MPI_COMM_WORLD, &numprocs);

    MPI_Comm_rank(MPI_COMM_WORLD, &myid);

    MPI_Comm_get_parent(&parent);

    MPI_Get_processor_name(processor_name, &namelen);

    pid = getpid();

    if (parent == MPI_COMM_NULL)

    {

	if (numprocs > 1)

	{

	    printf("Error: only one process allowed for the master.\n");

	    PrintUsage();

	    error_code = MPI_Abort(MPI_COMM_WORLD, -1);

	    exit(error_code);

	}

	printf("Welcome to the Mandelbrot/Julia set explorer.\n");

	master = 1;

	/* Get inputs-- region to view (must be within x/ymin to x/ymax, make sure

	xmax>xmin and ymax>ymin) and resolution (number of pixels along an edge,

	N x M, i.e. 256x256)

	*/

	read_mand_args(argc, argv, &imax_iterations, &ipixels_across, &ipixels_down,

	    &x_min, &x_max, &y_min, &y_max, &julia, &julia_constant.real,

	    &julia_constant.imaginary, &divergent_limit,

	    &alternate_equation, filename, &num_colors, &use_stdin, &save_image,

	    &num_children);

	check_mand_params(&imax_iterations, &ipixels_across, &ipixels_down,

	    &x_min, &x_max, &y_min, &y_max, &divergent_limit, &num_children);

	if (julia == 1) /* we're doing a julia figure */

	    check_julia_params(&julia_constant.real, &julia_constant.imaginary);

	/* spawn slaves */

	child_comm = (MPI_Comm*)malloc(num_children * sizeof(MPI_Comm));

	child_request = (MPI_Request*)malloc(num_children * sizeof(MPI_Request));

	result = (header_t*)malloc(num_children * sizeof(header_t));

	if (child_comm == NULL || child_request == NULL || result == NULL)

	{

	    printf("Error: unable to allocate an array of %d communicators, requests and work objects for the slaves.\n", num_children);

	    error_code = MPI_Abort(MPI_COMM_WORLD, -1);

	    exit(error_code);

	}

	printf("Spawning %d slaves.\n", num_children);

	for (i=0; i

	{

	    error = MPI_Comm_spawn(argv[0], MPI_ARGV_NULL, 1,

		MPI_INFO_NULL, 0, MPI_COMM_WORLD, &child_comm[i], &error_code);

	    if (error != MPI_SUCCESS)

	    {

		error_str[0] = '\0';

		length = MPI_MAX_ERROR_STRING;

		MPI_Error_string(error, error_str, &length);

		printf("Error: MPI_Comm_spawn failed: %s\n", error_str);

		error_code = MPI_Abort(MPI_COMM_WORLD, -1);

		exit(error_code);

	    }

	}

	/* send out parameters */

	for (i=0; i

	{

	    MPI_Send(&num_colors, 1, MPI_INT, 0, 0, child_comm[i]);

	    MPI_Send(&imax_iterations, 1, MPI_INT, 0, 0, child_comm[i]);

	    MPI_Send(&ipixels_across, 1, MPI_INT, 0, 0, child_comm[i]);

	    MPI_Send(&ipixels_down, 1, MPI_INT, 0, 0, child_comm[i]);

	    MPI_Send(&divergent_limit, 1, MPI_DOUBLE, 0, 0, child_comm[i]);

	    MPI_Send(&julia, 1, MPI_INT, 0, 0, child_comm[i]);

	    MPI_Send(&julia_constant.real, 1, MPI_DOUBLE, 0, 0, child_comm[i]);

	    MPI_Send(&julia_constant.imaginary, 1, MPI_DOUBLE, 0, 0, child_comm[i]);

	    MPI_Send(&alternate_equation, 1, MPI_INT, 0, 0, child_comm[i]);

	}

    }

    else

    {

	master = 0;

	MPI_Recv(&num_colors, 1, MPI_INT, 0, 0, parent, MPI_STATUS_IGNORE);

	MPI_Recv(&imax_iterations, 1, MPI_INT, 0, 0, parent, MPI_STATUS_IGNORE);

	MPI_Recv(&ipixels_across, 1, MPI_INT, 0, 0, parent, MPI_STATUS_IGNORE);

	MPI_Recv(&ipixels_down, 1, MPI_INT, 0, 0, parent, MPI_STATUS_IGNORE);

	MPI_Recv(&divergent_limit, 1, MPI_DOUBLE, 0, 0, parent, MPI_STATUS_IGNORE);

	MPI_Recv(&julia, 1, MPI_INT, 0, 0, parent, MPI_STATUS_IGNORE);

	MPI_Recv(&julia_constant.real, 1, MPI_DOUBLE, 0, 0, parent, MPI_STATUS_IGNORE);

	MPI_Recv(&julia_constant.imaginary, 1, MPI_DOUBLE, 0, 0, parent, MPI_STATUS_IGNORE);

	MPI_Recv(&alternate_equation, 1, MPI_INT, 0, 0, parent, MPI_STATUS_IGNORE);

    }

    if (master)

    {

	colors = malloc((num_colors+1)* sizeof(color_t));

	if (colors == NULL)

	{

	    MPI_Abort(MPI_COMM_WORLD, -1);

	    exit(-1);

	}

	Make_color_array(num_colors, colors);

	colors[num_colors] = 0; /* add one on the top to avoid edge errors */

    }

    /* allocate memory */

    if ( (in_grid_array = (int *)calloc(ipixels_across * ipixels_down, sizeof(int))) == NULL)

    {

	printf("Memory allocation failed for data array, aborting.\n");

	MPI_Abort(MPI_COMM_WORLD, -1);

	exit(-1);

    }

    if (master)

    {

	int istep, jstep;

	int i1[400], i2[400], j1[400], j2[400];

	int ii, jj;

	struct sockaddr_in addr;

	int len;

	char line[1024], *token;

	if ( (out_grid_array = (int *)calloc(ipixels_across * ipixels_down, sizeof(int))) == NULL)

	{

	    printf("Memory allocation failed for data array, aborting.\n");

	    MPI_Abort(MPI_COMM_WORLD, -1);

	    exit(-1);

	}

	srand(getpid());

	if (!use_stdin)

	{

	    addr.sin_family = AF_INET;

	    addr.sin_addr.s_addr = INADDR_ANY;

	    addr.sin_port = htons(DEFAULT_PORT);

	    listener = socket(AF_INET, SOCK_STREAM, 0);

	    if (listener == -1)

	    {

		printf("unable to create a listener socket.\n");

		MPI_Abort(MPI_COMM_WORLD, -1);

		exit(-1);

	    }

	    if (bind(listener, &addr, sizeof(addr)) == -1)

	    {

		addr.sin_port = 0;

		if (bind(listener, &addr, sizeof(addr)) == -1)

		{

		    printf("unable to create a listener socket.\n");

		    MPI_Abort(MPI_COMM_WORLD, -1);

		    exit(-1);

		}

	    }

	    if (listen(listener, 1) == -1)

	    {

		printf("unable to listen.\n");

		MPI_Abort(MPI_COMM_WORLD, -1);

		exit(-1);

	    }

	    len = sizeof(addr);

	    getsockname(listener, &addr, &len;;

	    printf("%s listening on port %d\n", processor_name, ntohs(addr.sin_port));

	    fflush(stdout);

	    sock = accept(listener, NULL, NULL);

	    if (sock == -1)

	    {

		printf("unable to accept a socket connection.\n");

		MPI_Abort(MPI_COMM_WORLD, -1);

		exit(-1);

	    }

	    printf("accepted connection from visualization program.\n");

	    fflush(stdout);

#ifdef HAVE_WINDOWS_H

	    optval = 1;

	    setsockopt(sock, IPPROTO_TCP, TCP_NODELAY, (char *)&optval, sizeof(optval));

#endif

	    printf("sending image size to visualizer.\n");

	    sock_write(sock, &ipixels_across, sizeof(int));

	    sock_write(sock, &ipixels_down, sizeof(int));

	    sock_write(sock, &num_colors, sizeof(int));

	    sock_write(sock, &imax_iterations, sizeof(int));

	}

	for (;;)

	{

	    /* get x_min, x_max, y_min, and y_max */

	    if (use_stdin)

	    {

		printf("input xmin ymin xmax ymax max_iter, (0 0 0 0 0 to quit):\n");fflush(stdout);

		fgets(line, 1024, stdin);

		printf("read <%s> from stdin\n", line);fflush(stdout);

		token = strtok(line, " \n");

		x_min = atof(token);

		token = strtok(NULL, " \n");

		y_min = atof(token);

		token = strtok(NULL, " \n");

		x_max = atof(token);

		token = strtok(NULL, " \n");

		y_max = atof(token);

		token = strtok(NULL, " \n");

		imax_iterations = atoi(token);

		/*sscanf(line, "%g %g %g %g", &x_min, &y_min, &x_max, &y_max);*/

		/*scanf("%g %g %g %g", &x_min, &y_min, &x_max, &y_max);*/

	    }

	    else

	    {

		printf("reading xmin,ymin,xmax,ymax.\n");fflush(stdout);

		sock_read(sock, &x_min, sizeof(double));

		sock_read(sock, &y_min, sizeof(double));

		sock_read(sock, &x_max, sizeof(double));

		sock_read(sock, &y_max, sizeof(double));

		sock_read(sock, &imax_iterations, sizeof(int));

	    }

	    printf("x0,y0 = (%f, %f) x1,y1 = (%f,%f) max_iter = %d\n", x_min, y_min, x_max, y_max, imax_iterations);fflush(stdout);

	    /*printf("sending the limits: (%f,%f)(%f,%f)\n", x_min, y_min, x_max, y_max);fflush(stdout);*/

	    /* let everyone know the limits */

	    for (i=0; i

	    {

		MPI_Send(&x_min, 1, MPI_DOUBLE, 0, 0, child_comm[i]);

		MPI_Send(&x_max, 1, MPI_DOUBLE, 0, 0, child_comm[i]);

		MPI_Send(&y_min, 1, MPI_DOUBLE, 0, 0, child_comm[i]);

		MPI_Send(&y_max, 1, MPI_DOUBLE, 0, 0, child_comm[i]);

		MPI_Send(&imax_iterations, 1, MPI_INT, 0, 0, child_comm[i]);

	    }

	    /* check for the end condition */

	    if (x_min == x_max && y_min == y_max)

	    {

		/*printf("root bailing.\n");fflush(stdout);*/

		break;

	    }

	    /* break the work up into 400 pieces */

	    istep = ipixels_across / 20;

	    jstep = ipixels_down / 20;

	    if (istep < 1)

		istep = 1;

	    if (jstep < 1)

		jstep = 1;

	    k = 0;

	    for (i=0; i<20; i++)

	    {

		for (j=0; j<20; j++)

		{

		    i1[k] = MIN(istep * i, ipixels_across - 1);

		    i2[k] = MIN((istep * (i+1)) - 1, ipixels_across - 1);

		    j1[k] = MIN(jstep * j, ipixels_down - 1);

		    j2[k] = MIN((jstep * (j+1)) - 1, ipixels_down - 1);

		    k++;

		}

	    }

	    /* shuffle the work */

	    for (i=0; i<500; i++)

	    {

		ii = rand() % 400;

		jj = rand() % 400;

		swap(&i1[ii], &i1[jj]);

		swap(&i2[ii], &i2[jj]);

		swap(&j1[ii], &j1[jj]);

		swap(&j2[ii], &j2[jj]);

	    }

	    /* send a piece of work to each worker (there must be more work than workers) */

	    k = 0;

	    for (i=0; i

	    {

		work[0] = k+1;

		work[1] = i1[k]; /* imin */

		work[2] = i2[k]; /* imax */

		work[3] = j1[k]; /* jmin */

		work[4] = j2[k]; /* jmax */

		/*printf("sending work(%d) to %d\n", k+1, i);fflush(stdout);*/

		MPI_Send(work, 5, MPI_INT, 0, 100, child_comm[i]);

		MPI_Irecv(&result[i], 5, MPI_INT, 0, 200, child_comm[i], &child_request[i]);

		k++;

	    }

	    /* receive the results and hand out more work until the image is complete */

	    while (k<400)

	    {

		MPI_Waitany(num_children, child_request, &index, &status);

		memcpy(work, &result[index], 5 * sizeof(int));

		/*printf("master receiving data in k<400 loop.\n");fflush(stdout);*/

		MPI_Recv(in_grid_array, (work[2] + 1 - work[1]) * (work[4] + 1 - work[3]), MPI_INT, 0, 201, child_comm[index], &status);

		/* draw data */

		output_data(in_grid_array, &work[1], out_grid_array, ipixels_across, ipixels_down);

		work[0] = k+1;

		work[1] = i1[k];

		work[2] = i2[k];

		work[3] = j1[k];

		work[4] = j2[k];

		/*printf("sending work(%d) to %d\n", k+1, index);fflush(stdout);*/

		MPI_Send(work, 5, MPI_INT, 0, 100, child_comm[index]);

		MPI_Irecv(&result[index], 5, MPI_INT, 0, 200, child_comm[index], &child_request[index]);

		k++;

	    }

	    /* receive the last pieces of work */

	    /* and tell everyone to stop */

	    for (i=0; i

	    {

		MPI_Wait(&child_request[i], &status);

		memcpy(work, &result[i], 5 * sizeof(int));

		/*printf("master receiving data in tail loop.\n");fflush(stdout);*/

		MPI_Recv(in_grid_array, (work[2] + 1 - work[1]) * (work[4] + 1 - work[3]), MPI_INT, 0, 201, child_comm[i], &status);

		/* draw data */

		output_data(in_grid_array, &work[1], out_grid_array, ipixels_across, ipixels_down);

		work[0] = 0;

		work[1] = 0;

		work[2] = 0;

		work[3] = 0;

		work[4] = 0;

		/*printf("sending %d to tell %d to stop\n", work[0], i);fflush(stdout);*/

		MPI_Send(work, 5, MPI_INT, 0, 100, child_comm[i]);

	    }

	    /* tell the visualizer the image is done */

	    if (!use_stdin)

	    {

		work[0] = 0;

		work[1] = 0;

		work[2] = 0;

		work[3] = 0;

		sock_write(sock, work, 4 * sizeof(int));

	    }

	}

    }

    else

    {

	for (;;)

	{

	    /*printf("slave[%d] receiveing bounds.\n", pid);fflush(stdout);*/

	    MPI_Recv(&x_min, 1, MPI_DOUBLE, 0, 0, parent, MPI_STATUS_IGNORE);

	    MPI_Recv(&x_max, 1, MPI_DOUBLE, 0, 0, parent, MPI_STATUS_IGNORE);

	    MPI_Recv(&y_min, 1, MPI_DOUBLE, 0, 0, parent, MPI_STATUS_IGNORE);

	    MPI_Recv(&y_max, 1, MPI_DOUBLE, 0, 0, parent, MPI_STATUS_IGNORE);

	    MPI_Recv(&imax_iterations, 1, MPI_INT, 0, 0, parent, MPI_STATUS_IGNORE);

	    /*printf("slave[%d] received bounding box: (%f,%f)(%f,%f)\n", pid, x_min, y_min, x_max, y_max);fflush(stdout);*/

	    /* check for the end condition */

	    if (x_min == x_max && y_min == y_max)

	    {

		/*printf("slave[%d] done.\n", pid);fflush(stdout);*/

		break;

	    }

	    x_resolution = (x_max-x_min)/ ((double)ipixels_across);

	    y_resolution = (y_max-y_min)/ ((double)ipixels_down);

	    MPI_Recv(work, 5, MPI_INT, 0, 100, parent, &status);

	    /*printf("slave[%d] received work: %d, (%d,%d)(%d,%d)\n", pid, work[0], work[1], work[2], work[3], work[4]);fflush(stdout);*/

	    while (work[0] != 0)

	    {

		imin = work[1];

		imax = work[2];

		jmin = work[3];

		jmax = work[4];

		k = 0;

		for (j=jmin; j<=jmax; ++j)

		{

		    coord_point.imaginary = y_max - j*y_resolution; /* go top to bottom */

		    for (i=imin; i<=imax; ++i)

		    {

			/* Call Mandelbrot routine for each code, fill array with number of iterations. */

			coord_point.real = x_min + i*x_resolution; /* go left to right */

			if (julia == 1)

			{

			    /* doing Julia set */

			    /* julia eq:  z = z^2 + c, z_0 = grid coordinate, c = constant */

			    icount = single_mandelbrot_point(coord_point, julia_constant, imax_iterations, divergent_limit);

			}

			else if (alternate_equation == 1)

			{

			    /* doing experimental form 1 */

			    icount = subtractive_mandelbrot_point(coord_point, julia_constant, imax_iterations, divergent_limit);

			}

			else if (alternate_equation == 2)

			{

			    /* doing experimental form 2 */

			    icount = additive_mandelbrot_point(coord_point, julia_constant, imax_iterations, divergent_limit);

			}

			else

			{

			    /* default to doing Mandelbrot set */

			    /* mandelbrot eq: z = z^2 + c, z_0 = c, c = grid coordinate */

			    icount = single_mandelbrot_point(coord_point, coord_point, imax_iterations, divergent_limit);

			}

			in_grid_array[k] = icount;

			++k;

		    }

		}

		/* send the result to the root */

		/*printf("slave[%d] sending work %d back.\n", pid, work[0]);fflush(stdout);*/

		MPI_Send(work, 5, MPI_INT, 0, 200, parent);

		/*printf("slave[%d] sending work %d data.\n", pid, work[0]);fflush(stdout);*/

		MPI_Send(in_grid_array, (work[2] + 1 - work[1]) * (work[4] + 1 - work[3]), MPI_INT, 0, 201, parent);

		/* get the next piece of work */

		/*printf("slave[%d] receiving new work.\n", pid);fflush(stdout);*/

		MPI_Recv(work, 5, MPI_INT, 0, 100, parent, &status);

		/*printf("slave[%d] received work: %d, (%d,%d)(%d,%d)\n", pid, work[0], work[1], work[2], work[3], work[4]);fflush(stdout);*/

	    }

	}

    }

    if (master && save_image)

    {

	imax_iterations = 0;

	for (i=0; i

	{

	    /* look for "brightest" pixel value, for image use */

	    if (out_grid_array[i] > imax_iterations)

		imax_iterations = out_grid_array[i];

	}

	if (julia == 0)

	    printf("Done calculating mandelbrot, now creating file\n");

	else

	    printf("Done calculating julia, now creating file\n");

	fflush(stdout);

	/* Print out the array in some appropriate form. */

	if (julia == 0)

	{

	    /* it's a mandelbrot */

	    sprintf(file_message, "Mandelbrot over (%lf-%lf,%lf-%lf), size %d x %d",

		x_min, x_max, y_min, y_max, ipixels_across, ipixels_down);

	}

	else

	{

	    /* it's a julia */

	    sprintf(file_message, "Julia over (%lf-%lf,%lf-%lf), size %d x %d, center (%lf, %lf)",

		x_min, x_max, y_min, y_max, ipixels_across, ipixels_down,

		julia_constant.real, julia_constant.imaginary);

	}

	dumpimage(filename, out_grid_array, ipixels_across, ipixels_down, imax_iterations, file_message, num_colors, colors);

    }

    if (master)

    {

	for (i=0; i

	{

	    MPI_Comm_disconnect(&child_comm[i]);

	}

	free(child_comm);

	free(child_request);

	free(colors);

    }

    MPI_Finalize();

    return 0;

}

void PrintUsage(void)

{

    printf("usage: mpiexec -n 1 pmandel [options]\n");

    printf("options:\n -n # slaves\n -xmin # -xmax #\n -ymin # -ymax #\n -depth #\n -xscale # -yscale #\n -out filename\n -i\n");

    printf("All options are optional.\n");

    printf("-i will allow you to input the min/max parameters from stdin and output the resulting image to a ppm file.");

    printf("  Otherwise the root process will listen for a separate visualizer program to connect to it.\n");

    printf("The defaults are:\n xmin %f, xmax %f\n ymin %f, ymax %f\n depth %d\n xscale %d, yscale %d\n %s\n",

	IDEAL_MIN_X_VALUE, IDEAL_MAX_X_VALUE, IDEAL_MIN_Y_VALUE, IDEAL_MAX_Y_VALUE, IDEAL_ITERATIONS,

	IDEAL_WIDTH, IDEAL_HEIGHT, default_filename);

    fflush(stdout);

}

int sock_write(int sock, void *buffer, int length)

{

    int result;

    int num_bytes;

    /*char err[100];*/

    int total = 0;

    struct timeval t;

    fd_set set;

    while (length)

    {

	num_bytes = send(sock, buffer, length, 0);

	if (num_bytes == -1)

	{

#ifdef HAVE_WINDOWS_H

	    result = WSAGetLastError();

	    if (result == WSAEWOULDBLOCK)

	    {

		FD_ZERO(&set);

		FD_SET(sock, &set);

		t.tv_sec = 1;

		t.tv_usec = 0;

		select(1, &set, NULL, NULL, &t);

		continue;

	    }

#else

	    if (errno == EWOULDBLOCK)

	    {

		FD_ZERO(&set);

		FD_SET(sock, &set);

		t.tv_sec = 1;

		t.tv_usec = 0;

		select(1, &set, NULL, NULL, &t);

		continue;

	    }

#endif

	    return total;

	}

	length -= num_bytes;

	buffer = (char*)buffer + num_bytes;

	total += num_bytes;

    }

    return total;

    /*return send(sock, buffer, length, 0);*/

}

int sock_read(int sock, void *buffer, int length)

{

    int result;

    int num_bytes;

    /*char err[100];*/

    int total = 0;

    struct timeval t;

    fd_set set;

    while (length)

    {

	num_bytes = recv(sock, buffer, length, 0);

	if (num_bytes == -1)

	{

#ifdef HAVE_WINDOWS_H

	    result = WSAGetLastError();

	    if (result == WSAEWOULDBLOCK)

	    {

		FD_ZERO(&set);

		FD_SET(sock, &set);

		t.tv_sec = 1;

		t.tv_usec = 0;

		select(1, &set, NULL, NULL, &t);

		continue;

	    }

#else

	    if (errno == EWOULDBLOCK)

	    {

		FD_ZERO(&set);

		FD_SET(sock, &set);

		t.tv_sec = 1;

		t.tv_usec = 0;

		select(1, &set, NULL, NULL, &t);

		continue;

	    }

#endif

	    return total;

	}

	length -= num_bytes;

	buffer = (char*)buffer + num_bytes;

	total += num_bytes;

    }

    return total;

    /*return recv(sock, buffer, length, 0);*/

}

color_t getColor(double fraction, double intensity)

{

    /* fraction is a part of the rainbow (0.0 - 1.0) = (Red-Yellow-Green-Cyan-Blue-Magenta-Red)

    intensity (0.0 - 1.0) 0 = black, 1 = full color, 2 = white

    */

    double red, green, blue;

    int r,g,b;

    double dtemp;

    fraction = fabs(modf(fraction, &dtemp));

    if (intensity > 2.0)

	intensity = 2.0;

    if (intensity < 0.0)

	intensity = 0.0;

    dtemp = 1.0/6.0;

    if (fraction < 1.0/6.0)

    {

	red = 1.0;

	green = fraction / dtemp;

	blue = 0.0;

    }

    else

    {

	if (fraction < 1.0/3.0)

	{

	    red = 1.0 - ((fraction - dtemp) / dtemp);

	    green = 1.0;

	    blue = 0.0;

	}

	else

	{

	    if (fraction < 0.5)

	    {

		red = 0.0;

		green = 1.0;

		blue = (fraction - (dtemp*2.0)) / dtemp;

	    }

	    else

	    {

		if (fraction < 2.0/3.0)

		{

		    red = 0.0;

		    green = 1.0 - ((fraction - (dtemp*3.0)) / dtemp);

		    blue = 1.0;

		}

		else

		{

		    if (fraction < 5.0/6.0)

		    {

			red = (fraction - (dtemp*4.0)) / dtemp;

			green = 0.0;

			blue = 1.0;

		    }

		    else

		    {

			red = 1.0;

			green = 0.0;

			blue = 1.0 - ((fraction - (dtemp*5.0)) / dtemp);

		    }

		}

	    }

	}

    }

    if (intensity > 1)

    {

	intensity = intensity - 1.0;

	red = red + ((1.0 - red) * intensity);

	green = green + ((1.0 - green) * intensity);

	blue = blue + ((1.0 - blue) * intensity);

    }

    else

    {

	red = red * intensity;

	green = green * intensity;

	blue = blue * intensity;

    }

    r = (int)(red * 255.0);

    g = (int)(green * 255.0);

    b = (int)(blue * 255.0);

    return RGBtocolor_t(r,g,b);

}

int Make_color_array(int num_colors, color_t colors[])

{

    double fraction, intensity;

    int i;

    intensity = 1.0;

    for (i=0; i

    {

	fraction = (double)i / (double)num_colors;

	colors[i] = getColor(fraction, intensity);

    }

    return 0;

}

void getRGB(color_t color, int *r, int *g, int *b)

{

    *r = getR(color);

    *g = getG(color);

    *b = getB(color);

}

void read_mand_args(int argc, char *argv[], int *o_max_iterations,

		    int *o_pixels_across, int *o_pixels_down,

		    double *o_x_min, double *o_x_max,

		    double *o_y_min, double *o_y_max, int *o_julia,

		    double *o_julia_real_x, double *o_julia_imaginary_y,

		    double *o_divergent_limit, int *o_alternate,

		    char *filename, int *o_num_colors, int *use_stdin,

		    int *save_image, int *num_children)

{	

    int i;

    *o_julia_real_x = NOVALUE;

    *o_julia_imaginary_y = NOVALUE;

    /* set defaults */

    *o_max_iterations = IDEAL_ITERATIONS;