working mpi convolution (2021-03-12)

MPI/03_convolution/convolution.jl

 using MPI
 
+"""
+    gaussianstencil(T)
+
+Create a Gaussian blur stencil of floating point type `T`
+
+For more information see
+<https://en.wikipedia.org/wiki/Kernel_(image_processing)>
+"""
+@inline function gaussianstencil(T)
+  SMatrix{5, 5, T, 25}(1,  4,  6,   4,  1,
+                       4, 16, 24,  16,  4,
+                       6, 24, 36,  24,  6,
+                       4, 16, 24,  16,  4,
+                       1,  4,  6,   4,  1) / 256
+end
+
 """
     split(N)
 
@@ -44,4 +60,192 @@ function partition(mpirank, mpisize, Nx, Ny)
     return (rank = (x = mpirank_x, y = mpirank_y),
             size = (x = mpisize_x, y = mpisize_y),
             rng  = (x = ix1:ix2  , y = iy1:iy2  ))
+end
+
+"""
+    mpi_applyblur!(B, A, mpicomm, part)
+
+Apply the Gaussian blur stencil to the matrix `A` and store the result in the
+matrix `B`
+"""
+function mpi_applyblur!(B, A, mpicomm, part)
+  # Get the floating point type we are using
+  T = eltype(A)
+
+  # Get the gaussianstencil
+  S = gaussianstencil(T)
+
+  # Get the size of the image
+  Lx = size(A, 1) - 4
+  Ly = size(A, 2) - 4
+
+  # Make sure that A and B are compatible sizes
+  @assert size(A) == size(B)
+
+  # Get the half width of the stencil
+  W = div(size(S, 1), 2)
+
+  # Get that stencil size is odd and same in both dimensions
+  @assert (2W+1, 2W+1) == size(S)
+
+  # Question: Which of our edges are on parallel boundaries and which are on
+  # physical boundaries?
+  I1 = part.rank.x == 0 ? 3 : 1
+  I2 = part.rank.x == part.size.x - 1 ? Lx + 2 : Lx + 4
+  J1 = part.rank.y == 0 ? 3 : 1
+  J2 = part.rank.y == part.size.y - 1 ? Ly + 2 : Ly + 4
+
+  # Communication
+  
+  send_req = fill(MPI.REQUEST_NULL, 8)
+  recv_req = fill(MPI.REQUEST_NULL, 8)
+
+  # West
+  if part.rank.x > 0
+    nabr_x = part.rank.x - 1
+    nabr_y = part.rank.y
+    nabr_rank = nabr_x + nabr_y * part.size.x
+    recv_req[1] = MPI.Irecv!(@view(A[1:2, 3:(Ly + 2)]),
+                             nabr_rank,
+                             777,
+                             mpicomm)
+    send_req[1] = MPI.Isend( @view(A[3:4, 3:(Ly + 2)]),
+                              nabr_rank,
+                              777,
+                              mpicomm)
+  end
+
+  # East
+  if part.rank.x < part.size.x - 1
+    nabr_x = part.rank.x + 1
+    nabr_y = part.rank.y
+    nabr_rank = nabr_x + nabr_y * part.size.x
+    recv_req[2] = MPI.Irecv!(@view(A[(Lx + 3):(Lx + 4), 3:(Ly + 2)]),
+                             nabr_rank,
+                             777,
+                             mpicomm)
+    send_req[2] = MPI.Isend( @view(A[(Lx + 1):(Lx + 2), 3:(Ly + 2)]),
+                              nabr_rank,
+                              777,
+                              mpicomm)
+  end
+
+  # South
+  if part.rank.y > 0
+    nabr_x = part.rank.x
+    nabr_y = part.rank.y - 1
+    nabr_rank = nabr_x + nabr_y * part.size.x
+    recv_req[3] = MPI.Irecv!(@view(A[3:(Lx + 2), 1:2]),
+                             nabr_rank,
+                             777,
+                             mpicomm)
+    send_req[3] = MPI.Isend( @view(A[3:(Lx + 2), 3:4]),
+                              nabr_rank,
+                              777,
+                              mpicomm)
+  end
+
+  # North
+  if part.rank.y < part.size.y - 1
+    nabr_x = part.rank.x
+    nabr_y = part.rank.y + 1
+    nabr_rank = nabr_x + nabr_y * part.size.x
+    recv_req[4] = MPI.Irecv!(@view(A[3:(Lx + 2), (Ly + 3):(Ly + 4)]),
+                             nabr_rank,
+                             777,
+                             mpicomm)
+    send_req[4] = MPI.Isend( @view(A[3:(Lx + 2), (Ly + 1):(Ly + 2)]),
+                              nabr_rank,
+                              777,
+                              mpicomm)
+  end
+
+  # South-West
+  if part.rank.x > 0 && part.rank.y > 0
+    nabr_x = part.rank.x - 1
+    nabr_y = part.rank.y - 1
+    nabr_rank = nabr_x + nabr_y * part.size.x
+    recv_req[5] = MPI.Irecv!(@view(A[1:2, 1:2]),
+                             nabr_rank,
+                             777,
+                             mpicomm)
+    send_req[5] = MPI.Isend( @view(A[3:4, 3:4]),
+                              nabr_rank,
+                              777,
+                              mpicomm)
+  end
+
+  # South-East
+  if part.rank.x < part.size.x - 1 && part.rank.y > 0
+    nabr_x = part.rank.x + 1
+    nabr_y = part.rank.y - 1
+    nabr_rank = nabr_x + nabr_y * part.size.x
+    recv_req[6] = MPI.Irecv!(@view(A[(Lx + 3):(Lx + 4), 1:2]),
+                             nabr_rank,
+                             777,
+                             mpicomm)
+    send_req[6] = MPI.Isend( @view(A[(Lx + 1):(Lx + 2), 3:4]),
+                              nabr_rank,
+                              777,
+                              mpicomm)
+  end
+
+  # North-West
+  if part.rank.x > 0 && part.rank.y < part.size.y - 1
+    nabr_x = part.rank.x - 1
+    nabr_y = part.rank.y + 1
+    nabr_rank = nabr_x + nabr_y * part.size.x
+    recv_req[7] = MPI.Irecv!(@view(A[1:2, (Ly + 3):(Ly + 4)]),
+                             nabr_rank,
+                             777,
+                             mpicomm)
+    send_req[7] = MPI.Isend( @view(A[3:4, (Ly + 1):(Ly + 2)]),
+                              nabr_rank,
+                              777,
+                              mpicomm)
+  end
+
+  # North-East
+  if part.rank.x < part.size.x - 1 && part.rank.y < part.size.y - 1
+    nabr_x = part.rank.x + 1
+    nabr_y = part.rank.y + 1
+    nabr_rank = nabr_x + nabr_y * part.size.x
+    recv_req[8] = MPI.Irecv!(@view(A[(Lx + 3):(Lx + 4), (Ly + 3):(Ly + 4)]),
+                             nabr_rank,
+                             777,
+                             mpicomm)
+    send_req[8] = MPI.Isend( @view(A[(Lx + 1):(Lx + 2), (Ly + 1):(Ly + 2)]),
+                              nabr_rank,
+                              777,
+                              mpicomm)
+  end
+
+  MPI.Waitall!(recv_req)
+
+  # Loop over (i,j) pixels of the output image
+  for j = 3:(Ly + 2)
+    for i = 3:(Lx + 2)
+      # storage for application of the stencil
+      b = -zero(T)
+
+      # Loop over the stencil indices (k + W + 1, l + W + 1)
+      for l = -W:W
+        # (I, J) are the input image pixels
+        # (max and min prevent going outside the image)
+        J = max(J1, min(J2, j + l))
+        for k = -W:W
+          I = max(I1, min(I2, i + k))
+
+          # apply the stencil to the pixel
+          @inbounds b += S[k + W + 1, l + W + 1] * A[I, J]
+        end
+      end
+
+      # Save the value to the output image
+      @inbounds B[i, j] = b
+    end
+  end
+
+  MPI.Waitall!(send_req)
+
 end
\ No newline at end of file

MPI/03_convolution/convolution_test.jl

 include("serial_convolution.jl")
 include("convolution.jl")
 
+using Random
+Random.seed!(7)
+
 MPI.Init()
 
 let
@@ -16,8 +19,20 @@ let
 
     part = partition(mpirank, mpisize, Nx, Ny)
 
-    # TODO: allocate our local data view
-    # TODO: Call a parallel applyblur!
-    # TODO: check result
+    # allocate our local data view
+    Lx = length(part.rng.x)
+    Ly = length(part.rng.y)
+    mpi_A = zeros(Lx + 4, Ly + 4)
+    mpi_B = similar(mpi_A)
+
+    # Copy the real data into our local data arrays
+    mpi_A[3:(Lx + 3 - 1), 3:(Ly + 3 - 1)] .= A[part.rng.x, part.rng.y]
+
+    # Call a parallel applyblur!
+    mpi_applyblur!(mpi_B, mpi_A, mpicomm, part)
+
+    # check result
+    @assert mpi_B[3:(Lx + 3 - 1), 3:(Ly + 3 - 1)] == B[part.rng.x, part.rng.y]
+
     # TODO: parallel timings
 end
\ No newline at end of file