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