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