seb/XM_40017
1
0 Fork 0
mirror of https://github.com/fverdugo/XM_40017.git synced 2026-11-15 06:57:17 +01:00
Code Issues Actions Packages Projects Releases Wiki Activity

Work in asp notebook

Browse Source
This commit is contained in:
Francesc Verdugo committed 2024-09-18 12:01:27 +02:00
1 parent e146513dfe
commit ab95b334eb
2 files changed
+1509 -1476
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

notebooks/asp.ipynb
+297 -268
View File
@@ -38,18 +38,10 @@
},
{
"cell_type": "code",
"execution_count": 17,
"execution_count": null,
"id": "1dc78750",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"🥳 Well done!\n"
]
}
],
"outputs": [],
"source": [
"using Printf\n",
"\n",
@@ -65,11 +57,19 @@
"function q1_answer(bool)\n",
" bool || return\n",
" msg = \"\"\"\n",
" It is not needed to start the loop over `j` with `j=k` as none of the entries of the matrix will be updated\n",
" in this particular case. Rememebr: `C[i,j] = min(C[i,j],C[i,k]+C[k,j])` if we substitute `j=k`, we get:\n",
" `C[i,k] = min(C[i,k],C[i,k]+C[k,k])`. Rememeber that `C[k,k]=0, thus, `C[i,k] = min(C[i,k],C[i,k])`, or\n",
" `C[i,k] = C[i,k]`. In other words, the new value of `C[i,k]` will correspond to the old value.\n",
" The same is true for `i=k`.\n",
" The we can change the loop order over i and j without changing the result. Rememebr:\n",
" \n",
" C[i,j] = min(C[i,j],C[i,k]+C[k,j])\n",
" \n",
" if we substitute j=k, we get\n",
" \n",
" C[i,k] = min(C[i,k],C[i,k]+C[k,k]).\n",
" \n",
" Since C[k,k]=0, thus, C[i,k] = min(C[i,k],C[i,k]), and C[i,k] = C[i,k].\n",
" \n",
" In other words, the value of C[i,k] will not be updated at iteration k.\n",
" \n",
" The same is true for i=k.\n",
" \"\"\"\n",
" println(msg)\n",
"end\n",
@@ -97,8 +97,8 @@
"id": "1faddbfa",
"metadata": {},
"source": [
"We represent the distance table as a matrix $C$, where $C_{ij}$ is the distance from node $i$ to node $j$ via a direct connection (a single hop in the graph). If there is no direct connection from $i$ to $j$, this is represented using a very large value in $C_{ij}$ (e.g. the largest possible floating point number, `inf`). \n",
"The next figure shows a simple directed graph with 4 nodes an its corresponding distance matrix (matrix labeled as \"input\")."
"We represent the distance table as a matrix $C$, where $C_{ij}$ is the distance from node $i$ to node $j$ via a direct connection (a single hop in the graph). If there is no direct connection from $i$ to $j$, this is represented using a large value in $C_{ij}$ representing infinity. \n",
"The next figure shows a simple directed graph with 4 nodes an its corresponding distance matrix (labeled as \"input\")."
]
},
{
@@ -121,7 +121,7 @@
"id": "ade31d26",
"metadata": {},
"source": [
"The ASP problem consists in computing the minimum distance between any two pair of nodes $i$ and $j$ in the graph. All these values can be also represented as a matrix (matrix labeled as \"output\" in the figure above). For instance, the minimum distance from node 2 to node 3 is 8 as highlighted in the figure."
"The ASP problem consists in computing the minimum distance between any two pair of nodes $i$ and $j$ in the graph. All these values can be also represented as a matrix (labeled as \"output\" in the figure above). For instance, the minimum distance from node 2 to node 3 is 8 as highlighted in the figure. You can understand both input and output matrices as distance tables. The key difference is that the input contains the distance using direct connections only, whereas the output contains the (minimum) distance allowing indirect connections."
]
},
{
@@ -136,21 +136,10 @@
},
{
"cell_type": "code",
"execution_count": 4,
"execution_count": null,
"id": "4fe447c5",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"floyd! (generic function with 1 method)"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"outputs": [],
"source": [
"function floyd!(C)\n",
" n = size(C,1)\n",
@@ -176,25 +165,10 @@
},
{
"cell_type": "code",
"execution_count": 6,
"execution_count": null,
"id": "860e537c",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4×4 Matrix{Int64}:\n",
" 0 9 6 1\n",
" 2 0 8 3\n",
" 5 3 0 6\n",
" 10 8 5 0"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"outputs": [],
"source": [
"inf = 1000\n",
"C = [\n",
@@ -228,7 +202,7 @@
"source": [
"### The algorithm explained\n",
"\n",
"The main idea of the algorithm is to perform as many iterations as nodes in the graph. At iteration $k$, we update the distance matrix $C$ by finding the shortest paths between each pair of nodes, allowing indirect paths via nodes from 1 to $k$. At the last iteration, it is allowed to visit all nodes and, thus, the distance table will contain the minimum possible distances, i.e. the solution of the ASP problem.\n",
"The main idea of the Floyd–Warshall algorithm is to perform as many iterations as nodes in the graph. At iteration $k$, we update the distance matrix $C$ by finding the shortest paths between each pair of nodes, allowing indirect paths via nodes from 1 to $k$. At the last iteration, it is allowed to visit all nodes and, thus, the distance table will contain the minimum possible distances, i.e. the solution of the ASP problem.\n",
"\n",
"This process is cleverly done with three nested loops:\n",
"\n",
@@ -243,7 +217,7 @@
"end\n",
"```\n",
"\n",
"At each outer iteration $k$, we do a loop over the distance matrix $C$. For each pair of nodes $i$ and $j$ we compare the current distance $C_{ij}$ against the distance via node $k$, namely $C_{ik}+C_{kj}$, and update $C_{ij}$ with the minimum."
"At each outer iteration $k$, we do a loop over the distance matrix $C$. For each pair of nodes $i$ and $j$ we compare the current distance $C_{ij}$ against the distance via node $k$, namely $C_{ik}+C_{kj}$, and update $C_{ij}$ with the minimum. I.e., at iteration $k$ one checks if it is beneficial to visit node $k$ to reduce the distance between nodes $i$ and $j$."
]
},
{
@@ -266,7 +240,7 @@
"id": "722e330c",
"metadata": {},
"source": [
"The update of the distance matrix at each iteration is illustrated in the next figure for a small ASP presented above. We highlight in green the distances that are updated in each iteration."
"The update of the distance matrix at each iteration is illustrated in the next figure for the small ASP presented above. We highlight in green the distances that are updated in each iteration. Note that some distances that were initially infinity (i.e., no connection) will be updated with a finite value in this process. You can understand this as adding new edges in the graph as illustrated in the figures below."
]
},
{
@@ -352,28 +326,17 @@
"source": [
"### Serial performance\n",
"\n",
"This algorithm is memory bound, meaning that the main cost is in getting and setting data from the input matrix `C`. In this situation, the order in which we traverse the entries of matrix `C` has a significant performance impact.\n",
"Before starting to parallelize this code, we want to make sure that we have an efficient sequential implementation. In this algorithm, the order in which we traverse the entries of matrix `C` has a significant performance impact.\n",
"\n",
"The following function computes the same result as for the previous function `floyd!`, but the nesting of loops over i and j is changed.\n"
]
},
{
"cell_type": "code",
"execution_count": 8,
"execution_count": null,
"id": "75cac17e",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"floyd2! (generic function with 1 method)"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"outputs": [],
"source": [
"function floyd2!(C)\n",
" n = size(C,1)\n",
@@ -399,19 +362,10 @@
},
{
"cell_type": "code",
"execution_count": 10,
"execution_count": null,
"id": "907bc8c9",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 1.777208 seconds\n",
" 2.967238 seconds\n"
]
}
],
"outputs": [],
"source": [
"n = 1000\n",
"C = rand(n,n)\n",
@@ -457,32 +411,19 @@
"\n",
"\n",
"<div class=\"alert alert-block alert-success\">\n",
"<b>Question (hard):</b> Can we really parallelize the loops over `i` and `j` ? To compute `C[i,j]` at iteration `k`, we first need to compute `C[i,k]` and `C[k,j]`. In order words, it seems that the order of the loops over `i` and `j` cannot be arbitrary. The first value of `i` and `j` in the loops should be `k`. However, this is not really necessary, why?\n",
"<b>Question (hard):</b> Can we really parallelize the loops over `i` and `j` ? To compute `C[i,j]` at iteration `k`, we first need to compute `C[i,k]` and `C[k,j]`. In order words, it seems that the order of the loops over `i` and `j` cannot be arbitrary. However, this is not really true, why?\n",
"</div>"
]
},
{
"cell_type": "code",
"execution_count": 18,
"execution_count": null,
"id": "8d05b686",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"It is not needed to start the loop over `j` with `j=k` as none of the entries of the matrix will be updated\n",
"in this particular case. Rememebr: `C[i,j] = min(C[i,j],C[i,k]+C[k,j])` if we substitute `j=k`, we get:\n",
"`C[i,k] = min(C[i,k],C[i,k]+C[k,k])`. Rememeber that `C[k,k]=0, thus, `C[i,k] = min(C[i,k],C[i,k])`, or\n",
"`C[i,k] = C[i,k]`. In other words, the new value of `C[i,k]` will correspond to the old value.\n",
"The same is true for `i=k`.\n",
"\n"
]
}
],
"outputs": [],
"source": [
"uncover = true\n",
"q1_answer(uncover)"
"uncover = false\n",
"q1_answer(true)"
]
},
{
@@ -616,7 +557,7 @@
"metadata": {},
"source": [
"<div class=\"alert alert-block alert-success\">\n",
"<b>Question:</b> How much data is communicated in each iteration in this parallel algorithm?\n",
"<b>Question:</b> How much data is send from the owner of row k in each iteration in this parallel algorithm?\n",
"</div>\n",
"\n",
" a) O(N²/P)\n",
@@ -627,18 +568,10 @@
},
{
"cell_type": "code",
"execution_count": 25,
"execution_count": null,
"id": "1bf4de56",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"It's not correct. Keep trying! 💪\n"
]
}
],
"outputs": [],
"source": [
"answer = \"x\" # replace x with a, b, c or d\n",
"floyd_check(answer)"
@@ -659,7 +592,7 @@
"source": [
"### Computation complexity\n",
"\n",
"Each process updates $N^2/P$ entries per iteration. The computation complexity is $O(N^2/P)$."
"Each process updates $N^2/P$ entries per iteration. The computation complexity per iteration is $O(N^2/P)$."
]
},
{
@@ -699,7 +632,7 @@
"metadata": {},
"source": [
"<div>\n",
"<img src=\"attachment:fig-asp-efficiency-comm-2.png\" align=\"left\" width=\"300\"/>\n",
"<img src=\"attachment:fig-asp-efficiency-comm-2.png\" align=\"left\" width=\"250\"/>\n",
"</div>"
]
},
@@ -714,7 +647,7 @@
"- On the receive side $O(N)/O(N^2/P) = O(P/N)$\n",
"\n",
"\n",
"In summary, the send/computation ratio is $O(P^2/N)$ and the receive/computation ratio is $O(P/N)$. The algorithm is potentially scalable if $P^2<<N$."
"In summary, the send/computation ratio is $O(P^2/N)$ and the receive/computation ratio is $O(P/N)$. The algorithm is potentially scalable if $P^2<<N$. Note that this is worse than for matrix-matrix multiplication, which is scalable for $P<<N$."
]
},
{
@@ -749,159 +682,288 @@
},
{
"cell_type": "markdown",
"id": "680e56cf",
"id": "bfede7cc",
"metadata": {},
"source": [
"### Code\n",
"We are going to implement the method in a function with the following signature\n",
"\n",
"We split the code into two functions. The first function is called on the main process (the process running this notebook). It splits the input matrix into blocks of rows. Then, we use a remotecall to compute Floyd's algorithm in each worker with its corresponding block of rows.\n"
"```julia\n",
" floyd_mpi!(C,comm)\n",
"```\n",
"\n",
"The signature is similar to the sequential function `floyd!(C)`, but there are some important differences\n",
"\n",
"- This function is going to be called on several MPI ranks, but only one rank (rank 0 in the implementation below) will receive the input matrix, the other ranks will receive an empty matrix. It would not be memory efficient if all ranks receive a copy of the input matrix as this one could be very large.\n",
"\n",
"- On output, only the value of C will be correct on rank 0. It will not be efficient to try to write the result on all ranks.\n",
"\n",
"- The parallel function takes an MPI communicator object. This allows the user to decide which communicator to use. For instance, `MPI.COMM_WORLD` directly or a duplicate of this one (which is the recommended approach).\n",
"\n",
"The implementation of the parallel function is given in the following code snipped:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "22234b2f",
"id": "544466a6",
"metadata": {},
"outputs": [],
"source": [
"function floyd_dist!(C)\n",
" m = size(C,1)\n",
" @assert mod(m,nworkers()) == 0\n",
" nrows_w = div(m,nworkers())\n",
" @sync for (i,w) in enumerate(workers())\n",
" rows_w = (1:nrows_w) .+ (i-1)*nrows_w\n",
" Cw = C[rows_w,:]\n",
" @async C[rows_w,:] = remotecall_fetch(floyd_worker!,w,Cw,rows_w)\n",
"code1 = quote\n",
" function floyd_mpi!(C,comm)\n",
" myC = distribute_input(C,comm)\n",
" floyd_iterations!(myC,comm)\n",
" collect_result!(C,myC,comm)\n",
" end\n",
" C\n",
"end"
"end;"
]
},
{
"cell_type": "markdown",
"id": "9fc3be11",
"id": "ce656413",
"metadata": {},
"source": [
"The second function is the one that runs on the workers. Note that we use MPI for communication in this case."
"In the following cells, we discuss the implementation of the helper functions `distribute_input`, `floyd_iterations!`, and `collect_result!`."
]
},
{
"cell_type": "markdown",
"id": "7beadf8e",
"metadata": {},
"source": [
"### Distributing the input matrix\n",
"\n",
"Since only rank 0 receives the input matrix C, this rank needs to split it row-wise and\n",
"send the pieces to all other ranks. This is done in the function below. We start by communicating the problem size N to all ranks (at the start only rank 0 knows the problem size). This is trivially done with an `MPI.Bcast!`. Once all ranks know the problem size, they can allocate space for their local part of `C` , called `myC` in the code. After this, rank 0 sends the pieces to all other ranks. We do it here with `MPI.Send` and `MPI.Recv!`. This can also be done with `MPI.Scatter!`, but it is more challenging since we are using a row partition and Julia stores the matrices in column major order. Note that this algorithm can also be implemented using a column-wise matrix partition. In this case, using `MPI.Scatter!` would provably be the best option.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d5a8d3a8",
"id": "d7ca1516",
"metadata": {},
"outputs": [],
"source": [
"@everywhere function floyd_worker!(Cw,rows_w)\n",
" comm = MPI.Comm_dup(MPI.COMM_WORLD)\n",
" rank = MPI.Comm_rank(comm)\n",
" nranks = MPI.Comm_size(comm)\n",
" m,n = size(Cw)\n",
" C_k = similar(Cw,n)\n",
" for k in 1:n\n",
" if k in rows_w\n",
" # Send row k to other workers if I have it\n",
" myk = (k-first(rows_w))+1\n",
" C_k .= view(Cw,myk,:)\n",
" for proc in 0:(nranks-1)\n",
" if rank == proc\n",
" continue\n",
" end\n",
" MPI.Send(C_k,comm;dest=proc,tag=0)\n",
"code2 = quote\n",
" function distribute_input(C,comm)\n",
" rank = MPI.Comm_rank(comm)\n",
" P = MPI.Comm_size(comm)\n",
" # Communicate problem size\n",
" if rank == 0\n",
" N = size(C,1)\n",
" if mod(N,P) !=0\n",
" println(\"N not multplie of P\")\n",
" MPI.Abort(comm,-1)\n",
" end\n",
" Nref = Ref(N)\n",
" else\n",
" Nref = Ref(0)\n",
" end\n",
" MPI.Bcast!(Nref,comm;root=0)\n",
" N = Nref[]\n",
" # Distribute C row-wise\n",
" L = div(N,P)\n",
" myC = similar(C,L,N)\n",
" if rank == 0\n",
" lb = L*rank+1\n",
" ub = L*(rank+1)\n",
" myC[:,:] = view(C,lb:ub,:)\n",
" for dest in 1:(P-1)\n",
" lb = L*dest+1\n",
" ub = L*(dest+1)\n",
" MPI.Send(view(C,lb:ub,:),comm;dest)\n",
" end\n",
" else\n",
" # Wait until row k is received\n",
" MPI.Recv!(C_k,comm,source=MPI.ANY_SOURCE,tag=0)\n",
" end\n",
" for j in 1:n\n",
" for i in 1:m\n",
" @inbounds Cw[i,j] = min(Cw[i,j],Cw[i,k]+C_k[j])\n",
" end\n",
" source = 0\n",
" MPI.Recv!(myC,comm;source)\n",
" end\n",
" return myC\n",
" end\n",
" Cw\n",
"end"
"end;"
]
},
{
"cell_type": "markdown",
"id": "bc81a83d",
"id": "661b5b3e",
"metadata": {},
"source": [
"<div class=\"alert alert-block alert-success\">\n",
"<b>Question:</b> Which of the following statements is true?\n",
"### Running Floyd's updates in parallel\n",
"\n",
"As discussed above, we need to communicate row $k$ of matrix $C$ in order to perform iteration $k$ of the algorithm. The function below is similar to the sequential function `floyd!`, but there is a key difference. At the start of iteration $k$, the owner of row $k$ sends it to the other processors. Once row $k$ is available, all ranks can do the Floyd update locally on their portion of matrix $C$."
]
},
{
"attachments": {
"fig-asp-efficiency-comm-2.png": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA1AAAAMbCAYAAABDhHSeAAAACXBIWXMAAB7CAAAewgFu0HU+AAAAGXRFWHRTb2Z0d2FyZQB3d3cuaW5rc2NhcGUub3Jnm+48GgAAIABJREFUeJzs3Xt81dWV9/HPSggJIJcqNxEUW8TitYpYFFCsYlWIM9MiUYeZWqVFGwUqGrk8WquDIIpVaFrQUeplShEfnWnUto/ai0FEW7x3ULxFCYpEBCn3QNbzxz4hh3BOrifnknzfrxevnvM7+7d/C6EkK3vvtczdERERERGRzGRmxwP/BIwC+gPdAQO+AN4ESoEH3X1dA+aaD+TWM2wHsC0y/9vAq+7+ST3zTgaOibxd6u5/rC+WqHv7AD+JunSlJziJMbPpwG2Rt3nuvivuWCVQIiIiIiKZx8y+CswBxhISprrsAe4HbnD3L+uYcyvQqQnh/BUocvc/x5n3D8C5kbeT3X1+Qyc2s2OBt6IuZbt7VRNijDd/NvC/wMDIpToTqKxEPVhERERERJLDzM4EVgEXUZM8OWFF6C/As8A7QHWi0Q6YCLxqZke1QEhDgGfMbFwLzN3S5lKTPNWrXQsGIiIiIiIiCWZmI4HfU7PVbiswG/hV7a10ZtYD+AEwg7CydCTwFzP7hrtvqOdRPwLejfNZe+AwYARwMZBDyC0eMLMX3X1tY39fyWZmQ4AbgfzG3KcESkREREQkQ5hZd+ARapKnD4Bz3P3DWOPdvQK4zcz+m7AqdWjkVzFh9aouL7v7qnrG3GdmDwJPE5KqTsD1wKQG/HaSzsxmA6cBXwd6NWUObeETEREREckctxFWfgC+BL4VL3mK5u7/C3w/6tJ3zOxriQjI3Z8DHoi69N1EzNtCzoz8ik6enm/MBEqgREREREQygJn1BP4t6tIN7v5RQ+939z8AqyNvs2jk1rV6/N+o133M7JAEzp1IfwUeBxYCUwiVAb/dmAm0hU9EREREJDNcCuRFXm8EHmzCHAuB4ZHXOxIRVETtEumdCDGmFXefXPuameXFGhuPEigRERERkcwwIur1/7j7zsZOECkf3uAS4o3Qo9b7tEueEkVb+EREREREMsPpUa9fSlkUsV0Y9foDd9+WskhamBIoEREREZE0F2n2Gl344J1UxVJbpKx6YdSl/05RKEmhLXwiIiIiIunvK9Q0zIVQga+ldTazr8S43gnoAgwC/olwNis78tkm4PYkxJYySqBERERERNJfx1rvE1kAIp4/NXL8VuDCBjTozWjawiciIiIikv4213rfJSVRxLYXeAI4wd2XpzqYlqYVKBERERGR9PcPYA81378no8/Sy8CWGNerCFsI1wGvAc+6e+0y5rXtiXqdHXdUbNE5i0d+pYwSKBERERGRGMysH3ADcCqwEri9AYlCi3B3N7N3CeeOAE4Eft/Cj/2Ru69K0Fxbo153auS9B0W93unuSqBERERERNKFmR0CXA9MAjpELg8BJprZr4Cb3P2zFIS2nJoEalhTJjCzY4A+kbefuvvfExFYA0QnUI1dPYsen4ziGXXSGSgREREREcDMupjZT4EPCStPHWoNaQ/8EFhjZjeZWeckh/hc1OvzzOzQJsxxN/BM5NelCYmqYT6Ken1MI++NHl+egFiaRQmUiIiIiLRpZtbBzKYC7wM3AfUlRl2AnwLvm9kUM8tt6RgjHiecOwLIAW5szM1m1hUYHnXpbwmKqyHeino91MxqJ6d1OSvOPCmhBEpERERE2iQzyzGzHwJrgDuB7rXHdOjQgTFjxpCXlxdrih7AzwgrUleYWYsej3H3ykic1a40swsbMcVN1KyqbQJ+l6jYGuBlaoo/dAEua8hNZnYscHbUpb8kNqzGUwIlIiIiIm2KBRcRVjMWAX1rj8nJyWH8+PGsWLGCe++9l5dffpnCwkLat28fa8rDgf8E3jKzfzezlvweez41/ZkMWGZm3zczi3eDmWWZWRHw46jLt7v7zhaMcz/uXg78IerSHWZ2QV33mNlRwDJqqvZtBv5vy0TYcJbiIhYiIiIiIkljZucAtwMnx/o8KyuL0aNHM336dPr373/A5+Xl5cyfP58lS5awd+/eeI95C7jF3ZclKOz9mFkf4I/A0VGXXwR+BTxP2OZnwGGELXsTgcFRY/8IfNvdo0uLV8+9lZoqeacksAofZvZ1wkpU9RZJJ5zrehx4k1BoIhcYAIwCLo68r/ZDd78vUfFExZXH/o2J89x9V9zxSqBEREREpLUzs9OA24CR8caMGDGCm266iWOPPbbe+d555x3uuusunnzySer4fnoFMMPdE77tLFIp8AlgRCNv/S1wibtvjzNviyVQkfnPIcR9UH1ja5nl7v8nkbFUa2wCpS18IiIiItJqmdlxZvYoIZkZGWvMkCFDePzxx1m6dGmDkieAo48+mkWLFvHb3/6W4cOHxxt2OvBnM3vGzGKueDWVu28k/H7GA+814JY3CInTP8VLnpLB3Z8FTgGeJjTkrc8bwOiWSp6aQitQIiIiItLqmFl/YDowgTiLBt/4xjf48Y9/zKhRo5r9vNLSUmbNmsUbb7wRb4gDjwEz3f3dZj+wFjM7DjgHOIJQ3MKAL4C3gRfc/bUGznMWNWeOXnb3LYmONepZfQkFIo4jFPDoBPwD+DTy63l3f7Olnh8VRxbwrahLf3T3uMmdEigRERERaTUi35TfCFwOxKyKN2DAAK6//nrGjBlDHbUXGs3defbZZ5kzZw6rV6+ON6wSWEw4I7Uu3iBJX0qgRERERCTjRc4EXQ9M4sAGuAD06dOHKVOmcPHFF9OuXctVHK+qquKpp57itttu46OPPoo3bDeh6MNN7v5ZiwUjCacESkREREQylpkdBBQStut1jTXm4IMP5qqrrmLChAnk5iar5y1UVlaydOlS7rzzTjZs2BBv2FagGJjt7l8mLThpMiVQIiIiIpJxzKw9oRnrLUCvWGMOOuggvve97zFp0iQ6d+4ca0hSbN++ncWLF7NgwQK2bIl7pGgjcAdwTzL7M0njKYESERERkYxhZu2AS4GbgSNjjenQoQOXXnopkydPpnv37skMr06bN2+muLiY+++/n5074+ZIa4H/AB6I1adJUk8JlIiIiIikPQvVHsYSkouBscbk5ORQUFDA1KlT6dUr5qJUWli/fj3FxcU8/PDD7N69O96wdwh9qx6pqyKcJJ8SKBERERFJa5Hmq7cDMXspZWVlMXr0aKZPn07//v2TGltzlJeXM3/+fJYsWcLevXvjDXuLULFvWRJDkzoogRIRERGRtGRmpxFWYUbGGzNixAhuuummBjfATUfvvPMOd911F08++SR1fG++Apjh7n9JYmgSgxIoEREREUkrkaawNwEXxRszZMgQpk+fztChQ5MXWAtbtWoVc+bM4YUXXqhr2LNAkbu/mqSwpBYlUCIiIiKSFsysP6Ec+QQgK9aYb3zjG/z4xz9m1KhRSYwsuUpLS5k1axZvvPFGvCEOPAbMdPd3kxeZgBIoEREREUkxM+sL3AhcDsTscDtgwACuv/56xowZQ6gn0bq5O88++yxz5sxh9erV8YZVAosJZ6TWJS+6tk0JlIiIiIikhJkdAlwPTAI6xBrTp08fpkyZwsUXX0y7djFzq1atqqqKp556ittuu42PPvoo3rDdwK+AG909bsdeSQwlUCIiIiKSVGZ2EFBI2K7XNdaYgw8+mKuuuooJEyaQm5ub1PjSUWVlJUuXLuXOO+9kw4a4OdJWoBi4zd3jduyV5lECJSIiIiJJYWbtgcuAW4CYjZo6derEZZddxqRJk+jcuXMyw8sI27dvZ/HixSxYsIAtW+LmSJ8DdwL3uHvcjr3SNEqgRERERKRFmVkW8F1CL6cjY43p0KEDl156KZMnT6Z79+5JjS8Tbd68meLiYu6//3527oybI60lNB5+wN33JC+61k0JlIiIiIi0CAvVHsYSvokfGGtMTk4OBQUFTJ06lV69Yi5KSR3Wr19PcXExDz/8MLt374437B1CP61H3L0qedG1TkqgRERERCThzOwcworTybE+z8rKYvTo0UyfPp3+/fsnNbbWaO3atSxYsIAlS5awd+/eeMPeBG5192VJDK3VUQIlIiIiIgljZqcRVjtGxhszYsQIbrrpJo499tikxdVWvPPOO9x1112UlJTUNWwFMMPd/5KksFoVJVAiIiIi0mxmdhxwE3BRvDFDhgxh+vTpDB06NHmBtVGrVq1izpw5vPDCC3UNexYocvdXkxRWq6AESkRERESazMz6E8qRTwCyYo0ZNGgQU6ZMIT8/P4mRCUBpaSmzZs3ijTfeiDfEgceAme7+bvIiy1xKoERERESk0cysL3AjcDkQs8PtgAEDuP766xkzZgyhnoSkgrvz5JNPMnfuXN5///14wyqBxcAt7r4uedFlHiVQIiIiItJgZnYIcD0wCegQa0yfPn2YMmUKF198Me3axcytJAWqqqp46qmnmDVrFh9//HG8YbuBXwE3unvcjr1tmRIoEREREamXmR0EFBK263WNNebggw/mqquuYsKECeTm5iY1Pmm4yspKli5dyp133smGDXFzpK1AMXCbu8ft2NsWKYESERERkbjMrD1wGXALELNRU6dOnbjsssuYNGkSnTt3TmZ40gzbt29n8eLFLFiwgC1b4uZInwN3Ave4e9yOvW2JEigREREROYCZZQHfJfRyOjLWmA4dOnDppZcyefJkunfvntT4JHE2b95McXEx999/Pzt3xs2R1hIaIj/g7nuSF136UQIlIiIiIvtYqPYwlvDN8sBYY3JycigoKGDq1Kn06hVzUUoy0Pr16ykuLuahhx6isrIy3rC3CeXqH/M2mkgogRIRERERAMzsHMKK08mxPs/KymL06NFMnz6d/v37JzU2SZ61a9eyYMEClixZwt69e+MNexO41d2XJTG0tKAESkRERKSNM7PTgFnAWfHGjBgxghtvvJHjjjsueYFJSr3zzjvcddddlJSULine truncated
}
},
"cell_type": "markdown",
"id": "18a9f60f",
"metadata": {},
"source": [
"<div>\n",
"<img src=\"attachment:fig-asp-efficiency-comm-2.png\" align=\"left\" width=\"250\"/>\n",
"</div>"
]
},
{
"cell_type": "markdown",
"id": "4ec6718c",
"metadata": {},
"source": [
" a) The processes are synchronized in each iteration due to the blocking send and receive of row k.\n",
" b) Receiving processes may overwrite the data in row k, which can lead to incorrect behavior.\n",
" c) The sending process can only continue the computation after the data are received in every other process.\n",
" d) The receiving process does not know the source of the received data."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4f4a57de",
"id": "138c77eb",
"metadata": {},
"outputs": [],
"source": [
"answer = \"x\" # replace x with a, b, c or d\n",
"floyd_impl_check(answer)"
]
},
{
"cell_type": "markdown",
"id": "c624722a",
"metadata": {},
"source": [
"### Testing the implementation"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "09937668",
"metadata": {},
"outputs": [],
"source": [
"function rand_distance_table(n)\n",
" threshold = 0.4\n",
" mincost = 3\n",
" maxcost = 10\n",
" infinity = 10000*maxcost\n",
" C = fill(infinity,n,n)\n",
" for j in 1:n\n",
" for i in 1:n\n",
" if rand() > threshold\n",
" C[i,j] = rand(mincost:maxcost)\n",
" end\n",
"code3 = quote\n",
" function floyd_iterations!(myC,comm)\n",
" L = size(myC,1)\n",
" N = size(myC,2)\n",
" rank = MPI.Comm_rank(comm)\n",
" P = MPI.Comm_size(comm)\n",
" lb = L*rank+1\n",
" ub = L*(rank+1)\n",
" C_k = similar(C,N)\n",
" for k in 1:N\n",
" if (lb<=k) && (k<=ub)\n",
" # Send row k to other workers if I have it\n",
" myk = (k-lb)+1\n",
" C_k[:] = view(myC,myk,:)\n",
" for dest in 0:(P-1)\n",
" if rank == dest\n",
" continue\n",
" end\n",
" MPI.Send(C_k,comm;dest)\n",
" end\n",
" else\n",
" # Wait until row k is received\n",
" MPI.Recv!(C_k,comm,source=MPI.ANY_SOURCE)\n",
" end\n",
" # Now, we have the data dependencies and\n",
" # we can do the updates locally\n",
" for j in 1:N\n",
" for i in 1:L\n",
" myC[i,j] = min(myC[i,j],myC[i,k]+C_k[j])\n",
" end\n",
" end\n",
" end\n",
" myC\n",
" end\n",
" C[j,j] = 0\n",
" end\n",
" C\n",
"end"
"end;"
]
},
{
"cell_type": "markdown",
"id": "ad1c3fca",
"metadata": {},
"source": [
"### Collecting back the results\n",
"\n",
"At this point, we have solved the ASP problem, but the solution is cut in different pieces, each one stored on a different MPI rank. It is often useful to gather the solution into a single matrix, e.g., to compare it against the sequential algorithm.The following function collects all pieces and stores them in $C$ on rank 0. Again, we implement this with `MPI.Send` and `MPI.Recv!` as it is easier as we are working with a row-partition. However, we could do it also with `MPI.Gather!`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "dd77ee3d",
"id": "c12fd15d",
"metadata": {},
"outputs": [],
"source": [
"using Test\n",
"load = 10\n",
"n = nworkers()*load\n",
"C = rand_distance_table(n)\n",
"C_seq = floyd!(copy(C))\n",
"C_par = floyd_dist!(copy(C))\n",
"@test C_seq == C_par"
"code4 = quote\n",
" function collect_result!(C,myC,comm)\n",
" L = size(myC,1)\n",
" rank = MPI.Comm_rank(comm)\n",
" P = MPI.Comm_size(comm)\n",
" if rank == 0\n",
" lb = L*rank+1\n",
" ub = L*(rank+1)\n",
" C[lb:ub,:] = myC\n",
" for source in 1:(P-1)\n",
" lb = L*source+1\n",
" ub = L*(source+1)\n",
" MPI.Recv!(view(C,lb:ub,:),comm;source)\n",
" end\n",
" else\n",
" dest = 0\n",
" MPI.Send(myC,comm;dest)\n",
" end\n",
" C\n",
" end\n",
"end;"
]
},
{
"cell_type": "markdown",
"id": "cbca58aa",
"metadata": {},
"source": [
"### Running and testing the code\n",
"\n",
"In the cell below, we run the parallel code and compare it against the sequential. Note that we can only compare both results on rank 0 since this is the only one that contains the result of the parallel code. We have also included a function that generates random distance tables of a given size $N$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2daef9be",
"metadata": {},
"outputs": [],
"source": [
"code = quote\n",
" using MPI\n",
" MPI.Init()\n",
" $code1\n",
" $code2\n",
" $code3\n",
" $code4\n",
" function input_distance_table(n)\n",
" threshold = 0.1\n",
" mincost = 3\n",
" maxcost = 9\n",
" inf = 10000\n",
" C = fill(inf,n,n)\n",
" for j in 1:n\n",
" for i in 1:n\n",
" if rand() > threshold\n",
" C[i,j] = rand(mincost:maxcost)\n",
" end\n",
" end\n",
" C[j,j] = 0\n",
" end\n",
" C\n",
" end\n",
" function floyd!(C)\n",
" n = size(C,1)\n",
" @assert size(C,2) == n\n",
" for k in 1:n\n",
" for j in 1:n\n",
" for i in 1:n\n",
" @inbounds C[i,j] = min(C[i,j],C[i,k]+C[k,j])\n",
" end\n",
" end\n",
" end\n",
" C\n",
" end\n",
" comm = MPI.Comm_dup(MPI.COMM_WORLD)\n",
" rank = MPI.Comm_rank(comm)\n",
" if rank == 0\n",
" N = 24\n",
" else\n",
" N = 0\n",
" end\n",
" C = input_distance_table(N)\n",
" C_par = copy(C)\n",
" floyd_mpi!(C_par,comm)\n",
" if rank == 0\n",
" C_seq = copy(C)\n",
" floyd!(C_seq)\n",
" if C_seq == C_par\n",
" println(\"Test passed 🥳\")\n",
" else\n",
" println(\"Test failed\")\n",
" end\n",
" end\n",
"end\n",
"run(`$(mpiexec()) -np 3 julia --project=. -e $code`);"
]
},
{
@@ -909,10 +971,12 @@
"id": "91a772df",
"metadata": {},
"source": [
"### Is this implementation correct?\n",
"## Is this implementation correct?\n",
"\n",
"Point-to-point messages are *non-overtaking* (i.e. FIFO order) between the specified sender and receiver according to section 3.5 of the MPI standard 4.0.\n",
"Unfortunately this is not enough in this case. The messages can still arrive in the wrong order if messages from different processes overtake each other."
"In the cell above, the result of the parallel code was provably the same as for the sequential code. However, is this sufficient to assert that the code is correct? Unfortunately no. In fact, the parallel code we implemented is not correct. There is no guarantee that this code computes the correct result. This is why:\n",
"\n",
"In MPI, point-to-point messages are *non-overtaking* between a given sender and receiver. Say that process 1 sends several messages to process 3. All these will arrive in FIFO order. This is according to section 3.5 of the MPI standard 4.0.\n",
"Unfortunately, this is not enough in our case. The messages could arrive in the wrong order *from different senders*. If process 1 sends messages to process 3, and then process 2 sends other messages to process 3, it is not granted that process 3 will receive first the messages from process 1 and then from process 2 (see figure below)."
]
},
{
@@ -920,7 +984,7 @@
"id": "1277772d",
"metadata": {},
"source": [
"If we are lucky all messages will arrive in order and we will process all rows in the right order in all processors."
"If we are lucky all messages will arrive in order. In our parallel code, all processors would receive first row one, then row 2, then row 3, etc. The computed result will be correct."
]
},
{
@@ -943,7 +1007,7 @@
"id": "df60e4e7",
"metadata": {},
"source": [
"However, FIFO ordering is not enough. In the next figure, communication between process 1 and process 3 is particularly slow. Note that process 3 receives messages from process 1 after it receives the messages from 2 even though FIFO ordering is satisfied between any two processors."
"However, FIFO ordering between pairs of processors is not enough to guarantee that rows arrive in consecutive order. The next figure shows a counter example. In this case, communication between process 1 and process 3 is particularly slow for some unknown reason. As result processor 3 receives first messages from processor 2, even though processor 1 sent the messages first. In our parallel code, the received rows would be first row 3, then row 4, then row 1, then row 2, which is not correct. Note however that process 3 received all messages from process 1 in the correct order (guaranteed by MPI). But this is not enough in our algorithm."
]
},
{
@@ -968,10 +1032,25 @@
"source": [
"### Possible solutions\n",
"\n",
"There are several solution to this synchronization problem:\n",
"\n",
"1. **Synchronous sends**: Use synchronous send MPI_SSEND. This is less efficient because we spend time waiting until each message is received. Note that the blocking send MPI_SEND used above does not guarantee that the message was received. \n",
"2. **MPI.Barrier**: Use a barrier at the end of each iteration over $k$. This is easy to implement, but we get a synchronization overhead.\n",
"3. **Order incoming messages**: The receiver orders the incoming messages, e.g. according to MPI.Status or the sender rank. This requires buffering and extra user code.\n",
"4. **MPI.Bcast!**: Communicate row k using `MPI.Bcast!`. One needs to know which are the rows owned by the other ranks."
"3. **Order incoming messages**: The receiver orders the incoming messages, e.g. using to MPI.Status to get the sender rank. This requires buffering and extra user code.\n",
"4. **MPI.Bcast!**: Communicate row k using `MPI.Bcast!`. One needs to know which are the rows owned by the other ranks since we cannot use `MPI.ANY_SOURCE` in `MPI.Bcast!`. This is trivial however if the number of rows is multiple of the number of ranks."
]
},
{
"cell_type": "markdown",
"id": "7c8a16e4",
"metadata": {},
"source": [
"## Summary\n",
"\n",
"- We learned how to parallelize Floyd's algorithm.\n",
"- The considered strategy based on a row-wise data partition has little communication overhead if the problem size is large enough.\n",
"- One needs to be careful in which order the messages are received to have a correct algorithm.\n",
"- There are several strategies to solve this synchronization problem, each one with pros and cons.\n"
]
},
{
@@ -979,61 +1058,11 @@
"id": "de96ad1b",
"metadata": {},
"source": [
"## Exercise \n",
"Rewrite the worker code of the parallel ASP algorithm so it runs correctly. Use the `MPI.Bcast!` to solve the problem of overtaking messages. Note: Only use `MPI.Bcast!`, do not use other MPI directives in addition. You can test your function with the following code cell. "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "31194529",
"metadata": {},
"outputs": [],
"source": [
"function floyd_par!(C,N)\n",
" comm = MPI.Comm_dup(MPI.COMM_WORLD)\n",
" nranks = MPI.Comm_size(comm)\n",
" rank = MPI.Comm_rank(comm)\n",
" T = eltype(C)\n",
" if rank == 0\n",
" buffer_root = Vector{T}(undef,N*N)\n",
" buffer_root[:] = transpose(C)[:]\n",
" else\n",
" buffer_root = Vector{T}(undef,0)\n",
" end \n",
" Nw = div(N,nranks)\n",
" buffer = Vector{T}(undef,Nw*N)\n",
" MPI.Scatter!(buffer_root,buffer,comm;root=0)\n",
" Cw = Matrix{T}(undef,Nw,N)\n",
" transpose(Cw)[:] = buffer\n",
" MPI.Barrier(comm)\n",
" floyd_worker_bcast!(Cw,comm)\n",
" buffer[:] = transpose(Cw)[:]\n",
" MPI.Gather!(buffer,buffer_root,comm;root=0)\n",
" if rank == 0\n",
" transpose(C)[:] = buffer_root[:]\n",
" end\n",
" C\n",
"end\n",
"## Exercise\n",
"\n",
"@everywhere function floyd_worker_bcast!(Cw,comm)\n",
" # Your implementation here\n",
"end\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1b7eb4c2",
"metadata": {},
"outputs": [],
"source": [
"load = 10\n",
"n = nworkers()*load\n",
"C = rand_distance_table(n)\n",
"C_seq = floyd!(copy(C))\n",
"C_par = floyd_par!(copy(C),n)\n",
"@test C_seq == C_par"
"### Exercise 1\n",
"\n",
"Modify the `floyd_iterations!` function so that it is guaranteed that the result is computed correctly. Use `MPI.Bcast!` to solve the synchronization problem. Note: only use `MPI.Bcast!`in `floyd_iterations!`, do not use other MPI directives. "
]
},
{
notebooks/figures/fig_jacobi.svg
+1212 -1208
View File
File diff suppressed because it is too large. Load diff
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 6.6 MiB

After

Width:  |  Height:  |  Size: 6.6 MiB