From 9ad56cd68572bb612c57a4ea9cd10628137ef15e Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Tue, 20 Aug 2024 15:23:15 +0200 Subject: [PATCH] Enhancements in Jacobi notebook --- docs/make.jl | 2 +- docs/src/solutions_for_all_notebooks.md | 106 +- notebooks/figures/fig_jacobi.svg | 2174 ++++++++++++++++++++--- notebooks/jacobi_method.ipynb | 298 ++-- 4 files changed, 2162 insertions(+), 418 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 1c23f39..14944b2 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -122,7 +122,7 @@ makedocs(; #"Distributed computing with MPI" => "mpi_tutorial.md", "Matrix-matrix multiplication"=>"matrix_matrix.md", "MPI (point-to-point)" => "julia_mpi.md", - #"Jacobi method" => "jacobi_method.md", + "Jacobi method" => "jacobi_method.md", #"All pairs of shortest paths" => "asp.md", #"Gaussian elimination" => "LEQ.md", #"Traveling salesperson problem" => "tsp.md", diff --git a/docs/src/solutions_for_all_notebooks.md b/docs/src/solutions_for_all_notebooks.md index ef41b4f..2409ffc 100644 --- a/docs/src/solutions_for_all_notebooks.md +++ b/docs/src/solutions_for_all_notebooks.md @@ -167,50 +167,70 @@ end ### Exercise 1 ```julia -@everywhere workers() begin - using MPI - comm = MPI.Comm_dup(MPI.COMM_WORLD) - function jacobi_mpi(n,niters) - nranks = MPI.Comm_size(comm) - rank = MPI.Comm_rank(comm) - if mod(n,nranks) != 0 - println("n must be a multiple of nranks") - MPI.Abort(comm,1) - end - n_own = div(n,nranks) - u = zeros(n_own+2) - u[1] = -1 - u[end] = 1 - u_new = copy(u) - for t in 1:niters - reqs = MPI.Request[] - if rank != 0 - neig_rank = rank-1 - req = MPI.Isend(view(u,2:2),comm,dest=neig_rank,tag=0) - push!(reqs,req) - req = MPI.Irecv!(view(u,1:1),comm,source=neig_rank,tag=0) - push!(reqs,req) - end - if rank != (nranks-1) - neig_rank = rank+1 - s = n_own+1 - r = n_own+2 - req = MPI.Isend(view(u,s:s),comm,dest=neig_rank,tag=0) - push!(reqs,req) - req = MPI.Irecv!(view(u,r:r),comm,source=neig_rank,tag=0) - push!(reqs,req) - end - for i in 3:n_own - u_new[i] = 0.5*(u[i-1]+u[i+1]) - end - MPI.Waitall(reqs) - for i in (2,n_own+1) - u_new[i] = 0.5*(u[i-1]+u[i+1]) - end - u, u_new = u_new, u - end - return u +function jacobi_mpi(n,niters) + comm = MPI.COMM_WORLD + nranks = MPI.Comm_size(comm) + rank = MPI.Comm_rank(comm) + if mod(n,nranks) != 0 + println("n must be a multiple of nranks") + MPI.Abort(comm,1) end + load = div(n,nranks) + u = zeros(load+2) + u[1] = -1 + u[end] = 1 + u_new = copy(u) + for t in 1:niters + reqs = MPI.Request[] + if rank != 0 + neig_rank = rank-1 + req = MPI.Isend(view(u,2:2),comm,dest=neig_rank,tag=0) + push!(reqs,req) + req = MPI.Irecv!(view(u,1:1),comm,source=neig_rank,tag=0) + push!(reqs,req) + end + if rank != (nranks-1) + neig_rank = rank+1 + s = load+1 + r = load+2 + req = MPI.Isend(view(u,s:s),comm,dest=neig_rank,tag=0) + push!(reqs,req) + req = MPI.Irecv!(view(u,r:r),comm,source=neig_rank,tag=0) + push!(reqs,req) + end + for i in 3:load + u_new[i] = 0.5*(u[i-1]+u[i+1]) + end + MPI.Waitall(reqs) + for i in (2,load+1) + u_new[i] = 0.5*(u[i-1]+u[i+1]) + end + u, u_new = u_new, u + + end + # Gather the results + if rank !=0 + lb = 2 + ub = load+1 + MPI.Send(view(u,lb:ub),comm,dest=0) + u_all = zeros(0) # This will nevel be used + else + u_all = zeros(n+2) + # Set boundary + u_all[1] = -1 + u_all[end] = 1 + # Set data for rank 0 + lb = 2 + ub = load+1 + u_all[lb:ub] = view(u,lb:ub) + # Set data for other ranks + for other_rank in 1:(nranks-1) + lb += load + ub += load + MPI.Recv!(view(u_all,lb:ub),comm;source=other_rank) + end + end + return u_all end ``` diff --git a/notebooks/figures/fig_jacobi.svg b/notebooks/figures/fig_jacobi.svg index 52b6518..14fa8e2 100644 --- a/notebooks/figures/fig_jacobi.svg +++ b/notebooks/figures/fig_jacobi.svg @@ -1,24 +1,24 @@ uu +u_newu_new +u_new[i] = u[i-1]+u[i]+u[1+1]u_new[i] = u[i-1]+u[i]+u[1+1] +ii +ii +i+1i+1 +i-1i-1 +ii +i+1i+1 +i-1i-1 +ii +ii +i+1i+1 +i-1i-1 +ii +uu +u_newu_new +uu +u_newu_new +uu +u_newu_new +uu +uu +u_newu_new +uu +u_newu_new +ii +i-1i-1 +i+1i+1 +ii +chnl_next_sndchnl_next_snd +chnl_next_rcvchnl_next_rcv +chnl_prev_rcvchnl_prev_rcv +chnl_prev_sndchnl_prev_snd +xx +uu +00 +LL +-1-1 +11 +00 +LL +n+2 pointsn+2 points +uu +?? +-1-1 +?? +?? +?? +?? +?? +?? +?? +?? +?? +11 +uu +u_newu_new +ii +i-1i-1 +i+1i+1 +ii +ii +i+1i+1 +i-1i-1 +ii +uu +u_newu_new +kk +ii +jj +uu +uu +uu +u_newu_new +ii +i-1i-1 +i+1i+1 +ii +kk +ii +jj +kk +ii +jj +uu +?? +-1-1 +?? +?? +?? +?? +?? +?? +?? +?? +?? +11 +uu +u_newu_new +?? +-1-1 +?? +?? +?? +?? +?? +?? +?? +?? +?? +11 +-1-1 +-1-1 +11 +11 +ss +rr +ss +uu +u_newu_new +rr +ss +rr +ss +uu +u_newu_new +rr +r-1r-1 +s+1s+1 +uu +uu +u_newu_new +communicationcommunication +computationcomputation +uu +u_newu_new +i-1i-1 +i+1i+1 +ii +i+1i+1 +-1-1 +-1-1 +-1-1 +11 +11 +11 +uu +uu +u_newu_new +communicationcommunication +computationcomputation +-1-1 +-1-1 +-1-1 +11 +11 +11 +(i+1,j)(i+1,j) +(i,j)(i,j) +(i-1,j)(i-1,j) +(i,j+1)(i,j+1) +(i,j-1)(i,j-1) + +uu +u_newu_new +locallocal +remoteremote +remoteremote +kk +ii +jj +Graph GGraph G +...... +...... +...... +...... +...... +...... +MasterMaster +WorkersWorkers +maxhopsmaxhops +uu +u_newu_new +rank 0rank 0 +rank 1rank 1 +rank 2rank 2 +rootroot +rand(1:10)rand(1:10) +uu +u_newu_new +locallocal +remoteremote +remoteremote +locallocal +remoteremote +Distance matrix CDistance matrix C + + id="path3" />i+1 +i-1 +i +u +u_new +i+1 +i-1 +i +u +u_new + \ No newline at end of file diff --git a/notebooks/jacobi_method.ipynb b/notebooks/jacobi_method.ipynb index ad8c6bb..1c91c5f 100644 --- a/notebooks/jacobi_method.ipynb +++ b/notebooks/jacobi_method.ipynb @@ -62,7 +62,8 @@ "gauss_seidel_1_check(answer) = answer_checker(answer,\"c\")\n", "jacobi_1_check(answer) = answer_checker(answer, \"d\")\n", "jacobi_2_check(answer) = answer_checker(answer, \"b\")\n", - "jacobi_3_check(answer) = answer_checker(answer, \"c\")" + "jacobi_3_check(answer) = answer_checker(answer, \"c\")\n", + "println(\"🥳 Well done! \")" ] }, { @@ -145,7 +146,7 @@ "metadata": {}, "source": [ "
\n", - "Note: In our version of the jacobi method, we return after a given number of iterations. Other stopping criteria are possible. For instance, iterate until the difference between u and u_new is below a tolerance. \n", + "Note: In our version of the jacobi method, we return after a given number of iterations. Other stopping criteria are possible. For instance, iterate until the difference between u and u_new is below a tolerance.\n", "
" ] }, @@ -178,7 +179,7 @@ "id": "798968b1", "metadata": {}, "source": [ - "### The Gauss-Seidel method\n", + "### The Gauss-Seidel method \n", "\n", "The usage of `u_new` seems a bit unnecessary at first sight, right?. If we remove it, we get another method called Gauss-Seidel.\n", "\n" @@ -275,7 +276,7 @@ "source": [ "### Parallelization strategy\n", "\n", - "Remember that a sufficiently large grain size is needed to achieve performance in a distributed algorithm. For Jacobi, one could update each entry of vector `u_new` in a different worker, but this would not be efficient. Instead, we use a parallelization strategy with a larger grain size that is analogous to the algorithm 3 we studied for the matrix-matrix multiplication:\n", + "Remember that a sufficiently large grain size is needed to achieve performance in a distributed algorithm. For Jacobi, one could update each entry of vector `u_new` in a different process, but this would not be efficient. Instead, we use a parallelization strategy with a larger grain size that is analogous to the algorithm 3 we studied for the matrix-matrix multiplication:\n", "\n", "- Each worker updates a consecutive section of the array `u_new` \n", "\n", @@ -299,16 +300,68 @@ }, { "cell_type": "markdown", - "id": "3f90d701", + "id": "22dc6e54", "metadata": {}, "source": [ "### Data dependencies\n", "\n", "Recall the Jacobi update:\n", "\n", - "`u_new[i] = 0.5*(u[i-1]+u[i+1])`\n", - "\n", - "Thus, in order to update the local entries in `u_new`, we also need remote entries of vector `u` located in neighboring processes. Figure below shows the entries of `u` needed to update the local entries of `u_new` in a particular process (CPU 2)." + "`u_new[i] = 0.5*(u[i-1]+u[i+1])`" + ] + }, + { + "cell_type": "markdown", + "id": "ba4113af", + "metadata": {}, + "source": [ + "Note that an entry in the interior of the locally stored vector can be updated using local data only. For this one, communication is not needed." + ] + }, + { + "attachments": { + "f1.png": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6gAAADyCAYAAABAvOgkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAewgAAHsIBbtB1PgAAABl0RVh0U29mdHdhcmUAd3d3Lmlua3NjYXBlLm9yZ5vuPBoAABgzSURBVHic7d19mF11YSfw750kBBIigWJiADUEBOXdRYvgCwtCX/VhH0TbLRVrdUlod/Gpq0+7a1uoL1316RarqzSpZXe77q61Pmu73bZYEFOtFUSs4BsCQhCEEEyUd4aQ3P3jnGHOTO6duTOZmfObyefzPPfJ79x7zr1fZm5O7pfzO+d2Nm7c2A0AAAC0bKjtAAAAAJAoqAAAABRi8bjln+50OlvaCAIAAMC+pdvtrk3y2ZHlMQW10+lsufjii2+b61AAAADsezZt2pRud/SySKb4AgAAUAQFFQAAgCIoqAAAABRBQQUAAKAICioAAABFUFABAAAogoIKAABAERRUAAAAiqCgAgAAUAQFFQAAgCIoqAAAABRBQQUAAKAICioAAABFUFABAAAowuK2AwAAPb05yZIkVyf5fstZAGBOKKgAUKY/SrIiyWsz/YK6JMmJSU5N8pIkB9f3/26SW/c2YGH+KMnLk/z3JB9p4fUXp/p9JcnjSYZbyAAw7ymoALAw/U6SdyVZ2uOxD89xlrlwdKoifu0cvd6aJOemKv6nJjklybL6sfVJNs1RDoAFRUEFgDK9LdUR0Fumuf2aVOW0m+T2JPcmOXtmopHkgizMog/QKgUVAMr0X/dy+08m+VSSryV5OMlLk3xlb0MV7Oenud2FqaZBfzrJV6e47ZYkN9W3byf5y2lmAKCmoAJAmdamutr+1lTnNE7VF2Y0zcL1r1IdDf1uplZQP5qx57oun8lQAPsqXzMDAGW6Jcn3YlruoN6c5P1JfnaOXm/3HL0OwD7FEVQAYCG4IMnP1eO/azMIANPnCCoAAABFUFABAAAogoIKAABAEZyDCgDzx7IkZ/R57LYk35/DLPPJLyZ5S5/HTqr/fEeSX+qzzvokd850KAD2pKACwPxxWJJr+jz29iRXzGGW+eTIJOdMss5x9a2XA2c2DgD9KKgAMH88nuTaPo85etrfZ5Lc0eextyd5WZKPJdncZ527ZyETAD0oqAAwf9yX5Ny2Q8xDt9a3Xt5Q//nVJH8xN3EA6MdFkgAAACiCggoAAEARTPEFgIVpbZINjeU1jfGGJK9pLH88/c/RBIA5o6ACwML03CS/2eexC8ctXxMFdaoOTvK9xnKnMb4iyfsby7+VZNNchAKY7xRUACjTFUmWZvrF8Z4kHxhwXVepnbpOqpLay7L6NmL/2Y8DsDAoqABQpsv2cvstqY7c7SvekqoU/niK2/3bVEeaH5zidj9OctSA626f4nMD7LMUVABgIdg6ze0emOZ2u5PcOc1tAejDVXwBAAAogoIKAABAERRUAAAAiqCgAgAAUAQFFQAAgCIoqAAAABRBQQUAAKAICioAAABFUFABAAAogoIKAABAERRUAAAAirC47QDQy5VXXnn40NDQG9rOMZGrrrrqVUk6SXLXXXfdsm3bth+1HGm8E5McUo/vSXJni1l6OSLJUfX40SQ3tZhlEFcmebLtEAzkDUkObzvEBJ6fZG09fijJ19uL0tOzkry4sfzFJLtbytLPKzP6P9m/kWRHi1km87Uk/9B2CAby/CTntx1iAp0kr2osfz3VPqQkJydZWY/vTrKlvSh7Wrt27fNXr169NklWrlz50Pnnn1/a/neM+++//yOXX375023nmGsKKkVatGjR0d1u9w/bzjGRG2+8Mbt3P/OZ7bw2sywQF7YdYBJ/FgV1vrg0ycvbDjEFF7UdYBKl799Kz/ehKKjzxYuSFP3ZY5zS3/vF2bJlS7Zs2ZIkWbduXVL4/nfNmjUbk+xzBdUUXwAAAIrgCCrzxS1tB+jhpJHB0qVLtw4PD29rM0wPx2X07/gTSW5vMUsva1NNJUySbqppeiXZP8kxbYdgrz1Q30pyVJLl1XBod7KssOn3wwclO5/duOObKW+K70mjw/13JIsLm+L7xGHJrmVtp2CvlfbZo5Pq9J0R9yX5YUtZ+jk+yaJ6/FiS77WYZQ9Lly5dNzw8fGCSdLvd3an2byVZluTotkO0TUFlvnjJ+vXrd7YdomloaKg7Mj799NM/unnz5ve2maeHB5McWo+/leSlLWbp5a+TvKYeP5nqvJWSvDjVuWPMb3+c5PK2Q4yzOcmZ1XD5E8nDhU1vP+P1yZff0bjj5UkebitNH93R4Wn/J9n8J+1F6WXVFcmDZ7Sdgr2yK+X9u7QsVekb8YEkH24pSz8/TnJQPf7nVOeLF+OEE0647qabbjorSXbu3PnY+vXri/odb9q06Yxut/ultnO0zRRfAAAAiqCgAgAAUAQFFQAAgCIoqAAAABRBQQUAAKAICioAAABFUFABAAAogoIKAABAERRUAAAAiqCgAgAAUAQFFQAAgCIoqAAAABRBQQUAAKAICioAAABFUFABAAAogoIKAABAERRUAAAAiqCgAgAAUAQFFQAAgCIoqAvHjUm+Wt+OG3CbzzS2+alZygUAADCQxW0HYMacmqRTj5cPuM0JSY6ux4fMeCIAAIApcAQVAACAIiioAAAAFEFBBQAAoAgKKgAAAEVQUAEAACiCggoAAEARFFQAAACKoKACAABQBAUVAACAIiio+7YlbQcAAAAYoaAuHMON8dIB1u8kOXSWsgAAAEyZgrpwPNwYP2eA9Y9PsnyWsgAAAEyZgrpw3N4Ynz7A+r86W0EAAACmQ0FdOG5ojC9M8qwJ1v3JJL8+u3EAAACmZnHbAZgxn0jy9nq8OsmfJ/nlJNsb63SS/OskH0t1gaThDHa+aus6nc5LN27cuKvtHE2XXHLJM+OtW7cekeS09tL01Pz7vTzl5Tu4MR5KefmOaTsAM+LwlPfeWjE63DWUvO349qL0sn38aSKnJnm8jSSD+eGzy/sZPnVg2wmYEaXtO8Z/Znt+ysu4qDE+MIXle+SRR57Z/3Y6nUUbN24sKl+32z2u7QwlUFAXjn9OVUp/oV7+mSR3Jdmc5P5UZeCMVB/WkuSyJG9McvScppymbrf7pbYzTOTWW29dn2R92zkm8KIk17cdYgJLU3Y+5q+31rdCPX5A8uGr2k4xievaDjCxb51X3WBGLUr5/y69PaMHJ0p0Sgr7Gd52223PjJcsWbIsheWjoqAuLBtS/d+0l9XLK5K8dtw63SQfTPKeVAUVAACgCArqwvLjJGcluTTVEYMXNB7bleTqJB9Kcm193z8l2VKPt85NRAAAgN4U1IXnyVRHSD+Y5Nmpzkd9PMm9SZ4at+6b5jba4O67774vrlu3ruhzeF73utcte/rppztJ8o1vfOPJO+64o6hzZFNNmx35O74ze/7+27Y4o+fT7Er13i3ZY20HYGDnZOx5UKVZkmS/elzie38oyQGN5RLf+8tSXVchqX5+pe1/m3a2HYCBXZPqvMmSNb8i8Ikku9sK0sf+Gd3/PpXC3v+nnHLKkiOPPHK/JFmxYsWuAw44oLT97xgXXXRRwef/z57Oxo0bu88sdDrHXnzxxbdNtAEAAADMhE2bNh3T7Xa/O7Lsa2YAAAAogoIKAABAERRUAAAAiqCgAgAAUAQFFQAAgCIoqAAAABRBQQUAAKAICioAAABFUFABAAAogoIKAABAERRUAAAAiqCgAgAAUAQFFQAAgCIoqAAAABRBQQUAAKAICioAAABFUFABAAAogoIKAABAERRUAAAAiqCgAgAAUAQFFQAAgCIoqAAAABRBQQUAAKAICioAAABFUFABAAAogoIKAABAERRUAAAAiqCgAgAAUAQFFQAAgCIoqAAAABRBQQUAAKAICioAAABFUFABAAAogoIKAABAERRUAAAAiqCgAgCz7VNJrktydttBACibggoAzLZXJDkryeq2gwBQNgUVAACAIixuOwAAsOCdkmRRkh+1HQSAsimoAMBs29Z2AADmBwUVAJhtb05yYJKrk9zechYACqagAgCz7X1J1iT5pSioAEzARZIAAAAogoIKAABAERRUAAAAiqCgAgAAUAQFFQAAgCIoqAAAABRBQQUAAKAICioAAABFUFABAAAogoIKAABAERRUAAAAiqCgAgAAUITFbQcAABa885IsTXJr20EAKJuCCgDMthvbDgDA/GCKLwAAAEVQUAEAACiCggoAAEARFFQAAACKoKACAABQBAUVAACAIiioAAAAFEFBBQAAoAiL2w4AvWzatOmUbrf7sbZzTOR973vf8bt37+4kyfbt27c88cQTj7adaZy1SQ6sxz9MsrW9KD0dmuQ59fjJJHe0mGUQP5Pk4bZDMJArk5zcdogJrKpvSfJ4kjtbzNLL8iRHNpa/laTbUpZ+jk/Sqcd3J3mkxSyT+WSSD7cdgoGcluSKtkNMoJPqvT/iriSPtZSlnyNT7UOS5MEkD7SYZQ8HHXTQqhUrVqxKkjVr1jz+1re+tbT97xjDw8NnXXrppcNt55hrCiqlOijJ6W2HmMi9996b3bt3jyye2GaWAazI2A+cpVmR5Nlth5jEkrYDMLATU/j+o2FFktVth5jEy9oOMIkT2g4wiRvaDsDADs782XckyUltB5jEiiTr2g7R9NBDD+Whhx5Kkuy3337F73+XLl26qO0MbTDFFwAAgCI4gsq80Ol0LtvdOFxZiPeMDFatWvW5bdu2bW4xSy//IcmyevyDJH/cYpZefjnJsfX46SS/12KWXtYk+bW2Q7DXPp/kurZDjPOWVFPwk2Q4yXvbi9LTaUleM7p4wlXJop2tpenp5vWNhetS/Z5L8sYkx7Qdgr2yO8llbYcYZ0mS320s/02S61vK0s+7kuxfj+9O8vEWs+xh9erVb37ggQfWJcmuXbue6na775lsm7nU6XSem+TitnO0TUFlXuh2u/9pw4YNRX1AGhoaemandtxxx23etm1baR8y35bRgnp/yvwQPFJQd6a8fC+OgroQfCHlvbfOSdkF9dczpqB+4hPJyYWd59ZpFtQvpryf4elRUOe7bsp7Xy3L2IL69ynv/OZ3ZLSg3pPCfoZHHHHE2Y2COrxhw4ai8m3atOmMbre7zxdUU3wBAAAogoIKAABAERRUAAAAiqCgAgAAUAQFFQAAgCIoqAAAABRBQQUAAKAICioAAABFUFABAAAogoIKAABAERRUAAAAirC47QDjvCPJonr8iSQ/GGCbtyQ5tB7/TZJvzkKuES9Icn49fjzJR+rxs5L8Yv3Y85KsTLItyQ1J/luSL0/jtY6un+/sJEck+YkkD6X6mVyT5H8l+f4E2/+LJOfW4/uS/I8J1n1hkvMay19N8rkJ1n91kpfU428n+esJ1gUAABhIaQX1vUmW1uN/zGAF9e1JjqvHWzO7BfW4JO+vx9tSFdRXJPnfqUpk05okJyf5N0muTPLvkuwe4DUOSvLBJL+aPX8/z0lybKrSelmSP6j/7PW8BzSyPprkz5M81ec135TktxrLX8jEBfV3kpxZj98ZBRUAAJgBpvjunVemKnIj5fSRJHemOro6opPk11KV78kckaqYX5yx5fSeJF+rn7tb37d/kt9OdSS11+/xK3WeJDkwyWkTvO6rxy2/LMnyPusuqx8fcd0EzwsAADAwBXX6lif5dJL9UpXEU1JN9T2q/vNnk2xprP/OJOsmeL79k/zfJCfUy8NJfi/J4ammDZ9aP/dzk2xsbPcL9XOPtzPVkdAR40voiINTTQdu2i9V+e7lFRk9yr09ydf7rAcAADAlCur0LU+yKsl/THJhkpsbj+1KcnWSczJ6NHVxkl+Z4PkuT/LievxYkrPq++4bt94PkmxIdfS0ue2qHs/ZnKbbr6CemdHzfr80wPrN+z+fwaYtAwAATEpB3TvXZfQ8z16+l+R/NpZf1We9lammAY94Zya/sNLvp7qYUVIdfb24T74Rp6Wa6jtes3C+O8mTPe7vt/5E56kCAABMiYK6dz6U0XNC+2mWxOP7rHNBkhX1eHuSjw/w2t0kf9JYPrfHOrekuphTkixJ74I8UjgfTbI5o0dRT87o1ZFHHJLRo7yJ808BAIAZpKDunS9MvsqY81BXprpo0nhnNsbXpjp/dBA3NMYv7fF4N9U03BHnjHv88CQvqsdfTHWV35GjokOpphk3nZXR98w9SW4bMCcAAMCkSvuamfnkkVTfSzqZxxrjxakuMPTkuHVOboxXJ/nNATMc1BgfkOoo7CPj1vlcqgspJXtO2z27Mb62sf6Ic5L8RWPZ9F4AAGDWKKjT98SA642fAtzrCOpPNMb/sr5Nx8r0LqgjTkx1MaWRab+9CudNSX5cP9f4QqugAgAAs8YU3zL0unjRdPT6fd6Z0WnGnYw9ajpSOLelOl81qa5AvLkeH5VkbT1+bpJjGts6/xQAAJhRC6GgLoT/huZU4XeluhjRdG739Hn+8dN2k+TYJEfU489n7JHeXl9P0zx6+p3s+fU3AAAAe6W0crerMV7Ud62xVs5GkDm2rTFek+RH07z1+07SyQrntRlrsvVN7wUAAGZcaQX14cb44AHWPyTVRYXmu+sb4/FX2p0J12X0COnaVFN3Jyqc30nyg3p8dqr3ydkTrA8AALDXSiuo9zfGJw6w/nnpfdGh+eazjfELM/2LJPXzQJJvNZbPbbzG95Lc1WObkXNMVyd5fZLD6uXmOaoAAAAzprSC+rXG+PWZuHyuSPLbsxtnzvxtklsbyx/JzF04aUTzqOdvpDr6PP7+fuu/uzEeucovAADAjCqtoP5lY3xSknf0We+QJH+VZN2sJ5obu5K8M6PTcE9I8vdJnjfAtsck+WiqQj+RZuE8ps/9g6zv6r0AAMCsKO17UK9ONRX1+Hr5g0lemeTTqa4auzLJGUl+JdU5qv+Y5DlJjp7roLPg/6U6UnlZvXx6ktuSfCrJNamm4T6R6sjq4UlenGqq7sn1+l+Z5Pn/IcnTGfs7353+hfPeJN9NdbXfJuefAgAAs6K0gvp0kotSlaaD6vteW9/G+06SC7KwjuhdnuSHSf4wyZIkS5O8sb7trYeT3Jiq+I64uX69fj6XsQX1ySRfmoEsAAAAeyhtim9SnYd6ZpIb+jz+VJI/TfKTqS7+s9D8l1QXiPqzVIVwIo+lOn/1olRHmSczvsxPdjR0/OP/lOooLgAAwIwr7QjqiJtTHek7KcnLU03jfTzJ3anOzdzeWPf0jH5n6mOznOtvM3pxoe5EKzZ8p7FNMnnpTKqptW9Ksj7VlOYXJjk0yQGpLlB0X73O11MV9kG9O8l/biw/Psn6f5Wx2Yen8FoAAABTUmpBTaoCeHN9m8jDkzw+k3Ym+dEUt9k1jW1GPJnqqOdMTWN+KlMrtHuTHQAAYEpKnOILAADAPkhBBQAAoAglT/GFpg9s3Lhxd9shmi655JJnxrfccstPp/oapJIsb4yfm+QP2grSx3GN8ZKUl29V2wGYET+V6uu5SnJUY7x/ynvvnzJ28YINSWdnO1EGcm7K+x2/qO0A7LWhlPd3c/zn9vOSPK+NIBPYvzE+MoX9DLds2fLMV1MuXrx46caNG4vKl+SwtgOUYCEX1BWZmf++XZnb81zp7TfaDjCRHTt2vCLJK9rOMYHVSf592yEmsDhl52P+Oj1jv16rNPul+Pf+HW9oO8EkzqhvMJM6Kf7vZs6ub6U6PIX9DLdvH73O6tDQUHH732530GuwLmwLuaB+JsmrZ+B5bkrykhl4HgAAACawkAsq89iuXbu2Dw0N/V3bOSZy5JFHntrtdjtJsm3bttsfffTRh9rONM4LkhxUj7cmubfFLL2sTjX1OKm+8ujbLWYZxFSugE27vpyyZ74cltFpXI8mubXFLL0cmOrrzUbclMG/Wm2unJrqCFeS3J6ktP1vU+n7NkY9mKTkzx6dVO/9Ed9N8khLWfo5NtUsxiS5P8kPWsyyh0MOOWTNypUrD0+SVatWlbj/HWN4eHhX2xnasJAL6seTXDMDz7N1Bp6DKbrkkku+meTn2s4BzEvvbDsAMC/dFJ89FrQdO3Zkx44dSZI777wz119/fcuJ6GUhF9RPth0AAACAwfmaGQAAAIqgoAIAAFAEBRUAAIAiKKgAAAAUQUEFAACgCAoqAAAARVBQAQAAKIKCCgAAQBEUVAAAAIqgoAIAAFAEBRUAAIAiKKgAAAAUQUEFAACgCAoqAAAARVBQAQAAKIKCCgAAQBEUVAAAAIqgoAIAAFAEBRUAAIAiKKgAAAAUQUEFAACgCAoqAAAARVBQAQAAKIKCCgAAQBEUVAAAAIqgoAIAAFAEBRUAAIAiKKgAAAAUQUEFAACgCAoqAAAARVjcXOh2u2s3bdrUVhYAAAD2Id1ud21zefG4xz/b7XbnLg0AAADUTPEFAACgCAoqAAAARfj/lUqiwpreHXAAAAAASUVORK5CYII=" + } + }, + "cell_type": "markdown", + "id": "97a5079d", + "metadata": {}, + "source": [ + "
\n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "343e60c3", + "metadata": {}, + "source": [ + "However, to update the entries on the boundary of the locally stored vector we need entries stored on other processors." + ] + }, + { + "attachments": { + "f2.png": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6gAAADyCAYAAABAvOgkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAewgAAHsIBbtB1PgAAABl0RVh0U29mdHdhcmUAd3d3Lmlua3NjYXBlLm9yZ5vuPBoAABgySURBVHic7d19lJ11YSfw70wSAoRIQJsY3kzCm4b3ovJSlQWxtj12PYtVu8eK3ZblRU/pqSvH3W13QcWu2pd1faPJWnZPX3brtud0u9bWV4wCLYixgq8gQigvgWgiIQEyhOTuH88zzDOTO3fuTGbm+c3k8znnnvk99/7uvV8ydy7znef3PHdg3bp1nQAAAEDLBtsOAAAAAImCCgAAQCEWjtl+7cDAwKY2ggAAAHBg6XQ6q5J8bnh7VEEdGBjYdPnll98z26EAAAA48Kxfvz6dzshpkSzxBQAAoAgKKgAAAEVQUAEAACiCggoAAEARFFQAAACKoKACAABQBAUVAACAIiioAAAAFEFBBQAAoAgKKgAAAEVQUAEAACiCggoAAEARFFQAAACKoKACAABQhIVtBwAAinR2fdmc5NMtZwHgAKGgAgDdvC7JdUk2ZP8K6lFJXpqq7J6cavXWHUl+b//iwYw6N8nHkjyT5PyWMixN9bv6niRPtJQBZp2CCgDMhGOSfC3Jyi63LZrlLDBZz0v1R5WhWXzO16b6Y87wH3SOra//bpJTZjEHtEpBBQC6+dtUy3s3T/H+izNSTrcn2ZhkbZIX7n80mJc+k2RB2yGgbQoqANDNxvoyVVuT/HL9GD9M0km1XFhBZS74fJKBKdxvVZIrU/1h579N8r47k9yVkZ+9Nyb5l1PIAHOaggoAdHNEfdmV5JEp3P/xJJ+a1kRQvmOSvDvJnZl8QT0yyd7G9qumKxTMJT5mBgDo5upUez7/vO0g0IITk3wgyfWz+Jx7J54C85+CCgAAo61OtSf0XW0HgQONggoAAEARFFQAAACKoKACAABQBGfxBQAm4/Qky7tcP5Tk5lnOAm1YnOpzgrtZVn89PskXxpnz2SR/MN2hYL5QUAGAybguyb/qcv3DqT5iA+a7wSQXTzDnsB5z/nl648D8oqACAJNxV5KlXa7/8WwHgZY8k+RN49z24iTvTfJAkmvGmXP/TISC+UJBBQAm47q2A0DL9iT5y3Fue0X99fEec4AenCQJAACAIiioAAAAFMESXwBgprw9yXGN7TX111OSfKBx/d1J/sdshQKgXAoqADBT3prk3C7Xn5jk3Y3tv4+CCjckeXNj+9D668lJtjWu35Tkp2cpE8w6BRUA6OaWJB9Mct9+PMafJvlKH/Pu3o/ngPni0CRHdLl+wZjrfzI7caAdCioA0M2X6sv++MR0BIEW3Jzk+CSdSd7v6/X9npnCc16T5D19zNs9hceGOUNBBQCA0Z7O1FYP7Jri/ZJkS32BA5qz+AIAAFAEBRUAAIAiKKgAAAAUQUEFAACgCAoqAAAARVBQAQAAKIKCCgAAQBEUVAAAAIqgoAIAAFAEBRUAAIAiKKgAAAAUYWHbAaCbG2644ejBwcE3tZ2jlxtvvPFVSQaS5P77779ry5YtP2k50linJTmyHj+Y5L4Ws3RzTJLj6/HOJBtbzNKPG5LsajsEfXlTkqPbDtHDi5Ksqsfbk3yzvShdPS/JWY3tm5PsbSnLeF6ZkT+yfyvJthazTOQbSb7Sdgj68qIkl7QdooeBJK9qbH8z1XtISc5IsqweP5BkU3tR9rVq1aoXrVixYlWSLFu2bPsll1xS2vvvKJs3b/7odddd92zbOWabgkqRFixYcEKn0/nDtnP0cscdd2Tv3ud+Z3t9m1nmibe0HWACfxIFda64OsnPtB1iEi5tO8AESn9/Kz3fh6OgzhUvSVL07x5jlP7aL86mTZuyadOmJMmaNWuSwt9/V65cuS7JAVdQLfEFAACgCPagMlfc1XaALk4fHixevPjRoaGhLW2G6WJtRn7Gn07ygxazdLMq1VLCJOmkWqZXkoOTnNR2CPbbY/WlJMcnWVKP9yT5TotZunl+Ri2RPuy+FLfEd+cJjY0Sv8erkyxtOwT767B7204wxkCy8/jG9iNJftxWmHGckmRBPX4yyQ9bzLKPxYsXrxkaGjosSTqdzt4k32450liHJjlhwlnznILKXPHSK664YnfbIZoGBwc7w+Pzzjvv4xs2bLi+zTxd/CjJC+rxd5K8rMUs3Xw6yevq8a5Ux62U5KxUx44xt/1RkuvaDjHGhiQX1OMnU95r/x1JPjayectlyRlPtpamq4HbGxslfo8/k+QX2g7B/hjYm+wo7NCT7x2crG0uF/9gko+0lWYcjyc5vB7/U6rjxYtx6qmn3rRx48YLk2T37t1PXnHFFUW9/65fv/78Tqdza9s52maJLwAAAEVQUAEAACiCggoAAEARFFQAAACKoKACAABQBAUVAACAIiioAAAAFEFBBQAAoAgKKgAAAEVQUAEAACiCggoAAEARFFQAAACKoKACAABQBAUVAACAIiioAAAAFEFBBQAAoAgKKgAAAEVQUAEAACiCggoAAEARFNT5444kX68va/u8z1837vOzM5QLAACgLwvbDsC0OTvJQD1e0ud9Tk1yQj0+ctoTAQAATII9qAAAABRBQQUAAKAICioAAABFUFABAAAogoIKAABAERRUAAAAiqCgAgAAUAQFFQAAgCIoqAAAABRBQT2wLWo7AAAAwDAFdf4YaowX9zF/IMkLZigLAADApCmo88cTjfEL+5h/SpIlM5QFAABg0hTU+eMHjfF5fcz/tZkKAgAAMBUK6vxxe2P8liTP6zH35UneMbNxAAAAJmdh2wGYNn+W5J31eEWSTyX5lSRbG3MGkvzrJJ9IdYKkofR3vGrrBgYGXrZu3bo9bedouuqqq54bP/roo8ckOae9NF01f76XpLx8RzTGgykv30ltB2BaHJ3yXltLG+MFKS/fqtGbH31JsuTpVpL0p8Tv8bK2AzAdfvOUthOMtuugMVe8KOW99hc0xoelsHw7dux47v13YGBgwbp164rK1+l01radoQQK6vzxT6lK6Zvr7Z9Lcn+SDUk2pyoD56f6H3mSXJvkrUlOmNWUU9TpdG5tO0Mv3//+969IckXbOXp4SZLb2g7Rw+KUnY+567L6UqolKf61/8cfbzvBBEr/HjMndQaTj9zYdooJvDMjOydKdGYKe3+75557nhsvWrTo0BSWj4qCOr9cmeqvaefW20uT/OKYOZ0kH0ryvlQFFQAAoAgK6vzyeJILk1yd6q/JJzZu25Pks0k+nOSL9XX/kGRTPX50diICAAB0p6DOP7tS7SH9UJKfSnU86lNJHkryzJi5b5vdaP175JFHbl6zZs1hbefo5Q1veMOhzz777ECSfOtb39p17733FnWMbKpls8M/47uz7/e/bQszcgz0nlSv3ZI92XYA+nZxRh8HVZpFSYaPJSvxtT+Y5JDGdomv/UNTnVchqf79Snv/bdrddgD69oVUx02WrPkRgU8n2dtWkHEcnJH332dS2Ov/zDPPXLR69eqDkmTp0qV7DjnkkNLef0e59NJLn2o7QxsG1q1b13luY2Dg5Msvv/yeXncAAACA6bB+/fqTOp3O3cPbPmYGAACAIiioAAAAFEFBBQAAoAgKKgAAAEVQUAEAACiCggoAAEARFFQAAACKoKACAABQBAUVAACAIiioAAAAFEFBBQAAoAgKKgAAAEVQUAEAACiCggoAAEARFFQAAACKoKACAABQBAUVAACAIiioAAAAFEFBBQAAoAgKKgAAAEVQUAEAACiCggoAAEARFFQAAACKoKACAABQBAUVAACAIiioAAAAFEFBBQAAoAgKKgAAAEVQUAEAACiCggoAAEARFFQAAACKoKACAABQBAUVAACAIiioAAAAFEFBBQAAoAgKKgBwIDkzyU1Jvth2EAD2tbDtAAAAs+iIJBcm6bQdBIB92YMKAABAEexBBQAOJLcmOSr2oAIUSUEFAA4kzyTZ3HYIALpTUAGAA8nRSS5JtQf1Yy1nAWAMBRUAOJCclOQjUVABiuQkSQAAABRBQQUAAKAICioAAABFUFABAAAogoIKAABAERRUAAAAiqCgAgAAUAQFFQAAgCIoqAAAABRBQQUAAKAICioAAABFUFABAAAowsK2AwAAzKJvJHllkk7bQQDYl4IKABxItie5pe0QAHRniS8AAABFUFABAAAogoIKAABAERRUAAAAiqCgAgAAUAQFFQAAgCIoqAAAABRBQQUAAKAIC9sOAN2sX7/+zE6n84m2c/Ty/ve//5S9e/cOJMnWrVs3Pf300zvbzjTGqiSH1eMfJ3m0vShdvSDJC+vxriT3tpilHz+X5Im2Q9CXG5Kc0XaIHpbXlyR5Ksl9LWbpZkmS1Y3t7yTptJRlPKckGajHDyTZ0WKWifxFko+0HYK+nJPkv7YdooeBVK/9YfcnebKlLONZneo9JEl+lOSxFrPs4/DDD1++dOnS5UmycuXKpy677LLS3n9HGRoauvDqq68eajvHbFNQKdXhSc5rO0QvDz30UPbu3Tu8eVqbWfqwNKN/4SzN0iQ/1XaICSxqOwB9Oy2Fv380LE2you0QEzi37QATOLXtABO4ve0A9O2IzJ33jiQ5ve0AE1iaZE3bIZq2b9+e7du3J0kOOuig4t9/Fy9evKDtDG2wxBcAAIAi2IPKnDAwMHDt3sbuykK8b3iwfPnyL23ZsmVDi1m6+Q9JDq3HDyf5oxazdPMrSU6ux88meU+LWbpZmeTtbYdgv305yU1thxjj11MtwU+SoSTXtxelq3OSvK6xfX2qnCV538hw5deT5Rvbi9LNPT+fPH1c2ynYL3uTXNt2iDEWJfnPje3PJLmtpSzj+e0kB9fjB5J8ssUs+1ixYsW/eeyxx9YkyZ49e57pdDrvm+g+s2lgYODYJJe3naNtCipzQqfT+S9XXnnl7rZzNA0ODj73prZ27doNW7ZsKe2XzN/MSEHdnDJ/CR4uqLtTXr6zoqDOB19Nea+ti1N2QX1HRhfU30t5x183fqk86ZvJhhvbi9LN8tMU1Dmvk/J+Ng/N6IL6+ZR3fPO7MlJQH0xh/4bHHHPMRY2COnTllVcWlW/9+vXndzqdA76gWuILAABAERRUAAAAiqCgAgAAUAQFFQAAgCIoqAAAABRBQQUAAKAICioAAABFUFABAAAogoIKAABAERRUAAAAiqCgAgAAUISFbQcY411JFtTjP0vycB/3+fUkL6jHn0ny7RnINezEJJfU46eSfLQePy/JL9e3HZdkWZItSW5P8j+T/OMUnuuE+vEuSnJMkucn2Z7q3+QLSf5Xkn/ucf+fTvKaevxIkj/tMffFSV7f2P56ki/1mP/qJC+tx99N8ukecwEAAPpSWkG9PsnienxL+iuo70yyth4/mpktqGuTfKAeb0lVUF+R5H+nKpFNK5OckeTfJrkhyW8k2dvHcxye5ENJfi37fn9emOTkVKX12iS/X3/t9riHNLLuTPKpJM+M85xvS/LvG9tfTe+C+p+SXFCPr4mCCgAATANLfPfPK1MVueFyuiPJfan2rg4bSPL2VOV7IsekKuaXZ3Q5fTDJN+rH7tTXHZzkd1LtSe32ffxanSdJDktyTo/nffWY7XOTLBln7qH17cNu6vG4AAAAfVNQp25Jkr9KclCqknhmqqW+x9dffz7Jpsb8a5Ks6fF4Byf5f0lOrbeHkrwnydGplg2fXT/2sUnWNe735vqxx9qdak/osLEldNgRqZYDNx2Uqnx384qM7OXemuSb48wDAACYFAV16pYkWZ7kPyZ5S5I7G7ftSfLZJBdnZG/qwiS/2uPxrktyVj1+MsmF9XWPjJn3cJIrU+09bd53eZfHbC7THa+gXpCR435v7WN+8/ovp79lywAAABNSUPfPTRk5zrObHyb588b2q8aZtyzVMuBh12TiEyv9bqqTGSXV3tfLx8k37JxUS33HahbO9ybZ1eX68eb3Ok4VAABgUhTU/fPhjBwTOp5mSTxlnDm/lGRpPd6a5JN9PHcnyX9vbL+my5y7Up3MKUkWpXtBHi6cO5NsyMhe1DMycnbkYUdmZC9v4vhTAABgGimo++erE08ZdRzqslQnTRrrgsb4i6mOH+3H7Y3xy7rc3km1DHfYxWNuPzrJS+rxzanO8ju8V3Qw1TLjpgsz8pp5MMk9feYEAACYUGkfMzOX7Ej1uaQTebIxXpjqBEO7xsw5ozFekeTdfWY4vDE+JNVe2B1j5nwp1YmUkn2X7V7UGH+xMX/YxUn+srFteS8AADBjFNSpe7rPeWOXAHfbg/r8xvhf1JepWJbuBXXYaalOpjS87Ldb4dyY5PH6scYWWgUVAACYMZb4lqHbyYumotv3876MLDMeyOi9psOFc0uq41WT6gzEG+rx8UlW1eNjk5zUuK/jTwEAgGk1HwrqfPhvaC4V/u1UJyOayuXBcR5/7LLdJDk5yTH1+MsZvae328fTNPeefi/7fvwNAADAfimt3O1pjBeMO2u0ZTMRZJZtaYxXJvnJFC/jfSbpRIXzixltovmW9wIAANOutIL6RGN8RB/zj0x1UqG57rbGeOyZdqfDTRnZQ7oq1dLdXoXze0kerscXpXqdXNRjPgAAwH4rraBuboxP62P+69P9pENzzeca4xdn6idJGs9jSb7T2H5N4zl+mOT+LvcZPsZ0RZI3Jjmq3m4eowoAADBtSiuo32iM35je5XNpkt+Z2Tiz5u+SfL+x/dFM34mThjX3ev5Wqr3PY68fb/57G+Phs/wCAABMq9IK6v9tjE9P8q5x5h2Z5G+SrJnxRLNjT5JrMrIM99Qkn09yXB/3PSnJx1MV+l6ahfOkca7vZ76z9wIAADOitM9B/Wyqpain1NsfSvLKJH+V6qyxy5Kcn+RXUx2jekuSFyY5YbaDzoC/TbWn8tp6+7wk9yT5P0m+kGoZ7tOp9qweneSsVEt1z6jnf22Cx/9Kkmcz+nu+N+MXzoeS3J3qbL9Njj8FAABmRGkF9dkkl6YqTYfX1/1ifRnre0l+KfNrj951SX6c5A+TLEqyOMlb68v+eiLJHamK77A76+cbz5cyuqDuSnLrNGQBAADYR2lLfJPqONQLktw+zu3PJPnjJC9PdfKf+eZjqU4Q9SepCmEvT6Y6fvXSVHuZJzK2zE+0N3Ts7f+Qai8uAADAtCttD+qwO1Pt6Ts9yc+kWsb7VJIHUh2bubUx97yMfGbqkzOc6+8ycnKhTq+JDd9r3CeZuHQm1dLatyW5ItWS5hcneUGSQ1KdoOiRes43UxX2fr03yR80tp+aYP7fZHT2oUk8FwAAwKSUWlCTqgDeWV96eWKC26fT7iQ/meR99kzhPsN2pdrrOV3LmJ/J5Art/mQHAACYlBKX+AIAAHAAUlABAAAoQslLfKHpg+vWrdvbdoimq6666rnxXXfd9dpUH4NUkiWN8bFJfr+tIONY2xgvSnn5lrcdgGnxs6k+nqskxzfGB6e81/6ZY7avz+QOD5ll3315ctIhbacYbceqthOw3wZT3s/m2N/bX5/kuDaC9HBwY7w6hf0bbtq06bmPply4cOHidevWFZUvyVFtByjBfC6oSzM9/317MrvHudLdb7UdoJdt27a9Iskr2s7Rw4ok/67tED0sTNn5mLvOy+iP1yrNQSn/tf8bbQfo7UenVxeYVgMp/2fzovpSqqNT2L/h1q0j51kdHBws7v230+n3HKzz23wuqH+d5NXT8Dgbk7x0Gh4HAACAHuZzQWUO27Nnz9bBwcG/bztHL6tXrz670+kMJMmWLVt+sHPnzu1tZxrjxCSH1+NHkzzUYpZuVqRaepxUH3n03Raz9KPgJY6M8Y8pe+XLURlZxrUzyfdbzNLNYak+3mzYxvT/0Wqz5exUe7iS5AdJSnv/bSr9vY0RP0pS8u8eA6le+8PuTrKjpSzjOTnVKsYk2Zzk4Raz7OPII49cuWzZsqOTZPny5SW+/44yNDS0p+0MbZjPBfWTSb4wDY/z6DQ8BpN01VVXfTvJL7SdA5iTrmk7ADAnbYzfPea1bdu2Zdu2bUmS++67L7fddlvLiehmPhfUv2g7AAAAAP3zMTMAAAAUQUEFAACgCAoqAAAARVBQAQAAKIKCCgAAQBEUVAAAAIqgoAIAAFAEBRUAAIAiKKgAAAAUQUEFAACgCAoqAAAARVBQAQAAKIKCCgAAQBEUVAAAAIqgoAIAAFAEBRUAAIAiKKgAAAAUQUEFAACgCAoqAAAARVBQAQAAKIKCCgAAQBEUVAAAAIqgoAIAAFAEBRUAAIAiKKgAAAAUQUEFAACgCAoqAAAARVBQAQAAKIKCCgAAQBEUVAAAAIqwsLnR6XRWrV+/vq0sAAAAHEA6nc6q5vbCMbd/rtPpzF4aAAAAqFniCwAAQBEUVAAAAIrw/wHmUqtWcFnIlwAAAABJRU5ErkJggg==" + } + }, + "cell_type": "markdown", + "id": "fce954fe", + "metadata": {}, + "source": [ + "
\n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "3f90d701", + "metadata": {}, + "source": [ + "Thus, in order to update the local entries in `u_new`, we also need some remote entries of vector `u` located in neighboring processes. Figure below shows the entries of `u` needed to update the local entries of `u_new` in a particular process (CPU 2)." ] }, { @@ -404,10 +457,9 @@ "id": "8ed4129c", "metadata": {}, "source": [ - "### Code\n", + "### MPI Code\n", "\n", - "\n", - "Take a look at the implementation below and try to understand it. Note that we have used MPIClustermanagers and Distributed just to run the MPI code on the notebook. When running it on a cluster, MPIClustermanagers and Distributed are not needed.\n" + "Take a look at the implementation below and try to understand it.\n" ] }, { @@ -417,7 +469,7 @@ "metadata": {}, "outputs": [], "source": [ - "] add MPI MPIClusterManagers" + "] add MPI" ] }, { @@ -427,30 +479,7 @@ "metadata": {}, "outputs": [], "source": [ - "using MPIClusterManagers \n", - "using Distributed" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "e0d63c6b", - "metadata": {}, - "outputs": [], - "source": [ - "if procs() == workers()\n", - " nw = 3\n", - " manager = MPIWorkerManager(nw)\n", - " addprocs(manager)\n", - "end" - ] - }, - { - "cell_type": "markdown", - "id": "d7fb9177", - "metadata": {}, - "source": [ - "First, we implement the function to be called on the MPI ranks." + "using MPI" ] }, { @@ -462,120 +491,99 @@ }, "outputs": [], "source": [ - "@everywhere workers() begin\n", + "code = quote\n", " using MPI\n", - " comm = MPI.Comm_dup(MPI.COMM_WORLD)\n", + " MPI.Init()\n", " function jacobi_mpi(n,niters)\n", + " comm = MPI.COMM_WORLD\n", " nranks = MPI.Comm_size(comm)\n", " rank = MPI.Comm_rank(comm)\n", " if mod(n,nranks) != 0\n", " println(\"n must be a multiple of nranks\")\n", " MPI.Abort(comm,1)\n", " end\n", - " n_own = div(n,nranks)\n", - " u = zeros(n_own+2)\n", + " load = div(n,nranks)\n", + " u = zeros(load+2)\n", " u[1] = -1\n", " u[end] = 1\n", " u_new = copy(u)\n", " for t in 1:niters\n", - " reqs = MPI.Request[]\n", + " # Communication\n", " if rank != 0\n", " neig_rank = rank-1\n", - " req = MPI.Isend(view(u,2:2),comm,dest=neig_rank,tag=0)\n", - " push!(reqs,req)\n", - " req = MPI.Irecv!(view(u,1:1),comm,source=neig_rank,tag=0)\n", - " push!(reqs,req)\n", + " s = 2\n", + " r = 1\n", + " MPI.Sendrecv!(view(u,s:s),view(u,r:r),comm;dest=neig_rank,source=neig_rank)\n", " end\n", " if rank != (nranks-1)\n", " neig_rank = rank+1\n", - " s = n_own+1\n", - " r = n_own+2\n", - " req = MPI.Isend(view(u,s:s),comm,dest=neig_rank,tag=0)\n", - " push!(reqs,req)\n", - " req = MPI.Irecv!(view(u,r:r),comm,source=neig_rank,tag=0)\n", - " push!(reqs,req)\n", + " s = load+1\n", + " r = load+2\n", + " MPI.Sendrecv!(view(u,s:s),view(u,r:r),comm;dest=neig_rank,source=neig_rank)\n", " end\n", - " MPI.Waitall(reqs)\n", - " for i in 2:(n_own+1)\n", + " # Local computation\n", + " for i in 2:(load+1)\n", " u_new[i] = 0.5*(u[i-1]+u[i+1])\n", " end\n", " u, u_new = u_new, u\n", " end\n", - " return u\n", - " end\n", - "end" - ] - }, - { - "cell_type": "markdown", - "id": "6eab32d0", - "metadata": {}, - "source": [ - "In order to check the result, we will compare it against the serial implementation. To this end, we need to define the serial implementation also in the workers." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "f1a8db8f", - "metadata": {}, - "outputs": [], - "source": [ - "@everywhere workers() function jacobi(n,niters)\n", - " u = zeros(n+2)\n", - " u[1] = -1\n", - " u[end] = 1\n", - " u_new = copy(u)\n", - " for t in 1:niters\n", - " for i in 2:(n+1)\n", - " u_new[i] = 0.5*(u[i-1]+u[i+1])\n", + " # Gather the results\n", + " if rank !=0\n", + " lb = 2\n", + " ub = load+1\n", + " MPI.Send(view(u,lb:ub),comm,dest=0)\n", + " u_all = zeros(0) # This will nevel be used\n", + " else\n", + " u_all = zeros(n+2)\n", + " # Set boundary\n", + " u_all[1] = -1\n", + " u_all[end] = 1\n", + " # Set data for rank 0\n", + " lb = 2\n", + " ub = load+1\n", + " u_all[lb:ub] = view(u,lb:ub)\n", + " # Set data for other ranks\n", + " for other_rank in 1:(nranks-1)\n", + " lb += load\n", + " ub += load\n", + " MPI.Recv!(view(u_all,lb:ub),comm;source=other_rank)\n", + " end\n", " end\n", - " u, u_new = u_new, u\n", + " return u_all\n", " end\n", - " u\n", - "end" - ] - }, - { - "cell_type": "markdown", - "id": "d2b04c67", - "metadata": {}, - "source": [ - "Finally, we call the parallel function on the workers, gather the results on the root rank, and compare against the sequential solution." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "68851107", - "metadata": { - "code_folding": [] - }, - "outputs": [], - "source": [ - "@everywhere workers() begin\n", - " # Call jacobi in parallel\n", - " niters = 10\n", - " load = 4\n", - " nranks = MPI.Comm_size(comm)\n", - " n = load*nranks\n", - " u = jacobi_mpi(n,niters)\n", - " # Gather results in root process and check\n", - " rank = MPI.Comm_rank(comm)\n", - " n_own = div(n,nranks)\n", - " if rank == 0\n", - " results = zeros(n+2)\n", - " results[1] = -1\n", - " results[n+2] = 1\n", - " rcv = view(results, 2:n+1)\n", - " else\n", - " rcv = nothing\n", + " function jacobi(n,niters)\n", + " u = zeros(n+2)\n", + " u[1] = -1\n", + " u[end] = 1\n", + " u_new = copy(u)\n", + " for t in 1:niters\n", + " for i in 2:(n+1)\n", + " u_new[i] = 0.5*(u[i-1]+u[i+1])\n", + " end\n", + " u, u_new = u_new, u\n", + " end\n", + " u\n", " end\n", - " MPI.Gather!(view(u,2:n_own+1),rcv,comm;root=0)\n", - " if rank == 0\n", - " @show results ≈ jacobi(n,niters)\n", - " end \n", - "end" + " function testit(load)\n", + " comm = MPI.COMM_WORLD\n", + " nranks = MPI.Comm_size(comm)\n", + " rank = MPI.Comm_rank(comm)\n", + " n = load*nranks\n", + " niters = 100\n", + " u_par = jacobi_mpi(n,niters)\n", + " if rank == 0\n", + " # Compare agains serial\n", + " u_seq = jacobi(n,niters)\n", + " if u_par ≈ u_seq\n", + " println(\"Test passed 🥳\")\n", + " else\n", + " println(\"Test failed\")\n", + " end\n", + " end\n", + " end\n", + " testit(3)\n", + "end\n", + "run(`$(mpiexec()) -np 4 julia --project=. -e $code`);" ] }, { @@ -584,7 +592,7 @@ "metadata": {}, "source": [ "
\n", - "Question: In function jacobi_mpi, how many messages per iteration are sent from a process away from the boundary?\n", + "Question: In function jacobi_mpi, how many messages per iteration are sent from a process away from the boundary (excluding the messages to collect the data on rank 0)?\n", "
\n", "\n", " a) 1\n", @@ -605,32 +613,6 @@ "jacobi_2_check(answer)" ] }, - { - "cell_type": "markdown", - "id": "075dd6d8", - "metadata": {}, - "source": [ - "
\n", - "Question: At the end of function jacobi_mpi ...\n", - "
\n", - "\n", - " a) each process holds the complete solution.\n", - " b) the complete solution is gathered in the root process. \n", - " c) each process contains the solution for the local partition. \n", - " d) the ghost cells of u contain the initial values -1 and 1 in all processes." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "c3b58002", - "metadata": {}, - "outputs": [], - "source": [ - "answer = \"x\" # replace x with a, b, c or d\n", - "jacobi_3_check(answer)" - ] - }, { "cell_type": "markdown", "id": "c9aa2901", @@ -842,7 +824,7 @@ "- We update $N^2/P$ items per iteration\n", "- We need data from 2 neighbors (2 messages per iteration)\n", "- We communicate $N$ items per message\n", - "- Communication/computation ratio is $O(P/N)$" + "- Communication/computation ratio is $2N/(N^2/P) = 2P/N =O(P/N)$" ] }, { @@ -876,7 +858,7 @@ "- We update $N^2/P$ items per iteration\n", "- We need data from 4 neighbors (4 messages per iteration)\n", "- We communicate $N/\\sqrt{P}$ items per message\n", - "- Communication/computation ratio is $O(\\sqrt{P}/N)$" + "- Communication/computation ratio is $ (4N/\\sqrt{P})/(N^2/P)= 4\\sqrt{P}/N =O(\\sqrt{P}/N)$" ] }, { @@ -939,7 +921,7 @@ "- Both 1d and 2d block partitions are potentially scalable if $P<