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