/* heat_equation_2d.c
 * PARALLEL [MPI] C PROGRAM TO SOLVE THE 2 DIMENSIONAL HEAT EQUATION.
 * THE INITIAL TEMPERATURE IS COMPUTED TO BE HIGH IN THE MIDDLE OF THE
 * DOMAIN UNDER CONSIDERATION AND ZERO AT THE BOUNDARIES. THE BOUNDARIES
 * ARE HELD AT ZERO TEMPERATURE THROUGH OUT THE SIMULATION.
 *
 * THE GIVEN GRID IS DECOMPOSED BY THE MASTER AND THEN DISTRIBUTED BY
 * COLUMNS TO WORKERS. A GRID POINT'S CURRENT TEMPERATURE DEPENDS UPON 
 * ITS PREVIOUS TIME STEP VALUE AS WELL AS THE VALUES OF THE NEIGHBORING
 * GRID POINTS. AS SUCH, AT EACH TIME STEP, WORKERS EXCHANGE BORDER/BOUNDARY
 * DATA WITH NEIGHBORS. ONCE ALL THE TIME STEPS ARE COMPLETED, WORKERS
 * SEND THEIR RESULTS TO MASTER.
 *
 * THE PROGRAM GENERATES A DATA FILE THAT CAN BE VISUALIZED WITH THIRD 
 * PARTY PROGRAMS [GNUPLOT, MATLAB, ETC.]
 *
 * TESTED SUCCESSFULLY WITH MPICH2 (1.3.1) COMPILED AGAINST GCC (4.1.2) 
 * IN A LINUX BOX WITH QUAD CORE INTEL XEON PROCESSOR (1.86 GHz) & 4GB OF RAM.
 *
 * ORIGINAL AUTHOR : BLAISE BARNEY; 2009/12/11
 *                   LAWRENCE LIVERMORE NATIONAL LABORATORY
 * FIRST COPIED    : GOWTHAM; Fri, 26 Nov 2010 21:30:30 -0500
 * LAST MODIFIED   : GOWTHAM; Fri, 10 Dec 2010 01:34:52 -0500
 *
 * URL:
 * http://sgowtham.net/blog/2010/12/11/mpi-c-2-dimensional-heat-equation/
 *
 * COMPILATION:
 * mpicc -g -Wall -lm heat_equation_2d.c -o heat_equation_2d.x 
 *
 * EXECUTION:
 * mpirun -np NPROC -machinefile MACHINEFILE ./heat_equation_2d.x
 *
 * NPROC       : NUMBER OF PROCESSORS ALLOCATED TO RUNNING THIS PROGRAM;
 *               NXY [GRID SIZE] MUST BE DIVISIBLE BY [NPROC - 1]
 * MACHINEFILE : FILE LISTING THE HOSTNAMES OF PROCESSORS ALLOCATED TO
 *               RUNNING THIS PROGRAM
 *
 * OUTPUT      : DATA FILE 'heat_equation_2d.dat' WITH THE FOLLOWING
 *               FORMAT
 *
 *               TSTEPS=#TSTEPS
 *               U(x, y) AT t = 0         [NXY x NXY]
 *               U(x, y) AT t = 1         [NXY x NXY]
 *               ...
 *               U(x, y) AT t = TSTEP     [NXY x NXY]
 *
 *               TOTAL NUMBER OF LINES IN 'heat_equation_2d.dat'
 *               (NXY * (TSTEPS + 1)) + 1
 *
*/
 
/* STANDARD HEADERS AND DEFINITIONS 
 * REFERENCE: http://en.wikipedia.org/wiki/C_standard_library
*/
#include <stdio.h>  /* Core input/output operations                         */
#include <stdlib.h> /* Conversions, random numbers, memory allocation, etc. */
#include <math.h>   /* Common mathematical functions                        */
#include <time.h>   /* Converting between various date/time formats         */
#include <mpi.h>    /* MPI functionality                                    */
 
#define NXY          20          /* Number of (x or y) grids in U(x,y) */
                                 /* NXY must be divisible by NPROC - 1 */
#define TSTEPS      500          /* Number of time steps               */
#define MASTER        0          /* Process ID for MASTER              */
#define BEGIN         1          /* Message tag                        */
#define DONE          2          /* Message tag                        */
#define LNEIGHBOR     3          /* Left neighbor tag                  */
#define RNEIGHBOR     4          /* Right neigbor tag                  */
#define NONE          5          /* Indicates no neighbor              */
 
struct ThermalDiffusivity {      /* Parameters used in difference      */
  double cx;                     /* equations - refer to class         */
  double cy;                     /* notes. Physically, these refer to  */
} parameters = {0.25, 0.25};     /* thermal diffusivity                */
 
 
/* MAIN PROGRAM BEGINS */
int main (int argc, char *argv[]) {
 
  /* VARIABLE DECLARATION & INITIALIZATION */
  double U[NXY][NXY],      /* 2D array to hold current values        */
         Unew[NXY][NXY],   /* 2D array to hold new values            */
         start_time,       /* Wall clock - start time                */
         end_time;         /* Wall clock - end time                  */
 
  int    proc_id,          /* Process ID                             */
         n_procs,          /* Number of processors                   */
         n_workers,        /* Number of WORKERS [nprocs - 1]         */
         cols_per_worker,  /* Number of columns of U(x,y) per WORKER */ 
         offset,           /* Offset                                 */
         source,           /* Source for a message                   */
         destination,      /* Destination for a message              */
         lneighbor,        /* Process ID of left neighbor            */
         rneighbor,        /* Process ID of right neighbor           */
         start,            /* Starting # for columns in WORKERS      */
         end,              /* Ending # for columns in WORKERS        */
         i,                /* Dummy/Running index                    */
         j,                /* Dummy/Running index                    */
         h;                /* Dummy/Running index for TSTEPS         */
 
  FILE   *output;          /* File pointer                           */
 
  MPI_Status status;       /* MPI structure containing return codes
                              for message passing operations         */
 
  MPI_Datatype cvector;    /* A custom MPI data type defined to 
                              contain column vectors                 */
 
 
  /* INITIALIZE MPI */
  MPI_Init(&argc, &argv);
 
  /* GET THE PROCESSOR ID AND NUMBER OF PROCESSORS */
  MPI_Comm_rank(MPI_COMM_WORLD, &proc_id);
  MPI_Comm_size(MPI_COMM_WORLD, &n_procs);
 
  /* NUMBER OF WORKERS [NPROC - 1] */
  n_workers = n_procs - 1;
 
  /* MPI DATA TYPE [CUSTOM] - TO SEND BLOCKS OF U(x,y) TO WORKERS AS COLUMNS
   * THANKS TO DR. STEVEN SEIDEL [http://www.cs.mtu.edu/~steve/] 
   * FOR INFORMING/TEACHING ME THE POSSIBILITY OF DEFINING CUSTOM DATA TYPES 
   * IN MPI, DURING A PREVIOUS DISCUSSION ABOUT MATRIX MULTIPILICATION
  */
  MPI_Type_vector(NXY, 1, NXY, MPI_DOUBLE, &cvector);
  MPI_Type_commit (&cvector);
 
  /* IF MASTER, THEN DO THE FOLLOWING:
   * PRINT GENERAL INFORMATION ABOUT THE PROGRAM
   * CHECK IF NPROCS MEETS THE PROGRAM REQUIREMENTS
   * POPULATE THE U(x,y) TO SIMULATE INITIAL TEMPERATURE DISTRIBUTION
   * SEND RELEVANT INFORMATION ABOUT NEIGHBORS, OF U(x,y) TO WORKERS
   * RECEIVE RELAVANT INFORMATION FROM WORKERS, AT EVERY TIME STEP
   * DISPLAY U(x,y) AT INITIAL AND FINAL STEPS
   * PRINT U(x,y) AT EVERY TIME STEP SO THAT IT CAN BE PLOTTED/ANIMATED
  */
  if (proc_id == MASTER) {
 
    printf("\n  2D Heat Equation [MPI/C]\n\n");
    printf("    Number of x-grids        : %d\n", NXY);
    printf("    Number of y-grids        : %d\n", NXY);
    printf("    Number of processors     : %d\n", n_procs);
    printf("    Number of workers        : %d\n", n_workers);
    printf("    Number of time steps     : %d\n", TSTEPS);
    printf("    Thermal diffusivity (Cx) : %4.3lf\n", parameters.cx);
    printf("    Thermal diffusivity (Cy) : %4.3lf\n\n", parameters.cy);
 
    /* CHECK IF NXY IS DIVISIBLE BY [NPROC + 1]
     * IF NOT, PRINT AN ERROR MESSAGE AND QUIT THE PROGRAM 
    */
    if (NXY%(n_procs - 1) != 0) {
      printf("    ERROR:\n");
      printf("    NXY not divisible by (NPROC + 1)\n\n");
 
      exit(1);
    }
 
    /* OPEN THE FILE FOR WRITING [OVER-WRITE, IF EXISTS]
     * PRINT ERROR MESSAGE, IF ANY & TERMINATE THE PROGRAM
    */
    if(!(output = fopen("heat_equation_2d.dat", "w"))) {
      fprintf(stderr, "\n\n    ERROR:\n");
      fprintf(stderr, "    Unable to open the output file\n\n");
 
      exit(1);
    }
 
    /* START THE CLOCK */
    start_time = MPI_Wtime();
 
    /* INITIALIZE THE GRID */
    printf("    Temperature distribution (at time t=0)\n\n");
 
    /* HEADER LINE INDICATING TOTAL NUMBER OF TIME STEPS
     * NEEDED BY MATLAB (&/OR OTHER) FOR PLOTTING
    */ 
    fprintf(output, "TSTEPS=%d\n", TSTEPS + 1);
 
    /* POPULATE THE ARRAY, U(x,y) 
     * BORDER ELEMENTS ARE ALL ZERO & HIGH VALUE IN THE MIDDLE
     * PRINT THE ARRAY - TO THE SCREEN AND TO THE FILE
    */
    for (i = 0; i < NXY; i++) {
      for (j = 0; j < NXY; j++) {
        U[i][j] = 100 * i * (NXY - i - 1) * j * (NXY - j - 1);
        printf("    %08.0lf", U[i][j]);
        fprintf(output, "%012.4lf  ", U[i][j]);
      }
      printf("\n");
      fprintf(output, "\n");
    }
 
    /* NUMBER OF COLUMNS PER WORKER 
     * IT'S IMPORTANT TO ALLOCATE EQUAL NUMBER OF COLUMNS TO
     * EACH WORKER - ELSE, SOME MIGHT FINISH FASTER WITHOUT
     * ADDITIONAL WORK
     */
    cols_per_worker = NXY/n_workers;
 
    /* INITIALIZE offset */
    offset = 0;
 
    for (i = 1; i <= n_workers; i++) {
 
      /* INFORM EACH WORKER OF ITS NEIGHBORS */
      if (i == 1) {
        lneighbor = NONE;
      } else {
        lneighbor = i - 1;
      }
 
      if (i == n_workers) {
        rneighbor = NONE;
      } else {
        rneighbor = i + 1;
      }
 
      /* DESTINATION PROCESSOR */
      destination = i;
 
      /* SEND INFORMATION TO EACH WORKER
       * OFFSET
       * NUMBER OF COLUMNS EACH WORKER HAS TO WORK WITH 
       * LEFT & RIGHT NEIGHBORS FOR EACH WORKER
       * RELEVANT PORTION OF U(x,y) TO EACH WORKER
      */
      MPI_Send(&offset,    1, MPI_INT, destination, BEGIN, MPI_COMM_WORLD);
      MPI_Send(&cols_per_worker, 1, MPI_INT, destination, BEGIN, MPI_COMM_WORLD);
      MPI_Send(&lneighbor, 1, MPI_INT, destination, BEGIN, MPI_COMM_WORLD);
      MPI_Send(&rneighbor, 1, MPI_INT, destination, BEGIN, MPI_COMM_WORLD);
 
      for(j=0; j < cols_per_worker; j++) {
        MPI_Send(&U[0][offset+j], 1, cvector, destination, BEGIN, MPI_COMM_WORLD);
      }
 
      /* INCREMENT THE offset */
      offset = offset + cols_per_worker;
    }
 
    /* RECEIVE RELEVANT INFORMATION FROM WORKERS  FOR EVERY TIME STEP */
    for(h = 1; h <= TSTEPS; h++) {
 
      for(i = 1; i <= n_workers; i++) {
 
        /* SOURCE PROCESSOR */
        source = i;
 
        MPI_Recv(&offset, 1, MPI_INT, source, DONE, MPI_COMM_WORLD, &status); 
        MPI_Recv(&cols_per_worker, 1, MPI_INT, source, DONE, MPI_COMM_WORLD, &status); 
 
        for(j=0; j < cols_per_worker; j++) {
          MPI_Recv(&U[0][offset+j], 1, cvector, source, h*DONE, MPI_COMM_WORLD, &status);
        }
 
      }
 
      /* PRINT U(x,y), THE TEMPERATURE DISTRIBUTION AFTER EVERY TIME STEP TO THE FILE */
      for (i = 0; i < NXY; i++) {
        for (j = 0; j < NXY; j++) {
          fprintf(output, "%012.4lf  ", U[i][j]);
        }
        fprintf(output, "\n");
      }
 
    }
 
    /* CLOSE THE FILE */
    fclose(output);
 
    /* PRINT U(x,y), THE TEMPERATURE DISTRIBUTION AFTER 'TSTEPS' STEPS TO THE SCREEN */
    printf("\n\n    Temperature distribution (after %d time steps)\n\n", h-1);
    for (i = 0; i < NXY; i++) {
      for (j = 0; j < NXY; j++) {
        printf("    %08.0lf", U[i][j]);
      }
      printf("\n");
    }
 
    printf("\n\n");
 
    /* STOP THE CLOCK */
    end_time = MPI_Wtime();
 
    /* PRINT TIMING INFORMATION */
    printf("    Time elapsed (sec)   : %6.4f\n\n", end_time - start_time);
 
  } /* MASTER LOOP ENDS */
 
 
  /* IF WORKER, THEN DO THE FOLLOWING:
   * INITIALIZE U(x,y) & Unew(x,y) TO ZERO
   * RECEIVE RELEVANT INFORMATION FROM MASTER
   * CALCULATE U(x, y, t+1) BY COMMUNICATING APPROPRIATELY WITH NEIGHBORING WORKERS
   * SEND RELEVANT INFORMATION TO MASTER 
  */
  if (proc_id != MASTER) {
 
    /* INITIALIZE ALL ELEMENTS OF U(x,y) & Unew(x,y) TO ZERO */
    for(i = 0; i < NXY; i++) {
      for(j = 0; j < NXY; j++) {
        U[i][j] = 0.0;
        Unew[i][j] = 0.0;
      }
    }
 
    /* RECEIVE THE COLUMNS, NEIGHBORS AND GRID INFORMATION FROM MASTER */
    MPI_Recv(&offset,    1, MPI_INT, MASTER, BEGIN, MPI_COMM_WORLD, &status);
    MPI_Recv(&cols_per_worker, 1, MPI_INT, MASTER, BEGIN, MPI_COMM_WORLD, &status);
    MPI_Recv(&lneighbor, 1, MPI_INT, MASTER, BEGIN, MPI_COMM_WORLD, &status);
    MPI_Recv(&rneighbor, 1, MPI_INT, MASTER, BEGIN, MPI_COMM_WORLD, &status);
 
    for(j=0; j < cols_per_worker; j++) {
      MPI_Recv(&U[0][offset+j], 1, cvector, MASTER, BEGIN, MPI_COMM_WORLD, &status);
    }
 
    /* FOR DEBUGGING PURPOSES ONLY 
     * RUN WITH NPROC = 3 WITH NXY = 16
     * THERE SHOULD BE TWO BLOCKS OF COMPLEMENTARY MATRICES
 
    if (proc_id == 1) {
    printf("\n");
      for (i = 0; i < NXY; i++) {
        for (j = 0; j < cols_per_worker; j++) {
          printf(" %08.0lf", U[0][i][j]);
        }
        printf("\n");
      }
    }
 
    if (proc_id == 2) {
    printf("\n");
      for (i = 0; i < NXY; i++) {
        for (j = cols_per_worker; j < 2*cols_per_worker; j++) {
          printf(" %08.0lf", U[0][i][j]);
        }
        printf("\n");
      }
    }
    */ 
 
 
    /* DETERMINE BORDER ELEMENTS - FIRST & LAST COLUMNS NEED EXTRA CARE:
     * ROW 0 CANNOT EXCHANGE INFORMATION WITH 0-1
     * ROW N CANNOT EXCHANGE INFORMATION WITH N+1
    */
    if (offset == 0) {
      /* DO NOT INCLUDE ZEROTH COLUMN */
      start = 1;
    } else {
      start = offset;
    }
 
    if ((offset + cols_per_worker) == NXY) {
      /* DO NOT INCLUDE THE LAST COLUMN */
      end = offset + cols_per_worker - 2;
    } else {
      end = offset + cols_per_worker - 1;
    }
 
    /* PERFORM TIME STEP [TSTEPS] ITERATIONS.
     * MUST COMMUNICATE BORDER COLUMNS WITH NEIGHBORS.
     * IF IT'S THE FIRST OR THE LAST COLUMN, THEN THE COMMUNICATION TAKES 
     * PLACE ONLY WITH ONE NEIGHBOR
    */
    for(h = 1; h <= TSTEPS; h++) {
      if (lneighbor != NONE) {
        MPI_Send(&U[0][offset], 1, cvector, lneighbor, RNEIGHBOR, MPI_COMM_WORLD);
 
        MPI_Recv(&U[0][offset-1], 1, cvector, lneighbor, LNEIGHBOR, MPI_COMM_WORLD, 
                 &status);
      }
 
      if (rneighbor != NONE) {
        MPI_Send(&U[0][offset+cols_per_worker-1], 1, cvector, rneighbor,
                 LNEIGHBOR, MPI_COMM_WORLD);
 
        MPI_Recv(&U[0][offset+cols_per_worker], 1, cvector, rneighbor,
                 RNEIGHBOR, MPI_COMM_WORLD, &status);
      }
 
      /* UPDATE THE VALUE OF GRID POINTS 
       * REFER TO THE SLIDES FROM WEEK #14 [PART #02]
      */
      for(i = 1; i <= NXY-2; i++) {
        for(j = start; j <= end; j++) {
          Unew[i][j] = U[i][j] + 
                       parameters.cx * (U[i-1][j] + U[i+1][j] - 2 * (U[i][j])) +
                       parameters.cy * (U[i][j-1] + U[i][j+1] - 2 * (U[i][j]));
        }
      }
 
      /* SET U(x,y) TO Unew(x,y) */
      for(i = 1; i <= NXY-2; i++) {
        for(j = start; j <= end; j++) {
          U[i][j] = Unew[i][j];
        }
      }
 
      /* SEND RELEVANT INFORMATION TO MASTER */
      MPI_Send(&offset, 1, MPI_INT, MASTER, DONE, MPI_COMM_WORLD);
      MPI_Send(&cols_per_worker, 1, MPI_INT, MASTER, DONE, MPI_COMM_WORLD);
 
      for(j=0; j < cols_per_worker; j++) {
        MPI_Send(&U[0][offset+j], 1, cvector, MASTER, h*DONE, MPI_COMM_WORLD);
      }
    }
 
  } /* WORKER LOOP ENDS */
 
 
  /* FINALIZE MPI */
  MPI_Finalize();
 
  /* INDICATE THE TERMINATION OF THE PROGRAM */
  return 0;
 
} /* MAIN PROGRAM ENDS */

[guest@dirac mpi_samples]$ which mpicc
alias mpicc='mpicc -g -Wall -lm'
	~/mpich2/1.3.1/gcc/4.1.2/bin/mpicc
 
[guest@dirac mpi_samples]$ which mpirun
alias mpirun='mpirun -machinefile $HOME/machinefile'
	~/mpich2/1.3.1/gcc/4.1.2/bin/mpirun
 
[guest@dirac mpi_samples]$ mpicc heat_equation_2d.c -o heat_equation_2d.x
 
[guest@dirac mpi_samples]$ mpirun -np 5 ./heat_equation_2d.x
 
  2D Heat Equation [MPI/C]
 
    Number of x-grids        : 20
    Number of y-grids        : 20
    Number of processors     : 5
    Number of workers        : 4
    Number of time steps     : 5000
    Thermal diffusivity (Cx) : 0.050
    Thermal diffusivity (Cy) : 0.050
 
    Temperature distribution (at time t=0)
 
    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000
    00000000    00032400    00061200    00086400    00108000    00126000    00140400    00151200    00158400    00162000    00162000    00158400    00151200    00140400    00126000    00108000    00086400    00061200    00032400    00000000
    00000000    00061200    00115600    00163200    00204000    00238000    00265200    00285600    00299200    00306000    00306000    00299200    00285600    00265200    00238000    00204000    00163200    00115600    00061200    00000000
    00000000    00086400    00163200    00230400    00288000    00336000    00374400    00403200    00422400    00432000    00432000    00422400    00403200    00374400    00336000    00288000    00230400    00163200    00086400    00000000
    00000000    00108000    00204000    00288000    00360000    00420000    00468000    00504000    00528000    00540000    00540000    00528000    00504000    00468000    00420000    00360000    00288000    00204000    00108000    00000000
    00000000    00126000    00238000    00336000    00420000    00490000    00546000    00588000    00616000    00630000    00630000    00616000    00588000    00546000    00490000    00420000    00336000    00238000    00126000    00000000
    00000000    00140400    00265200    00374400    00468000    00546000    00608400    00655200    00686400    00702000    00702000    00686400    00655200    00608400    00546000    00468000    00374400    00265200    00140400    00000000
    00000000    00151200    00285600    00403200    00504000    00588000    00655200    00705600    00739200    00756000    00756000    00739200    00705600    00655200    00588000    00504000    00403200    00285600    00151200    00000000
    00000000    00158400    00299200    00422400    00528000    00616000    00686400    00739200    00774400    00792000    00792000    00774400    00739200    00686400    00616000    00528000    00422400    00299200    00158400    00000000
    00000000    00162000    00306000    00432000    00540000    00630000    00702000    00756000    00792000    00810000    00810000    00792000    00756000    00702000    00630000    00540000    00432000    00306000    00162000    00000000
    00000000    00162000    00306000    00432000    00540000    00630000    00702000    00756000    00792000    00810000    00810000    00792000    00756000    00702000    00630000    00540000    00432000    00306000    00162000    00000000
    00000000    00158400    00299200    00422400    00528000    00616000    00686400    00739200    00774400    00792000    00792000    00774400    00739200    00686400    00616000    00528000    00422400    00299200    00158400    00000000
    00000000    00151200    00285600    00403200    00504000    00588000    00655200    00705600    00739200    00756000    00756000    00739200    00705600    00655200    00588000    00504000    00403200    00285600    00151200    00000000
    00000000    00140400    00265200    00374400    00468000    00546000    00608400    00655200    00686400    00702000    00702000    00686400    00655200    00608400    00546000    00468000    00374400    00265200    00140400    00000000
    00000000    00126000    00238000    00336000    00420000    00490000    00546000    00588000    00616000    00630000    00630000    00616000    00588000    00546000    00490000    00420000    00336000    00238000    00126000    00000000
    00000000    00108000    00204000    00288000    00360000    00420000    00468000    00504000    00528000    00540000    00540000    00528000    00504000    00468000    00420000    00360000    00288000    00204000    00108000    00000000
    00000000    00086400    00163200    00230400    00288000    00336000    00374400    00403200    00422400    00432000    00432000    00422400    00403200    00374400    00336000    00288000    00230400    00163200    00086400    00000000
    00000000    00061200    00115600    00163200    00204000    00238000    00265200    00285600    00299200    00306000    00306000    00299200    00285600    00265200    00238000    00204000    00163200    00115600    00061200    00000000
    00000000    00032400    00061200    00086400    00108000    00126000    00140400    00151200    00158400    00162000    00162000    00158400    00151200    00140400    00126000    00108000    00086400    00061200    00032400    00000000
    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000
 
 
    Temperature distribution (after 5000 time steps)
 
    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000000    00000000    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000000    00000000    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000000    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000000    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000000    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000000    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000000    00000000    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000000    00000000    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000
 
 
    Time elapsed (sec)   : 2.5108
 
[guest@dirac mpi_samples]$

% heat_equation_2d.m
% MATLAB 'SCRIPT m' FILE TO READ THE FILE 'heat_equation_2d.dat', 
% ONE DATA SET AT A TIME AND PLOT IT REPEATEDLY TO SIMULATE AN ANIMATION.
%
% TESTED SUCCESSFULLY WITH MATLAB (v7.8.0.347 R2009a; STUDENT VERSION) IN A MACBOOK PRO 
% (OS X 10.6.5) WITH INTEL CORE 2 DUO PROCESSOR (2.4 GHz) & 4GB OF RAM.
%
% FIRST WRITTEN: GOWTHAM; Fri, 10 Dec 2010 23:12:18 -0500
% LAST MODIFIED: GOWTHAM; Sun, 12 Dec 2010 09:34:23 -0500
%
% INPUT:
% 'heat_equation_2d.dat' PRODUCED BY RUNNING 'heat_equation_2d.c'
% AND IT HAS THE FOLLOWING FORMAT:
% 
% TSTEPS=# 
% U(x,y) AT t = 0         [NXY x NXY]
% U(x,y) AT t = 1         [NXY x NXY]
% U(x,y) AT t = 2         [NXY x NXY]
% ...
% U(x,y) AT t = TSTEPS    [NXY x NXY]
%
% FOR AN EXPLANATION OF 'NXY', 'TSTEPS', ETC., PLEASE REFER TO 'heat_equation_2d.c'.
% 
% OUTPUT:
% IF 'mencoder' [AND OTHER REQUIRED UTILITIES] IS AVAILABLE, THIS SCRIP WILL PRODUCE 
% 'heat_equation_2d.avi'. IF NOT, 'New Screen Recording' OPTION MAY BE USED IN
% 'QuickTime Player' TO CREATE THE ANIMATION.
%
 
% CLEAR THE PREVIOUSLY SET VARIABLES
clear
 
% CLEAR THE COMMAND WINDOW
clc
 
% PAUSE FOR 15 SECONDS
% USEFUL WHEN 'mencoder' IS NOT AVAILABLE AS THIS WILL GIVE JUST ENOUGH TIME TO HIDE 
% UNWANTED SCREENS / WINDOWS AND START RECORDING VIA QuickTime Player.
% COMMENT IT OUT IF NOT NEEDED.
pause(15);
 
% OPEN THE DATA FILE FOR READING. 
% SAME BASENAME WILL BE USED FOR THE ANIMATION FILE
filename   = 'heat_equation_2d';
filehandle = fopen(filename '.dat', 'r');
 
% COUNT THE NUMBER OF TIME STEPS - GIVEN IN THE FIRST LINE
TSTEPS = fscanf(filehandle, 'TSTEPS=%d\n', 1);
 
% VARIABLE DECLARATION AND INITIALIZATION
% THIS IS THE NUMBER OF ROWS & COLUMNS IN EACH BLOCK OF DATA
NXY = 20;
 
% x & y GRIDS FOR THE MESH [3D] PLOT
x = 1:1:NXY;
y = 1:1:NXY;
 
% GENERATE X & Y ARRAYS FOR 3D PLOT
% TRANSFORMS THE DOMAIN SPECIFIED BY VECTORS x & y INTO ARRAYS X and Y, 
% WHICH CAN BE USED TO EVALUATE FUNCTIONS OF TWO VARIABLES AND 3D
% MESH/SURFACE PLOTS. THE ROWS OF THE OUTPUT ARRAY X ARE COPIES OF THE
% VECTOR x; COLUMNS OF THE OUTPUT ARRAY Y ARE COPIES OF THE VECTOR y.
% http://www.mathworks.com/help/techdoc/ref/meshgrid.html
[X, Y] = meshgrid(x,y);
 
% LOOP THROUGH THE DATA FILE TO READ 'NXY x NXY' ELEMENTS AT ONCE
% AND IMPORT THEM INTO 'TSTEPS' NUMBER OF MATRICES OF ORDER NXY.
for i = 1:TSTEPS
 
  % INITIALIZE ALL ELEMENTS OF MATRIX [OF ORDER NXY] TO ZERO
  Temperature(i).Z = zeros(NXY, NXY);
 
  % SCAN THE DATA FILE, READ 'NXY x NXY' ELEMENTS, TREAT THEM
  % AS DOUBLE PRECISION AND STORE THEM IN THE MATRIX.
  % SYNTAX:
  % Z = fscanf(filehandle, 'FORMAT', [COLS, ROWS]);
  Temperature(i).Z = fscanf(filehandle, '%f', [NXY, NXY]);
 
  % PLOT THE DATA IN THE ABOVE MATRIX AS A FUNCTION OF X & Y
  figure1 = figure(1);
  mesh(X, Y, Temperature(i).Z)
 
  % HANDLE FOR AXIS PROPERTY SO THAT IT CAN BE MODIFIED TO
  % FIT OUR REQUIREMENS
  ap = gca;
 
  % MODIFY CERTAIN AXIS PROPERTIES APPROPRIATELY
  set(ap, 'XLimMode', 'manual', 'YLimMode', 'manual', ...
          'ZLimMode', 'manual', ...
          'Xlim', [0 20], 'Ylim', [0 20], 'Zlim', [0 1000000], ...
          'FontSize', 14)
 
  % CHANGE THE COLORMAP FROM 'default' TO 'hsv'
  % THE hsv COLORMAP IS CREATED WITH 65536 COLORS
  % http://www.mathworks.com/help/techdoc/ref/colormap.html
  colormap(hsv(65536))
 
  % DISPLAY A COLOR BAR THAT WILL SERVE AS A LEGEND
  colorbar
 
  % AXES LABELS FOR THE PLOT
  xlabel('Grid points along x-axis');
  ylabel('Grid points along y-axis');
  zlabel('U(x, y), Temperature distribution');
 
  % TITLE FOR THE PLOT
  title(['TSTEP= ' num2str(i)]);
 
  % CAPTURE MOVIE FRAME AS AN ARRAY FOR LATER USE.
  % COMMENT THIS LINE IF 'mencoder' IS NOT AVAILABLE.
  % http://www.mathworks.com/help/techdoc/ref/getframe.html
  M(i) = getframe;
 
  % TIME LAG BETWEEN SUCCESSIVE PLOTS
  pause(0.00001);
 
end 
% LOOP THROUGH DATA FILE ENDS
 
% CLOSE THE DATA FILE
fclose(datafilename);
 
 
% CREATE AN ANIMATION USING THE MOVIE FRAMES.
% COMMENT THIS SECTION IF 'mencoder' IS NOT AVAILABLE
 
% CREATE AN AVI [AUDIO/VIDEO INTERLEAVED] FILE FROM THE MOVIE FRAMES ARRAY
% 'None' IS THE ONLY COMPRESSION AVAILABLE OPTION ON UNIX (LIKE) OS.
% http://www.mathworks.com/help/techdoc/ref/movie2avi.html
movie2avi(M, [filename '.avi'], 'compression', 'None', 'fps', 25);
 
% USE 'mencoder' TO COMPRESS THE MOVIE SO THAT IT CAN BE PLAYED ON 
% QuickTime Player OR SOME OTHER SUCH APPLICATION.
% USE eval TO RUN THE MAC/LINUX COMMAND
% http://www.mplayerhq.hu/DOCS/HTML/en/menc-feat-selecting-codec.html
avicompression = ['!mencoder -oac copy -ovc xvid -xvidencopts bitrate=1000 ' filename '.avi -o ' filename '_compressed.avi'];
eval(avicompression);
 
% RENAME THE COMPRESSED AVI
% USE eval TO RUN THE MAC/LINUX COMMAND
avirename = ['!mv ' filename '_compressed.avi ' filename '.avi'];
eval(avirename);