git.gag.com Git - debian/atlc/blob - src/non_gui/finite_difference_mpi.c

   1 /* atlc - arbitrary transmission line calculator, for the analysis of
   2 transmission lines are directional couplers.
   3
   4 Copyright (C) 2002. Dr. David Kirkby, PhD (G8WRB).
   5
   6 This program is free software; you can redistribute it and/or
   7 modify it under the terms of the GNU General Public License
   8 as published by the Free Software Foundation; either package_version 2
   9 of the License, or (at your option) any later package_version.
  10
  11 This program is distributed in the hope that it will be useful,
  12 but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 GNU General Public License for more details.
  15
  16 You should have received a copy of the GNU General Public License
  17 along with this program; if not, write to the Free Software
  18 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307,
  19 USA.
  20
  21 Dr. David Kirkby, e-mail drkirkby at ntlworld.com
  22
  23 */
  24 #include "config.h"
  25
  26 #ifdef ENABLE_MPI /* file only needed on MPI systems. */
  27
  28 #include <mpi.h>
  29
  30 #ifdef HAVE_STDLIB_H
  31 #include <stdlib.h>
  32 #endif
  33
  34 #include "definitions.h"
  35
  36
  37 extern int coupler;
  38 extern int width, height;
  39 extern double **Vij, **Er;
  40 extern signed char **oddity;
  41 extern int num_pes;
  42 extern int dielectrics_to_consider_just_now;
  43 extern struct strip strip_map[MAX_PES+1];
  44
  45 /* these variables are only needed on PE 0,
  46    but are declared here for all for convenience.
  47    by switching to a dynamic allocation scheme,
  48    we could make them local to finite_difference,
  49    and only allocate them on PE 0.
  50 */
  51 MPI_Request energy_rcv_requests[MAX_PES];
  52 double energies[MAX_PES];
  53
  54 /*
  55   the job of a worker PE is to work on a columnar strip
  56   of the voltage array for a given number of iterations,
  57   then send off the results to PE 0.
  58
  59   in order for a worker PE to calculate the edges of
  60   its strip, it needs two additional columns of data
  61   to the left and two additional columns of data to
  62   the right of its strip.  initially, PE 0 supplies
  63   the extra columns.  after the first iteration, the
  64   the PE needs to update the columns adjacent to its
  65   strip edges, and to do that, it needs the columns
  66   that are adjacent to those (hence each worker PE
  67   maintains 4 columns of data in addtion to its strip).
  68   the outermost columns are sent from the neighbor PEs
  69   after they are done computing them for the current
  70   iteration.  after they are received, the columns
  71   adjacent to the strip edges are updated, and the PE
  72   moves on to the next iteration.
  73
  74   after the prescribed number of iterations are completed,
  75   each worker PE sends its strip of the voltage matrix
  76   to PE 0, so that PE 0 has a complete up-to-date copy
  77   of the voltage array after each call to finite_difference.
  78   while that data is being sent, the worker PE computes
  79   the energy_per_metre ot its strip and then sends that
  80   off to PE 0, so that PE 0 can compute the overall
  81   capacitance_per_metre.
  82
  83 */
  84
  85
  86 /*
  87   this routine is only run on the worker PEs
  88 */
  89
  90 void mpi_worker(int rank) {
  91
  92   int width_height[2];
  93   int strip_params[2];
  94   int start_col, num_cols;
  95   MPI_Status status;
  96   int control, done;
  97   int i,j,iterations;
  98   MPI_Request send_requests[2];
  99   MPI_Status send_statuses[2];
 100   MPI_Request rcv_requests[2];
 101   double energy_per_metre;
 102   int index;
 103
 104   /* get the total width and height of the voltage
 105          matrix.  the worker PE needs to know the column
 106          height in order to run its calculations.  the
 107          width is also currently needed for the electric
 108          field subroutines.
 109   */
 110   MPI_Recv(&width_height,
 111                    2,
 112                    MPI_INT,
 113                    0,
 114                    MSG_TAG_WIDTH_HEIGHT,
 115                    MPI_COMM_WORLD,
 116                    &status);
 117
 118   width = width_height[0];
 119   height = width_height[1];
 120
 121   /* get the location and size of the strip of the
 122          voltage matrix that has been assigned to this
 123          PE.  strictly speaking, the PE does not need
 124          to know the starting column number since it
 125          uses its own local indexing scheme, but it
 126          is sent anyway as a debugging aid if needed.
 127   */
 128   MPI_Recv(&strip_params,
 129                    sizeof(strip_params)/sizeof(int),
 130                    MPI_INT,
 131                    0,
 132                    MSG_TAG_STRIP_PARAMS,
 133                    MPI_COMM_WORLD,
 134                    &status);
 135
 136   /* this is the starting column in the global voltage matrix
 137          of the columnar strip this PE has been assigned*/
 138   start_col = strip_params[0];
 139
 140   /* this is the width of the columnar strip that this PE
 141          has been assigned from the global voltage matrix */
 142   num_cols = strip_params[1];
 143
 144   /* allocate matrixes big enough to contain the
 145          assigned strip and supporting data */
 146   oddity=cmatrix(0,num_cols+4,0,height);
 147   Vij=dmatrix(0,num_cols+4,0,height);
 148   Er=dmatrix(0,num_cols+4,0,height);
 149
 150   /* get the oddity data to use in computing
 151          the assigned strip */
 152   MPI_Recv(oddity[0],
 153                    (num_cols+4)*height,
 154                    MPI_DOUBLE,
 155                    0,
 156                    MSG_TAG_NODE_TYPE,
 157                    MPI_COMM_WORLD,
 158                    &status);
 159
 160   /* get the Er data to use in computing
 161          the assigned strip */
 162   MPI_Recv(Er[0],
 163                    (num_cols+4)*height,
 164                    MPI_DOUBLE,
 165                    0,
 166                    MSG_TAG_ER,
 167                    MPI_COMM_WORLD,
 168                    &status);
 169
 170   /*************************************************
 171    * all of the data received above this point is  *
 172    * sent only once in the lifetime of the program *
 173    *************************************************/
 174
 175   done = 0;
 176   do {
 177
 178         /* recieve a control word that tells
 179            the PE whether to set off on another
 180            set of iterations, or to exit */
 181
 182         MPI_Recv(&control,
 183                          1,
 184                          MPI_INT,
 185                          0,
 186                          MSG_TAG_CONTROL,
 187                          MPI_COMM_WORLD,
 188                          &status);
 189
 190         switch (control) {
 191         case CONTROL_VALUE_RECEIVE_DATA:
 192           /* receive the strip of the voltage matrix
 193                  that we are to update. this is sent by PE 0. */
 194           MPI_Recv(Vij[1],
 195                            (num_cols+2)*height,
 196                            MPI_DOUBLE,
 197                            0,
 198                            MSG_TAG_VIJ,
 199                            MPI_COMM_WORLD,
 200                            &status);
 201
 202           /* receive the current value of
 203                  dielectrics_to_consider_just_now.
 204                  this is sent by PE 0. */
 205           MPI_Recv(&dielectrics_to_consider_just_now,
 206                            1,
 207                            MPI_INT,
 208                            0,
 209                            MSG_TAG_ORDINARY_INTERIOR_POINTS,
 210                            MPI_COMM_WORLD,
 211                            &status);
 212
 213           break;
 214         case CONTROL_VALUE_SEND_DATA:
 215           /* send our strip to PE 0 */
 216           MPI_Send(Vij[2],
 217                            num_cols*height,
 218                            MPI_DOUBLE,
 219                            0,
 220                            MSG_TAG_VIJ,
 221                            MPI_COMM_WORLD);
 222           break;
 223         case CONTROL_VALUE_DO_ITERATIONS:
 224           /* receive the number of iterations we are
 225                  to compute for. this is sent by PE 0 at
 226                  the beginning of finite_difference. */
 227           MPI_Recv(&iterations,
 228                            1,
 229                            MPI_INT,
 230                            0,
 231                            MSG_TAG_ITERATIONS,
 232                            MPI_COMM_WORLD,
 233                            &status);
 234
 235           i=0;
 236           do {
 237
 238                 /* update our strip of the voltage matrix */
 239                 do_columns(2, num_cols, 0);
 240
 241                 /* send the columns that the neighbor PEs
 242                    require to the nieghbor PEs */
 243                 MPI_Isend(Vij[num_cols+1],
 244                                   height,
 245                                   MPI_DOUBLE,
 246                                   (rank+1)%num_pes,
 247                                   MSG_TAG_VIJ_RBORDER,
 248                                   MPI_COMM_WORLD,
 249                                   &send_requests[1]);
 250
 251                 MPI_Isend(Vij[2],
 252                                   height,
 253                                   MPI_DOUBLE,
 254                                   rank-1,
 255                                   MSG_TAG_VIJ_LBORDER,
 256                                   MPI_COMM_WORLD,
 257                                   &send_requests[0]);
 258
 259                 /* receive the columns that we need
 260                    to update the columns adjacent
 261                    to our strip edges from the neighbor
 262                    PEs */
 263                 MPI_Irecv(Vij[num_cols + 3],
 264                                   height,
 265                                   MPI_DOUBLE,
 266                                   (rank+1)%num_pes,
 267                                   MSG_TAG_VIJ_LBORDER,
 268                                   MPI_COMM_WORLD,
 269                                   &rcv_requests[0]);
 270
 271                 MPI_Irecv(Vij[0],
 272                                   height,
 273                                   MPI_DOUBLE,
 274                                   rank-1,
 275                                   MSG_TAG_VIJ_RBORDER,
 276                                   MPI_COMM_WORLD,
 277                                   &rcv_requests[1]);
 278
 279                 /* update the columns adjacent to our strip
 280                    edges */
 281                 MPI_Waitany(2,
 282                                         rcv_requests,
 283                                         &index,
 284                                         &status);
 285
 286                 if (0 == index) {
 287                   update_voltage_array(num_cols + 2,0,Vij_TO_Vij);
 288                 } else {
 289                   update_voltage_array(1,0,Vij_TO_Vij);
 290                 }
 291
 292                 MPI_Waitany(2,
 293                                         rcv_requests,
 294                                         &index,
 295                                         &status);
 296
 297                 if (0 == index) {
 298                   update_voltage_array(num_cols + 2,0,Vij_TO_Vij);
 299                 } else {
 300                   update_voltage_array(1,0,Vij_TO_Vij);
 301                 }
 302
 303
 304                 MPI_Waitall(2, send_requests, send_statuses);
 305
 306                 i++;
 307           } while(i < iterations);
 308
 309
 310           energy_per_metre=0.0;
 311           for(i=2;i<2+num_cols;++i)
 312                 for(j=0;j<height;++j) {
 313                   energy_per_metre+=find_energy_per_metre(i,j);
 314                 }
 315
 316           MPI_Send(&energy_per_metre,
 317                            1,
 318                            MPI_DOUBLE,
 319                            0,
 320                            MSG_TAG_ENERGY,
 321                            MPI_COMM_WORLD);
 322           break;
 323         case CONTROL_VALUE_EXIT:
 324           done = 1;
 325           break;
 326         }
 327
 328   } while (0 == done);
 329
 330 }
 331
 332
 333
 334 void do_columns(int start_col, int num_cols, int calculate_edges)
 335 {
 336   int i;
 337   int end_col;
 338
 339   end_col = start_col + num_cols - 1;
 340
 341   for(i=start_col; i<=end_col; ++i)
 342   {
 343         update_voltage_array(i, calculate_edges,Vij_TO_Vij);
 344   }
 345 }
 346
 347
 348
 349 /* There are three versions of finite_difference. They all take the same
 350 arguments and return the same data. One is single threaded, one is
 351 multi-threaded, and this one uses MPI. */
 352
 353 double finite_difference(int accuracy)
 354 {
 355   int i, j, a;
 356   double capacitance_per_metre=0.0, energy_per_metre;
 357   int left_start_col, left_num_cols;
 358   int right_start_col, right_num_cols;
 359   MPI_Status status;
 360   int control;
 361   MPI_Request send_requests[2];
 362   MPI_Status send_statuses[2];
 363   MPI_Request rcv_requests[2];
 364   MPI_Status rcv_statuses[2];
 365   int index;
 366
 367   for(i=1; i<num_pes; i++) {
 368
 369         /* arm worker */
 370         control = CONTROL_VALUE_DO_ITERATIONS;
 371         MPI_Send(&control,
 372                          1,
 373                          MPI_INT,
 374                          i,
 375                          MSG_TAG_CONTROL,
 376                          MPI_COMM_WORLD);
 377
 378         /* send the iteration count */
 379         MPI_Send(&accuracy,
 380                          1,
 381                          MPI_INT,
 382                          i,
 383                          MSG_TAG_ITERATIONS,
 384                          MPI_COMM_WORLD);
 385   }
 386
 387   /* these are the parameters of the left and right
 388          half-strips that PE 0 is responsible for */
 389   left_start_col = strip_map[0].start_col;
 390   left_num_cols = strip_map[0].num_cols;
 391
 392   right_start_col = strip_map[num_pes].start_col;
 393   right_num_cols = strip_map[num_pes].num_cols;
 394
 395   for(a=1; a<=accuracy; ++a) {
 396
 397         /* update the left and right half-strips */
 398         do_columns(left_start_col+1, left_num_cols-1, 1);
 399         do_columns(right_start_col, right_num_cols-1, 1);
 400
 401         /* send the columns that PE 1 and PE N-1
 402            need to those PEs */
 403         MPI_Isend(Vij[left_num_cols - 2],
 404                           height,
 405                           MPI_DOUBLE,
 406                           1,
 407                           MSG_TAG_VIJ_RBORDER,
 408                           MPI_COMM_WORLD,
 409                           &send_requests[0]);
 410
 411         MPI_Isend(Vij[right_start_col+1],
 412                           height,
 413                           MPI_DOUBLE,
 414                           num_pes-1,
 415                           MSG_TAG_VIJ_LBORDER,
 416                           MPI_COMM_WORLD,
 417                           &send_requests[1]);
 418
 419         /* receive the columns that we need
 420            to update the columns adjacent
 421            to our strip edges from PE 1 and
 422            PE N-1 */
 423         MPI_Irecv(Vij[right_start_col-2],
 424                           height,
 425                           MPI_DOUBLE,
 426                           num_pes-1,
 427                           MSG_TAG_VIJ_RBORDER,
 428                           MPI_COMM_WORLD,
 429                           &rcv_requests[0]);
 430
 431         MPI_Irecv(Vij[left_num_cols+1],
 432                           height,
 433                           MPI_DOUBLE,
 434                           1,
 435                           MSG_TAG_VIJ_LBORDER,
 436                           MPI_COMM_WORLD,
 437                           &rcv_requests[1]);
 438
 439         MPI_Waitall(2, rcv_requests, rcv_statuses);
 440
 441         /* update the columns adjacent to our strip
 442            edges */
 443         update_voltage_array(left_num_cols, 1,Vij_TO_Vij);
 444         update_voltage_array(right_start_col - 1, 1,Vij_TO_Vij);
 445
 446         MPI_Waitall(2, send_requests, send_statuses);
 447
 448   }  /* end of accuracy loop */
 449
 450   /* set up receives for the energies from the worker PEs */
 451   for(i=1; i<num_pes; i++) {
 452         MPI_Irecv(&energies[i-1],
 453                           1,
 454                           MPI_DOUBLE,
 455                           i,
 456                           MSG_TAG_ENERGY,
 457                           MPI_COMM_WORLD,
 458                           &energy_rcv_requests[i-1]);
 459   }
 460
 461
 462   /* The energy in the matrix has now been minimised a number
 463          (accuracy) times, so we now calcuate the capacitance to see if it is
 464          converging */
 465
 466   energy_per_metre=0.0;
 467
 468   /* sum up all of the PE energies.
 469      they can arrive in any order. */
 470   for(i=1; i<num_pes; i++) {
 471         MPI_Waitany(num_pes-1,
 472                                 energy_rcv_requests,
 473                                 &index,
 474                                 &status);
 475         energy_per_metre += energies[index];
 476   }
 477
 478   /* add in the left half-strip energy */
 479   for(i=left_start_col;i<left_start_col+left_num_cols;++i)
 480         for(j=0;j<height;++j)
 481           {
 482                 energy_per_metre+=find_energy_per_metre(i,j);
 483           }
 484
 485   /* add in the right half-strip energy */
 486   for(i=right_start_col;i<right_start_col+right_num_cols;++i)
 487         for(j=0;j<height;++j)
 488           {
 489                 energy_per_metre+=find_energy_per_metre(i,j);
 490           }
 491
 492   if(coupler==FALSE)
 493         capacitance_per_metre=2*energy_per_metre;
 494   else if (coupler==TRUE)
 495         capacitance_per_metre=energy_per_metre;
 496
 497   return(capacitance_per_metre);
 498 }
 499 #endif