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

Changes in jacobi method

Browse Source
This commit is contained in:
Francesc Verdugo committed 2024-09-09 13:15:34 +02:00
1 parent 86e2f88e7b
commit 4a9dc6f2c4
2 files changed
+130360 -129868
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

notebooks/figures/fig_jacobi.svg
+129779 -129766
View File
File diff suppressed because it is too large. Load diff
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 6.5 MiB

After

Width:  |  Height:  |  Size: 6.1 MiB

notebooks/jacobi_method.ipynb
+581 -102
View File
@@ -45,10 +45,18 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 104,
"id": "1dc78750",
"metadata": {},
"outputs": [],
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"🥳 Well done! \n"
]
}
],
"source": [
"using Printf\n",
"\n",
@@ -63,6 +71,24 @@
"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\")\n",
"sndrcv_check(answer) = answer_checker(answer,\"b\")\n",
"function sndrcv_fix_answer()\n",
" msg = \"\"\"\n",
" One needs to carefully order the sends and the receives to avoid cyclic dependencies\n",
" that might result in deadlocks. The actual implementation is left as an exercise. \n",
" \"\"\"\n",
" println(msg)\n",
"end\n",
"jacobitest_check(answer) = answer_checker(answer,\"a\")\n",
"function jacobitest_why()\n",
" msg = \"\"\"\n",
" The test will pass. The parallel implementation does exactly the same operations\n",
" in exactly the same order than the sequential one. Thus, the result should be\n",
" exactly the same. Note however this is often not true in other parallel algorithms.\n",
" Specially, when one does parallel reductions.\n",
" \"\"\"\n",
" println(msg)\n",
"end\n",
"println(\"🥳 Well done! \")"
]
},
@@ -526,45 +552,73 @@
"source": [
"## MPI implementation\n",
"\n",
"We consider the implementation using MPI. The programming model of MPI is generally better suited for data-parallel algorithms like this one than the task-based model provided by Distributed.jl. In any case, one can also implement it using Distributed.jl, but it requires some extra effort to setup the remote channels right for the communication between neighbor processes.\n",
" \n",
"Take a look at the implementation below and try to understand it.\n"
"We consider the implementation of the Jacobi method using MPI. The programming model of MPI is generally better suited for data-parallel algorithms like this one than the task-based model provided by Distributed.jl. In any case, one can also implement it using Distributed.jl, but it requires some extra effort to setup the remote channels right for the communication between neighbor processes.\n",
"\n",
"We are going to implement the method in a function with the following signature\n",
"\n",
"```julia\n",
" u = jacobi_mpi(n,niters,comm)\n",
"```\n",
"\n",
"The signature will be very similar to the sequential function `jacobi`, but there are some important differences:\n",
"\n",
"1. The parallel one 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",
"2. Function `jacobi_mpi` will be called on multiple MPI ranks. Thus, the values of `n` and `niters` should be the same in all ranks. The caller of this function has to make sure that this is true. Otherwise, the parallel code will not work.\n",
"3. The result `u` is NOT the same as the result of the sequential function `jacobi`. It will contained only the local portion of `u` stored at each rank. To recover the same vector as in the sequential case, one needs to gather these pieces in a single rank. We will do this later.\n",
"\n",
"The implementation of the parallel function is given in the following code snipped:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e15082fb",
"execution_count": 53,
"id": "77088426",
"metadata": {},
"outputs": [],
"source": [
"] add MPI"
"code1 = quote\n",
" function jacobi_mpi(n,niters,comm)\n",
" # Initialization\n",
" u, u_new = init(n,comm)\n",
" for t in 1:niters\n",
" # Communication\n",
" ghost_exchange!(u,comm)\n",
" # Local computation\n",
" local_update!(u,u_new)\n",
" u, u_new = u_new, u\n",
" end\n",
" return u\n",
" end\n",
"end;"
]
},
{
"cell_type": "markdown",
"id": "cae5af50",
"metadata": {},
"source": [
"In the following cells, we discuss the implementation of the helper functions `init`, `ghost_exchange`, and `local_update`."
]
},
{
"cell_type": "markdown",
"id": "5da3e049",
"metadata": {},
"source": [
"### Initialization\n",
"\n",
"Let us start with function `init`. This is its implementation:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a66cbf9a",
"execution_count": 52,
"id": "ee17f3cb",
"metadata": {},
"outputs": [],
"source": [
"using MPI"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d39c7bb2",
"metadata": {
"code_folding": []
},
"outputs": [],
"source": [
"code = quote\n",
" using MPI\n",
" MPI.Init()\n",
" function jacobi_mpi(n,niters)\n",
" comm = MPI.COMM_WORLD\n",
"code2 = quote\n",
" function init(n,comm)\n",
" nranks = MPI.Comm_size(comm)\n",
" rank = MPI.Comm_rank(comm)\n",
" if mod(n,nranks) != 0\n",
@@ -573,53 +627,409 @@
" end\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",
" # Communication\n",
" if rank != 0\n",
" neig_rank = rank-1\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 = 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",
" # 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",
" if rank == 0\n",
" u[1] = -1\n",
" end\n",
" # Gather the results\n",
" if rank == nranks-1\n",
" u[end] = 1\n",
" end\n",
" u_new = copy(u)\n",
" u, u_new\n",
" end\n",
"end;"
]
},
{
"cell_type": "markdown",
"id": "dee2b3b9",
"metadata": {},
"source": [
"This function crates and initializes the vector `u` and the auxiliary vector `u_new` and fills in the boundary values. Note that we are not creating the full arrays like in the sequential case. We are only creating the parts to be managed by the current rank. To this end, we start by computing the number of entries to be updated in this rank, i.e., variable `load`. We have assumed that `n` is a multiple of the number of ranks for simplicity. If this is not the case, we stop the computation with `MPI.Abort`. Note that we are allocating two extra elements in `u` (and `u_new`) for the ghost cells or boundary conditions. The following figure displays the arrays created for `n==9` and `nranks==3` (thus `load == 3`). Note that the first and last elements of the arrays are displayed with gray edges denoting that they are the extra elements allocated for ghost cells or boundary conditions."
]
},
{
"attachments": {
"fig.png": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA6gAAAFzCAYAAAAzJyTtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAewgAAHsIBbtB1PgAAABl0RVh0U29mdHdhcmUAd3d3Lmlua3NjYXBlLm9yZ5vuPBoAACAASURBVHic7d15mBxVvf/xd08mmewkISEEIoTFsCMSBWRT2d0uCIr7DoGrXnC/qCigV38oCoKKTBCv4Mp1AwHZjIIiyBIBZQmIAUNICISELJBJMjP9++N0p2tmumeyzPQ5PfV+PU8/nOruTH/o7uqqb9Wpcwqtra1FJEmSJEmKrCl2AEmSJEmSwAJVkiRJkpSI5m7LRxcKhSdiBJEkSZIk5UuxWJwG3Fhe7lKgFgqFJ2bOnPlovUNJkiRJkvJn1qxZFIuVYZHs4itJkiRJSoIFqiRJkiQpCRaokiRJkqQkWKBKkiRJkpJggSpJkiRJSoIFqiRJkiQpCRaokiRJkqQkWKBKkiRJkpJggSpJkiRJSoIFqiRJkiQpCRaokiRJkqQkWKBKkiRJkpJggSpJkiRJSoIFqiRJkiQpCRaokiRJkqQkWKBKkiRJkpJggSpJkiRJSoIFqiRJkiQpCRaokiRJkqQkWKBKkiRJkpJggSpJkiRJSoIFqiRJkiQpCRaokiRJkqQkWKBKkiRJkpJggSpJkiRJSoIFqiRJkiQpCRaokiRJkqQkWKBKkiRJkpJggSpJkiRJSkJz7ACS6ueyyy4b097e/ufYOXpz8cUXb7948eIWgFWrVi1ZtWrV0tiZupkITCi1XwQWRMxSzUhgamb50VhBejE9015AeB9TMpXwPgIsBZZEzFLNBML3EKANmB8xSzXDge0yy48BnZGy1LIzlYP0C4FVEbNUsw0wutR+HngmYpZqxgFbldprgSfiRempubl52MSJE6eVlz/+8Y/PGzduXHvESL0qFotfOPXUU6+JnUNKhQWqlCNNTU1DgJfFztGbZcuW8fTTT5cXX1K6pWoEsGXsEH1I+vMGXho7QB+2Ld1SNQIYHztEH/aKHaAPO8UO0IcRwJTYIXoxgsR+Z9rb27PbEYrF4h4R4/SpUChM6PtZUn7YxVeSJEmSlATPoEo5ViwWL2xqakqqi2pnZ+dXgaEAU6ZMeWTRokXfjxypu/+i0n3xReCsiFmqeS3w+szyp2MFqaEJ+Fpm+Vrg1khZavkS4awQwL+B70TMUs1MKmee1wKfj5ilmv2At2aWvwysiJSllvMy7T8Cv4sVpIYvAmNK7UXA+RGzVPMeYO9SuwM4I2KWHiZMmLDr0qVLP1ReXrly5UUTJkx4Mmam7orF4heAsbFzSCmyQJXy7cczZ868J3aIrClTppxDqUCdOnXqvEWLFn0jcqTu3kKlQG0DUss3hK4Famr5mulaoN4CfDNOlJo+R6VAXUh67+HRVArUdaSX74N0LVC/RyiyUvJ1oFBq30V67+EnqRSoS0gv3350LVCTyjdjxozX33zzzesL1Dlz5vz4K1/5yt0xM3XX2tr6cSxQpars4itJkiRJSoIFqiRJkiQpCRaokiRJkqQkWKBKkiRJkpJggSpJkiRJSoIFqiRJkiQpCRaokjZXM7AX8AHCfJF/Bu4p3T5Z5ywFwjygrcADwEqgE3gauAb4jzrnqWYL4GTgt8ACYA2wGngE+DawQ7xo6+0BnAvcASwFisAywmf7n8CweNGAMA3RccCPgccI89GuA54Efgq8Kl609aYAnwJuBhYD7cALwH2EeVYnxou23kGEdfZ+wjylReAZ4HrgRCrTsMQyGngf8GtgPmFapzbCZ34JsEu8aOtNJ3yetxGmgykCy4HbgY9Tma4oliHAG4AfEn5jXiB8F58CfkH4vYxtEnA64Xv3NGFdfhH4B+F3aOtIuYYC+xJ+ry8hfKblbdsHImWScsF5UCVtjv8kzH83ssbjf65jlh2AK4FXVnlsMvDG0u3HwPsJc/fV2/uB84HxVR6bXrp9AHgXcHX9Yq3XAvwEOJ6exck44ODS7YPA6wg75PX2KuBHwE5VHpsKvAN4O3A2oXCI4RzgvwnvZ9ZI4GWl28mEwuFv9Y0GwDaEdeXgKo9NAo4p3d5JKFTX1i/aem8FvlvK091Opdv7gA8RDkrUW4HwPXwHPQ/2jyV8T19FyHcUYT7detuH8N7sVuWxbQhzOr+FMA/xp+qYK+tThHV1VLf7m4E9S7eTCAekbqtjrjNLt+7rcNmUOmaRcscCVdLmmEzP4vR5QjFTb7tRKU4fBq4lnGkZCxxGKKgA3l26/5x6BwTeRChO1wG/B24hnKHcHngvsB1hR+3nhB2zf9U530jghFL7aUKR/I/S8isIhfPQUvtHVN7TejqISnF6J3Aj4czpJMJO7H6E4uEcwvfgFxEynkjYsV0F/I5wJno1lQMQ4wlnhX4L7Fp6Xj3tQKU4/Rfhc36U8N07BDiW8B4eSziD9Yk654NQIE8iHEj6Q+m2hHAQ4l3AzsBw4HLgQcJZ4HoqlHJQynV1KUMH4QDE+wjfgT0I38GDCWdX6+mVVIrTewnfxSeALQkHRw4pPfZJYC7w/Trng/B7M4qwflwP/IWwPuxIOBA2iZD3KsK6Uq+DYlPoWpwWCb0MtqjT60u5ZoEqaXM8Q+h+N6d0uwd4DfDLSHn+AXwGuKHb/d8APkLozgihO9lXCF3d6qkIXAr8D6HLYtZXCUXrgYQd748QpzB4hnDm4HJ6njlrJRTVLYQCYk9CV+p6mw2cQfi+ZZ0LXAR8tLT8SeIUqGsIn/EFhAMQWV8H7iIcjNiWcLY3RmHwKGFduYbQDb7sfELh9ePS8qnAFwhdQ+vtx4Sza90P1HyVkPsown7M6YRipt6WAmcBlxEKrKzvEoqt0YR1+sDScr39hbCudD/7+DXCb+DnSsufJM73sJ3w+3we4bcn61zgr4Su3FsSeqB8o065FhJ6GZS3bX8jnC2/uE6vL+Wa16BK2hwXE46Af5VwJuu5iFn+CsygZ3Fa9l3g36X2eMIR+nr7CDCTnsUphB3cszLLM+qSqKsXCGdcLqV6t86/Eg5IlO1bj1DdXAEcQc/iFMIBgM9RKbheTpzrKI8mFHXdi1MI16N+M7Mc43N+kHCW72q6FqdlP6FyRnIE1buIDrTPAe+hei+CtYSDKGUx3sMi4X35Dj2LU4C/UynyIc66cjVwKLW7xp5NuNYTwtnJ7t1s6+EtwKfpWZxC6I3z/zLL9fycv0I4eHQe4ez983V8bSn3LFAlDRZLCV1ne/NUPYL0YnEfj2fzxSis1lK9qMqKnbHajmzWSkJXvLIYGfv6nBdk2jHyPU8YbKg3sT/n1NeV8oBSvYmd8RmqH4AoW0fX/wc/Z0lJsECVlBfbUDkC/yzweMQstbwp074zWorahgGvzyzfFStILw6hcg30PfS+gx5LdjTpFD/niVSuUV0FPBQxSy2prysFwqBsZSlm3Idw/TuEs+r1vhZ6Q6T+OUsaABaokvKghdDdrjzoxdfp+2xrvb0C+GKpvYLK9bIp+Tawe6n9a8IgRCmZRLgesOyrsYL04q2EAbEgdF+9MmKWaoYC/0cYXAzgW8S5/rQ3u1Pp+tlGuG42NWcC+5fas0mvuJpAuD67fFbyKxGz1HI08OFSeyFhqhxJOWCBKmmwawF+RWW+v6tIb4d2X8I1vKMIo4B+gMr1sqn4EuH6WQgjfp4SMUs1kwhzjr60tPw14Lp4cao6lnCgpEA4CPE2KtcApqCZMDpzeV35PfDleHGq2o2Qazyhm+2ppHeg5JNUpjj6N2FE35RsQVg3di4tXwL8LF6cqg4nHARrJlzj+3bCXMyScsBRfCUNZkMJOznlbqk3EkZiTKnb58uo7HB3EOZN/HWv/6L+ziYM+gPhrN9RxJkDtZYJhLNUe5WWLwY+Gy9OVW8inLEaSuhK+SbC6KCpGEIont9WWr4NeDNx5kCtZTrwR8L0VkXC6L2XR03U08epjDS7gDCgV+xr37NGE34Hy2d3f0wYvC0lryWM0jyCMCL2W6jvnNqSIrNAlTSYtVIpTq8HjqfvwWHqaRvCqMPjCUXz++k68mcKTqYyuvBjhJ3HBbWfXnfNhJ3ZcnH6HeA06j/nZG9eSejKO5TQXfZ11B5ZNZZvUClO/0SYJzOlaxInADcRilOA/yKMzJ2St1EZoXkBYcqtes9l3JsC4SBJuTi9gjA9T0oH7HYnjD5cLk6PJRTUknLELr6SBquDCV1lAe4mHIVPqTiFMFfm1qX2p0ivOB1H5WzQEsI1YSkVpxB2sA8sta8kveIUwrW7IwiFwFtIrzjdm/C+QRgs5z9IqziFcAa/PKDPl0ivOB1O+JzL3bePIa3iFML1z8eU2tcTemt0xItT1QXAmFL7fVicSrnkGVRJWV8ijNRazSLgwjpm2VzHZtqfJa1r/cqOK/13Hmm+t4dRGSznQkLO1GQ/50+QXnE6hcoZqxuoPU9vTP9B5YD1WcDyiFlqKX/OzxIO7KTmIMJ10BDmEX4wYpZajsu0Pwm0xwpSwxjCtacAt5PeAGKS6sQCVVLWJ6g9Wfv9pFlE1fKSTDu1QVQgvM/jS+1HSKubXdnUTDvF9xAqGZ8njPSZmkZ6DyHNKWUKVDI+RnojcEP67yFUMnYSBjpLzdaEa6Eh3XVFUh1YoErKupfQFbGaf9YzSD9oybTXREtR2/BMO7Wux2WNlDHFzxga4z3MrispZmwiXL8LaeaDxvqc15JeTwNojPdQUh1YoErKOiR2gAGS4s5YVur5GkEjvIeNkDF1voebrxHew0bIKGmAWKBKGqxOJVxnBaH7Z2qeB3YqtV+IGaQXs4BfltpPxwzSiyMIZ9dSG+yl7G4qn3OK30OAz1CZ7/TJmEFq6KDyHq6OGaQXPyPMwwvwTMwgvTiecBY11eLvESqf84qYQSTFZYEqabBaHDtAHzpIc9ChrOdJt6gqS7Ggymoj/c/52dItZam/hytIv6hKaT7WataS/ucsqQ6cZkbSLine truncated
}
},
"cell_type": "markdown",
"id": "5885e6e7",
"metadata": {},
"source": [
"<div>\n",
"<img src=\"attachment:fig.png\" align=\"left\" width=\"450\"/>\n",
"</div>"
]
},
{
"cell_type": "markdown",
"id": "0efe5f7a",
"metadata": {},
"source": [
"### Communication\n",
"\n",
"Once the working arrays have been created, we can start with the iterations of the Jacobi method. Remember that this is implemented in parallel by updating the ghost cells with data from neighboring processes and then performing the sequential Jacobi update. See figure:"
]
},
{
"attachments": {
"fig15.png": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA6gAAAFACAYAAACxyVHuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAewgAAHsIBbtB1PgAAABl0RVh0U29mdHdhcmUAd3d3Lmlua3NjYXBlLm9yZ5vuPBoAACAASURBVHic7N15fFT19f/x12cmIQFkE1wANxBFrTu4b1itimQC2mKrtfrtRtXazbbWrmJba+1m7S5fv9Wqrbb5VWBugFatpe4buGvdFVcU2cISSOZ+fn+cO2QyTCAsmXszeT8fj/vIvXPv3Dm5M3PnnvvZ3DXXXOMRERERERERiVkq7gBEREREREREQAmqiIiIiIiIJERV0fLJzrlX4whEREREREREehbv/W7AP/PL7RJU59yrU6ZMeb7cQYmIiIiIiEjPM23aNLxv6xZJVXxFREREREQkEZSgioiIiIiISCIoQRUREREREZFEUIIqIiIiIiIiiaAEVURERERERBJBCaqIiIiIiIgkghJUERERERERSQQlqCIiIiIiIpIISlBFREREREQkEZSgioiIiIiISCIoQRUREREREZFEUIIqIiIiIiIiiaAEVURERERERBKhKu4ARDbGe+8aGxsvizuODbnuuusOWLhw4SCAhQsXvvbiiy++EndMRXYBRkbzy4H5McZSSn/g4ILlu4FcTLGUdPjhhx+dTqerAE455ZSnDjzwwEVxx9QR59yDdXV1s8rxWtls9nPOueHleK3NMXv27BFPPPHErgArV65c8thjjz0ed0xFBgEHFCzPjSmODRlXMP84sCSmODpyAHYcARYAL8cYy3pGjBixy7Bhw0YCbLvttsunTJmStPNvO+l0+vJTTz11TVe/ThAEewFndfXrbK7m5ub0L3/5y2Pyy88+++z8xYsXL48zphIOxn4/wT73C2KMpZQRwK7R/FLgsRhjWc/gwYMH7rXXXgfmly+55JK5MYazUWEY3jBx4sQX446jHJSgSuJddtllbsyYMd+NO44NeeGFF3jyySfzi8dsaNuEmBB3ABuRuGP4wAMPrJsfN25c4uIr5L3/LVCWBBX4DDC2TK+1yd577z3uvffewofq44qlkxL92SL58SXOK6+8wiuv2D3LUaNGQcLPv9XV1T8HujxB9d6Pds4l+re96NyR9M9+0uMDyMQdQKH333+/3XschuExqVRyK5em0+l7gR6RoCb3XRAREREREZEeRSWo0h09jlUVSYwwDI8i+j717j1w5erVPBdzSEVW7AGt/Ww+HUK/RFWzgTVDYPUubcsDnwRaYgunBOeWHey9B2DNmjULgKRV494XGBxzDK+QsCpma9euPYioClxNTZ/cmjW9ElbFt3lHaB7WtjwwgdU/lxZUv+/7JlQvjC+WUpYdAD5t871WQp9EnX/79HGjVq1a0h/Ae58D7ok5pGL9aN/EIg4h1rQjMcIw7AUckV+uqdl2wZo1YcKadiw7EHxU2FTTBL1fiDeeYk17Qm4bm0/noF+izr+1tentm5vf3ym/7L2/C/AxhlTKcXEHEAclqNLteO8vqq+vvzPuOAqNGjVqCTAQYOTII196+ulZ58ccUpHBv4bFh9p87RpYkrD4jq+Hud9uW/7zxXDq4vjiWV91de8H165tBuCpp576y89+9rNvxhxSO42NjTO89xNjDuN/M5nMFTHH0M6RRx55P3A4wHbbjVj1xhtPJeyzP+YTMP/CtuWkfTcB3INt80dOh9tujC+WUqpvh9aoHd6Or8BriTqGu+/+4auffPLWwwFaWlqaM5nMuJhDaiebzR7inHso5jASd1zOPPPMHYB38sujR3/xxieeuDQbY0glVM2FXG+b3+m/8OKFG9y87AZNg6VRG/u+q5N2fttrr0vOfOyxK7+cX/7zn//8oYaGhrVxxlQsCIIcPbDGa4/7h0VERERERCSZlKCKiIiIiIhIIihBFRERERERkURQgioiIiIiIiKJoARVREREREREEkEJqoiIiIiIiCSChpmRnqwa2A/YPVp+Bng6vnBERERERHo2JajS0+wBXASMAfYHagrWfR+4NI6gRERERKTbGhlNAK8Dz8UYS7enBFV6mn2B8wqWffTXxRCLiIiIiHQ/w4AvYAUeY4FBBeuuBr4cR1CVQgmq9DSLgZuBedE0H3iV9icWEREREZGO7AVcEncQlUoJqvQ0/4kmEREREZHNsRKYQ1uBxzzgLmC3GGOqGEpQRUREREREOu9B4NS4g6hUGmZGREREREREEkEJqoiIiIiIiCSCElQRERERERFJBCWoIiIiIiIikghKUEVERERERCQRlKBKpZiHjXFaPP0tzqBERERERKTzNMxM5doD+EfB8geA5o08pwZ4pmB5PPD8Vo6rqwwABpV4fJtyByIiIiIiIptHCWrl6gWMLFh2nXiOK3pOr60aUdcaQ+kaAS3lDkRERERERDaPElSpFMviDkBERERERLaM2qCKiIiIiIhIIqgEVXqig2h/cyYd/R2GVRXOewd4s1xBiYiIiIj0dEpQpSe6B+hT4vHPRFPelcAlZYlIRERERLqTXWkr5IC2vGoA7ft0WQ4sKldQlUAJqvRErwC9O7Hd4q4ORERERES6pQeAHUs8/j/RlPcH4PwyxFMxlKBKT7Rv3AGIiIiISLe2DBuicWNWdXUglUYJqoiIiIiIyKbZK+4AKpV68RUREREREZFEUIIqIiIiIiIiiaAEVQp1ph69iIiIiIhIl1CCWrlWFy13ptfaUj2RiYiIiIiIlIUS1Mq1rGh5eCeec0RXBCIiIiIiItIZSlAr1/u0HxT46I1s74DPdV04IiIiIiIiG6YEtbLdVzD/eaB6A9teBBzeteGIiIiIiIh0TAlqZftTwfwHgFuAIUXbDAGuBn4KLC1TXCIiIiIiIuupijsA6VIzgLuAY6Pl04EJwCNYG9XtgQOxz0ELcBYwu/xhbhrn3KFBENTGHUehr33ta+u+S8uXv90PPnZknPGsb+3AtvlcKnnxvT6q/fLvDoEbmuKJpbQwDNfNb7/99iODIDg1xnDW473fIe4YgNFJOy5XXXXVus9+c3NTOnmf/UW7tl9OWnzFXt81eTH6gmuZVYk7/y5b9vqg/HwqlUon7TsC7BF3AEDijsvChQsH3nzzzeuWFy26f1TSPlvgCwqaVg5IXnxr+7fNtybu/Lto0ZIRhcuTJ08+5ZxzzmmNKx5powS1soXAGcA/gQOix2qAo4q2WwKcC/y7fKFtkSviDqBYr1691s2//vqjI+DRq2IMZyOaa+CvCY4PYNb3446gWGvBT9bOO+98BvbdkvbOjabE2HbbbdfNL1q0oA8sSPhnP+nfzf9mbEqqRbsm7RguWNA2X1VVVQvMii2Y5KohYcdlwIAB7ZbfeuufHwU+Gk80nfHOnkn77Le3qnfS4nvjjfbLNTU1M+OJRIopQa18C4HDsDaoZwIHY1W7Q2AB8Dfgt9F8NTCt4LmLyxqpiIiIiIj0aEpQe4Y1wC+iKQUMBJqwar2FWkhgT76XXnqpb2xsXBN3HBuSTqerq6urUwC5XC4XhmHSqoikafu+h6z/3sfNAb0KlhP3fldXV/fC4gQ7fuEGNo+V975s728qlVrjvU/c+5WXSqWqqqur0wDe+7C1tTVpn/0U7TuwS+KxrCmYT+Jnv5q2PjVyQKLOv6lUKp1Op6sAqqqqknj+baelpcWX6aVyJPPzvk51dfW6z35ra+ta7325jk1nFX72W7FjmiSJvvZwzqWqqqqqC5YT/Xkkee9vl3HXXHPNui+bc270lClTno8zIBEREREREekZpk2btqf3/rn8snrxFRERERERkURQgioiIiIiIiKJoARVRCqZ2/gmIiLr6JwhIhIzJagiUomqgG8Ac9B5TrrGsLgDkK3uKOBRYP+4A5GKNAjoHXcQIt2BLtxEpBLtBnwXOBn4WryhSAUaDTwP3AD0jTkW2XouwcYMvwmojTkWqSwO+1w9jg39JyIboARVRCrRi8AXo/kfAofHGItUlhrgZiwx3QlYHW84shV9Cngb2A/4acyxSGX5CnAqsDM6Z4hslBJUEalUfwRuwcaJuwVVyZSt4yfAQcBC4CySNyaobL73gHOw9/TzwP/EGo1UijHAFdH8F4EnYoxFpFtQgioilew84FlgV2A2MCDecKSbm4JdYHrgM8A78YYjXeAO4AdYlcxpwPh4w5FubjgwHegF/B3433jDEekelKCKSCVbhl1gvoW1LbsVq6IpsqkmAr+L5r8LNMYYi3Sty4DrsNoXfwMOiTcc6aa2BW7DqvU+AXw63nBEug8lqCJS6V7DOktaCnwQK0ntF2tE0t0cg7U7TQN/AC6PNxzpYh74LDAT2AZLMo6ONSLpbvpin599gAVY+9NlsUYk0o0oQRWRnuApYBKwAktSb8PubotszAnYTY3eWKnpF+INR8okB5wN3A0MBP6JqvtK5wyi7abGMiADvBlrRCLdjBJUEekp/oMlp+9jvfreBewea0SSdGdiY+lugyWnZwCtsUYk5bQCq33RCPTBSsQ+GWtEknTDsN+WI4FFwIdQp0gim0wJqoj0JA9j1TVfBz4AzAc+EmtEkkQO+DY2bmF19Pd0NDxET7Qaq32Rb5P6R2z82z5xBiWJNBa4H9gXeAM4FvvNEZFNpARVRHqaZ4EjgHuB/lgnKFdivSyKbAsE2Pi5KeBq4FygJc6gJFY5rIOb70TznwDuAfaIMyhJlPOwz8QuwH+x6r3PxhqRSDemBFVEeqI3gXFYYgpwMVaaenhcAUkiHAM8AkzASs4+CXwZjXUq1nHS5Vib5HewsXAfB6ZiJavSM20L/Bn4PdZD/N+Bw7DO+URkMylBFZGeqhW4BPgwsBCr8nsv8Es0XmpPMwArKZ0LjABexG5WXB9fSJJQ/8Gqct6GdZx1KfBA9Jj0LBmsfelZtP2eTAaWxxmUSCVQgioiPd10YC9gGtb28EvAy8A3gNoY45Ku57COkJ4Bvhgt/wEYgzo2Line truncated
}
},
"cell_type": "markdown",
"id": "f7045b08",
"metadata": {},
"source": [
"<div>\n",
"<img src=\"attachment:fig15.png\" align=\"left\" width=\"450\"/>\n",
"</div>"
]
},
{
"cell_type": "markdown",
"id": "e340cfd0",
"metadata": {},
"source": [
"The communication step happens in function `ghost_exchange!`. This one modifies the ghost cells in the input vector `u` by importing data using MPI point-to-point communication. See the following code:"
]
},
{
"cell_type": "code",
"execution_count": 59,
"id": "f3221768",
"metadata": {},
"outputs": [],
"source": [
"code3 = quote\n",
" function ghost_exchange!(u,comm)\n",
" load = length(u)-2\n",
" rank = MPI.Comm_rank(comm)\n",
" nranks = MPI.Comm_size(comm)\n",
" if rank != 0\n",
" neig_rank = rank-1\n",
" u_snd = view(u,2:2)\n",
" u_rcv = view(u,1:1)\n",
" dest = neig_rank\n",
" source = neig_rank\n",
" MPI.Sendrecv!(u_snd,u_rcv,comm;dest,source)\n",
" end\n",
" if rank != (nranks-1)\n",
" neig_rank = rank+1\n",
" u_snd = view(u,(load+1):(load+1))\n",
" u_rcv = view(u,(load+2):(load+2))\n",
" dest = neig_rank\n",
" source = neig_rank\n",
" MPI.Sendrecv!(u_snd,u_rcv,comm;dest,source)\n",
" end\n",
" end\n",
"end;"
]
},
{
"cell_type": "markdown",
"id": "8948edc4",
"metadata": {},
"source": [
"Note that we have used `MPI.Sendrecv!` to send and receive values."
]
},
{
"cell_type": "markdown",
"id": "186254b6",
"metadata": {},
"source": [
"<div class=\"alert alert-block alert-success\">\n",
"<b>Question:</b> Is this other implementation based on `MPI.Send` and `MPI.Recv!` correct?\n",
"</div>\n",
"\n",
"```julia\n",
" function ghost_exchange!(u,comm)\n",
" load = length(u)-2\n",
" rank = MPI.Comm_rank(comm)\n",
" nranks = MPI.Comm_size(comm)\n",
" if rank != 0\n",
" neig_rank = rank-1\n",
" u_snd = view(u,2:2)\n",
" u_rcv = view(u,1:1)\n",
" dest = neig_rank\n",
" source = neig_rank\n",
" MPI.Send(u_snd,comm;dest)\n",
" MPI.Recv!(u_rcv,comm;source)\n",
" end\n",
" if rank != (nranks-1)\n",
" neig_rank = rank+1\n",
" u_snd = view(u,(load+1):(load+1))\n",
" u_rcv = view(u,(load+2):(load+2))\n",
" dest = neig_rank\n",
" source = neig_rank\n",
" MPI.Send(u_snd,comm;dest)\n",
" MPI.Recv!(u_rcv,comm;source)\n",
" end\n",
" end\n",
"```\n",
"\n",
" a) It is correct.\n",
" c) It is incorrect, but it might provide the right result depending on the MPI implementation.\n",
" b) It is incorrect, and it is guaranteed that it will result in a dead lock.\n",
" d) This implementation does not work when distributing over just a single MPI rank.\n"
]
},
{
"cell_type": "code",
"execution_count": 62,
"id": "77e94001",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"It's not correct. Keep trying! 💪\n"
]
}
],
"source": [
"answer = \"x\" # replace x with a, b, c or d\n",
"sndrcv_check(answer)"
]
},
{
"cell_type": "markdown",
"id": "fd3f87e4",
"metadata": {},
"source": [
"<div class=\"alert alert-block alert-success\">\n",
"<b>Question:</b> How would you fix the implementation above, while still using `MPI.Send` and `MPI.Recv!` instead of `MPI.Sendrecv!` ?\n",
"</div>\n",
"\n",
"Think about the answer. To uncover the the solution run next cell."
]
},
{
"cell_type": "code",
"execution_count": 75,
"id": "16a76ea1",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"One needs to carefully order the sends and the receives to avoid cyclic dependencies\n",
"that might result in deadlocks. The actual implementation is left as an exercise. \n",
"\n"
]
}
],
"source": [
"sndrcv_fix_answer()"
]
},
{
"cell_type": "markdown",
"id": "0bc9ed9d",
"metadata": {},
"source": [
"### Local computation\n",
"\n",
"Once the ghost values have the right values, we can perform the Jacobi update locally at each process. This is done in function `local_update!`. Note that here we only update the data *owned* by the current MPI rank, i.e. we do not modify the ghost values. There is no need to modify the ghost values since they will updated by another rank, i.e. the rank that own the value. In the code this is reflected in the loop over `i`. We do not visit the first nor the last entry in `u_new`."
]
},
{
"cell_type": "code",
"execution_count": 55,
"id": "cfce8c20",
"metadata": {},
"outputs": [],
"source": [
"code4 = quote\n",
" function local_update!(u,u_new)\n",
" load = length(u)-2\n",
" for i in 2:(load+1)\n",
" u_new[i] = 0.5*(u[i-1]+u[i+1])\n",
" end\n",
" end\n",
"end;"
]
},
{
"cell_type": "markdown",
"id": "ffc1b539",
"metadata": {},
"source": [
"### Running the code\n",
"\n",
"Not let us put all pieces together and run the code. If not done yet, install MPI."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "d990bd1c",
"metadata": {},
"outputs": [],
"source": [
"] add MPI"
]
},
{
"cell_type": "markdown",
"id": "9ac0d914",
"metadata": {},
"source": [
"The following cells includes all previous code snippets into a final one. Note that we are eventually calling function `jacobi_mpi` and showing the result vector `u`. Run the following code for 1 MPI rank, then for 2 and 3 MPI ranks. Look into the values of `u`. Does it make sense?"
]
},
{
"cell_type": "code",
"execution_count": 92,
"id": "0926b746",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"u = [-1.0, -0.7999999999999999, -0.5999999999999999, -0.3999999999999999, -0.19999999999999996, 0.0, 0.19999999999999996, 0.3999999999999999, 0.5999999999999999, 0.7999999999999999, 1.0]\n"
]
}
],
"source": [
"using MPI\n",
"code = quote\n",
" using MPI\n",
" MPI.Init()\n",
" $code1\n",
" $code2\n",
" $code3\n",
" $code4\n",
" n = 9\n",
" niters = 200\n",
" comm = MPI.Comm_dup(MPI.COMM_WORLD)\n",
" u = jacobi_mpi(n,niters,comm)\n",
" @show u\n",
"end\n",
"run(`$(mpiexec()) -np 1 julia --project=. -e $code`);"
]
},
{
"cell_type": "markdown",
"id": "d5d5bb40",
"metadata": {},
"source": [
"### Checking the result\n",
"\n",
"Checking the result visually is not enough in general. To check the parallel implementation we want to compare the result against the sequential implementation. The way we do the computations (either in parallel or sequential) should not affect the result. However, how can we compare the sequential and the parallel result? The parallel version gives a distributed vector. We cannot compare this one directly with the result of the sequential function. A possible solution is to gather all the pieces of the parallel result in a single rank and compare there against the parallel implementation.\n",
"\n",
"\n",
"The following function gather the distributed vector in rank 0."
]
},
{
"cell_type": "code",
"execution_count": 39,
"id": "ec9dcdec",
"metadata": {},
"outputs": [],
"source": [
"code5 = quote\n",
" function gather_final_result(u,comm)\n",
" load = length(u)-2\n",
" rank = MPI.Comm_rank(comm)\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",
" # If I am not rank 0\n",
" # I send my own data to rank 0\n",
" u_snd = view(u,2:(load+1))\n",
" MPI.Send(u_snd,comm,dest=0)\n",
" u_root = zeros(0) # This will nevel be used\n",
" else\n",
" u_all = zeros(n+2)\n",
" # If I am rank 0\n",
" nranks = MPI.Comm_size(comm)\n",
" n = load*nranks\n",
" u_root = zeros(n+2)\n",
" # Set boundary\n",
" u_all[1] = -1\n",
" u_all[end] = 1\n",
" u_root[1] = -1\n",
" u_root[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",
" u_root[lb:ub] = view(u,lb:ub)\n",
" # Receive and set data from 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",
" u_rcv = view(u_root,lb:ub)\n",
" MPI.Recv!(u_rcv,comm;source=other_rank)\n",
" end\n",
" end\n",
" return u_all\n",
" # NB onle the root (rank 0) will\n",
" # contain meaningfull data\n",
" return u_root\n",
" end\n",
"end;"
]
},
{
"cell_type": "markdown",
"id": "93fb1bd6",
"metadata": {},
"source": [
"Run the following cell to see the result. Is the result as expected?"
]
},
{
"cell_type": "code",
"execution_count": 97,
"id": "8e173a5f",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"u_root = [-1.0, -0.7999999999999999, -0.5999999999999999, -0.3999999999999999, -0.19999999999999996, 0.0, 0.19999999999999996, 0.3999999999999999, 0.5999999999999999, 0.7999999999999999, 1.0]\n",
"u_root = Float64[]\n",
"u_root = Float64[]\n"
]
}
],
"source": [
"code = quote\n",
" using MPI\n",
" MPI.Init()\n",
" $code1\n",
" $code2\n",
" $code3\n",
" $code4\n",
" $code5\n",
" n = 9\n",
" niters = 200\n",
" comm = MPI.Comm_dup(MPI.COMM_WORLD)\n",
" u = jacobi_mpi(n,niters,comm)\n",
" u_root = gather_final_result(u,comm)\n",
" @show u_root\n",
"end\n",
"run(`$(mpiexec()) -np 3 julia --project=. -e $code`);"
]
},
{
"cell_type": "markdown",
"id": "d0edbad5",
"metadata": {},
"source": [
"Now that we have collected the parallel vector into a single array, we can compare it against the sequential implementation. This is done in the following cells."
]
},
{
"cell_type": "code",
"execution_count": 98,
"id": "6abd3a76",
"metadata": {},
"outputs": [],
"source": [
"code6 = quote\n",
" function jacobi(n,niters)\n",
" u = zeros(n+2)\n",
" u[1] = -1\n",
@@ -633,53 +1043,122 @@
" end\n",
" u\n",
" end\n",
" 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`);"
]
},
{
"cell_type": "markdown",
"id": "eff25246",
"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 (excluding the messages to collect the data on rank 0)?\n",
"</div>\n",
"\n",
" a) 1\n",
" b) 2\n",
" c) 3\n",
" d) 4\n",
"\n"
"end;"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "98bd9b5e",
"execution_count": 102,
"id": "58c2f44e",
"metadata": {},
"outputs": [],
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Test passed 🥳\n"
]
}
],
"source": [
"code = quote\n",
" using MPI\n",
" MPI.Init()\n",
" $code1\n",
" $code2\n",
" $code3\n",
" $code4\n",
" $code5\n",
" $code6\n",
" n = 12\n",
" niters = 100\n",
" comm = MPI.Comm_dup(MPI.COMM_WORLD)\n",
" u = jacobi_mpi(n,niters,comm)\n",
" u_root = gather_final_result(u,comm)\n",
" rank = MPI.Comm_rank(comm)\n",
" if rank == 0\n",
" # Compare agains serial\n",
" u_seq = jacobi(n,niters)\n",
" if isapprox(u_root,u_seq)\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`);"
]
},
{
"cell_type": "markdown",
"id": "1bd7a876",
"metadata": {},
"source": [
"Note that we have used function `isapprox` to compare the results. This function checks if two values are the same within machine precision. Using `==` is generally discouraged when working with floating point numbers as they can be affected by rounding-off errors."
]
},
{
"cell_type": "markdown",
"id": "a6602ada",
"metadata": {},
"source": [
"<div class=\"alert alert-block alert-success\">\n",
"<b>Question:</b> What happens if we use `u_root == u_seq` to compare the parallel and the sequential result?\n",
"</div>\n",
"\n",
" a) The test will still pass.\n",
" c) The test will fail due to rounding-off errors.\n",
" b) The test might pass or fail depending on `n`.\n",
" d) The test might pass or fail depending on the number of MPI ranks."
]
},
{
"cell_type": "code",
"execution_count": 105,
"id": "6b1c9f68",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"It's not correct. Keep trying! 💪\n"
]
}
],
"source": [
"answer = \"x\" # replace x with a, b, c or d\n",
"jacobi_2_check(answer)"
"jacobitest_check(answer)"
]
},
{
"cell_type": "markdown",
"id": "0825df69",
"metadata": {},
"source": [
"Run cell below for an explanation of the correct answer."
]
},
{
"cell_type": "code",
"execution_count": 106,
"id": "7a1a1342",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The test will pass. The parallel implementation does exactly the same operations\n",
"in exactly the same order than the sequential one. Thus, the result should be\n",
"exactly the same. Note however this is often not true in other parallel algorithms.\n",
"Specially, when one does parallel reductions.\n",
"\n"
]
}
],
"source": [
"jacobitest_why()"
]
},
{
@@ -687,7 +1166,7 @@
"id": "c9aa2901",
"metadata": {},
"source": [
"### Latency hiding\n",
"## Latency hiding\n",
"\n",
"Can our implementation above be improved? Note that we only need communications to update the values at the boundary of the portion owned by each process. The other values (the one in green in the figure below) can be updated without communications. This provides the opportunity of overlapping the computation of the interior values (green cells in the figure) with the communication of the ghost values. This technique is called latency hiding, since we are hiding communication latency by overlapping it with computation that we need to do anyway.\n",
"\n",