content/programming_strategies/src/legacy_codes/mandelbrot_3.c

   1
   2                 /**********************************/
   3                 /*                                */
   4                 /*     The Mandelbrot Set         */
   5                 /*                                */
   6                 /*      z = z*z + c               */
   7                 /*                                */
   8                 /* Adapted to calculate the Set   */
   9                 /* in stripes parallel to Img axis*/
  10                 /*                                */
  11                 /* Improved routines structure    */
  12                 /* to facilitate parallel         */
  13                 /*     implementation             */
  14                 /**********************************/
  15
  16                                         /* Michel Vallieres */
  17
  18 #include <stdio.h>
  19 #include <time.h>
  20 #include <mpi.h>
  21                                         /* large distance; if |z| reaches */
  22                                         /* it, point C is not a member    */
  23                                         /* of the set                     */
  24 #define LARGE_RADIUS 2.0
  25                                         /* Maximum number of iterations */
  26                                         /* before declaring a point in  */
  27                                         /* the Mandelbrot set           */
  28 #define MAX_ITERATIONS 80
  29 #define NX 600                          /* image size */
  30 #define NY 500
  31 #define N_STRIPE 37                     /* number of stripes */
  32 #define MAX_PROCS 100
  33
  34                                         /* prototypes   */
  35 void compute_stripe( double crmin, double crmax, double cimin, double cimax,
  36                      int i_min, int i_max, double store_iter[],
  37                      int *store_n);
  38 void grid_stripes( double *crmin, double *crmax, double *cimin, double *cimax,
  39                       int begin_stripe[], int last_stripe[] );
  40
  41
  42 int main( int argc, char *argv[] )
  43 {
  44   double crmin, crmax, cimin, cimax;
  45   double store_iter[NY*(NX/N_STRIPE+1)];
  46   int begin_stripe[N_STRIPE], last_stripe[N_STRIPE];
  47   int stripe;
  48   int store_n, i_point;
  49   int count;
  50   int my_id, numprocs, ip;
  51   int sent_stripe, recv_stripe, stripe_in_proc[MAX_PROCS];
  52   int from_node, number_from_node;
  53   int end_signal;
  54
  55   MPI_Status  mesgStatus;
  56
  57   MPI_Init(&argc, &argv);               /* start MPI virtual machine */
  58
  59   MPI_Comm_rank(MPI_COMM_WORLD, &my_id);
  60   MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
  61
  62
  63                                         /* initialize grid & stripes */
  64   grid_stripes ( &crmin, &crmax, &cimin, &cimax, begin_stripe, last_stripe );
  65
  66   count = 0;
  67                                         /**********/
  68                                         /* master */
  69                                         /**********/
  70   if(my_id == 0)
  71   {
  72      sent_stripe = 0;
  73      recv_stripe = 0;                   /* logic of stripes */
  74                                         /* send original stripe to slaves */
  75      for( ip = 1; ip < numprocs; ip++)
  76      {
  77         MPI_Ssend(&sent_stripe, 1, MPI_INT, ip, 5, MPI_COMM_WORLD);
  78         stripe_in_proc[ip] = sent_stripe;
  79         sent_stripe++;
  80      }
  81
  82      for( recv_stripe = 0; recv_stripe < N_STRIPE; recv_stripe++)
  83      {
  84                                         /* Recieve result back from a node */
  85         MPI_Recv(store_iter, NY*(NX/N_STRIPE+1), MPI_DOUBLE, MPI_ANY_SOURCE, 7,
  86                                          MPI_COMM_WORLD, &mesgStatus);
  87         from_node = mesgStatus.MPI_SOURCE;
  88         MPI_Get_count( &mesgStatus, MPI_DOUBLE, &number_from_node );
  89         fprintf( stderr, "From node %d -- stripe & # %d %d \n", from_node,
  90                              stripe_in_proc[from_node], number_from_node );
  91                                         /* Termination logic. */
  92         end_signal = -1;
  93         fprintf( stderr, "I'm sending a termination signal to %d \n", from_node);
  94         fflush(stderr);
  95         MPI_Ssend( &end_signal, 1, MPI_INT, from_node, 5, MPI_COMM_WORLD);
  96         fprintf( stderr, "I have sent a termination signal to %d \n", from_node);
  97         fflush(stderr);
  98      }
  99
 100   }
 101                                         /**********/
 102                                         /* slave  */
 103                                         /**********/
 104   else
 105   {
 106      for ( ; ; )
 107        {
 108                                         /* I'm getting a stripe to do, YAHOO! */
 109         MPI_Recv(&stripe, 1, MPI_INT, 0, 5, MPI_COMM_WORLD, &mesgStatus);
 110         fprintf(stderr, "my ID: %d -- recieved stripe: %d \n", my_id, stripe);
 111
 112         if( stripe < 0)
 113         {
 114           fprintf(stderr, "Node %d is quitting. \n", my_id);
 115           MPI_Finalize();
 116           exit(0);
 117         }
 118
 119                                         /* find the points in the set */
 120         compute_stripe( crmin, crmax, cimin, cimax, begin_stripe[stripe],
 121                      last_stripe[stripe], store_iter, &store_n);
 122
 123                                         /* send stripe results back to master */
 124        MPI_Ssend(store_iter, store_n, MPI_DOUBLE, 0, 7, MPI_COMM_WORLD);
 125        }
 126  }
 127
 128   MPI_Finalize();
 129 }
 130
 131
 132
 133
 134 void grid_stripes( double *crmin, double *crmax, double *cimin, double *cimax,
 135                       int begin_stripe[], int last_stripe[] )
 136 {
 137   int  w_stripe, stripe, left_over;
 138                                         /* define grid spanning   */
 139                                         /* the complex plane for  */
 140                                         /* the C variable         */
 141   *crmin=-1.7;
 142   *crmax = 0.8;
 143   *cimin = -1.0;
 144   *cimax = 1.0;
 145                                         /* stripe definition            */
 146                                         /* stripe # does not have to be */
 147                                         /* commensurate with grid size  */
 148                                         /* in which case add 1 to some  */
 149                                         /* stripes */
 150   w_stripe= NX/N_STRIPE;
 151   left_over = NX % N_STRIPE;
 152                                         /* fist stripe */
 153   begin_stripe[0] = 0;
 154   last_stripe[0] = begin_stripe[0] + w_stripe - 1;
 155   if ( left_over > 0 )
 156     {
 157       last_stripe[0] = last_stripe[0] + 1;
 158       left_over--;
 159     }
 160                                         /* other stripes */
 161   for ( stripe = 1; stripe < N_STRIPE ; stripe++ )
 162     {
 163       begin_stripe[stripe] = last_stripe[stripe-1] + 1 ;
 164       last_stripe[stripe] = begin_stripe[stripe] + w_stripe - 1;
 165       if ( left_over > 0 )
 166         {
 167       last_stripe[stripe] = last_stripe[stripe] + 1;
 168       left_over--;
 169         }
 170     }
 171 }
 172
 173
 174
 175
 176 void compute_stripe( double crmin, double crmax, double cimin, double cimax,
 177                      int i_min, int i_max, double store_iter[],
 178                      int *store_n)
 179 {
 180    double z_real, z_img, z_magnitude;
 181    double c_real, c_img, z_current_real;
 182    double dcr, dci;
 183    int    i, j;
 184    int    counter;
 185
 186                                         /* scan C values on complex plane */
 187    dcr = (crmax - crmin)/(double)(NX-1);
 188    dci = (cimax - cimin)/(double)(NY-1);
 189                                         /* scan stripe */
 190    *store_n=0;
 191    for ( i=i_min ; i<=i_max; i++ )
 192    {
 193       for ( j=0 ; j<NY; j++ )
 194       {
 195            c_real = crmin + i*dcr;
 196            c_img = cimin + j*dci;
 197                                         /* z starts a origin */
 198            z_real = 0.0;
 199            z_img  = 0.0;
 200                                         /* set counter to zero */
 201            counter = 0;
 202                                         /* iterate map for             */
 203                                         /*  MAX_ITERATIONS times or ...*/
 204            while ( counter < MAX_ITERATIONS )
 205            {
 206              z_current_real = z_real;
 207              z_real = z_real*z_real - z_img*z_img + c_real;
 208              z_img  = 2.0*z_current_real*z_img + c_img;
 209              counter = counter + 1;
 210                                         /* ... until magnitude of z */
 211                                         /* is larger than some      */
 212                                         /* large radius             */
 213              z_magnitude = z_real*z_real+z_img*z_img;
 214              if ( z_magnitude > LARGE_RADIUS ) break;
 215            }
 216                                         /* store iteration numbers */
 217                                         /* to plot the set         */
 218            store_iter[*store_n] =
 219                     (double)(255*counter) / (double)MAX_ITERATIONS ;
 220            (*store_n)++;
 221        }
 222     }
 223 }
 224