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