added MPI examples

examples/05_mpi/CMakeLists.txt

+cmake_minimum_required(VERSION 2.8)
+find_package(MPI)
+find_package(BLAS)
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11")
+include_directories(${MPI_CXX_INCLUDE_PATH})
+set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${MPI_CXX_LINK_FLAGS}")
+set (FILES deadlock1 deadlock2 deadlock3 deadlock4 deadlock5 diffusion1d diffusion1d_nonblocking errors example1 example2 series1 series2 series3 simpson1 simpson2 simpson3 simpson4 pdaxpy pddot pdscal)
+foreach(f ${FILES})
+  add_executable(${f} ${f}.cpp)
+  target_link_libraries(${f} ${MPI_CXX_LIBRARIES} ${BLAS_LIBRARIES})
+endforeach(f ${FILES})  

examples/05_mpi/README.txt

+To build the examples using cmake, you can type at the command prompt: mkdir build; cd build; cmake -DCMAKE_BUILD_TYPE=Release ..
+Required packages: cmake, MPI, BLAS 
+For Euler: module load gcc; module load mvapich2 (or open_mpi); module load mkl
+
+
+To build the examples manually, following this example: mpic++ -o example1 example1.cpp 

examples/05_mpi/aligned_allocator.hpp

+// Example codes for HPC course
+// (c) 2012 Andreas Hehn, ETH Zurich
+
+#ifndef HPC12_ALIGNED_ALLOCATOR_HPP
+#define HPC12_ALIGNED_ALLOCATOR_HPP
+
+#ifdef _WIN32
+#include <malloc.h>
+#else
+#include <cstdlib>
+#endif
+
+#if __cplusplus >= 201103L 
+#define NOEXCEPT_SPEC noexcept
+#else
+#define NOEXCEPT_SPEC
+#endif
+
+namespace hpc12 {
+
+// Alignment must be a power of 2 !
+template <typename T, unsigned int Alignment>
+class aligned_allocator {
+  public:
+    typedef T*              pointer;
+    typedef T const*        const_pointer;
+    typedef T&              reference;
+    typedef T const&        const_reference;
+    typedef T               value_type;
+    typedef std::size_t     size_type;
+    typedef std::ptrdiff_t  difference_type;
+
+    template <typename U>
+    struct rebind {
+        typedef aligned_allocator<U,Alignment> other;
+    };
+
+    aligned_allocator() NOEXCEPT_SPEC {
+    }
+
+    aligned_allocator(aligned_allocator const& a) NOEXCEPT_SPEC {
+    }
+
+    template <typename U>
+    aligned_allocator(aligned_allocator<U,Alignment> const& b) NOEXCEPT_SPEC {
+    }
+
+    pointer allocate(size_type n) {
+        pointer p;
+#ifdef _WIN32
+        p = _aligned_malloc(n*sizeof(T), Alignment);
+        if(p == 0)
+            throw std::bad_alloc();
+#else
+        if(posix_memalign(reinterpret_cast<void**>(&p), Alignment, n * sizeof(T) ))
+            throw std::bad_alloc();
+#endif
+        return p;
+    }
+
+    void deallocate(pointer p, size_type n) NOEXCEPT_SPEC {
+        std::free(p);
+    }
+
+    size_type max_size() const NOEXCEPT_SPEC {
+        std::allocator<T> a;
+        return a.max_size();
+    }
+
+#if __cplusplus >= 201103L 
+    template <typename C, class... Args>
+    void construct(C* c, Args&&... args) {
+        new ((void*)c) C(std::forward<Args>(args)...);
+    }
+#else
+    void construct(pointer p, const_reference t) {
+        new((void *)p) T(t);
+    }
+#endif
+
+    template <typename C>
+    void destroy(C* c) {
+        c->~C();
+    }
+
+    bool operator == (aligned_allocator const& a2) const NOEXCEPT_SPEC {
+        return true;
+    }
+
+    bool operator != (aligned_allocator const& a2) const NOEXCEPT_SPEC {
+        return false;
+    }
+
+    template <typename U, unsigned int UAlignment>
+    bool operator == (aligned_allocator<U,UAlignment> const& b) const NOEXCEPT_SPEC {
+        return false;
+    }
+
+    template <typename U, unsigned int UAlignment>
+    bool operator != (aligned_allocator<U,UAlignment> const& b) const NOEXCEPT_SPEC {
+        return true;
+    }
+};
+
+}
+
+#undef NOEXPECT_SPEC
+
+#endif //HPC12_ALIGNED_ALLOCATOR_HPP

examples/05_mpi/deadlock1.cpp

+// Example codes for HPC course
+// (c) 2012 Matthias Troyer, ETH Zurich
+
+#include <iostream>
+#include <string>
+#include <mpi.h>
+
+int main(int argc, char** argv) {
+  MPI_Status status;
+  int num;
+  
+  MPI_Init(&argc, &argv);
+  MPI_Comm_rank(MPI_COMM_WORLD,&num);
+  
+  double ds=3.1415927; // to send
+  double dr;           // to receive
+  int tag=99;
+  
+  if(num==0) {
+    MPI_Ssend(&ds,1,MPI_DOUBLE,1,tag,MPI_COMM_WORLD);
+    MPI_Recv (&dr,1,MPI_DOUBLE,1,tag,MPI_COMM_WORLD,&status);
+  }
+  else {
+    MPI_Ssend(&ds,1,MPI_DOUBLE,0,tag,MPI_COMM_WORLD);
+    MPI_Recv (&dr,1,MPI_DOUBLE,0,tag,MPI_COMM_WORLD,&status);
+  }
+  
+  MPI_Finalize();
+  return 0;
+}

examples/05_mpi/deadlock2.cpp

+// Example codes for HPC course
+// (c) 2012 Matthias Troyer, ETH Zurich
+
+#include <iostream>
+#include <string>
+#include <mpi.h>
+
+int main(int argc, char** argv) {
+  MPI_Status status;
+  int num;
+  
+  MPI_Init(&argc, &argv);
+  MPI_Comm_rank(MPI_COMM_WORLD,&num);
+  
+  double ds=3.1415927; // to send
+  double dr;           // to receive
+  int tag=99;
+  
+  if(num==0) {
+    MPI_Ssend(&ds,1,MPI_DOUBLE,1,tag,MPI_COMM_WORLD);
+    MPI_Recv (&dr,1,MPI_DOUBLE,1,tag,MPI_COMM_WORLD,&status);
+  }
+  else {
+    MPI_Recv (&dr,1,MPI_DOUBLE,0,tag,MPI_COMM_WORLD,&status);
+    MPI_Ssend(&ds,1,MPI_DOUBLE,0,tag,MPI_COMM_WORLD);
+  }
+  
+  MPI_Finalize();
+  return 0;
+}

examples/05_mpi/deadlock3.cpp

+// Example codes for HPC course
+// (c) 2012 Matthias Troyer, ETH Zurich
+
+#include <iostream>
+#include <string>
+#include <mpi.h>
+
+int main(int argc, char** argv) {
+  MPI_Status status;
+  int num;
+  
+  MPI_Init(&argc, &argv);
+  MPI_Comm_rank(MPI_COMM_WORLD,&num);
+  
+  double ds=3.1415927; // to send
+  double dr;           // to receive
+  int tag=99;
+  
+  if(num==0) {
+    MPI_Sendrecv(&ds,1,MPI_DOUBLE,1,tag,
+                 &dr,1,MPI_DOUBLE,1,tag,
+                 MPI_COMM_WORLD,&status);
+  }
+  else {
+    MPI_Sendrecv(&ds,1,MPI_DOUBLE,0,tag,
+                 &dr,1,MPI_DOUBLE,0,tag,
+                 MPI_COMM_WORLD,&status);
+  }
+  
+  MPI_Finalize();
+  return 0;
+}

examples/05_mpi/deadlock4.cpp

+// Example codes for HPC course
+// (c) 2012 Matthias Troyer, ETH Zurich
+
+#include <iostream>
+#include <string>
+#include <mpi.h>
+
+int main(int argc, char** argv) {
+  MPI_Status status;
+  int num;
+  
+  MPI_Init(&argc, &argv);
+  MPI_Comm_rank(MPI_COMM_WORLD,&num);
+  
+  double ds=3.1415927; // to send
+  double dr;           // to receive
+  int tag=99;
+  
+  // allocate a buffer and attach it to MPI
+  int buffer_size = sizeof(double) + MPI_BSEND_OVERHEAD;
+  char* buffer = new char[buffer_size];
+  MPI_Buffer_attach(buffer,buffer_size);
+  
+  
+  if(num==0) {
+    MPI_Bsend(&ds,1,MPI_DOUBLE,1,tag,MPI_COMM_WORLD);
+    MPI_Recv (&dr,1,MPI_DOUBLE,1,tag,MPI_COMM_WORLD,&status);
+  }
+  else {
+    MPI_Bsend(&ds,1,MPI_DOUBLE,0,tag,MPI_COMM_WORLD);
+    MPI_Recv (&dr,1,MPI_DOUBLE,0,tag,MPI_COMM_WORLD,&status);
+  }
+  
+  // detach the buffer, making sure all sends are done
+  MPI_Buffer_detach(buffer,&buffer_size);
+  delete[] buffer;
+  
+
+  MPI_Finalize();
+  return 0;
+}

examples/05_mpi/deadlock5.cpp

+// Example codes for HPC course
+// (c) 2012 Matthias Troyer, ETH Zurich
+
+#include <iostream>
+#include <string>
+#include <mpi.h>
+
+int main(int argc, char** argv) {
+  MPI_Status status;
+  int num;
+  
+  MPI_Init(&argc, &argv);
+  MPI_Comm_rank(MPI_COMM_WORLD,&num);
+  
+  double ds=3.1415927; // to send
+  double dr;           // to receive
+  int tag=99;
+  
+  if(num==0) {
+    MPI_Send(&ds,1,MPI_DOUBLE,1,tag,MPI_COMM_WORLD);
+    MPI_Recv (&dr,1,MPI_DOUBLE,1,tag,MPI_COMM_WORLD,&status);
+  }
+  else {
+    MPI_Send(&ds,1,MPI_DOUBLE,0,tag,MPI_COMM_WORLD);
+    MPI_Recv (&dr,1,MPI_DOUBLE,0,tag,MPI_COMM_WORLD,&status);
+  }
+  
+  MPI_Finalize();
+  return 0;
+}

examples/05_mpi/diffusion1d.cpp

+// Example codes for HPC course
+// (c) 2012 Matthias Troyer, ETH Zurich
+
+
+#include <vector>
+#include <algorithm>
+#include <iostream>
+#include <iterator>
+#include <cassert>
+#include <mpi.h>
+
+int main(int argc, char** argv)
+{
+  double L  = 10000.;  // extent of the domain
+  double dx = 0.1;  // spatial resolution
+  double dt = 0.005;  // time step
+  double c  = 1.;    // diffusion constant
+
+  double coefficient = c*dt/dx/dx;
+  
+  unsigned iterations = 10000; // number of iterations
+  unsigned N          = L/dx;  // number of mesh points
+  
+  // now initialize MPI and get information about the number of processes
+  
+  MPI_Init(&argc,&argv);
+  
+  int size;
+  int rank;
+  MPI_Status status;
+  MPI_Comm_size(MPI_COMM_WORLD,&size);
+  MPI_Comm_rank(MPI_COMM_WORLD,&rank);
+ 
+  // get neighbors and assume periodic boundary conditions;
+  int left = (rank > 0 ? rank-1 : size-1);
+  int right = (rank < size-1 ? rank+1 : 0);
+
+  // get number of points to be used on this process
+  int points_per_rank =  (N+size-1)/size;
+  int local_N = points_per_rank;
+  // take care of roundoff issues for the last segment
+  // if N is not a multiple of size
+  if (rank==size-1)
+    local_N = N-(size-1)*points_per_rank;
+  
+  // add space for the two ghost cells on the left and right end
+  local_N+=2;
+  std::vector<double> density(local_N);
+  std::vector<double> newdensity(local_N);
+  
+  std::cout << "Rank " << rank << " has " << local_N << " points\n";
+  
+  // set the density to 1 in the middle
+  // we need to shift and limit to our range
+  int mystart = rank*points_per_rank;
+  int start1 = N/4-mystart+1; // take care of ghost cell offset
+  int end1 = N/4-mystart+1;   // take care of ghost cell offset
+  
+  if (start1 < local_N && end1 >= 0)
+    std::fill(density.begin()+std::max(0,start1),
+              density.begin()+std::min(end1,local_N),1.);
+
+  for (int t=0; t<iterations; ++t) {
+    // first get the ghost cells and send our boundary values to
+    // the neighbor for their ghost cells
+    
+    // avoid deadlocks by a clear ordering who sends and receives first
+    // even ones send left and odd ones receive right
+    // even ones receive left and odd ones send right
+    // even ones send right and odd ones receive left
+    // even ones receive right and odd ones send left
+    
+    // make sure we have an even number of ranks for this to work
+    assert(size %2 == 0);
+    
+    if (rank % 2 == 0) {
+      MPI_Send(&density[1],1,MPI_DOUBLE,left,0,MPI_COMM_WORLD);
+      MPI_Recv(&density[0],1,MPI_DOUBLE,left,0,MPI_COMM_WORLD,&status);
+      MPI_Send(&density[local_N-2],1,MPI_DOUBLE,right,0,MPI_COMM_WORLD);
+      MPI_Recv(&density[local_N-1],1,MPI_DOUBLE,right,0,MPI_COMM_WORLD,&status);
+    }
+    else {
+      MPI_Recv(&density[local_N-1],1,MPI_DOUBLE,right,0,MPI_COMM_WORLD,&status);
+      MPI_Send(&density[local_N-2],1,MPI_DOUBLE,right,0,MPI_COMM_WORLD);
+      MPI_Recv(&density[0],1,MPI_DOUBLE,left,0,MPI_COMM_WORLD,&status);
+      MPI_Send(&density[1],1,MPI_DOUBLE,left,0,MPI_COMM_WORLD);
+    }
+
+    // do calculation 
+    for (int i=1; i<local_N-1;++i)
+      newdensity[i] = density[i] + coefficient * (density[i+1]+density[i-1]-2.*density[i]);
+    
+    // and swap
+    density.swap(newdensity);
+  }
+  
+  MPI_Finalize();
+ }
+

examples/05_mpi/diffusion1d_nonblocking.cpp

+// Example codes for HPC course
+// (c) 2012 Matthias Troyer, ETH Zurich
+
+
+#include <vector>
+#include <algorithm>
+#include <iostream>
+#include <iterator>
+#include <cassert>
+#include <mpi.h>
+
+int main(int argc, char** argv)
+{
+  double L  = 10000.;  // extent of the domain
+  double dx = 0.1;  // spatial resolution
+  double dt = 0.005;  // time step
+  double c  = 1.;    // diffusion constant
+
+  double coefficient = c*dt/dx/dx;
+  
+  unsigned iterations = 10000; // number of iterations
+  unsigned N          = L/dx;  // number of mesh points
+  
+  // now initialize MPI and get information about the number of processes
+  
+  MPI_Init(&argc,&argv);
+  
+  int size;
+  int rank;
+  MPI_Status status[4];
+  MPI_Request reqs[4];
+  MPI_Comm_size(MPI_COMM_WORLD,&size);
+  MPI_Comm_rank(MPI_COMM_WORLD,&rank);
+ 
+  // get neighbors and assume periodic boundary conditions;
+  int left = (rank > 0 ? rank-1 : size-1);
+  int right = (rank < size-1 ? rank+1 : 0);
+
+  // get number of points to be used on this process
+  int points_per_rank =  N/size;
+  int local_N = points_per_rank;
+  // take care of roundoff issues for the last segment
+  // if N is not a multiple of size
+  if (rank==size-1)
+    local_N = N-(size-1)*points_per_rank;
+  
+  // add space for the two ghost cells on the left and right end
+  local_N+=2;
+  std::vector<double> density(local_N);
+  std::vector<double> newdensity(local_N);
+  
+  std::cout << "Rank " << rank << " has " << local_N << " points\n";
+  
+  // set the density to 1 in the middle
+  // we need to shift and limit to our range
+  int mystart = rank*points_per_rank;
+  int start1 = N/4-mystart+1; // take care of ghost cell offset
+  int end1 = N/4-mystart+1;   // take care of ghost cell offset
+  
+  if (start1 < local_N && end1 >= 0)
+    std::fill(density.begin()+std::max(0,start1),
+              density.begin()+std::min(end1,local_N),1.);
+
+  for (int t=0; t<iterations; ++t) {
+    // first start the communication to get  the ghost cells and send
+    // our boundary values to the neighbor for their ghost cells
+    
+    // avoid deadlocks by a clear ordering who sends and receives first
+    // even ones send left and odd ones receive right
+    // even ones receive left and odd ones send right
+    // even ones send right and odd ones receive left
+    // even ones receive right and odd ones send left
+    
+    // make sure we have an even number of ranks for this to work
+    assert(size %2 == 0);
+    if (rank % 2 == 0) {
+      MPI_Isend(&density[1],1,MPI_DOUBLE,left,0,MPI_COMM_WORLD,&reqs[0]);
+      MPI_Irecv(&density[0],1,MPI_DOUBLE,left,0,MPI_COMM_WORLD,&reqs[1]);
+      MPI_Isend(&density[local_N-2],1,MPI_DOUBLE,right,0,MPI_COMM_WORLD,&reqs[2]);
+      MPI_Irecv(&density[local_N-1],1,MPI_DOUBLE,right,0,MPI_COMM_WORLD,&reqs[3]);
+    }
+    else {
+      MPI_Irecv(&density[local_N-1],1,MPI_DOUBLE,right,0,MPI_COMM_WORLD,&reqs[0]);
+      MPI_Isend(&density[local_N-2],1,MPI_DOUBLE,right,0,MPI_COMM_WORLD,&reqs[1]);
+      MPI_Irecv(&density[0],1,MPI_DOUBLE,left,0,MPI_COMM_WORLD,&reqs[2]);
+      MPI_Isend(&density[1],1,MPI_DOUBLE,left,0,MPI_COMM_WORLD,&reqs[3]);
+    }
+
+    // do calculation  of the interior
+    for (int i=2; i<local_N-2;++i)
+      newdensity[i] = density[i] + coefficient * (density[i+1]+density[i-1]-2.*density[i]);
+    
+    // wait for the ghost cells to arrive
+    MPI_Waitall(4,reqs,status);
+    
+    // do the boundaries
+    newdensity[1] = density[1] + coefficient * (density[2]+density[0]-2.*density[1]);
+    newdensity[local_N-2] = density[local_N-2] + coefficient * (
+                               density[local_N-1]+density[local_N-3]-2.*density[local_N]);
+    
+    // and swap
+    density.swap(newdensity);
+  }
+  
+  MPI_Finalize();
+ }
+

examples/05_mpi/dvector.hpp

+// Example codes for HPC course
+// (c) 2012 Andreas Hehn and Matthias Troyer, ETH Zurich
+
+// a stripped-down distributed vector class for the examples of the HPC course
+
+#ifndef HPC12_DVECTOR_HPP
+#define HPC12_DVECTOR_HPP
+
+#include <cassert>
+#include <iostream>
+#include <vector>
+#include "aligned_allocator.hpp"
+
+#include <mpi.h>
+
+
+template <typename T, typename Allocator = hpc12::aligned_allocator<T,64> >
+class dvector : public std::vector<T,Allocator>
+{
+  public:
+    typedef T value_type;
+  
+  dvector(std::size_t n, MPI_Comm c = MPI_COMM_WORLD)
+  : comm_(c)
+  , global_size_(n)
+  {
+    int s;
+    MPI_Comm_rank(comm_,&rank_);
+    MPI_Comm_size(comm_,&s);
+    
+    // calculate the local block size
+    block_size_ = (global_size_+s-1)/s;
+    if (rank_*block_size_ < global_size_);
+      this->resize(std::min(block_size_,global_size_-rank_*block_size_));
+  }
+
+  value_type const* data() const {
+      return this->empty() ? 0 : &this->front();
+  }
+
+  value_type* data() {
+      return this->empty() ? 0 : &this->front();
+  }
+
+  // a swap function should be added
+  
+  std::size_t global_size() const { return global_size_;}
+  std::size_t offset() const { return rank_ * block_size_;}
+
+  std::size_t block_size() const { return block_size_;}
+  
+  MPI_Comm& communicator() const { return comm_;}
+
+  private:
+    mutable MPI_Comm comm_;
+    int rank_;
+    std::size_t global_size_;
+    std::size_t block_size_;
+};
+
+
+
+
+#endif 
\ No newline at end of file

examples/05_mpi/errors.cpp

+// Example codes for HPC course
+// (c) 2012 Matthias Troyer, ETH Zurich
+
+#include <iostream>
+#include <string>
+#include <stdexcept>
+#include <mpi.h>
+
+// throw a std::runtime_error if there was an MPI error
+void check_error(int err)
+{
+  if (err != MPI_SUCCESS) {
+    // get the error text for the error code
+    int len = MPI_MAX_ERROR_STRING;
+    char txt[MPI_MAX_ERROR_STRING];
+    MPI_Error_string(err, txt, &len);
+    // and throw an exception
+    throw std::runtime_error(txt);
+  }
+}
+
+int main(int argc, char** argv) {
+  
+  MPI_Init(&argc, &argv);
+  int num;
+  
+  MPI_Comm_rank(MPI_COMM_WORLD,&num);
+  
+  // tell MPI to return an error code instead of aborting
+  MPI_Errhandler_set(MPI_COMM_WORLD, MPI_ERRORS_RETURN);