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