![](\"https://upload.wikimedia.org/wikipedia/commons/thumb/3/39/VU_logo.png/800px-VU_logo.png?20161029201021\"){kind=logo}
\n",
"\n",
- "## Serial Algorithm\n",
- "To demonstrate the algorithm, we will consider a simple system of linear equations $Ax = b$:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 66,
- "id": "2b369b73",
- "metadata": {},
- "outputs": [],
- "source": [
- "A = [1.0 4.0 5.0 8.0 1.0; \n",
- " 2.0 -1.0 4.0 3.0 0.0; \n",
- " 7.0 6.0 3.0 -4.0 5.0; \n",
- " -3.0 4.0 2.0 2.0 2.0; \n",
- " 0.0 -4.0 2.0 1.0 2.0]\n",
+ "### Programming large-scale parallel systems\n",
"\n",
- "b = [61.0; 24.0; 37.0; 29.0; 12.0];"
+ "# Gaussian elimination"
]
},
{
"cell_type": "markdown",
- "id": "53124eb8",
+ "id": "ff0fbd76",
"metadata": {},
"source": [
- "The code in the following cell converts the general problem $Ax=b$ to the upper triangular equation system $Ux=y$. Note that this function assumes that the pivots are all nonzero. This function will be erroneos if any of the diagonal entries are zero!"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 67,
- "id": "7a7b926a",
- "metadata": {},
- "outputs": [
- {
- "data": {
- "text/plain": [
- "convert_to_upper_triangular! (generic function with 1 method)"
- ]
- },
- "execution_count": 67,
- "metadata": {},
- "output_type": "execute_result"
- }
- ],
- "source": [
- "function convert_to_upper_triangular!(A,b)\n",
- " n = size(A,1)\n",
- " # Upper Triangularization: convert Ax=b to Ux=y\n",
- " for k in 1:n\n",
- " for j in k+1:n\n",
- " # Divide by pivot\n",
- " A[k,j] = A[k,j] / A[k,k]\n",
- " end\n",
- " b[k] = b[k] / A[k,k]\n",
- " A[k,k] = 1\n",
- " # Substract lower rows\n",
- " for i in k+1:n \n",
- " for j in k+1:n\n",
- " A[i,j]=A[i,j] - A[i,k] * A[k,j]\n",
- " end\n",
- " b[i] = b[i] - A[i,k] * b[k]\n",
- " A[i,k] = 0\n",
- " end\n",
- " end\n",
- " return A, b #U,y\n",
- "end\n"
+ "## Contents\n",
+ "\n",
+ "In this notebook, we will learn\n",
+ "\n",
+ "- How to parallelize Gaussian elimination\n",
+ "- How the data partition can create (or solve) load imbalances"
]
},
{
"cell_type": "markdown",
- "id": "78ef1849",
+ "id": "8dcee319",
"metadata": {},
"source": [
- "The function in the following cell solves the upper triangular equation system using backwards substitution. Note that the function alters the input values. "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 68,
- "id": "5a134433",
- "metadata": {},
- "outputs": [
- {
- "data": {
- "text/plain": [
- "solve_upper_triangular! (generic function with 1 method)"
- ]
- },
- "execution_count": 68,
- "metadata": {},
- "output_type": "execute_result"
- }
- ],
- "source": [
- "function solve_upper_triangular!(U,y)\n",
- " n = size(U,1)\n",
- " for step in reverse(1:n)\n",
- " if U[step,step] == 0\n",
- " if y[step] != 0\n",
- " return \"No solution\"\n",
- " else\n",
- " return \"Infinity solutions\"\n",
- " end\n",
- " else\n",
- " # Backwards substitution\n",
- " y[step] = y[step] / U[step,step]\n",
- " end\n",
- " for row in reverse(1:step-1)\n",
- " y[row] -= U[row,step] * y[step]\n",
- " end\n",
- " end\n",
- " return y \n",
- "end"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 69,
- "id": "b92332f7",
- "metadata": {},
- "outputs": [
- {
- "data": {
- "text/plain": [
- "5-element Vector{Float64}:\n",
- " 1.0000000000000009\n",
- " 1.999999999999999\n",
- " 2.9999999999999964\n",
- " 4.000000000000002\n",
- " 5.000000000000005"
- ]
- },
- "execution_count": 69,
- "metadata": {},
- "output_type": "execute_result"
- }
- ],
- "source": [
- "U,y = convert_to_upper_triangular!(A,b)\n",
- "sol = solve_upper_triangular!(U,y)"
+ "## Gaussian elimination\n",
+ "\n",
+ "[Gaussian elimination](https://en.wikipedia.org/wiki/Gaussian_elimination) is provably one of the first algorithms you have learned to solve linear equations like this one:\n",
+ "\n",
+ "$$\n",
+ "\\left[\n",
+ "\\begin{matrix}\n",
+ "1 & 3 & 1 \\\\\n",
+ "1 & 2 & -1 \\\\\n",
+ "3 & 11 & 5 \\\\\n",
+ "\\end{matrix}\n",
+ "\\right]\n",
+ "\\left[\n",
+ "\\begin{matrix}\n",
+ "x \\\\\n",
+ "y \\\\\n",
+ "z \\\\\n",
+ "\\end{matrix}\n",
+ "\\right]\n",
+ "=\n",
+ "\\left[\n",
+ "\\begin{matrix}\n",
+ "9 \\\\\n",
+ "1 \\\\\n",
+ "35 \\\\\n",
+ "\\end{matrix}\n",
+ "\\right]\n",
+ "$$\n",
+ "\n",
+ "\n",
+ "The value of $x$, $y$, and $z$ can be found by creating an *augmented matrix* and applying Gaussian elimination to it:\n",
+ "\n",
+ "$$\n",
+ "\\left[\n",
+ "\\begin{matrix}\n",
+ "1 & 3 & 1 & 9 \\\\\n",
+ "1 & 2 & -1 & 1 \\\\\n",
+ "3 & 11 & 5 & 35 \\\\\n",
+ "\\end{matrix}\n",
+ "\\right]\n",
+ "\\rightarrow\n",
+ "\\left[\n",
+ "\\begin{matrix}\n",
+ "1 & 3 & 1 & 9 \\\\\n",
+ "0 & -1 & -2 & -8 \\\\\n",
+ "0 & 2 & 2 & 8 \\\\\n",
+ "\\end{matrix}\n",
+ "\\right]\n",
+ "\\rightarrow\n",
+ "\\left[\n",
+ "\\begin{matrix}\n",
+ "1 & 3 & 1 & 9 \\\\\n",
+ "0 & 1 & 2 & 8 \\\\\n",
+ "0 & 0 & 1 & 4 \\\\\n",
+ "\\end{matrix}\n",
+ "\\right]\n",
+ "$$\n",
+ "\n",
+ "The result is an upper diagonal matrix with ones on the diagonal. From this matrix the values $x$, $y$, and $z$ can be found via [backward substitution](https://en.wikipedia.org/wiki/Triangular_matrix).\n",
+ "\n"
]
},
{
"cell_type": "markdown",
- "id": "962ec4a9",
+ "id": "94c10106",
"metadata": {},
"source": [
- "We can test if the obtained solution is correct using `@test`:"
+ "### Serial implementation\n",
+ "\n",
+ "Gaussian elimination can be implemented as shown in the following function. Note that the result is overwritten in-place on the input matrix."
]
},
{
"cell_type": "code",
"execution_count": 70,
- "id": "0e336c85",
+ "id": "e4070214",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
- "\u001b[32m\u001b[1mTest Passed\u001b[22m\u001b[39m"
+ "gaussian_elimination! (generic function with 2 methods)"
]
},
"execution_count": 70,
@@ -176,22 +120,261 @@
}
],
"source": [
- "using Test\n",
- "@test sol ≈ [1.0; 2.0; 3.0; 4.0; 5.0]"
+ "function gaussian_elimination!(B)\n",
+ " n,m = size(B)\n",
+ " @inbounds for k in 1:n\n",
+ " for t in (k+1):m\n",
+ " B[k,t] = B[k,t]/B[k,k]\n",
+ " end\n",
+ " B[k,k] = 1\n",
+ " for i in (k+1):n \n",
+ " for j in (k+1):m\n",
+ " B[i,j] = B[i,j] - B[i,k]*B[k,j]\n",
+ " end\n",
+ " B[i,k] = 0\n",
+ " end\n",
+ " end\n",
+ " B\n",
+ "end"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "cb7f89e5",
+ "metadata": {},
+ "source": [
+ "Let us test the function with the example above:"
]
},
{
"cell_type": "code",
- "execution_count": null,
- "id": "28dff449",
+ "execution_count": 68,
+ "id": "eb30df0d",
"metadata": {},
- "outputs": [],
- "source": []
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "3×4 Matrix{Float64}:\n",
+ " 1.0 3.0 1.0 9.0\n",
+ " 1.0 2.0 -1.0 1.0\n",
+ " 3.0 11.0 5.0 35.0"
+ ]
+ },
+ "execution_count": 68,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "A = Float64[1 3 1; 1 2 -1; 3 11 5]\n",
+ "b = Float64[9,1,35]\n",
+ "B = [A b]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 71,
+ "id": "52bfada2",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "3×4 Matrix{Float64}:\n",
+ " 1.0 3.0 1.0 9.0\n",
+ " 0.0 1.0 2.0 8.0\n",
+ " 0.0 0.0 1.0 4.0"
+ ]
+ },
+ "execution_count": 71,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "gaussian_elimination!(B)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "563aa0d0",
+ "metadata": {},
+ "source": [
+ "We get the same result as shown before (as expected)."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "39f2e8ef",
+ "metadata": {},
+ "source": [
+ "## Parallelization\n",
+ "\n",
+ "Gaussian elimination is expensive and thus it makes sense to try to seed-up this computation by using several processors. It requires $(2n^3)/3$ operation, where $n$ is the number of rows in the input matrix. Thus, the time complexity is $O(n^3)$, which rapidly grows with $n$.\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "1b1a6469",
+ "metadata": {},
+ "source": [
+ "### Where can we extract parallelism?\n",
+ "\n",
+ "```julia\n",
+ "n,m = size(B)\n",
+ "for k in 1:n\n",
+ " for t in (k+1):m\n",
+ " B[k,t] = B[k,t]/B[k,k]\n",
+ " end\n",
+ " B[k,k] = 1\n",
+ " for i in (k+1):n \n",
+ " for j in (k+1):m\n",
+ " B[i,j] = B[i,j] - B[i,k]*B[k,j]\n",
+ " end\n",
+ " B[i,k] = 0\n",
+ " end\n",
+ "end\n",
+ "```\n",
+ "\n",
+ "- The loop over k cannot be parallelized (the state of B at iteration k depends on the state at iteration k-1)\n",
+ "- Loops over t, i, and j can be parallelized\n",
+ "- BUT the loop over t needs to be done before loop over i and j"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "a67e0aad",
+ "metadata": {},
+ "source": [
+ "### Data partition\n",
+ "\n",
+ "Let us start considering a simple 1D block partition. We assume that each process contains only a portion of the input matrix consisting in a block of consecutive rows. In this algorithm, the data stored in a process does not correspond with the data updated at a given iteration over the outer loop over k. The data updated is only a subset of the data stored, which leads to load imbalances (as we will see in a second).\n",
+ "\n",
+ "Consider next figure. Let's find out which is the data updated and the data used by CPU 3 at iteration $k$.\n"
+ ]
+ },
+ {
+ "attachments": {
+ "g20491.png": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAABS8AAAJNCAYAAADQ04CWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAewgAAHsIBbtB1PgAAABl0RVh0U29mdHdhcmUAd3d3Lmlua3NjYXBlLm9yZ5vuPBoAACAASURBVHic7N13uCRlmbDx+0xgYAYYcnLIGQQERIwY1xyWNec1YN5VV1fd/VDXtLrmHBFRMeKacA0YMICCogSVICBIziBxcn9/PKc5VW9XV1f16e6qOdy/6+oLuqfCc6qrqt963gSSNH4fBjqZ1yOaDecu6WZmjv8VE9zvEuBs4Ibp14sGLP8n8ufKwrFGp3Hakvx3+atmw1mnPIn8sXv3GPbxVWauy8+NYfuSpPbaADg489plwvs/ILPvrUqWm08+zr3HH5rGaBPy3+eyZsNZpywjf+w2GfH210+2v3jE25+1BU0HILYEPlZhuTuI5MMtwI3AX4DTgUvGF5okzdqRzBQ0LwI+32Askma8C3gyMA94LnAM8IsmA5IkTcyuwGmZ9/9LVJxNwsOBH03//0pgH+CaPssuIR/nacAh4wtNY3YY8J3M+/cAr2solnXNq4DXZN4/EfjmCLe/EvgMcOD0+/cn+2ucycvmLSYeHoZ1HXAKceJ+g0huttVGwP/LvP8b8ImGYpEADgcOzbz/LHB+Q7HMRXsA/5Z5/yZgxYj38d9E8gXgJiIho9F4C7Bo+v9vA97WYCwavTOArwFPB6aAjxAF1jVNBiXVdCDw0wHLLGemAcBNRCOA84A/TL/+Ns4AJeWsT77hzkeBCxuKRdKMtUSu5vvT7/+V6Jnzp8YiUuvsSL5r2GxeNxEPl0sm+hdUty35eE9uNhxNUFu7jX+SfFwPbzacsWqi2/hXM/u8kGoVZnW7ja/OLOsD6Ghlz5lrh1jfbuPDm0S3cYD9icJqdz/PGdN+pHE5hNmXny8nyil3n3Dsde0C/Djzekuz4Wgdd3fy18E3JrTfN2T2eQuw+YDlNyYf5+8GLL8n+evk/5Uvrhp2In9sh6nUfjyTKd/MRe8lf+z+aUz7OSmzjx+PaR9DseVlO50FrEo+mwcsBTad/u+8dKXpz48kHj6eBvxmjDFKUpl7kG9V/h4i0TjIZcQ4TF2dUQaliVoD/DXzfpJjraqas4D/Ax47/f7NRKXDysYikiZvO+BfgFcA3yV6DPy1dI1mbAQ8LPP+1qYCkYa0Cfkuwp8Crh/xPpaSv05Gvf27sg3JH9s7mgpEY/VO4HvT//8w4CHAz5oLZ4bJy3Z6NFEL3M8SYhDVewPPJFpOZO0A/JxIHHx3DPFJ0iD/wUwly9+BL1Zc75HjCUcNuIEYU0vt9jFmkpe7EJWfX2guHGnWjiJaFKeWTr+2J8ZiThsCTAFPIB7WjgC+MsYYpbui1xINcSAa6nygwjqrgOMy7+1ivm67gvz3eVZTgaiv7xM94bq9Ed4B3Ke5cNQmRd3G71ZzGw8kxq5Kt7Ocdg1obLfxuy67jTdvkt3GdyIKm939fXSM+7Lb+PjMttu4hjepbuMQCZyLMvs6fYz7kkatqNv4oOFGILqiPoyoWFtesI21wMvGEO9sHEA+xm81G47WcZPuNr4RMd5sd3/HlS8+tHuR/7u+Oqb93BWl54yNpCZrUt3GIX7/svt64Bj3VVlR12Ote35BtML8XPL5IuBLzEy4IEmT8DLyLfuPaSgOSYOtJd8y+h7EbKDSXHYz8BPg2USF23eSf+9OYvUwJI3C84lu412fbioQSQN9Cbg98/7VTQWSZbfxuWM58EJga6LbedfuRCKhSrP8rG2AfYhWoBsRY5fdAFxKtMoY9YzBs7ElsC+wjIi1Q8R6GRHrXXk8jl2IGZ+3IwriVxEDXV8zou1vBxxEtKpdRbToO4/Rt4BbjxiAe1dgK6JVxS3E93wG8V23zTwiCbCMOEcXAlcDFzPTUnq2pojjvwuwBXE8Lp/eflNjYc0HnpV5/zfg9w3FMlsbE8d3G6Kb0+1EC8TfE9/lqC0khgG5G3FtrQLOZXAr9c2ILpA7EeMRzSNaN1w5Heu6Pi7afOCexDHZhvj9uY4Yk25cMyBuSFy/2xMPW38H/kLMTFzUHXVd9w3gjZn3zwV+2VAs0qRdBRxOTOyRnYBiHnA0sBf5hzhJ9WVbMl8J/LSpQCQN9HfgeOCp0+8fR5SJL20sIkxezjVriaTB+eRnbnsl0W13Tcm6i4kxrx4LPIg4OftZTgza+h5ibM1BTiMe/Ocnnx9E+bglPyPGHEqtDzyKmK3sgcDOJdtYQTyAvR/4YYVYZ+vLwKGZ9/ehXqLwHCJRB5GQK2ui/Rjie+36JPGdQIxZ9p/AfgXrrSEKDEcyeMa+fg4jZrk8jOIW3L8nHgDSlgx17Et0nXwo0QWkrAXxRcQDxseIpE2ZZwP/Nf3/WyT/dgzlye5/BP44YPt7ETMpPppIWha5mhhL678ZrjvuQmJCgZcSQ0+kVgBfB95EJEsn6aFEMrvrW9RL1P6QqHTp2pPeiX6+RyTrIH9f2Y7ye8o5zIzt18884vp5EXBfirsfdogk1geJa75qMusUZs6J5cQ5zvRn/w/4Z2I8tqzfE4m7rIVEa6B/JO4Re5bsczUxedsHgW9XiPU44t4MMb5y16aUH9uL6G2htBn5e8xpzBSCqtiDuE89lpkxslKXEOf6O4nkfRXbkE8I/474ziHOoXcCTyF+a1LXEJWBHyS+wzZ6CfDvyWf/DXy2ZJ2zgAuA3abfP4mYvOSuXPmnu5YO8Hbi9+c5mc+3J66FOkM37ETct/ciKjA3JH6XbyDuWb8iyhJl5fKs7v1v4+Tz9eh/b2R6n4OSrsuIMtZexN+6IVFx1m0EcBJRIVplwr1R2Jj87/qgMl3WFPmWfauJiu4q5hNl9oOI++BGxLPRTcBtRFnqXOL3ok5MWfOIOQseQPzWbElMjnY1cDZwAlExN1v7Ec9Jy4jjcSWRcDgR+PMItj+MQ4nf9K5vUb3slH6va4iW0/2W2Sj5t1FcJ10bE+Xc/Yjvb1PiWrmauK5/Tb1rpTshb1d6zs4H7k+cl8uIRkrdpFLZc+2O5O9BGxF/543M3IPOYvb3oIWUH9uVxPWTrrNh5v1y6pc15hO9Th9AlOm2II7LdcS99cfT7+tYj3y5N41rHlHOPYRoZLBken9nAD9idA2DxiU97lD8/WR9g5ly+zxirpV3jT40rUtGMeZl6m0F2yzrArY+caNM16nyOpbiB7ysa4fc9vF9tnfdkNv7Fr0/aKN2YrLP7Wqunx136aIByz412dc7iQTf16h2PFYzk8Sragp4H1HgqLKPzxKVJHXHvHxjxe2nr2uJ5HuZVwy57Q69SaSs9YguZqsK1uv3+jv5Gbmr2J5oUVxl+7cwk6yb1JiXH0tiGJQsTP0pWb8oeVj1709fg8byuzuRlKyzzVPJJ2vLXJlZr1tQfiDl97Q/FGzn7CH//p/Rm7BPnTzkts8r2NaWyTK/GrDvrilicPCVNfZ/A/D0ittflqzbbWH4D8QDSJX9/YbyAvtsDTPm5RTxO5BdbxWRiK/iE8m6jy5fXGqFYce87Gcx8RCa3d7FlA+1NY8oZ39setkq95C/Af/KTIV12baHuSd3KK6wmCKSdB8kKiyqbOcKokJkgwGxjsKfk32nDR/KbJKse2qFdRYSs19fRrVjsXZ6u0+pEdcSooJy0O/LaqICcffizQx0CFFJWLaPE5lp2DDJMS9nM/b9xsm6RQ0vNqDa91f0+niFGPYBvsngcsl1xLUy6Nm4a+dk/W5r1CmiEuWvffaTjiU/j0jkfZT8GNZlr0uIBk5Vhper+tyXvo4p2Nbjk2XqVAxtQEwIOigXsJJoJFLWwCn15GQb/5P5txcSx6vf/lYQ99Rx5xmGHfNyK6IBRXbdvzHTiKKfJcTf1l1nUAMe3QXsSO8FMNvk5bb03mTeWbL84oIYOkQm/m9ETePF9L9h/5jyQt2ok5dFcdxO3FTOJW7aKwqW6RAPnIMKirNxYrK/SScv08TRbURB8I9E7XHRMXlbuuESH+mzjeuJGrxziFqq7L+9i/qFlncX7GMt8UBxPtF9s98P10rgISXbHkfycilR4Cha53ai8HHe9P+n/74GePGA49G17fTfXnRs/kYUWi8lf/3fQXR/nVTyMlvQWk1vS8JBmkpePpj8QPLZ19XE93cZ+QmCstdqlQJSmrw8mOjSXTfmSwuWWz4d37lEC8miCSg6xDWa1rxmNZ28nE/MdN3vWrpg+u8vqiRYSyQDBilKXt6f3t+NbuujfhUS42zNXzd5uR5RmZhd51bqVR48JVl/nBNtSaMy6uQlRMVuus2ymVbTBESd16+I1lT9jDp5udkstvd7yntmjcIkk5ebTS8zzLGoOhP9vYnhfOpsewXVK+O6Xk5x+aTfb+ljmWzyMls2XEG+hdsgTSYvp4jnpKrHtvs6iyhrDFKUvJwPfH7A9j+WbGeHWfz9JxOtF8u0IXm5B1G+rbP/5eRb0pcpSl4uInpYVd3fKeRbCY/aMMnL3emtqDqd6jmKXybr1kkIaw7akd4Tf7bJS4gH1Ow2f1KybDd5eSvRau8ZxNiCU8ly6xPNpb9TEPN/lGx/R2JMvnsn6/x++vN+r3430pXED+83iRvSHvQmT9cjWjQVtUIcZ3PnE5N9TTJ5+cfM/19KdIPM1vzNBx5Jb3KoQ7UH3GcXrHcq0dIxe/yXAM9jpuXC2oJ9Vk1e/pF4iLgfxbVZOxI1nGkriSuIQmmRjZk5x9IfpOdQfk4W1U5OAd9PtrMK+AzxUJUteC8gkmTpebKS/HADRaaI6zi73hqiC+sOybK7EIWxbmHjfPJJ5XElL3dK4hume1KV5OUyZr6TbIHycsq/v34Fyb3oTSJeQiTC0vvx5sRDQtqC4tQ+sWZlk5fLybeg/AHwRCJBvYBI/D2AaBGSupQo/P8fURu8L73DwCwgzqnPEudJNtajS2LcjpnjlT0mN1B+bIseaIdJXv5nsk4HOJO4R2WP71JiWJH0e1jL4Ak20uTlX4gx7zpEa+W3Ew91XQuI+1xagOtQr/VNHXWSl5sQrWqzy19N3H/q2CnZxrjGE5VGaRzJy+0LtvmmkuWLkpdXEL/ZXycmJjmOeGAsqgw5jf4V6/OIB/YL6a24ujXzb0WvovJuUfLyGuIe0o3169MxFTUWOId6iae6JpW8LCpTdX/rTiSSk58gjse3iIrBbHmjSvLysUQjgnQf1xLlxi9Mb6c7lnL6W/a8Kn/09HJFyaUzp/fxSaJBSLaX3c1E4iO7/LiSl3sn+zmp5vpVkpfrM3Pep61ob6H8Onlrn/3OJyazS49r97nmf4ny1HfIl++6r78xuGdOUfLyfclnK4nz7zSibLqa3oRrUfKyO65o9rrudw/6A+UtMC+g+B50G+XH9j0F2xomebkvxS2XbyW6bH+e+B7+VrDMWqLhyiBFycujks8uI4bM+z5xfRVdd8dW2New6iYv70NvI7IT6O3+X+atyfovK19cc92O9J70o0heHp1ss2xw1UVEoaxfsqfIEeQfiG8gkqBltk1iGjQJRT9vJZo/V/U08jfq26j3t9ZxIvm/cZLJy+7rT5T/fevTW2C7hPLC6Ob03vy+Tvm4udvRv7vDoOTlU4gkX1VbEcnw7D7eUGG9TybrPLzGPrtem2zjeiLpVGaKSDpm1zuH8hbMz0mWX01MMFDmmfQmrjqML3n5xGQ/VVsmZFVJXmZlHyaGmShqEVEAye7zBwyuOd2e3kqidJzBVFHhdjUxOUod/0m91i8PJ9/qdzXVak6zrXWHGZu1bvJyP3oflr9J+TmwDVGoTwuXZa1L0+Rl93UhM2M+FllIJIyz65w44G8aVtXk5fbkK646RCvYXYbY5xT5FvprGH8XKGm2xpG8hN6H4LIxvLsJiJOIB7uyFlfbEGPQpve6Kr1gDkjW+VaFdVLd5OVviS6jO5UsuwXxfJD2qPnIEPutalLJy0cny95IlLPK9rcp0SLyR0Tld5k96R2S63SiEUHaOATiHPpmsvzt5CvSiuxKb4L0jxRXiC8mnqG6z0RpuX5cycsjkv18uHzxHlWSl1n3Spb/as39db052c5aIqG9Y8GyU8Q45GnX4kG9E9Pk5VXMlNsvI2ZoTxNN2xGNOrK6yctfExXsZWXErYlK2vQeVNZLsyttrfvdCuuk6iYv1ydasmbXWU40nEqfXaeIa/viZPmVDK7QTZOX2W38gOLedzvR24ClQzTYGoc6ycvD6e3xdwz1fx/T72uYZzvNITvSe8KPInmZtl4ZxyyJH0z28c8Dlh9V8nIYb0n2/aox7efEZD+TTl7eQbWH1qX0JlLKxkV7XbLseVQbz2V/imv46ox1U9Uy8i3FBh0/mH3ycmPyD/urKR9fNmseUfuV3f8TSpZPxzl8S8X9FHXBH1fy8u3Jft48xDYmnbx8frK/31N9XK+9yRcMLqd8WIqi5OWghOeovCzZ739XWGfSycu0m9SfqXaf2ZXeAtrLS5YvSl7eyswkUGWWkS/wr2U8Y19WSV4eQG8Lk18zeFzTMr9Otnf/WWxLmoRxJS/TnjtnlSy7Cb3JhEEeTn6oihsY3KJxFMnLDSgfWqfI/ciXr24nPzHoKE0qefnZZNmy8le/ffUzj97hbb7B4N+zKXqHWSrrOQe9iZOzGDxcz3MpbjE2ruTlp5L9FE3GWqaJ5OV9yJcvVxG/y4NsR28X3aeVLN9vyIkz6D/pZ5Gl1P+9fij5iombGFxh2UTy8r+S5VcyeEzurYleZ9n1/kR5IjlNXnZf76O4wqFrIb09cz4zIL5hVU1e/gu9Qx28jfK/o5+dku2UTd6pu4Ad6b1IRpG8fGnBdkc90HbaRL2sKyI0m7zcjHwS7Ztj2s+J5P/GSScv6wx6/KJk3dP6LDdF7w/AoBZ/WZ8uiHMcyUvoLYzuOGD52SYv/y1Zf9A1kDo4Wf/bfZa7f7LclVS/njciWoNm1x9X8jLthv+CIbYxyeTlVMH+7ltzGx9P1i8bgiFNXp5PvQez2ViPmCCqu+8q3bYmmbzcnN5xOutcj+lEdWeXLFuUvPyfkuVTP0jWrZsIqGJQ8vIfyH+fHSKRMdvf+a8n26w6VpTUlHElL9Oum5eNYJupdyX7eOaA5UeRvBzWG5J9j6vr4KSSl9lkw62UJzXqOjyJ40yqT+CygN6eRPfos+ye5JOQq4iJZaooGlNxXMnLtFLsgTXXbyJ5+b1kG3Uqmu9F/nv5bcmyRcnLWxmu98Qw0kYHg37zJ528XJ/eocH+q+J+DqQ3gfeokuWLkpc/odq94dBkvUsqxljXoOTlPKKrfnaZVdSvMMiaIp/kXkuDvXJGeaNWu9xY8Nmou0pfQrQ06qo7vtYk3UCMa9bV5lhn45gay36FuBl1HUxxd/xdyXelvI74Ua+qbkJvNk5J3o/7e35i8v6TNdf/PZHA6rofxbVij0zef438d1fmFiIhMQk7Ju8vL1yqPfYmP9PeH4lCdh1fS97Xqfk+hugeNAkrifOt6yAmlzit4jDy4y1dRHS3qurTRIGqa2/qda2vc59KH4iHnRl2WM8luq9nu5J9lEh4Vr0v9JNWbOw0y+1J66q0HD2O4YaOSt6XTQrUtKPJ32PbHGsV2Vau8xjt7+Grk/dvICrnqlhNb/fdJ/dZ9nnky4zHUl5xl3Uk+e9znNJhatpeNtyTfMu+i4EP1Vj/t8TQAl2HUG+Sk08Rw25NQtvvQf9IvgXqtVSfu+J04EvJZ2W9DItUvU5OJd+AYntm1wtmGN3JhV6b+ew24hjOpiVoh3zZcIrBjYPGpmy8Oq3bin6E6z4kzyNaqGxM/MgXZdlXZv6/bivDUZoiH2vRQLTZm8820+t0xh/axFxK9UILRFLrZPKTW9yL3sRkmgD8GVGLU9WpxEPAKLpWLiV+xJYSY9qlrSvSVsuDBsqejQ3Ij39yI4Nrg4ucxkzyYwsiWXBRskz6HZxQcx8/Al5SO7L60u/4lgnsczbSsUl/VLhUubTFcp2E+S+H2F+ZuxHXxmKKu7Rlf/M3mF7m+hHHMKx0fKBuV7iqLiVatxyYbLNsvOeu68hXblXZV9Y4Z5ZMvZEYMqL7wNoBXk/xoPjDuDl5P44u8dK64Kbk/QZEK6CqSagqLiC6YHfHjD9ohNsetWuICTO65ao2x1rF1Zn/34DozTSKiTaWku/BcRX1yxY/IMrZ3TJuv0rRByXv00RNmUuJiUfG0XMgaz2iC2/WlWPe52w9inxS+Mvkn3erOJ58w4P7UW04KxjvhC+pi4myevcZv23XdToU11epdw/+HPnWpIdR/fn/YnobxZQ5nXxSbweifDkJmxK997LH62qiN1i/npV1XEm+NfAONDSpo8nLuavogeOGCuvtQQxG/WiiafigSXiyljLZhODORKyPJSZ6KJugIbWAuFGnD2rrsjOHXCebvNyX3uTlvsn7uvvpEGPw1O0mAlGL9DhiAp97Ub+mZ5xJhf3Jj294JfGj3z23qkq7eW5FbwEn7QZU9zsY5twYRnq/mG0rsHFLB99eSbRAXkL52JWpFcy0GqwzRtEws7FnbUeMpfQEojvhoHGuUpvSnuRlOt7k6UNs43Tyycuq3ecup97v1q3J+zq/PcNaSLSQyA7FsIJoeTPKwdPT8bHHOauw1Gbpb/Ma4pqrajFxD9qZuNcWVbh2t9s16ZY6XYuIst7OxBAe/X4Ds8+NTcU6Kr8g34X0M8RzzweJhOOw7ke+AcnJ1G/heBtRDtxj+v1+BcusR747+RrqD8f1C8afvNyIfCKww3jmYRiltGJ70HjdRf6YvN+/4nq3Uz6+bh0bEPegXSi/B2XPz7Zd12nF9s9qrn8SkezsDtuwGdGy9twK65Z19y9ydfK+bpl8WDsRFf7ZcvR5RC5nVC1403JvndnKR8rk5dyVJi9vpbzWaCkx5tcRDD+cwHzi5jDuhMVi4B3AK5jdObyEuZW8HKYmM+0iWNQtKj2XRrGfKh5JdIXcdYh1u8b54J0mqfZhNLVbg76DDvW/g3GNcZlK7x2T6pI0rLSQ9p/Tr9mo2rWwQ2/LnqoWEl1ZXkf1cbSKtCkxld5nhulWlraIrNpqsG4L4fS8HmYA9LpeQv67vokYV+3nI95P+rc5vJDuqtL7x40MruRYj5i88hlEa7m6XZEn2Yp7ATHG5jOBB1O/PD3JWMfhKKLVevd7Xn/6/WuJIVZOIMax/y29D+5l0tnBt57ebl3Z5PkmxLmUTXQvIz/UyvnUf/5KE2zjkFZqdyeparM0WfwPRAVxHdsk76uWDf/K7IYTWkgMLfNMIgnb5ntQFWl3+7qNMVYTDQUOzny2C9WSl3XHem+i8vdA4lk529Pw18S4oqNsnJDeWxp7fjB5OXelreXKmi1vRvxIH9zn36+dfhX9KO5NvdaZs7UR0TKw34zO1xOx3lbwb3vQ4ACzEzBMIjZdp+hHK605GqYrcN3YnkfUghf96K4gasW7k9BkbcLskp11jKs75aLk/ULyPxK3UT8puIJ8zeO4pD/co54kbNTGMYZZ+v31s5zhkruLgOOIFslFbiJqf4setrq1722U3ntGcT+r+re2/UEKeq/da6jX1b2q9Pe86LdUuitIK7cG9V46BPgi0apnWOP+je7aj+iaWrU1WJH1iMqNtldS9nM90avnW+Rbz88nevrci6gkXE1UTP+Q6D58PuXScsX9qT8LdGqK3mFe0t/MYRIVw0zEV1cTlX2zlX6H/zaCbVYtj/x9Fvs4mLgHpT1Z6pjUPaiKefS28BvmnE1bRFb9Luq2EG6iLHlk8v54YgiMUTcka03jFJOXc1fazLqse+KHyScu7yDGiDiOGMOv7OHlDOrXRs3GO8knLlcCXyAmzTiV8sTaL+if9JwLhrme03VWFyyTfjbMoOZ11tmHGKw6u84ZRM3SL4gxovo5nPHNJJ9Kb+RXUG/M0X6uSd6vIX4kuvsb9r49ift9mjxqU8u+Iul3eBrDt4bsmk3Bs4rXk09criXGAPoSUdtaFv9xxKQubZTeZ0ZxP6szNm/bfZcY5qA7tvQexJipDyU/SPxspdfsXOqdINVxYPI+fQDOOpRoBJA+aK8kxgW7jJmuyOlEQC9nMkNPdO1HtCjcPPm820LpEiLWbu+A7AP5CwvWW5f9hHj+eS8xBFVRYm0B8Ux1b+DNwHeImaf7lUXHUSnajSMrbYxRp3XobNapK32GXES7k97dRPGoVS3T1BmaIusQYpLDtMHJKuIedCn970Evo52NezYiX07vMFyFapobqNqde12o2O6Qv2/tTtyDRj0pVmsqtk1ezk270zs2YL8BZ/churd03UIMAP2Hivua1HgOEDN3ZScdWU405T+p4vqTjHVYs5ntcJjxJ9JjUjRLfZoMGcV+yryZ/JgsnycKzEWJ1dQkx+BIa7lPJwq/o7aWuC67x3B94vjUScwsYTL3+0vJT1gzzgmTRiH9Dt9BDHjdVhsRs5Z2rSVmETy+4vqNjVFTQXrvGdf9bF11HtH646fM/L7vykwCs6xSp470mr1kRNuV1iWL6a2YP6PPsguJCv/sPevPwJuIiVcGtYB5HpNLXs4DPks+AXkBMRHY8Qx+ID2cuZW8hGjB/niie+pTiDEg70vxdzJF/OY+BPgn4n6cSofoOoHhxnBOpd9Nel4N02JuEq3sbiXK79ky6EaMv6J3WB2ifJ19Hnsvs+vKDXDOLNcv070HZctAZxP3oO8z+B70bNqZvExbPk4Rye+6k6alvcDaPuZqHUcSyefuhLV7EQ19Rl2xnZbJGytfm7ycm15Z8Nn/9Vn2CeQz9v9D9cTlfOpNTjFbjyP/Y/JRqicuoXf8kXFIf9zqXGMb1Vw+tcvgRQauk7b6g94m+sPsp2pX7kXkB1D/O9EqoUriEibzHXelx2W3Me8rWyjZhUhmVDXMdzaMi5P36ezvbZOe7+P8DkfhEeQLYcdSPXEJk70+6kqHNhnmu0jXmUSXuEm6kOg98FNm/tYdiILqwxjNw9Gy5P3FI9imtK55JL0TW/Qrbz6efDfNM4lJW6q2TJnkGHMPIV/BeAHRortqIqlt4+Gl6ky0l7qIeAb6H6IsfjBxJwIqbAAAIABJREFUX30YkczMbntj4H+JIQLSFrlppeipRBJp1NKGBcM00JjEMDKriVZg2UY1d6O9yUuI7zBbfv0I7a7Iewz54eL+BNyH6i1r23pdryL+hmxFwibUn1ArbQ09lyq2z2WmXLjT9GfjqNjeLnl/8Yi2W5sDsc89OxCDhWedRf9av3Rg6R/U2Ne+TLZbaBrr92usuyMxaPa4pV3s6tRk7TF4kVIHUL/lZtotqihxnZ47B9XcxyKqz/q7I/lj9gvqNU0/tMays3Um+dq/Peh98B+V2X4HdZcfVjpD4u4T2u+w0pkEH9pIFNXN5h64IdWvwybM9hyH3tnjR9HapW0uIQqq2SEqtiMm7hnFEC7p79CoZj2V1iXpGHcrgB/1WfYxyfs3U73csi2zS7jVlcb6NqonkTZkMq0u05aLVceRhtHFt5pIOr6DmMxoVyKBle3qvJSoXE+lEyTuNaKYUlcl8exB/WeAccWWSmc8HldZeVTSSTEndZyGVXQPqpq43JLJzl1RV5qoHGY8z7TsO6lJTCflr0S5MDse7w5EAnMU5f4F5Bs/rKHBZL7Jy7llEfANehOK7ylZJ205Waelyj/VWBZ6m63X7RY4yViHlQ7oXichOdvEySbAA2ssvwP5sU5XUZy8/B35cT8eSL1aukdRvWtK+h2XTTSV2pioIa8j7XpQJ9m8nBhjsGsKeHrN/VeVJtmeUHP9w0cVyAC/T97PZjKAqrL3lbr3lLTL14Npd1f32Vwfj6d+y+7ssR13l6J0aJPHUu+h9Z7kWxivpfe6mSuuJIZ3yXZj3Qr4GflWVXVtSf78v4Tycf6kuegJRMvJrK/TP8mX9myo0yPoPjWWhd4x2Oo+x6Wxnlxj3UOZ3dBGVVWZSLKfcZU5LgP+FXhr8nlRuf1XBcuMo6fjzeQnbVtM72Stg8zm96KOdMiF2UwoU8Vsr5P0O3zELGKZhNlc15O+B9WVluPSOT0G2ZN8y8vlzM1K2UuJ5/PsHCfbMpqK7T3I38P+zPBjs86aycu5Y0tirLb0h+g3xEQO/aTJm6rdCpeSH3+yilvI3/TqdjlPk59VY11MFDom4U/J+wdXXG894KUj2P/Lai6bvQd8n+JxQK4lWkB2LQJeUGM/r6ix7LDfMUQNeN0ES1pIrntOfjF5/++Mp/vFN8nXsD+eGAO2il3Jd8Ufp3PIJ/D3ZfytSrIDcS+lt6tfmYvIF/IWMZ7uXaMy7PWxAHjtEPvLHttFjHfc4JPJJ8o2BZ5bY/20pdRPaHe3tNm6lugCmi3Yb0YM2J8mXqo6OHlfJwkjzQW7Acckn60muhL3k3ZJrNNb5Nk1li3adt1uv7OJ9Vk19zWsdDiXOgm5tAXaqP1v8r6oBeE55Ft2bcb4GlD8Jnn/5BrrLiY/+d84nZa832/M+5vtdZJWbD+DdrdOXJfuQXWl5/jTaq7/jOT972kw8TZmVxJ5h2yvoy2JCdpmU1GRVgr9bhbbmjWTl+u+RcQD3unEGD1Z1xMXbdlsWRcn759YYZ9TwMep3w17DfnBY7egXhePi5P3VWIFeB8z40CM26nJ+2dQbebBdzGaGJ9ItRrCfYFXJZ99qmT5zyTvjyRabg7yNOq1KL2E3laeW1RY7wCGSzpdmLyv2zXkWPLdYbYEvsJwtezz6H89/JV8YWoD4IMUz4yZbvPDTK5b2hpicPquDYB7jXmf2e9wiqhlrSNtSfFi6heOujZn8HcyG+ng21XvgUfSO0REFen1UffY1rES+HTy2dup1r3sUfR+Zx8ZRVAtdyMxaV22lchSonvrQ4bY3oOS93WGkZHWdQ8jEvZpBeQHybdmSaWVJFXH+L4n9ZNH6b6qlI/K1q8a616Mr2dJKm2lV7XydRfgqSOOJZUmVosmQukAn0w+exvjmZTp2OT9i6leyTjJWe5/Qb5sP+4Wn7O9Tn5Evmy/DfD6WUU0XsNe1wcSE1DNZl91j21dXyefbNyfaMBRxeb0Du2QNjqZa64lnruz+YhNiQr9+w+5zfQ57udDbkdzxI7EDT37KpvkYinR3P65RALx6oL1O0TiMm1FUeRxyXp3UN5acDFwdJ99duid0Sv1tWT5T1O9a+ADknVXUV7wW0QkbvrFOo7uoVPE4LjZ/fwf/btNrw+8v098Fw3Y11P7rPd3YuyLfvYhmpdn1/k55UmXBUSX8uw6Z1OecH0McT4VxViWYP19sux3KD9HHkDcrIv2866S9SCORXb5G6k/1t4jiZYZ2e38iuqT5GwKvIgYdLms9e0hBfv5MP1bGi4iks5Fx2Wc4708K9nXG2uu/6dk/UEtKdPr5zvUL5Cn96XVRNxV7013J76LWwfEe2VmH8PMdrhXEudaesc4zlpAjH3U7x44qCvJ25Plf0i91pdbJuunXbFSmxPnZnadcyifvOeRxHFP72dllbPLkuV/OSCu1JOT9d9ec/0qnpTs490lyy4hWlxml7+d+i2uT82sv4rJTsgnDesQeu9tdVrg3xP4PFH5lm7nFxW2lf7Ovr/CPrck7m3p/tLeIEWuySy/kmoV5F3vS/Z3dIV1ltJbLuu+xtEI5rBkHzcxuLy+ITH0SBpf2qAga5ihmp6TbP+7fZbbhChPZpf9NsPN7r0//e/FU0TZMbufzzG4EnUf4lkhPV7fGCK+qs7I7GcN9RqvbEw+zkEtv6aI86a7/B3UH1boBck+11CvZWvXAsqH9No52U/RDPaDfDLZxocrrLM5USmTngNVWm1my7KrqJ/AfHyyz7LyDcAXkuUvpncCmdR84nzOrnc95XN1pGW7shb3RdL762Nrrl/Fe5N99GvVvREzlQbd160Md987PbONNVg2vMvbkd4bx01E18v0tbZg2aLXyVRvxTef3h++VcAniJvttkSN0z2BNxAtf7rLnUTvjW9Q8jJNlnb/3lOIH/avZ15vSNadYmb8xexFdDRxMW5HtAY9kOgieX5mudOJ7nXZdcc1tt2L6f0bL5yO6YFEQfvRxANv9nh+nOjG331/0YD9pMnL7xI/Ot3j8jmiVc4y4m+9P/ABehOKt1GtRdXBRO1Xdt2bgXcSYyFtQ5zPjwa+ysz5+neiS3p2vbLkZZr86hDn2QuJxP2mxEQwhxNJp+5+1gBHJesNSl5CdGdJ93cuUUv19eTVrzbz9QXbWAUcBxxB1FrtMv06iLgO3kQkg1Zm1hk0dED6w9ghCoTdY7M5USh9Kflr83hmhm3oMN7k5Ybkk0lpd6FB6iYvD6T3mNw6vd/vkv/++hWQlpAvWHdfVxJjBj+O6Oa0M3GtHEokDT9E7wPoOJOXEJUhaZxfJxJVy4ixD/cH/oUY1yd7D/pJst6g5OVe9P7u3EY8yB5P/th+qGD9uslLiPtH0T4/RNzndyfO9cOJ6ytd9iYG//7NteQlxEPx95J1VlC9VcV25JM3/1c/ZKkRRcnLPZj5ze2+9iLu3Y8gKgyPovf3Jvs6kWoJj39M1lsDPL9k+YOJyt8Ocf9alVm3SvLyW8n+vkfcE6uMR/ngZN219PbCybo7+YrrbHmlw3iSl0WNAE6jf8OOezCTXF2erFeWvFw7vd4RVEukHUYkP7LbL+tK/0R6f59Oo9qwHkuApxCt39cSf2M/j0z20SGSWf26Od+PGMOzQ2+F+DiTl29N9vWkGuvWTV5C73PHt4jyW9VxW+cR5Zz0enk/1ZI3uxLPsRcPiHcUycv02XoNcV73cxAzzwjpPahK8jJNCn6fePao2uusbvJyGb2VAefQv0fRZkS5NL0uyu7LMLeSlxD3gBOS5e8gytlVbUm+bOhwQipMXg77OptokVm32+J96d9Crt/rMiIhdnry+aDk5RSRUKiyj+ML1j+AfBKmyuta4mH258nn40peTtF7sxj0+imR9JhN8vKdREKlzn7vAB5e4297IvkfuUGvtUQBJW0BW5a8nEcksuue/0fS+xBRJXl5b3oL5P1e6WzGWS+jN7lb9zUoebmASNjU2ebFxA/5zZnPxj3T3ucz+1pLFM6qqpu8hEj8VzkWZbNPb0o8KMzm+xsU7yiSlzuTb3lT5XUrce9MW5hWGcQ7LSj1e51XsO4wyUuI1i1Vr8ns6xqqdUebi8lLiOEh0oeKVVQbBuEVyXrj7n4pjUpR8nI2r+XAW6j+ID4POLNgOycRY/E+gSibvJr4jck+CP43+dbmVZKXj+oT92rgcqKiqvsqKgOdVLDu74DXTcf5BGKc+OPJl/c+zkzStfsa1/BjzyyI8Vbgy0RDgH8hZgL/JTPHcy29reUGJS+z98lfEUMEvIS4/z6GSHa8md7nhw4xYeOgRNibCtbrEJWlHyDKfE8lxh18xfRnPyfKB9nly5KXAJ8t2MelRPLlmUQS5WVEojt7vNL4xpm83J38Mf9yjXWHSV6mzwNl10k6fFDXUvKVwN3X7cSx/H/A84gywQuIhgxfoLdhz7iTl/Po7R3XPUez96BXEYnG7D3o3cwM2dWhWvLy4QX76kxv9wryx7Zo0uC6yUuI6yRtHb+KqGh9NVHOeTExBNoNBbEdV2Efcy15CdGDLM27rKiwXtcLk3VfXD9kzTU7MnwB6w6iNeGHiBrB2Yy19g9UfyD+zXTcUD95CXEhvYPim0v2VZS8hHyt4aDX6czM+P3z5N/GOavwYqonaT/KTPfU2SYvIW40acGn6PVX6s1O3vVIqh3/m5kZKLlO8hKi1vlLFfbRmf5bXzm93jDJS4iWkSdX2FdZ8rK7nRMrxp19rSLO9yrJpIXEj1daY170Oo2Za3WSyct7J3EcWWPdYZKX84jJkrLJwX73g0HbeWWF7RS9riIePMoe5kaRvIQ4T7Ity8tef2FmKIRhkpdTxIPsoGt+lMlLiBZCRV2a+r1+SHn38qy5mryESLgcm6y7mvLhBSDf5fIKJjdOrjRbo0pe3kIk6KoO+ZK1N70tgwa9vkMkv+omLyF6R1XZx2cL1t2J/kNO9Xv9lLgnTCp5CTE2XdX41hKt3DZJPq+avKz7OoPqE+Y9l/oNRNLXoOTlAmIyoTrbfDfRsjb72TiTl5AvZ99M9W70wyQvoXyYs+zr4yXb2IiYNHM239+4k5cQrUoHPVenr+OJe1Dd5CXEs2uVfRxTsO4wyUuIZ8phKraPpdqzxFxMXkL87WlL1FX0TmZU5IeZde5g/BM0DTTMpBIarb9T7eJYSdzo/07cnM4nChGrRxTHj4nC18uJh6a7k0+G3kDU1n6R/MzHx5BPAlaJZwVRW/U2otC5H/GQmyY+z+2z/snEhDMvJbpVHEC+AHUTkWA9lrhYuzF9hXhA67q1QqzDup24OT+OqKV4MPkuHDcRBdfPkJ/t+L3MXJc3Drnvo4jv83nEAOu7MXN8ut1kvkrUTtWZka7rh8TxfwlR03UAM+fKGiKJ8SXi3OgmyX5CPllz0YB93EbUFn+OOCcfRH4Q/ZVEt6LvEt1juhOZnE/+esrOkl7mt0RSfCfinNyZ4nGkrqqwnQdPb+twYtKM/emtnb+DOL/PII7Nj4kWwlWsIloefJk4tw4n34XldiKB+kWiMNs9/z/ATJI8O4v0OJxCnNfd7lEvJFqYrO27xozPkp8MrMo6a4na3Q8QCeYDiO7T6dg2V1bYzoeIsXifRlTqPJjih5SriXvwyURh8yQG3/8+ShSEIb7HYZ1J/I0vIK6Tg8n/nt9CnItfJe6Dy6c//w75ay+dfKBIh6h8+Nj0fg4gjkd6bK8rWPd28tfjxRX213Xi9L6eQlRK/AO9E2lcRrRk+hr1Cvy3zCIuiGs3u36dpGxV5yX7qJpgXU20XP0DcQ107UF8b0X3sP2J7rRdHyPusdJcdhNRofUHooXSDyiegKWKc4heTF9hcKXQSuLafitRZhrGS4n7zsuJysI6ScSLgfsQZYhDyxdlNXH//w8mf0/4Z+J+9UrKEw/XE+Xd4+n9jSjzeaIVa53JR+8gWme+g+rl588TZaIjiTJ51W7LHeK8/ALFlYNZq4kGDf9OtKYsSwreTjyDfZB41pukDxPXCURZ6J+o1wKzrucTZYmXE2X7YZLttxC9zp5MfId1Zkq/hXhmPmqI/dZ1HnFdf5XBye5VRMLwvxj+HvQKotz7CuIeVPW8no0vE3/n+6jW+OZyYvz6Y4jr6a5qFXHvWc7M7PILiHvL+vQf+/hu5MfIPJbhcxMjM85ZUbVuW8LMALxXM/Pw20aLiVjnEbEOW/gcp3lEgmkJkaQaRfLoqcSPVNe7iAJmagEzyZermV3SpMh609tfPb39YX8Iy8wnxiTqjqdYJenSJpsz0yLgJmZaaIzKxsQ1cCMt+GGZ9ngiWdb1OKKbzbpoPeL4bkBUvtzEeCs/6lqfiG8hcW0MUymxLlhCJORWE/fRNv8urUs+zsyQFTcTLc+uby4cqZbNqDYJwc3MNAK4maj8GLV5xG/fPxIVaVsQlWJXTu/vl8QDeLYi7XHMJJtWE2Pz1bGEqFDehN4J1S6i/7jTU0RPmicSvUa6laBXE12OTyYqoi/NrPMI8uOAdoepGKfdiQfuBwE7EH/vVcTf9r/TMXR/8xaSH+f3BsortuYRCbxDie9rN2B74lguIH7rryMqDE8ivpuqrWOLbE90uX3w9N+15fS+7iCeCy4iEuGnAj+jeqV21jIiUfGo6f0tJb7TvxFlsK8xU05cSn7oqMuIxh/jsoBoeLDj9PtfEN/rIOn3eiNR8V/HhsR1spTe6+RCohJjkHlEou6hROX8dkT5fgPiXLmJSK79mejtdyqDn7mWkB+D8BqqN7zoF+NjiYYN9yTOsbXENZO9B2V7YD2GmQY2a4iEax2Lieuo6B50Mb0tT7cjP/breUT3/DruTcT9AKICYkviPnAl0YPr+0Rjmzq9nJYRCeBh4zqAmd6eEPeMQY0m6kr38WsiSVvFPOL3JtuzptsCd0XB8m9kZkiFDpG4/3OdYCUpq1+3caktfs3M+TnOArGk4WxFfpiRNzcbjiRpDnsZM783a4mW/5LaZT2i4qp7rdatVJOkHiYv1XYPIX+OHtZsOJIS72Tm+ryGarMrS5I0jIVES8fu786Xmg1HUoHnM3ONrqHaGPmSVMrkpdYF2UHkf4PDlUhtsRX5ibyOaDYcSdJdwLOY+d1ZTfXJ9iSN33xibPfuNfr5ZsORNFeYvNS6YHtifMjuefqEZsORNC07Y/FpjHfmYEmSICqxf87M789XS5eWNEkvZObavJkYo1SSZs3kpdYVuxEzVR9MzOQuqXkHMHNd1plxV5Kk2bg7MYN9d+zL+5QvLmkCNgAuYSa38Opmw+m1oOkAJElz3gVNByCpx5lNByBJukv6E/lZjyU17w5gh6aDKGMXIUmSJEmSJEmtZMtLad11JvCGzPtTmgpEkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJ0tw31XQAkjQGLwTmTf//pcAPGoxFkiRJmoQjmHnGvwT4YYOxSJIkqcQKoDP9MnEpSZKku4LVzJSBv9dwLJI0MvMGLyJJkiRJkiRJk2fyUpIkSZIkSVIrmbyUJEmSJEmS1EomLyVJkiRJkiS1kslLSZIkSZIkSa1k8lKSJEmSJElSK5m8lCRJkiRJktRKC5oOQJJaYh7wNuBemc+uAV4B3NhIRJIkSdJ4zQPeDhyS+exqogx8UyMRSZIk3QWsADrTrx9UWH4R8JXMOh3gAmD3cQUoSZIkjdhqZsqy36uw/CLga1gGliRJmrg6ycvNgF+SL7SdAmw5zgAlSZKkEauTvNwM+BX5MvBvsAwsSZI0EVWTlzsD55AvtH0bWDzuACVJkqQRq5q83Bk4l3wZ+FvABuMOUJIkSaFK8vIQ4CryhbbP4FjAkiRJWjdVSV4eQoxpmS0Dfwgn85UkSZqoQcnLhwM3Z5ZZC/zXpIKTJEmSxmBQ8vIRWAaWJElqhbLk5fOBVZl/XwE8Y6LRSZIkSaNXlrx8Afky8HLg6RONTpIkSXcqSl5OETXL2S4yNwIPmnx4kiRJ0sgVJS+LysA3AA9sID5JkiRNS5OX6wFfJF9ouxy4R1MBSpIkSSOWJi/XA44lXwa+DDigqQAlSZIUssnLk4GfkC+0nQncrbHoJEmSpNHLJi9/BfyMfBn4dGC7xqKTJEnSnbLJy/T1e2Dj5kKTJEmSxiKbvExfvwM2ai40SZIkZWWTl7eRL7itAA5vLjRJkiRpLLLJy7QMvBx4QnOhSZIkKSsd8/J15Atvq4HnNBadJEmSNHrpmJevp7cM/OzGopMkSdKdimYbfymwlnzh7fmNRCdJkiSNXtFs4y+ntwz8vEaikyRJ0p2KkpcARwBrMv+2FnjlxKOTJEmSRq8oeQnwInrLwP8y8egkSZJ0p37JS4CnA6vId6E5cqLRSZIkSaPXL3kJ8AzyZeC1wL9NNDpJkiTdqSx5CfBkYCX5BOa7JhadJEmSNHplyUuAp9BbBn7zxKKTJEnSnQYlLwEeDdyBCUxJkiTNDYOSlwCPwTKwJElS46okLwEeBNxCvvD2MWBqzPFJkiRJo1YleQnwSOB28mXgj2IZWJIkaWKqJi8BHgD8nXzh7dPAvHEGKEmSJI1Y1eQlwGHAzeTLwJ/EMrAkSdJE1EleAtwTuI584e1LwIJxBShJkiSNWJ3kJcAhwPXky8DHYhlYkiRp7OomLwEOBK4hX3j7GrBwHAFKkiRJI1Y3eQlwEHAt+TLwV7EMLKlF5jcdgCSNwT7AucDZwG+BX1VY5yrg+8AW0+udTYz7sxA4azxhSpIkSSOzD3AOUY49FTipwjpXAj8ENidfBp4P/HE8YUqSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEnS0KaaDkCdFwEPbjoKSZLUWhfA1BubDkIarc5LgAc2HYUkSWqtv8DUm8HkZQt0Pgs8v+koJElSa/0Wpg5tOghptDrHAM9tOgpJktRav4Gp+4LJywZ1ngk3bgcXvRUOWr/paCRJUltdvgou/yDc61KY+kjT0Uiz03kWXL8dXPYWOMAysCRJ6uPSVXDN++Hgy0xeNqZzNrB301FIkqR1ynyYWtt0ENLwOucBezQdhSRJWncsaDoA3el64G9NB9HrjHvA7ufBkjuajqTYOXvCBrfDTpc2HUmxK7aGGzeHfc9uOpJia+fBmQfAgac3HUl/Z+4PO10ES29pOpJif9kNFqyGXS5uOpJiV28B12wD+/2p6Uj6O/3Adp+DZ+0H210GW9zYdCTFLtgl/rvbX5uNo58bNoHLtof9/9h0JP219hzcCljWdBDSmF0HXNJ0EL3OuAfscS4sXt50JMXO2RMW3wY7XtZ0JMUu3wZu2qy9ZeDV8+GP+7f03j/tjANglwth41ubjqTYebvDeith5xY+wwJctSVctxXc/c9NR9Jfa8sf087aD5ZdCpvd1HQkxc7fFaY67S0DX7cpXLHMMvBQtgbu1nQQAqLlZaeTeR3VdER93ALs33QQJX4NtLkL3euBc5oOosT6wMqmgxjgeuChTQdR4gTgC00HUeII4OKmgxhgVdMBDHAF8OSmgyjxzelXWz0JuLLpIAZo6TnYeVVSVuhAZ17TUUmz0zkvOac/1XREfdwK7NV0ECVOAT7QdBAljgRanDRiY2BF00EMcANwWNNBlPgpcHTTQZR4OXBh00EM0NLyx52uAg5vOogS3waOazqIEk8HLm86iAFaeg52XpOWgS0AS5IkSZIkSWolk5eSJEmSJEmSWsnkpSRJkiRJkqRWMnkpSZIkSZIkqZVMXkqSJEmSJElqJZOXkiRJkiRJklrJ5KUkSZIkSZKkVjJ5KUmSJEmSJKmVTF5KkiRJkiRJaiWTl5IkSZIkSZJayeSlJEmSJEmSpFYyeSlJkiRJkiSplUxeSpIkSZIkSWolk5eSJEmSJEmSWsnkpSRJkiRJkqRWMnkpSZIkSZIkqZVMXkqSJEmSJElqJZOXkiRJkiRJklrJ5KUkSZIkSZKkVjJ5KUmSJEmSJKmVTF5KkiRJkiRJaiWTl5IkSZIkSZJayeSlJEmSJEmSpFYyeSlJkiRJkiSplUxeSpIkSZIkSWolk5eSJEmSJEmSWsnkpSRJkiRJkqRWMnkpSZIkSZIkqZVMXkqSJEmSJElqJZOXkiRJkiRJklrJ5KUkSZIkSZKkVjJ5KUmSJEmSJKmVTF5KkiRJkiRJaiWTl5IkSZIkSZJayeSlJEmSJEmSpFYyeSlJkiRJkiSplaaaDuCuq3M2sPfM+6/eCE+/orFw+tsHuAi4o+lA+tgNWA5c1nQgfWwNLAX+0nQgfcwD9gLObjqQEnsDlwM3Nx1IH7sAa4C/NR1IH1tMv85tOpAS+9Duc3Av4GrgxqYD6WPn6f9e1GgU/W1K3As9B2s7cnN42zbJh/Nham0j4Ugj0TkP2GPm/ZduhGdd2Vg4/e1N3FeXNx1IH7sR5fPLmw6kj62BjYHzmw6kj3nAnsA5TQdSYi/iGefWpgPpYxdgFXBp04H0sSVRBmnrcxi0tvxxp72Aq4Cbmg6kj52ADu19DtsU2Ao4r+lASrT0HHzT5vCWrbOfLGgqFKXudhHw7aajKPBG4ASgjYlVgJcSN9RvNR1IHw8G7gF8relA+lgAHEl74wP4T+DntLfw+wKiUNnWY3hv4P60Nz6At9Du+N4AnAKc0XQgfTx7+r9tPYYHAI+ivfFBa8/Bne9DHDtpDlt2IfCdpqMo8Cbg+0TlVRu9jEhctvHYATwMuDvwlaYD6WMR8B+0Nz6IMvrPgAubDqSPFxFJra83HUgf9wPuS7u/47fS7vjlPNWPAAAgAElEQVT+AzgZOKvpQPp4LtGIpK3H8CDg4bQ3PmjtObjL/YBHNh2FgGh52elkXkc1HVEftwD7Nx1EiV8DH2k6iBKvp901uusDK5sOYoDrgYc2HUSJE4AvNB1EiSOAi5sOYoBVTQcwwBXAk5sOosQ3p19t9SSgja2qslp6DnZelZQVOtBxyB+t4zrnJef0p5qOqI9biVZHbXUK8IGmgyhxJPDnpoMosTGwoukgBrgBOKzpIEr8FDi66SBKvJz2Jn67Wlr+uNNVwOFNB1Hi28BxTQdR4um0t3V8V0vPwc5r0jKwBWBJkiRJkiRJrWTyUpIkSZIkSVIrmbyUJEmSJEmS1EomLyVJkiRJkiS1kslLSZIkSZIkSa1k8lKSJEmSJElSK5m8lCRJkiRJktRKJi8lSZIkSZIktZLJS0mSJEmSJEmtZPJSkiRJkiRJUiuZvJQkSZIkSZLUSiYvJUmSJEmSJLWSyUtJkiRJkiRJrWTyUpIkSZIkSVIrmbyUJEmSJEmS1EomLyVJkiRJkiS1kslLSZIkSZIkSa1k8lKSJEmSJElSK5m8lCRJkiRJktRKJi8lSZIkSZIktZLJS0mSJEmSJEmtZPJSkiRJkiRJUiuZvJQkSZIkSZLUSiYvJUmSJEmSJLWSyUtJkiRJkiRJrWTyUpIkSZIkSVIrmbyUJEmSJEmS1EomLyVJkiRJkiS1kslLSZIkSZIkSa1k8lKSJEmSJElSK5m8lCRJkiRJktRKJi8lSZIkSZIktZLJS0mSJEmSJEmtZPJSkiRJkiRJUiuZvJQkSZIkSZLUSiYvJUmSJEmSJLWSyUtJkiRJkiRJrbSg6QDU9bN7Ae9tOooC6wOvA65qOpA+dgM2op3HDuB+wLa0N74FwHzaGx/AEuAVwKOaDqSPfYGdaO8xPAjYgvbGB+0/BzcBng8c2nQgfRwEdGjvMdwHWEp744PWnoOfOBhe2nQQ0pj95N608vpjEfAfwLVNB9LHrsBi2nnsAB4A3I32xreIKAe3NT6IMvC/Ao9vOpA+9gaW0d5jeE9gK9obH7S2/HGnpcALiWfaNroHsJb2HsP9gE1pb3zQ2nPw0/eEF+U+MXnZGtdvAxzWdBQF5gEHArc1HUgfGxIxtvHYAWwPbEB745sHTNHe+AAWAncnCsBttJR4eGjrMdyWKKC3Nb6uNse3HrAXsHnTgfSx2fR/23oMtyCOYVvj62phfFdv3XQE0vhdty2tvP6YDxwM3N50IH1sSPzGt/HYQfvLwPNpfxl4AZH82KHpQPpYSru/4+2wDDxbC4kk9ZZNB9LHZkQFfluP4ZbEMWxrfF0tjO/qbZqOQHfqnA2dTuZ1VNMR9XELsH/TQZT4NfCRpoMo8XrgnKaDKLE+sLLpIAa4Hnho00GUOAH4QtNBlDgCuLjpIAZY1XQAA1wBPLnpIEp8c/rVVk8Crmw6iAFaeg52XpWUFTrQccgfreM65yXn9KeajqiPW4mKq7Y6BfhA00GUOBL4c9NBlNgYWNF0EAPcQCuTCnf6KXB000GUeDlwYdNBDNDS8sedrgIObzqIEt8Gjms6iBJPBy5vOogBWnoOdl6TloEtAEuSJEmSJElqJZOXkiRJkiRJklrJ5KUkSZIkSZKkVjJ5KUmSJEmSJKmVTF5KkiRJkiRJaiWTl5IkSZIkSZJayeSlJEmSJEmSpFYyeSlJkiRJkiSplUxeSpIkSZIkSWolk5eSJEmSJEmSWsnkpSRJkiRJkqRWMnkpSZIkSZIkqZVMXkqSJEmSJElqJZOXkiRJkiRJklrJ5KUkSZIkSZKkVjJ5KUmSJEmSJKmVTF5KkiRJkiRJaiWTl5IkSZIkSZJayeSlJEmSJEmSpFYyeSlJkiRJkiSplUxeSpIkSZIkSWolk5eSJEmSJEmSWsnkpSRJkiRJkqRWMnkpSZIkSZIkqZVMXkqSJEmSJElqJZOXkiRJkiRJklrJ5KUkSZIkSZKkVjJ5KUmSJEmSJKmVTF5KkiRJkiRJaiWTl5IkSZIkSZJayeSlJEmSJEmSpFYyeSlJkiRJkiSplUxeSpIkSZIkSWolk5eSJEmSJEmSWsnkpSRJkiRJkqRWMnkpSZIkSZIkqZVMXkqSJEmSJElqpQVNB6CuazcGdm06igLzgB2A25oOpI8NgE1p57ED2AJYRHvjW2/6v22ND2A+sB3tjXEJ0NbrF2BrYCHtja+rzfEtALalvTFuOP3ftsa3LXEM2xpfVwvjO28L2LPpIKQxa3UZeEdgVdOB9LGYdpeBN6fdZeAl0/9ta3wQZeC70d4YlwBLaW98WxHPOm2Nr6vN8S0AtqG9MW5E3KPbGp9l4KFdsDnslvtkqqFIROdsYO+Z959bDc9f2Vg4/S0GlgNrmw6kjw2I2FY0HUgf6xEFjzuaDqTEYuD2poMosZj4ftc0HUgf6wMd2nsOLiR+ND0Hh7cBUTBa3XQgfaw//d/ljUbR3wLiPPQcrO01C+C96yUfzoeptv4mSxV0zgP2mHl/1Go4wjJwfRsQZaM2HjuwDDwKloFnxzLw7FkGnh3LwEP794Xw7oVNRyEgkpedTuZ1VNMR9XELsH/TQZT4NfCRpoMo8XrgnKaDKLE+7S30dl0PPLTpIEqcAHyh6SBKHAFc3HQQA7S1VUvXFcCTmw6ixDenX231JODKpoMYoKXnYOdVSVmhAx2H/NE6rnNeck5/qumI+rgV2KvpIEqcAnyg6SBKHAn8uekgSmxMe5NuXTcAhzUdRImfAkc3HUSJlwMXNh3EAC0tf9zpKuDwpoMo8W3guKaDKPF04PKmgxigpedg5zVpGdgCsCRJkiRJkqRWMnkpSZIkSZIkqZVMXkqSJEmSJElqJZOXkiRJkiRJklrJ5KUkSZIkSZKkVjJ5KUmSJEmSJKmVTF5KkiRJkiRJaiWTl5IkSZIkSZJayeSlJEmSJEmSpFYyeSlJkiRJkiSplRY0HYAkSZIkSZorOhsAdwO2BZYB2wCbAQuBDaf/uxGRj9gIWAGsBG5L/nstcAVw6fR/L4epFZP8SyS1g8lLSZIkSZJUw8IpOHQT+NVTgP2AuwO7AdsRicox6VwDnA2cCZwx/TobplaOb5+SmmbyUpIkSZIklejsBTwIuC9wd1h1ACycBxw+4UC2mn49KPPZKuj8ATgB+CFwKkxNOCxJ42TyUpIkSZIkZXT2JBKE3dc2+X9fOOF4Si0EDp1+vRG4Ef50MXx0I/jEhjB1a6PRSZo1k5eSJEmSJN2ldeYB9weeRLSmXNZsPLOyKey7KXwC4ErofAM4GjgJpjrNhiZpGCYvJUmSJEm6y+nMBw5jJmG5bbPxjMWGwD9Pv86HzqeAT8LUbU0GJameeU0HIEmSJEmSJqWzJ3TeT8zg/TPgZczNxGVqd+C9wF+h81roLGk6IEnVmLyUJEmSJGlO6yyEzpOh81PgHODVxMQ3d0VbAe9hJom5qOmAJJUzeSlJkiRJ0pzU2QE6bwMuAb4OPASn4u7qJjHPhM4Dmw5GUn8mLyVJkiRJmlM6O0+P73gBcCQ9s4UrY0/gROh8ATpbNB2MpF4mLyVJkiRJmhM6+0Lny8D5wIuAhQ0HtK6YAp4N/Ak6j206GEl5Ji8lSZIkSVqndfaPloOcCTwdmN9wQOuqrYHvQudd0DFfIrWEF6MkSZIkSeukzs7Q+RpwBtFy0KTl7E0BrweOh86mTQcjyeSlJEmSJEnrmM7G0HkncDbwFJyEZxweDfwOOns1HYh0V2fyUpIkSZKkdUJnHnSeA5wHvAFYv+GA5rpd/z97dx4mV1Xnf/x9uzv7QkjYEvYtssgmO7JJUHFGRkFQcQER0BlxFAdHUNFxGQdwQ0VHkVV0REQQUUd/DKgo26igMpIQFkmALCxZIIEkJOnz++NU0p2mq7q7uuqe09Xv1/Pcp/pW+tb55Pbt29/61l2A2yDslTqINJzZvJQkSZIkKXthb+Be4Dt49/AybQbcCmG/1EGk4crmpSRJkiRJ+ZsBeARgGlOAWyAckjqINBzZvJQkSZIkSaptI+AXEF6eOog03Ni8lCRJkiRJ6ttE4CYIm6YOIg0nNi8lSZIkSZL6Z3vgBgijUgeRhgubl5IkSZIkSf13KHBx6hDScGHzUpIkSZIkaWDOgHBc6hDScNCROoAkSZIkSWopq4BlwHPA0srXyyrPbwQUwCRgM2BzYGSamIP2DQi/gWJJ6iBSK7N5KUmSJEmS6vEccB8wC5gNzAQeAOZC0dn/lwmbAzsCewJ7AfsBryD/s0WnAhcB70qcQ2ppNi8lSZIkSVJ/PA/cDl/YCX43C256IxRrB/+yxZPAk8CdXc+FKcBRwN8DxwMTBj9OU5wC4QdQ/DJ1EKlV5f4phiRJkiRJSmcOcAHxJjUbQ3EMfGQu/PTpxjQuqykWQXEdFO8iHuH4TuCO5o03KBdAKFKHkFqVR17m4xAI/5k6xEt9exS87uOw9aLUSXr3/R1g0gT4uwzXHcBt+8CcLeCUTPOtbofL2+EfM80HcOk4OOyDsMubUifp3Q93hVHbwhsyXYd37wazJsOpmeYD+FZb3tvg5RvBvqfD3q9KnaR31+8dH9+U6Tr8845wz0Q4LdN8kPE2uGfqAFIJDs20Bh4Jx34Spi5NnaR3398OJo+BYzJcdwC/3hcenwonZ5rvhZFwdUem+/6KS8fCEf8C09/a7cm9SgzwAvAI8BDxiMiNgLfHKQDXvgzGToNjy16H9wHPAPsCW5U8di17Ab+A8Lf+L5Jt/VFx2UQ44D2w56tTJ+nd9XtCWyccl+k6vHdn+PNG8O5M80HG2+BL9nV+MpBMmAnsmjqFJEkaUtoHdg0xKTdhNjA9dQpJkjR0eNq4JEmSJEmSpCzZvJQkSZIkSZKUJa95mY9lxGt3ZGbutrDFAhj1YuokvZs3DUatgk0yvSbn0o1g2QTY+onUSXoXivgz3m5O6iTVzd0WNn0Kxq5InaR3C7aA9rWw2dOpk/TuuQnw7CTY+vHUSaqbs13e2+Bj28DkRTD++dRJerdw8/i4xZNpc1SzfBwsngLbPJY6SXXZboMTgSmpQ0hNlmkNPGc7mDYPRq5OnaR3udfASybB8vH51sCdbfHve5b7/oq528JmT8KYld2e3AiY3IAXXwksBQZRXy/YAjrWwKY5/f6OALYgfZ/jCaAf+45s64+Kx7aBKc/AuBdSJ+ndws2hCLD5U6mT9G7ZeFgy2Rq4Lo3a12nwwkwIodt0WepEVSwj7xsG3AlcnDpEDecAs1KHqGE0kGljer1FwIzUIWq4Gbg6dYgaziDeITJnmb4xXG8+cGLqEDXcUJlydQKwIHWIPmS6DYazetQKAYJnzWiIC7N7bNOXpE5UxXJgl9QhargbuCh1iBrOA+5PHaKGicCq1CH6sBg4fMOnwtm9/F0YyPQXCI2qaW4FrmjQazVQmFr5fw5mPQ12+lQ/w2Zaf6y3EDgudYgabgSuSx2ihpOAealD9CHTbfCl+zoLYEmSJEmSWtdfgDcAe0ORc7OnAYoFwAxYkvKsqJPjGW6SGsXmpSRJkiRJrWcu8FZgHyhuiqfYDgfFM/DWb8DitYkCbA8clGhsqSXZvJQkSZIkqXWsAD4N7AbFtcOnadndzc/ASU8CnYkCHJVoXKkl2byUJEmSJKk1/AzYHYpPQZHpjVbKcvNK4EuJBj8y0bhSS7J5KUmSJEnS0PZXYAYUx0LxaOowGfkEMDvBuIdAGJlgXKkl2byUJEmSJGloWgl8lHhdy1+lDpOfYhXw8QQDjwUOTDCu1JJsXkqSJEmSNPTcRWxaXgDFmtRhMnYD8McE4x6QYEypJdm8lCRJkiRp6FgBnAscBsUDqcPkrwjAVxMMvFOCMaWWZPNSkiRJkqSh4Q5gbyguhGJt6jBDyI+AJSWPuXPJ40kty+alJEmSJEn5+wHxaMsHUwcZeoqVwHUlD2rzUmoQm5eSJEmSJGWvmFc5BVr1uaXk8baCMKbkMaWWZPNSkiRJkiS1ul8BnSWO1wZMLXE8qWXZvJQkSZIkSS2uWAQ8UvKg40seT2pJNi8lSZIkSdJwUPb1QseVPJ7UkmxeSpIkSZKk4aDs5uWEkseTWpLNS0mSJEmSNBw8VfJ4HnkpNYDNS0mSJEmSNBwsK3k87zYuNYDNS0mSJEmSNByU3bx8seTxpJZk81KSJEmSJA0Ho0seb2XJ40ktyealJEmSJEkaDiaVPN6zJY8ntSSbl5IkSZIkaTiYXPJ4C0seT2pJNi8lSZIkSdJwsFXJ4z1Z8nhSS7J5KUmSJEmShoP9ShzrKSieK3E8qWXZvJQkSZIkSS0uTAR2LnHA2SWOJbU0m5eSJEmSJKnVHUi5PZC/ljiW1NJsXkqSJEmSpFZ3Qsnj3V3yeFLLsnkpSZIkSZJaWOgAji950P8teTypZdm8lCRJkiRJrezVwCYljvcYFF7zUmoQm5eSJEmSJKmVfajk8X5R8nhSS7N5KUmSJEmSWlTYn3jkZZl+UvJ4UkuzeSlJkiRJklrVx0oe7xnglpLHlFqazUtJkiRJktSCwmHAG0oe9IdQrC55TKml2byUJEmSJEktZkQBXAwUJQ98ScnjSS2vI3UASZIkSZKkxjp3IrB9yYP+Dor7Sh5TankeeSlJkiRJklrIsZvDuZMTDPzFBGNKLc8jLyVJkiRJUosIo+GpU2Fs2aeL3w/8rOQxpWGh7F9mrRdmArt2zf/oKTjxkWRxqjsQ+CvwfOogVexNzPZQ6iBVbANsAtybOkgVBfFnfHfqIDUcQPz5LkkdpIo9gNXAA6mDVDGtMv0xdZAaDgbuSh2ihv2BucBTqYNUsVvlcWbSFNVtBmwL/CF1kBoy3QY/NRX+bbseT7ZD0ZkijdQYYTYwvWv+2ifhrX9LFqe6A4H/A15IHaSKvYHlwMOpg1SxHTCZfGvgDmA/8q6BDwRmA0tTB6liT2AVMWNmrt0B3rx5+eOeOhuuWjyABTKtP9bbH5gDPJ04RzW7A4F8a+DNif0Aa+AB++w0OG/b7s945GU2Jj0D3JM6RS/2JzZlcn3TvjNxZ5rjuoP4OzaBvPMdQL75APYlNi/npg5SxXbACvJdh2uATck3H8BB5J1vH+BR8m1Qb1V5zHUdvgzYknzzQbbb4JS9iPsYqYVt/DR5fsC2P/EN8TOpg1QxnVif57juAEYCY8k33yhijZlrPoj5ZgOPpw5SxQ7AMrJbh5/eJU3jctY8uOrmAS50ENmtvw3sA/yNLBvUAGwNrCXfdbgbMJV880G22+CUfYgHPwxr/0j84aybdk4TI8yEELpNl6XJ0adlxE/VcnUn8Q5yuToHmJU6RA2jgRdTh+jDImBG6hA13AxcnTpEDWcQPzHN2erUAfowHzgxdYgabqhMuToBWJA6RB8y3QbDWT1qhQDB65WrXu9jwxp4pzQxwuwe23Sud+VdDuySOkQNdwMXpQ5Rw3nEU2hzNZF41GDOFgOHpw5Rw63AFalDbCi8CsKqXv52NntaDWHvOgJnWn+stxA4LnWIGm4ErksdooaTgHmpQ/Qh020wnN3z92w4Hnk5lfgp1jpjUwWRJEmSSjKNDWvgMamCSFLjhd2JH+aOTDD4l6D4c4JxpWHDT+8lSZIkSdIQFbYEfgFMSjD4Q8CnE4wrDSs2LyVJkiRJ0hAUphFPYd86weCrgZOhWJFgbGlYsXkpSZIkSZKGmLAFcAvx5oApnAtFznesl1qGzUtJkiRJkjSEhK2B3wG7Jgrwc/K+aZbUUobjDXskSZIkSdKQFLYBfg3skCjAXOLp4iHR+NKw45GXkiRJkiRpCEjeuFwNvA2KxYnGl4Ylj7wcmElA0W1+DbAsURZJkiSpDNbAkjIQtic2LrdNGOL9UNyZcHxpWPLIy/4pgE8BS4DFlekRYJ+EmSRJkqRm6q0GfhjYO2EmScNS2JZ4V/GUjcvPQfHthONLw5ZHXvZtJHA58I5uz80BXgc8kCKQJEmS1GQjgSuBt3V77lHg77AGllSqsB3xiMvtEoa4BvhEwvGlYc3mZW0TgOuA13Z77j5i0TYvSSJJkiSpuSYBPwaO7PbcH4BjgSdTBJI0XIWdiY3LLROG+A1wqjfokdLxtPHqpgG/ZcPG5f8Ah2HjUpIkSa1pS+Ib9SO7PXczMAMbl5JKFaaTvnF5P3AcFKsSZpCGPZuXvXs5cDcbXs/nKuDvgedSBJIkSZKabA9iDbxXt+euBF6PN+iRVKqwK3AbaRuXc4FjoFiaMIMkbF72ZgZwO7B1t+cuBE4FVidJJEmSJDXX0cQaeKvKfAA+Dbwba2BJpQq7EG/Os0XCEI8DM6B4ImEGSRVe83JDJwOXASMq82uAMwHvKCZJkqRWdQpwKRvWwO+rPCdJJQq7AL8CpiYM8QTwKigeSZhBUjc2L7t8ELgIKCrzy4G3AP+dLJEkSZLUXOcA57NhDfxm4BfJEkkapsKuxMZlyiMubVxKGbJ5Ce3AN4F/7PbcQuK1fe5JkkiSJElqrnbgW8B7uz23gFgD35skkaRhLOwF3AJskjDEQuA1UDycMIOkXti8hKuB3bvNzwL+DpiTJI0kSZLUfN9jwxp4JvA64LE0cSQNX2Ef4H+AKQlDzAeOhOKhhBkkVeENezYs2gA+gI1LSZIktbbuNXAA3o+NS0mlC3uRvnE5j3iquI1LKVM2L+GuHvPXAQelCCJJkiSVpHsNXADXAwckyiJpWAp7EE8VT9m4fJx4xOWDCTNI6oPNy3gnxa93m59E3IHOSBNHkiRJarp/Ai7vNr8x8UYZR6WJI2l4CbsAN5P2GpdziUdceo1LKXM2L+NpMh8AvtLtuXHATcBrkiSSJEmSmqsTOIOX1sA/BV6dJJGkYSK8jPR3FZ+DdxWXhgybl1EAPgR8pttzY4nF23FJEkmSJEnNta4G/my359bVwG9MkkhSiwvTiY3LqQlDPAgcCsWjCTNIGgCblxv6N+DcbvMjgR8C70gTR5IkSWq6T7JhDTwKuBZ4U5o4klpT2B74DTAtYYiZwBFQzEuYQdIA2bx8qQuBfyV+Eg3QAVwFvCtRHkmSJKnZLgQ+0m1+JLGBeUqaOJJaS5gI/IS0R1zOAo6GYmHCDJLqYPOyd18kXsS8szLfDlwBvD9ZIkmSJKm5vkDvNfBpyRJJagFhBPAjYI+EIf4MHA7FgoQZJNXJ5mV1lwAnA2sq8wXwNeCsZIkkSZKk5voW8F66GphtwKXAB5MlkjTUfY20NwK7l3jE5TMJM0gaBJuXtf0X8XqXqyvzBXAR8bpAkiRJUiu6DHg7G36I/xXgvGSJJA1R4RzgHxMG+CPwaigWJcwgaZBsXvbtWuB4YGW35z4NXJAmjiRJktR0PyDWwKu6PfdZrIEl9Vs4HviPdOP/fiVwFBSL02WQ1Ag2L/vnZ8BxwIpuz51DvC6QJEmS1Ip+Su818IVp4kgaOsJOwJUk6zk89gi8ZiEUy9KML6mROlIHSOBq4Hfd5h/p53K/BF4BbNXj+QmAO0RJkiTl7ErgN93mH+3ncr8A9gW27PG8NbCkKsIo4IfAxEQBfgcH/Bie9Ya7UosYjs3LR+h/w7KnByqTJEmSNJQMpgaeVZkkqT8uAPZJNPbvgdfDk+9MNL6kJvC0cUmSJEmS1ADhYOCfEw1+H/A6KJ5LNL6kJrF5KUmSJEmSBimMBq4C2hMMPgs42pvzSK3J5qUkSZIkSRqsDwPTE4z7CLFx+XSCsSWVwOalJEmSJEkahDANOCfBwI8DM6CYn2BsSSWxeSlJkiRJkgbjE8D4ksdcDrweirkljyupZMPxbuOZunNP4COpU/RiJHAaMC91kCq2IW7HOa47gBnAZuSbbwTxmjS55gMYC7wN2Dd1kCp2ArYg33V4CLAx+eaD+EFazvkmAm8Ctk8dpIrdKo+5rsNXABPINx9kux/8zqFwSuoQUpPdvhdZ/v4xAjgDeDJ1kCq2BgJ5rjuAI8m7Bh5Ntvv+9cYA7wAOSh2kih2BKcBHYK8JsOb0ctsLIcCnrofPHAMc08s3HNaVL1u518DjgROBnVMHqWJXoJN81+H+xPcRueaDbLfB7x0ed39dikRJRJhJ/GWruOZZeNvCZHGqm048FH9F6iBVbA+sAnI9TWBT4g7rkdRBqmgj/jGanTpIDTsDC4ifrOZoW2At8ETqIFVMrkwPpw5Sw8vIexvcCXgGWJo6SBVbE/+eP5Y6SBUbEd/APpQ6SA2ZboMf3xj+fbMeT7ZD0ZkkjtQQYTYbXBPve0vhnTk2CKcT96srUwepYntitgWpg1SxGfGDq5xr4J2AB1MHqWFn4nuc51MHqWJbYA0wD762Kfzz5HKH/+zT8MlaN+eZQvwA3xq4fjsTP8DJ9e7t2xA/xHk8dZAqJgGb4DZYh09Mhs9smjqFgNi8DKHbdFnqRFUsA/ZMHaKGO4GLU4eo4Rzine9yNRp4MXWIPiwiHsGaq5uBq1OHqOEMYE7qEH1YnTpAH+YTP3XO1Q2VKVcnkO+b63Uy3QbDWT1qhQDBS/5oiAuze2zTl6ROVMVyYJfUIWq4G7godYgazgPuTx2ihonEAyBythg4PHWIGm4FroAwAsKCXv5eNXP6Tj/ynUm+zfN1Mq0/1lsIHJc6RA03AtelDlHDSeR7Bus6mW6D4eyev/cWwJIkSZIkqR7/QLyEUlnuA95T4niSMmDzUpIkSZIk1aPMs1NeBE6BIvejZiU1mM1LSZIkSZI0QOPagdeVOOBXoPhzieNJyoTNS0mSJEmSNEBv2px4/dAyPAV8rqSxJGXG5qUkSZIkSRqgo6eWONiXoMj1rteSmszmpSRJkiRJGqD9y7pRzzLgmyWNJSlDNi8lSZIkSdIAjChghyklDXY9FMtKGktShmxeSpIkSZKkAdh3LIxsL2mwn5Q0jqRM2byUJEmSJEkDsP+EkgYKwJ0ljSUpUzYvJUmSJEnSAOwwuqSBHoLiqZLGkpQpm5eSJEmSJGkAtiqreXlPSeNIyr1+4jAAACAASURBVJjNS0mSJEmSNABbl9W8fLikcSRlzOalJEmSJEkagGllNS8fLWkcSRnrSB1AkiRJkiQNFaGAVWU1L/9W0jiSMuaRl5IkSZIkqb8mw6iyegnzSxpHUsZsXkqSJEmSpP4aX+JYi0scS1KmbF5KkiRJkqT+GlPSOAF4tqSxJGXM5qUkSZIkSeqvspqXy6FYU9JYkjJm81KSJEmSJPVXWc3LpSWNIylzNi8lSZIkSVJ/ldW8fKGkcSRlzualJEmSJEnqr7Kal54yLgmweSlJkiRJkvpvREnjrC1pHEmZs3kpSZIkSZJyszp1AEl5sHkpSZIkSZJy45GXkgCbl5IkSZIkKT9e81ISYPNSkiRJkiTlxyMvJQE2LyVJkiRJkiRlyualJEmSJEmSpCzZvJQkSZIkSZKUJZuXkiRJkiRJkrJk81KSJEmSJElSlmxeSpIkSZIkScqSzUtJkiRJkiRJWbJ5KUmSJEmSJClLNi8lSZIkSZIkZcnmpSRJkiRJkqQs2byUJEmSJEmSlKWO1AG0zto28v15tJNvNohN+FzzrfuAINd87ZXHXPOtk/PPuCDvfEPlZ5x7vpx/xkXlMdd8ue8H18kw35q2LGNJDWUNPAg5r7vc/zYNlfoo021wRTuMKWGcUFD//9/6ozFy3s+0Efc1uebLfT+4Tob5XloDF1W+U00XZgK7ds1f0QmnrU0Wp7oRwBogpA5SxQigE8hx3UEsONqA1amD1DCC/POtJf6cc7Rur7omaYrq3AYHz21wcNwG6/ahNvhye48n26HIdVuU+iHMBqZ3zV/aCe/JsY6zBh6c3Pf96xoeueaDrOuP49vg+p5/n5rg9gCH1Vvf5L4NQrb1x3oZb4OANXAjZLoNfrgNvlDCPkb9EGZCCN2my1InqmIZsGfqEDXcCVycOkQN5wCzUoeoYTTwYuoQfVgEzEgdooabgatTh6jhDGBO6hB9yPAP5gbmAyemDlHDDZUpVycAC1KH6EOm22A4q0etECB4yR8NcWF2j236ktSJqlgO7JI6RA13AxelDlHDecD9qUPUMBFYlTpEHxYDh6cO0btwXC9/n5ox/XYQIc8EHmnU/7hJMq0/1lsIHJc6RA03AtelDlHDScC81CH6kOk2GM7uuT+wAJYkSZIkSZKUJZuXkiRJkiRJkrKU4YU5JUmSJElSpn4Gm/8WOh6Hef/cxHFyvZahpJLZvJQkSZIkSf1UrCY2Fl+EYknqNJJan6eNS5IkSZIkScqSzUtJkiRJkiRJWbJ5KUmSJEmSJClLNi8lSZIkSZIkZcnmpSRJkiRJkqQs2byUJEmSJEmSlCWbl5IkSZIkSZKyZPNSkiRJkiRJUpZsXkqSJEmSJEnKks1LSZIkSZIkSVmyeSlJkiRJkiQpSzYvJUmSJEmSJGXJ5qUkSZIkSZKkLNm8lCRJkiRJkpQlm5eSJEmSJEmSsmTzUpIkSZIkSVKWbF5KkiRJkiRJypLNS0mSJEmSJElZ6kgdQJIkSZIkDQVhG+DXMH8qdOwPmx3R5AG/C8WnmjyGpMzZvJQkSZIkSf0xAtgBpq2bn9Dk8aY0+fUlDQGeNi5JkiRJkiQpSzYvJUmSJEmSJGXJ5qUkSZIkSZKkLNm8lCRJkiRJkpQlm5eSJEmSJEmSsmTzUpIkSZIkSVKWbF5KkiRJkiRJypLNS0mSJEmSJElZsnkpSZIkSZIkKUs2LyVJkiRJkiRlyealJEmSJEmSpCzZvJQkSZIkSZKUpY7UAbTefhA+lzrES104Et70z7DTU6mT9O4/t4HJI+CtGa47gJ8dDI9sCh/MNN+qDvhyO3w003wAXxgDR78b9jkqdZLefXsnGD0VTs50Hd6yN9w3Cf4l03wA57flvQ1+aQIceBIcunfqJL27Yrf4+O5M1+Edu8BdE+DDmeaDjLfBA1IHkEqwf7418FvOgu0WpU7Su29sDZsW8OYM1x3ATa+ERzfLtwZeNgq+0Q7nZpoP4PNj4LVnwF6vTZ2kh41LHu+g+vYR/29fmDkZPpTxzzjb+qPii+PhkLfDIfulTtK7y3eFtk44NdN1+Nvd4Q8T4exM80HG2+CBPZ8oUqQQQJgJ7Jo6hSRJGlLaoehMHUKqX5gNTE+dQpIkDR2eNi5JkiRJkiQpSzYvJUmSJEmSJGXJa17m40XghdQhXurZjWD8cmhfmzpJ75aPj9nGrEidpHerRsGLI2HCstRJqls6CSYtTZ2iumc3gnHPQ8ea1El6t3x8vNbK2Ax/fyFufytHw8TnUiepLvdt8LmJcR8zYnXqJL17flx8HPd82hzVrB4BK8a4DdZlFDAmdQipyayB67JsQqyNcq2BV46O+/9ca+BQxJ9xlvv+imxr4DZgYonjrQLq2M5XjYqT9Uf9nt0ovsexBq6P78MGwRo4H2EmhNBtuix1oiqWAXumDlHDncDFqUPUcA4wK3WIGkYT3zTkbBEwI3WIGm4Grk4dooYzgDmpQ/Qh04JovfnAialD1HBDZcrVCcCC1CH6kOk2GM7qUSsECJ41oyEuzO6xTV+SOlEVy4FdUoeo4W7gotQhajgPuD91iBomEptiOVsMHJ46xEuFHXv529TMqd73emcCjzTyf94EmdYf6y0EjksdooYbgetSh6jhJGBe6hB9yHQbDGf33BdYAEuSJEmSJEnKks1LSZIkSZIkSVmyeSlJkiRJkiQpSzYvJUmSJEmSJGXJ5qUkSZIkSZKkLNm8lCRJkiRJkpQlm5eSJEmSJEmSsmTzUpIkSZIkSVKWbF5KkiRJkiRJypLNS0mSJEmSJElZsnkpSZIkSZIkKUs2LyVJkiRJkiRlyealJEmSJEmSpCzZvJQkSZIkSZKUpY7UASRJkiRJ0pCwCrgHHpgOI16EHec0ebzHm/z6koYAm5eSJEmSJKkfiieA/YBbgbnAu9PmkTQceNq4JEmSJEmSpCzZvJQkSZIkSZKUJZuXkiRJkiRJkrJk81KSJEmSJElSlmxeSpIkSZIkScqSzUtJkiRJkiRJWbJ5KUmSJEmSJClLHakDSFLrCKOArYGplcfxwMaVfxwNrKxMS4H5wDzgYSheLD+rJEmSJEn5s3kpSXUJI4ED4Cvbw77T4bBZwE4MfL+6GsJM4B7gf4BboHimwWElSZIkSRqSbF5KUr+FzYHjKtOhwFg4a7AvOgLYqzK9G+iEcCtwJfBjKFYOdgBJkiRJkoYqm5eSVFMYBZwAnA4cBrQ3ecA24NWVaT6EC4FLoVjR5HElSZIkScqON+yRpF6FaZXG4RPA94AjaX7jsqdpwFeBmRBeV/LYkiRJkiQlZ/NSkjYQtoXwDeAR4CPAJokDAWwH/DeEKyGMSR1GkiRJkqSy2LyUJADCJAhfAh4E3ke8O3hu3gXcGRuskiRJkiS1PpuXkkQ4FXgI+BdgZOIwfdkb+B2E6amDSJIkSZLUbDYvJQk+QB6nh/fX1sBtELZPHUSSJEmSpGayeSlJQ9MWwM8hbJQ6iCRJkiRJzWLzUpKGrl2Bb6YOIUmSJElSs3SkDiBJQ9+ylTDhv4E/AQ8ADwPPAkuJHxJNAjYGtgL2B14JHAA04s7hJ0H4IRQ3NuC1JEmSJEnKis1LSarPX4EfwKF/D3c9DJ0n1/jeRZXHPwKVJmOYBJxKvLP5ToPM8mUIP4di9SBfR5IkSZKkrNi8zMZPXgP8LHWKXowhnpa6JHWQKl4ObAPkeuOSXYBp5PmzhXhUYAf55gOYCFwAPNm8IR7eoX/9wxVr4CdPwMWPwZ3LgIOJp27vwIDXYVF5bH8Q/mM1nPUyGFnvpTy2h0/eAzzWy7/tSLw+Zs4/43byzrcJcB5wSuogVRxUecx1HW4NTCHffJDtNnj+9vDR1CGkJrvhGLL8/WM08G3gudRBqtgN2BLYOXWQKnYFppLnzxZgRGXKNR/EGvgLwNOpg1SxL7A7+a7DnYi/I7nmg2zrj/UmA58ETksdpIqDgbXkuw63ATYl33yQ7Tb4+R3gIxs8Y/MyGx2rgeWpU1TxPPlmWwvkvO5WAZ3km6+98phrPoAArKCpGds6a//7ohVw5QPw+Ufh6Z5HN64F1lB3vrXAOffBT+fCD14JW06o73VO3xY+O7OXf1hJXIc5/4wh73ydxPWYa8Y1lcdc863AbbBOI1alTiA1X9Y18Avkm62TvGvgF8m7Bl73PjjXfFBKDTwovg9rjJzzBfKugVeTd425ErfBOlkDZyTMhBC6TZelTlTFMmDP1CFquBO4OHWIGs4BZqUOUcNoYnGZs0XAjOYOEf7U4/dx3fQUhHMg1Lo25c3A1Q3KsTmEWVWy9Gfar5cXPQOY05h8TZP76e7zgRNTh6jhhsqUqxOABalD9CHTbTCc1cvvuTdb1BAXZvfYpi9JnaiK5cQzaHJ1N3BR6hA1nAfcnzpEDROJza2cLQYOTx2ihluBK1KHqOFM4JHUIfqQaf2x3kLguNQhargRuC51iBpOAualDtGHTLfBcHbPGtgCWJJeaiXwOWB7KC6EYkU5wxZPAq+l/lPkj29gGEmSJEmSkrN5KUkb+gmwOxTnQfF8+cMXjxFv4lOPIxsYRJIkSZKk5GxeSlL0JHAiFG+E4m9poxQ3AHfVseB+EEY2Oo0kSZIkSanYvJQk+CHxaMsfpQ7STT3X0RwBbNfgHJIkSZIkJWPzUpIozodiUeoUPdxIvDvdQO3Y6CCSJEmSJKVi81KSslQsBB6uY8HJjU4iSZIkSVIqNi8lKV/z61hmfMNTSJIkSZKUiM1LScrXgjqWaW94CkmSJEmSErF5KUn56qhjmaUNTyFJkiRJUiI2LyUpX5vVscyShqeQJEmSJCkRm5eSlKVQAHvVseCsRieRJEmSJCkVm5eSlKeDgUkDXGYpMLcJWSRJkiRJSsLmpSTl6e11LHMbFKHhSSRJkiRJSsTmpSRlJ0wBTqljwR83OokkSZIkSSnZvJSk/HwWGDfAZZ4HbmpCFkmSJEmSkrF5KUlZCQcB761jwaug8E7jkiRJkqSWYvNSkrIRNgauYeD75heBixqfR5IkSZKktGxeSlIWwghi43K7Ohb+EhSPNDaPJEmSJEnp2byUpORCAVwKvLaOhR8F/r2xeSRJkiRJyoPNS0lKKhTA16nv7uJrgHdC8UJjM0mSJEmSlIeO1AEkafgKbcAlwOl1vsC5UNzRwECSJEmSJGXF5qUkJRHaiKeKv7vOF7gG+HLj8kiSJEmSlB+bl5JUutAOXE59p4oD/BY4FYrQuEySJEmSJOXH5qUklSq0A1cC76zzBf4KvAGKVY3LJEmSJElSnrxhjySVZtCNy9nAa6BY2rhMkiRJkiTlq8qRl2EscATwKmAqMAVYBjwB/BH4JRRLelmug+rXb1sKPAP8BYpFVcYdRddplLOg+F3f/4XwWmBboBOKy/r+/l5fYypwEPAKYBqwCRCAp4C7gR/3/v+VpP4K7cBVwDvqfIEHgaOgWNCwSJKkHsJY4EhiDbwFG9bAfyDWwL18gNRnDfw0sQZeXGXc0cDJlZmZUNzej6zHANsAa6G4vO/v7/U1etbAmwKdxBr4LmIN7AdmkiQpqR7NyzAKOAv4CDC5xnIrIVwDfAKKed2eH0G8c24tAcKvgI9C8Yce/zau2/LfBvrRvOSfgDcAq4E6mpfh/cDFNb7hDOBrEM6D4isDf31JCh3AfwFvrvMFHsDGpSQ1URgNfAj4MLVr4BUQvg98Eor53Z4fRf9q4FuJNfAfe/zbhG7LfxPoR/OSM4HXA6uI11EeoPBBoFZtewZwMYSPQlGrVpYkSWqqbs3LMAX4MXBY5Ym1wO8r0zPAWGAHYAbxyMRTiZ/U7lbltRcBc7rNbw5sBRSV17gTwpuguKkh/5P6ja48LiP+X+cBTxIL15cDBxKbqhdBGAPF+UlSShqi1p8qPpjG5Qwbl5LULGET4EbglZUn1gL/SzzS8hliHbg9cDTxSMzTiDXwy6u84DPA3G7zWwBbEmvgo4EjIBwPxc8a+/8YsHU18HN01cBPEWvgPYADiP/3r1Vq4M8nSSlJkoa9SvMyjAR+RizEAH4JfACKh166SGgnvgn/LLVv+PNTKE7tsewk4KPEIzs7gCsh7Jj4dJRfAYcC/wvFmpf+czgCuJ5YrH4KwuVQPFVmQElDVWgDvgO8vc4X+AvwaiieblwmSVKXMBL4ObFRB/ALYg38cC/f2w68hb5r4JugOK3HshsDHwfOJp6ptK4Gfm5w+QflFmrXwEcSa+DJwGcgXAHFM2UGlCRJgq4b9nyarsbltcDre29cAhRrobgG2Bv40cCGK5ZCcQ7w3coTk4HjBvYajVbcC8UdvRdtAMVtwGcqMyOBo0oKJmno+yL1Ny7/BBxt41KSmurf6Wpcfh84tvfGJVRq4O8Ta+DrBzZMsQSKD1fGgHgW0xsHHreRinv6qIF/Q1w/EE+Lf1U5uSRJkjbUBmEy8P7K/OPAGbE460uxHIqP1Tlu96bnvnW+Rpn+1O3rjZKlkDSEhNOI10+rxz3ExqVHuEhS04RNiNeNhHia93v7WQMvg+LjdQ5qDSxJkjRAbcBJwPjK/H/Ggqzpul8HaFIJ4w3WTt2+rvJpvCStE44AvlXnwncSr3FZ5Y60kqQGeRvxmu4A34gfzDfdnG5fWwNLkiT1QxtwRLf5H5c07rbdvk55vct+CBsBH6zMPADcljCMpOyFjYjXuax1PbRq/gT8PRTPNjaTJKkX1sA1hUnAByozM+nfHdAlSZIargPYp/L1MuDBksY9vtvX95Q0Zj+EA4FtKjOTgF2AdxGvzfkY8Jbq1wWSJAC+zoZvTvvrL8QjLjN/MytJLWNdDbwUeKSkMd/U7eucauCDgK0rMxsTa+BTiDXwXOCt1sCSJCmVDuJdtAGegiI0d7gwATiXWAwBLAFubO6YA/IB4ilE3XUSL1b+RY+GklRbOID6btAzG3hNvKGDJKkkZdbAE4GPAe+oPLEY+ElzxxyQDwJv7fHcWuKd1b9kDSxJklLqoOt6l883+LWPhfDHbvNTiEc1rrvD+VrgPZm9WX8K+BvQDkwl3l28jdhw3QXCmVA8lTCfpLx9DigGuMwi4PXuWySpdM2qgf+hHzXw6Zk1BHurgduJNfDLILwfiqcT5pMkScNYG/F0cegq4BplCvEuiuum7egq2n4LHA7Fj3pfNJXiQ1DsCMV2wBjiRco/T2xGnAD8CsK4ZPEkZSzsAhw9wIXWAidC4U0QJKl8zaqBN6F6DXwbcBgUZV1js5+KD/ZSA3+RmPvNwK0QxlZdXJIkqYk6gGeI17PZDEIbFJ0Neu27gKu6zT9PPEXmT1As7MfybX1/y4C+b4CKTuL1j86B8CjwTWB34qnl5zd+vD/vRNepRDkZARwL7Jk6SBWbE38uOa47gAOIv1+55htBPLIi13wAo4HXEI8EydFWwET49QXwqgEuetfNcMiWNHf9HwpMaPIYg9VG3vnGEX+4o1IHqWL7ymOu6/Bg4h2dc80H2W6DP94PjksdopU9A2wEbA6haOCp43cAV3ebX1cD3wvFk/1YPpca+F8rNfA3gD2A9xM/1G+we3cmy98/OoA3Ak+kDlLFFsSfS47rDmA/4vVTc803hqFRAx9D1z0ZcrMlMWOu6/AQ4j4+13yQbf2x3ljgKGItnKMdgDXkuw4PI667XPNBttvgTfvCP2zwTAHhGrqucbM7FDPrHyCMAV6ozFwFxakDXH40sKIy830o+nHtuHAr8RfqWSgmDWy8fudqA+YTG2W3QPHqBrzmTGDXrvnvPQ/vzOkU+nW2BJ4GXkwdpIotgNXEU29zNIlYHC1IHaSKgvgzzrUwh5hvMV37htxsDqyF/2uHl4/o/2KLO2G7BbCsyddZY0Jlmt/kcQZjK/LeBqcBzwHLUwepYlPi73Kulx4YD0zEbbAO54yHC3rWNu0N/KB5mAs/BE6szOwCxexBvNY4uvYRV0Bx2iCW/y4UJ/djmd8Q75i+BIrJAxuv37naiDXMZsD/g+KYBrzmbGB61/x3nod3WQMP3BbEbItTB6ki9xq4jfj3PcN9/3pDpAbmmdRBqphIrEGsP+q3JfGmdo2+vEqjbAYE4r46R9bAdfvYBPjcRj2eDO+FECrTJwY3QBjT7bWurPM1nq8s/+t+fv8Dle9v8mmX4Y7KOA26G2WY2W1dBQiXNeZ1G24Z+R51CXAncHHqEDWcA8xKHaKG0eRblK+zCJiROkQNN8Mm10BY0+N3uq/pcyXlOwOYU9JY9VqdOkAf5tPV4MjRDZUpVyeQ75vXdTLdBsNZvew7mnS03XAU3tdtvX5skK81rttrXV7na6ysLH9LP7//ocr3P1jfeP3OdVdjxwmze2zTlzTmdRtuOfGu67m6G7godYgazgPuTx2ihonAqtQh+rAYODx1iBpuBa5IHaKGM4lHkecs0/pjvYXkfQrIjcB1qUPUcBIwL3WIPmS6DYaze9bAbcA1wLoLhv8ThI0TJgT4c+VxHwgja39rmEK8Jk/35Zpli8rjsprfJWkYOnoy8dSjgfh5M5JIkvrt+3TVdWdC6PkJf9nW1bKvgNDHkfxhU+Lpct2Xa5Z1l22xBpYkSUm0QfEc8OXK/FTgir6bhgBhMoSvNiHTLyuPG9F1Ons1p9PVMPhlrW8cnPAqugrE3zdvHElD064TB7hAAP7ajCSSpP4qltJ15No0Yg3cj8t/hClNroE3Bt7Sx/eeQdc1L5tZAx8NbFuZsQaWJElJrCt6/oN490OIF6b+Hwh79b5IGA3hPcD/Aa9rQqZLiNcWA/gKhCOr5Hgj8MnKzDzgvwY+VBgbT1UJe1T59wLCscAPK090VvJJUjfbDPQi1nMqHxxJktL6LHB75evjiTVwlcvlhNHxckvcB7y2CVm+Rdd1L78KocrpouF44im5EK9T9f2BDxXGQ/hWHzXwPwA/qDyxFmtgSZKUSEd8KNZAeANwLbEYOxz4M4T/A/5AvAjvWOIp2uvumARNOU2leArCacRiaWPg1/FaO/wvsam5MfBK4BWVBVYCb4eingsZtwHviVOYC/yFeF2HNuLFaV8ObN3t+z8GxT11jCOppW0xdoAL/K0pMSRJA1Ss6fZB9auJN8D5S48aeBxdNfC6/f29TciyEMLpxA/kJwO3QbiTWAMvI9bAhwL7VBZYAbwNipV1DNYGvDdOYS6xIbuAeEbTlsDubFgDnwtFs09PlyRJ6lVH15fFsxBeTzwN5Tzi6TN7VKaelgCXAuc3J1bxIwjHAF8HXgYcXJl6+hPwHij+WOdAa4EHiXc83Jau02J6mgP8a8wlST2N6ej7ezaQ612rJWkYKpZC+DtiM+/jxMsoVauBFxNr4AualOVaCEuINfDOwCGVqad7iTVwvR+qr6F/NfCjwIehyPmmYJIkqcX1eMNdrAG+CeFS4EDgVcQb1WxCfLP9BPF6N7+q8invKuKn1jDoO4sWt0DYjfgJ+KHANsTrYC4lHrX0W+AuKMIgxlgBvAzCNsBRwPbEpu3EbuPcBdwORWf940hqbaMHerOeeo4UlyQ1TbEG+EblztcHAUcSa+BNiUc9Pk6sgX9dpQZeSVcNPH+QWW6GsCtdNfDWxBp4CbGZeBsUdw1yjBfoqoFnANsRj7icQFcNfCdwhzWwJElKrcrRQsUa4I7KNABFJ3DLIDP1fL1fV6YmKh4DrmruGJJa10ibl5LUEoo1xGtg3t7Xd/ZYbi2NrYHXAr+qTE1UPAZc2dwxJEmSBqet72+RJNVWDHQBm5eSJEmSJPWDzUtJKl89N1eQJEmSJGnYsXkpSeXz+mGSJEmSJPWDzUtJkiRJkiRJWbJ5KUmSJEmSJClLNi8lSZIkSZIkZcnmpSRJkiRJkqQs2byUJEmSJEmSlCWbl5IkSZIkSZKy1JE6gCQNfd99GPb+6gAW+H3TokiSJEmS1EJsXkrSoH15Dnz5wtQpJEmSJElqNZ42LkmSJEmSJClLNi8lSZIkSZIkZcnmpSRJkiRJkqQs2byUJEmSJEmSlCWbl5IkSZIkSZKyZPNSkiRJkiRJUpZsXkqSJEmSJEnKUkfqAJI0tLUXsPlImLfxABYKUCxtWiRJkiRJklqEzUtJGpTpo2DmW4C3DGChJcDkJgWSJEmSJKlleNq4JEmSJEmSpCzZvJQkSZIkSZKUJZuXkiRJkiRJkrJk81KSJEmSJElSlmxeSpIkSZIkScqSzUtJkiRJkiRJWbJ5KUmSJEmSJClLNi8lSZIkSZIkZcnmpSRJkiRJkqQs2byUJEmSJEmSlCWbl5IkSZIkSZKyVKQOMHyFmcCuXfNXdMJpncniVNcBrAVC6iBVtFce1yZNUV078fdsTeogNXSQf76Mt8Ft22HOAPelS4DJZa3ztsqU+88493ydlSlHue8H3Qbr9qE2+HLPD5rboch1W5T6IcwGpnfNX9oJ78lxm868/qCDmC3Xfb818ODlvg1afwye2+DguA0OXqbb4Ifb4AsebJmHMBNC6DZdljpRFcuAPVOHqOFO4OLUIWo4B5iVOkQNo4EXU4fowyJgRuoQ1e362x6/y/2ZFpcY8AxgTonj1WN16gB9mA+cmDpEDTdUplydACxIHaIPmW6D4axe9h8Wchriwuwe2/QlqRNVsRzYJXWIGu4GLkodoobzgPtTh6hhIrAqdYg+LAYOTx2ihluBK1KHqOFM4JHUIfqQaf2x3kLguNQhargRuC51iBpOAualDtGHTLfBcHbPGtgCWJIkSZIkSVKWbF5KkiRJkiRJylJH6gCSNLQ99iKcfAtcfeEAFsr08HxJkiRJkvJi81KSBuX5TvjuQvjuLamTSJIkSZLUajxtXJIkSZIkSVKWbF5KkiRJkiRJypLNS0mSJEmSJElZsnkpSZIkSZIkKUs2LyVJkiRJkiRlyealJEmSJEmSpCzZvJQkSZIkSZKUJZuXkiRJkiRJkrJk81KSJEmSJElSlmxeSpIkSZIkScqSzUtJkiRJkiRJWbJ5oCJfSQAAAqFJREFUKUmSJEmSJClLNi8lSZIkSZIkZcnmpSRJkiRJkqQs2byUJEmSJEmSlCWbl5IkSZIkSZKyZPNSkiRJkiRJUpZsXkqSJEmSJEnKks1LSZIkSZIkSVmyeSlJkiRJkiQpSzYvJUmSJEmSJGXJ5qUkSZIkSZKkLNm8lCRJkiRJkpQlm5eSJEmSJEmSsmTzUpIkSZIkSVKWbF5KkiRJkiRJypLNS0mSJEmSJElZsnkpSZIkSZIkKUs2LyVJkiRJkiRlyealJEmSJEmSpCzZvJQkSZIkSZKUJZuXkiRJkiRJkrLUkTqA1tsTwjmpQ7zUx0fCSe+Gly9InaR3F2wJUwKckeG6A/jBEfDQFPhEpvlWdsCn2+D8TPMBfGI0HPNWeOV+qZP07ovbwZgpcGam6/DGA+AvE+HfMs0H8NHMt8FPjodD3wCv2SF1kt59ZXp8PCvTdXjzHnD7OPhMpvkg423w0NQBpBLsnWkNPAJOOR2mP506Se/Onwab7QenZbjuAP7rMPjbJvnWwM+NgvPbM933V5w3Go59Gxx4cOokvfvCtjBuIrwv03V4/cHw10nWwIPxyXFw+HFw9PTUSXp30c7QvhY+kOk6/OVecNd4+HSm+SDjbfCwnk8UKVIIIMwEdk2dQpIkDSntUHSmDiHVL8wGMn0jLEmScuSRl+n8CdYsg6X7wSaevi9Jkqp4oRNenPn/27l7FISBIIDCD7Q2EfuQO1jpFcSzxpPYSGqxtdLsWhl/TmClMIu87wSvXGaGhToBr+ga6Ut7uF9hWMLCN7AkSfrg9oSxhyp5eRlrAscO2m10iCRJKtXQQ7UCUnSJ9CNTOO2g2USHSJKkUl0OMF8D2W1nrAe05+gISZJUslnGwaX+ywhNoX9JSpKkMtQJyABvolBu6KD0PbAAAAAASUVORK5CYII="
+ }
+ },
+ "cell_type": "markdown",
+ "id": "6d8b79ba",
+ "metadata": {},
+ "source": [
+ "-
@@ -14,4 +14,4 @@
var myIframe = document.getElementById("notebook");
iFrameResize({log:true}, myIframe);
});
-
Settings
This document was generated with Documenter.jl version 0.27.25 on Friday 8 September 2023. Using Julia version 1.9.3.
Settings
This document was generated with Documenter.jl version 0.27.25 on Friday 8 September 2023. Using Julia version 1.9.3.
+The code in the following cell converts the general problem $Ax=b$ to the upper triangular equation system $Ux=y$. Note that this function assumes that the pivots are all nonzero. This function will be erroneos if any of the diagonal entries are zero!
-function convert_to_upper_triangular!(A,b)
- n = size(A,1)
- # Upper Triangularization: convert Ax=b to Ux=y
- for k in 1:n
- for j in k+1:n
- # Divide by pivot
- A[k,j] = A[k,j] / A[k,k]
- end
- b[k] = b[k] / A[k,k]
- A[k,k] = 1
- # Substract lower rows
- for i in k+1:n
- for j in k+1:n
- A[i,j]=A[i,j] - A[i,k] * A[k,j]
- end
- b[i] = b[i] - A[i,k] * b[k]
- A[i,k] = 0
- end
- end
- return A, b #U,y
-end
-Contents¶
In this notebook, we will learn
+-
+
- How to parallelize Gaussian elimination +
- How the data partition can create (or solve) load imbalances +
convert_to_upper_triangular! (generic function with 1 method)-
The function in the following cell solves the upper triangular equation system using backwards substitution. Note that the function alters the input values.
-function solve_upper_triangular!(U,y)
- n = size(U,1)
- for step in reverse(1:n)
- if U[step,step] == 0
- if y[step] != 0
- return "No solution"
- else
- return "Infinity solutions"
- end
- else
- # Backwards substitution
- y[step] = y[step] / U[step,step]
- end
- for row in reverse(1:step-1)
- y[row] -= U[row,step] * y[step]
- end
- end
- return y
-end
-Gaussian elimination¶
Gaussian elimination is provably one of the first algorithms you have learned to solve linear equations like this one:
+$$ +\left[ +\begin{matrix} +1 & 3 & 1 \\ +1 & 2 & -1 \\ +3 & 11 & 5 \\ +\end{matrix} +\right] +\left[ +\begin{matrix} +x \\ +y \\ +z \\ +\end{matrix} +\right] += +\left[ +\begin{matrix} +9 \\ +1 \\ +35 \\ +\end{matrix} +\right] +$$
+The value of $x$, $y$, and $z$ can be found by creating an augmented matrix and applying Gaussian elimination to it:
+$$ +\left[ +\begin{matrix} +1 & 3 & 1 & 9 \\ +1 & 2 & -1 & 1 \\ +3 & 11 & 5 & 35 \\ +\end{matrix} +\right] +\rightarrow +\left[ +\begin{matrix} +1 & 3 & 1 & 9 \\ +0 & -1 & -2 & -8 \\ +0 & 2 & 2 & 8 \\ +\end{matrix} +\right] +\rightarrow +\left[ +\begin{matrix} +1 & 3 & 1 & 9 \\ +0 & 1 & 2 & 8 \\ +0 & 0 & 1 & 4 \\ +\end{matrix} +\right] +$$
+The result is an upper diagonal matrix with ones on the diagonal. From this matrix the values $x$, $y$, and $z$ can be found via backward substitution.
solve_upper_triangular! (generic function with 1 method)-
U,y = convert_to_upper_triangular!(A,b)
-sol = solve_upper_triangular!(U,y)
-5-element Vector{Float64}:
- 1.0000000000000009
- 1.999999999999999
- 2.9999999999999964
- 4.000000000000002
- 5.000000000000005
-We can test if the obtained solution is correct using @test:
Serial implementation¶
Gaussian elimination can be implemented as shown in the following function. Note that the result is overwritten in-place on the input matrix.
using Test
-@test sol ≈ [1.0; 2.0; 3.0; 4.0; 5.0]
+function gaussian_elimination!(B)
+ n,m = size(B)
+ @inbounds for k in 1:n
+ for t in (k+1):m
+ B[k,t] = B[k,t]/B[k,k]
+ end
+ B[k,k] = 1
+ for i in (k+1):n
+ for j in (k+1):m
+ B[i,j] = B[i,j] - B[i,k]*B[k,j]
+ end
+ B[i,k] = 0
+ end
+ end
+ B
+end
Test Passed
+gaussian_elimination! (generic function with 2 methods)
Let us test the function with the example above:
+
+A = Float64[1 3 1; 1 2 -1; 3 11 5]
+b = Float64[9,1,35]
+B = [A b]
3×4 Matrix{Float64}:
+ 1.0 3.0 1.0 9.0
+ 1.0 2.0 -1.0 1.0
+ 3.0 11.0 5.0 35.0
+gaussian_elimination!(B)
+3×4 Matrix{Float64}:
+ 1.0 3.0 1.0 9.0
+ 0.0 1.0 2.0 8.0
+ 0.0 0.0 1.0 4.0
+We get the same result as shown before (as expected).
+Parallelization¶
Gaussian elimination is expensive and thus it makes sense to try to seed-up this computation by using several processors. It requires $(2n^3)/3$ operation, where $n$ is the number of rows in the input matrix. Thus, the time complexity is $O(n^3)$, which rapidly grows with $n$.
+Where can we extract parallelism?¶
n,m = size(B)
+for k in 1:n
+ for t in (k+1):m
+ B[k,t] = B[k,t]/B[k,k]
+ end
+ B[k,k] = 1
+ for i in (k+1):n
+ for j in (k+1):m
+ B[i,j] = B[i,j] - B[i,k]*B[k,j]
+ end
+ B[i,k] = 0
+ end
+end
+-
+
- The loop over k cannot be parallelized (the state of B at iteration k depends on the state at iteration k-1) +
- Loops over t, i, and j can be parallelized +
- BUT the loop over t needs to be done before loop over i and j +
Data partition¶
Let us start considering a simple 1D block partition. We assume that each process contains only a portion of the input matrix consisting in a block of consecutive rows. In this algorithm, the data stored in a process does not correspond with the data updated at a given iteration over the outer loop over k. The data updated is only a subset of the data stored, which leads to load imbalances (as we will see in a second).
+Consider next figure. Let's find out which is the data updated and the data used by CPU 3 at iteration $k$.
+
+By looking into the code above, you can see that B[i,j] is updated at iteration k, if and only if i and j are greater or equal than k. Thus, from all entries owned by CPU3, only the ones fulfilling this condition will be updated. See the entries highlighted in the next figure.
+Data dependencies¶
At iteration k, CPU3 needs part of row k to update its entries. The process owning row k needs to broadcast this entries to all other processes that require them. This communication pattern is almost identical to the one previously studied in Floyd's algorithm.
+
+Load imbalance¶
Do all processes process the same amount of data at a given iteration? This answer is no. To understand this let us find out the data updated by CPU1 in the figure below. At iteration k, CPU 1 has not data to process since k is larger than all row ids that CPU 1 owns. This is in contrast to CPU 2 and CPU 3 which have data to process. As a result, CPU 1 is idle (doing nothing), which is waisting computational resources that could be otherwise used.
+
+Addressing the load imbalance¶
For this algorithm, a cyclic distribution mitigates the load imbalance.
+-
+
- Cyclic partitions are often useful in algorithms with predictable load imbalance +
- A cyclic partition is a form of static load balancing +
- Cyclic partition are not suitable for all communication patters. E.g., algorithms that require communication between nearest-neighbors like Jacobi. +
+-
@@ -14,4 +14,4 @@
var myIframe = document.getElementById("notebook");
iFrameResize({log:true}, myIframe);
});
-
Settings
This document was generated with Documenter.jl version 0.27.25 on Friday 8 September 2023. Using Julia version 1.9.3.
Settings
This document was generated with Documenter.jl version 0.27.25 on Friday 8 September 2023. Using Julia version 1.9.3.
Getting started
Introduction
The programming of this course will be done using the Julia programming language. Thus, we start by explaining how to get up and running with Julia. After studying this page, you will be able to:
- Use the Julia REPL,
- Run serial and parallel code,
- Install and manage Julia packages.
Why Julia?
Courses related with high-performance computing (HPC) often use languages such as C, C++, or Fortran. We use Julia instead to make the course accessible to a wider set of students, including the ones that have no experience with C/C++ or Fortran, but are willing to learn parallel programming. Julia is a relatively new programming language specifically designed for scientific computing. It combines a high-level syntax close to interpreted languages like Python with the performance of compiled languages like C, C++, or Fortran. Thus, Julia will allow us to write efficient parallel algorithms with a syntax that is convenient in a teaching setting. In addition, Julia provides easy access to different programming models to write distributed algorithms, which will be useful to learn and experiment with them.
You can run the code in this link to learn how Julia compares to other languages (C and Python) in terms of performance.
Installing Julia
This is a tutorial-like page. Follow these steps before you continue reading the document.
- Download and install Julia from julialang.org;
- Follow the specific instructions for your operating system: Windows, MacOS, or Linux
- Download and install VSCode and its Julia extension;
The Julia REPL
Starting Julia
There are several ways of opening Julia depending on your operating system and your IDE, but it is usually as simple as launching the Julia app. With VSCode, open a folder (File > Open Folder). Then, press Ctrl+Shift+P to open the command bar, and execute Julia: Start REPL. If this does not work, make sure you have the Julia extension for VSCode installed. Independently of the method you use, opening Julia results in a window with some text ending with:
julia>You have just opened the Julia read-evaluate-print loop, or simply the Julia REPL. Congrats! You will spend most of time using the REPL, when working in Julia. The REPL is a console waiting for user input. Just as in other consoles, the string of text right before the input area (julia> in the case) is called the command prompt or simply the prompt.
Basic usage
The usage of the REPL is as follows:
- You write some input
- press enter
- you get the output
For instance, try this
julia> 1 + 1A "Hello world" example looks like this in Julia
julia> println("Hello, world!")Try to run it in the REPL.
Help mode
Curious about what the function println does? Enter into help mode to look into the documentation. This is done by typing a question mark (?) into the input field:
julia> ?After typing ?, the command prompt changes to help?>. It means we are in help mode. Now, we can type a function name to see its documentation.
help?> printlnPackage and shell modes
The REPL comes with two more modes, namely package and shell modes. To enter package mode type
julia> ]Package mode is used to install and manage packages. We are going to discuss the package mode in greater detail later. To return back to normal mode press the backspace key several times.
To enter shell mode type semicolon (;)
julia> ;The prompt should have changed to shell> indicating that we are in shell mode. Now you can type commands that you would normally do on your system command line. For instance,
shell> lswill display the contents of the current folder in Mac or Linux. Using shell mode in Windows is not straightforward, and thus not recommended for beginners.
Running Julia code
Running more complex code
Real-world Julia programs are not typed in the REPL in practice. They are written in one or more files and included in the REPL. To try this, create a new file called hello.jl, write the code of the "Hello world" example above, and save it. If you are using VSCode, you can create the file using File > New File > Julia File. Once the file is saved with the name hello.jl, execute it as follows
julia> include("hello.jl")Make sure that the file "hello.jl" is located in the current working directory of your Julia session. You can query the current directory with function pwd(). You can change to another directory with function cd() if needed. Also, make sure that the file extension is .jl.
The recommended way of running Julia code is using the REPL as we did. But it is also possible to run code directly from the system command line. To this end, open a terminal and call Julia followed by the path to the file containing the code you want to execute.
$ julia hello.jlThe previous line assumes that you have Julia properly installed in the system and that it's usable from the terminal. In UNIX systems (Linux and Mac), the Julia binary needs to be in one of the directories listed in the PATH environment variable. To check that Julia is properly installed, you can use
$ julia --versionIf this runs without error and you see a version number, you are good to go!
In this tutorial, when a code snipped starts with $, it should be run in the terminal. Otherwise, the code is to be run in the Julia REPL.
Avoid calling Julia code from the terminal, use the Julia REPL instead! Each time you call Julia from the terminal, you start a fresh Julia session and Julia will need to compile your code from scratch. This can be time consuming for large projects. In contrast, if you execute code in the REPL, Julia will compile code incrementally, which is much faster. Running code in a cluster (like in DAS-5 for the Julia assignment) is among the few situations you need to run Julia code from the terminal.
Running parallel code
Since we are in a parallel computing course, let's run a parallel "Hello world" example in Julia. Open a Julia REPL and write
julia> using Distributed
+Getting started · XM_40017 Getting started
Introduction
The programming of this course will be done using the Julia programming language. Thus, we start by explaining how to get up and running with Julia. After studying this page, you will be able to:
- Use the Julia REPL,
- Run serial and parallel code,
- Install and manage Julia packages.
Why Julia?
Courses related with high-performance computing (HPC) often use languages such as C, C++, or Fortran. We use Julia instead to make the course accessible to a wider set of students, including the ones that have no experience with C/C++ or Fortran, but are willing to learn parallel programming. Julia is a relatively new programming language specifically designed for scientific computing. It combines a high-level syntax close to interpreted languages like Python with the performance of compiled languages like C, C++, or Fortran. Thus, Julia will allow us to write efficient parallel algorithms with a syntax that is convenient in a teaching setting. In addition, Julia provides easy access to different programming models to write distributed algorithms, which will be useful to learn and experiment with them.
Tip You can run the code in this link to learn how Julia compares to other languages (C and Python) in terms of performance.
Installing Julia
This is a tutorial-like page. Follow these steps before you continue reading the document.
- Download and install Julia from julialang.org;
- Follow the specific instructions for your operating system: Windows, MacOS, or Linux
- Download and install VSCode and its Julia extension;
The Julia REPL
Starting Julia
There are several ways of opening Julia depending on your operating system and your IDE, but it is usually as simple as launching the Julia app. With VSCode, open a folder (File > Open Folder). Then, press Ctrl+Shift+P to open the command bar, and execute Julia: Start REPL. If this does not work, make sure you have the Julia extension for VSCode installed. Independently of the method you use, opening Julia results in a window with some text ending with:
julia>
You have just opened the Julia read-evaluate-print loop, or simply the Julia REPL. Congrats! You will spend most of time using the REPL, when working in Julia. The REPL is a console waiting for user input. Just as in other consoles, the string of text right before the input area (julia> in the case) is called the command prompt or simply the prompt.
Basic usage
The usage of the REPL is as follows:
- You write some input
- press enter
- you get the output
For instance, try this
julia> 1 + 1
A "Hello world" example looks like this in Julia
julia> println("Hello, world!")
Try to run it in the REPL.
Help mode
Curious about what the function println does? Enter into help mode to look into the documentation. This is done by typing a question mark (?) into the input field:
julia> ?
After typing ?, the command prompt changes to help?>. It means we are in help mode. Now, we can type a function name to see its documentation.
help?> println
Package and shell modes
The REPL comes with two more modes, namely package and shell modes. To enter package mode type
julia> ]
Package mode is used to install and manage packages. We are going to discuss the package mode in greater detail later. To return back to normal mode press the backspace key several times.
To enter shell mode type semicolon (;)
julia> ;
The prompt should have changed to shell> indicating that we are in shell mode. Now you can type commands that you would normally do on your system command line. For instance,
shell> ls
will display the contents of the current folder in Mac or Linux. Using shell mode in Windows is not straightforward, and thus not recommended for beginners.
Running Julia code
Running more complex code
Real-world Julia programs are not typed in the REPL in practice. They are written in one or more files and included in the REPL. To try this, create a new file called hello.jl, write the code of the "Hello world" example above, and save it. If you are using VSCode, you can create the file using File > New File > Julia File. Once the file is saved with the name hello.jl, execute it as follows
julia> include("hello.jl")
Warning Make sure that the file "hello.jl" is located in the current working directory of your Julia session. You can query the current directory with function pwd(). You can change to another directory with function cd() if needed. Also, make sure that the file extension is .jl.
The recommended way of running Julia code is using the REPL as we did. But it is also possible to run code directly from the system command line. To this end, open a terminal and call Julia followed by the path to the file containing the code you want to execute.
$ julia hello.jl
The previous line assumes that you have Julia properly installed in the system and that it's usable from the terminal. In UNIX systems (Linux and Mac), the Julia binary needs to be in one of the directories listed in the PATH environment variable. To check that Julia is properly installed, you can use
$ julia --version
If this runs without error and you see a version number, you are good to go!
Note In this tutorial, when a code snipped starts with $, it should be run in the terminal. Otherwise, the code is to be run in the Julia REPL.
Tip Avoid calling Julia code from the terminal, use the Julia REPL instead! Each time you call Julia from the terminal, you start a fresh Julia session and Julia will need to compile your code from scratch. This can be time consuming for large projects. In contrast, if you execute code in the REPL, Julia will compile code incrementally, which is much faster. Running code in a cluster (like in DAS-5 for the Julia assignment) is among the few situations you need to run Julia code from the terminal.
Running parallel code
Since we are in a parallel computing course, let's run a parallel "Hello world" example in Julia. Open a Julia REPL and write
julia> using Distributed
julia> @everywhere println("Hello, world! I am proc $(myid()) from $(nprocs())")
Here, we are using the Distributed package, which is part of the Julia standard library that provides distributed memory parallel support. The code prints the process id and the number of processes in the current Julia session.
You will probably only see output from 1 process. We need to add more processes to run the example in parallel. This is done with the addprocs function.
julia> addprocs(3)
We have added 3 new processes. Plus the old one, we have 4 processes. Run the code again.
julia> @everywhere println("Hello, world! I am proc $(myid()) from $(nprocs())")
Now, you should see output from 4 processes.
It is possible to specify the number of processes when starting Julia from the terminal with the -p argument (useful, e.g., when running in a cluster). If you launch Julia from the terminal as
$ julia -p 3
and then run
julia> @everywhere println("Hello, world! I am proc $(myid()) from $(nprocs())")
You should get output from 4 processes as before.
Installing packages
One of the most useful features of Julia is its package manager. It allows one to install Julia packages in a straightforward and platform independent way. To illustrate this, let us consider the following parallel "Hello world" example. This example uses the Message Passing Interface (MPI). We will learn more about MPI later in the course.
Copy the following block of code into a new file named "hello_mpi.jl"
# file hello_mpi.jl
using MPI
MPI.Init()
@@ -14,4 +14,4 @@ julia> DataFrame(a=[1,2],b=[3,4])
You should get an error or a
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
Copy the contents of previous code block into a file called Project.toml and place it in an empty folder named newproject. It is important that the file is named Project.toml. You can create a new folder from the REPL with
julia> mkdir("newproject")
To install all the packages registered in this file you need to activate the folder containing your Project.toml file
(@v1.8) pkg> activate newproject
and then instantiating it
(newproject) pkg> instantiate
The instantiate command will download and install all listed packages and their dependencies in just one click.
Getting help in package mode
You can get help about a particular package operator by writing help in front of it
(@v1.8) pkg> help activate
You can get an overview of all package commands by typing help alone
(@v1.8) pkg> help
Package operations in Julia code
In some situations it is required to use package commands in Julia code, e.g., to automatize installation and deployment of Julia applications. This can be done using the Pkg package. For instance
julia> using Pkg
-julia> Pkg.status()
is equivalent to calling status in package mode.
(@v1.8) pkg> status
Conclusion
We have learned the basics of how to work with Julia. If you want to further dig into the topics we have covered here, you can take a look at the following links:
Settings
This document was generated with Documenter.jl version 0.27.25 on Friday 8 September 2023. Using Julia version 1.9.3.
is equivalent to calling status in package mode.
(@v1.8) pkg> statusConclusion
We have learned the basics of how to work with Julia. If you want to further dig into the topics we have covered here, you can take a look at the following links:
Settings
This document was generated with Documenter.jl version 0.27.25 on Friday 8 September 2023. Using Julia version 1.9.3.
Programming Large-Scale Parallel Systems (XM_40017)
Welcome to the interactive lecture notes of the Programming Large-Scale Parallel Systems course at VU Amsterdam!
What
This page contains part of the course material of the Programming Large-Scale Parallel Systems course at VU Amsterdam. We provide several lecture notes in jupyter notebook format, which will help you to learn how to design, analyze, and program parallel algorithms on multi-node computing systems. Further information about the course is found in the study guide (click here) and our Canvas page (for registered students).
Material will be added incrementally to the website as the course advances.
This page will eventually contain only a part of the course material. The rest will be available on Canvas. In particular, the material in this public webpage does not fully cover all topics in the final exam.
How to use this page
You have two main ways of studying the notebooks:
- Download the notebooks and run them locally on your computer (recommended). At each notebook page you will find a green box with links to download the notebook.
- You also have the static version of the notebooks displayed in this webpage for quick reference.
How to run the notebooks locally
To run a notebook locally follow these steps:
- Install Julia (if not done already). More information in Getting started.
- Download the notebook.
- Launch Julia. More information in Getting started.
- Execute these commands in the Julia command line:
julia> using Pkg
+Home · XM_40017 Programming Large-Scale Parallel Systems (XM_40017)
Welcome to the interactive lecture notes of the Programming Large-Scale Parallel Systems course at VU Amsterdam!
What
This page contains part of the course material of the Programming Large-Scale Parallel Systems course at VU Amsterdam. We provide several lecture notes in jupyter notebook format, which will help you to learn how to design, analyze, and program parallel algorithms on multi-node computing systems. Further information about the course is found in the study guide (click here) and our Canvas page (for registered students).
Note Material will be added incrementally to the website as the course advances.
Warning This page will eventually contain only a part of the course material. The rest will be available on Canvas. In particular, the material in this public webpage does not fully cover all topics in the final exam.
How to use this page
You have two main ways of studying the notebooks:
- Download the notebooks and run them locally on your computer (recommended). At each notebook page you will find a green box with links to download the notebook.
- You also have the static version of the notebooks displayed in this webpage for quick reference.
How to run the notebooks locally
To run a notebook locally follow these steps:
- Install Julia (if not done already). More information in Getting started.
- Download the notebook.
- Launch Julia. More information in Getting started.
- Execute these commands in the Julia command line:
julia> using Pkg
julia> Pkg.add("IJulia")
julia> using IJulia
-julia> notebook()
- These commands will open a jupyter in your web browser. Navigate in jupyter to the notebook file you have downloaded and open it.
Authors
This material is created by Francesc Verdugo with the help of Gelieza Kötterheinrich. Part of the notebooks are based on the course slides by Henri Bal.
License
All material on this page that is original to this course may be used under a CC BY 4.0 license.
Acknowledgment
This page was created with the support of the Faculty of Science of Vrije Universiteit Amsterdam in the framework of the project "Interactive lecture notes and exercises for the Programming Large-Scale Parallel Systems course" funded by the "Innovation budget BETA 2023 Studievoorschotmiddelen (SVM) towards Activated Blended Learning".
Settings
This document was generated with Documenter.jl version 0.27.25 on Friday 8 September 2023. Using Julia version 1.9.3.
- These commands will open a jupyter in your web browser. Navigate in jupyter to the notebook file you have downloaded and open it.
Authors
This material is created by Francesc Verdugo with the help of Gelieza Kötterheinrich. Part of the notebooks are based on the course slides by Henri Bal.
License
All material on this page that is original to this course may be used under a CC BY 4.0 license.
Acknowledgment
This page was created with the support of the Faculty of Science of Vrije Universiteit Amsterdam in the framework of the project "Interactive lecture notes and exercises for the Programming Large-Scale Parallel Systems course" funded by the "Innovation budget BETA 2023 Studievoorschotmiddelen (SVM) towards Activated Blended Learning".
Settings
This document was generated with Documenter.jl version 0.27.25 on Friday 8 September 2023. Using Julia version 1.9.3.
-
@@ -14,4 +14,4 @@
var myIframe = document.getElementById("notebook");
iFrameResize({log:true}, myIframe);
});
-
Settings
This document was generated with Documenter.jl version 0.27.25 on Friday 8 September 2023. Using Julia version 1.9.3.
Settings
This document was generated with Documenter.jl version 0.27.25 on Friday 8 September 2023. Using Julia version 1.9.3.
-
@@ -14,4 +14,4 @@
var myIframe = document.getElementById("notebook");
iFrameResize({log:true}, myIframe);
});
-
Settings
This document was generated with Documenter.jl version 0.27.25 on Friday 8 September 2023. Using Julia version 1.9.3.
Settings
This document was generated with Documenter.jl version 0.27.25 on Friday 8 September 2023. Using Julia version 1.9.3.
-
@@ -14,4 +14,4 @@
var myIframe = document.getElementById("notebook");
iFrameResize({log:true}, myIframe);
});
-
Settings
This document was generated with Documenter.jl version 0.27.25 on Friday 8 September 2023. Using Julia version 1.9.3.
Settings
This document was generated with Documenter.jl version 0.27.25 on Friday 8 September 2023. Using Julia version 1.9.3.
Tasks¶
Creating a task¶
A task is a piece of computation work that can be run asynchronously (i.e., that can be run in the background). To create a task, we first need to create a function that represents the work to be done in the task. In next cell, we generate a task that generates and sums two matrices.
+Tasks¶
Creating a task¶
Technically, a task in Julia is a symmetric co-routine. More informally, a task is a piece of computation work that can be started (scheduled) at some point in the future, and that can be interrupted and resumed. To create a task, we first need to create a function that represents the work to be done in the task. In next cell, we generate a task that generates and sums two matrices.
Tasks do not run in parallel¶
It is also important to note that tasks do not run in parallel. We were able to run code while previous tasks was running because the task was idling most of the time. If the task does actual work, the current process will be busy running this task and it is likely that we cannot run other code at the same time. Let's illustrate this with an example. The following code computes an approximation of $\pi$ using Leibniz formula. The quality of the approximation increases with the value of n.
Tasks do not run in parallel¶
It is also important to note that tasks do not run in parallel. We were able to run code while previous tasks was running because the task was idling most of the time in the sleep function. If the task does actual work, the current process will be busy running this task and preventing to run other tasks. Let's illustrate this with an example. The following code computes an approximation of $\pi$ using Leibniz formula. The quality of the approximation increases with the value of n.
yield¶
If tasks do not run in parallel, what is the purpose of tasks? Tasks are handy since they can be interrupted and to switch control to other tasks. This is achieved via function yield. When we call yield, we provide the opportunity to switch to another task. The function below is a variation of function compute_π in which we yield every 1000 iterations. At the call to yield we allow other tasks to take over. Without this call to yield, once we start function compute_π we cannot start any other tasks until this function finishes.
function compute_π_yield(n)
+ s = 1.0
+ for i in 1:n
+ s += (isodd(i) ? -1 : 1) / (i*2+1)
+ if mod(i,1000) == 0
+ yield()
+ end
+ end
+ 4*s
+end
+You can check this behavior experimentally with the two following cells. The next one creates and schedules a task that computes pi with the function compute_π_yield. Note that you can run the 2nd cell bellow while this task is running since we call to yield often inside compute_π_yield.
fun = () -> compute_π_yield(3_000_000_000)
+t = Task(fun)
+schedule(t)
+1+1
+Example: Implementing function sleep¶
Using yield, we can implement our own version of the sleep function as follows:
+function mysleep(secs)
+ final_time = time() + secs
+ while time() < final_time
+ yield()
+ end
+ nothing
+end
+You can check that it behaves as expected.
+@time mysleep(3)
+a) 10*t
b) t
@@ -8274,7 +8405,7 @@ d) near 0*t a) 10*t
b) t
@@ -8307,7 +8438,7 @@ d) near 0*t a) 10*t
b) t
@@ -8340,7 +8471,7 @@ d) near 0*t a) infinity
b) 1 second
@@ -8390,7 +8521,7 @@ d) 3 secondsa) infinity
b) 1 second
diff --git a/dev/julia_basics.ipynb b/dev/julia_basics.ipynb
index 781b196..ff73688 100644
--- a/dev/julia_basics.ipynb
+++ b/dev/julia_basics.ipynb
@@ -449,7 +449,7 @@
"source": [
"\n",
"\n",
- "Question: What will be the value of `x` in the last line ? (Think your answer before executing next cell to find out the result) \n",
+ "Question (NB1-Q1): What will be the value of `x` in the last line ? (Think your answer before executing next cell to find out the result) \n",
"\n",
"\n"
]
@@ -557,7 +557,7 @@
"metadata": {},
"source": [
"\n",
- "Question: What will be the value of `x` in the last line ?\n",
+ "Question (NB1-Q2): What will be the value of `x` in the last line ?\n",
"\n"
]
},
@@ -1030,7 +1030,7 @@
"metadata": {},
"source": [
"\n",
- "Question: Which will be the value of `x` below? \n",
+ "Question (NB1-Q3): Which will be the value of `x` below? \n",
"\n",
"\n"
]
@@ -1645,14 +1645,6 @@
"\n",
"This notebook is part of the course [Programming Large Scale Parallel Systems](https://www.francescverdugo.com/XM_40017/) at Vrije Universiteit Amsterdam and may be used under a [CC BY 4.0](https://creativecommons.org/licenses/by/4.0/) license."
]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "f8d92f25",
- "metadata": {},
- "outputs": [],
- "source": []
}
],
"metadata": {
diff --git a/dev/julia_basics/index.html b/dev/julia_basics/index.html
index 4b9c0ab..169bd33 100644
--- a/dev/julia_basics/index.html
+++ b/dev/julia_basics/index.html
@@ -1,5 +1,5 @@
-Julia Basics · XM_40017
+Julia Basics · XM_40017
Tip
@@ -14,4 +14,4 @@
var myIframe = document.getElementById("notebook");
iFrameResize({log:true}, myIframe);
});
-
Settings
This document was generated with Documenter.jl version 0.27.25 on Friday 8 September 2023. Using Julia version 1.9.3.
Settings
This document was generated with Documenter.jl version 0.27.25 on Friday 8 September 2023. Using Julia version 1.9.3.
-
-