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