[MPICH] VIsualization ??? FYI - mandel viewer - Just use MPI_Send/MPI_Recv
Jayesh Krishna
jayesh at mcs.anl.gov
Tue Apr 24 15:29:24 CDT 2007
Hi,
Once you run your programs as shown below (using mpiexec and separating
different executables using ":") you can use any MPI function
(MPI_Send/MPI_Recv/MPI_Bcast etc) to exchange information across the
processes. You do not need to use MPI_Open_port()/MPI_Comm_accept() etc to
communicate across the MPI processes (the gui and process).
Regards,
Jayesh
-----Original Message-----
From: Ralph Pass [mailto:rppass at rppass.com]
Sent: Monday, April 23, 2007 3:11 PM
To: Jayesh Krishna
Subject: Re: [MPICH] VIsualization ??? FYI - mandel viewer
Jayesh:
I have looked at your code and it looks very similar. I did add the
broadcast of the port_name from the root process to the other processes
and my visualization program could see it. However, the MPI_Comm_accept
and MPI_Comm_connect pair do not seem to want to talk, even if I start
both programs with something like:
mpiexec -n 4 process : -localonly 1 gui
All the processes report 5 processes, and gui reports its number as 4.
I am about ready to try this with a Unix/Linux box to see if there is
something different (than the behavior on my Windows machines).
Thanks,
Ralph
Jayesh Krishna wrote:
> Hi,
> I have attached the code for MandelViewer, which is a program similar to
> your requirements. The program opens a port and broadcasts this port
number
> to all visualizers in MPI_COMM_WORLD. The visualizers connect to this port
> and the program sends the appropriate information to the visualizer.
>
>
> Regards,
> Jayesh
>
> -----Original Message-----
> From: rppass at rppass.com [mailto:rppass at rppass.com]
> Sent: Tuesday, April 17, 2007 10:54 AM
> To: Jayesh Krishna
> Cc: rppass at rppass.com; 'Trach-Minh Tran'; mpich-discuss at mcs.anl.gov
> Subject: RE: [MPICH] VIsualization ???
>
> Thanks you for your efforts.
>
> SO I take it that if I start the two jobs (computation and visualization)
> with the same mpiexec command, then I can publish names?
>
> To try to work around this, I start the computation and have it print out
> its port-name (looking like tag=0 port..........) I then enter that
> port-name into the visualization program and try to connect to that
explicit
> port name. This does not connect. So I have two problems! The first
being
> the name lookup and then the second, if I get the port-name from the name
> lookup, using it to connect (unless the format or structure of the
port-name
> from the name lookup is different than that returned by MPI_OPEN_PORT
> command.
>
> Ralph
>
>
>
>> Hi,
>> MPICH2 on windows only supports the smpd process manager. I am
>> currently looking into the code to find out the extend of support
>> available in smpd process manager for publishing a name (You can
>> publish names across MPI processes run using the same mpiexec command
>> using smpd. MPD allows you to specify the machines in the MPD ring
>> during mpdboot time and the published name is published to all the
>>
> machines in the MPD ring.).
>
>> (Note: There is a -sethosts option in smpd which could be of help. I
>> will dig into the code and let you know soon.)
>>
>> Regards,
>> Jayesh
>>
>> -----Original Message-----
>> From: owner-mpich-discuss at mcs.anl.gov
>> [mailto:owner-mpich-discuss at mcs.anl.gov] On Behalf Of
>> rppass at rppass.com
>> Sent: Tuesday, April 17, 2007 7:25 AM
>> To: Trach-Minh Tran
>> Cc: Ralph Pass; mpich-discuss at mcs.anl.gov
>> Subject: Re: [MPICH] VIsualization ???
>>
>> Naive me, I just got MPICH2 from the MPICH2 web site at Argonne.
>>
>> I am running on Windows XP.
>>
>> How do I tell the version of MPICH2 that I got?
>>
>> Thank you,
>> Ralph Pass
>>
>>
>>
>>> On 04/13/2007 09:26 PM, Ralph Pass wrote:
>>>
>>>> I am trying to get a separate process (visualizer) to connect to a
>>>> server for computation of results to be visualized). I am tyring to
>>>> emulate the code in Using MPI-2 (I downloaded the version yesterday).
>>>>
>>>> I try to publish the name of the computation program and then accept
>>>> a comm.
>>>>
>>>> On the visualizer side I try to lookup the name.
>>>>
>>>> I am using C++ on a Windows XP system.
>>>>
>>>> I run the computation program (using mpiexec to start it) and it is
>>>> fat and dumb, awaiting a connection and then a data transfer.
>>>> On the visualization side, when I try the lookup call [looks like:
>>>> MPI_Lookup_name("ImageProcessing", MPI_INFO_NULL, mMPIPortName);]
>>>>
>>>> I get an error return: 536999969. When I convert this to a string,
>>>> I
>>>> get:
>>>>
>>>> Error message Invalid service name (see MPI_Publish_name), error stack:
>>>> MPID_NS_Lookup(70): Lookup failed for service name ImageProcessing
>>>>
>>>>
>>>>
>>>> I start the computation with mpiexec -n 4 processing
>>>>
>>>>
>>>> and the visualizer with
>>>>
>>>> mpiexec -localonly 1 visualizer
>>>>
>>>> Any help or comments? (Like you need to start the two programs in
>>>> the same mpiexec command)
>>>>
>>>>
>>> Which mpiexec did you use? The publish_name/lookup_name does not work
>>> with the OSC's mpiexec. However with the MPICH2 MPD version, it does.
>>> After setting up the MPD with "mpdboot ...", I just run
>>>
>>> mpiexec -n 1 ./viz < /dev/null &
>>> mpiexec -n 4 ./comp
>>>
>>> Hope this help.
>>>
>>> -Minh.
>>>
>>>
>>>
>>
>>
>>
>
>
>
> ------------------------------------------------------------------------
>
> /* -*- Mode: C; c-basic-offset:4 ; -*- */
> /*
> * (C) 2001 by Argonne National Laboratory.
> * See COPYRIGHT in top-level directory.
> */
> #include <stdio.h>
> #include <stdlib.h>
> #include <math.h>
> #include <string.h>
> #include "mpi.h"
> #ifdef HAVE_WINDOWS_H
> #include <process.h>
> #endif
>
> /* 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 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);
> 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);
> 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);
>
> #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;
> MPI_Comm comm;
>
> void swap(int *i, int *j)
> {
> int x;
> x = *i;
> *i = *j;
> *j = x;
> }
>
> 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 temp[5];
> /* 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 save_image = 0;
> int result;
> char mpi_port[MPI_MAX_PORT_NAME];
> MPI_Info info;
>
> MPI_Init(&argc, &argv);
> MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
> MPI_Comm_rank(MPI_COMM_WORLD, &myid);
> MPI_Get_processor_name(processor_name, &namelen);
>
> if (numprocs == 1)
> {
> PrintUsage();
> MPI_Finalize();
> exit(0);
> }
>
> if (myid == 0)
> {
> printf("Welcome to the Mandelbrot/Julia set explorer.\n");
>
> /* 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);
> check_mand_params(&imax_iterations, &ipixels_across, &ipixels_down,
> &x_min, &x_max, &y_min, &y_max, &divergent_limit);
>
> if (julia == 1) /* we're doing a julia figure */
> check_julia_params(&julia_constant.real,
&julia_constant.imaginary);
>
> MPI_Bcast(&use_stdin, 1, MPI_INT, 0, MPI_COMM_WORLD);
> MPI_Bcast(&num_colors, 1, MPI_INT, 0, MPI_COMM_WORLD);
> MPI_Bcast(&imax_iterations, 1, MPI_INT, 0, MPI_COMM_WORLD);
> MPI_Bcast(&ipixels_across, 1, MPI_INT, 0, MPI_COMM_WORLD);
> MPI_Bcast(&ipixels_down, 1, MPI_INT, 0, MPI_COMM_WORLD);
> MPI_Bcast(&divergent_limit, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
> MPI_Bcast(&julia, 1, MPI_INT, 0, MPI_COMM_WORLD);
> MPI_Bcast(&julia_constant.real, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
> MPI_Bcast(&julia_constant.imaginary, 1, MPI_DOUBLE, 0,
MPI_COMM_WORLD);
> MPI_Bcast(&alternate_equation, 1, MPI_INT, 0, MPI_COMM_WORLD);
> }
> else
> {
> MPI_Bcast(&use_stdin, 1, MPI_INT, 0, MPI_COMM_WORLD);
> MPI_Bcast(&num_colors, 1, MPI_INT, 0, MPI_COMM_WORLD);
> MPI_Bcast(&imax_iterations, 1, MPI_INT, 0, MPI_COMM_WORLD);
> MPI_Bcast(&ipixels_across, 1, MPI_INT, 0, MPI_COMM_WORLD);
> MPI_Bcast(&ipixels_down, 1, MPI_INT, 0, MPI_COMM_WORLD);
> MPI_Bcast(&divergent_limit, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
> MPI_Bcast(&julia, 1, MPI_INT, 0, MPI_COMM_WORLD);
> MPI_Bcast(&julia_constant.real, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
> MPI_Bcast(&julia_constant.imaginary, 1, MPI_DOUBLE, 0,
MPI_COMM_WORLD);
> MPI_Bcast(&alternate_equation, 1, MPI_INT, 0, MPI_COMM_WORLD);
> }
>
> if (myid == 0)
> {
> 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 (myid == 0)
> {
> int istep, jstep, cur_proc;
> int i1[400], i2[400], j1[400], j2[400];
> int ii, jj;
> 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)
> {
> result = MPI_Info_create(&info);
> if (result != MPI_SUCCESS)
> {
> printf("Unable to create an MPI_Info to be used in
MPI_Open_port.\n");
> MPI_Abort(MPI_COMM_WORLD, -1);
> }
> result = MPI_Open_port(/*info*/MPI_INFO_NULL, mpi_port);
> if (result != MPI_SUCCESS)
> {
> printf("MPI_Open_port failed, aborting.\n");
> MPI_Abort(MPI_COMM_WORLD, -1);
> exit(-1);
> }
> result = MPI_Bcast(mpi_port, MPI_MAX_PORT_NAME, MPI_CHAR, 0,
MPI_COMM_WORLD);
> if (result != MPI_SUCCESS)
> {
> printf("Unable to broadcast the port name to
COMM_WORLD.\n");
> MPI_Abort(MPI_COMM_WORLD, -1);
> }
> printf("%s listening on port:\n%s\n", processor_name, mpi_port);
> fflush(stdout);
>
> result = MPI_Comm_accept(mpi_port, info, 0, MPI_COMM_WORLD,
&comm);
> if (result != MPI_SUCCESS)
> {
> printf("MPI_Comm_accept failed, aborting.\n");
> MPI_Abort(MPI_COMM_WORLD, -1);
> exit(-1);
> }
> printf("accepted connection from visualization program.\n");
> fflush(stdout);
>
> printf("sending image size to visualizer.\n");
> /*
> printf("sending pixels_across.\n");
> fflush(stdout);
> */
> result = MPI_Send(&ipixels_across, 1, MPI_INT, 0, 0, comm);
> if (result != MPI_SUCCESS)
> {
> printf("Unable to send pixel width to visualizer,
aborting.\n");
> MPI_Abort(MPI_COMM_WORLD, -1);
> }
> /*
> printf("sending pixels_down.\n");
> fflush(stdout);
> */
> result = MPI_Send(&ipixels_down, 1, MPI_INT, 0, 0, comm);
> if (result != MPI_SUCCESS)
> {
> printf("Unable to send pixel height to visualizer,
aborting.\n");
> MPI_Abort(MPI_COMM_WORLD, -1);
> }
> /*
> printf("sending num_colors.\n");
> fflush(stdout);
> */
> result = MPI_Send(&num_colors, 1, MPI_INT, 0, 0, comm);
> if (result != MPI_SUCCESS)
> {
> printf("Unable to send num_colors to visualizer,
aborting.\n");
> MPI_Abort(MPI_COMM_WORLD, -1);
> }
> result = MPI_Send(&imax_iterations, 1, MPI_INT, 0, 0, comm);
> if (result != MPI_SUCCESS)
> {
> printf("Unable to send imax_iterations to visualizer,
aborting.\n");
> MPI_Abort(MPI_COMM_WORLD, -1);
> }
> }
>
> 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);
> result = MPI_Recv(&x_min, 1, MPI_DOUBLE, 0, 0, comm,
&status);
> if (result != MPI_SUCCESS)
> {
> printf("Unable to receive x_min from the visualizer,
aborting.\n");
> MPI_Abort(MPI_COMM_WORLD, -1);
> }
> result = MPI_Recv(&y_min, 1, MPI_DOUBLE, 0, 0, comm,
&status);
> if (result != MPI_SUCCESS)
> {
> printf("Unable to receive y_min from the visualizer,
aborting.\n");
> MPI_Abort(MPI_COMM_WORLD, -1);
> }
> result = MPI_Recv(&x_max, 1, MPI_DOUBLE, 0, 0, comm,
&status);
> if (result != MPI_SUCCESS)
> {
> printf("Unable to receive x_max from the visualizer,
aborting.\n");
> MPI_Abort(MPI_COMM_WORLD, -1);
> }
> result = MPI_Recv(&y_max, 1, MPI_DOUBLE, 0, 0, comm,
&status);
> if (result != MPI_SUCCESS)
> {
> printf("Unable to receive y_max from the visualizer,
aborting.\n");
> MPI_Abort(MPI_COMM_WORLD, -1);
> }
> result = MPI_Recv(&imax_iterations, 1, MPI_INT, 0, 0, comm,
&status);
> if (result != MPI_SUCCESS)
> {
> printf("Unable to receive imax_iterations from the
visualizer, aborting.\n");
> MPI_Abort(MPI_COMM_WORLD, -1);
> }
> }
> 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);
>
> /* 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]);
> }
>
> /*printf("bcasting the limits: (%f,%f)(%f,%f)\n", x_min, y_min,
x_max, y_max);fflush(stdout);*/
> /* let everyone know the limits */
> MPI_Bcast(&x_min, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
> MPI_Bcast(&x_max, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
> MPI_Bcast(&y_min, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
> MPI_Bcast(&y_max, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
> MPI_Bcast(&imax_iterations, 1, MPI_INT, 0, MPI_COMM_WORLD);
>
> /* check for the end condition */
> if (x_min == x_max && y_min == y_max)
> {
> break;
> }
>
> /* send a piece of work to each worker (there must be more work
than workers) */
> cur_proc = 1;
> k = 0;
> for (i=1; i<numprocs; i++)
> {
> temp[0] = k+1;
> temp[1] = i1[k]; /* imin */
> temp[2] = i2[k]; /* imax */
> temp[3] = j1[k]; /* jmin */
> temp[4] = j2[k]; /* jmax */
>
> /*printf("sending work(%d) to %d\n", k+1,
cur_proc);fflush(stdout);*/
> MPI_Send(temp, 5, MPI_INT, cur_proc, 100, MPI_COMM_WORLD);
> cur_proc++;
> k++;
> }
> /* receive the results and hand out more work until the image is
complete */
> while (k<400)
> {
> MPI_Recv(temp, 5, MPI_INT, MPI_ANY_SOURCE, 200,
MPI_COMM_WORLD, &status);
> cur_proc = temp[0];
> MPI_Recv(in_grid_array, (temp[2] + 1 - temp[1]) * (temp[4] +
1 - temp[3]), MPI_INT, cur_proc, 201, MPI_COMM_WORLD, &status);
> /* draw data */
> output_data(in_grid_array, &temp[1], out_grid_array,
ipixels_across, ipixels_down);
> temp[0] = k+1;
> temp[1] = i1[k];
> temp[2] = i2[k];
> temp[3] = j1[k];
> temp[4] = j2[k];
> /*printf("sending work(%d) to %d\n", k+1,
cur_proc);fflush(stdout);*/
> MPI_Send(temp, 5, MPI_INT, cur_proc, 100, MPI_COMM_WORLD);
> k++;
> }
> /* receive the last pieces of work */
> /* and tell everyone to stop */
> for (i=1; i<numprocs; i++)
> {
> MPI_Recv(temp, 5, MPI_INT, MPI_ANY_SOURCE, 200,
MPI_COMM_WORLD, &status);
> cur_proc = temp[0];
> MPI_Recv(in_grid_array, (temp[2] + 1 - temp[1]) * (temp[4] +
1 - temp[3]), MPI_INT, cur_proc, 201, MPI_COMM_WORLD, &status);
> /* draw data */
> output_data(in_grid_array, &temp[1], out_grid_array,
ipixels_across, ipixels_down);
> temp[0] = 0;
> temp[1] = 0;
> temp[2] = 0;
> temp[3] = 0;
> temp[4] = 0;
> /*printf("sending %d to tell %d to stop\n", temp[0],
cur_proc);fflush(stdout);*/
> MPI_Send(temp, 5, MPI_INT, cur_proc, 100, MPI_COMM_WORLD);
> }
>
> /* tell the visualizer the image is done */
> if (!use_stdin)
> {
> temp[0] = 0;
> temp[1] = 0;
> temp[2] = 0;
> temp[3] = 0;
> result = MPI_Send(temp, 4, MPI_INT, 0, 0, comm);
> if (result != MPI_SUCCESS)
> {
> printf("Unable to send a done command to the visualizer,
aborting.\n");
> MPI_Abort(MPI_COMM_WORLD, -1);
> }
> }
> }
> }
> else
> {
> if (!use_stdin)
> {
> result = MPI_Info_create(&info);
> if (result != MPI_SUCCESS)
> {
> printf("Unable to create an MPI_Info to be used in
MPI_Open_port.\n");
> MPI_Abort(MPI_COMM_WORLD, -1);
> }
> result = MPI_Bcast(mpi_port, MPI_MAX_PORT_NAME, MPI_CHAR, 0,
MPI_COMM_WORLD);
> if (result != MPI_SUCCESS)
> {
> printf("Unable to broadcast the port name on
COMM_WORLD.\n");
> MPI_Abort(MPI_COMM_WORLD, -1);
> }
> /*
> printf("%s listening on port: %s\n", processor_name, mpi_port);
> fflush(stdout);
> */
>
> result = MPI_Comm_accept(mpi_port, info, 0, MPI_COMM_WORLD,
&comm);
> if (result != MPI_SUCCESS)
> {
> printf("MPI_Comm_accept failed, aborting.\n");
> MPI_Abort(MPI_COMM_WORLD, -1);
> exit(-1);
> }
> }
> for (;;)
> {
> MPI_Bcast(&x_min, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
> MPI_Bcast(&x_max, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
> MPI_Bcast(&y_min, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
> MPI_Bcast(&y_max, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
> MPI_Bcast(&imax_iterations, 1, MPI_INT, 0, MPI_COMM_WORLD);
>
> /* check for the end condition */
> if (x_min == x_max && y_min == y_max)
> {
> break;
> }
>
> x_resolution = (x_max-x_min)/ ((double)ipixels_across);
> y_resolution = (y_max-y_min)/ ((double)ipixels_down);
>
> MPI_Recv(temp, 5, MPI_INT, 0, 100, MPI_COMM_WORLD, &status);
> while (temp[0] != 0)
> {
> imin = temp[1];
> imax = temp[2];
> jmin = temp[3];
> jmax = temp[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 */
> temp[0] = myid;
> MPI_Send(temp, 5, MPI_INT, 0, 200, MPI_COMM_WORLD);
> MPI_Send(in_grid_array, (temp[2] + 1 - temp[1]) * (temp[4] +
1 - temp[3]), MPI_INT, 0, 201, MPI_COMM_WORLD);
> /* get the next piece of work */
> MPI_Recv(temp, 5, MPI_INT, 0, 100, MPI_COMM_WORLD, &status);
> }
> }
> }
>
> if (myid == 0 && save_image)
> {
> imax_iterations = 0;
> for (i=0; i<ipixels_across * ipixels_down; ++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);
> }
>
> MPI_Comm_disconnect(&comm);
>
> MPI_Finalize();
> free(colors);
> return 0;
> } /* end of main */
>
> void PrintUsage(void)
> {
> printf("usage: mpiexec -n x pmandel [options]\n");
> printf("options:\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);
> }
>
> 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<num_colors; 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 i;
>
> *o_julia_real_x = NOVALUE;
> *o_julia_imaginary_y = NOVALUE;
>
> /* set defaults */
> *o_max_iterations = IDEAL_ITERATIONS;
> *o_pixels_across = IDEAL_WIDTH;
> *o_pixels_down = IDEAL_HEIGHT;
> *o_x_min = IDEAL_MIN_X_VALUE;
> *o_x_max = IDEAL_MAX_X_VALUE;
> *o_y_min = IDEAL_MIN_Y_VALUE;
> *o_y_max = IDEAL_MAX_Y_VALUE;
> *o_divergent_limit = INFINITE_LIMIT;
> strcpy(filename, default_filename);
> *o_num_colors = IDEAL_ITERATIONS;
> *use_stdin = 0; /* default is to listen for a controller */
> *save_image = NOVALUE;
>
> *o_julia = 0; /* default is "generate Mandelbrot" */
> *o_alternate = 0; /* default is still "generate Mandelbrot" */
> *o_divergent_limit = INFINITE_LIMIT; /* default total range is assumed
> if not explicitly overwritten */
>
> /* We just cycle through all given parameters, matching what we can.
> Note that we force casting, because we expect that a nonsensical
> value is better than a poorly formatted one (fewer crashes), and
> we'll later sort out the good from the bad
> */
>
> for (i = 0; i < argc; ++i) /* grab command line arguments */
> if (strcmp(argv[i], "-xmin\0") == 0 && argv[i+1] != NULL)
> sscanf(argv[i+1], "%lf", &*o_x_min);
> else if (strcmp(argv[i], "-xmax\0") == 0 && argv[i+1] != NULL)
> sscanf(argv[i+1], "%lf", &*o_x_max);
> else if (strcmp(argv[i], "-ymin\0") == 0 && argv[i+1] != NULL)
> sscanf(argv[i+1], "%lf", &*o_y_min);
> else if (strcmp(argv[i], "-ymax\0") == 0 && argv[i+1] != NULL)
> sscanf(argv[i+1], "%lf", &*o_y_max);
> else if (strcmp(argv[i], "-depth\0") == 0 && argv[i+1] != NULL)
> sscanf(argv[i+1], "%d", &*o_max_iterations);
> else if (strcmp(argv[i], "-xscale\0") == 0 && argv[i+1] != NULL)
> sscanf(argv[i+1], "%d", &*o_pixels_across);
> else if (strcmp(argv[i], "-yscale\0") == 0 && argv[i+1] != NULL)
> sscanf(argv[i+1], "%d", &*o_pixels_down);
> else if (strcmp(argv[i], "-diverge\0") == 0 && argv[i+1] != NULL)
> sscanf(argv[i+1], "%lf", &*o_divergent_limit);
> else if (strcmp(argv[i], "-jx\0") == 0 && argv[i+1] != NULL)
> sscanf(argv[i+1], "%lf", &*o_julia_real_x);
> else if (strcmp(argv[i], "-jy\0") == 0 && argv[i+1] != NULL)
> sscanf(argv[i+1], "%lf", &*o_julia_imaginary_y);
> else if (strcmp(argv[i], "-alternate\0") == 0 && argv[i+1] != NULL)
> sscanf(argv[i+1], "%d", &*o_alternate);
> else if (strcmp(argv[i], "-julia\0") == 0)
> *o_julia = 1;
> else if (strcmp(argv[i], "-out\0") == 0 && argv[i+1] != NULL)
> strcpy(filename, argv[i+1]);
> else if (strcmp(argv[i], "-colors\0") == 0 && argv[i+1] != NULL)
> sscanf(argv[i+1], "%d", &*o_num_colors);
> else if (strcmp(argv[i], "-i\0") == 0)
> *use_stdin = 1;
> else if (strcmp(argv[i], "-save\0") == 0)
> *save_image = 1;
> else if (strcmp(argv[i], "-nosave\0") == 0)
> *save_image = 0;
> else if (strcmp(argv[i], "-help\0") == 0)
> {
> PrintUsage();
> MPI_Finalize();
> exit(0);
> }
>
> if (*save_image == NOVALUE)
> {
> if (*use_stdin == 1)
> *save_image = 1;
> else
> *save_image = 0;
> }
> #if 0
> if (myid == 0)
> {
> printf("Doing %d iterations over (%.2lf:%.2lf,%.2lf:%.2lf) on a
%d x %d grid\n",
> *o_max_iterations, *o_x_min,*o_x_max,*o_y_min,*o_y_max,
> *o_pixels_across, *o_pixels_down);
> }
> #endif
> }
>
> 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)
> {
> /* Get first batch of limits */
> #if PROMPT
> if (*m_x_min == NOVALUE || *m_x_max == NOVALUE)
> {
> /* grab unspecified limits */
> printf("Enter lower and higher limits on x (within -1 to 1): ");
> scanf("%lf %lf", &*m_x_min, &*m_x_max);
> }
>
> if (*m_y_min == NOVALUE || *m_y_max == NOVALUE)
> {
> /* grab unspecified limits */
> printf("Enter lower and higher limits on y (within -1 to 1): ");
> scanf("%lf %lf", &*m_y_min, &*m_y_max);
> }
> #endif
>
> /* Verify limits are reasonable */
>
> if (*m_x_min < MIN_X_VALUE || *m_x_min > *m_x_max)
> {
> printf("Chosen lower x limit is too low, resetting to
%lf\n",MIN_X_VALUE);
> *m_x_min = MIN_X_VALUE;
> }
> if (*m_x_max > MAX_X_VALUE || *m_x_max <= *m_x_min)
> {
> printf("Chosen upper x limit is too high, resetting to
%lf\n",MAX_X_VALUE);
> *m_x_max = MAX_X_VALUE;
> }
> if (*m_y_min < MIN_Y_VALUE || *m_y_min > *m_y_max)
> {
> printf("Chosen lower y limit is too low, resetting to
%lf\n",MIN_Y_VALUE);
> *m_y_min = MIN_Y_VALUE;
> }
> if (*m_y_max > MAX_Y_VALUE || *m_y_max <= *m_y_min)
> {
> printf("Chosen upper y limit is too high, resetting to
%lf\n",MAX_Y_VALUE);
> *m_y_max = MAX_Y_VALUE;
> }
>
> /* Get second set of limits */
> #if PROMPT
>
> if (*m_max_iterations == NOVALUE)
> {
> /* grab unspecified limit */
> printf("Enter max interations to do: ");
> scanf("%d",&*m_max_iterations);
> }
> #endif
>
> /* Verify second set of limits */
>
> if (*m_max_iterations > MAX_ITERATIONS || *m_max_iterations < 0)
> {
> printf("Warning, unreasonable number of iterations, setting to
%d\n", MAX_ITERATIONS);
> *m_max_iterations = MAX_ITERATIONS;
> }
>
> /* Get third set of limits */
> #if PROMPT
> if (*m_pixels_across == NOVALUE || *m_pixels_down == NOVALUE)
> {
> /* grab unspecified limits */
> printf("Enter pixel size (horizonal by vertical, i.e. 256 256): ");
> scanf(" %d %d", &*m_pixels_across, &*m_pixels_down);
> }
> #endif
>
> /* Verify third set of limits */
>
> if (*m_pixels_across > MAX_WIDTH)
> {
> printf("Warning, image requested is too wide, setting to %d pixel
width\n", MAX_WIDTH);
> *m_pixels_across = MAX_WIDTH;
> }
> if (*m_pixels_down > MAX_HEIGHT)
> {
> printf("Warning, image requested is too tall, setting to %d pixel
height\n", MAX_HEIGHT);
> *m_pixels_down = MAX_HEIGHT;
> }
>
> #if 0
> printf("%d iterations over (%.2lf:%.2lf,%.2lf:%.2lf), %d x %d grid,
diverge value %lf\n",
> *m_max_iterations, *m_x_min,*m_x_max,*m_y_min,*m_y_max,
> *m_pixels_across, *m_pixels_down, *m_divergent_limit);
> #endif
> }
>
> void check_julia_params(double *julia_x_constant, double
*julia_y_constant)
> {
>
> /* Hey, we're not stupid-- if they must Prompt, we will also be Noisy
*/
> #if PROMPT
> if (*julia_x_constant == NOVALUE || *julia_y_constant == NOVALUE)
> {
> /* grab unspecified limits */
> printf("Enter Coordinates for Julia expansion: ");
> scanf("%lf %lf", &*julia_x_constant, &*julia_y_constant);
> }
> #endif
>
> #if 0
> /* In theory, any point can be investigated, although it's not much
point
> if it's outside of the area viewed. But, hey, that's not our problem
*/
> printf("Exploring julia set around (%lf, %lf)\n", *julia_x_constant,
*julia_y_constant);
> #endif
> }
>
> /* This is a Mandelbrot variant, solving the code for the equation:
> z = z(1-z)
> */
>
> /* This routine takes a complex coordinate point (x+iy) and a value
stating
> what the upper limit to the number of iterations is. It eventually
> returns an integer of how many counts the code iterated for within
> the given point/region until the exit condition ( abs(x+iy) > limit) was
met.
> This value is returned as an integer.
> */
>
> int subtractive_mandelbrot_point(complex_t coord_point,
> complex_t c_constant,
> int Max_iterations, double divergent_limit)
> {
> complex_t z_point, a_point; /* we need 2 pts to use in our calculation
*/
> int num_iterations; /* a counter to track the number of iterations
done */
>
> num_iterations = 0; /* zero our counter */
> z_point = coord_point; /* initialize to the given start coordinate */
>
> /* loop while the absolute value of the complex coordinate is < our
limit
> (for a mandelbrot) or until we've done our specified maximum number of
> iterations (both julia and mandelbrot) */
>
> while (absolute_complex(z_point) < divergent_limit &&
> num_iterations < Max_iterations)
> {
> /* z = z(1-z) */
> a_point.real = 1.0; a_point.imaginary = 0.0; /* make "1" */
> a_point = subtract_complex(a_point,z_point);
> z_point = multiply_complex(z_point,a_point);
>
> ++num_iterations;
> } /* done iterating for one point */
>
> return num_iterations;
> }
>
> /* This is a Mandelbrot variant, solving the code for the equation:
> z = z(z+c)
> */
>
> /* This routine takes a complex coordinate point (x+iy) and a value
stating
> what the upper limit to the number of iterations is. It eventually
> returns an integer of how many counts the code iterated for within
> the given point/region until the exit condition ( abs(x+iy) > limit) was
met.
> This value is returned as an integer.
> */
>
> int additive_mandelbrot_point(complex_t coord_point,
> complex_t c_constant,
> int Max_iterations, double divergent_limit)
> {
> complex_t z_point, a_point; /* we need 2 pts to use in our calculation
*/
> int num_iterations; /* a counter to track the number of iterations
done */
>
> num_iterations = 0; /* zero our counter */
> z_point = coord_point; /* initialize to the given start coordinate */
>
> /* loop while the absolute value of the complex coordinate is < our
limit
> (for a mandelbrot) or until we've done our specified maximum number of
> iterations (both julia and mandelbrot) */
>
> while (absolute_complex(z_point) < divergent_limit &&
> num_iterations < Max_iterations)
> {
> /* z = z(z+C) */
> a_point = add_complex(z_point,coord_point);
> z_point = multiply_complex(z_point,a_point);
>
> ++num_iterations;
> } /* done iterating for one point */
>
> return num_iterations;
> }
>
> /* This is a Mandelbrot variant, solving the code for the equation:
> z = z(z+1)
> */
>
> /* This routine takes a complex coordinate point (x+iy) and a value
stating
> what the upper limit to the number of iterations is. It eventually
> returns an integer of how many counts the code iterated for within
> the given point/region until the exit condition ( abs(x+iy) > limit) was
met.
> This value is returned as an integer.
> */
>
> int exponential_mandelbrot_point(complex_t coord_point,
> complex_t c_constant,
> int Max_iterations, double divergent_limit)
> {
> complex_t z_point, a_point; /* we need 2 pts to use in our calculation
*/
> int num_iterations; /* a counter to track the number of iterations
done */
>
> num_iterations = 0; /* zero our counter */
> z_point = coord_point; /* initialize to the given start coordinate */
>
> /* loop while the absolute value of the complex coordinate is < our
limit
> (for a mandelbrot) or until we've done our specified maximum number of
> iterations (both julia and mandelbrot) */
>
> while (absolute_complex(z_point) < divergent_limit &&
> num_iterations < Max_iterations)
> {
> /* z = z(1-z) */
> a_point.real = 1.0; a_point.imaginary = 0.0; /* make "1" */
> a_point = subtract_complex(a_point,z_point);
> z_point = multiply_complex(z_point,a_point);
>
> ++num_iterations;
> } /* done iterating for one point */
>
> return num_iterations;
> }
>
> void complex_points_to_image(complex_t in_julia_coord_set[],
> int in_julia_set_size,
> int i_pixels_across,int i_pixels_down,
> int *out_image_final, int *out_max_colors)
> {
> int i, i_temp_quantize;
> double x_resolution_element, y_resolution_element, temp_quantize;
> double x_max, x_min, y_max, y_min;
>
> /* Convert the complex points into an image--
>
> first, find the min and max for each axis. */
>
> x_min = in_julia_coord_set[0].real;
> x_max = x_min;
> y_min = in_julia_coord_set[0].imaginary;
> y_max = y_min;
>
> for (i=0;i<in_julia_set_size;++i)
> {
> if (in_julia_coord_set[i].real < x_min)
> x_min = in_julia_coord_set[i].real;
> if (in_julia_coord_set[i].real > x_max)
> x_max = in_julia_coord_set[i].real;
> if (in_julia_coord_set[i].imaginary < y_min)
> y_min = in_julia_coord_set[i].imaginary;
> if (in_julia_coord_set[i].imaginary > y_max)
> y_max = in_julia_coord_set[i].imaginary;
> }
>
> /* convert to fit image grid size: */
>
> x_resolution_element = (x_max - x_min)/(double)i_pixels_across;
> y_resolution_element = (y_max - y_min)/(double)i_pixels_down;
>
> #ifdef DEBUG
> printf("%lf %lf\n",x_resolution_element, y_resolution_element);
> #endif
>
>
> /* now we can put each point into a grid space, and up the count of
> said grid. Since calloc initialized it to zero, this is safe */
> /* A point x,y goes to grid region i,j, where
> xi = (x-xmin)/x_resolution (position relative to far left)
> yj = (ymax-y)/y_resolution (position relative to top)
> This gets mapped to our 1-d array as xi + yj*i_pixels_across,
> taken as an integer (rounding down) and checking array limits
> */
>
> for (i=0; i<in_julia_set_size; ++i)
> {
> temp_quantize =
> (in_julia_coord_set[i].real - x_min)/x_resolution_element +
> (y_max - in_julia_coord_set[i].imaginary)/y_resolution_element
*
> (double)i_pixels_across;
> i_temp_quantize = (int)temp_quantize;
> if (i_temp_quantize > i_pixels_across*i_pixels_down)
> i_temp_quantize = i_pixels_across*i_pixels_down;
>
> #ifdef DEBUG
> printf("%d %lf %lf %lf %lf %lf %lf %lf\n",
> i_temp_quantize, temp_quantize, x_min, x_resolution_element,
> in_julia_coord_set[i].real, y_max, y_resolution_element,
> in_julia_coord_set[i].imaginary);
> #endif
>
> ++out_image_final[i_temp_quantize];
> }
>
> /* find the highest value (most occupied bin) */
> *out_max_colors = 0;
> for (i=0;i<i_pixels_across*i_pixels_down;++i)
> {
> if (out_image_final[i] > *out_max_colors)
> {
> *out_max_colors = out_image_final[i];
> }
> }
> }
>
> complex_t exponential_complex(complex_t in_complex)
> {
> complex_t out_complex;
> double e_to_x;
> /* taking the exponential, e^(x+iy) = e^xe^iy = e^x(cos(y)+isin(y)
*/
> e_to_x = exp(in_complex.real);
> out_complex.real = e_to_x * cos(in_complex.imaginary);
> out_complex.imaginary = e_to_x * sin(in_complex.imaginary);
> return out_complex;
> }
>
> complex_t multiply_complex(complex_t in_one_complex, complex_t
in_two_complex)
> {
> complex_t out_complex;
> /* multiplying complex variables-- (x+iy) * (a+ib) = xa - yb + i(xb +
ya) */
> out_complex.real = in_one_complex.real * in_two_complex.real -
> in_one_complex.imaginary * in_two_complex.imaginary;
> out_complex.imaginary = in_one_complex.real * in_two_complex.imaginary
+
> in_one_complex.imaginary * in_two_complex.real;
> return out_complex;
> }
>
> complex_t divide_complex(complex_t in_one_complex, complex_t
in_two_complex)
> {
> complex_t out_complex;
> double divider;
> /* dividing complex variables-- (x+iy)/(a+ib) = (xa - yb)/(a^2+b^2)
> + i (y*a-x*b)/(a^2+b^2) */
> divider = (in_two_complex.real*in_two_complex.real -
> in_two_complex.imaginary*in_two_complex.imaginary);
> out_complex.real = (in_one_complex.real * in_two_complex.real -
> in_one_complex.imaginary * in_two_complex.imaginary) / divider;
> out_complex.imaginary = (in_one_complex.imaginary *
in_two_complex.real -
> in_one_complex.real * in_two_complex.imaginary) / divider;
> return out_complex;
> }
>
> complex_t inverse_complex(complex_t in_complex)
> {
> complex_t out_complex;
> double divider;
> /* inverting a complex variable: 1/(x+iy) */
> divider = (in_complex.real*in_complex.real -
> in_complex.imaginary*in_complex.imaginary);
> out_complex.real = (in_complex.real - in_complex.imaginary) / divider;
> out_complex.imaginary = (in_complex.real - in_complex.imaginary) /
divider;
> return out_complex;
> }
>
> complex_t add_complex(complex_t in_one_complex, complex_t in_two_complex)
> {
> complex_t out_complex;
> /* adding complex variables-- do by component */
> out_complex.real = in_one_complex.real + in_two_complex.real;
> out_complex.imaginary = in_one_complex.imaginary +
in_two_complex.imaginary;
> return out_complex;
> }
>
> complex_t subtract_complex(complex_t in_one_complex, complex_t
in_two_complex)
> {
> complex_t out_complex;
> /* subtracting complex variables-- do by component */
> out_complex.real = in_one_complex.real - in_two_complex.real;
> out_complex.imaginary = in_one_complex.imaginary -
in_two_complex.imaginary;
> return out_complex;
> }
>
> double absolute_complex(complex_t in_complex)
> {
> double out_double;
> /* absolute value is (for x+yi) sqrt( x^2 + y^2) */
> out_double = sqrt(in_complex.real*in_complex.real +
in_complex.imaginary*in_complex.imaginary);
> return out_double;
> }
>
> /* This routine takes a complex coordinate point (x+iy) and a value
stating
> what the upper limit to the number of iterations is. It eventually
> returns an integer of how many counts the code iterated for within
> the given point/region until the exit condition ( abs(x+iy) > limit) was
met.
> This value is returned as an integer.
> */
>
> int single_mandelbrot_point(complex_t coord_point,
> complex_t c_constant,
> int Max_iterations, double divergent_limit)
> {
> complex_t z_point; /* we need a point to use in our calculation */
> int num_iterations; /* a counter to track the number of iterations
done */
>
> num_iterations = 0; /* zero our counter */
> z_point = coord_point; /* initialize to the given start coordinate */
>
> /* loop while the absolute value of the complex coordinate is < our
limit
> (for a mandelbrot) or until we've done our specified maximum number of
> iterations (both julia and mandelbrot) */
>
> while (absolute_complex(z_point) < divergent_limit &&
> num_iterations < Max_iterations)
> {
> /* z = z*z + c */
> z_point = multiply_complex(z_point,z_point);
> z_point = add_complex(z_point,c_constant);
>
> ++num_iterations;
> } /* done iterating for one point */
>
> return num_iterations;
> }
>
> void output_data(int *in_grid_array, int coord[4], int *out_grid_array,
int width, int height)
> {
> int i, j, k;
> int result;
>
> k = 0;
> for (j=coord[2]; j<=coord[3]; j++)
> {
> for (i=coord[0]; i<=coord[1]; i++)
> {
> out_grid_array[(j * height) + i] = in_grid_array[k];
> k++;
> }
> }
> if (!use_stdin)
> {
> result = MPI_Send(coord, 4, MPI_INT, 0, 0, comm);
> if (result != MPI_SUCCESS)
> {
> printf("Unable to send coordinates to the visualizer,
aborting.\n");
> MPI_Abort(MPI_COMM_WORLD, -1);
> }
> result = MPI_Send(in_grid_array, (coord[1] + 1 - coord[0]) *
(coord[3] + 1 - coord[2]), MPI_INT, 0, 0, comm);
> if (result != MPI_SUCCESS)
> {
> printf("Unable to send coordinate data to the visualizer,
aborting.\n");
> MPI_Abort(MPI_COMM_WORLD, -1);
> }
> }
> }
>
> /* Writes out a linear integer array of grayscale pixel values to a
> pgm-format file (with header). The exact format details are given
> at the end of this routine in case you don't have the man pages
> installed locally. Essentially, it's a 2-dimensional integer array
> of pixel grayscale values. Very easy to manipulate externally.
> */
>
> /* You need the following inputs:
> A linear integer array with the actual pixel values (read in as
> consecutive rows),
> The width and height of the grid, and
> The maximum pixel value (to set greyscale range). We are assuming
> that the lowest value is "0".
> */
>
> 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[])
> {
> FILE *ifp;
> int i, j, k;
> #ifdef USE_PPM
> int r, g, b;
> #endif
>
> printf("%s\nwidth: %d\nheight: %d\ncolors: %d\nstr: %s\n",
> filename, in_pixels_across, in_pixels_down, num_colors,
input_string);
> fflush(stdout);
>
> if ( (ifp=fopen(filename, "w")) == NULL)
> {
> printf("Error, could not open output file\n");
> MPI_Abort(MPI_COMM_WORLD, -1);
> exit(-1);
> }
>
> #ifdef USE_PPM
> fprintf(ifp, "P3\n"); /* specifies type of file, in this case ppm */
> fprintf(ifp, "# %s\n", input_string); /* an arbitrary file identifier
*/
> /* now give the file size in pixels by pixels */
> fprintf(ifp, "%d %d\n", in_pixels_across, in_pixels_down);
> /* give the max r,g,b level */
> fprintf(ifp, "255\n");
>
> k=0; /* counter for the linear array of the final image */
> /* assumes first point is upper left corner (element 0 of array) */
>
> if (in_max_pixel_value < 1)
> {
> for (j=0; j<in_pixels_down; ++j) /* start at the top row and work
down */
> {
> for (i=0; i<in_pixels_across; ++i) /* go along the row */
> {
> fprintf(ifp, "0 0 0 ");
> }
> fprintf(ifp, "\n"); /* done writing one row, begin next line */
> }
> }
> else
> {
> for (j=0; j<in_pixels_down; ++j) /* start at the top row and work
down */
> {
> for (i=0; i<in_pixels_across; ++i) /* go along the row */
> {
> getRGB(colors[(in_grid_array[k] * num_colors) /
in_max_pixel_value], &r, &g, &b);
> fprintf(ifp, "%d %d %d ", r, g, b); /* +1 since 0 = first
color */
> ++k;
> }
> fprintf(ifp, "\n"); /* done writing one row, begin next line */
> }
> }
> #else
> fprintf(ifp, "P2\n"); /* specifies type of file, in this case pgm */
> fprintf(ifp, "# %s\n", input_string); /* an arbitrary file identifier
*/
> /* now give the file size in pixels by pixels */
> fprintf(ifp, "%d %d\n", in_pixels_across, in_pixels_down);
> /* gives max number of grayscale levels */
> fprintf(ifp, "%d\n", in_max_pixel_value+1); /* plus 1 because 0=first
color */
>
> k=0; /* counter for the linear array of the final image */
> /* assumes first point is upper left corner (element 0 of array) */
>
> for (j=0;j<in_pixels_down;++j) /* start at the top row and work down
*/
> {
> for (i=0;i<in_pixels_across;++i) /* go along the row */
> {
> fprintf(ifp, "%d ", in_grid_array[k]+1); /* +1 since 0 = first
color */
> ++k;
> }
> fprintf(ifp, "\n"); /* done writing one row, begin next line */
> }
> #endif
>
> fclose(ifp);
> }
>
More information about the mpich-discuss
mailing list