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