hwloc: null check for cpuset
[mpich-dev.git] / examples / pmandel_service.c
1 /* -*- Mode: C; c-basic-offset:4 ; indent-tabs-mode:nil ; -*- */
2 /*
3  *  (C) 2001 by Argonne National Laboratory.
4  *      See COPYRIGHT in top-level directory.
5  */
6 #include <stdio.h>
7 #include <stdlib.h>
8 #include <math.h>
9 #include <string.h>
10 #include "mpi.h"
11 #ifdef HAVE_WINDOWS_H
12 #include <process.h>
13 #endif
14
15 /* definitions */
16
17 #define PROMPT 1
18 #define USE_PPM /* define to output a color image, otherwise output is a grayscale pgm image */
19
20 #ifndef MIN
21 #define MIN(x, y)       (((x) < (y))?(x):(y))
22 #endif
23 #ifndef MAX
24 #define MAX(x, y)       (((x) > (y))?(x):(y))
25 #endif
26
27 #define NOVALUE 99999       /* indicates a parameter is as of yet unspecified */
28 #define MAX_ITERATIONS 10000 /* maximum 'depth' to look for mandelbrot value */
29 #define INFINITE_LIMIT 4.0  /* evalue when fractal is considered diverging */
30 #define MAX_X_VALUE 2.0     /* the highest x can go for the mandelbrot, usually 1.0 */
31 #define MIN_X_VALUE -2.0    /* the lowest x can go for the mandelbrot, usually -1.0 */
32 #define MAX_Y_VALUE 2.0     /* the highest y can go for the mandelbrot, usually 1.0 */
33 #define MIN_Y_VALUE -2.0    /* the lowest y can go for the mandelbrot, usually -1.0 */
34 #define IDEAL_ITERATIONS 100 /* my favorite 'depth' to default to */
35 /* the ideal default x and y are currently the same as the max area */
36 #define IDEAL_MAX_X_VALUE 1.0
37 #define IDEAL_MIN_X_VALUE -1.0
38 #define IDEAL_MAX_Y_VALUE 1.0
39 #define IDEAL_MIN_Y_VALUE -1.0
40 #define MAX_WIDTH 2048       /* maximum size of image, across, in pixels */
41 #define MAX_HEIGHT 2048      /* maximum size of image, down, in pixels */
42 #define IDEAL_WIDTH 768       /* my favorite size of image, across, in pixels */
43 #define IDEAL_HEIGHT 768      /* my favorite size of image, down, in pixels */
44
45 #define RGBtocolor_t(r,g,b) ((color_t)(((unsigned char)(r)|((unsigned short)((unsigned char)(g))<<8))|(((unsigned long)(unsigned char)(b))<<16)))
46 #define getR(r) ((int)((r) & 0xFF))
47 #define getG(g) ((int)((g>>8) & 0xFF))
48 #define getB(b) ((int)((b>>16) & 0xFF))
49
50 typedef int color_t;
51
52 typedef struct complex_t
53 {
54   /* real + imaginary * sqrt(-1) i.e.   x + iy */
55   double real, imaginary;
56 } complex_t;
57
58 /* prototypes */
59
60 void read_mand_args(int argc, char *argv[], int *o_max_iterations,
61                     int *o_pixels_across, int *o_pixels_down,
62                     double *o_x_min, double *o_x_max,
63                     double *o_y_min, double *o_y_max, int *o_julia,
64                     double *o_julia_real_x, double *o_julia_imaginary_y,
65                     double *o_divergent_limit, int *o_alternate,
66                     char *filename, int *num_colors, int *use_stdin,
67                     int *save_image);
68 void check_mand_params(int *m_max_iterations,
69                        int *m_pixels_across, int *m_pixels_down,
70                        double *m_x_min, double *m_x_max,
71                        double *m_y_min, double *m_y_max,
72                        double *m_divergent_limit);
73 void check_julia_params(double *julia_x_constant, double *julia_y_constant);
74 int additive_mandelbrot_point(complex_t coord_point, 
75                             complex_t c_constant,
76                             int Max_iterations, double divergent_limit);
77 int additive_mandelbrot_point(complex_t coord_point, 
78                             complex_t c_constant,
79                             int Max_iterations, double divergent_limit);
80 int subtractive_mandelbrot_point(complex_t coord_point, 
81                             complex_t c_constant,
82                             int Max_iterations, double divergent_limit);
83 complex_t exponential_complex(complex_t in_complex);
84 complex_t multiply_complex(complex_t in_one_complex, complex_t in_two_complex);
85 complex_t divide_complex(complex_t in_one_complex, complex_t in_two_complex);
86 complex_t inverse_complex(complex_t in_complex);
87 complex_t add_complex(complex_t in_one_complex, complex_t in_two_complex);
88 complex_t subtract_complex(complex_t in_one_complex, complex_t in_two_complex);
89 double absolute_complex(complex_t in_complex);
90 void dumpimage(char *filename, int in_grid_array[], int in_pixels_across, int in_pixels_down,
91           int in_max_pixel_value, char input_string[], int num_colors, color_t colors[]);
92 int single_mandelbrot_point(complex_t coord_point, 
93                             complex_t c_constant,
94                             int Max_iterations, double divergent_limit);
95 color_t getColor(double fraction, double intensity);
96 int Make_color_array(int num_colors, color_t colors[]);
97 void output_data(int *in_grid_array, int coord[4], int *out_grid_array, int width, int height);
98 void PrintUsage(void);
99
100 #ifdef USE_PPM
101 const char *default_filename = "pmandel.ppm";
102 #else
103 const char *default_filename = "pmandel.pgm";
104 #endif
105
106 /* Solving the Mandelbrot set
107  
108    Set a maximum number of iterations to calculate before forcing a terminate
109    in the Mandelbrot loop.
110 */
111
112 /* Command-line parameters are all optional, and include:
113    -xmin # -xmax #      Limits for the x-axis for calculating and plotting
114    -ymin # -ymax #      Limits for the y-axis for calculating and plotting
115    -xscale # -yscale #  output pixel width and height
116    -depth #             Number of iterations before we give up on a given
117                         point and move on to calculate the next
118    -diverge #           Criteria for assuming the calculation has been
119                         "solved"-- for each point z, when
120                              abs(z=z^2+c) > diverge_value
121    -julia        Plot a Julia set instead of a Mandelbrot set
122    -jx # -jy #   The Julia point that defines the set
123    -alternate #    Plot an alternate Mandelbrot equation (varient 1 or 2 so far)
124 */
125
126 int myid;
127 int use_stdin = 0;
128 MPI_Comm comm;
129
130 void swap(int *i, int *j)
131 {
132     int x;
133     x = *i;
134     *i = *j;
135     *j = x;
136 }
137
138 int main(int argc, char *argv[])
139 {
140     complex_t coord_point, julia_constant;
141     double x_max, x_min, y_max, y_min, x_resolution, y_resolution;
142     double divergent_limit;
143     char file_message[160];
144     char filename[100];
145     int icount, imax_iterations;
146     int ipixels_across, ipixels_down;
147     int i, j, k, julia, alternate_equation;
148     int imin, imax, jmin, jmax;
149     int temp[5];
150     /* make an integer array of size [N x M] to hold answers. */
151     int *in_grid_array, *out_grid_array = NULL;
152     int numprocs;
153     int  namelen;
154     char processor_name[MPI_MAX_PROCESSOR_NAME];
155     int num_colors;
156     color_t *colors = NULL;
157     MPI_Status status;
158     int save_image = 0;
159     int result;
160     char mpi_port[MPI_MAX_PORT_NAME];
161     MPI_Info info;
162
163     MPI_Init(&argc, &argv);
164     MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
165     MPI_Comm_rank(MPI_COMM_WORLD, &myid);
166     MPI_Get_processor_name(processor_name, &namelen);
167
168     if (numprocs == 1)
169     {
170         PrintUsage();
171         MPI_Finalize();
172         exit(0);
173     }
174
175     if (myid == 0)
176     {
177         printf("Welcome to the Mandelbrot/Julia set explorer.\n");
178
179         /* Get inputs-- region to view (must be within x/ymin to x/ymax, make sure
180         xmax>xmin and ymax>ymin) and resolution (number of pixels along an edge,
181         N x M, i.e. 256x256)
182         */
183
184         read_mand_args(argc, argv, &imax_iterations, &ipixels_across, &ipixels_down,
185             &x_min, &x_max, &y_min, &y_max, &julia, &julia_constant.real,
186             &julia_constant.imaginary, &divergent_limit,
187             &alternate_equation, filename, &num_colors, &use_stdin, &save_image);
188         check_mand_params(&imax_iterations, &ipixels_across, &ipixels_down,
189             &x_min, &x_max, &y_min, &y_max, &divergent_limit);
190
191         if (julia == 1) /* we're doing a julia figure */
192             check_julia_params(&julia_constant.real, &julia_constant.imaginary);
193
194         MPI_Bcast(&use_stdin, 1, MPI_INT, 0, MPI_COMM_WORLD);
195         MPI_Bcast(&num_colors, 1, MPI_INT, 0, MPI_COMM_WORLD);
196         MPI_Bcast(&imax_iterations, 1, MPI_INT, 0, MPI_COMM_WORLD);
197         MPI_Bcast(&ipixels_across, 1, MPI_INT, 0, MPI_COMM_WORLD);
198         MPI_Bcast(&ipixels_down, 1, MPI_INT, 0, MPI_COMM_WORLD);
199         MPI_Bcast(&divergent_limit, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
200         MPI_Bcast(&julia, 1, MPI_INT, 0, MPI_COMM_WORLD);
201         MPI_Bcast(&julia_constant.real, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
202         MPI_Bcast(&julia_constant.imaginary, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
203         MPI_Bcast(&alternate_equation, 1, MPI_INT, 0, MPI_COMM_WORLD);
204     }
205     else
206     {
207         MPI_Bcast(&use_stdin, 1, MPI_INT, 0, MPI_COMM_WORLD);
208         MPI_Bcast(&num_colors, 1, MPI_INT, 0, MPI_COMM_WORLD);
209         MPI_Bcast(&imax_iterations, 1, MPI_INT, 0, MPI_COMM_WORLD);
210         MPI_Bcast(&ipixels_across, 1, MPI_INT, 0, MPI_COMM_WORLD);
211         MPI_Bcast(&ipixels_down, 1, MPI_INT, 0, MPI_COMM_WORLD);
212         MPI_Bcast(&divergent_limit, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
213         MPI_Bcast(&julia, 1, MPI_INT, 0, MPI_COMM_WORLD);
214         MPI_Bcast(&julia_constant.real, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
215         MPI_Bcast(&julia_constant.imaginary, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
216         MPI_Bcast(&alternate_equation, 1, MPI_INT, 0, MPI_COMM_WORLD);
217     }
218
219     if (myid == 0)
220     {
221         colors = malloc((num_colors+1)* sizeof(color_t));
222         if (colors == NULL)
223         {
224             MPI_Abort(MPI_COMM_WORLD, -1);
225             exit(-1);
226         }
227         Make_color_array(num_colors, colors);
228         colors[num_colors] = 0; /* add one on the top to avoid edge errors */
229     }
230
231     /* allocate memory */
232     if ( (in_grid_array = (int *)calloc(ipixels_across * ipixels_down, sizeof(int))) == NULL)
233     {
234         printf("Memory allocation failed for data array, aborting.\n");
235         MPI_Abort(MPI_COMM_WORLD, -1);
236         exit(-1);
237     }
238
239     if (myid == 0)
240     {
241         int istep, jstep, cur_proc;
242         int i1[400], i2[400], j1[400], j2[400];
243         int ii, jj;
244         char line[1024], *token;
245
246         if ( (out_grid_array = (int *)calloc(ipixels_across * ipixels_down, sizeof(int))) == NULL)
247         {
248             printf("Memory allocation failed for data array, aborting.\n");
249             MPI_Abort(MPI_COMM_WORLD, -1);
250             exit(-1);
251         }
252
253         srand(getpid());
254
255         if (!use_stdin)
256         {
257             result = MPI_Info_create(&info);
258             if (result != MPI_SUCCESS)
259             {
260                 printf("Unable to create an MPI_Info to be used in MPI_Open_port.\n");
261                 MPI_Abort(MPI_COMM_WORLD, -1);
262             }
263             result = MPI_Open_port(/*info*/MPI_INFO_NULL, mpi_port);
264             if (result != MPI_SUCCESS)
265             {
266                 printf("MPI_Open_port failed, aborting.\n");
267                 MPI_Abort(MPI_COMM_WORLD, -1);
268                 exit(-1);
269             }
270             result = MPI_Bcast(mpi_port, MPI_MAX_PORT_NAME, MPI_CHAR, 0, MPI_COMM_WORLD);
271             if (result != MPI_SUCCESS)
272             {
273                 printf("Unable to broadcast the port name to COMM_WORLD.\n");
274                 MPI_Abort(MPI_COMM_WORLD, -1);
275             }
276             printf("%s listening on port:\n%s\n", processor_name, mpi_port);
277             fflush(stdout);
278
279             result = MPI_Comm_accept(mpi_port, info, 0, MPI_COMM_WORLD, &comm);
280             if (result != MPI_SUCCESS)
281             {
282                 printf("MPI_Comm_accept failed, aborting.\n");
283                 MPI_Abort(MPI_COMM_WORLD, -1);
284                 exit(-1);
285             }
286             printf("accepted connection from visualization program.\n");
287             fflush(stdout);
288
289             printf("sending image size to visualizer.\n");
290             /*
291             printf("sending pixels_across.\n");
292             fflush(stdout);
293             */
294             result = MPI_Send(&ipixels_across, 1, MPI_INT, 0, 0, comm);
295             if (result != MPI_SUCCESS)
296             {
297                 printf("Unable to send pixel width to visualizer, aborting.\n");
298                 MPI_Abort(MPI_COMM_WORLD, -1);
299             }
300             /*
301             printf("sending pixels_down.\n");
302             fflush(stdout);
303             */
304             result = MPI_Send(&ipixels_down, 1, MPI_INT, 0, 0, comm);
305             if (result != MPI_SUCCESS)
306             {
307                 printf("Unable to send pixel height to visualizer, aborting.\n");
308                 MPI_Abort(MPI_COMM_WORLD, -1);
309             }
310             /*
311             printf("sending num_colors.\n");
312             fflush(stdout);
313             */
314             result = MPI_Send(&num_colors, 1, MPI_INT, 0, 0, comm);
315             if (result != MPI_SUCCESS)
316             {
317                 printf("Unable to send num_colors to visualizer, aborting.\n");
318                 MPI_Abort(MPI_COMM_WORLD, -1);
319             }
320             result = MPI_Send(&imax_iterations, 1, MPI_INT, 0, 0, comm);
321             if (result != MPI_SUCCESS)
322             {
323                 printf("Unable to send imax_iterations to visualizer, aborting.\n");
324                 MPI_Abort(MPI_COMM_WORLD, -1);
325             }
326         }
327
328         for (;;)
329         {
330             /* get x_min, x_max, y_min, and y_max */
331             if (use_stdin)
332             {
333                 printf("input xmin ymin xmax ymax max_iter, (0 0 0 0 0 to quit):\n");fflush(stdout);
334                 fgets(line, 1024, stdin);
335                 /*printf("read <%s> from stdin\n", line);fflush(stdout);*/
336                 token = strtok(line, " \n");
337                 x_min = atof(token);
338                 token = strtok(NULL, " \n");
339                 y_min = atof(token);
340                 token = strtok(NULL, " \n");
341                 x_max = atof(token);
342                 token = strtok(NULL, " \n");
343                 y_max = atof(token);
344                 token = strtok(NULL, " \n");
345                 imax_iterations = atoi(token);
346                 /*sscanf(line, "%g %g %g %g", &x_min, &y_min, &x_max, &y_max);*/
347                 /*scanf("%g %g %g %g", &x_min, &y_min, &x_max, &y_max);*/
348             }
349             else
350             {
351                 printf("reading xmin,ymin,xmax,ymax.\n");fflush(stdout);
352                 result = MPI_Recv(&x_min, 1, MPI_DOUBLE, 0, 0, comm, &status);
353                 if (result != MPI_SUCCESS)
354                 {
355                     printf("Unable to receive x_min from the visualizer, aborting.\n");
356                     MPI_Abort(MPI_COMM_WORLD, -1);
357                 }
358                 result = MPI_Recv(&y_min, 1, MPI_DOUBLE, 0, 0, comm, &status);
359                 if (result != MPI_SUCCESS)
360                 {
361                     printf("Unable to receive y_min from the visualizer, aborting.\n");
362                     MPI_Abort(MPI_COMM_WORLD, -1);
363                 }
364                 result = MPI_Recv(&x_max, 1, MPI_DOUBLE, 0, 0, comm, &status);
365                 if (result != MPI_SUCCESS)
366                 {
367                     printf("Unable to receive x_max from the visualizer, aborting.\n");
368                     MPI_Abort(MPI_COMM_WORLD, -1);
369                 }
370                 result = MPI_Recv(&y_max, 1, MPI_DOUBLE, 0, 0, comm, &status);
371                 if (result != MPI_SUCCESS)
372                 {
373                     printf("Unable to receive y_max from the visualizer, aborting.\n");
374                     MPI_Abort(MPI_COMM_WORLD, -1);
375                 }
376                 result = MPI_Recv(&imax_iterations, 1, MPI_INT, 0, 0, comm, &status);
377                 if (result != MPI_SUCCESS)
378                 {
379                     printf("Unable to receive imax_iterations from the visualizer, aborting.\n");
380                     MPI_Abort(MPI_COMM_WORLD, -1);
381                 }
382             }
383             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);
384
385             /* break the work up into 400 pieces */
386             istep = ipixels_across / 20;
387             jstep = ipixels_down / 20;
388             if (istep < 1)
389                 istep = 1;
390             if (jstep < 1)
391                 jstep = 1;
392             k = 0;
393             for (i=0; i<20; i++)
394             {
395                 for (j=0; j<20; j++)
396                 {
397                     i1[k] = MIN(istep * i, ipixels_across - 1);
398                     i2[k] = MIN((istep * (i+1)) - 1, ipixels_across - 1);
399                     j1[k] = MIN(jstep * j, ipixels_down - 1);
400                     j2[k] = MIN((jstep * (j+1)) - 1, ipixels_down - 1);
401                     k++;
402                 }
403             }
404
405             /* shuffle the work */
406             for (i=0; i<500; i++)
407             {
408                 ii = rand() % 400;
409                 jj = rand() % 400;
410                 swap(&i1[ii], &i1[jj]);
411                 swap(&i2[ii], &i2[jj]);
412                 swap(&j1[ii], &j1[jj]);
413                 swap(&j2[ii], &j2[jj]);
414             }
415
416             /*printf("bcasting the limits: (%f,%f)(%f,%f)\n", x_min, y_min, x_max, y_max);fflush(stdout);*/
417             /* let everyone know the limits */
418             MPI_Bcast(&x_min, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
419             MPI_Bcast(&x_max, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
420             MPI_Bcast(&y_min, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
421             MPI_Bcast(&y_max, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
422             MPI_Bcast(&imax_iterations, 1, MPI_INT, 0, MPI_COMM_WORLD);
423
424             /* check for the end condition */
425             if (x_min == x_max && y_min == y_max)
426             {
427                 break;
428             }
429
430             /* send a piece of work to each worker (there must be more work than workers) */
431             cur_proc = 1;
432             k = 0;
433             for (i=1; i<numprocs; i++)
434             {
435                 temp[0] = k+1;
436                 temp[1] = i1[k]; /* imin */
437                 temp[2] = i2[k]; /* imax */
438                 temp[3] = j1[k]; /* jmin */
439                 temp[4] = j2[k]; /* jmax */
440
441                 /*printf("sending work(%d) to %d\n", k+1, cur_proc);fflush(stdout);*/
442                 MPI_Send(temp, 5, MPI_INT, cur_proc, 100, MPI_COMM_WORLD);
443                 cur_proc++;
444                 k++;
445             }
446             /* receive the results and hand out more work until the image is complete */
447             while (k<400)
448             {
449                 MPI_Recv(temp, 5, MPI_INT, MPI_ANY_SOURCE, 200, MPI_COMM_WORLD, &status);
450                 cur_proc = temp[0];
451                 MPI_Recv(in_grid_array, (temp[2] + 1 - temp[1]) * (temp[4] + 1 - temp[3]), MPI_INT, cur_proc, 201, MPI_COMM_WORLD, &status);
452                 /* draw data */
453                 output_data(in_grid_array, &temp[1], out_grid_array, ipixels_across, ipixels_down);
454                 temp[0] = k+1;
455                 temp[1] = i1[k];
456                 temp[2] = i2[k];
457                 temp[3] = j1[k];
458                 temp[4] = j2[k];
459                 /*printf("sending work(%d) to %d\n", k+1, cur_proc);fflush(stdout);*/
460                 MPI_Send(temp, 5, MPI_INT, cur_proc, 100, MPI_COMM_WORLD);
461                 k++;
462             }
463             /* receive the last pieces of work */
464             /* and tell everyone to stop */
465             for (i=1; i<numprocs; i++)
466             {
467                 MPI_Recv(temp, 5, MPI_INT, MPI_ANY_SOURCE, 200, MPI_COMM_WORLD, &status);
468                 cur_proc = temp[0];
469                 MPI_Recv(in_grid_array, (temp[2] + 1 - temp[1]) * (temp[4] + 1 - temp[3]), MPI_INT, cur_proc, 201, MPI_COMM_WORLD, &status);
470                 /* draw data */
471                 output_data(in_grid_array, &temp[1], out_grid_array, ipixels_across, ipixels_down);
472                 temp[0] = 0;
473                 temp[1] = 0;
474                 temp[2] = 0;
475                 temp[3] = 0;
476                 temp[4] = 0;
477                 /*printf("sending %d to tell %d to stop\n", temp[0], cur_proc);fflush(stdout);*/
478                 MPI_Send(temp, 5, MPI_INT, cur_proc, 100, MPI_COMM_WORLD);
479             }
480
481             /* tell the visualizer the image is done */
482             if (!use_stdin)
483             {
484                 temp[0] = 0;
485                 temp[1] = 0;
486                 temp[2] = 0;
487                 temp[3] = 0;
488                 result = MPI_Send(temp, 4, MPI_INT, 0, 0, comm);
489                 if (result != MPI_SUCCESS)
490                 {
491                     printf("Unable to send a done command to the visualizer, aborting.\n");
492                     MPI_Abort(MPI_COMM_WORLD, -1);
493                 }
494             }
495         }
496     }
497     else
498     {
499         if (!use_stdin)
500         {
501             result = MPI_Info_create(&info);
502             if (result != MPI_SUCCESS)
503             {
504                 printf("Unable to create an MPI_Info to be used in MPI_Open_port.\n");
505                 MPI_Abort(MPI_COMM_WORLD, -1);
506             }
507             result = MPI_Bcast(mpi_port, MPI_MAX_PORT_NAME, MPI_CHAR, 0, MPI_COMM_WORLD);
508             if (result != MPI_SUCCESS)
509             {
510                 printf("Unable to broadcast the port name on COMM_WORLD.\n");
511                 MPI_Abort(MPI_COMM_WORLD, -1);
512             }
513             /*
514             printf("%s listening on port: %s\n", processor_name, mpi_port);
515             fflush(stdout);
516             */
517
518             result = MPI_Comm_accept(mpi_port, info, 0, MPI_COMM_WORLD, &comm);
519             if (result != MPI_SUCCESS)
520             {
521                 printf("MPI_Comm_accept failed, aborting.\n");
522                 MPI_Abort(MPI_COMM_WORLD, -1);
523                 exit(-1);
524             }
525         }
526         for (;;)
527         {
528             MPI_Bcast(&x_min, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
529             MPI_Bcast(&x_max, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
530             MPI_Bcast(&y_min, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
531             MPI_Bcast(&y_max, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
532             MPI_Bcast(&imax_iterations, 1, MPI_INT, 0, MPI_COMM_WORLD);
533
534             /* check for the end condition */
535             if (x_min == x_max && y_min == y_max)
536             {
537                 break;
538             }
539
540             x_resolution = (x_max-x_min)/ ((double)ipixels_across);
541             y_resolution = (y_max-y_min)/ ((double)ipixels_down);
542
543             MPI_Recv(temp, 5, MPI_INT, 0, 100, MPI_COMM_WORLD, &status);
544             while (temp[0] != 0)
545             {
546                 imin = temp[1];
547                 imax = temp[2];
548                 jmin = temp[3];
549                 jmax = temp[4];
550
551                 k = 0;
552                 for (j=jmin; j<=jmax; ++j)
553                 {
554                     coord_point.imaginary = y_max - j*y_resolution; /* go top to bottom */
555
556                     for (i=imin; i<=imax; ++i)
557                     {
558                         /* Call Mandelbrot routine for each code, fill array with number of iterations. */
559
560                         coord_point.real = x_min + i*x_resolution; /* go left to right */
561                         if (julia == 1)
562                         {
563                             /* doing Julia set */
564                             /* julia eq:  z = z^2 + c, z_0 = grid coordinate, c = constant */
565                             icount = single_mandelbrot_point(coord_point, julia_constant, imax_iterations, divergent_limit);
566                         }
567                         else if (alternate_equation == 1)
568                         {
569                             /* doing experimental form 1 */
570                             icount = subtractive_mandelbrot_point(coord_point, julia_constant, imax_iterations, divergent_limit);
571                         }
572                         else if (alternate_equation == 2)
573                         {
574                             /* doing experimental form 2 */
575                             icount = additive_mandelbrot_point(coord_point, julia_constant, imax_iterations, divergent_limit);
576                         }
577                         else
578                         {
579                             /* default to doing Mandelbrot set */
580                             /* mandelbrot eq: z = z^2 + c, z_0 = c, c = grid coordinate */
581                             icount = single_mandelbrot_point(coord_point, coord_point, imax_iterations, divergent_limit);
582                         }
583                         in_grid_array[k] = icount;
584                         ++k;
585                     }
586                 }
587                 /* send the result to the root */
588                 temp[0] = myid;
589                 MPI_Send(temp, 5, MPI_INT, 0, 200, MPI_COMM_WORLD);
590                 MPI_Send(in_grid_array, (temp[2] + 1 - temp[1]) * (temp[4] + 1 - temp[3]), MPI_INT, 0, 201, MPI_COMM_WORLD);
591                 /* get the next piece of work */
592                 MPI_Recv(temp, 5, MPI_INT, 0, 100, MPI_COMM_WORLD, &status);
593             }
594         }
595     }
596
597     if (myid == 0 && save_image)
598     {
599         imax_iterations = 0;
600         for (i=0; i<ipixels_across * ipixels_down; ++i)
601         {
602             /* look for "brightest" pixel value, for image use */
603             if (out_grid_array[i] > imax_iterations)
604                 imax_iterations = out_grid_array[i];
605         }
606
607         if (julia == 0)
608             printf("Done calculating mandelbrot, now creating file\n");
609         else
610             printf("Done calculating julia, now creating file\n");
611         fflush(stdout);
612
613         /* Print out the array in some appropriate form. */
614         if (julia == 0)
615         {
616             /* it's a mandelbrot */
617             sprintf(file_message, "Mandelbrot over (%lf-%lf,%lf-%lf), size %d x %d",
618                 x_min, x_max, y_min, y_max, ipixels_across, ipixels_down);
619         }
620         else
621         {
622             /* it's a julia */
623             sprintf(file_message, "Julia over (%lf-%lf,%lf-%lf), size %d x %d, center (%lf, %lf)",
624                 x_min, x_max, y_min, y_max, ipixels_across, ipixels_down,
625                 julia_constant.real, julia_constant.imaginary);
626         }
627
628         dumpimage(filename, out_grid_array, ipixels_across, ipixels_down, imax_iterations, file_message, num_colors, colors);
629     }
630
631     MPI_Comm_disconnect(&comm);
632
633     MPI_Finalize();
634     free(colors);
635     return 0;
636 } /* end of main */
637
638 void PrintUsage(void)
639 {
640     printf("usage: mpiexec -n x pmandel [options]\n");
641     printf("options:\n -xmin # -xmax #\n -ymin # -ymax #\n -depth #\n -xscale # -yscale #\n -out filename\n -i\n");
642     printf("All options are optional.\n");
643     printf("-i will allow you to input the min/max parameters from stdin and output the resulting image to a ppm file.");
644     printf("  Otherwise the root process will listen for a separate visualizer program to connect to it.\n");
645     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",
646         IDEAL_MIN_X_VALUE, IDEAL_MAX_X_VALUE, IDEAL_MIN_Y_VALUE, IDEAL_MAX_Y_VALUE, IDEAL_ITERATIONS,
647         IDEAL_WIDTH, IDEAL_HEIGHT, default_filename);
648     fflush(stdout);
649 }
650
651 color_t getColor(double fraction, double intensity)
652 {
653     /* fraction is a part of the rainbow (0.0 - 1.0) = (Red-Yellow-Green-Cyan-Blue-Magenta-Red)
654     intensity (0.0 - 1.0) 0 = black, 1 = full color, 2 = white
655     */
656     double red, green, blue;
657     int r,g,b;
658     double dtemp;
659
660     fraction = fabs(modf(fraction, &dtemp));
661
662     if (intensity > 2.0)
663         intensity = 2.0;
664     if (intensity < 0.0)
665         intensity = 0.0;
666
667     dtemp = 1.0/6.0;
668
669     if (fraction < 1.0/6.0)
670     {
671         red = 1.0;
672         green = fraction / dtemp;
673         blue = 0.0;
674     }
675     else
676     {
677         if (fraction < 1.0/3.0)
678         {
679             red = 1.0 - ((fraction - dtemp) / dtemp);
680             green = 1.0;
681             blue = 0.0;
682         }
683         else
684         {
685             if (fraction < 0.5)
686             {
687                 red = 0.0;
688                 green = 1.0;
689                 blue = (fraction - (dtemp*2.0)) / dtemp;
690             }
691             else
692             {
693                 if (fraction < 2.0/3.0)
694                 {
695                     red = 0.0;
696                     green = 1.0 - ((fraction - (dtemp*3.0)) / dtemp);
697                     blue = 1.0;
698                 }
699                 else
700                 {
701                     if (fraction < 5.0/6.0)
702                     {
703                         red = (fraction - (dtemp*4.0)) / dtemp;
704                         green = 0.0;
705                         blue = 1.0;
706                     }
707                     else
708                     {
709                         red = 1.0;
710                         green = 0.0;
711                         blue = 1.0 - ((fraction - (dtemp*5.0)) / dtemp);
712                     }
713                 }
714             }
715         }
716     }
717
718     if (intensity > 1)
719     {
720         intensity = intensity - 1.0;
721         red = red + ((1.0 - red) * intensity);
722         green = green + ((1.0 - green) * intensity);
723         blue = blue + ((1.0 - blue) * intensity);
724     }
725     else
726     {
727         red = red * intensity;
728         green = green * intensity;
729         blue = blue * intensity;
730     }
731
732     r = (int)(red * 255.0);
733     g = (int)(green * 255.0);
734     b = (int)(blue * 255.0);
735
736     return RGBtocolor_t(r,g,b);
737 }
738
739 int Make_color_array(int num_colors, color_t colors[])
740 {
741     double fraction, intensity;
742     int i;
743
744     intensity = 1.0;
745     for (i=0; i<num_colors; i++)
746     {
747         fraction = (double)i / (double)num_colors;
748         colors[i] = getColor(fraction, intensity);
749     }
750     return 0;
751 }
752
753 void getRGB(color_t color, int *r, int *g, int *b)
754 {
755     *r = getR(color);
756     *g = getG(color);
757     *b = getB(color);
758 }
759
760 void read_mand_args(int argc, char *argv[], int *o_max_iterations,
761                     int *o_pixels_across, int *o_pixels_down,
762                     double *o_x_min, double *o_x_max,
763                     double *o_y_min, double *o_y_max, int *o_julia,
764                     double *o_julia_real_x, double *o_julia_imaginary_y,
765                     double *o_divergent_limit, int *o_alternate,
766                     char *filename, int *o_num_colors, int *use_stdin,
767                     int *save_image)
768 {       
769     int i;
770
771     *o_julia_real_x = NOVALUE;
772     *o_julia_imaginary_y = NOVALUE;
773
774     /* set defaults */
775     *o_max_iterations = IDEAL_ITERATIONS;
776     *o_pixels_across = IDEAL_WIDTH;
777     *o_pixels_down = IDEAL_HEIGHT;
778     *o_x_min = IDEAL_MIN_X_VALUE;
779     *o_x_max = IDEAL_MAX_X_VALUE;
780     *o_y_min = IDEAL_MIN_Y_VALUE;
781     *o_y_max = IDEAL_MAX_Y_VALUE;
782     *o_divergent_limit = INFINITE_LIMIT;
783     strcpy(filename, default_filename);
784     *o_num_colors = IDEAL_ITERATIONS;
785     *use_stdin = 0; /* default is to listen for a controller */
786     *save_image = NOVALUE;
787
788     *o_julia = 0; /* default is "generate Mandelbrot" */
789     *o_alternate = 0; /* default is still "generate Mandelbrot" */
790     *o_divergent_limit = INFINITE_LIMIT; /* default total range is assumed
791                                          if not explicitly overwritten */
792
793     /* We just cycle through all given parameters, matching what we can.
794     Note that we force casting, because we expect that a nonsensical
795     value is better than a poorly formatted one (fewer crashes), and
796     we'll later sort out the good from the bad
797     */
798
799     for (i = 0; i < argc; ++i) /* grab command line arguments */
800         if (strcmp(argv[i], "-xmin\0") == 0 && argv[i+1] != NULL) 
801             sscanf(argv[i+1], "%lf", &*o_x_min);
802         else if (strcmp(argv[i], "-xmax\0") == 0 && argv[i+1] != NULL) 
803             sscanf(argv[i+1], "%lf", &*o_x_max);
804         else if (strcmp(argv[i], "-ymin\0") == 0 && argv[i+1] != NULL) 
805             sscanf(argv[i+1], "%lf", &*o_y_min);
806         else if (strcmp(argv[i], "-ymax\0") == 0 && argv[i+1] != NULL) 
807             sscanf(argv[i+1], "%lf", &*o_y_max);
808         else if (strcmp(argv[i], "-depth\0") == 0 && argv[i+1] != NULL) 
809             sscanf(argv[i+1], "%d", &*o_max_iterations);
810         else if (strcmp(argv[i], "-xscale\0") == 0 && argv[i+1] != NULL) 
811             sscanf(argv[i+1], "%d", &*o_pixels_across);
812         else if (strcmp(argv[i], "-yscale\0") == 0 && argv[i+1] != NULL) 
813             sscanf(argv[i+1], "%d", &*o_pixels_down);
814         else if (strcmp(argv[i], "-diverge\0") == 0 && argv[i+1] != NULL) 
815             sscanf(argv[i+1], "%lf", &*o_divergent_limit);
816         else if (strcmp(argv[i], "-jx\0") == 0 && argv[i+1] != NULL) 
817             sscanf(argv[i+1], "%lf", &*o_julia_real_x);
818         else if (strcmp(argv[i], "-jy\0") == 0 && argv[i+1] != NULL) 
819             sscanf(argv[i+1], "%lf", &*o_julia_imaginary_y);
820         else if (strcmp(argv[i], "-alternate\0") == 0 && argv[i+1] != NULL)
821             sscanf(argv[i+1], "%d", &*o_alternate);
822         else if (strcmp(argv[i], "-julia\0") == 0)
823             *o_julia = 1;
824         else if (strcmp(argv[i], "-out\0") == 0 && argv[i+1] != NULL)
825             strcpy(filename, argv[i+1]);
826         else if (strcmp(argv[i], "-colors\0") == 0 && argv[i+1] != NULL)
827             sscanf(argv[i+1], "%d", &*o_num_colors);
828         else if (strcmp(argv[i], "-i\0") == 0)
829             *use_stdin = 1;
830         else if (strcmp(argv[i], "-save\0") == 0)
831             *save_image = 1;
832         else if (strcmp(argv[i], "-nosave\0") == 0)
833             *save_image = 0;
834         else if (strcmp(argv[i], "-help\0") == 0)
835         {
836             PrintUsage();
837             MPI_Finalize();
838             exit(0);
839         }
840
841         if (*save_image == NOVALUE)
842         {
843             if (*use_stdin == 1)
844                 *save_image = 1;
845             else
846                 *save_image = 0;
847         }
848 #if DEBUG2
849         if (myid == 0)
850         {
851             printf("Doing %d iterations over (%.2lf:%.2lf,%.2lf:%.2lf) on a %d x %d grid\n",
852                 *o_max_iterations, *o_x_min,*o_x_max,*o_y_min,*o_y_max,
853                 *o_pixels_across, *o_pixels_down);
854         }
855 #endif
856 }
857
858 void check_mand_params(int *m_max_iterations,
859                        int *m_pixels_across, int *m_pixels_down,
860                        double *m_x_min, double *m_x_max,
861                        double *m_y_min, double *m_y_max,
862                        double *m_divergent_limit)
863 {
864     /* Get first batch of limits */
865 #if PROMPT
866     if (*m_x_min == NOVALUE || *m_x_max == NOVALUE)
867     {
868         /* grab unspecified limits */
869         printf("Enter lower and higher limits on x (within -1 to 1): ");
870         scanf("%lf %lf", &*m_x_min, &*m_x_max);
871     }
872
873     if (*m_y_min == NOVALUE || *m_y_max == NOVALUE)
874     {
875         /* grab unspecified limits */
876         printf("Enter lower and higher limits on y (within -1 to 1): ");
877         scanf("%lf %lf", &*m_y_min, &*m_y_max);
878     }
879 #endif
880
881     /* Verify limits are reasonable */
882
883     if (*m_x_min < MIN_X_VALUE || *m_x_min > *m_x_max)
884     {
885         printf("Chosen lower x limit is too low, resetting to %lf\n",MIN_X_VALUE);
886         *m_x_min = MIN_X_VALUE;
887     }
888     if (*m_x_max > MAX_X_VALUE || *m_x_max <= *m_x_min)
889     {
890         printf("Chosen upper x limit is too high, resetting to %lf\n",MAX_X_VALUE);
891         *m_x_max = MAX_X_VALUE;
892     }
893     if (*m_y_min < MIN_Y_VALUE || *m_y_min > *m_y_max)
894     {
895         printf("Chosen lower y limit is too low, resetting to %lf\n",MIN_Y_VALUE);
896         *m_y_min = MIN_Y_VALUE;
897     }
898     if (*m_y_max > MAX_Y_VALUE || *m_y_max <= *m_y_min)
899     {
900         printf("Chosen upper y limit is too high, resetting to %lf\n",MAX_Y_VALUE);
901         *m_y_max = MAX_Y_VALUE;
902     }
903
904     /* Get second set of limits */
905 #if PROMPT
906
907     if (*m_max_iterations == NOVALUE)
908     {
909         /* grab unspecified limit */
910         printf("Enter max interations to do: ");
911         scanf("%d",&*m_max_iterations);
912     }
913 #endif
914
915     /* Verify second set of limits */
916
917     if (*m_max_iterations > MAX_ITERATIONS || *m_max_iterations < 0)
918     {
919         printf("Warning, unreasonable number of iterations, setting to %d\n", MAX_ITERATIONS);
920         *m_max_iterations = MAX_ITERATIONS;
921     }
922
923     /* Get third set of limits */
924 #if PROMPT
925     if (*m_pixels_across == NOVALUE || *m_pixels_down == NOVALUE)
926     {
927         /* grab unspecified limits */
928         printf("Enter pixel size (horizonal by vertical, i.e. 256 256): ");
929         scanf(" %d %d", &*m_pixels_across, &*m_pixels_down);
930     }
931 #endif
932
933     /* Verify third set of limits */
934
935     if (*m_pixels_across > MAX_WIDTH)
936     {
937         printf("Warning, image requested is too wide, setting to %d pixel width\n", MAX_WIDTH);
938         *m_pixels_across = MAX_WIDTH;
939     }
940     if (*m_pixels_down > MAX_HEIGHT)
941     {
942         printf("Warning, image requested is too tall, setting to %d pixel height\n", MAX_HEIGHT);
943         *m_pixels_down = MAX_HEIGHT;
944     }
945
946 #if DEBUG2
947     printf("%d iterations over (%.2lf:%.2lf,%.2lf:%.2lf), %d x %d grid, diverge value %lf\n",
948         *m_max_iterations, *m_x_min,*m_x_max,*m_y_min,*m_y_max,
949         *m_pixels_across, *m_pixels_down, *m_divergent_limit);
950 #endif
951 }
952
953 void check_julia_params(double *julia_x_constant, double *julia_y_constant)
954 {       
955
956     /* Hey, we're not stupid-- if they must Prompt, we will also be Noisy */
957 #if PROMPT
958     if (*julia_x_constant == NOVALUE || *julia_y_constant == NOVALUE)
959     {
960         /* grab unspecified limits */
961         printf("Enter Coordinates for Julia expansion: ");
962         scanf("%lf %lf", &*julia_x_constant, &*julia_y_constant);
963     }
964 #endif
965
966 #if DEBUG2
967     /* In theory, any point can be investigated, although it's not much point
968     if it's outside of the area viewed.  But, hey, that's not our problem */
969     printf("Exploring julia set around (%lf, %lf)\n", *julia_x_constant, *julia_y_constant);
970 #endif
971 }
972
973 /* This is a Mandelbrot variant, solving the code for the equation:
974 z = z(1-z)
975 */
976
977 /* This routine takes a complex coordinate point (x+iy) and a value stating
978 what the upper limit to the number of iterations is.  It eventually
979 returns an integer of how many counts the code iterated for within
980 the given point/region until the exit condition ( abs(x+iy) > limit) was met.
981 This value is returned as an integer.
982 */
983
984 int subtractive_mandelbrot_point(complex_t coord_point, 
985                                  complex_t c_constant,
986                                  int Max_iterations, double divergent_limit)
987 {
988     complex_t z_point, a_point; /* we need 2 pts to use in our calculation */
989     int num_iterations;  /* a counter to track the number of iterations done */
990
991     num_iterations = 0; /* zero our counter */
992     z_point = coord_point; /* initialize to the given start coordinate */
993
994     /* loop while the absolute value of the complex coordinate is < our limit
995     (for a mandelbrot) or until we've done our specified maximum number of
996     iterations (both julia and mandelbrot) */
997
998     while (absolute_complex(z_point) < divergent_limit &&
999         num_iterations < Max_iterations)
1000     {
1001         /* z = z(1-z) */
1002         a_point.real = 1.0; a_point.imaginary = 0.0; /* make "1" */
1003         a_point = subtract_complex(a_point,z_point);
1004         z_point = multiply_complex(z_point,a_point);
1005
1006         ++num_iterations;
1007     } /* done iterating for one point */
1008
1009     return num_iterations;
1010 }
1011
1012 /* This is a Mandelbrot variant, solving the code for the equation:
1013 z = z(z+c)
1014 */
1015
1016 /* This routine takes a complex coordinate point (x+iy) and a value stating
1017 what the upper limit to the number of iterations is.  It eventually
1018 returns an integer of how many counts the code iterated for within
1019 the given point/region until the exit condition ( abs(x+iy) > limit) was met.
1020 This value is returned as an integer.
1021 */
1022
1023 int additive_mandelbrot_point(complex_t coord_point, 
1024                               complex_t c_constant,
1025                               int Max_iterations, double divergent_limit)
1026 {
1027     complex_t z_point, a_point; /* we need 2 pts to use in our calculation */
1028     int num_iterations;  /* a counter to track the number of iterations done */
1029
1030     num_iterations = 0; /* zero our counter */
1031     z_point = coord_point; /* initialize to the given start coordinate */
1032
1033     /* loop while the absolute value of the complex coordinate is < our limit
1034     (for a mandelbrot) or until we've done our specified maximum number of
1035     iterations (both julia and mandelbrot) */
1036
1037     while (absolute_complex(z_point) < divergent_limit &&
1038         num_iterations < Max_iterations)
1039     {
1040         /* z = z(z+C) */
1041         a_point = add_complex(z_point,coord_point);
1042         z_point = multiply_complex(z_point,a_point);
1043
1044         ++num_iterations;
1045     } /* done iterating for one point */
1046
1047     return num_iterations;
1048 }
1049
1050 /* This is a Mandelbrot variant, solving the code for the equation:
1051 z = z(z+1)
1052 */
1053
1054 /* This routine takes a complex coordinate point (x+iy) and a value stating
1055 what the upper limit to the number of iterations is.  It eventually
1056 returns an integer of how many counts the code iterated for within
1057 the given point/region until the exit condition ( abs(x+iy) > limit) was met.
1058 This value is returned as an integer.
1059 */
1060
1061 int exponential_mandelbrot_point(complex_t coord_point, 
1062                                  complex_t c_constant,
1063                                  int Max_iterations, double divergent_limit)
1064 {
1065     complex_t z_point, a_point; /* we need 2 pts to use in our calculation */
1066     int num_iterations;  /* a counter to track the number of iterations done */
1067
1068     num_iterations = 0; /* zero our counter */
1069     z_point = coord_point; /* initialize to the given start coordinate */
1070
1071     /* loop while the absolute value of the complex coordinate is < our limit
1072     (for a mandelbrot) or until we've done our specified maximum number of
1073     iterations (both julia and mandelbrot) */
1074
1075     while (absolute_complex(z_point) < divergent_limit &&
1076         num_iterations < Max_iterations)
1077     {
1078         /* z = z(1-z) */
1079         a_point.real = 1.0; a_point.imaginary = 0.0; /* make "1" */
1080         a_point = subtract_complex(a_point,z_point);
1081         z_point = multiply_complex(z_point,a_point);
1082
1083         ++num_iterations;
1084     } /* done iterating for one point */
1085
1086     return num_iterations;
1087 }
1088
1089 void complex_points_to_image(complex_t in_julia_coord_set[],
1090                              int in_julia_set_size,
1091                              int i_pixels_across,int i_pixels_down,
1092                              int *out_image_final, int *out_max_colors)
1093 {
1094     int i, i_temp_quantize;
1095     double x_resolution_element, y_resolution_element, temp_quantize;
1096     double x_max, x_min, y_max, y_min;
1097
1098     /* Convert the complex points into an image--
1099
1100     first, find the min and max for each axis. */
1101
1102     x_min = in_julia_coord_set[0].real;
1103     x_max = x_min;
1104     y_min = in_julia_coord_set[0].imaginary;
1105     y_max = y_min;
1106
1107     for (i=0;i<in_julia_set_size;++i)
1108     {
1109         if (in_julia_coord_set[i].real < x_min)
1110             x_min = in_julia_coord_set[i].real;
1111         if (in_julia_coord_set[i].real > x_max)
1112             x_max = in_julia_coord_set[i].real;
1113         if (in_julia_coord_set[i].imaginary < y_min)
1114             y_min = in_julia_coord_set[i].imaginary;
1115         if (in_julia_coord_set[i].imaginary > y_max)
1116             y_max = in_julia_coord_set[i].imaginary;
1117     }
1118
1119     /* convert to fit image grid size: */
1120
1121     x_resolution_element =  (x_max - x_min)/(double)i_pixels_across;
1122     y_resolution_element =  (y_max - y_min)/(double)i_pixels_down;
1123
1124 #ifdef DEBUG
1125     printf("%lf %lf\n",x_resolution_element, y_resolution_element);
1126 #endif
1127
1128
1129     /* now we can put each point into a grid space, and up the count of
1130     said grid.  Since calloc initialized it to zero, this is safe */
1131     /* A point x,y goes to grid region i,j, where
1132     xi =  (x-xmin)/x_resolution   (position relative to far left)
1133     yj =  (ymax-y)/y_resolution   (position relative to top)
1134     This gets mapped to our 1-d array as  xi + yj*i_pixels_across,
1135     taken as an integer (rounding down) and checking array limits
1136     */
1137
1138     for (i=0; i<in_julia_set_size; ++i)
1139     {
1140         temp_quantize =
1141             (in_julia_coord_set[i].real - x_min)/x_resolution_element +
1142             (y_max -  in_julia_coord_set[i].imaginary)/y_resolution_element *
1143             (double)i_pixels_across;
1144         i_temp_quantize = (int)temp_quantize;
1145         if (i_temp_quantize > i_pixels_across*i_pixels_down)
1146             i_temp_quantize = i_pixels_across*i_pixels_down;
1147
1148 #ifdef DEBUG
1149         printf("%d %lf %lf %lf %lf %lf %lf %lf\n",
1150             i_temp_quantize, temp_quantize, x_min, x_resolution_element,
1151             in_julia_coord_set[i].real, y_max, y_resolution_element,
1152             in_julia_coord_set[i].imaginary);
1153 #endif
1154
1155         ++out_image_final[i_temp_quantize];
1156     }
1157
1158     /* find the highest value (most occupied bin) */
1159     *out_max_colors = 0;
1160     for (i=0;i<i_pixels_across*i_pixels_down;++i)
1161     {
1162         if (out_image_final[i] > *out_max_colors)
1163         {
1164             *out_max_colors = out_image_final[i];
1165         }
1166     }
1167 }
1168
1169 complex_t exponential_complex(complex_t in_complex)
1170 {
1171     complex_t out_complex;
1172     double e_to_x;
1173     /* taking the exponential,   e^(x+iy) = e^xe^iy = e^x(cos(y)+isin(y) */
1174     e_to_x = exp(in_complex.real);
1175     out_complex.real = e_to_x * cos(in_complex.imaginary);
1176     out_complex.imaginary = e_to_x * sin(in_complex.imaginary);
1177     return out_complex;
1178 }
1179
1180 complex_t multiply_complex(complex_t in_one_complex, complex_t in_two_complex)
1181 {
1182     complex_t out_complex;
1183     /* multiplying complex variables-- (x+iy) * (a+ib) = xa - yb + i(xb + ya) */
1184     out_complex.real = in_one_complex.real * in_two_complex.real -
1185         in_one_complex.imaginary * in_two_complex.imaginary;
1186     out_complex.imaginary = in_one_complex.real * in_two_complex.imaginary +
1187         in_one_complex.imaginary * in_two_complex.real;
1188     return out_complex;
1189 }
1190
1191 complex_t divide_complex(complex_t in_one_complex, complex_t in_two_complex)
1192 {
1193     complex_t out_complex;
1194     double divider;
1195     /* dividing complex variables-- (x+iy)/(a+ib) = (xa - yb)/(a^2+b^2)
1196     + i (y*a-x*b)/(a^2+b^2) */
1197     divider = (in_two_complex.real*in_two_complex.real -
1198         in_two_complex.imaginary*in_two_complex.imaginary);
1199     out_complex.real = (in_one_complex.real * in_two_complex.real -
1200         in_one_complex.imaginary * in_two_complex.imaginary) / divider;
1201     out_complex.imaginary = (in_one_complex.imaginary * in_two_complex.real -
1202         in_one_complex.real * in_two_complex.imaginary) / divider;
1203     return out_complex;
1204 }
1205
1206 complex_t inverse_complex(complex_t in_complex)
1207 {
1208     complex_t out_complex;
1209     double divider;
1210     /* inverting a complex variable: 1/(x+iy) */
1211     divider = (in_complex.real*in_complex.real -
1212         in_complex.imaginary*in_complex.imaginary);
1213     out_complex.real = (in_complex.real - in_complex.imaginary) / divider;
1214     out_complex.imaginary = (in_complex.real - in_complex.imaginary) / divider;
1215     return out_complex;
1216 }
1217
1218 complex_t add_complex(complex_t in_one_complex, complex_t in_two_complex)
1219 {
1220     complex_t out_complex;
1221     /* adding complex variables-- do by component */
1222     out_complex.real = in_one_complex.real + in_two_complex.real;
1223     out_complex.imaginary = in_one_complex.imaginary + in_two_complex.imaginary;
1224     return out_complex;
1225 }
1226
1227 complex_t subtract_complex(complex_t in_one_complex, complex_t in_two_complex)
1228 {
1229     complex_t out_complex;
1230     /* subtracting complex variables-- do by component */
1231     out_complex.real = in_one_complex.real - in_two_complex.real;
1232     out_complex.imaginary = in_one_complex.imaginary - in_two_complex.imaginary;
1233     return out_complex;
1234 }
1235
1236 double absolute_complex(complex_t in_complex)
1237 {
1238     double out_double;
1239     /* absolute value is (for x+yi)  sqrt( x^2 + y^2) */
1240     out_double = sqrt(in_complex.real*in_complex.real + in_complex.imaginary*in_complex.imaginary);
1241     return out_double;
1242 }
1243
1244 /* This routine takes a complex coordinate point (x+iy) and a value stating
1245 what the upper limit to the number of iterations is.  It eventually
1246 returns an integer of how many counts the code iterated for within
1247 the given point/region until the exit condition ( abs(x+iy) > limit) was met.
1248 This value is returned as an integer.
1249 */
1250
1251 int single_mandelbrot_point(complex_t coord_point, 
1252                             complex_t c_constant,
1253                             int Max_iterations, double divergent_limit)
1254 {
1255     complex_t z_point;  /* we need a point to use in our calculation */
1256     int num_iterations;  /* a counter to track the number of iterations done */
1257
1258     num_iterations = 0; /* zero our counter */
1259     z_point = coord_point; /* initialize to the given start coordinate */
1260
1261     /* loop while the absolute value of the complex coordinate is < our limit
1262     (for a mandelbrot) or until we've done our specified maximum number of
1263     iterations (both julia and mandelbrot) */
1264
1265     while (absolute_complex(z_point) < divergent_limit &&
1266         num_iterations < Max_iterations)
1267     {
1268         /* z = z*z + c */
1269         z_point = multiply_complex(z_point,z_point);
1270         z_point = add_complex(z_point,c_constant);
1271
1272         ++num_iterations;
1273     } /* done iterating for one point */
1274
1275     return num_iterations;
1276 }
1277
1278 void output_data(int *in_grid_array, int coord[4], int *out_grid_array, int width, int height)
1279 {
1280     int i, j, k;
1281     int result;
1282
1283     k = 0;
1284     for (j=coord[2]; j<=coord[3]; j++)
1285     {
1286         for (i=coord[0]; i<=coord[1]; i++)
1287         {
1288             out_grid_array[(j * height) + i] = in_grid_array[k];
1289             k++;
1290         }
1291     }
1292     if (!use_stdin)
1293     {
1294         result = MPI_Send(coord, 4, MPI_INT, 0, 0, comm);
1295         if (result != MPI_SUCCESS)
1296         {
1297             printf("Unable to send coordinates to the visualizer, aborting.\n");
1298             MPI_Abort(MPI_COMM_WORLD, -1);
1299         }
1300         result = MPI_Send(in_grid_array, (coord[1] + 1 - coord[0]) * (coord[3] + 1 - coord[2]), MPI_INT, 0, 0, comm);
1301         if (result != MPI_SUCCESS)
1302         {
1303             printf("Unable to send coordinate data to the visualizer, aborting.\n");
1304             MPI_Abort(MPI_COMM_WORLD, -1);
1305         }
1306     }
1307 }
1308
1309 /* Writes out a linear integer array of grayscale pixel values to a
1310 pgm-format file (with header).  The exact format details are given
1311 at the end of this routine in case you don't have the man pages
1312 installed locally.  Essentially, it's a 2-dimensional integer array
1313 of pixel grayscale values.  Very easy to manipulate externally.
1314 */
1315
1316 /* You need the following inputs:
1317 A linear integer array with the actual pixel values (read in as
1318 consecutive rows),
1319 The width and height of the grid, and
1320 The maximum pixel value (to set greyscale range).  We are assuming
1321 that the lowest value is "0".
1322 */
1323
1324 void dumpimage(char *filename, int in_grid_array[], int in_pixels_across, int in_pixels_down,
1325                int in_max_pixel_value, char input_string[], int num_colors, color_t colors[])
1326 {
1327     FILE *ifp;
1328     int i, j, k;
1329 #ifdef USE_PPM
1330     int r, g, b;
1331 #endif
1332
1333     printf("%s\nwidth: %d\nheight: %d\ncolors: %d\nstr: %s\n",
1334         filename, in_pixels_across, in_pixels_down, num_colors, input_string);
1335     fflush(stdout);
1336
1337     if ( (ifp=fopen(filename, "w")) == NULL)
1338     {
1339         printf("Error, could not open output file\n");
1340         MPI_Abort(MPI_COMM_WORLD, -1);
1341         exit(-1);
1342     }
1343
1344 #ifdef USE_PPM
1345     fprintf(ifp, "P3\n");  /* specifies type of file, in this case ppm */
1346     fprintf(ifp, "# %s\n", input_string);  /* an arbitrary file identifier */
1347     /* now give the file size in pixels by pixels */
1348     fprintf(ifp, "%d %d\n", in_pixels_across, in_pixels_down);
1349     /* give the max r,g,b level */
1350     fprintf(ifp, "255\n");
1351
1352     k=0; /* counter for the linear array of the final image */
1353     /* assumes first point is upper left corner (element 0 of array) */
1354
1355     if (in_max_pixel_value < 1)
1356     {
1357         for (j=0; j<in_pixels_down; ++j) /* start at the top row and work down */
1358         {
1359             for (i=0; i<in_pixels_across; ++i) /* go along the row */
1360             {
1361                 fprintf(ifp, "0 0 0 ");
1362             }
1363             fprintf(ifp, "\n"); /* done writing one row, begin next line */
1364         }
1365     }
1366     else
1367     {
1368         for (j=0; j<in_pixels_down; ++j) /* start at the top row and work down */
1369         {
1370             for (i=0; i<in_pixels_across; ++i) /* go along the row */
1371             {
1372                 getRGB(colors[(in_grid_array[k] * num_colors) / in_max_pixel_value], &r, &g, &b);
1373                 fprintf(ifp, "%d %d %d ", r, g, b); /* +1 since 0 = first color */
1374                 ++k;
1375             }
1376             fprintf(ifp, "\n"); /* done writing one row, begin next line */
1377         }
1378     }
1379 #else
1380     fprintf(ifp, "P2\n");  /* specifies type of file, in this case pgm */
1381     fprintf(ifp, "# %s\n", input_string);  /* an arbitrary file identifier */
1382     /* now give the file size in pixels by pixels */
1383     fprintf(ifp, "%d %d\n", in_pixels_across, in_pixels_down);
1384     /* gives max number of grayscale levels */
1385     fprintf(ifp, "%d\n", in_max_pixel_value+1); /* plus 1 because 0=first color */
1386
1387     k=0; /* counter for the linear array of the final image */
1388     /* assumes first point is upper left corner (element 0 of array) */
1389
1390     for (j=0;j<in_pixels_down;++j) /* start at the top row and work down */
1391     {
1392         for (i=0;i<in_pixels_across;++i) /* go along the row */
1393         {
1394             fprintf(ifp, "%d ", in_grid_array[k]+1); /* +1 since 0 = first color */
1395             ++k;
1396         }
1397         fprintf(ifp, "\n"); /* done writing one row, begin next line */
1398     }
1399 #endif
1400
1401     fclose(ifp);
1402 }