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