Merge pull request #36 from fverdugo/francesc

Enhancements in Jacobi notebook
2026-11-15 06:17:17 +01:00 · 2024-08-20 15:30:14 +02:00
commit 4f0746e0a4
@@ -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",
@@ -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
 ```

@@ -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": [
    "<div class=\"alert alert-block alert-info\">\n",
-    "<b>Note:</b> 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",
+    "<b>Note:</b> 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",
    "</div>"
   ]
  },
@@ -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/9pw4YNRX1AGhoaemandtxxx23etm1baR8yLine truncated
+    }
+   },
+   "cell_type": "markdown",
+   "id": "97a5079d",
+   "metadata": {},
+   "source": [
+    "<div>\n",
+    "<img src=\"attachment:f1.png\"  align=\"left\" width=\"450\"/>\n",
+    "</div>"
+   ]
+  },
+  {
+   "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/MLine truncated
+    }
+   },
+   "cell_type": "markdown",
+   "id": "fce954fe",
+   "metadata": {},
+   "source": [
+    "<div>\n",
+    "<img src=\"attachment:f2.png\"  align=\"left\" width=\"450\"/>\n",
+    "</div>"
+   ]
+  },
+  {
+   "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": [
    "<div class=\"alert alert-block alert-success\">\n",
-    "<b>Question:</b>  In function jacobi_mpi, how many messages per iteration are sent from a process away from the boundary?\n",
+    "<b>Question:</b>  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",
    "</div>\n",
    "\n",
    "    a) 1\n",
@@ -605,32 +613,6 @@
    "jacobi_2_check(answer)"
   ]
  },
-  {
-   "cell_type": "markdown",
-   "id": "075dd6d8",
-   "metadata": {},
-   "source": [
-    "<div class=\"alert alert-block alert-success\">\n",
-    "<b>Question:</b>  At the end of function jacobi_mpi ...\n",
-    "</div>\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<<N$\n",
    "- The 2d block partition has the lowest communication complexity\n",
    "- The 1d block partition requires to send less messages (It can be useful if the fixed cost of sending a message is high)\n",
-    "- The best strategy for a given problem size will thus depend on the machine, but the 2d block partition is the one used in practice since it has the lowest communication complexity.\n",
+    "- The best strategy for a given problem size will thus depend on the machine.\n",
    "- Cyclic partitions are impractical for this application (but they are useful in others) \n",
    "\n"
   ]
@@ -977,15 +959,15 @@
 ],
 "metadata": {
  "kernelspec": {
-   "display_name": "Julia 1.9.1",
+   "display_name": "Julia 1.10.0",
   "language": "julia",
-   "name": "julia-1.9"
+   "name": "julia-1.10"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
-   "version": "1.9.1"
+   "version": "1.10.0"
  }
 },
 "nbformat": 4,