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