src/mandelbrot_MPE/MS3_MPE.c

   1
   2
   3                         /*   MS3  */
   4
   5                /* Adapted for MPE display  */
   6
   7 #include <math.h>
   8 #include <stdio.h>
   9 #include <stdlib.h>
  10 #include <string.h>
  11 #include <mpi.h>
  12 #include <mpe.h>
  13 #include <mpe_graphics.h>
  14
  15
  16                 /********************************/
  17                 /*                              */
  18                 /*     The Mandelbrot Set       */
  19                 /*                              */
  20                 /*      z = z*z + c             */
  21                 /*                              */
  22                 /********************************/
  23
  24 /*  Simplified main() via two functions   */
  25 /*  iterate_map() and set_grid()          */
  26
  27 /*  Slice image - strips parallel to real C */
  28 /*  via functions                           */
  29 /*  set_slices() and calculate_slice()      */
  30
  31 /*  Parallel version -- Master-Slave model  */
  32 /*                                          */
  33 /*  master()  running on node zero          */
  34 /*  slave()   running on all other nodes    */
  35 /*                                          */
  36
  37
  38 /* use:                  */
  39 /*                       */
  40 /* MS3                   */
  41
  42
  43                                         /* Michel Vallieres */
  44
  45
  46                                         /* large distance; if |z| goes  */
  47                                         /* beyond it, point C is not a  */
  48                                         /* member of the set            */
  49 #define LARGE_RADIUS 2.0
  50                                         /* Maximum number of iterations */
  51                                         /* before declaring a point in  */
  52                                         /* the Mandelbrot set           */
  53 #define MAX_ITERATIONS 120
  54                                         /* image size */
  55 #define NX 700
  56 #define NY 500
  57                                         /* # of slices parallel to C real */
  58 #define N_slices 64
  59                                         /* Maximum # processes */
  60 #define MAX_PROCESSES 1000
  61                                         /* # of colors in MPE palette */
  62 #define NC 128
  63
  64
  65
  66 /* Iterate the Mandelbrot map -- decide if the (input) point  */
  67 /* in complex C plane belongs or not to the set               */
  68 /* prints out the number of iterations pro-rated to           */
  69 /* MAX_ITERATIONS                                             */
  70 double iterate_map( double c_real, double c_img, int *count )
  71 {
  72    double z_real, z_img, z_current_real, z_magnitude;
  73    double color;
  74    int    counter;
  75
  76                                         /* z starts a origin */
  77    z_real = 0.0;
  78    z_img  = 0.0;
  79                                         /* set counter to zero */
  80    counter = 0;
  81                                         /* iterate map for             */
  82                                         /*  MAX_ITERATIONS times or ...*/
  83    while ( counter < MAX_ITERATIONS )
  84    {
  85      z_current_real = z_real;
  86      z_real = z_real*z_real - z_img*z_img + c_real;
  87      z_img  = 2.0*z_current_real*z_img + c_img;
  88      counter = counter + 1;
  89                                         /* ... until magnitude of z */
  90                                         /* is larger than some      */
  91                                         /* large radius             */
  92      z_magnitude = z_real*z_real+z_img*z_img;
  93      if ( z_magnitude > LARGE_RADIUS ) break;
  94    }
  95                                         /* write C real and       */
  96                                         /* imaginary parts if     */
  97                                         /* C belongs to the set   */
  98          /*
  99            if ( z_magnitude < LARGE_RADIUS )
 100            printf ( " %f %f  \n ", c_real, c_img );
 101          */
 102                                         /* or iteration number      */
 103                                         /* for all scanned C values */
 104                                         /* for color picture        */
 105    count++;
 106
 107    color = (double)(counter) / (double)MAX_ITERATIONS;
 108    return color;
 109
 110 }
 111
 112
 113 /* Set up the grid in complex C plane to generate image  */
 114 void set_grid( double *crmin, double *crmax, double *cimin,
 115                double *cimax, double *dcr, double *dci )
 116 {
 117                                         /* define grid spanning   */
 118                                         /* the complex plane for  */
 119                                         /* the C variable         */
 120    *crmin = -1.7;
 121    *crmax = 0.8;
 122    *cimin = -1.1;
 123    *cimax = 1.1;
 124                                         /* C spacing on complex plane */
 125    *dcr = (*crmax - *crmin)/(NX-1);
 126    *dci = (*cimax - *cimin)/(NY-1);
 127 }
 128
 129
 130 /* Calculate the Mandelbrot set in slice  */
 131 void calculate_slice( int slice, int *begin_slice,
 132             int *end_slice, int *count,
 133             double crmin, double cimin, double dcr, double dci,
 134                       double *store_slice, MPE_Point store_points[] )
 135 {
 136   int    i, j, k;
 137   double c_real, c_img, color;
 138   int    color_pointer;
 139                                         /* store pointer */
 140   k = 0;
 141                                         /* scan over slice */
 142   for ( j=begin_slice[slice] ; j<=end_slice[slice] ; j++ )
 143   {
 144     for ( i=0 ; i<NX ; i++ )
 145     {
 146        c_real = crmin + i*dcr;
 147        c_img = cimin + j*dci;
 148                                         /* iterate Mandelbrot map */
 149        color = iterate_map( c_real, c_img, count );
 150
 151                                         /* store iteration # for */
 152                                         /* use in MPE display    */
 153                                         /* color is NC scaled by */
 154                                         /* iterations/max_iterations */
 155        store_points[k].x = i;
 156        store_points[k].y = j;
 157        color_pointer = NC*color;
 158        store_points[k].c = color_pointer;
 159
 160                                         /* store resulting color */
 161        store_slice[k] = color;
 162        k++;
 163     }
 164   }
 165 }
 166
 167
 168 /* Setup the slices through which compute the Mandelbrot set  */
 169 int set_slices( int *begin_slice, int *end_slice, int *offset )
 170 {
 171   int slice, average, leftover, temp, temp_offset;
 172   MPI_Status status;
 173
 174   average = NY / N_slices;
 175   leftover = NY % N_slices;
 176   temp = -1;
 177   temp_offset = 0;
 178   for ( slice=0 ; slice<N_slices ; slice++ )
 179     {
 180       begin_slice[slice] = temp + 1;
 181       end_slice[slice] = begin_slice[slice] + average - 1;
 182       if ( leftover >0 )
 183         {
 184           end_slice[slice]++;
 185           leftover--;
 186         }
 187       temp = end_slice[slice];
 188       offset[slice] = temp_offset;
 189       temp_offset = temp_offset + NX *
 190               ( end_slice[slice] - begin_slice[slice] + 1 );
 191     }
 192   return average;
 193 }
 194
 195
 196 /* terminates graph interactively */
 197 int quit_by_touch( MPE_XGraph graph, int WINDOW_SIZE_X,
 198                    int WINDOW_SIZE_Y  )
 199 {
 200   int x, y, button, xc, yc, w, h, xs, ys;
 201
 202   xc = 0.95*WINDOW_SIZE_X;              /* draw filled rectangle     */
 203   yc = 0.95*WINDOW_SIZE_Y;              /* lower left bottom         */
 204   w = 0.05*WINDOW_SIZE_X;               /* corner of graphic window  */
 205   h = WINDOW_SIZE_Y - yc;
 206                                         /* write Q in small rectangle */
 207   MPE_Fill_rectangle( graph, xc, yc, w, h, MPE_BLACK );
 208   xs = xc + 0.4*w;
 209   ys = yc + 0.8*h;
 210   MPE_Draw_string( graph, xs, ys, MPE_WHITE, "Q" );
 211
 212                                         /* infinite loop            */
 213                                         /* wait for different touch */
 214   for ( ; ; )
 215     {
 216                                         /* wait for mouse to touch graph */
 217                                         /* can only be done by master    */
 218       MPE_Get_mouse_press( graph, &x, &y, &button );
 219                                         /* result is touch point */
 220       fprintf( stderr, "\n Graph coordinates of the touch: %d %d %d \n",
 221                        x, y, button );
 222       if ( x < xc+w && y > yc )
 223         {
 224           fprintf( stderr,
 225                    "\n Last touch was in the -- quit -- rectangle \n\n" );
 226           return;
 227         }
 228     }
 229 }
 230
 231
 232
 233 /* slave processes -- calculate a specific slice   */
 234 /* specified by "slice" as received from process 0 */
 235 void slave( int my_rank, int *begin_slice, int *end_slice,
 236             int average, double crmin, double cimin, double dcr,
 237             double dci, double *store_slice , MPE_Point store_points[],
 238             MPE_XGraph graph )
 239 {
 240   static int count;
 241   int        number, store_n, slice;
 242   int        done;
 243   MPI_Status status;
 244                                          /* signal ->  slice is done */
 245   done = 1;
 246
 247   for ( ; ; )
 248     {
 249                                         /* a new slice to calculate */
 250       MPI_Recv( &slice, 1, MPI_INT, 0, 325, MPI_COMM_WORLD, &status );
 251
 252                                         /* suicide signal */
 253       if ( slice < 0 )
 254         {
 255           free( store_slice );
 256           free( store_points );
 257           MPE_Close_graphics( &graph );
 258           MPI_Finalize();
 259           exit(0);
 260         }
 261                                        /* calculate requested slice */
 262       calculate_slice( slice, begin_slice, end_slice, &count,
 263                        crmin, cimin, dcr, dci, store_slice,
 264                        store_points );
 265
 266                                        /* send results back to master */
 267       //  number = NX * (average+1) ;
 268       //  MPI_Ssend( store_slice, number, MPI_DOUBLE, 0, 327,
 269       //                MPI_COMM_WORLD );
 270
 271                                        /* send signal to node 0 */
 272       MPI_Ssend( &done, 1, MPI_INT, 0, 327, MPI_COMM_WORLD );
 273
 274                                        /* draw the slice in MPE graph */
 275       store_n = NX * ( end_slice[slice]-begin_slice[slice]+1 );
 276       MPE_Draw_points( graph, store_points, store_n);
 277
 278                                         /* update the graph */
 279       MPE_Update( graph );
 280
 281       fprintf( stderr, "slave %d  ->  slice %d calculated & drawn\n",
 282                             my_rank, slice );
 283     }
 284 }
 285
 286
 287
 288 /* Master node -- control the calculation */
 289 void master(int numprocs, int *begin_slice, int *end_slice,
 290             int average, double crmin, double cimin, double dcr,
 291             double dci, double *store_slice, int *offset,
 292             MPE_XGraph graph )
 293 {
 294   int    k, k_slice, slice, process;
 295   int    number, source;
 296   int    done;
 297   static int count;
 298   static int list_slices[MAX_PROCESSES];
 299   int    kill, recv_slice;
 300   double mandelbrot[NX*NY];
 301   MPI_Status status;
 302                                         /* # points in set */
 303    count = 0;
 304                                         /* scan over slices */
 305    slice = 0;
 306                                         /* send a slice to each slave */
 307    for ( process=1 ; process<numprocs ; process++ )
 308      {
 309                                         /* generate Mandelbrot set */
 310                                         /*  in each process        */
 311        MPI_Ssend( &slice, 1, MPI_INT, process, 325,
 312                 MPI_COMM_WORLD );
 313        list_slices[process] = slice;
 314        slice++;
 315      }
 316
 317                                         /* scan over received slices */
 318    for ( k_slice=0 ; k_slice<N_slices ; k_slice++ )
 319      {
 320                                         /* receive a slice */
 321        // number = NX * ( average+1);
 322        // MPI_Recv( store_slice, number, MPI_DOUBLE, MPI_ANY_SOURCE, 327,
 323        //               MPI_COMM_WORLD, &status );
 324
 325                                         /* watch for done signal */
 326        MPI_Recv( &done, 1, MPI_INT, MPI_ANY_SOURCE, 327,
 327                         MPI_COMM_WORLD, &status );
 328
 329                                         /* source of slice */
 330        source = status.MPI_SOURCE;
 331                                         /* received slice */
 332        recv_slice = list_slices[source];
 333                                         /* send a slice to this slave */
 334        if ( slice < N_slices )
 335          {
 336             MPI_Ssend( &slice, 1, MPI_INT, source, 325,
 337                                         MPI_COMM_WORLD );
 338             list_slices[source] = slice;
 339             slice++;
 340          }
 341                                         /* actual number */
 342        // number = NX * ( end_slice[recv_slice]-begin_slice[recv_slice]+1 );
 343
 344                                         /* store set in matrix */
 345        // for ( k=0 ; k<number ; k++ )
 346        //   mandelbrot[offset[recv_slice]+k] = store_slice[k];
 347
 348      }
 349
 350                                         /* ouput Mandelbrot set */
 351    // fprintf( stderr, " The Mandelbrot Set will be written out \n" );
 352    // for ( k=0; k<NX*NY ; k++ )
 353    //        printf( " %f \n", mandelbrot[k] );
 354
 355    // fprintf( stderr, " The Mandelbrot Set has been written out \n" );
 356
 357                                         /* calculation is done!    */
 358                                         /* all stripes done        */
 359                                         /* wait for a touch - quit */
 360    sleep( 2 );
 361    quit_by_touch( graph, NX, NY );
 362
 363                                         /* end all slave processes */
 364    kill = -1;
 365    for ( process=1 ; process<numprocs ; process++ )
 366      {
 367        MPI_Ssend( &kill, 1, MPI_INT, process, 325, MPI_COMM_WORLD );
 368      }
 369                                         /* end itself */
 370    free( store_slice );
 371    MPE_Close_graphics( &graph );
 372    MPI_Finalize();
 373    exit( 0 );
 374 }
 375
 376
 377
 378 int main( int argc, char *argv[] )
 379 {
 380    double crmin, crmax, cimin, cimax;
 381    double dcr, dci;
 382    int    average;
 383    int    begin_slice[1000], end_slice[1000], offset[1000];
 384    int    my_rank, numprocs;
 385    double *store_slice;
 386    MPE_Point *store_points;
 387    MPE_Color colors[500];
 388    MPE_XGraph graph;
 389    int    ierr;
 390
 391    char displayname[MPI_MAX_PROCESSOR_NAME+4] = "";
 392
 393                             /* MPI exposure  */
 394    MPI_Init( &argc, &argv );
 395    MPI_Comm_size( MPI_COMM_WORLD, &numprocs );
 396    MPI_Comm_rank( MPI_COMM_WORLD, &my_rank );
 397
 398    if ( my_rank == 0 )
 399         strcpy( displayname, getenv( "DISPLAY" ) );
 400
 401    MPI_Bcast( displayname, MPI_MAX_PROCESSOR_NAME+4, MPI_CHAR,
 402                0, MPI_COMM_WORLD );
 403    fprintf( stdout, "%d : $DISPLAY at process 0 = %s\n",
 404              my_rank, displayname );
 405    fflush( stdout );
 406
 407    ierr = MPE_Open_graphics( &graph, MPI_COMM_WORLD, displayname,
 408                               -1, -1, 700, 500, 0 );
 409
 410    if ( ierr != MPE_SUCCESS )
 411      {
 412        fprintf( stderr, "%d : MPE_Open_graphics() fails\n", my_rank );
 413        ierr = MPI_Abort( MPI_COMM_WORLD, 1 );
 414        exit(1);
 415      }
 416
 417    fprintf( stderr, "After Open Graphics call \n" );
 418
 419                                         /* make a new color palette */
 420    MPE_Make_color_array( graph, NC, colors );
 421
 422                                         /* set slicing */
 423    average = set_slices( begin_slice, end_slice, offset );
 424
 425                                         /* setup C grid */
 426    set_grid( &crmin, &crmax, &cimin, &cimax, &dcr, &dci );
 427
 428                                         /* color storage */
 429    store_slice = (double *)malloc( NX*(average+1)*sizeof(double) );
 430
 431                                         /* color storage for MPE graph */
 432    store_points = (MPE_Point *)malloc( NX*(average+1)*sizeof(MPE_Point) );
 433
 434                                         /* master node */
 435    if ( my_rank == 0 )
 436      {
 437        master( numprocs, begin_slice, end_slice, average,
 438                crmin, cimin, dcr, dci, store_slice, offset, graph );
 439      }
 440                                         /* slave nodes  */
 441    else
 442      {
 443        slave( my_rank, begin_slice, end_slice, average,
 444               crmin, cimin, dcr, dci, store_slice, store_points, graph );
 445      }
 446
 447 }
 448
 449