{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 11: From dense to sparse linear algebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Syllabus\n",
    "**Week 1:** Matrices, vectors, matrix/vector norms, scalar products & unitary matrices  \n",
    "**Week 2:** TAs-week (Strassen, FFT, a bit of SVD)  \n",
    "**Week 3:** Matrix ranks, singular value decomposition, linear systems, eigenvalues  \n",
    "**Week 4:** Matrix decompositions: QR, LU, SVD + test + structured matrices start"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Recap of the previous lecture\n",
    "- Algorithms for the symmetric eigenvalue problems (QR-algorithm, Divide-and-Conquer, bisection)\n",
    "- SVD and its applications (collaborative filtering, integral equations, latent semantic indexing)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Today lecture\n",
    "\n",
    "Today we will talk about **sparse matrices**, where they arise, how we store them, how we operate with them.\n",
    "\n",
    "\n",
    "- Formats: list of lists and compressed sparse row format, relation to graphs\n",
    "- Fast matrix-by-vector product\n",
    "- Fast direct solvers for Gaussian elimination\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Sparse matrices are ubiqitous in PDE\n",
    "Consider the simplest partial differential equation (PDE), called **Laplace equation**:  \n",
    "$$\n",
    "   \\Delta T = \\frac{\\partial^2 T}{\\partial x^2} + \\frac{\\partial^2 T}{\\partial y^2} = f.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Discretization\n",
    "$$\\frac{\\partial^2 T}{\\partial x^2} \\approx \\frac{T(x+h) + T(x-h) - 2T(x)}{h^2} + \\mathcal{O}(h^4),$$\n",
    "same for $\\frac{\\partial^2 T}{\\partial y^2},$\n",
    "and we get a linear system.  \n",
    "    First, let us do **one-dimensional case**:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "After the discretization of the one-dimensional Laplace equation with Dirichlet boundary conditions we have\n",
    "$$\\frac{u_{i+1} + u_{i-1} - 2u_i}{h^2} = f_i,$$\n",
    "\n",
    "or in the matrix form\n",
    "\n",
    "$$\n",
    "   A u = f,\n",
    "$$\n",
    "                   and ($n = 5$ illustration)\n",
    "\\begin{equation}\n",
    "A=-\\frac{1}{h^2}\\begin{bmatrix}\n",
    "   2& -1 & 0 & 0 & 0\\\\\n",
    "   -1 & 2 & -1 & 0 &0 \\\\\n",
    "   0 & -1 & 2& -1 & 0 \\\\\n",
    "   0 & 0 & -1 & 2  &-1\\\\\n",
    "   0 & 0 & 0 & -1 & 2\n",
    "   \\end{bmatrix}\n",
    "   \\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The matrix is **triadiagonal** and **sparse**  \n",
    "(and also **Toeplitz**: all elements on the diagonal are the same)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Block structure in 2D\n",
    "In two dimensions, we get equation of the form  \n",
    "$$\\frac{4u_{ij} -u_{(i-1)j} - u_{(i+1)j} - u_{i(j-1)}-u_{i(j+1)}}{h^2} = f_{ij},$$\n",
    "\n",
    "or in the **Kronecker product form**  \n",
    "\n",
    "$$\\Delta_2 = \\Delta_1 \\otimes I + I \\otimes \\Delta_1,$$\n",
    "\n",
    "where $\\Delta_1$ is a 1D Laplacian, and $\\otimes$ is a **Kronecker product** of matrices. \n",
    "\n",
    "For matrices $A$ and $B$ its Kronecker product is defined as a block matrix of the form \n",
    "$$\n",
    "   [a_{ij} B]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "In the block matrix form the 2D-Laplace matrix can be written in the following form:  \n",
    "\n",
    "$$\n",
    "   A = -\\frac{1}{h^2}\\begin{bmatrix}\n",
    "   \\Delta_1 + 2I & -I & 0 & 0 & 0\\\\\n",
    "   -I & \\Delta_1 + 2I  & -I & 0 &0 \\\\\n",
    "   0 & -I & \\Delta_1 + 2I & -I & 0 \\\\\n",
    "   0 & 0 & -I & \\Delta_1 + 2I   &-I\\\\\n",
    "   0 & 0 & 0 & -I & \\Delta_1 + 2I \n",
    "   \\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We can create this matrix using **scipy.sparse** package (actually this is **not the best** sparse matrix package)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x10eccc390>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQwAAAEACAYAAABGTkjoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF3pJREFUeJzt3X+QXWV9x/H3h5+K0GIGJqDBkolBpEPJTpUpynZ3O5WC\nwwCC/LDFiYU6ndGQuE4tP2ZKYvwDpMWVNZXp1DiDRKkgkQY7CEE3mE75UdskIBFJaDJjLCxRiSVj\nS6F8+8c5N3uz7P19zr3n3Pt5zdzZc8+ee549s9kn3+d7nvN8FRGYmTXjkF7/AGZWHu4wzKxp7jDM\nrGnuMMysae4wzKxp7jDMrGld7zAknSvpGUk7JF2bUxu7JT0paYukJ9J98yRtlPSspIckHdvB+b8q\naVrSU1X7ap5f0vXp9T4j6ZwM2lolaU96fVsknZdFW+nnT5I0JelpST+StDzn66vVXi7XKOlNkh6X\ntFXSdkk35Xx9tdrL7XeYnuPQ9Lz3Z3p9EdG1F3AosBM4GTgc2Aq8O4d2dgHzZu27BfjLdPta4OYO\nzj8MDAFPNTo/cFp6nYen170TOKTDtlYCn57j2I7aSs9xArAk3T4a+Anw7hyvr1Z7eV7jUenXw4DH\ngLPzur467eV2fel5Pg18HdiQ5b/PbkcYZwI7I2J3RLwK/ANwYU5tadb7C4A70u07gIvaPXFEbAZe\navL8FwJ3RcSrEbGb5BdyZodtwRuvr+O20vZeiIit6fZ+4MfA28nv+mq1B/ld46/TzSNI/hN7iZyu\nr057kNP1SVoAfBD4SlUbmVxftzuMtwM/rXq/h5l/HFkK4GFJP5T08XTf/IiYTrengfkZt1nr/G8j\nuc6KrK75GknbJK2tCi8zbUvSySTRzeN04fqq2nss3ZXLNUo6RNJWkuuYioinyfH6arQH+f0OJ4DP\nAK9X7cvk+rrdYXRrHvr7I2IIOA/4pKThg36IJBbL7Wdp4vydtn07sBBYAjwP3Jp1W5KOBu4FVkTE\nywedMIfrS9v7VtrefnK8xoh4PSKWAAuA35c0Nuv7mV7fHO2NktP1STofeDEitjB3BNPR9XW7w/gZ\ncFLV+5M4uHfLREQ8n37dC3ybJMSalnQCgKQTgRczbrbW+Wdf84J0X9si4sVIkYSdlRAyk7YkHU7S\nWdwZEfelu3O7vqr21lXay/sa0zZ+BfwT8Lt04fdX1d57cry+9wEXSNoF3AX8gaQ7yej6ut1h/BBY\nLOlkSUcAlwMbsmxA0lGSjkm33wKcAzyVtrM0PWwpcN/cZ2hbrfNvAK6QdISkhcBi4IlOGkp/4RUf\nIrm+TNqSJGAtsD0ivlj1rVyur1Z7eV2jpOMq4b+kNwMfALbkeH1ztlf54836+iLihog4KSIWAlcA\n34+Ij2Z2fa1mXzt9kQwTfkKSXLk+h/MvJMn6bgV+VGkDmAc8DDwLPAQc20EbdwH/CfwvSU7mT+ud\nH7ghvd5ngD/qsK2rgK8BTwLb0l/8/CzaSj9/NsnYdyvJH9IW4Nwcr2+u9s7L6xqB04F/T9t7EvhM\no38fObWX2++w6jwjzNwlyeT6lH7AzKyhXIYk6sLkLDPrvswjDEmHkgw5/pAkefKvwEci4seZNmRm\nXZdHhNHNyVlm1kV5dBjdmpxlZl12WA7nbDjGkeRMq1mPRMScE7qakUeE0dTkrGOOWcn8+Sv52MdW\nMjU1VfM2zpe+FHz4w8nXdm8vrVy5suPbtW7P7ZWxvampKVauXHng1ak8IowDk7NI5g9cDnxk9kGL\nFq1i1y7YuhWuuqr2yZYtS74+8sjB782ssdHRUUZHRw+8/+xnP9vR+TKPMCLiNWAZ8CCwHfhmzHGH\nZHISFi6EvXth/fr651y2DEZG4LnnYPPmrH9iM2tWLvMwIuKBiHhXRLwzIm6a65jh4aTTOP54uOce\nGB+vf85ly2DRoqRzabXTqO5hu8Htub0it9eJnsz0lBSVdsfHYd06OPJIuPRSuPjipDOZy+bNsG1b\nEmksWuThiVmrJBEdJD173mFA0mls2pQMT44/Pok8anUaAGvWJDmNkRF3Gmat6IsOA5LoYfly2LUr\nyW240zDLXqcdRmFWDa/kNJwINSuuwkQYFZVIY+/eJKcxMVH/XJs3J52LcxpmjfVNhFExPAyjo/DK\nK8ndkzVrGh+/aFEyPGl0rJl1pmcRxg9+EHVzFOPjyaSuffvg2GNh9WrnNMw6Vdqk56c+FU0NI8bG\nYMsWJ0LNslDaIcmiRc0lLFevdiLUrCh6mvRcsyb54643WQucCDXLSmkjDIAzzki+NpruPTsROj7e\n+HgnQs2y1/PbqpVoABpHGp4RataZ0iY9q9tt5RkRzwg1a1+phyQVw8MzT6M2Slh6RqhZ7xQiwqjW\naiJ0/344//zmjnci1AZdX0QY1VpNhB59dJIIXb7ciVCzvBUuwoDWEqHVOY2hIZiaqt+2cxo2yPoi\n6TmXVhOhN94IO3bA4sWeRm5WS98NSSpaTYROTSWdxY4dToSa5aWjCEPSbuC/gP8DXo2IMyXNA74J\n/BawG7gsIvbN+lzDCKOaE6Fm2eh1hBHAaEQMRcSZ6b7rgI0RcQrwvfR9R5wINSuGTiOMXcB7IuIX\nVfueAUYiYlrSCcCmiDh11udaijCg/USoJ3eZzehp0lPSfwC/IhmS/F1E/L2klyLiren3Bfyy8r7q\ncy13GOAZoWad6rTD6LTy2fsj4nlJxwMb0+jigIiIWnVUV61adWB7dnWmWoaHk1clp7F5c+1OoDIj\ntPKU6/r19TuMSifR6LxmZbJp0yY2bdqU2fkyu60qaSWwH/g4SV7jBUknAlNZDElmcyLUrHU9S3pK\nOkrSMen2W4BzgKeADcDS9LClwH3ttlGPE6Fm3dd2hCFpIfDt9O1hwNcj4qb0turdwDvI6LZqLc5p\nmLWmb2d6tiKvlbuaPa9ZWfR6HkYhNFuouZIIXbw4WYin0ZCjkwLQZv2oLzoMaC2nccklSRL0ttsa\nL/d3xhkznYZzGjbo+mJIUuHl/szqcw5jFidCzWpzDmMWL/dnlp++izCque6J2cEcYdThuidm2err\nCAOcCDWr5qRnE5wINUt4SNIEJ0LNsjEQEUY1J0JtkDnCaFGriVCYmUbuRKgNuoGLMKC1RGglR7F1\na/KIvHMaVmZOerbJiVAbRB6StMmJULPWDWyEUc3L/dmgcISRAS/3Z9YcRxgp1z2xQeCkZ4acCLV+\nl/uQRNJXJU1Leqpq3zxJGyU9K+khScdWfe96STskPSPpnHZ/sF5wItSsvoYRhqRhknojX4uI09N9\ntwA/j4hbJF0LvDUirpN0GvAN4L3A24GHgVMi4vVZ5yxkhFHNiVDrR7lHGBGxGXhp1u4LgDvS7TuA\ni9LtC4G7IuLViNgN7ATOpIScCDV7o3ZLJc6PiOl0exqYn26/DXis6rg9JJFG6VSihPXrZ4YbtSKH\niYmDnz2ZnKx/fCWyeOSRg9+bFV2ntVXr1k+tHDLXznZqq3Zb5Q9+27bkVb1vrmMnJ5PXzp2u5WrF\n0JPaqpJOBu6vymE8wxz1UyVdBxARN6fHfRdYGRGPzzpf4XMYs7lYkvWDXk3cqlU/dQNwhaQj0lKK\ni4En2v3hiiSv5f5c98TKpJm7JHcBI8BxJPmKG4F/pEb9VEk3AFcBrwErIuLBOc5ZuggDvNyflZ8n\nbnWZJ3dZmflZki7z5C4bZI4wOuDl/qxsHGH0kOue2KBxhNEhJ0KtTJz0LAAnQq0sPCQpgHYTofv3\nz8wgrcWJUCsSRxgZayURWplGPjrqRKh1hyOMgmklEToykuQznAi1snCEkQMnQq2onPQsKCdCrYg8\nJCkozwi1fuQIowu83J8VhSOMEvByf9YvHGF0ieueWBE46VkiToRar3lIUiJOhFrZOcLokXYSoStW\nNI4enAi1ehxhlFSridAFC2DtWhgbcyLUescRRg+1kgiFpLPYssU5DWtfr2qrrpK0R9KW9HVe1fdK\nW1u124aHk45i0aIkGdoo97B6tXMa1lvt1lZdCbwcEV+YdWzf1FbtNtc9sW7oVW1VgLka7Zvaqt3m\nuidWBp0kPa+RtE3SWknHpvveRlJPtaK0tVW7rTI8gcadxsQEXHll8nRroxmh1bdynQi1TrVbW/V2\nYHW6/TngVuDqGseWtrZqt1XXcl2/PvlaK2FZXQB6167ka71EqAtAD6ZC1Fat9b1+rq3abc5pWB56\nMg8jLcBc8SGgcgelb2urdltlGNFMTmNycmZ4Mj7e+LwXX+ychrWnmduqdwH/ArxL0k8lXQV8XtKT\nkraR1F0dB4iI7SQ1V7cDDwCfcCjRPtc9saLxxK2Ca3e5v8WLYWqq/rk9uWvw+GnVAdDqU6433gg7\ndiSdxurVnhFqM/wsyQBo9SnXqamks9ixwzNCLVuOMErGBaCtE44wBowTodZLjjBKyHVPrF1Oeg4o\nL/dn7fCQZEB5uT/rBUcYfcB1T6xZjjDMdU+saxxh9AnXPbFmOOlpB7SbCB0a8jTyQeEhiR3QTiJ0\naAj27Ws85HAi1MARRt9qZd2LsTHYs8eJ0EHgCMPm1GwiFGDJEidCrTmOMPqYE6E2myMMq6mVuiee\n3GXNcIQxILxGqIEjDGuS655YFhxhDBA/5WrdqK16kqQpSU9L+pGk5en+eZI2SnpW0kNVxYxcX7Wg\nqnMajSKCiYmZ1cgrdU/qRRqVnIbvnvS3ZoYkrwLjEfHbwO8Bn5T0buA6YGNEnAJ8L31fqa96OXAa\ncC7wZUke+hSEn3K1TjRTW/WFiNiabu8HfkxS/vAC4I70sDuAi9Jt11ctAdc9sXa09D9/WuVsCHgc\nmB8R0+m3poH56bbrq5ZEq4lQSPIaa9Z4ctegarq2qqSjgXuBFRHxsjSTN4mIkFQvi/mG77m2au9V\nkpjr188MN2olNicmZjqB225LZoa6lmvx9aq26uHAd4AHIuKL6b5ngNGIeCEtnTgVEac2U1/Vd0mK\nxcv9DY5u3CURsBbYXuksUhuApen2UuC+qv2ur1oiToRasxpGGJLOBn4APMnM0OJ6kk7gbuAdwG7g\nsojYl37mBuAq4DWSIcyDs87pCKOgXPekv3kBHctUq5O71q2DI49MOo1Gx3t40nvuMCxznhHav9xh\nWC6cCO1PfvjMcuFEqM3FEYY15Lon/cMRhuXOdU+soumZnja4Wp0RWn3LdXKy/vGeEVouHpJY01pN\nhE5OwqOP+u5JkXhIYl3TaiJ0+fKks3AitH84wrC2OBFaTo4wrCecCB1MjjCsba57Uj6e6Wk95Rmh\n5eIhifWUZ4QOFkcYlhknQovPEYYVhhOh/c8RhmXKOY1ic4RhhdJuTmP//qSjqcc5jd5zhGG5aSWn\nMTkJO3cmQxUXgM6PIwwrrFZyGiMjyZ0TF4Autk5qq66StEfSlvR1XtVnXFvVDtRyhcadxrJlybqg\nlQpr9RKh1cMeJ0K7q5lVw08AToiIrWkxo38jKYt4GfByRHxh1vGnAd8A3ktS8exh4JSIeL3qGA9J\nBogTocWR+5CkTm1VgLkadm1VO4gnd/WPdmurPpbuukbSNklrJR2b7nNtVZuTC0CXX6u1Vb9FUpho\nv6TbgdXptz8H3ApcXePjrq1qQJKwfO65xit3VSZ3rVuXdBpQ/67I8HAy7PHKXQcrTG3VWd8/Gbg/\nIk53bVVrxHVPeif3p1XT2qp3AL+IiPGq/SdGxPPp9jjw3oj446qk55nMJD3fWd1DuMMwJ0J7oxvz\nMN4PXAmMzbqF+nlJT0raBowA4wARsZ2k5up24AHgE+4dbDYnQsvJMz2t59p5ynXFisbRg59yfSPP\n9LTSa/UpV4C1axvPCPVTrtlzhGGF0EoidM0auPde2LHDidBWeYk+6xtOhObPQxLrG06EFp8jDCsk\nL/eXD0cY1pe83F8xOcKwwnLdk+w56Wl9zYnQbHlIYn3NidBicYRhpeFEaOccYdjAaDcReuON9c/r\nRGjz3GFYabSyRujERDI8Wbw4mRE6NtZ4TdGREXcajXhIYqVTSYRCEnU0KjUwNpZ0Gpde6hIGHpLY\nwKkkQqFxpAGwenVry/01s4zgoHKHYaXVak7jlVdc96RTHpJYqXm5v9Z44pYNPE/uap5zGDbwPLmr\nexxhWF9pdXLX3r3N3T3pl8lduUYYkt4k6XFJWyVtl3RTun+epI2SnpX0UFURI9dVtZ7KKxHqyV2J\nZsoMHBURv5Z0GPDPwF8AFwA/j4hbJF0LvDUirmumrmp6TkcYlhsnQmvrWtJT0lHAI8DHgHuBkYiY\nTos1b4qIUyVdD7weEZ9PP/NdYFVEPDbrXO4wLFftJkKHhmBqqv65y9xp5J70lHSIpK3ANDAVEU8D\n8yNiOj1kGpifbruuqhVCO4nQoSHYt6/xkGOQE6ENa6umw4klkn4TeFDS2Kzvh6R64cKc33NtVeuG\nZcuSDqCZWq5TU8k08ptvTjqDeonQZctmhj7bthU30uhJbdUDB0t/Bfw38GfAaES8IOlEksjj1Gbq\nqqb7PSSxrmk1p7FuHRx5ZHL3pJmSB2UanuR9l+S4yh0QSW8GPgBsATYAS9PDlgL3pdsbgCskHSFp\nIbAYeKLdH84sC60+5XrllTPPnjRaI3TQnnKtG2FIOp2kEPMh6evOiPhrSfNI6qe+A9gNXBYR+9LP\n3ABcBbwGrIiIB+c4ryMM6zrPCPXUcLOW5TW5qwyPxrvDMGvDoC7352dJzNrguiftcYRhA2sQ6554\nSGLWgUFLhHpIYtYBPxrfGkcYZqlWEqGV0gWXXNJ4IeIiJUIdYZhlpJVE6CWXwHHHwW23DVYi1BGG\nWZV+T4Q6wjDLUGUa+aJFSTLUOY2DOcIwq6EfZ4Q6wjDLieuevJEjDLM6+m25P0/cMstZP03u8pDE\nLGee3DXDEYZZC8pe98QRhlkXDXrdk4aLAJvZjEpUsX5944WFK1HFpk3Ja8+e+sdXIotHHjn4fZF4\nSGLWhlYToZOT8Oijvb974iGJWQ+0mghdvjzpLMqeCG23tuoqSXskbUlf51V9xrVVbWBUOo1mchqT\nkzOrkY+PNz7vxRcXb3JX3Q4jIv4HGIuIJcDvAGOSziYpTvSFiBhKXw8ApLVVLwdOA84FvizJUYz1\ntUFKhDZT+ezX6eYRwKHAS+n7ucZBFwJ3RcSrwG5JO4EzgcfmONasL7SbCL3nnuRrvZxG0RKh7dZW\nBbhG0jZJayvFjnBtVRtQ1U+5NhpGTEzMDE927SpXsaR2aquOArcDq9NDPgfcClxd6xRz7XRtVes3\nw8PJqzK5a/Pm+rVcJydnJnetX1//zkklsmh03tkKUVs1Iv6mat/JwP0Rcbprq5olilr3pCe1VSWd\nUHXYh4Cn0m3XVjWjf+uetFtb9WvAEpLhxi7gzyNiOv2Ma6uaUczl/vx4u1mBtfNo/N69cNZZyXbW\nnYZnepoVWDuPxp91FuzcWcwZoY4wzLqkCIlQRxhmJdEPiVBHGGZd1OtEqCMMsxIpe90TRxhmPdKL\nuieOMMxKqox1TxxhmPVQt+ueeOKWWcl1s+6JhyRmJVemuifuMMwKop3l/r7zncZ5iurl/jrlDsOs\nQFpNhAKsXdv8cn+dcg7DrGBaSYSuWQP33gs7djSXCHXS06wP5ZUIddLTrA/lmQjthDsMswLLq+5J\nu9xhmBVcXjNC2+EchlkJZDUj1ElPswGRRSK0K0lPSYemNVTvT9/Pk7RR0rOSHqoqZFTI2qpZ1mVw\ne26vV+0VIRHabA5jBbCdmaJE1wEbI+IU4Hvp+8LWVi3qPwC35/baaa/dGaFZJEKbKZW4APgg8BVm\n6qleQFJ+gPTrRen2gdqqEbEbqNRWNbMMtbvcX6ea+d9/AvgM8HrVvvmVOiQkNVfnp9uurWrWBdUr\ndzVSqeW6eHEGDUdEzRdwPvC36fYoSUlEgJdmHffL9OuXgD+p2v8V4OI5zht++eVXb171/uYbvRoV\nY34fcIGkDwJvAn5D0p3AtKQTIuIFSScCL6bH/ww4qerzC9J9B+kkS2tmvVN3SBIRN0TESRGxELgC\n+H5EfJSkhurS9LClwH3ptmurmvWxRhHGbJF+vRm4W9LVwG7gMoCI2C7pbpI7Kq8Bn/CEC7P+0ZOJ\nW2ZWTj2fI2Fm5eEOw8ya5g7DzJrmDsPMmuYOw8ya5g7DzJrmDsPMmvb/l3ireRMKSK8AAAAASUVO\nRK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10cedbc10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.sparse\n",
    "from scipy.sparse import csc_matrix\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "n = 20\n",
    "ex = np.ones(n);\n",
    "lp1 = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr'); \n",
    "e = sp.sparse.eye(n)\n",
    "A = sp.sparse.kron(lp1, e) + sp.sparse.kron(e, lp1)\n",
    "A = csc_matrix(A)\n",
    "plt.spy(A, aspect='equal', marker='.', markersize=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Sparsity pattern \n",
    "\n",
    "The ``spy`` command plots the sparsity pattern of the matrix: the $(i, j)$ pixel is drawn, if the corresponding matrix element is non-zero.\n",
    "\n",
    "Sparsity pattern is really important for the understanding the complexity of the sparse linear algebra algorithms.  \n",
    "\n",
    "Often, only the sparsity pattern is needed for the analysis of \"how complex\" the matrix is."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Sparse matrix: formal definition\n",
    "\n",
    "The formal definition of \"sparse matrix\" is that the number of **non-zero elements** is much less than the total number of\n",
    "\n",
    "elements, so you can do the basic linear algebra operations (like solving linear systems at the first place) can be done faster, \n",
    "\n",
    "than if working for with the full matrix.\n",
    "\n",
    "The **scipy.sparse** package has tools for solving sparse linear systems:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.contour.QuadContourSet instance at 0x1020a6f38>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD+CAYAAADBCEVaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGttJREFUeJzt3WuoXXeZx/FvW413cRIxPW0yNYHTJ4kodmSCWG8tVjpe\n0r5qKyhhKr6pQ4uCmPhiYF5Mp/pCxRn6QrwQy5hpQIxxROxNxWGwpdCCmOQxJRMwnfTESa0XREjb\nMy/22nV1Z9/O2eu//rffB0p39jknWWfttb/72f+zzt4Xra6uIiIiZbk49gaIiEj3FHcRkQIp7iIi\nBVLcRUQKpLiLiBRIcRcRKdBL5vkkMzsF/B54Djjv7rvNbCNwL3AFcAq4yd2faT5/P3Br8/m3u/t9\nk/7ujRs36lxMEZE1evrppy+a9vF5J/dV4L3ufpW7726u2wfc7+5XAg82f8bMdgE3A7uA64G7zUzP\nEEREerSW6I4+SuwBDjSXDwA3NpdvAA66+3l3PwU8AexGRER6s5bJ/QEze9TMPtFct9ndV5rLK8Dm\n5vJlwOnW154GLl94S0VEZG7zxv1qd78K+Dvgk2b2rvYH3X2VwQPAJFpXFxHp0Vxxd/czzf9/A3yX\nwTLLipldCmBmS8DZ5tOfBLa2vnxLc52IiPRk5tkyZvZK4BJ3/4OZvQp4P/BPwBFgL/D55v+Hmy85\nAnzbzL7IYDlmGXgkwLZLAXbu3Bl7E7J17Nix2JsgCbto1qtCmtk2BtM6DB4M/t3d/6U5FfIQ8Ndc\neCrk5xicCvkscIe7/2jS369TIcNTQKVveuAJb9apkDPjHpriPp6CLDKgB4rxFPeEKNgiYdX0QKC4\n90jxFklbSfFX3ANQxEXKkmP0FfcFKeQidUo9+Ir7Oijo9dmxY8fcn3v8+PGAWyIpSjH0ivucFPRu\nrSWWEp4ekLqTSugV9ykU9AspyjKJHiAuFDP0ivsYtUVdwZYYanowiBF5xX1EaWFXuCVXpcW/78Ar\n7o2co66AS21yDn9fkVfcySvsCrnIeDkFv4/AVx/3lMOukIssJuXghw581XFPLeyKuUhYqcU+ZOCr\njXsqYVfQ57d9+/bYm5C8kydPxt6EbKQS+lCBrzLuscNec9AV6HTV/MAQO/QhAl9d3GOGveSoK9r1\nKPlBIGbkuw78rLjPfJs9ma6UoCveMjTrWMg5/u37a+xpPrSiJvc+p/Zco66ISyi5Rr/PyHc5vWty\nDyCHsCvi0rdJx1zq0d+xY0eRU3wxk3sfU3vKUVfMJRcpx76PyHc1vWty70hqYVfMJVejx25KsS9p\nii9icg85tacUdQVdSpdS6ENGvovpXZP7AlIIu4I+3vLycuxN6MyJEydib0Iy2sd77NDnPsUr7hPE\nDHttQS8p1Oux3u+/9AeFFEKfc+CzX5YJsSQTK+ylRb32aKegxAeAGKEPEfhFl2a0LLNGfYc996Ar\n4GmbdvvkGv7hfabPyOc4wSvuLX2GPceoK+RlGXd75hT8viOfW+CzjnvsFwhbrxzCrpDXKcfgx5jk\nu7Bz586gLwmcddy71MfUnmrUFXKZJpfgb9++PXjgc5reFXfChz3FqCvosoj28ZNS6PuY4nMJvOIe\nWCphV8wllNFjK4XYh57icwh8tqdCdrXeHnJqjx12BV1iix36kIHvKu7rXXfXqZBThAq7oi4yMDwW\nY0U+5DJN6tN71XEPIVbYaw567AfT9cjtzI5FxV6j7+OHrampdlkmxNQeIzIlRz3HaHet5CDFiHyI\n/bno9B51WcbMLgEeBU67+4fNbCNwL3AFcAq4yd2faT53P3Ar8Bxwu7vft64tz0zfISol6gr4dNP2\nT+7hj7FkU9MEP9fkbmafBt4GvMbd95jZF4D/c/cvmNlngb9y931mtgv4NvC3wOXAA8CV7v78pL97\nPZN7alN7n4HKNeqKeH9yjVefke96H8WY3hee3M1sC/AB4J+BTzdX7wHe01w+APwE2AfcABx09/PA\nKTN7AtgN/HzNWy4XyCXsCnlc4/Z/DsHvc5KvYYKfZ1nmS8BngNe2rtvs7ivN5RVgc3P5Ml4c8tMM\nJvhi9RGy1KOumKcv5Xc/GrW8vBz9FMoSTI27mX0IOOvuj5nZe8d9jruvmtm0pZW4P7Ed0eWSTOio\npRp1xTx/qce+jym+y+k9xdMiZ03u7wD2mNkHgJcDrzWze4AVM7vU3Z8ysyXgbPP5TwJbW1+/pbmu\nOLWFXUEvWwpvjDFO6Cm+5OWZi6d90N0/5+5b3X0bcAvwkLt/DDgC7G0+bS9wuLl8BLjFzDaY2TZg\nGXgkzKaXK5Wwb9++/YX/pB6p3e6h7w9dfZ8pvC1n21p/iWm4xHIXcMjMPk5zKiSAux81s0PAUeBZ\n4DZ3T2ZZpqudH/Kgjx32VO7QkoZUJvrYv+kaWoiX/83ul5gWOQ2yi7iHip+iLrmIvYwRKvBdfF+L\nrLuvNe6zToWcuiwj/YgZ9pSefkseYh8zsQehXFQT91Sn9lgHauw7qOQv5jEU4n7TxfeS0rq7Xjgs\nohhhrz3oS0tLwf7uM2fOBPu7Uxbrbe50Pvx0ivucuo5i32EvOeohg70Wa9mOEh8IYkS+68CXdGqk\n4l6BUsKeSsS7MO17yT38JQVyPVL5haYq4r7oOliuU3uuUS8p4usx7vvPLfh9TvGa3serIu4pUdgv\nVHvM5zG6j3KJfV+h1Pr7hRT3GbqMZB9hzyHqivnicop9X1N8l4EvYXqv5lTIGqQc9qWlpRf+k+7l\nsH9TPj5TsOj7VIzKKu7r+eZTOe80l9fH6FIOwSlRyvu9lhfcS6E7WcW9b10diLWFPdWw1CjF2yKX\nwKd2v1orrblnLpUDMLWAyIu1b58U1udLWNNOnSb3wEJO7SmEPcXJUKZL5TYr+dVVU6C4T5BCOKeJ\nvX2pBELWL4XbMPZxPEvq2zeN4h5QqOkh5gGXQhCkW7Fv01JfRjv2D1WLjnvsnRtC7LBLuUoMfM2K\njvt6dXGgpfqSpOsRe7KT/sS8rVN9Se1cH3gUd5lKUa+Tbvf8Ke4BlDC1a1qXGMdAqtN7jhT3DMQI\nu8hQCYGvkeI+ovYDS2GXcWo/LnLsgn5DtWNdPwXs86Cq7Q68adOmhf+Oc+fOdbAleVhaWurtt1u7\n/g3WGl8SWHFPmMK+mC7ivei/UVr8cw58DDHflanYuJd4jnsopYS9j5iv1eg2lRD7PgNfm507d3Ls\n2LFO/i6tuXeoyyWZvqb2nMO+adOmF/2Xgxy3eZy+jpvc3iwnJdnEvesXsh8nxx+aLCLHsJcQxrac\nv58cj59F5NaHbOJekz4OopzumDkHcC1y/D77OI5yi2oqFPeO5PSUL5ew5xa6LuX0vedyPEFe99NF\nKe6JCT2l5HBHzClsoeUyzYc+rjS9r53iLsnIIWIxad/IWijuFUl1alfU55fyvkr1+KqV4t6BHN6Q\nN8U7XsqhSl2q+y7kcZbCG9bntDykuDdyutFKkGKYcqT9KJMo7hVIbWpXkLqV2hSf2vEWW6zfli8y\n7jm+9EANzxxSi1Bpati3NdxPujL1tWXM7OXAT4GXARuA77n7fjPbCNwLXAGcAm5y92ear9kP3Ao8\nB9zu7veF23yZJZUpqobwpGDTpk1JvH6NXn8mvqmTu7v/GbjG3d8KvAW4xszeCewD7nf3K4EHmz9j\nZruAm4FdwPXA3WZW5LODoZp+KWK9FPZ+aX/PVsP9dmZ43f1PzcUNwCXAb4E9wIHm+gPAjc3lG4CD\n7n7e3U8BTwC7u9zgEoV6qpnC1K7QxJHCfg91/GlpZj4z425mF5vZ48AK8GN3/yWw2d1Xmk9ZATY3\nly8DTre+/DRweYfbKxlJITA10/6v2zyT+/PNsswW4N1mds3Ix1eB1Sl/xbSPSaEUljTEvh1SePZY\nq7nXw939d8APgLcBK2Z2KYCZLQFnm097Etja+rItzXXSM92pRMLIZVloatzN7PVm9rrm8iuA64DH\ngCPA3ubT9gKHm8tHgFvMbIOZbQOWgUdCbHiXcrmxchF7WpQX0+1Rp1mT+xLwULPm/jDwfXd/ELgL\nuM7MfgVc2/wZdz8KHAKOAj8EbmuWbWSC0h5YFJI0xbxdQjyLLO1+E8LU89zd/RfA34y5/mngfRO+\n5k7gzk62LnGpnk6lJRmR2ZaXlzlx4kTszQim6HPQpV+a2tOm26cuiruISIEU98LEWpLRVJiHWLeT\nlgr7p7iLiBRIcY+olJ/4a2rPSym3Vyn3n1AUdxGRAinuIiIFUtxFKlTK0oxMprgXJMYZCYqEzEtn\nzPRLcRcRKZDiLiJSIMVdRKRAiruISIEUdxGRAinuIiIFUtxFRAqkuIuIFEhxFxEpkOIuIlIgxV1E\npECKe0HOnDkTexNEJBGKuyzk3LlzsTdB1iHG7abho1+Ku4hIgRR3EZECKe4RnTx5MvYmdEJLM3nR\n7VUHxV1EslTKcBSK4l4Y/dBKUqTjsn+Ku3RCT/XzoNupHoq7dEbhEEmH4r6AEydOxN4Ekbnpwbcu\nintkIX4oFHN9UwFJU8zbJdX19tKHM8Ud/dRdJDe6z86muEvnNL2nRbdHt3J5YFHcCxX7qbCCkobY\nt0Ps47BmL5n1CWa2FfgW8AZgFfiqu3/FzDYC9wJXAKeAm9z9meZr9gO3As8Bt7v7fWE2X1J27tw5\nNm3aFHszqhU77BLXPJP7eeBT7v4m4O3AJ81sJ7APuN/drwQebP6Mme0CbgZ2AdcDd5uZniFMEepp\nXgpTkwITRwr7PdTxl8uySGwzo+vuT7n7483lPwLHgMuBPcCB5tMOADc2l28ADrr7eXc/BTwB7O54\nu5NR+k/cu5BCaGqi/S2wxjV3M3sjcBXwMLDZ3VeaD60Am5vLlwGnW192msGDgUSQwvQOCk5fUtnP\nqRx3k9QwlM0ddzN7NfAd4A53/0P7Y+6+ymA9fpJpHxPqeKqZSnhKpf0rbXPF3cxeyiDs97j74ebq\nFTO7tPn4EnC2uf5JYGvry7c01/Xm+PHjff5zyUtpilKAwkhpv4Y83moYgroyM+5mdhHwdeCou3+5\n9aEjwN7m8l7gcOv6W8xsg5ltA5aBR7rbZMnduXPnkopRzrQv0xdr2Jxncr8a+ChwjZk91vx3PXAX\ncJ2Z/Qq4tvkz7n4UOAQcBX4I3NYs2yRtkYmgq/W7kFNJStP7kKK0mBT3X4rHWZdyeuYw8zx3d/8v\nJj8IvG/C19wJ3LnAdkklhoHS+fDzSzHqfegqrDX8MBX0G6pVSXmq0vLCbKnvo5SPrxop7okJ/bQv\n9Ttg6gGLIYd9Evq4ymk5JBWKe0dyeqqXeuAhj6CFlss+yOF4qpHinqA+ppRc7pDDwOUQuS7U9v3O\no8v7Q05D2KKyifuxY8eC/xt66pe2ksOX6/eVy5DQhdz6MPNsGYnj5MmTbN++Pei/cebMGZaWloL+\nG6G0Q5jjmTY5hnxUH2HPLagpUdw7dOLECZaXl2NvxprkHPih0VCmGPsSYt6W48Sew5JMlysUxcb9\n+PHj7NixI/ZmLKSP6R3KCHzbpJD2Ef3SIj5OX2EvYWqP+VIoxcZd1qa0wI9TQ3hDy3Fir1U2P1Dt\ny6LTQtdP/fqcXs6cOaM7r4zV97HR9XG/6P0yx2cRinsG+j6wFHhp6/t4yDGkKVLcA8jhBzezKPAC\nZRwHJdwf10Nxz0SMaUbLNPWKddtrau+O4j5GFwdYiGkh1oGvwNcl1u2dathT3a5Zio57ie/IFDPw\ninzZYt7GoY7rWpdkoPC4xxbqwIo5SSjyZdJt2r3Yw6XiLuuiGJQhhQdrTe1hKO4TpP6uLymsAw7D\nEDsOsnap3G6phz2F+9l66TdUM9bXyxPMYxiK0n/LNWcpxLwt53DmQJN7D0I+PUztDpLKRCh/keJt\nEvK4rX05Zqj4yX2RFxBLaTKeZnhHSWlb2zHRNN+/1GI+lNowMk1O2zpOVpN7H2/YEUof00SqB6PW\n5vuR+n7u4/hMZWqPfaYMZBb3GHJ7i69UAz/UDlCqEcpJLvsyt7DHuB91PbwWvyyTmj7e0COX5SS4\ncPlASziTpR7wSVIfOEqluBcqxXX4eSj2f5FrzIf6jHoqyzEpuWh1dTXqBmzcuHFNG7Bz5851/TuL\nvitT15Hs++34cov8LCVFP/eIj5Nz2Bfd9vWut691Webpp5++aNrHNblH0vf7rea0VDOPeYKYwgNA\nieGepu8lGE3skynucwoRxxiBh/Km+ElqC2tMMdbVS3rl1RB0tkxkMSaPkydPFnUQSzyxjiVN7LNV\nE/cuzjtN/XUw1kqRl/WKeeyk/HpNKZzfPlRN3LtSWuBBkZf5xT5WUg57arJbcz927Ni6z5hJXd9r\n8KPaB3gt6/IyWyrhK3kpJsRv32d3KiSs/3RIWPyUyKHQ8YsZ+VEKfX1SCfpQDi++t8iSzHrirlMh\nAwl9amHsKb5NE30dUgs6hJ/WU/yeu6K4JyylwA8p9GVJOW4lL8P0YeayjJl9A/ggcNbd39xctxG4\nF7gCOAXc5O7PNB/bD9wKPAfc7u73Tfv7+16Wge6WZqC/wKUW+XEU+/SlHPOhvqLe5b5Y9CyZWMsy\n3wT+FfhW67p9wP3u/gUz+2zz531mtgu4GdgFXA48YGZXuvvza97yTPT1m5/DAz7lyI+7syj48eQQ\n8lE5hj1VM+Pu7j8zszeOXL0HeE9z+QDwEwaBvwE46O7ngVNm9gSwG/h5VxsMi58xs8gbeIzT56/2\np7hUM42C34/cY9XnEkzX+yqlc9vb1rvmvtndV5rLK8Dm5vJlvDjkpxlM8MXrO/CQ9hQ/zbQ7l8I/\nWe4BH6fvdfUU92GoNyFa+Aeq7r5qZtPWzeOeazlB19M79P/iXLlHfpxZd76S459ieEKJ9bIbXUt1\naof1x33FzC5196fMbAk421z/JLC19XlbmuuqEePVF9t3lJJCP84id9A+bpeaAr1WsX8LuzbrjfsR\nYC/w+eb/h1vXf9vMvshgOWYZeGTRjQwlxPQOcV9et8Rpvis13sFTEPuUxlC3e8pTO8wRdzM7yOCH\np683s18D/wjcBRwys4/TnAoJ4O5HzewQcBR4FrjN3ZNclgkt9uun1zTNS3piB32o5gf0LF9+YKir\n15gJMb23pbROrNBLKKkEHcJHvaupfZEfpurlB+YQanlmKPYU3zZ6B1TsZb1SinlbLmEPTXHvSarv\ngqTYy7xSjflQzUsw4yjujdDT+1CqkR8adwdW8OuUesyH+ox6LlM7ZL7mDt2tuw/1Efi2VCM/i4Jf\njlwiPk7OYV/0l5e05r5GfU3wQ6lP8pNMCoKin66cI94WY/klp4l9SHEfo+/AQ76RHzVPQPQA0L1S\nwj1NrDX1HMMOBSzLQPdLM0N9B35U7qHvQs0PBDUEex4xf1AaKuxdvJ6MlmUWEGOCbytlml9EiMCF\neMBQiLuVwpkvuU7sQ0VM7hBueof4E/yommMvZUoh5m0hw97Vq0Bqcu/A8IZOJfJ6qzspQWpBH8p9\nYh8qZnKHsNP7UCqBn0bBl9SkGvK2vqKuyT1RqU3x4+jdjyS2HGI+1OekHuqNOcYpanKHfqb3tpQj\nP4uCL4vKKeKj+l5+6TrsmtwDy2GSn6Tmdz2S+eQc70lKWVOfpbjJHfqf3ttyjPwi9ACQrxLDPU3M\nqIdYjpk1uRcZd4gb+KHaQj8PPRiEU1us55HClB5qnb3auEMagR9S6MPJ8QFDIQ4nhaC3xYq71tx7\nMnrAKfbdUSgltaAP9Xl2zKiiJ/ehlCb4SRR7kfmkGvJRocNe9bJMWw6BH6XgS+1yCXlbX9O64t6S\nY+DHUfSlNDlGfJw+l2EU9zFKifwkir+kppR4TxJjbV1xn6L0yM+iBwFZVOnRniXmD0wV9znUHvn1\n0ANDOWoP9HrEjPqQ4r5GCn069ACi8KYmhagPKe4LUOhF6pZSzEcp7h1S7EXKlnLMRynugSn4InnK\nKeTjKO4RKfwiceUe8GkU9wzoQUBkbUqO9rwU98LpgUFypUAvRnGXzuiBpHwKbj4UdxGRAs2K+8V9\nbYiIiPRHcRcRKVCQd2Iys+uBLwOXAF9z98+H+HdERGS8zid3M7sE+DfgemAX8BEz00/iRER6FGJZ\nZjfwhLufcvfzwH8ANwT4d0REZIIQcb8c+HXrz6eb60REpCch1tzXdGrjrNN5RERk7UJM7k8CW1t/\n3spgehcRkZ6EmNwfBZbN7I3A/wI3Ax8J8O+IiMgEnU/u7v4s8A/Aj4CjwL3urt9pFhHpUfSXHxAR\nke7pN1RFRAqkuIuIFEhxFxEpUJDXlplXLa9BY2bfAD4InHX3NzfXbQTuBa4ATgE3ufszzcf2A7cC\nzwG3u/t9Mba7S2a2FfgW8AYGvwvxVXf/SoX74eXAT4GXARuA77n7/tr2A7zwUiWPAqfd/cOV7oNT\nwO8ZfF/n3X13V/sh2uRe2WvQfJPB99m2D7jf3a8EHmz+jJntYnD66K7ma+42sxKeYZ0HPuXubwLe\nDnyyub2r2g/u/mfgGnd/K/AW4BozeyeV7YfGHQzOqBue1VHjPlgF3uvuV7n77ua6TvZDzB1UzWvQ\nuPvPgN+OXL0HONBcPgDc2Fy+ATjo7ufd/RTwBIN9lTV3f8rdH28u/xE4xuBlKaraDwDu/qfm4gYG\nz1p/S2X7wcy2AB8AvgYMf0u9qn3QMvpb+p3sh5hxr/01aDa7+0pzeQXY3Fy+jBf/Rm9x+6X5Bber\ngIepcD+Y2cVm9jiD7/fH7v5L6tsPXwI+Azzfuq62fQCDyf0BM3vUzD7RXNfJfogZd51g33D3Vabv\nj2L2lZm9GvgOcIe7/6H9sVr2g7s/3yzLbAHebWbXjHy86P1gZh9i8POnx7hwagXK3wctV7v7VcDf\nMViqfFf7g4vsh5hxr/01aFbM7FIAM1sCzjbXj+6XLc112TOzlzII+z3ufri5urr9MOTuvwN+ALyN\nuvbDO4A9ZvY/wEHgWjO7h7r2AQDufqb5/2+A7zJYZulkP8SM+wuvQWNmGxj8oOBIxO3p2xFgb3N5\nL3C4df0tZrbBzLYBy8AjEbavU2Z2EfB14Ki7f7n1odr2w+vN7HXN5VcA1wGPUdF+cPfPuftWd98G\n3AI85O4fo6J9AGBmrzSz1zSXXwW8H/gFHe2HaHGv6TVozOwg8N+Di/ZrM/t74C7gOjP7FXBt82fc\n/ShwiME++SFwW/PULHdXAx9lcHbIY81/11PfflgCHmrW3B8Gvu/uD1Lffmgbfj+17YPNwM9ax8J/\nNqc2drIf9NoyIiIFKuVcURERaVHcRUQKpLiLiBRIcRcRKZDiLiJSIMVdRKRAiruISIH+H+O6NUen\nDwcmAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10ed31110>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.sparse\n",
    "import scipy.sparse.linalg\n",
    "import seaborn as sns\n",
    "from scipy.sparse import csc_matrix, csr_matrix\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "n = 512\n",
    "ex = np.ones(n);\n",
    "lp1 = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr'); \n",
    "e = sp.sparse.eye(n)\n",
    "A = sp.sparse.kron(lp1, e) + sp.sparse.kron(e, lp1)\n",
    "A = csr_matrix(A)\n",
    "rhs = np.ones(n * n)\n",
    "sol = sp.sparse.linalg.spsolve(A, rhs)\n",
    "plt.contourf(sol.reshape((n, n)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## What we need to find out to see how it actually works\n",
    "\n",
    "**Question 1:** How to store the sparse matrix in memory?\n",
    "\n",
    "\n",
    "**Question 2:** How to solve linear systems with sparse matrices fast?\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Compressed sparse row (CSR)\n",
    "A matrix is stored as 3 different arrays: \n",
    "```python\n",
    "ia, ja, sa\n",
    "```\n",
    "where:\n",
    "\n",
    "- **ia** is an integer array of length $n+1$ \n",
    "- **ja** is an integer array of length **nnz** \n",
    "-  **sa** is an real-value array of length **nnz**\n",
    "- **nnz** is the total number of non-zeros for the matrix\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Idea behind CSR.\n",
    "- For each row $i$ we store the column number of the non-zeros (and their) values\n",
    "- We stack this all together into **ja** and **sa** arrays\n",
    "- We save the location of the **first non-zero element** in each row"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## CSR helps for matrix-by-vector product as well\n",
    "```python\n",
    "   for i in range(n):\n",
    "       for k in range(ia[i]:ia[i+1]):\n",
    "           y[i] += sa[k] * x[ja[k]]\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Let us do a short timing test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10000 loops, best of 3: 46.7 µs per loop\n",
      "10000 loops, best of 3: 56.8 µs per loop\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.sparse\n",
    "import scipy.sparse.linalg\n",
    "from scipy.sparse import csc_matrix, csr_matrix, coo_matrix\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "n = 60\n",
    "ex = np.ones(n);\n",
    "lp1 = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr'); \n",
    "e = sp.sparse.eye(n)\n",
    "A = sp.sparse.kron(lp1, e) + sp.sparse.kron(e, lp1)\n",
    "A = csr_matrix(A)\n",
    "rhs = np.ones(n * n)\n",
    "B = coo_matrix(A)\n",
    "%timeit A.dot(rhs)\n",
    "%timeit B.dot(rhs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "As you see, **CSR** is faster, and for more **unstructured patterns** the gain will be larger."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Sparse matrices and efficiency\n",
    "Sparse matrix give complexity reduction.  \n",
    "\n",
    "But they are **not very good** for parallel implementation.  \n",
    "\n",
    "And they do not give maximal efficiency due to **random data access**.  \n",
    "\n",
    "Typically, peak efficiency of $10\\%-15\\%$ is considered good.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## How we measure efficiency of linear algebra operations\n",
    "\n",
    "The standard way to measure the efficiency of linear algebra operations on a particular computing architecture is to \n",
    "\n",
    "use **flops** (number of floating point operations per second)\n",
    "\n",
    "The peak performance is determined as \n",
    "\n",
    " **frequency x number of cores x pipeline size x 2**\n",
    "\n",
    "We can measure peak efficiency of an ordinary matrix-by-vector product."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time:  8.5e-04, Efficiency:  2.3e+00 Gflops\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import time\n",
    "n = 1000\n",
    "a = np.random.randn(n, n)\n",
    "v = np.random.randn(n)\n",
    "t = time.time()\n",
    "np.dot(a, v)\n",
    "t = time.time() - t\n",
    "print('Time: {0: 3.1e}, Efficiency: {1: 3.1e} Gflops'.\\\n",
    "      format(t,  ((2 * n ** 2)/t) / 10 ** 9))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time:  1.4e-04, Efficiency:  2.2e-02 Gflops\n"
     ]
    }
   ],
   "source": [
    "n = 1000\n",
    "ex = np.ones(n);\n",
    "a = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr'); \n",
    "rhs = np.random.randn(n)\n",
    "t = time.time()\n",
    "a.dot(rhs)\n",
    "t = time.time() - t\n",
    "print('Time: {0: 3.1e}, Efficiency: {1: 3.1e} Gflops'.\\\n",
    "      format(t,  (3 * n) / t / 10 ** 9))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## How to make things more efficient\n",
    "\n",
    "Sparse matrix computations dominate the linear algebra computations nowdays.\n",
    "\n",
    "They allow us to work with much larger matrices, but they utilize only $10\\%-15\\%$ percent of the peak computer performance.\n",
    "\n",
    "It means, that our computer architecture **is not well suited** for standard sparse matrix algorithms.\n",
    "\n",
    "There are many possible solutions of the problem, for example:\n",
    "\n",
    "1. Use block sparse format\n",
    "2. Reorder equations to make them more \"block\"\n",
    "3. Instead of vectors, use \"block vectors\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time:  2.9e-04, Efficiency:  1.0e-01 Gflops\n"
     ]
    }
   ],
   "source": [
    "\n",
    "n = 1000\n",
    "k = 10\n",
    "ex = np.ones(n);\n",
    "a = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr'); \n",
    "rhs = np.random.randn(*(n, k))\n",
    "t = time.time()\n",
    "a.dot(rhs)\n",
    "t = time.time() - t\n",
    "print('Time: {0: 3.1e}, Efficiency: {1: 3.1e} Gflops'.\\\n",
    "      format(t,  (3 * n * k) / t / 10 ** 9))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Florida sparse matrix collection\n",
    "\n",
    "There are many other types of matrices besides tridiagonal!\n",
    "\n",
    "[Florida sparse matrix collection](http://www.cise.ufl.edu/research/sparse/matrices/) which contains all sorts of matrices for different applications.\n",
    "\n",
    "It also allows for finding test matrices as well! \n",
    "\n",
    "We can have a look."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<iframe src=http://yifanhu.net/GALLERY/GRAPHS/search.html width=700 height=450></iframe>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import HTML\n",
    "HTML('<iframe src=http://yifanhu.net/GALLERY/GRAPHS/search.html width=700 height=450></iframe>')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Visualization of sparse matrices and graph\n",
    "\n",
    "Sparse matrices and fast algorithms (especially for linear systems) have deep connection with graph theory.\n",
    "\n",
    "First of all, sparse matrix can be treated as an **adjacency matrix** of a certain graph:\n",
    "\n",
    "The vertices $(i, j)$ are connected, if the corresponding matrix element is non-zero.\n",
    "\n",
    "Can you guess why it is important."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Graph structure is important for LU decomposition\n",
    "\n",
    "Why sparse linear systems can be solved fast? \n",
    "\n",
    "Because the LU-decomposition is often also sparse (i.e., matrices $L$ and $U$ are sparse as well).\n",
    "\n",
    "And how sparse they are, and which variables/equations select for the elimination at each step, is governed by the sparse structure.\n",
    "\n",
    "And solving linear systems with **sparse** triangular matrices is very easy.  \n",
    "\n",
    "\n",
    "\n",
    "<font color='red'> Note that the inverse matrix is not sparse! </font>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-2.  1.  0.  0.  0.  0.  0.]\n",
      " [ 1. -2.  1.  0.  0.  0.  0.]\n",
      " [ 0.  1. -2.  1.  0.  0.  0.]\n",
      " [ 0.  0.  1. -2.  1.  0.  0.]\n",
      " [ 0.  0.  0.  1. -2.  1.  0.]\n",
      " [ 0.  0.  0.  0.  1. -2.  1.]\n",
      " [ 0.  0.  0.  0.  0.  1. -2.]]\n",
      "[[-0.875 -0.75  -0.625 -0.5   -0.375 -0.25  -0.125]\n",
      " [-0.75  -1.5   -1.25  -1.    -0.75  -0.5   -0.25 ]\n",
      " [-0.625 -1.25  -1.875 -1.5   -1.125 -0.75  -0.375]\n",
      " [-0.5   -1.    -1.5   -2.    -1.5   -1.    -0.5  ]\n",
      " [-0.375 -0.75  -1.125 -1.5   -1.875 -1.25  -0.625]\n",
      " [-0.25  -0.5   -0.75  -1.    -1.25  -1.5   -0.75 ]\n",
      " [-0.125 -0.25  -0.375 -0.5   -0.625 -0.75  -0.875]]\n"
     ]
    }
   ],
   "source": [
    "#Indeed, it is not sparse\n",
    "n = 7\n",
    "ex = np.ones(n);\n",
    "a = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr'); \n",
    "a = a.todense()\n",
    "b = np.array(np.linalg.inv(a))\n",
    "print a\n",
    "print b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## And the factors...\n",
    "\n",
    "$L$ and $U$ are typically much better. In the tridiagonal case they are even bidiagonal!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.          0.          0.          0.          0.          0.          0.        ]\n",
      " [-0.5         1.          0.          0.          0.          0.          0.        ]\n",
      " [-0.         -0.66666667  1.          0.          0.          0.          0.        ]\n",
      " [-0.         -0.         -0.75        1.          0.          0.          0.        ]\n",
      " [-0.         -0.         -0.         -0.8         1.          0.          0.        ]\n",
      " [-0.         -0.         -0.         -0.         -0.83333333  1.          0.        ]\n",
      " [-0.         -0.         -0.         -0.         -0.         -0.85714286\n",
      "   1.        ]]\n"
     ]
    }
   ],
   "source": [
    "p, l, u = scipy.linalg.lu(a)\n",
    "print l"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 2D-case\n",
    "However, for a 2D case everything is much worser:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x1159538d0>"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQ8AAAEDCAYAAAAr7YFFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHjdJREFUeJztnX+IJOd55z+9ThTHiTQ9g3dHiqSVxeF+Lg4OEgkiJL6T\nZKygOIlsQ1AsuGS52MZwOcUkISDpj7OxOcfHYcWEw4IgJchK0MkknCxfEI4s24nDYeeUkxxhmWcw\naL1WYvUIdrrXwhjkTN0fXTWqqa3fXdVV1f39wLLV3dX1VvV0fft5n/d53+8oCAKEEKIqJ7o+ASHE\nMJF4CCFqIfEQQtRC4iGEqIXEQwhRC4mHEKIWP9Rl42Z2K/AJ4DXA/e7+31po4yxwAfhX4BV3v8HM\ndoBHgGuAs8Dt7j6refw/BX4Z2Hf3N4fPZR7fzO4Gfis8n99x979Zsq0PAe8FXgp3u8fdH1+2rfD9\nVwOfAk4BAfAn7v7HLV5fVnutXKOZvRb4W+BHgEuAz7j73S1eX1Z7rVxfrN3XAE8BL7j7rzZ1fZ1F\nHuEF/Q/gVuBNwB1m9pMtNBUAN7n79e5+Q/jcXcAT7j4Bngwf1+XPWFxDnNTjm9mbgF9ncb23Ap80\nsyp/g7S2AuDe8Pquj33plm0L4BXgd939p4CfA347/Bu1dX1Z7bVyje7+feBmd78O+GngZjN7S1vX\nl9Nem39DgA8Az4Xt0NT1ddltuQH4prufdfdXgP8JvKOltkaJx7cBD4bbDwLvrHtgd/8ycFDy+O8A\nHnb3V9z9LPBNFp/DMm3Bxde3dFthey+6+zPh9svAN4Arae/6stqD9q7xe+HmJSwi4ANaur6c9qCl\n6zOzq4C3A/fH2mjk+roUjyuBb8cev8CrX5QmCYDPm9lTZva+8Lldd5+G21Ngt+E2s47/EyyuM6Kp\na77TzL5mZg+Y2biNtszsDcD1wFdZwfXF2vtK+FQr12hmJ8zsGRbX8UV3/zotXl9Ge9De3/CPgD8A\nDmPPNXJ9XYrHqurif8Hdrwd+iUUY/O/iL7p70Oa5lDj+sm3fB1wLXAd8B/h4022Z2Y8DfwV8wN2/\nG3+tjesL2/vLsL2XafEa3f0w7EZcBfx7M7s58Xqj15fS3k20dH1m9iss8mNPkx7ZLHV9XYrHPwNX\nxx5fzXHVawR3/074/0vA/2IRhk3N7HIAM7sC2G+42azjJ6/5qvC52rj7vrsH4Zfgfl4NMxtpy8x+\nmIVwPOTuj4ZPt3Z9sfb+PGqv7WsM25gDfw38DCv4+8Xa+9kWr+/ngdvM7HngYeCtZvYQDV1fl+Lx\nFPBGM3uDmV3CIlHzWJMNmNnrzOzScPvHgF8Eng3bORPudgZ4NP0Itck6/mPAu83sEjO7Fngj8A/L\nNBT+8SPexeL6GmnLzEbAA8Bz7v6J2EutXF9We21do5m9PuoimNmPArcAT7d4fantRTdy09fn7ve4\n+9Xufi3wbuAL7v4bTV3fqMtZtWb2S7w6VPuAu/9hw8e/lkW0AYth6b9w9z8Mh6o+DZxm+aHah4Eb\ngdez6D/+F+AzWcc3s3tYDIX9gEVY/rkl2vogcBOLcDcAngfeH/Vnl2krfP9bgL8D/olXw9e7WXyh\n2ri+tPbuAe5o4xrN7M0sEoYnwn8Puft/z/t+tNTep9q4vkTbNwK/7+63NXV9nYqHEGK4qMJUCFEL\niYcQohYSDyFELVqZ22IrmLMihOiWxhOm4ZwVB97GYoz4/wJ3uPs3Gm1ICNEpbXRbVjlnRQjREW2I\nx6rmrAghOqSNnEelftBoNDraf3//QvNnk2B7+3UcHHyveEe1p/bWvL2TJy9Nne9SljbEo/aclVOn\nLmMVRWsnT17aehtqT+0Npb26tCEeR3NWgH9hMWfljrJvHo1GrUYgJ09eyksvfbd4R7Wn9ta8vWVF\nqvGch7v/APjPwOdYrF70SNWRlt3draZPSwjRMK3UeYTLqD1eZt/t7W0ODo4vjhUEAbu7W0yn8zZO\nTwjRAJ1XmJ4/f57xeHzR80EQcOrUZR2ckRCiDJ2LRxESECH6SS/EY2/vXGr0ESEBEaJ/9EI8YCEg\n+/sXMkVESVQh+kVvxCMiEpEkURJVCNEPeiceeSiJKkR/6K145BWKSUCE6J7eigcsBGQ0Si+/l4AI\n0S29Fg+A6XSuJKoQPaT34gFKogrRRwYhHnkoiSpENwxKPJREFaI/DEo8AFWiCtETBicee3vnMkdg\nQElUIVbF4MQD8kdglAMRYjUMUjwgewQmQgIiRLsMVjwilAMRohsGLx5FORAJiBDtMHjxgPwcCCiJ\nKkQbLLWGqZmdBS4A/wq84u43mNkO8AhwDXAWuN3dZ8u0M5mcZj5frGe6tbXF3t65i/aJnkuLNKIk\n6ip8YYTYFJaNPALgJne/3t1vCJ+7C3jC3SfAk+Hj2pw6dRmz2YwgCAiCgNlsxmRy+qL9JpPThRGG\nujBCNEcTq6cnEw63ATeG2w8CX6KmgGSJwXw+P4pGqppEte0LI8Sm0ETk8Xkze8rM3hc+t+vu03B7\nCuxWPehkcjrXPS6KQOq6yykCEWJ5RsvYO5rZFe7+HTM7CTwB3Ak85u7bsX3Ou/tOzmGOTmBnZ+ci\nD5fUkx6NlralHI1GHB4eLnUMIQbOUl61S4lHHDP7IPAy8D4WeZAXzewK4Ivu/m9z3hq89NJ3mUxO\nM5vl51VHo9FRwrRutyVJ212YdbYrVHvDbm9Zo+va3RYze52ZXRpu/xjwi8CzwGPAmXC3M8Cjy5xg\nlKPY37/AdDo/GlXZ2ztXOERbBnVhhKjHMjmPXeDLZvYM8FXgf7v73wAfA24xsz3greHj/APtbhVG\nHWmUiVbKIAERojq1R1vc/XngupTnzwNvK3ucEydOLN31aAJ54wpRjc4rTIuEY2sru3ajyGmu6nmo\nElWI8jRR57EUyZGTKCma1h2JisPiFaZp1abRvlW7NKpEFaI8nUceyeHSaDRlPB4zHo+PxCESg6wK\n0yRROXsdlAMRopjOxQM48qiNi8Xe3rnMqKIMed2dMkhAhMinF+IBxWKRFo0UHS8rH5I3hT+OciBC\nZNN5zqMKVSORtP3zinCSeZIoiapRGCEupjeRRx9IExutiSpEOhKPkkhAhDiOxCOBljQUohwSjwTT\n6Vy+MEKUQOKRwnQ6P5qMlxyxUSWqEAskHgUoiSpEOhKPEmTVi0hAxCYj8VgSCYjYVCQeJSiavasc\niNhEJB4libxx02bcKgciNhGJRw2ypuxLQMQmIfGoiQREbDoSjyWQgIhNpnBWrZn9KfDLwL67vzl8\nLtOP1szuBn6LhX/t74SLIq8t4/E4dcUyzcYV606ZyOPPgFsTz6X60ZrZm4BfB94UvueTZrbW0U3W\nMgFKoop1p/DGdvcvA0kbt9tY+NAS/v/OcPsdwMPu/oq7nwW+CdzAmpM3jFt24SEhhkbdqCDLj/Yn\ngBdi+70AXFmzjcGwt3dOs3HFxrH0SmLuHphZnn9CoSnLyZOXLnsalWijvcPDw1wPmt3drZV5467D\n56n2+k9d8Zia2eUxP9r98Pl/Bq6O7XdV+Fwu6+IFOp3OMy0fgiDgxIkTrSdR19lbVe0139Yy1O22\nZPnRPga828wuMbNrgTcC/7DUGQ6MqBI1DSVRxTpRKB5m9jDwfxab9m0z+49k+NG6+3PAp4HngMeB\n/+Tu3XtJdkCecZQERKwDhd0Wd78j46VUP1p3/yjw0WVOal3IqgEB5EwnBs9a12B0zd7eOba3tzNf\n12xcMWQkHi1z/vz51OUMQUsaimEj8VgRWWuCKIkqhsqgHOOGymRymvl8nlkDAsqBiOGhyKNldnZ2\nmM1mucIRoQhEDAmJR0tMJqc5deoyDg6S04LyUQ5EDAV1Wxpmd3erVJSRhcy1xVBQ5FGTKLI4deoy\nJpPTwPLCEaEkqhgCEo8aJOevzGazIwFpEgmI6DMSj4aYz+eFPrd1kICIviLxqEHa+h1bW4tEZ9zn\ndn//AkEQsL9/4aL9R6NR6vNpKIkq+ojEowKTyWl2d7c4deqyyrmNSFyyHuehHIjoIxKPkkR5jroJ\n0ajCdDQaMR6Pj9Y+rSIiEhDRJyQeKyCKVtLEp2gJwyQSENEXJB4tk7Y0YXx0ZjI5XTmakYCIPiDx\nKEHW0oIR8W5IRFFuZD6fH/u/Kkqiiq6ReBRQJBxpNFUsloeSqKJrJB4ViZKeceoUiUWJ0ioJ0zQk\nIKIrJB4Jou5GJAbJbkWZbkZRsVi8m5O1zkcVJCCiC+p61X4IeC/wUrjbPe7+ePjaYL1q492N2WyW\n2f2YTufHXhuNRhflPKKJbWWW0k/Ll1Tt9mgynVg1db1qA+Bed78+/BcJx6C9apM3bPLxaDRKvUGD\nIGhsbkud0ZfoHJREFaukrlctQFpcPmiv2ryuRlI4kjd43VGTJMscJwgCeeOKlbFMVHCnmX3NzB4w\ns6jTPmiv2rzkZVoU0hTx6f1NjNIoByJWQd3FgO4DPhxufwT4OPCejH0H41V7cHCQ6zcbf9/h4eEx\nARmPx5nHzbu+aJnCpmlKiMrQl7+f2lsttcTD3SNvWszsfuCz4cPBe9VOp/PMX+7t7e1jyc39/QtH\nuQ73b6Uet6i9NoQjQt64aq+orWWo1W0Jza0j3gU8G26vhVdtVpckrZ5jb+/cRaMlZambHC2Lkqii\nTcoM1T4M3Ai83sy+DXwQuMnMrmPRJXkeeD8svGrNLPKq/QED9apNDsUOmagSVbYOomlGPbhBgiGE\nhVHEUTXKKGqvjKdLVES2bBenDQFZ57B+3ds7efLSpbL+Wj29JHW7Jqs+7smTl7K9vZ0qNIpARJMM\npoBLlCev5F05ENEUEo81JavYTElU0RQSjzWkKNmr6fyiCSQeG4wERCyDxGMNqeIfIwERddFoy5qS\nrCzN68poOr+ogyKPDSFuRpUciVEORNRB4rGBZBWbSUBEFSQe4hgSEFEWiccGEq8yTUusSkBEGSQe\nG0qU/8hKlKqQTBQh8RCp0YeSqKIIiYfIHaaVgIgsJB4CyJ+uLwERaUg8xBFpNSARyoGIJBIPcYy9\nvXOpUYhm44okEg9RGiVRRRyJh0hFORBRRJkFkK8GPgWcYrHg8Z+4+x+b2Q7wCHANcBa43d1n4XsG\n61crXmU8HueWsmtJw82mTOTxCvC77v5TwM8Bv21mPwncBTzh7hPgyfDx4P1qxavs7Z3LndqvHMhm\nU8ar9kV3fybcfhn4BgsLyduAB8PdHgTeGW4P2q9WHGc6nWeOwARBwIkT+l3YVCr95c3sDcD1wFeB\nXXefhi9Ngd1we9B+teJiskZgQEnUTab0YkBm9uPAXwEfcPfvmtnRa+4emFmeAUyuOcy6e4GuS3tB\nEGR2Y+SNO9z26lJKPMzsh1kIx0Pu/mj49NTMLnf3F0P7yci/trJf7bqa6gypveRKY6PRiK2trWO+\nMkU5jtFo1HoSdSif5xDaa92r1sxGwAPAc+7+idhLjwFnwu0zwKOx5wfvV7tJpC1RGARBqjdvmWOJ\nzaBM5PELwH8A/snMng6fuxv4GPBpM3sP4VAtrI9frVgwn89LWWJGRJWoWhN1/SkUD3f/e7IjlLdl\nvOejwEeXOC/RIpPJaWazGaPRqNRNXtUjV+bam4HG2TaMSDjg1Zs8iiiSydAo71EXjcKsNxIPcURa\n3gMWlaZlfWCSSEDWF/m2bBh7e+eORR9l35NFkbVltI9yIOuHIo8NZG/vXGrVaNXoYjI5XTqJqghk\n/ZB4bCiRgESCESVPs0rRm0ACsl5IPDaYvb1zR05yUbciLirj8Ti3y7K3d47t7e1KbUpA1gflPMRF\n5AlGkvPnz+dWRKblVzSMux4o8hCtkiVEqkQdPhIP0QlKog4fiYdonbxRHAnIcJF4iNaZTucSkDVE\n4iFWQjSqk+UNIwEZHhIPsXKyitSURB0WEg/RCWmjMEqiDguJh+iMvCUNRf+ReIjOyJssJwHpPxIP\n0Sl5c2mUA+k3Eg/RKVnJU5C5dt+ReIjOKRIQdWH6yTJetR8C3gu8FO56j7s/Hr5HXrWiEnt753IX\nFtJkuv5RZlZt5FX7TGj89I9m9gQLIbnX3e+N75zwqr0S+LyZTdz9sOFzF2vGdDrPXeVslcZSophl\nvGoB0sba5FUrapNnbQnIG7dH1PWq/Ur41J1m9jUze8DMok6rvGrF0uR54yqJ2g+qetX+JQuv2pfN\n7D7gw+HLHwE+Drwn4+3yqlV7lcnyxk1aRrTNunyeTVPVq/bPI69ad9+PvX4/8Nnwobxq1V5j7O9f\nyBxtkTfu8m0tQ22v2tDcOuJdwLPhtrxqRaNkzcQFVaJ2SV2v2nuAO8zsOhZdkueB94O8asXqkS9M\nNyzjVft4znvkVStqU9WUSuba3aBxL9ErqgpHhCpRV4/EQ/SGusIRRwKyOiQeohc0IRwREpDVIPEQ\nnbC7u8WpU5e1VvClQrL2kXiIlROfABclO/Nm1iYZjUaFptzKgbSP7CZF50RCkucul6wmnU7nubNw\nIzQbtz0UeYiVU3VIdWtrK/dxEYpA2kHiIToh3kUp6q5EXZrRaMT29vZRhFJFRCQgzSPxEJ0QCcJ4\nPM7sriT3n07nnD9/PvUYRTkQUBK1aZTzEJ1RRjSWOUZy+DdKoioH0gyKPMTakiUs6sI0g8RDrDWa\njdseEg+xsSgHshwSD7HWxEdqkmhJw+VQwlSsPfHcR7K7oiRqfRR5iI1COZDmkHgIESIBqYbEQ2wU\n83l+abxyIOVRzkOsPZPJaebzeSmrBi1pWJ5c8TCz1wJ/C/wIcAnwGXe/28x2gEeAa4CzwO3uPgvf\nI59a0RvqdEWURC1HbrfF3b8P3Ozu1wE/DdxsZm8B7gKecPcJ8GT4OOlTeyvwSTNT10h0QlEXpGhd\nEOVA8injVfu9cPMS4DXAAXAb8GD4/IPAO8Nt+dSK3lCmmzKdznMFpMyEu02ljOnTCTN7BpgCX3T3\nrwO77j4Nd5kCu+G2fGpFbyia6r+1tcVkcrpQZJRETaeMb8shcJ2ZbQGfM7ObE68HZpb36RfK/7p7\ngaq9+uzs7HBwcHD0eHt7+9i0/J2dHWazWa4AjEaj1Nfn83mppQ+jJOrh4WHFs6/HWnnVArj73Mz+\nGvgZYGpml7v7i6HtZORbW9mnFuRVq/bSSVtR/eDg4GhBoGVXXA+CAPdvlTpOZLq9Tt64rXrVmtnr\nzWwcbv8ocAvwNAs/2jPhbmeAR8Nt+dSKUiyzenpUq1FUsxERBEFqhBHlM6qsK6Ik6qsU5TyuAL4Q\n5jy+CnzW3Z8EPgbcYmZ7wFvDx7j7c0DkU/s48qkVKaStnp5G1orq0fKDVZYhnM1mFyU/4++vkhiV\ngCwYlclIt0ywLmG22itH2s0X7w7krYo+Go2OFXCVWUE9IhKi+XzO1tbWRRFHlWMlz6MpVtxtWWoo\nSTUYonUmk9Ps7i5GNtKI/+oX3cBBEBwdp8xISZJoLdS0rsp0Wi6BGp3HpkcgEg/RKlEyMggCZrPZ\n2g17brKASDxEqySTmmmRQvy5oqKtuPXC3t650rmKKqu0l40+IjZVQDQxTrTK1tZW5eHUvFxCPCdQ\npdtSZUQl2reKKGziXBhFHqJVykQHVX/pI8oO1dalamn6unXJipB4iNbJiw7KdieSVIk66s5PqWpr\nuWlJVHVbRKsU/Ro3YfxURDRCE7W1bGVqEZvShZF4iMpUqYfIY5kZq1F5elSzUbTYz2w2OxribVM4\nIjZBQCQeohLLCMd4PD66yZsosioTtbQdZeSx7iuSSTw2kOiGim7gqGvR5hd9NBqtpIuSJIpQom1Y\nnaCs+5KGSphuGPEbJ0rwBUFQ2gCpqA6jj+ztnTsmXHt751bWpVjnJKoijw0jb3izbHek7tySLohu\n3LRuUtY6H22dx7rlQBR5bBh5w49FEcVkcppTpy479q/szVd12LMJ4r/4aZHVqrsT6xaBSDw2jGT5\ndfzXMO9m6jLx2Cb7+xdW2g1bp0IydVs2gHjXYjweHxOB+K9hm6F1F8Kzv38ht9sSkXy+za7YOiVR\nJR5rTvJGqHsTL7vsX1dJ1jpiGL+xq6yvUVZ01sUXRuKx5lT5BS26wdOGWpM3V5bAxKs8k0OnQ6SJ\nbtzQBUQ5D7EyovU8ZrPZsYrPodFk/mfISVSJxwYTd0xralm9ovUw4pFQ27Ni26Lp8x5qErWuV+2H\ngPcCL4W73uPuj4fvkVdtTyj6ZU9bx7MJ4sfsex1IHxhqEjVXPNz9+2Z2s7t/z8x+CPj70Ks2AO51\n93vj+ye8aq8EPm9mk9A4SqyYol/I2WyWORpRNS8RiUSVwqsuaj+aoM4CR0XEq32HQhnHuDSvWoC0\n7NqRVy1w1swir9qvNHCunRHdGHXXnmiasn3uKiMc8V+/+PHjU9mzSFopZFG0evlQaMJwKotVGEs1\nRaF4hC73/w/4N8B97v51M/s14E4z+03gKeD33X3Gwqs2LhSD96qN3xhRkq/LL32VL20UCZT9pWz7\nV2+oYpFG1rU0kQAdyihMHa/am4D7gA+Hu3wE+DjwnoxDrJVX7Wg0qvz+Jq+var3Esov5xn1iszg8\nPCw8r+3t7cY+hz5/X5qaL7NKb9y61PGq/Vl3/1L0vJndD3w2fLh2XrXT6fxYXsD9W5Xe37SJj/u3\nSicho/Mt8kvJO9aJEycKFyQuI2iz2ayRz6HvplZNRW+r8MbtxKvWzC6P7fYu4Nlwe+28auPDaHHD\noa6osnZnvDBrPB4zGo0Yj8dH/6bTOdPpPHd+R5M3g6hOn+tAiiKPK4AHw7zHCeAhd3/SzD5lZtex\n6JI8D7wfFl61ZhZ51f6ANfCqTX7p5/P5yhOoywx3ls2PRIsCJdsp0+0pE6rXXSF9aLQxzb+vORB5\n1RZQdOMWCciyYXYTdRLJyXDx57POvexQbXR98aFaOC66TYps37st0F5tS9MCsqxXrea2JEhGFfGR\nivganBFtV0ku+yWsOyGt6s0+tAKnNkn7LJoY2u1bIZnK02Mkh2WjeRjw6q9n8o/XdqHTsrNRgyBg\nPh/e0oHiYvq2pKHEI0byVz4rwtjfv3CUdGw759HEmqHRGqVJhjq3ZIg0+Vn3RUDUbYmRl+xKU/2u\nkoBlhljFetOHJGovIo/J5DS7u1udDIPu7m5VXo8zYhXTytOSb1mRRFWGOrdkiLTxWXc9G7dz8djZ\n2WE2mxEEwcrXeMjLiq9TjiAqNopqPaJ6j3UqF+878Vqbpihrl9EWnXdbkhnovL5h3RWo4iJRdq5H\nmV/2VeU8kiJXpdsSny0rseiWvM+/bh6jyyUNOxeP8Xh8bP5EVngXv4GqTE5L3nhRhFOHrvqYWcNz\nfUmcie7pQkA6F4/z58+zvb2dO007WZJdJXPd96RiZNZcZip7JHrr1KUSzbFqAelcPKB6OF0l+dRU\nuXAbN2zZSsRkpNR3QRTdscpCss4TpmVYtedHlFBsen3PiMh5TSIgmmaVSdReRB5ZZP0yt91tCYKg\nteTiujqvdUHU5dva2srs+i0r/PHRwLKs0gM3jVUlUXsbeWjh3GZYx1qOqDYnPsSfV9xX95d4MjnN\nwcFB5e9hX763bSfUeyseeVTNeVSlzYRkkTVBGlEpfETcMqHofX0ano0X5NWt5zlx4sTKbuZ1iBDb\nFJBed1uSRDUaVW6IrHUqit7TJm3e0Kuesl6WNNvLrteDLaLr7kdTtJVE7Z14JHMCTSQro/eXFZG+\nTX1ehihkj1ZFj/JFq169vKmJeYeHhxdFH0U3ed2+f50fnj7Sli9M5+IR/yKkfQma/MNVWb5vHYh/\n8ZPh66p/+dMWJKqbj1mlsB8eHl4UyQ0x6d1GErXznEey+jNJk/kHFVd1h6b/94MmcyCdi0cS3eDN\nEf+FTi5yPBqNBtlt6QNDPe+IpgSkVLfFzF7DwtzpBXf/VTPbAR4BrgHOAreHpk9Le9X2oZ+5LgIW\nH6JMfmHiK6uvgia7LV3Tht3kqmnCF6Zs5PEBFiuiR3f0XcAT7j4BngwfJ71qbwU+Ga68Xppk5WWT\n/duyx1qXnEefSA5R920YuQpFw+1D+PEJF6te6otexm7yKuDtwH8Ffi98+jbgxnD7QeBLLASkca/a\nJpM8mzYLNYriou0uR1tgvZYEqHMtZSZBDoky3ZY/Av4AiN95u+4+DbenwG64vXZetUMnHm2t0807\nRMp8/mXqdPryI5grHmb2K8C+uz8detRehLsHZpYnpbkyGwTBCGA0Gp0HtrNeb4KyoVqTbUK/vVXV\n3vDaa6LL0QRFkcfPA7eZ2duB1wKXmdlDwNTMLnf3F83sCmA/3L+WVy1AEAQ71U69Hk0LgxBd0Ifv\ncW4y093vcfer3f1a4N3AF9z9N1h40p4JdzsDPBpur51XrRAinap1HlGo9DHgFjPbA94aPsbdnwMi\nr9rHWQOvWiFEOn3wqhVCDJDeVZgKIYaBxEMIUQuJhxCiFhIPIUQtJB5CiFpIPIQQtZB4CCFq8f8B\nJ5epf1mO5iUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11583af10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n = 20\n",
    "ex = np.ones(n);\n",
    "lp1 = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr'); \n",
    "e = sp.sparse.eye(n)\n",
    "A = sp.sparse.kron(lp1, e) + sp.sparse.kron(e, lp1)\n",
    "A = csc_matrix(A)\n",
    "T = scipy.sparse.linalg.spilu(A)\n",
    "plt.spy(T.L, marker='.', color='k', markersize=8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "For two-dimensional case the number of nonzeros in the L factor grows as $\\mathcal{O}(N^{3/2})$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Sparse matrices and graph ordering\n",
    "The number of non-zeros in LU decomposition has a deep connection to the graph theory.\n",
    "\n",
    "(I.e., there is an edge between $(i, j)$ if $a_{ij} \\ne 0$.\n",
    "\n",
    "``networkx package`` can be used to visualize graphs, given only the adjacency matrix. \n",
    "\n",
    "It may even recover to some extend the graph structure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeIAAAFECAYAAAD7rD1QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XdYFNf7NvAb6XUpuzTpGhQUsaNRQYy9VxQQWzAa/UWM\nJka/RhNNrJHYu8aGqFiiidhj7BWNYu8KRlAEAens7v3+YdhXVBQQWMTzua653N2ZOfOMwD5zZk7R\nIEkIgiAIgqAWldQdgCAIgiB8zEQiFgRBEAQ1EolYEARBENRIJGJBEARBUCORiAVBEARBjUQiFgRB\nEAQ1EolYEARBENRIJGJBEARBUCORiAVBEARBjUQiFgRBEAQ1EolYEARBENRIJGJBEARBUCORiAVB\nEARBjUQiFgRBEAQ1EolYEARBENRIJGJBEARBUCORiAVBEARBjUQiFgRBEAQ1EolYEARBENRIJGJB\nEARBUCORiAVBEARBjUQiFgRBEAQ1EolYEARBENRIJGJBEARBUCORiAVBEARBjUQiFgRBEAQ1EolY\nEARBENRIJGJBEARBUCORiAVBEARBjUQiFgRBEAQ1EolYEARBENRIJGJBEARBUCORiAVBEARBjUQi\nFgRBEAQ1EolYEARBENRIJGJBEARBUCORiAVBEARBjUQiFgRBEAQ1EolYEARBENRIJGJBEARBUCOR\niAVBEARBjUQiFgRBEAQ1EolYEARBENRIJGJBEARBUCORiAVBEARBjUQiFgRBEAQ1EolYEARBENRI\nJGJBEARBUCORiAVBEARBjUQiFgRBEAQ1EolYED5A8TEx2BwYiM2BgXgcG6v2cgRBKD4tdQcgCB+L\n+JgYHB03DgDgPX06rOzti1xGVlYW7ty5gz1ffonRR48CAEKOHYNR375wdnaGRCKBsbExTExM8v1r\nbGwMbW3tfGXl5OTgwOjR6LtlCwBgM4Be69e/30m+oiTOWRAqOg2SVHcQgvAxWNOtG/pv3w4A+D87\nOzxp3Bi6urrQ0dGBjo4OdHV1Ve+zsrKQmpqKxMREJCQkICEhAY8fP0ZqairMzMzQLjsba1JTAQB+\nOjrYTiI3Nxc6OjrQ1taGhoYGlEolFAoFFAoF5HI5NDQ0oKGhAQBQKpXQ0NBAFwC///cVEGRkhKdN\nm8LF3h6VL16EtrY27L74AvYuLq8ldx0dnUKd8+bAQPQKDwcAbOjdG/4bN4rkLAivEIlYEIrh5WTS\ndOpUGFtYID09HWlpaUhPT0dsbCwuXLiAU6dO4eLFi4iPj0e73Fzs+G9/Px0d7NXTQ25uripZKpVK\nvOvPUUNDAyShBaD9f5/tAiB/R7x5SVhDQwOampqoVKnSi9dKJVrm5kKpVGK/piaySXRUKlVxdtPQ\nwB5dXWhqaoIkFAoFcnJyoKGhAX19fejp6UFfXx8GBgaqxcjISLVYnjqFX27eBAB0AXDQyAg9tLSw\nOjkZALDJ3x+9/0vUgvCxEolYEIpAqVTi4sWL2D5oECZduADgRYLZ+1KyysnJgVKpBABUqlQJEokE\nFhYWMNbVhcfDh8jJycFOpRJyDQ3Y2trC3t4ezs7OqFq1KqpVqwY3NzdYWlrCyMgIurq6qlpsYWJ7\n+vQpHjx4gIcPH+LRo0d49OgRbt++jXv37iE+Ph7JycnIyMiAQqEAAGhra0NPTw9yuRw6OjpwdnZG\n3dhYrExMBAD00NTEdlJ1Pi/HoqWlparN6+joQEtLS3WhkBcPcnPRLDUVShJ7KlVCplyODgqFKtF3\nAXDIxAQ1a9aEr68vevTogdq1axf6nAWhIhCJWBDeISEhAfv378eePXuwZ88e6OjooG12NlY8fQrg\nRbI6aWmJzMxMpKenw9LSEoaGhnj+/DmSkpJQpUoVuLq6qpZq1arB1dUVMpnsnQknJycHCQkJePLk\nCR4/fownT57kW179TFdXF5aWlrCysoKlpWW+5eXPDAwM8PDhQ1y+fBnR0dE4d+4czp49C6lUijo1\na6Le48cwMzOD16RJ+NTHB6mpqTh06BAiIyNx7Ngx3Lt3D+bm5tDV1UV2djaSk5OhqakJY2Nj6Onp\noVKlSpDL5cjIyEBKSgoMDAxUxzc3MsInt2+DJE6bm+PGvXt49uwZNDQ0oFAooKGhAUtLS1Vy7tmz\nJ6pVqyZuaQsVlkjEgvAKuVyO06dPY9u2bfjzzz8RExMDc3NzZGVlITMzE1ZWVmB2Nho+fYpcuRy7\nAOibmMDV1RX16tWDm5ubKtk6OjpCU1NTVTZJJCcnF5hIX/0sLS0NMpmsUMlVJpNBX1+/2Oe9fft2\nhISEYObMmbhx4waio6Nx8eJF/Pvvv3Bzc0OtWrVQq1YteHp6wsXFBaePHUP09Ol48uQJNj1/jqpu\nbnBzc4O1tTV0dHQQFxeHu3fv4s6dO0hMTISNjQ2srKxgamoKIyMjaGm9aCuanZ2NhIQEPHjwAAkJ\nCcjOzgbw4m4CAFVtvKeWFjbLX9yE3xwQUOINywRBXUQiFgQA0dHRWLVqFfbu3YubN2+qbrGamJhA\nU1MTqampqoSgr68PNzc3+Pj4oEuXLnB3d0daWtpba6p5nyUkJMDAwKDAZPrqZ6ampqqEVBb8/f1R\nuXJlzJo1S/VZWloarly5gosXLyI6Olq1tMnMxKacHADArKZN4fL117hy5QqOHDmCU6dOoUqVKvD2\n9oa3tzfq16+PjIwM3L17V5Wc8/69d+8eTExMUKVKFbi4uMDR0RF6enqIjY1VXRDI5XK0y8lBxH/H\nE4lYqEhEIhY+Ci/f1nT/9lucjIpCREQELly4gJSEBLQFoAFgNwDT/5KhUqnE48ePkZGRgapVq8LB\nwQESiQSZmZn5EmxGRsZribSgBGtpaQldXd1ye5s1ISEBHh4e+OOPP9CwYcN86x4+fIj9+/erltZZ\nWViflgYA+NLGBpuysmBiYgIvLy/Ur18fpqamePLkCU6ePIljx47BxsZGlZi9vb1h/985K5VKxMfH\nq5Lzy4n67t27SE5ORuXKlaGnqYnajx4hKzMTF+3s4NOyJXybNIHmvn3Q0tQsV/+PglAUIhELH4VV\n3bph4H9dh7oA+AOArq4ubG1t0So9HUufPAEAdAVUDYm0tLRUrYLzuhfp6emp/tXX14e+vj50dXVR\nqVKlfEtey+SCPtfcuRNz7t8HAHxXrRqy2rRRNXIiqVre9b44+7zrfUxMDK5fvw5vb28kJSWpavdZ\nWVmq2+RSqRRaJJyvXQMA3KlWDZp6ekhLS0NKSgqePXuGlJQUpKWlwdDQEBKJBDo6OlAqlUhPT1c9\nU5ZIJDA1NYWJiQn09PTeGI9CoUBWVpZqyczMREZGBrKysvI1/BK1ZOFDJQb0ECqcvNpmTk4Orjo7\nY+PWrfC4excD/1uvq6MDa3NzVKpUCQ8ePED8f7ecAcCjZk18OWsWjIyMVF2K3rQUZ93Tp09x9uxZ\nnD17Fo3/S/wAYGhoCBtn53z9fF/ublTU98UtgyQePHgAbW1tXL16FXv27IGrqyu8vLzg6emJKlWq\nQFNTs0jl5+Tk4M6dO7h+/bpqycrKQu3atWFjYwNNTU0kJibi6tWrUCqV8PT0RJ06deDp6QknJ6d3\nHu/chAnAvn0l9asjCGohasRChbPR3x99Nm4E8P9rv1oAOgAgXu93a2lmhs6amtDS1kZ83bqQ2tjA\nxMQEEokEBgYGqv6yecvL7wtal9ftKCkpCVu3bsX69etx6dIldO/eHQEBAajm7Izj48cDUO+t6bzb\nzfv27cOBAwcglUrRunVr1K1bF6NGjcLBgwfh6elZoseMi4vD6dOnVUtUVBSsra1Ro0YNGBkZISUl\nBVeuXEFKSgqaNWsGb29v+Pj4wNPTM1/DNwB4HBuLI2PHAihft/gFoShEIhYqjLi4OCxZsgTXZ8zA\npv9a3vrr6aka+FhZWUEmk+Hhw4fQ1taGh4cH0tPTcf36daSlpcHU1BTa2trIzMxEZmYmsrOzoa2t\nrRqtSktLK99gGMCLWqRSqYRcLlct2dnZyM3NVTX40tPTg7GxMczMzGBoaPjOJF6UhP+m92/re5ye\nno7Dhw9j37592LdvH548eYKWLVuiVatWaNWqFRwcHFTbzp4xA5dmzkSbNm3QfMaMUktyCoUCV69e\nxenTp3Hq1CmcPn0a9+7dg5ubG2QyGXJzc3H//n08efIETZo0ydcArLAjfAlCeSYSsfDBO3PmDKZP\nn47du3fDyMgI6c+eobVCAT09PeS2aoXBw4ZBLpdj8uTJyM7Oxk8//YROnTrlS1bnz5/HjBkzcPDg\nQXz55Zf46quvYGFhgefPnyM5OVm1pKSk5HuftyQlJeHevXt4+PCh6vmnXC6Hvr6+amhIIyMjGBoa\nwsDAIF/S1NbWViX6vGfJeXJycl57PlrQ+5y0NPhkZoJKJQ7o6EDbwAC6/w00kpubq3q2amJiAktL\nS1SuXBnW1tb54nk5oWdu3Ijvo6IAAGE9e6Lv5s1l9jNNTU1FVFSUqtZ86tQpKJVKODs7Q0tLC0+e\nPEFcXBy8vLxUidnLywsGBgZlFqMglBSRiIUPUnp6OiZOnIjffvtN1bVIX18f9evXh7+/P/z8/GBh\nYYFjx45h/PjxiI+Px+TJk9GrV6+3dge6ffs2Zs2ahYiICAQGBmL06NFwcnJ647YkcerUKaxfvx6b\nN2+Gi4sLAgIC4Ofn96KvMYm0tLRCJfKCttHW1oapqWm+Ja+B06tL/NKlCPn7bwDAD56e+MfBAceP\nH4eZmRkaN26M+vXrw93dHZUqVXprQs97r7lzJ365dQsA0A3A/dq10aFDB7Rq1QqNGzcu09poXiOy\nvBrz6dOnceHCBchkMkgkEqSlpSEuLg516tSBj48PvL298emnn8LExKTMYhSE4hKJWPhgXL9+HaGh\nodi2bRuSkpIAALa2tvDz88MXX3yB6tWrq2q5UVFRmDBhAq5du4YffvgBQUFBqgEkCiMuLg5z587F\n8uXL0b59e4wZMwYeHh4AgCtXriA8PBzh4eHQ09NDYGAg/P39UaVKlRI9X5LIyMh4a6JOTk5GQkIC\nbt++Dfvz5xGekQHgReKM1NaGXC5XTdSQ99w773VBS942uRkZuLtgASpVqoRzlpZYs2EDGjVqhISE\nBNy6dQtNmzZV3dKuUaNGmQ9LmZOTg+joaFWN+dSpU3j06FG+rmfVqlVDy5Yt4e3tjaZNm8LCwqJM\nYxSEwhCJWCi30tPTsWrVKmzcuBHnzp1DVlYWNDQ04OjoiG+++QZDhgx5LblevnwZEydOxOnTpzF+\n/HgEBwe/V80tJSUFS5YsQWhoKCwsLJCbm4usrCz4+/sjMDAQnp6eZZ6AlEol/vnnH9Vz3rNnz6JB\ngwZo0rAh7KKjoampiVtVq2L7rl2Qy+Xo1q0b2rZtCxsbG6SmpiIlJQWpqalvXV7dJj09/cWt/txc\naGlpwcXFBZqamnj+/DkSEhKgVCpRpUoVuLu7o27duqo+168meWNj4yJdEBVVYmIizpw5g9OnT+Pk\nyZM4efIkNDU1oauri+Tk5Bfd1Vq1wmeffYZqzs64OWcOANHQS1AvkYgFtXl1UAtLOzscOnQIy5Yt\nw6FDh/D48WNVYyktLS2MGDECI0aMeGOt5vbt2/jhhx+wf/9+jBkzBsOGDXvv54WJiYnYvHkzwsPD\ncfXqVbi7u+P27dtwdnbGuHHj0KFDhzJLwrGxsflaN1taWqJVq1Zo3bo1fHx8YGRk9No+JHH+/Hls\n2LABmzZtgqmpKfz9/dGnTx+4uLgU6fhKpRJpaWl49uwZli1bhoULF6Jz585o166dasSsixcv4ubN\nm3jw4AH09PRgZmYGfX19aGhoIC0tDampqXj+/Dn09PSKVUN/eclIScGpCRMAvDmJkoRcLkdubi5u\n3Lihup19/Phx3Lt3D/r6+vgsLU01BeSMxo3RZ8MGODo6FvMnJAjFJxKxoDYvz1Xb38QE4RkZUCqV\ncHBwgLW1Ne7du4dPPvkEISEh6Nq16xtrUjExMfjpp5/w+++/IyQkBCEhIe/1XDA9PR1//PEHwsPD\nceTIEbRr1w4BAQFo27YtdHR0IJfLsXXrVsyYMQO5ubn47rvv0Lt3b2hraxf7mG+Slpamat28f/9+\nVevm1q1bo1WrVqpRqQpLqVTi+PHj2LBhA7Zs2QIXFxfVs3QbG5sixxcbG4vhw4fjzp07WL58OT79\n9FPVOoVCgXPnzuHAgQPYv38/zp49i7p166pqou7u7sjIyCiw9l2YGrrPs2f4/b/j9dLRwV5dXdW8\ny3K5HEqlUtXKXUtLK9+S94zcNzUV2/7rQz7CwQER2dnQ1dXNN/qXq6urmAlKKH0UBDVQKBScVLcu\nCZAAv7S1ZWhoKAcNGkRTU1P279+f586dK3D/uLg4jhgxgubm5hw7diwTExOLHUtOTg4jIyMZEBBA\niUTCtm3bct26dUxNTS1wH6VSyb1799LX15eOjo6cP38+09PTix2DQqHg2bNnOWXKFDZv3pyGhoZs\n3rw5p06dyqioKCoUimKX/aqcnBzu3r2b/fr1o6mpKVu0aMHly5czKSmpSOUolUpGRETQxsaGw4YN\nY0pKyhu3S0tL465duzhq1Ch6eHjQ1NSUXbt25YIFC3jjxg0qlcoin8OmgADV7860Ro2YnJzMtLQ0\nZmVlUS6Xv1ZmTk4O9+7dy8GDB1Mmk7F+/fqcMHYsV3TqxIiAAMbHxFCpVPLGjRtcvnw5g4KC6Ojo\nSCsrK/bq1Yvz58/nxYsXS/TnIAh5RCIWytxff/1FBwcH6mpo8FtXV4Z6e9O7cWPa2Nhw8uTJfPz4\ncYH7JiYm8rvvvqOZmRlDQkIYHx9frBgUCgWPHj3KL7/8klKplJ9++ikXLFjw1mMX5NSpU+zWrRst\nLS05efLkQl8UxMTEcOXKlezduzctLCzo5ubGkJAQRkZG8vnz50WOozgyMjK4ZcsW9ujRgyYmJuzU\nqRPXr19fpOMnJSUxODiYdnZ23L59+zu3j4uLY1hYGPv3709bW1s6ODjw888/54YNG/jkyZNCHTM+\nJoYRAQGc8emn9G3a9I3bZGdnMzIykgMHDqSFhQUbNWrEWbNm8d69e4U+t/v373Pt2rUMDg6mq6sr\nzc3N2blzZ86aNYtnzpxhbm4u4x48YERAgCqhC0JRiUQslJmLFy/S19eX+vr6dHJy4rfffktnZ2d6\neXlx/fr1zM7OLnDflJQUTpo0iRYWFhw8eDBjivmFFx0dzbFjx9LR0ZHu7u6cMmUK7969W9xTyufa\ntWscNGgQzc3NOWrUKJ47dSrfF/Tz58+5c+dOjhgxgtWrV6eFhQX79OnD3377jbGxsSUSw/tISUnh\nmjVr2LZtW5qYmLB3797cvn07s7KyCrX/oUOH6Orqyh49evDRo0eF2kepVPLq1aucO3cuO3bsSBMT\nE9apU4fffvst9+3bx4yMjLfun5WVRZlMxps3b5IkMzMzuWPHDgYFBdHMzIxNmzblnDlziv378qq4\nuDhGRERw+PDh9PDwoLGxMYdaW6tq5xEBASVyHOHjIhKxUOoePHjA/v3708TEhIaGhqxbty4lEgkD\nAgJ46tSpt+6bnp7OmTNnUiaTMTAwkLdu3Sry8e/du8epU6eyZs2atLe355gxY3jx4sVi3RItjNjY\nWI4aNYq9dHRUX9DBMhkNDQ3p6+tbKrebS9qTJ0+4aNEient709zcnIMGDeL+/fspl8vful9mZibH\njx9PmUzGpUuXFvkcc3JyePToUU6cOJGNGzemoaEhP/vsM06fPp3nzp17Y3lff/01O3fuTH9/f0ok\nEjZv3pwLFizgv//+W6RjF0diYiJ/9fZW/ZwXtGxZ6scUKh6RiIVSk5SUxG+++YampqZ0cnKivr4+\nzczMOGHChHd+SWZlZXHBggW0tbVl9+7defny5SIdOyEhgQsXLmSTJk1oYWHBoUOH8siRI2WW/HJy\ncji+Zk3VF3S3SpXYsGFDjho1ilu2bCl0jbE8iI2N5axZs1ivXj1aWVnxq6++4okTJ956IRMdHU0v\nLy82a9aM165dK/axk5OTuX37dg4fPpzVqlWjhYUF/fz8OH/+fC5YsIC9evWikZERtbW1OX/+/GI/\nqngf8TExDDIx4c8NGtBWKuX333/PnJycMo9D+HCJRCyUuMzMTM6cOZMWFhasUaMGtbS0KJVKuXTp\nUmZmZr5139zcXK5cuZKOjo5s164do6KiCn3c58+fMywsjO3bt6dEIqG/vz///PPPt97yLg0nT56k\nh4cHW3p787fOnRkREMA716/z77//5tSpU9mxY0eam5vTycmJ/v7+nD9/Ps+dO8fc3NwyjbM4bty4\nwUmTJrF69ep0cnLi2LFjC7y7IJfLOW/ePFpYWHDy5Mnv/XNISUnh3LlzWbduXWpra1NbW5uWlpYc\nMGAAa9euzSVLlrxX+e/D1dWV169fZ1xcHFu3bs1GjRrxzp07aotH+LCIRCyUGLlcztWrV9PW1pbO\nzs7U0dGhnp4eZ8yY8c7bwAqFghs2bKCrqyt9fHx49OjRQh0zOzubf/75J/v06UMTExO2b9+eYWFh\nZdbY6WUpKSkcPnw4ra2tGR4e/tZzViqVvH79OletWsXBgwezRo0aNDIyYvPmzfm///2PO3fu5NOn\nT8sw+qJRKpX8559/OGbMGDo4ONDNzY2TJ09+46ODmJgYduzYkTVq1OCJEyeKdJxnz55xzZo17NSp\nE42NjdmxY0euXr2aSUlJVCgU/Oeffzhz5kx6enqyUqVK9PLy4vfff8/Dhw+X6QVY1apVVc+pFQoF\nf/31V0qlUoaFhZVZDMKHSyRi4b0plUpGRkbS2dmZZmZmNDIyoqWlJVu1avXOVshKpZLbt2+nh4cH\nGzZsyH379hUqaR8+fJhDhgyhhYUFmzRpwkWLFhW6xW1JUyqV3Lp1KytXrszg4OBid6V69uwZ9+zZ\nw4kTJ7Jly5Y0MTFhtWrVOHDgQC5btoyXL18ul8+VFQoFjx07xuHDh9PS0pL169dnaGgoHz58qNpG\nqVRy06ZNtLGx4fDhwwvs6kSST58+5cqVK9muXTuamJiwa9euDAsLY3JycoH75OTk0MbGhsuWLePY\nsWNZr149Ghsbs3379pw9ezYvXbpUam0CSNLFxeW1i5Dz58+zevXq7Nu371vPVxBEIhaK7OXuGjs2\nb6arqyt1dHRob29PPz8/SqVSrly58p01wn379rFhw4asVasWd+zY8c7tL1y4wDFjxtDe3p4eHh6c\nNm1akbqilIaYmBh27tyZ1atX5+HDh0u0bLlczosXL3Lx4sUMCgpi1apVaWpqyjZt2nDSpEncv39/\nufuCz83N5b59+zhw4ECamZnRx8eHS5YsUdXuk5KS+Pnnn9Pe3p47duxQ7ff48WMuXbqUrVq1oomJ\nCXv27MmNGze+tS/3qyZOnMjhw4er3ickJDAiIoKDBw+mk5MTbWxsGBQUxLVr15Z4Qy5nZ+c33opO\nT0/nkCFD6OzszJMnT5boMYWKQyRiocgiXhpMoQvA2rVrMzw8nL6+vvz000/f+Wzs6NGj9Pb2pqur\nKzds2PDWWt7du3c5ZcoUuru708HBgWPHjmV0dHRJn1KRyeVyzpkzhxYWFpw0aVKhu/i8r8ePH3P7\n9u0cM2YMmzVrRkNDQ9aqVYtDhw7l2rVreevWrVKt+RVFZmYmf//9d/r5+akeG+QNlPL333/T2dmZ\ntWvXZpMmTSiRSNinTx9u2bKFaWlpxTpebGwszczM3vhYQqlU8tatW1y8eDG7d+9OMzMz1qhRgyNH\njuTOnTvf+1GGo6PjWy8Kt23bRktLS/7000/vbHkufHxEIhaKbEnbtqpEPM7dnStWrKBMJuOUKVPe\n+iUTFRXFtm3b0tHRkb/99luBjZMeP37M+fPns3HjxpRKpRw2bBiPHTtWbm7Lnj9/nvXr16ePjw+v\nX7+u1liys7N5+vRpzp49m35+fqxcuTJlMhm7dOnC6dOn88iRI+/si1sWUlNTGRYWxhYtWlBXV5dS\nqZQGBgZ0d3ensbEx58+fXyI/3y5dunD58uXv3E4ul/P06dP8+eef6ePjQ0NDQ3p7e3Py5Mk8efJk\nkRvO2dvb88GDB2/d5uHDh/T19WWzZs3eua3wcRGJWCiSS5cu0UYq5axmzTineXM62NiwUqVKDAgI\nKHBgjMuXL7N79+60sbHhggUL3lh7TE1N5dq1a9m2bVtKJBIGBgYyMjKyXHUDSUtL4+jRoymTyfjb\nb7+Vm5rnq2JiYrhp0yaGhISwQYMGNDAwYMOGDRkSEsJNmzaV+eAh9+/fZ2hoKBs3bkxzc3P6+/vz\n//7v/+jt7U1TU1N27tyZrq6ubNas2Xtf2OzevZv16tUr8n55w3B+/fXXqmE4u3XrxoULF/LmzZvv\n/FlXrly5UIOGyOVyTp8+nTKZjBEREUWOU6iYRCIWCu369eu0sbHhhg0buH//ftrZ2XHEiBG8fv06\nv/nmG5qbm7N79+48evSo6lZgYGAgZTIZf/nll3xjMcc9eMANffrwV29vdu3QgSYmJuzYsSPDw8OL\nfWuyNEVGRtLR0ZGBgYHFGgZTndLT03nkyBFOnz6dnTt3pkwmo52dHf38/DhnzhyeOXOmxC947ty5\nwxkzZrBBgwaUSqUMDg7mnj17XmvJ/PDhQ/76669s0KABjYyMaKCtzVEuLgzv3btYw0UqFAo6OTnx\n7Nmz7xV/XFwc161bx379+tHW1paOjo78/PPPuXHjRiYkJLy2va2tbb7Gae9y5swZVq1alYMGDVJL\nC3+hfBGJWCiUO3fu0M7OjsuWLWNISAjt7Oy4b9++fNukpqZy/vz5dHR0pFQqpZGRESdOnJivQVFa\nWhq3bt3KYZUrq25v/9ygQbntqhMXF0c/Pz+6uLi8dr4fqryLpDVr1nDIkCGsVasWDQ0N2axZM373\n3Xfcvn17sS42bty4wSlTprBOnTq0tLTk0KFDeeDAgULf5r19+zbHVq+u+r34ysGhWHFMmzaNgwYN\nKvJ+BVEqlbxy5QrnzJmTbxjOMWPGcP/+/bx74wZ76+lxdbduRbp4eP78OQcOHMhPPvnkvS8chA+b\nSMTCO8Xpfa1MAAAgAElEQVTExNDJyYn/+9//6O7uTj8/vzd20YmPj+eIESNoZmbGHj16sFmzZrS1\nteX48eM5f/58dunShSYmJmzUqBF7amuX6/F5FQoFly5dSqlUyrFjx77XzEofgpSUFO7bt48//vgj\n27RpQ4lEwqpVqzIoKIiLFy/mxYsXKZfLX5vg4MqVK5w0aRI9PDxUXZP+/vvvYjdIerkhYG89PWpo\naLBJkyY8f/58oct4/PgxTU1N+ezZs2LF8C7Z2dk8cuSIahjOHlpa7/W7vHHjRspkMs6YMaPctIMQ\nypZIxMJbPXr0iFWqVGGHDh0ok8kYFhb22vOyxMREjh07lmZmZhwxYgTj4uIYExPDefPmsUGDBqpR\nkJo3b85Vq1bRysqKSxcsKLcz1ly5coVNmjRho0aNykULbXVQKBS8dOkSly5dygEDBtDV1ZUmJiYc\n8tIEBwFGRqxcuTJDQkJ49OjREkkiebMq5f1e7Nmzh25ubtTQ0GDVqlUZERFRqGfzffr04bx58947\nnsII69VL9X/ySwEzQb3L/fv32bRpU3722WdFusUtVAwiEQsFevLkCatWrUpHR0f6+vq+1tLz1RmR\nDhw4wClTprB+/fo0Nzdnv379+PvvvzM9PZ3x8fEMDg5mpUqV6Onpyd27d5e7q//MzEx+//33lEql\nXLhwoehm8oqEhARObdjw/ydiQ0Pq6urS0tKSjRo1YkBAAL///nv+9ttvPHToEGNiYkrs//DixYts\n2rQpNTU1aW5uzqlTp761LcGhQ4fo7u5eJg3q4mNi6Keryzm+vnS2s+Po0aOLNVxpbm4uJ0+eTCsr\nq0JNJylUHBokCUF4RWJiIurUqYPExET89NNPGDlyJCpVqgQAyMjIwKJFi/DLL7+gbt26cHR0xKFD\nh5CWloauXbuiW7du8Pb2hra2tqq8o0ePokePHli2bBmSk5Mxe/Zs5OTkYOTIkQgKCoKBgYG6ThUA\ncPDgQQwdOhS1atXCvHnzYGtrq9Z4yqOrV6+i3Wef4Qs7O7i6usJ7+nTIKldGfHw87t69+9py7949\nJCYmwtHRES4uLq8tzs7OMDExKVIMsbGxGDVqFHbs2IFKlSqhb9++mDBhAhwdHfNtRxLu7u5YtmwZ\nmjVrVpL/DW9kZWWF6OhoaGlpITAwEFlZWdi0aROsrKyKXNaJEycQGBiItm3bIjQ0VO1/G0IZUPOF\ngFAO3b59m6amppTJZDy4Z4/qVmHMnTucO3cuLSws6OLiQktLS1arVo1jx47l6dOnC6zhHjhwgDKZ\njPv371d9plQqefDgQXbq1IlSqZTjxo1Tyy25hIQE9u/fn/b29vzjjz/K/PgfinPnztHa2ppr164t\n0n4ZGRm8evUqd+7cyXnz5nHkyJHs3Lkza9asSQMDA1pYWLBBgwbs3bs3x40bx+XLl/Ovv/7ivXv3\n3lqrfPr0KUeMGEE9PT3q6OiwRYsWPHz4cL4a8Jw5cxhQRu0PLC0tGRcXR/JFF6UJEybQzs6uyGNr\n50lOTmZAQADd3Nx44cKFkgxVKIdEIhby2bJlC7W1tVmnTh1mZmbmazzTXVOT2tradHNz45QpU3j1\n6tV3lhcZGUmZTMYjR44UuM3Nmzf51Vdf0czMjAEBATxz5kxJntIbKZVKrlmzhlZWVhw5cmSRhlL8\n2Bw9epQymYzbtm0r0XKVSiXj4+N54sQJhoWFcfLkyRwwYAB9fHxob29PHR0dVqlSha1ateKQIUM4\nY8YMRkREMCoqiklJSSRftDyeOnUqJRIJDQwMWLVqVa5cuZKZmZlMSkqiqalpmYxBbmVlpUrEef74\n4w/KZDIuWLCg2LfI161bR6lUytmzZ5e7RzlCyRGJWCD5oltRcHAw9fT02K5dOyoUihfPgOvWVSXi\nMdWqFWrQgjx5w/oVdozdZ8+ecdasWXRwcOCnn37KzZs3l8rUgDdv3uRnn33GOnXqiG4j77B3715K\npVLu3bu3zI+dlZXFGzducPfu3Vy4cCFHjx7Nbt260dPTk0ZGRjQ1NWXdunXZs2dPjh49mgEBATQ3\nN6exsTElEglH/d//8SsHB/5Yu3apNwh8UyImyVu3brFWrVrs27dvsVve3759m15eXmzbtq1a5lsW\nSp9IxAJPnTrFqlWrsnLlyuzZsycfPHjAb7/9lubm5qxTowb9DQy4tkePIn2ZbdiwgVZWVjx37lyR\n48nNzeXmzZvZpEkTOjo6ctasWSXSFSU7O5s///wzLSwsGBoa+kHM/6tO27Zto0wmK/SUlGVJqVQy\nISGBp0+f5oYNGzhlyhR+/vnnbN68OWUyGTU0NNjlvwvIsugiZ2VlxUePHr1xXXp6Ovv27ctatWq9\ncZrIwsjJyeH48eNpY2PDyMjI9wlVKIdEIv6I5eTk8IcffqClpSUbNGjAzz77jH379qWZmRlDQkK4\nePFi2traFjh0ZUFWr15NGxsbXrp06b1jPH36NP39/WlmZsavvvqq2F9kx44dY40aNdi+fXu1z9j0\nIVi7di2tra2LdSFVHmRnZ3Oql9f/HzSmYUNmZmaW2vGsra0LTMTkiwuHBQsWUCaTvVdbhEOHDtHB\nwYEjRowo1fMRypZIxB+pGzdusEGDBmzdurVqcgUbGxtOmzaNSUlJPHLkCGUyWZEGUiDJxYsX087O\njteuXSvReGNjYzl27FhKpVJ27tyZBw8eLNRzt2fPnnHIkCG0sbHhpk2byu340OXJwoULaWdnV6g2\nAOVZfEwMRzo5cYyrK1v7+tLGxoYzZ84slfYA70rEeY4fP047Ozt+//33xe7alZSUxJ49e9LDw4OX\nL18uVhlC+SIS8UdGqVRy0aJFNDc3Z1BQEM3MzGhoaMhly5apJmO4fPkyLS0t87VyLozZs2fT0dGR\nt2/fLo3QSb64zbdkyRJWr16dnp6eXLVq1RsnkVAqlYyIiKCtrS2HDBlSaqMsVTTTpk2ji4tLke+C\nlEepqamUSCSqYTIvXLjAPn36UCqV8vvvv3/jmNHFZW1tXeg5juPj4+nj48M2bdoUe2hXpVLJFStW\nUCqVctGiReIC8wMnEvFH5NGjR2zZsiXt7OxoZWXFypUrs0aNGvkGnY+NjaWDgwPDwsKKVPbUqVNZ\npUqVMpveTaFQcPfu3WzdujWtrKz4448/MjoqihEBAfytSxe2adGC7u7uPHbsWJnE86FTKpUcN24c\n3dzcKszITqtXr2anTp1e+/zWrVv84osvVI9gitIAsSA2NjaFTsTki3YQ33zzDZ2cnBgVFVXs416/\nfp1169Zl+5YtubZHj3I5Up3wbiIRV3B5YwNPadiQxnp61NfXZ58+fRgQEEAvL698EzI8e/aMNWvW\n5MyZMwtdvlKp5MSJE1m9evUifRGVpMuXL3Pw4MH5xq8eX7PmazP9CG+mUCg4fPhw1q1bt0Rrierm\n6+vLLVu2FLj+33//5ejRo2lmZsaBAwe+1xSMRU3EeTZv3kypVMqVK1cW+9jZ2dkc6+5ersduF95O\nJOIK7uV+wCMcHfngwQP+73//Y+3atVV9MckXwzv6+PgwJCSk0Le5lEolv/32W3p4eJSLqQHXdO+e\nr8+zl5cXx44dy3379lX4SRuKKzc3l/369WPTpk2ZnJys7nBKzP3792lhYfHGxxavSkxM5KRJkyiT\nydijR49i1VBtbGyKfSfh6tWrrF69OgcPHlzsBlgv/50vbd++WGUI6lNJ3SN7CaXrxs2bqtdNmzTB\nunXrsH37duzbtw9mZmYAAKVSiX79+sHKygq//vorNDQ03lmuUqlESEgIDh48iL///huWlpaldg6F\nZdG3LwaYmmJzQAB+vX4d06dPh7a2NiZPngxLS0t4e3vjhx9+wOHDh5Gdna3ucNUuOzsbvXv3Rnx8\nPPbu3QuJRKLukEpMWFgY/Pz8oKur+85tzc3NMXHiRNy9exdNmjRBly5d0KZNGxw6dAgs5AjAhfmb\nKYibmxvOnDmDpKQkNGvWDDExMUUuw3v6dIQ2bYoBpqZYGhMDuVxe7HiEsifGmq7AUlNTYWNhgRFV\nq6Ju3bq47uKCsIgIHD58GNbW1gBejMk7cuRIREdHY8+ePYX64lIqlRg6dCguX76MXbt2wdTUtLRP\npVB+/PFHZGRkYObMma+tS0tLw/Hjx/H333/j4MGDuHbtGry8vODr64sWLVqgfv36+cbGrugyMjLQ\nvXt3GBoaIjw8vFA/9w8FSVSvXh1r1qxBo0aNirx/dnY2wsLCMGPGDFhYWGDcuHHo2LGjaqz1N6lc\nuTLOnDmDypUrv1fcoaGhmDVrFsLCwtCyZcsi7T927Fjo6+vj6NGjaNWqFb777rtixyKUMXVWx4XS\n1aJFC8pkMiqVSi5cuJDOzs6vNUyZMWMGPTw8Cn1bMjc3l0FBQfT29i53w0I2atSIBw4cKNS2ycnJ\n/OOPP/j111/T09OTxsbGbNeuHX/55RdGRUVV6JmXkpOT2bRpU/br169CDmpy8uRJurq6vndLYrlc\nzoiICNapU4c1a9ZkWFhYgf9ftra2JdbI7eDBg7S2tubUqVOLNKxlixYtGBkZybt379LCwuK9nnkL\nZUsk4goqMjKSlSpV4uHDh7ly5Ura29u/1iVl7dq1dHBwKPQXSE5ODnv16sVWrVqVu2euSUlJNDY2\nLvYztoSEBG7ZsoXDhw+nm5sbzczM2KVLF86dO5fR0dEVZpzfhIQE1q1bl8OHD68w5/SqL7/8kj//\n/HOJladUKrlnzx56e3vT2dmZixYteu33zNbWlrGxsSV2zNjYWDZq1IhdunQp1EWyQqGgRCJRjas9\nb948NmnSpML+jCsakYgroKysLJqamrJDhw5cv349bW1teePGjXzb7N27l5aWlrxy5Uqhy+zSpQs7\nduxYLkf02bx5M9u1a1di5T169Ijh4eEMDg5mlSpVKJPJ2KtXLy5evJg3btz4IPtt/vvvv3R3d+e4\nceM+yPgLIysrixYWFqXWje748ePs2LEjra2tOX36dFWvg8qVK5doIiZfnMuwYcP4ySefvHOUuhs3\nbtDR0VH1XqFQsEmTJpw3b16JxiSUDpGIK6ABAwZQT0+Pa9eupZWV1Wt/xOfOnSvSGMIZGRls27Yt\ne/ToUW67BAUHB3POnDmlVv6DBw+4evVq9uvXj3Z2drS1tWXfvn25cuXKD2LIzLt379LFxYXTpk1T\ndyilasuWLfT19S3140RHRzMgIIAWFhYc9X//x956elzdtWup9OFds2YNpVIpw8PDC9xm/fr17NGj\nR77Prl+/TgsLiwoxOEtFJxJxBRMVFUUtLS1++eWXtLS0fG2Iyjt37tDW1rbQU9o9f/6cvr6+DAgI\nKLfPE5VKJe3t7ctsSEalUslbt25x2bJl7NOnDy0tLenk5MRBgwZx3bp15W5AjKtXr9LOzo4LFy5U\ndyilrlOnTly9enWZHW/Xrl3sb2pa6n14//nnH7q4uHDEiBHMycl5bf3IkSM5ffr01z6fMWMGP/vs\nswp7B6SiEIm4AsnNzaW9vT3t7OwolUp56tSpfOufPHnCTz75hIsWLSpUecnJyWzSpAk///zzct14\n6dq1a7S3t1fbl41SqeSVK1c4f/58du/enebm5nR1deXQoUMZERFRJvPhFuTcuXO0trbmmjVr1BZD\nWXn8+DElEkmpNyJMTEzk/PnzWa9ePdrZ2fF/NWqoEnG4n1+pHTcpKYkdOnRgkyZNXhvXumnTpm9s\nqJibm8t69epx+fLlpRaX8P5EIq5Axo8fTy0tLUokEh4+fDjfurS0NDZs2JDjx48vVFmJiYls0KAB\nhw0bVu4bfMyZM4fBwcHqDkNFoVDwn3/+YWhoKDt27EgTExN6eHhwxIgR3L59e76BVErTsWPHKJPJ\nuHXr1jI5nrrNmTOHQUFBpVJ2bm4uIyMj2atXL0okEvr7+3Pfvn2Uy+WMj4lhREAARzo7s5+/f6kc\nP49CoeDkyZNpa2vLI0eOqGIzNDQscDz16OhoSqXScnenRvj/RCKuAOIePOCKjh3ZVUODBtrar03W\nkJubyw4dOnDAgAGFqjU+efKEnp6eHD169AdxS6t9+/aMiIhQdxgFys3N5alTpzht2jS2atWKRkZG\nrFevHr/55hvu2rWrVGpw+/bto1Qq5Z49e0q87PKqbt26RZ6o5F2uXr3KMWPG0MbGhl5eXlyyZEmB\nCS8lJYUuLi6FfuzzPnbv3k1LS0v+NGECl7Zvz75GRm99Pv3DDz+wY8eOH8Tf88dIJOIKYONLw9v9\n0rRpvnVKpZKDBg1iu3bt3vhs6VWPHj2iu7s7v//++w/ijzYrK4vGxsZlVsssCVlZWTxy5AgnTZpE\nHx8fGhoasnHjxhw/fjwPHDjAjIyM9yp/27ZtRWqMVxFcunSJdnZ2JfII5dmzZ1yyZAm9vLxoY2PD\nMWPGFLr9wYkTJ2hlZVUm467fvXuXA8zMCvV8Ojs7W9UXWih/RCKuAH6qX1/1xxhkYsLQ0FDV4P0T\nJkxg/fr1882wVJCYmBh+8sknJdoHs7QdOHCAjRo1UncY7yUjI4MHDhzg+PHj2bhxYxoaGtLHx4eT\nJk3ikSNHitRSfd26dbS2tua5c+dKMeLy59tvv+XYsWOLvb9cLue+ffvo7+9PiUTCnj17MjIyslgN\nFCdNmsRWrVqVySOd8N69C91Q7MyZM7S0tGR8fHypxyUUjUjEH7j79+9TZmrKOc2bs4++Prdv3syg\noCBKJBLWq1ePtra2jIuLe2c5d+/epbOzM0NDQ8sg6pIzZswYTpw4Ud1hlKjU1FTu2rWL3377LevV\nq0cjIyO2atWK06ZN4+nTpwtMDosWLaKdnV2h+4ZXFLm5ubSxsSlWq/mbN29y/PjxtLe3Z7169Th/\n/vxizxH8cjyffvopf/311/cqpzBOHDrEXjo69Dcw4Oa3dG/KM2bMGPbq1avU4xKKRiTiD5hSqWTr\n1q05ZcoUKpVKOjk58fLlyyRf9Cs0MTFh9erV6ezszJ9//rnAxho3btygvb39B9m9xdPTk8ePH1d3\nGKUqKSmJ27dv54gRI+jh4UGJRMKOHTsyNDSUB3bt4saAAP5QuzZd7O15584ddYdb5vbs2cMGDRoU\nevvU1FSuWLGCTZs2paWlJb/++mtevHixRGO6e/cuZTIZL1y4UKLlvmrw4MGcOHEid+zYwWrVqr1z\ntqmMjAy6urp+NA34PhQiEX/AVq9ezdq1a6ue/X7xxRcMDQ3l8ePHKZPJePbsWSqVSp49e5ZDhgyh\nmZkZO3XqxB07dqhqVVeuXKGtre17zYeqLnFxcTQ1NS23/ZtLy+PHjxkREcGhQ4eyr7Gx6tbkUBsb\nDhkyhOPHj+fs2bMZFhbGPXv2MCoqivfv3y93w5KWFH9/fy5YsOCt2ygUCh48eJD9+vWjRCJh165d\nuX379kK1myiudevW0d3d/b2f+Rfk33//pZmZmeoxVKdOnQr1WOno0aO0sbFhYmJiqcQlFJ2YfekD\nFR8fj1q1amHv3r2oU6cOAGDbtm349ddfcfv2baxevRpt27bNt09aWho2b96M5cuX4/79+2jfvj3+\n+OMPzJ49G4GBgeo4jfeSN6Xj1q1b1R2K2mwODESv8HAAQGjTpjAMDERCQgKePn2Kp0+fvvZaQ0MD\nUqkUMpkMUqn0tdevvrewsICWlpaaz7JgqampcHBwwJ07d2BhYfHa+nv37mHNmjVYs2YNjI2NMXDg\nQAQGBpbJtJ0kERgYCAsLC8yfP7/Eyx8zZgyys7Mxd+5cAMD9+/dRv359nD17Fs7Ozm/dd8SIEUhN\nTcXq1atLPC6h6EQi/kD16NED1apVw9SpU1WfXbt2DTVq1MDSpUsxePDgt+6/ceNGfP7559DU1ETD\nhg0RHByMbt26fVDT4fXt2xfe3t744osv1B2K2jyOjUVoy5YwNjbGF7//Dit7+wK3JYmMjIzXEvSb\nEnbe62fPnsHY2PidCfvl1yYmJu81P29RrFy5EpGRkdi2bZvqs/T0dGzduhWrVq3C5cuX4e/vjwED\nBqBOnTplFlee5ORk1K5dG4sWLUL79u1LrNxnz56hatWqOH/+PBwdHVWfT5s2DcePH8eff/751nNN\nS0uDh4cHFi1ahHbt2pVYXELxiET8Adq6dSvGjx+PCxcuQE9PDwCQkpICb29vJCcnY8WKFWjVqtVr\n+8XHxODouHF4kpCAn86dw/LVq9G6dWts374dy5cvx8WLFxEUFITg4GC4u7uX9WkViVKphI2NDU6f\nPg0nJyd1h6NWwcHB8PLyeufFV3EoFAokJye/M2G//DorK6tQte2XXxf3AtDHxwdff/01unTpguPH\nj2PVqlXYtm0bmjRpggEDBqBTp05qv7g8evQoevfujX/++QdWVlYlUubUqVNx48YNrFmzJt/nOTk5\n8PT0xLRp09C1a9e3lnHgwAF8/vnnuHTpEkxMTEokLqF4RCL+wCQlJaFmzZrYvHkzmjRpAuDFRObt\n2rWDu7s7ZDIZnj9/jlmzZr2278u3Mef6+iLk4MF86+/evYuVK1di1apVcHZ2RnBwMPz8/GBoaFj6\nJ1ZE//zzD/r06YMbN26oOxS169u3L9q0aYOgoCB1hwLgxe9jYWrbL7/W09MrVMLOe539/Dn2hoRg\n9549cB42DNt27oSWlhYGDhyIoKAg2NjYqPu/IZ+8C+edO3e+d608MzMTzs7O+Ouvv1CjRo3X1v/9\n998YMGAArly5AiMjo7eWFdS7N6zPn0fDhg3hPX36W++oCKVHJOIPTP/+/SGRSDBv3jwAL2qGgYGB\nyM3NxaZNm3D27Fl88cUXiI6OVu2TmZmJP/74A3+PHIkl8fEAgM0BAei1fv0bjyGXy7Fr1y6sWLEC\nx44dg5+fH4KDg1GvXr0yv7VXkOnTp+PRo0eq/4ePWY8ePeDv74+ePXuqO5RiIYnnz58XKmHnvfZJ\nTsbv/311fVu1KnqGhaFhw4bl5vfzVbm5ufj0008xYMAADB8+/L3KWrRoEfbu3YsdO3YUuE1QUBBs\nbW0xY8aMt5a1vlcvBG7ZAuDt3wlCKVNLEzGhWHbv3k0nJ6d8g3OMGjWKzZo1U80RLJfLaWZmxtjY\nWB46dIiDBg2imZkZW7duzf4BAfzS1pYRAQGFnq7t4cOH/Pnnn+ns7MzatWtzwYIFBQ7xV5Z8fX25\nc+dOdYdRLrRv3/6j+7/Y6O9f6jMelbSbN29SKpWquhgWR25uLp2cnHjixIm3bhcfH0+pVPrOeYwj\nXhqVb3nHjsWOS3g/IhF/IFJTU+ng4MC9e/eqPgsNDaW7u3u+4R2vX7/O6tWr08LCgh4eHvzll19U\nw+21bduWGzZsKNbxFQoF9+/fTz8/P0okEgYFBfHIkSNqGQbz+fPnNDIyKtRoYR8DX1/fN868U5HF\nx8Rwcdu27GdqWipzAJeWFStWsFatWu/s71uQ9evX09vbu1DbLly4kM2aNXvr32h8TAznt2zJrgDN\nDA0/iGFtKyKRiD8Qw4YN48CBA1Xvw8PDaW9vz5iYGCYkJHD+/Pls0KABra2t2bp1a7Zt2zbf/k+e\nPKFEImFaWtp7x5KQkMDQ0FC6ubnR1dWVM2fO5OPHj9+73ML6888/y2Ty9w9F48aNK/ygJm+Snp5O\nPT29D6ofuVKpZPfu3Tlq1Khi7evh4cFdu3YVanu5XM569eq9c37m1NRU6unpEQB79+5d5LiE9ycS\n8Qfg8OHDtLW15bXoaEYEBHBeixa0MjdnaGgoO3XqRIlEwsDAQO7Zs4e5ubl88OABpVJpvrFuFy9e\nzD59+pRoXEqlksePH+eAAQMokUjYo0cP7tmzp9TnLv7qq684bdq0Uj3Gh6R27dof3djSeVxcXHjt\n2jV1h1EkT58+pZ2dHfft21ek/SIjI1mrVq0i1VrPnj1LKyurdw7e4e7uzk6dOhGAanpFoeyIRFzO\nZWRk8JNPPuG2bdvyPc/pqaXFFi1acPXq1W+cRq969eqMiopSvff29ub27dtLLc7k5GQuXryY9erV\no6OjIydNmsTY2NhSOZarqyvPnz9fKmV/iKpVq1ascZYrgk6dOnHLli3qDqPIDhw4wMqVK6tGxSqM\nZs2aMbwQ40m/atiwYRw6dOhbt+nfvz+XLFlCU1NT6unplcidM6HwKqm7sZjwdj/++CNq166Nbt26\nIf6/Fs8A0KFjR/z111/o378/jI2NX9uvTZs22Lt3LwDg4cOHuHTp0msjbZUkiUSCoUOHIioqCr//\n/jseP34MT09PdOjQAb///jtyc3NL5Dj3799HcnIyPD09S6S8iiArK0vVn/xjU6NGDVy5ckXdYRTZ\nZ599Bn9/fwQHB4OF6Lhy/PhxPHz4EL169SrysaZMmYIdO3bgzJkzBW5Tv359REVF4dChQ8jKykKH\nDh2KfBzhPaj7SkAo2NmzZ1XTlq1bt45W5ubsb2rKuS1avLOBSmRkpKpRR2hoaL7ny2UlPT2da9as\nYdOmTWltbc3vvvuOt27deq8yly5dysDAwBKKsGKwtLQs1AxbFdG6devo5+en7jCKJSsri7Vr1+ay\nZcveuW2nTp24aNGiYh9r3bp1rFOnToGPjU6dOsXatWuTfFE71gY4rVGjIvWwEIpPJOJyKjs7mx4e\nHlyzZg2nTp1KBwcHXr58mVKplI8ePXrn/mlpaTQ0NGRqaiobNGhQ5OdRJe3atWscPXo0ZTIZmzdv\nzvXr16u6XBVF9+7duXbt2lKI8MNlYmJSLrqUqcO5c+dYs2ZNdYdRbFevXqVUKuX169cL3ObSpUu0\nsrJ6r8kjlEolmzdvznnz5r1xfWZmJvX19ZmRkUGFQsGe2tofXPewD5m4NV1OzZgxA7a2tjhx4gQ2\nbdqEkydPQiaTQaFQwNra+p37GxoawsvLC+Hh4Xjw4AF8fX3LIOqCVa9eHbNmzcLDhw8xbNgwrFmz\nBnZ2dggJCcGlS5cKVYZcLsfBgwffOHznxywzMxP6+vrqDkMtqlevjtu3b5fYo4+y5ubmhsmTJyMw\nMCqjh+YAACAASURBVBA5OTlv3GbmzJkICQl5r5+xhoYGFi1ahMmTJyMuLu619Xp6enBzc8OFCxdQ\nqVIltX9ffGxEIi6Hrly5gjlz5kChUODOnTs4cuQIbG1tceXKFdSoUaPQowe1adMGq1atQq9evcrN\nDDo6Ojro1asX9u7di6ioKEgkErRr1w5eXl5YsWIFnj9/XuC+Z86cgaOjY6EuRD4WCoUCcrkcOjo6\n6g5FLQwMDFC5cmXcvn1b3aEU29ChQ2FjY4MffvjhtXUPHjxAZGQkvvzyy/c+jpubG4KDgzF69Og3\nrs97TgwAPVasQB89PXQF0Oinn9772MLbiURczigUCvTr1w8mJiawtbVFZGSkakD2y5cvv3Fs2YK0\nadNGNSZzeeTk5ITJkyfj/v37mDhxIiIjI+Hg4IDBgwfj9OnTrzVi2bt3L9q0aaOmaMunvIZa5XVo\nx7LwoTbYyqOhoYGVK1dizZo1OHToUL51oaGhCA4OhqmpaYkca8KECThx4gT++uuv19Y1aNAAZ8+e\nBQBY2dvD9dtvEQngp8aNsTkwEI9jY0skBuF1IhGXE/ExMdgcGIiJtWrhxuXL6Nu3L1avXp2vpnPl\nyhXUrFmz0GVqaGhALpeX+xqklpaWqnX11atXUaVKFQQGBsLT0xPz5s3DtehobA4MxMOFC9Hwv7mX\nhReysrI+2tvSeT70RAwAlpaWWLlyJfr3749nz54BABISEhAWFoaRI0eW2HEMDAwwb948DBs2DNnZ\n2fnWvZyIAcDV1RW9DAyw7MkT9AoPx5GxY0ssDiE/kYjLiaPjxqFXeDimXL2K4U5O+Omnn16r5eTd\nmi6sTZs2oXr16jhw4EBJh1tqbGxsMHbsWNy8eRNz587FqVOnMLF+ffQKD8fKxERkbd2q7hDLlczM\nzI+261KeipCIAaBdu3bo0qULhg4dCpKYP38+evXqVeIzSXXu3FnVZuNlNWrUQGxsLFJTUwEAVapU\nKVTXKuH9iURcDt2/fx/jxo1DRkaG6jOSRUrEJLFx40b4+/ur+hN/SPIajISHh6Nzly6qz7f9/jt8\nfHzw448/4vDhw69d1X9sRI244iRi4EUjzcv//IOpjRrhxowZGBgYWCrHmTt3LmbPno27d++qPtPS\n0kKtWrVw/vx5AICZmRm2Zmaiv4kJ5vr6wnv69FKJRRCJuNzwnj4dc5s3x5fW1hh//DgePHgAd3d3\n7Ny5EwAQHx8PTU1NWFpaFqq8s2fPQlNTE4MHD8bff//9wbYqBYDWv/6KzQEBCJZK0X7BAowbNw6Z\nmf+PvTOPqzH74/in0mJr3/dFokiSdQwqImRJSdllLdsY69i3kcEYjN1YxhKR3TBki5CQNYZos0Si\nbGm7n98fpvvTtGi5t3vjvl8vL+5zznPO56p7v88557tkYMKECdDW1oarqyvmz5+PCxcuVOr3WRZk\nK+JPntMPHz4s0uu4MlG1alUEWlhg6uXL2JWVhaS1a8Uyj7m5OcaPH49Ro0blW/V+vj29fPlysEoV\n2E2diihDQ1mtYjEiM8RSgp6JCTQHDcJbV1fYOzlhx44d2LBhA8aNG4fu3bvj9OnTpdqWzlsN6+np\nwdLSEpGRkWJUL170TEzgvX07msyfj6NhYejQoQMWLlyIy5cv4/Hjx/jhhx/w6tUrjBw5ElpaWvna\nc3JyJC1frHzLoUt5qKiowNTUFA8ePJC0lDKTnZ2Nw4cPw9fXF6dPnaqQOceNG4e4uDjs379feC3P\nEF+/fh2hoaGoXbs2GjVqhMOHDxcb0SCjfMgMsRTx8uVLaGtrC1+3bdsWN2/ehIODA4YMGYKsrKwS\nrfhyc3Oxa9cuobe0m5tbpdye/i89e/bEiRMnkJqaKrympqaGzp07Y8mSJbh27Rri4+MxbNgwPHny\nBP7+/tDW1s7XnpubK8F3IHq+5fSWn1MZt6dJIiIiAgEBATA0NERQUBBatWqF2VeuYGOXLvBRUhJr\n6JCSkhJWrVqFMWPG4N27dwD+b4hHjRqFuXPnwsbGBqmpqWjdujVCZf4ZYkNmiKWI1NRUaGlp5bum\noqKCmTNnolOnTnj9+jUaNWqECxcuFDvO+fPnoaOjg7p16wL4FMZ0/PhxsemuKNTV1dGpUycEBwcX\n2UdTUxPdu3fH8uXLcevWLTx48AD9+/fHw4cP0bt3b+jo6KBbt25YtmwZbt26BYFAUIHvQPTIVsSf\nqEyG+O7du5g2bRqsrKwwZMgQGBkZ4fLlyzh//jxGjBgB2wYNMOjAAbxo0QKXo6PFqqVNmzZo3bo1\n5syZAwCwtrbG8+fP8ebNG/j7+8PKygqPHj1Cv3798Oeff4pVy7eMzBBLEf9dEX/O48ePsXr1akyd\nOhXe3t4YMmRIvpXh5wQHB+eLHW7RogXu3r1bZP/KxIABA7B58+YS99fR0YG3tzdWrVqFu3fv4s6d\nO/Dx8cHt27fRvXt36Onp5WuvbF6ishXxJ6TdED99+hRLliyBo6MjXF1d8fHjR4SGhuLOnTuYOnUq\nLCwsCtwzaNAgbNy4UezaFi9ejE2bNuH27dt4//49BAIBBgwYAAUFBVhZWeHhw4fo3Lkzbty4gYSE\nBLHr+RaRGWIpoihDnOcxXa9ePfj4+CAmJgZVq1aFnZ0dtmzZks94ZGdnIzQ0NJ8hVlZWRqtWrQoN\n4q9suLq6Ijk5ucRpMf+LgYEBfH19sX79esTGxuLatWvw8PBAVFQU3N3dYWhoCF9fX6xbtw4PHjyQ\nesMsWxF/wtbWVuoMcXp6OjZt2oS2bduiXr16uHPnDhYtWoSkpCQsXrwYDRs2LDYRS48ePXDx4kU8\nefJErDr19PQwfswYLHNxwc9Nm8LK1BTv378HAFhaWuLhw4dQVlZGz549sX37drFq+VaRGWIpIjU1\ntVBD/OTJE1StWlXYpqamhuXLl+Pw4cNYsWIF2rRpg5iYGABAWFgYatWqBXNz83xjfF4WsTKjoKCA\nfv36YcuWLSIZz8TEBP369cOmTZsQHx+PCxcuoF27djh37hycnZ1hYmKCvn37CtulDdmK+BM2NjaI\ni4uTuOd0VlYWDhw4gJ49e8LU1BQHDhzA8OHD8eTJE2zcuBGurq5QUFAo0VjVqlWDt7d3hWwJm8fE\nYH1KChbcvYs+GhpCz+m8FTEA4fa0tD+cVkZkhliKePnyZYEzYqDoRB5OTk6IjIxEz5490bp1a0yZ\nMgXbtm2Dr69vgb5ubm44fvz4V/EhGjBgALZt2yaWUCULCwsMGjQIW7duRVJSEk6fPo3vv/8ex48f\nR7NmzfK1P378WOTzlxbZivgTKioqMDMzw/379yt8boFAgHPnzmHYsGEwNDTEr7/+irZt2wo9kr28\nvMr8Mxo4cCA2bdok/s/tZytzfQMDoSE2NTVFcnIyMjMz0axZM+Tm5ubLviVDNMgMsRRR1NZ0cYk8\nFBQUEBgYiJs3byI2NhY7d+4sNC9t7dq1oaCggLt374pcd0VTu3ZtWFlZiX2FLycnB2trawwdOhTB\nwcF49uwZjhw5gkaNGuHAgQNwcHDI156cnCxWPYUhS+jxfyr6nPj27duYMmUKLCwsEBAQAAsLC1y9\nehVnz57F0KFDoampWe45mjZtiipVqiAiIkIEiovmib09/DU14auiAqtRo5CdnY0nT55AUVERxsbG\nSEhIgJycHPr27YutW7eKVcu3iMwQSwkkC/WaBkqW2jLv7NPe3h5z586Fp6cnkj5L0i4nJ/fVhDEB\npXfaEgVycnKwtbVFYGAg9uzZgxcvXiA0NBR2dnbYtWsXbG1tUbduXWH7y5cvxa5JltDj/1SEIX78\n+DEWLVoEBwcHuLu7QyAQ4NChQ7h16xYmT54MMzMzkc4nJycndqetp0+fYsHixRh75gwMAwJw4swZ\nNG7cWFiJ6fPt6T59+mDnzp0SPwL42pAZYinh7du3UFFRgbKycoG2klZdCg4ORkBAAG7duoUGDRqg\nYcOGWLJkiXAL92sJYwI+xRSHhYVJ1BNcXl4e9vb2GDNmDPbv34+UlBRs374dlpaW2LRpE6ysrPK1\n5yXzFyWyFfH/EZchTktLw4YNG+Ds7Ax7e3vcv38fv/32GxISErBw4ULY29uLfM7P6du3L/bt2yeM\n9RU1o0aNwvDhw1G/fn34+fkhODgYjRo1Em5B5zls5f27bt26OHr0qFi0fKvIDLGUUNT5MEnExMR8\n0RC/ffsWx48fR48ePYSxxxcvXsTff/8NJycnXLx4Ea6urjh//jw+fvworrdRYaipqX0xpriiUVBQ\ngKOjI3788UccOXIEqampWL9+PQwMDLB69WqYmZkJ2w8fPixMrl8eZCvi/yNKQ5yZmYm9e/eiR48e\nMDMzw9GjRzFq1Cg8ffoU69evR5s2bSAvXzFfn3p6emjVqhV2794t8rH37t2LmJgYTJs2DQDg6OgI\neXl5qKurF+qwBUAWUywOKEMqiIyMpJOTU4Hr8fHxNDAw+OL9W7duZadOnQpcFwgEDA4OpqGhIfv1\n6sWBmppc5uLC5MREkeiWJMePH6ejo6OkZZSYzMxMnj9/nnPnzqWLiwtr1KjBxo0bc9KkSTx27Bjf\nvn1b6jEDAwO5fPlyMaitfHz8+JHKysr8+PFjme7Pzc3l6dOnOXjwYGpqarJNmzbcsGEDX79+LWKl\npWf//v1s2bKlSMd8/fo1DQ0Nee7cuXzXZ86cySFDhlBTU5MCgYChoaH08PAQtqelpVFVVZWpqaki\n1fMtI1sRSwllcdT6nLzc0v9FTk4OvXr1QkxMDPSuXcPGV68w+tQpnBw/XiS6JYmLiwtevHiBmzdv\nSlpKiVBSUsJ3332HadOm4eTJk0hJScEvv/wCZWVlzJ8/H/r6+vnaMzIyvjimLHzp/ygrK8PCwgL/\n/PNPqe67efMmJk6cCDMzM4wdOxa1a9fGjRs3cPr0afj7+xfq/FjRdOzYEQ8ePBCpV/jEiRPRtWtX\ntGzZMt91X19fHDp0CFWrVkVcXFyBFbGamhrc3d2xa9cukWn51pEZYimhOEetevXqffHec+fOoUuX\nLkX2UVNTQ2MnJ+Hrffv2ISAgADdu3Ci7aAkj6pjiikZFRQVt2rTB7NmzER4ejhcvXmD27NkgienT\np0NHR0dY8jE8PLzQko+y8KX8lHR7OjExEUFBQahfvz48PDygoKCAo0eP4vr165gwYQKMjY0rQG3J\nUVRUFMazi4IzZ87g6NGjWLBgQYE2GxsbGBoawtzcHFFRUbC0tERcXFy+EKq+ffvKtqdFiMwQSwnl\nWRHv3bsX7du3R82aNYvt1yooCLv9/PCjpSXUfX2hr6+Pzp07o3nz5ti8eXO++seVhf79+2P79u1f\nRfnDatWqoW3btsKSjsnJycKSj+PHj4e2tna+9qRHj6AZHo7na9fi+Wce8t8yxRniV69eYd26dWjd\nujUcHR0RHx+PVatWIS4uDgsWLPjiA6+kGThwIP78889yVxTLyMjAkCFDsHLlSqipqRXax9fXFx8/\nfsSVK1dQs2ZN1KxZE8+ePRO229etC9PoaPzh4SH73RMBVSQtQMYnijPEw4YNK/be4OBgjBw58otz\n5JUTbJaUBAcHB9y6dQs//fQTjh49irVr1+LHH39E7969MWzYsFKVXJQktWvXRq1atXDs2DF4eHhI\nWo5IqVGjBjp06IAOHToA+JQy8dy5czh16hQCAwNhdfs29uTkAI8fY4KLC7I7dYKSkhKUlZWL/Lu4\ntuL+Li4VozRha2uLnTt3Cl9nZGTg8OHD2L59O06fPo327dtj3Lhx6NChQ6ERCtKMra0tTExMcPz4\ncXTs2LHM48yZMweOjo7F7qD16tULc+bMQbVq1QD833Pa0NAQAHBp+nTsyswEDh/GblVVeMtSX5YL\nmSGWElJTU+Hg4JDvmkAgQExMDGxtbYu879mzZ4iOji7VB9PExASDBg3CnDlzsGbNGnh4eMDDwwMJ\nCQnYsGED2rVrBysrKwwbNgxeXl5SfwaZF1P8tRni/5JX8rFz584AgG3e3sCePQA+VaaqZmaGrKws\nZGZmIjMzE2/evBG+Ls/fWVlZUFRULJMBL4/xL+nfn6eM1NfURI2wMKxwdcUlLS0cDQuDo6Mjevfu\njS1bthS5Aqws5MUUl9UQX79+HX/88ccXc7UbGxujfv36iIqKQm5urvCc+Pvvvy/TvDKKR478CnIe\nfgV4eXnBx8cH3t7ewmuPHj1C69at8yXm+C/Lli3DtWvXSn1O+urVK9jY2CAiIgK1a9fO15ZXpHzN\nmjW4du0a+vbti2HDhsHGxqZ0b6qCSE9Ph5mZGWJjY4usXvU18jwpCYcCAnAyLAyLY2JgVEgFH1FA\nEtnZ2fkMdHmNe3n+/u81eXl5oWF2y8j4tFIDMMfREYMPHRKu4r4GyvO7npOTg2bNmiEwMBADBw78\nYv+8XbLLly9j165dEAgEmPtvfeTnSUn4tV07JCQkYNn9+9AzMSnT+5HxLxL12ZYhpHXr1jx16lS+\nawcPHmT79u2Lva9Zs2Y8evRomeb8+eef6e3tXWyfhw8fcvLkydTT02ObNm0YHBxc5vAQcdK7d+9v\nNoynSZMmZf4dqOwIBAJmZWXx3bt3TE1N5RZPTxIgAYb4+Ulanljo06cPf/vtt1Lft2jRIrq6ulIg\nEJSo/8uXL6moqMg1a9Zwy5Yt9PvP/2fnzp1pZ2dXah0yCiIzxFKCnZ0db968me/aggULOG7cuCLv\nefToEbW1tZmVlVWmOd+9e0cDAwNGRUV9sW9mZiZDQkLo4uJCXV1dTpgwgQ8ePCjTvOLgxIkTlSqm\nWJSsXLmSvXr1krQMqSA5MZELW7TgCEPDryJWvjBOnTrF+vXrl9igkmRsbCy1tLQYGxtb4nueJSSw\nr6oqh+nr81BoKJs2bZqv3c7Ojs7OziUeT0bRyLympYTCwpeK8phOTkzE7t69sdnTEx7t20NRUbFM\nc1avXh0zZszA5MmTv9hXSUkJ3t7eOHnyJM6fPw+SaNGiBdq1a4c9e/ZI3GvZ2dkZKSkplSamWJT4\n+Pjg6NGjSE9Pl7QUiaNnYoIRx45h+7t3UKxeXdJyxELr1q3x7t07XLt2rUT9SWLYsGGYPHkyrKys\nSjzPuSlT8OebN1iTnIxX27fniyUGPuXdNjIyKpV2GYUjM8QSJjkxESG9e6PZ8+fI+U8Ch6IM8bkp\nU+C9YwdmX78Ox+fPyzW/v78/EhMTERYWVuJ7rK2thQXOBw0ahBUrVsDU1BQ//fQT4uLiyqWnrOTF\nFFd0IQhpQEtLC66urmJJgVgZqVmzJtq1a4d9+/ZJWopYkJeXx8CBA0tcCGLLli14/fo1xo4dW+Y5\n5eTk8OHDB2Fa1rS0NGRkZMDAwKDMY8r4PzJDLGHOTZmCnjt2YB+Ji9OnC6/n5ubi3r17hXpMCz7z\nr9Mpp3OSoqIi5s2bh8mTJ0MgEJTqXmVlZfj6+uLs2bM4deoUMjIy0LhxY7i7u2P//v3ljncsLV9T\nTHFpkeX/zU+vXr3yhTF9bfTv3x+7du36Yt7458+fY9KkSfjjjz9QpUrpgmRaBQUhwNAQgzQ18bJJ\nk3zFH2JiYqClpQVVVdUyvwcZ/0dmiKWIAwcOIDw8HMAnj2ldXd0CSTpSUlKwLikJQ3V10V9VFdW8\nvMo9r9e/Y+z5NxSmLNStWxdLly5FUlISfH19sWjRIpiZmWHmzJnFen2LEmtra1hbW3+TlWHc3d1x\n7949PHr0SNJSpIKOHTsiKioKz8u5YyStmJqawtHREfv37y+23+jRozFo0KACoZEl4X12NnZnZcF1\n+XIcDQvLl+ryzp07UFNT+2ISIRklQ2aIJUxetqtAY2McIeHn54c+ffrg3LlzBbalL1y4gEaNGqFJ\ny5ZY9eQJzMeOxbnIyHJrkJeXR1BQEKZOnVru1WTVqlXRr18/RERE4NixY3j16hUcHBzg4eGBI0eO\nIDc3t9x6i2PgwIHf5Pa0kpISevXqJSva/i/VqlVD586dERoaKmkpYuNLdYoPHjyIa9euYcaMGWUa\nf82aNRgwYAC6d++Oy5cvw8DAQPigFxMTgxo1asgMsaiQtLeYjE/ExcVRSUmJY8eO5ZQpU1i1alW2\nadOGWVlZFAgEXLp0KXV0dHjw4EHhPVeuXKGNjY3INLi6unL16tUiGy+Pd+/e8Y8//mCTJk1oamrK\nOXPm8MmTJyKfhyTT09OppqbGFy9eiGV8aebKlSu0tLQslTft18zBgwf5/fffS1qG2MjIyKCWlhYT\nEhIKtKWlpdHY2LhASGRJ+fDhA7W1tYVe1r169WKvXr04dOhQkqSbmxtbtGjBkJCQsr8BGUJkK2Ip\nwdzcHC4uLvjjjz/g6+sLV1dXvHz5EvXq1UObNm2wdetWREZG5sse5ejoiHfv3pW62kxRBAUFYc6c\nOXj//r1IxsujevXqGDRoECIjI3HgwAE8ffoU9erVQ/fu3fH333+X+my6OFRVVeHh4SFVdYorCkdH\nR6ioqCAiIkLSUqQCNzc33LlzB48fP5a0FLGgoqICHx+fQpP5TJkyBR06dICzs3OZxg4JCYGTk5PQ\ny7pnz564e/duvq1pBQUF2YpYRMgMsRQxYcIEqKqo4DdnZ2idP4+Rw4YJwxQsLCwKhCnJycmhS5cu\nOHjwoEjmd3JyQsuWLbF8+XKRjFcYDg4OWL16NRISEuDu7o4pU6agVq1aWLBggcjO8/JSXn5ryMnJ\noX///pW2GpWoUVZWRrdu3b5qb/JBgwZh06ZN+R5mz58/jwMHDuCXX34p87irVq1CQECA8LW7uzse\nPXqEf/75B+np6UhLS0N2drbMEIsKSS/JZfwfgUDAfqqqwsxA3kpK/PPPP/n+/XvOmDGDWlpaXLBg\nQb7MVkePHuV3330nMg3//PMPtbW1K6zot0Ag4OXLl+nv7091dXV6eXkxLCyMubm5ZR4zNzeXpqam\nvH79ugiVVg4eP35MDQ0NfvjwQdJSpIK///6bTZo0kbQMsSEQCGhvby/cgs7IyGCdOnUYGhpa5jGj\noqJoZmbGnJycfNd9fX2poKDAs2fP0snJifXq1eONGzfKpV/GJ2QrYilCTk4ONnXqCF87Ozujb9++\nqFatGmbPno3IyEhERETA3t4ef//9t7DP7du3kZKSIhINtWvXRo8ePQqtUyoO5OTk0LhxY2zYsAHx\n8fFwdnbGDz/8ABsbGyxatKhM70teXr5S1ykuD0ZGRmjcuDEOHDggaSlSgYuLC+Li4r5ab3I5Obl8\nTlvz589H3bp14enpWeYxV61aheHDh+crpgF8Ko1YpUoVnD9/Hra2tnjz5o0sfElUSPpJQEZ+Ht67\nR28lJQ7R1eWKJUsK7XPo0CFaWlqye/fujI+Pp5eXFzdt2iQyDU+ePKGmpiYTJZQiUCAQ8MKFC+zf\nvz/V1NTYq1cvnjlzplROSA8ePKCurm6Z039WZrZt20Z3d3dJy5AaRowYwQULFkhahti4HR1NL0VF\nrnRzo66GRrkcIVNTU4t0dox/8IDd5OQ43MCAM6ZMoYaGBl++fFke6TL+RbYiljIsbWxgPmYMXjRv\njh1FxPV27twZd+7cQcOGDdGoUSOQxN69e0WmwdDQEEOHDsWsWbNENmZpkJOTQ/PmzbF582bExcWh\nefPmGDFiBGxtbbF06VK8evXqi2PUqlULtWvX/iZjirt3746LFy8iOTlZ0lKkAh8fn686uUfMokXY\nnZ2NgOPHMcbaulzVpvLKiero6BRouzxzJvaRWP3sGdQvXMDbt29lZ8QiQmaIpZCAgACEh4cjPj6+\nyLqhKioqmD59Oq5cuYKMjAwcOXJEpMZ40qRJOHjwIGJiYkQ2ZlnQ0NDA6NGjcefOHaxbtw5Xr16F\npaWlMFaZxVTxHDhwIDZt2lSBaqWDatWqoXv37tguK9YOAGjZsiVSUlJw9+5dSUsROenp6fjn/n3h\na6tatco8lkAgKOCk9Tl56S2BT2lE88pPyig/MkMshZibm6N169Zo2LAh1q9f/8W+R44cgZ2dHUaN\nGgUPDw+RnIepq6tj4sSJmDZtWrnHEgVycnL4/vvvsW3bNjx8+BAODg4YNGgQ7O3tsWLFCqSlpRW4\nx9vbG6dPnxbZ+XllQpby8v8oKCigZ8+e2LVrl6SliIzo6GgMHToU5ubmuKqnh2Vt2qCHggIaTZ1a\n5jFPnDiBmjVrolmzZoW2H69SBT7KyughL49ms2fLVsMiRGaIpZTRo0fj/v372L59Oz58+PDF/gMG\nDECHDh3w3XffoUmTJpg5cyYy/lNEorSMHDkSUVFRuHTpUrnGETVaWloYN24c7t27h+XLlyMiIgLm\n5ubCWOW8VXLNmjXRpUsX7NixQ8KKK55WrVohPT0d169fl7QUqcDHxwe7du0qdgdF2snIyMCWLVvQ\nrFkzdO3aFaampoiJicG+w4cx5vRpvHF2xp0HD8o8ft5qWE5OrkDbP//8g1Pnz+MvRUXE1quHaurq\nMkMsSiR7RC2jKAQCAevVq8cmTZpw8+bNX+wfGxtLfX195ubmMjExkT179qS5uTn37dtXrkxLGzZs\nYOvWraU+W9Pz588ZFBRES0tLNmjQgKtWrWJ6ejpPnjxJBwcHScuTCNOmTeMPP/wgaRlSgUAgoJmZ\nWaUMabt//z5//PFHamtrs0OHDjxw4ACzs7ML9Fu0aBFHjBhRpjni4+OpqanJd+/eFdru5+fHadOm\nUUlJiX369OH169dZv379Ms0loyAyQyzFrFu3jo0bN2aLFi1K1N/W1paXLl0Svg4LC2PdunXZoUMH\n/vPPP2XSkJ2dzTp16vCvv/4q0/0VTW5uLo8fP84ePXpQXV2dfX186KOiwjXu7l9tofiiuH//PvX0\n9L5Jz/HCmDRpEidPnixpGSUiOzube/fuZbt27aijo8OJEyfy4cOHxd5z8+ZNWlpalmm+KVOmcMyY\nMYW2xcTEUEdHh5cuXaKmpiYXLFjA8PDwEn8vyfgyMkMsxbx//56amprU09PjrVu3vth/8uTJh7xB\n+AAAIABJREFU/Omnn/Jdy8zM5KJFi6ilpcUpU6YU+cRbHHv37mWDBg3KlWRDEjx9+pQz7O2FCVJ+\nrFWL27dvZ1JSkqSlVRjNmzfn4cOHJS1DKoiOjqaFhYVU7+48fvyYs2bNopGREVu0aMFt27YxIyOj\nRPcKBAIaGhrywYMHpZrz48eP1NXV5b179wpt9/X15c8//8y//vqLWlpaPHjwII8cOcL27duXah4Z\nRSM7I5ZiqlWrBn9/f5ibm2PdunVf7N+1a9cC6S6VlJQwfvx43Lx5EwkJCahbty52795dqrOybt26\nQUVFpdLlbzYwMEC9evWEr6tVq4bQ0FA0bNgQVlZWwkpNjx49qtRnh8UhS3n5fxo0aABFRUVERUVJ\nWko+BAIBwsLC0KNHD9SrVw/Jycn466+/EBERgd69e0NFRaVE48jJycHNzU2Y7Kek7NmzB/b29rCx\nsSnQFhMTg7CwMIwcORJJSUn4+PEjdNTUcOfnn2EdE4PnFVTi9KtH0k8CMoonLi6Oampq1NDQ4Pv3\n74vtm5ubSz09vWK3sM6ePcv69evT1dWVMTExJdZx+vRpWlhYMDMzs8T3SAPJiYkM8fPjeGtrjgkI\nIPlp5XDnzh2uWrWKvXr1ooGBAY2Njenn58e1a9fy7t27Ur1qKg2vXr2impoaX716JWkpUsGMGTOk\n5tw8NTWVS5YsobW1NevXr89Vq1bxzZs35RozODiYHh4epbqnRYsW3Lt3b6FtPj4+wmQo48ePp6Ki\nInf6+gp3mUL8/MqlV8YnZIa4EtC9e3fa2dmVyGnL39+fS5cuLbZPdnY2f/vtN2pra3P8+PEl/vB3\n6NCBK1asKFFfaSPPGeX169cF2gQCAR88eMANGzawX79+NDMzo66uLr28vLh8+XLeuHGj0m3Lf463\ntzfXrFkjaRlSwZ07d2hkZCSxn6dAIOClS5eEWeN69+7N8+fPi+zB7+XLl1RVVS3xA3N0dDSNjY0L\ndf66ffs2dXV1+fbtW5Kku7s7TU1NPxlfmSEWKTJDXAk4ffo0jY2NS+QccfDgQTo7O5do3OTkZPbv\n359GRkbcsWPHF78MoqOjqa+vL/xgVjb69+/PuXPnlqhvfHw8t27dysGDB7N27drU1NRkly5duHjx\nYkZFRRX6xSWtHDp0SOZY8xn169dneHh4hc757t07rlu3jg0bNqSlpSUXLlwotprZjRs35unTp0vU\nd8iQIUV+Jry9vblw4ULh6zp16rBt27bcvG4duwL808vrm3OAFBcyQ1wJyAtl0tLS+qLT1vv371mz\nZs1SVU+KiIigg4MDW7du/cXx/fz8OHv27BKPLU3keX+WxWHt6dOn3LVrFwMCAlivXj2qqqqyQ4cO\n/PnnnxkRESHVW/ZZWVnU1dXl/fv3JS1FKpg/fz4D/j2mEDd37tzhyJEjhQ9yR48eFftqfNq0aSXy\nDn/9+jXV1dX57NmzAm23bt2inp5evs+Kmpoax40bR3Nzc5mjloiRGeJKwrp162htbc1Ro0Z9sW+X\nLl24bdu2Uo2fk5PDlStXUkdHh2PGjGFaWlqh/R4+fEhNTU2xPc2Lmx49enxx674kpKSkcN++fRw7\ndiwdHR1Zo0YNuri4cNasWTx9+rTUlSEcO3Ysp0+fLmkZUkFsbCx1dXXFtquRmZnJ4OBgtmrVivr6\n+pw2bRoTEhLEMldhnDt3jg0bNvxiv2XLltHHx6fQNi8vLy5atEj4Ojc3l/Ly8pwxYwbl5OQq9P18\nC8gMsYR5lpDAED8/hvj5FbvN8/79e2poaFBdXf2LTlt//PEHvb29y6TnxYsXHDx4MA0MDLhly5ZC\nn94DAwOLjDmUdq5cuUIjI6N8NZ1FQVpaGo8cOcKJEyeyWbNmrF69Olu2bMmffvqJx44dK7cTTnmJ\njo6mmZlZpT7rFiVOTk48ceKESMeMi4vjlClTqKenR2dnZ4aEhEgkhjsrK4vq6upMTk4uso9AIKCN\njQ3Pnj1boO3GjRsFVsPPnj2jvLw8raysZMccYkBmiCVMaRwfJk6cSHNzc27ZsqXYfsnJyVRTUyuX\nsYmMjKSTkxNbtGjB6OjoAuNramoyLi6uzONLkg4dOnD9+vVinePdu3c8ceIEp02bxlatWrF69eps\n3Lgxx48fz4MHD1a4F7NAIGD9+vV55syZCp1XWlm8eDH9/f3LPU5OTg4PHz7MTp06UVNTk2PGjOHd\nu3dFoLB8dO/enVu3bi2yPSwsjPXq1SvUL8TT05OLFy/Od+306dOUk5MjgFJFW8goGTJDLGE+N8TL\nXVyK7RsfH88aNWqwadOmXxy3efPm/Pvvv8ulLScnh+vWraOuri4DAwPzGY/p06ezX79+5RpfUoSH\nh9PKyqpCHa4yMjJ49uxZzp07l23btmWNGjXYoEEDjh49mnv27KmQrf7Fixdz4MCBYp+nMpCYmEhN\nTc0yn+0/f/6cCxYsoLm5OZ2cnPjHH398caeqIlm9ejX79OlTZLunpydXrVpV4HqeQ+Z/30tQUBDl\n5eW/2XSx4kZmiCVMcmIid/r6spu8PGsoK/Pp06fF9u/WrRvV1NS+6FQVFBQkMoeUly9fcvjw4dTT\n0+OGDRuYm5vL9PR06urq8ubNmyKZo6Jp2bIld+zYIbH5s7KyePHiRS5cuJAdO3akmpoa69aty+HD\nh3PHjh3lKu5eFE+fPi3R0ca3wnfffVeqrGMCgYDh4eH09fWluro6Bw0axKioKDEqLDuPHj2irq5u\noUcRSUlJ1NDQKPS4pHv37vz111/zXXuWkMARhobsAvDYoUNi0/wtIzPEUkL//v2poqJCGxubYh19\nTp8+TS0trS86bcXExNDY2FikiSmuXr3KZs2asUmTJoyKiuLSpUvZuXNnkY1fkfz111+sV6+e1JyZ\n5uTk8OrVq1y6dCm7detGLS0tWllZcdCgQdy8eTPj4uJE8rN0d3fn9u3bRaC48rNixYpiV415pKen\n8/fff6ednR1tbGz422+/VYoEKdbW1rx27VqB69OnT2dgYGCB69HR0TQwMCjw/SOLGxY/MkMsJWzd\nuJH9VFXZo0oVduvUqcgv3TwnC1VV1WINtkAgYK1atQr9IJaH3Nxcbty4kXp6evT396eJvj5/b9fu\ni85m0oZAIKCDgwMPHjwoaSmFkpuby1u3bnHlypX08fGhvr4+TUxM2Lt3b65bt4737t0rk2EODg6m\nm5ubGBRXPp49e0Z1dfUiP0fR0dEcOnQo1dXV6e3tzVOnTlWqjGsjR44UZsXKIzMzk/r6+rx9+3aB\n/l27di00ouC3Nm1khljMyAyxlLDDx0f4y95LRaXYxBPr16+njo7OF522fvzxR86cOVPESj/x+vVr\njho1ij2qVKm0H9KQkBA2bdq0Uny5CgQC3r9/n+vXr2ffvn1pampKPT09ent78/fff+fNmzdLtLr/\n8OEDNTQ0+Pjx4wpQLf24uroyNDRU+DojI4Nbtmxhs2bNaGJiwrlz537xuEhaOXToENu0aZPv2q5d\nu9i6desCfa9evUpDQ8MCDyVpaWmsVqUKu8nJVbqH7cqEzBBLCZ9v/3grKVFfX5979uwptG9e0o5G\njRoVO+bZs2dLFE9YHtZ06CDUPdvRUeJhOqUhJyeHNjY2PHXqlKSllIn4+Hhu2bKF/v7+rFWrFjU1\nNdm1a1f++uuvvHLlCnNycgq9b/DgwfkyJn3LLFmwgKNMTbnBw4OBQ4Z8seZvZeLt27esUaNGvkx4\nrVq14q5duwr09fDw4LJly/JdEwgEdHd3J4BKm8SnsiBHfqVlZyoZz5OScHDECBw5cgRyHTsi7skT\nPH78GMePH4ejo2OB/hMmTMDq1asRGRkJOzu7QsfMycmBvr4+rl27BlNTU7Hp3t6nD6KvXcO71q1x\n7tIlDB06FKNHj4a+vr5Y5hQlmzdvxrZt2xAWFiZpKeXm6dOnCA8PR3h4OM6ePYsnT56gRYsWaN26\nNVq1agUnJycoKiriwJ492Ovvj06dO6N1UBD0TEwkLV1ITk4OsrKyCvzJzMws9HpxbSW5rnX+PFY/\newYAmGJriyGHDsHS0lLC/wuiw8XFBT/88AM8PDxw+/ZtuLm5ISEhAYqKisI+V69eRdeuXREbG5uv\n0lNQUBDmzZuHzMxMZGdnS0L+N4PMEEsZM2fOxPz582Fvb48mTZrgyJEjiIyMhKGhYb5+CQkJqFOn\nDgYOHIhVq1YVOV7//v3RpEkTBAYGik0zSTRr1gw//PADmjRpgl9//RU7duxAjx498OOPP6JOnTpi\nm7u8ZGdno1atWggJCUHTpk0lLUekpKSk4Pz58zh79izCw8Px4MEDNG3aFM1fvMDcW7cAAH907oym\nCxaI3eCV1KiShLKyMpSVlaGkpFTgT2mvf+mexytXYlx4OABgt58fvLdvl+SPTOQsXLgQSUlJ+P33\n3xEYGAhtbW3Mnj07Xx8PDw+0b98eI0eOFF4LCwtD3759kZycjLFjx2Lp0qUVLf2bQmaIpRBbW1uk\npaUhKysLgwcPxqlTp3D27FlUrVo1X78OHTrg/PnzSElJKdCWR2hoKNatW1fqGqWl5eTJkxg+fDhi\nYmKgqKiIly9fYtWqVVi5ciWaN2+OCRMm4LvvvhOrhrLy+++/IywsDPv375e0FLGSlpaGiIgInJ8w\nAQvu3gUA9K1RA9FmZmI3eCW9rqCgUKH/J8+TkrCtd2/8888/mHvlilTtDoiC69evw9vbG1evXoW5\nuTlu3boFIyMjYXtUVBS6d++ebzWckJCApk2bwsDAALdu3UJWVhbk5WWl68WK5HbFZRTFs2fPqKio\nyIYNG9LHx4d+fn708fEp4FR05swZVq9evdjyiG/fvmXNmjWLzB0tSlxcXLhu3bp8196/f89Vq1bR\n0tKSzZs35969e6UmZCiPDx8+UF9f/4ux2V8LyYmJ3NytG72UlPhQCrJASZrMzExqa2sXW8e7spJX\no3z27Nns0aNHgfaOHTvy999/F77OyMhgw3r1OMPBgV0B+nl5VaTcbxaZIZZS1q1bR3l5eerr63P/\n/v1s2rRpAYcJgUBAMzMz2tnZFTuWu7t7oQ4aoubSpUs0NjZmRkZGgbacnBzu3r2bjRs3Zu3atbl2\n7dpC+0mKoKAg9u7dW9IyKpRu3bpx7dq1kpYhFYwaNUpsEQaS5FlCAseam7OXigr3BAfna8v7vH6e\nCtff35+jTE2FDpg7fX0rWvI3icwQSzEtmzalp4ICfatVY+S5c5+KcoeE5OuzZs0aKisrFxoXmMfq\n1asrzMjkee0WhUAg4JkzZ9ipUyfq6+tz3rx5pSrZKC7S09OppaXF2NhYSUupMI4dO0YHB4dKEb4l\nbq5evUpzc3Op260pLzt9fYVGddd/wgvd3d3zpblct24dbW1tubFbt0obklhZkRliKWabl5fwAzG5\nTh1GR0dTW1s7X1q9R//8Q095eY61sCgyxu/x48fU0NCokEowt27doq6ubonCmG7fvs0BAwZQQ0OD\nY8aMYXx8vNj1Fce0adM4dOhQiWqoSHJzc2llZcWLFy9KWorEySuKcfr0aUlLERmRkZEcoK5eqFG9\nePEiTUxMhKvhyMhI6ujo8N69e7SrXZvd5OS409dXFjdcQcgMsRTzeWyxp7w8o6OjuXfvXhobGwsT\nMnzeZ0cRtUXJT2XfKipetk+fPqWKO3z8+DEnTJhATU1N+vn5Faj2VFGkpKR8c8kufvnll0pbvEPU\nLFmyhAMGDJC0jHKTlpbGkSNHUl9fn6uXL+dkW1vOsLfPZ1Tbt2/P1atXk/xU+tTU1JT79u3j7t27\nKScnx3HjxklK/jeJzBBLMcmJiQzx82NQs2ZUBGhpacmcnBz+/PPPbNSoEd+/f5/PEA/S0ipym3fO\nnDkVVkM4NjaWWlpafPnyZanuS0tL4y+//EIjIyO2bduWx48fr/Bt07Fjx35TX0IpKSlUU1OTiuMB\nSZOcnEx1dfV8CTAqEwKBgLt376aRkRGHDBki/Jl6eXkx+LPz4QsXLtDU1JSZmZnMzs6mi4sLp0yZ\nwmfPnrF69epUUlKq9MlMKhsyQ1xJGDlyJOXk5DhmzBgKBAL26dOH3t7efBofzxA/P/pVr04lgDo6\nOly5cmWBD9KNGzdoYWFRYYZtxIgRHD9+fJnuzczM5KZNm2hnZ0cHBwdu3769wgqs51WmKe1DRGWm\nT58+XLJkiaRlSAUeHh7FRiFIK3FxcezUqRNtbW157ty5fG2Ojo6MjIwUvm7Xrp3QSW/ixIls164d\ns7Oz2bFjRyopKX2xoIwM0SMzxJUEgUDABg0aEACPHz/OjIwMNmvWTOjpGRUVRQD09fVlmzZtWK9e\nPZ48eTLf/WZmZhUWovPkyRNqamqWa5s3NzeXhw8fZuvWrWlqasrffvutQlYrQ4YM4fTp08U+j7Rw\n/vx5Wltbf3WOSmVhz549BfIzSzNZWVn85ZdfqKWlxZ9//rnQ+srq6upMSUkh+elnbW5uzszMTO7Z\ns4dmZmZMSUnhunXraGRkREVFxTLXaJZRdmSGuBKRmprKqlWrUlFRkXfv3mVycjLNzMy4c+dOkqSZ\nmRmrVKnC69evMzQ0lObm5vT09BTGR44aNYrz5s2rML0TJkzg8OHDRTJWZGQkvby8qK2tzalTpzI5\nOVkk4xbGgwcPqKWlxfT0dLHNIU3kOSqdOHFC0lIkzsePH6mlpcVHjx5JWsoXuXjxIu3t7enm5lak\nt39qaipVVVWFO2Ft27bl+vXrGRMTI3T8zPt9V1VVpb+/f0W+BRn/IjPElYwTJ05QTk6OGhoafPTo\nEa9fv04dHR1evnyZ61etYleA/VRVmRAbyw8fPnDevHnU1NTkTz/9xEOHDrFJkyYVpvXly5ciDwmK\njY1lQEAANTQ0OGTIEN67d09kY3+Or6/vN1UYYdWqVfT09JS0DKlg5MiRnDVrlqRlFElaWhoDAgKo\nr6/PHTt2FHvcFBUVRQcHB5LkuXPnaG5uztTUVNapU4cbNmxgdnY2mzdvTm9vb1apUoXv37+vqLch\n4zNkhrgSMmDAAFapUoVmZmZMSkri/v37aWRkxM2fxf+NMTcXbuM+fvyYffv2paGhIatVq1ahXsFz\n5swRSwzzixcvOHPmTOro6LBbt26MiIgQ6fg3btygvr5+sTWfvybevHlDDQ0NPnnyRNJSJM6VK1do\nYWEhdVv1AoGAISEhNDQ05NChQ/nq1asv3rNz505hRi0XFxeuX7+ePXr04JAhQ0iSP40fz6F6euyp\nrMweXbqIVb+MopElEK2EbNiwAVpaWiCJtm3bonnz5hg1ahTOnD0r7JMQHw9tbW3Y2tpi8uTJcHR0\nxJiAAHTMycFEGxscPXiwQrSOHTsWYWFhuPVvkQFRoaOjg1mzZiE+Ph5t27ZFnz590LJlSxw4cAAC\ngaDc49vb26Nx48bYtGmTCNRKPzVr1oSPjw82bNggaSkSx9HREdWrV8f58+clLUVIXFwcOnXqhDlz\n5iAkJARr166FhobGF+979OgRLC0tER4ejvj4eKSkpCAxMRErVqzAtWvXELtiBdY+f45dmZnorqxc\nAe9ERqFI+klARtm4ceMG5eXl6eHhQXt7e758+ZK9PD050sSEv7ZqRa2aNWlsbMyIiAhu2LCBI0aM\n4CBNTeGKuZeKCvv161chK6ClS5eyi5iftrOzs7lr1y42atSINjY2XL9+fblTaF68eJFmZmYV5rEt\naa5fv05jY2NZ6ArJxYsXc+DAgZKWwaysLC5cuJBaWloMCgoq9e/i4MGDuXr1ajo7O3PChAnU19dn\nYmIiP3z4wLp163Ji7dqyLFpSgMwQV2LGjBlDRUVFBgQEsHHjxnzx4gVbtGjB6dOnc8uWLaxevTr7\n9u0r7P95zPGCZs04efJkobelOPM+Z2Rk0MTEpEIyOAkEAp46dYru7u7U19fnzz//XKItvKJwcXGp\nlOEsZaV58+bct2+fpGVInGfPnkk8pvjChQusX78+27dvX+aCFM7Ozly8eDHNzMyop6cnjKQYM2YM\nGzduzKoKChxrYcEQPz9ZFi0JIjPElRiBQEAjIyPa29tz2LBh/P777xkXF0czMzPu2LGDc+fOZZUq\nVYQFH/IShPzerh1N9PT45s0bxsbGslu3brSwsODevXvFFme8fv16uri4iGXsorh58yb79etHDQ0N\njh07lgkJCaUeIywsjDY2NszJyRGDQunjzz//pJubm6RlSAWdOnXili1bKnze169fc/jw4TQwMODO\nnTvL9Zk0MzOjk5MTLSwshM6HJ06coKqqKqtWrcoOHTrIco1LATJDXMm5desWFRQUOH/+fPbr149u\nbm6MioqitrY2L126xJ49e1JRUZGJ/3naHTBgQL4MUidOnKCdnR1dXV3FEmucnZ1Na2triYTIJCYm\n8scff6SmpiZ79+7N69evl/hegUDAJk2acM+ePWJUKD1kZGRQW1ubDx48kLQUibN79246OztX2HwC\ngYA7d+6kgYEBhw8fztevX5drvMzMTCoqKrJmzZrs3r07BQIBU1JSWKNGDaqpqbFhw4ZSVQHtW0Zm\niL8CRo8eTUVFRd64cYNeXl7s0qULQ0NDaWhoyISEBNatW5cGBgb5zpeeP39OHR0d3rx5U3gtOzub\nv//+O3V0dBgYGCjy7FI7d+5k48aNJfYE/vr1awYFBdHQ0JBubm4MCwsrkZb9+/fT0dHxm1k5jB8/\nnhMmTJC0DImTF1NcEcVIHj58yPbt27N+/fq8cOFCucd7lpDADZ07s5ucHI10dPjmzRvhEZGGhgZN\nTEzEGosvo3TIDPFXQGZmJg0MDGhmZsb09HR26tSJvXr1YlBQEB0cHIQ5ZFu2bJkvJGPVqlVs2bJl\nAQPz8uVLjhw5kjo6OlyxYoXInHdyc3Pp4OAg8TPIjx8/cuPGjaxbty4bNmzI4ODgYt9jbm4u7ezs\nePTo0QpUKTkePHhAHR0d2WqJZGBgIOfMmSO28bOysrhgwQJqaWlx4cKFInMM/NwfZF3Hjnz16hXr\n1KnDqlWrUktLq8Iy7MkoGTJD/JVw7tw5qqiocPDgwczIyKCLiwsHDhzIfv360dPTk2fOnKGioiKH\nDBkiNLw5OTl0cnIq0hnp1q1bdHFxoZ2dnci2lI8cOUJbW1upOHPNzc3loUOH+P3339Pc3JzLli3j\nu3fvCu27fft2fv/99xWsUHK4ublx69atkpYhcaKiomhpaSmW3ZCIiAjWq1ePHTp0EGkmr4yMDAYY\nGwsN8aZu3Vi7dm2qqKhQU1Pzm3mgrEzIDPFXhK+XF3tUqcIlrVrx4d27bNGiBQMCAtiiRQtOnTqV\nEydOZI0aNbhgwQLhPVFRUdTX1y/Ss1ggEHDfvn20tLRk165dy50lSyAQ8LvvvuOff/5ZrnFEzcWL\nF+np6UltbW1OmzaNz58/z9eenZ1NS0tLhoeHS0hhxbJ37162aNFC0jIkjkAgoJ2dnUh/7q9eveKw\nYcNoaGjIXbt2idTIx8fH09HRkSry8hysq8u1HTvSVF+fVlZW1NbW5u+//y6yuWSIDpkh/orY5uUl\nfAre0Lkz09LS6OjoyNGjR9PMzIybN2+mvb09tbS0uGHDBuF9w4cP54gRI4odOyMjQ7iFNmnSJL55\n86bMOsPDw2lhYSGVyeXv37/PYcOGUV1dncOGDeP9+/eFbWvXrqW7u7sE1VUc2dnZNDIy4o0bNyQt\nReIsWrSIgwYNKvc4AoGAO3bsoIGBAQMCApiWliYCdf/nxIkT1NPTo4GBAXV0dHj8+HHq6Oiwf//+\nVFVVZWBgoEjnkyE6ZIb4K+Lzc6GuAPfv38+UlBTa2dlx5MiR1NbWZnBwMDU0NKijo8P9+/eT/PSE\nrqenx6ioqC/O8eTJE/bv35+GhobcvHlzmdMAdujQgStXrizTvRVBcnIyp02bRm1tbXp6evLixYv8\n+PEjDQ0Nee3aNUnLqxBmzZolsqIdlZmnT59SXV29yGOLkhAbG0s3Nzfa29vz0qVLIlT3ycAvWLCA\n+vr67Nq1K6tUqcJly5ZRR0eHmzZtorKyMlu1aiVL1CLFyAzxV0RenHBQs2asAlBJSYk9evTgtWvX\naG1tTX9/fxoaGnLmzJm0t7enjo4Oz549S5LctGkTGzduXOKz28jISDZt2pSNGzcuU6KOq1ev0sDA\nQOqTzL97947Lly+nubk5W7Zsyd49e3Kkick3kQDh8ePH1NDQKNfux9dCx44dy3RmnpmZyfnz51NL\nS4uLFi0SeZa29PR0du/enU2aNOGaNWsoJyfHvn370tjYmFevXqWxsTENDQ2/mUpilRWZIf5KmTVr\nFgHQ0NCQWlpanDNnDs3MzOjl5UV7e3u2atWK/v7+1NHR4fXr15mbm8vvvvuOq1evLvEcubm53Lp1\nK42MjNinT59SF5Pw9vZmUFBQad+aRMjOzmZwcDD7q6sLdx2WtGrFp0+fSlqaWPH09CzV78TXSkhI\nCF1dXUt1z7lz52hra8uOHTsyLi5O5Jru3LlDGxsb4RGKoqIi69WrR1tra27s1o2jzcxYTVFRLHPL\nEC0yQ/wV08XdXVgWsbGDA+vXr09dXV22atWKbm5u1NbW5sKFC2loaMiHDx/yxo0b1NHR4YsXL0o1\nz9u3bzl16lRqaWlx3rx5Ja5YdPfuXWpra5c7cUFFsuuz7f/hBgbCmExvb28uXryY586dk/pVfmk4\nfvw47e3tv5kY6qLIyMgocUxxamoqhwwZQkNDQ+7evVss/3chISHU1tbmxo0bmZOTQ319faqqqrJl\ny5bc+pmvyK9t2oh8bhmiR1Z96Sumt4YG9gPY8uYNrGNj0bdvX2RnZ+PKlStISEhA06ZNsW3bNkya\nNAnt27eHnp4e+vTpg0mTJpVqnho1amDevHm4fPkyrl27BltbW4SGhoJksffVqVMHHh4eWLJkSTne\nZcXSOigIWz094a2khMAjR5CamoqTJ0/Cw8MDcXFxGDt2LLS1tdGoUSMEBATgzz//xD8lWcTmAAAg\nAElEQVT//COSilCSwNXVFR8+fMDFixclLUWiqKiowMfHB1u3bi2yD0ns2LEDdnZ2UFJSQkxMDLy8\nvCAnJycyHTk5OZgwYQImTJiAY8eOYeDAgXBzc8OLFy/g7OyMEydO4OnTp8L+xoaGIptbhhiR9JOA\nDPHxX+ctVVVVnj59mh07dqS8vDzV1NTYqFEjTpo0iTNnzmTDhg2ZlJREIyMjnj9/vszznjx5kvXr\n12ebNm2+6HUbHx9PTU3NSpflZ+jQoZw5c2ahbR8+fGBERASXLFnCnj170szMjOrq6mzfvj1nzJjB\nI0eOiDxrmThZvHgx+/TpI2kZEufy5cusVatWoSvcBw8esG3btmzQoAEjIyPFMv/z58/p7OzMdu3a\nMSUlhST5008/EQC9vb35/v17jhgxgrVMTbnG3f2b8GP4WpAZ4q+YPOetED8/zp0xg/Ly8lRQUGBo\naCiXLVsmfK2hocHw8HCOGDGCzs7O3LJlC+3t7cvlZZmdnc1Vq1ZRV1eXw4cPF35xFMbo0aM5ZsyY\nMs8lCe7du0cdHZ0Sb0M/e/aM+/fv55QpU+js7MyaNWuyVq1a7N27N1esWMHLly9LZTgX+SnTmpqa\nWrE/w28BgUBAW1vbfA+pmZmZnDdvHrW0tLh48WKxeSZfunSJJiYm/Omnn4QOldu2bSMAdu7cmQ8f\nPmSjRo3o6ekp8rAoGeJHZoi/Ie7fv08tLS0C4PDhw/n3339TWVmZAKihocHU1FR6eXnR09OTLi4u\nXLp0abnnTE1N5ejRo6mjo8Nly5YV6jWanJxMTU3NMlVHkiRdunQpsyNTTk4Ob926xfXr13Pw4MGs\nX78+q1WrxubNm3Ps2LEMDg5mXFyc1JzN9uvXj4sWLZK0DIkzffJkjrOyYoifHw/s3s26deuyc+fO\nYstHLRAIuGbNGmpra+dLDXvkyBECYIMGDXjgwAHq6Ohw6dKlUvP7IqN0yAzxN0ZWVhY7depEALS2\ntmZwcDCrVatGZTk59lRS4io3N7Zu0YI+Pj7U0tLikydPRDLvnTt32K5dO9atW5d///13gfapU6eK\nJGlCRXL27FlaW1uXOZb6v7x584anTp3iggUL2LVrV+rr61NXV5ddunTh/PnzGRYWJrEwlIsXL9LK\nykpk77Wysrl7d+FxT6+qVRkaGio24/fhwwcOHDiQtra2vHfvHslPxRwWtWzJbgB11NU5fvx4mpiY\nMCIiQiwaZFQMcuQXPGpkfJWsWbMGAQEBUFRUxOjRo/Ho118R+q9D0cTatXGiWjVUr14dJiYmCA4O\nFsmcJHH48GH88MMPsLW1xZIlS2BtbQ0ASEtLg7W1Nc6fPw8bGxuRzCduSKJp06aYOnUqunbtKpbx\nk5KScOnSJURGRiIyMhLR0dGwsLBA06ZN0axZMzRt2hR2dnZQUFAQ+fz/1dKwYUP88ssvcHNzE+tc\nkkQgECA5ORkJCQlITEwU/sl7bX33LnZnZQEAtnt7o3dIiFh0xMfHo0ePHrCyssLGjRtRo0YNAMDc\nRo0w/do1AIC/piYeOzlh+/bt0NbWFosOGRWDzBB/w9y8eRPNmjVDVlYW+qqqYtPr1wCAfqqqOFmj\nBgQCAbKyshASEgJXV1eRzZuZmYnly5dj4cKFGDRoEKZNmwZVVVUEBQUhOjoau3btEtlc4iYkJAQr\nVqzAuXPnKmS+7Oxs3Lx5E5GRkUID/fTpUzg5OaFp06ZCA21gYCDyudeuXYtjx45h3759Ih+7onj/\n/j2SkpKKNLRPnjyBpqYmTE1NhX/MzMyE/1ZRUMAfnp7Q0taG/5490DMxEbnG48ePo2/fvpg4cSLG\njRsHOTk5CAQCDB06FKkbN2Lfv1/Z0+vXx6zoaLE/hMkQPzJD/I2Tnp6Ohg0b4ml8PDrJyyMrNxdV\nAEBODuHVq+NtZiY0NTWRmJgIJSUlkc6dnJyMqVOn4ujRo5g3bx68vb1hY2ODw4cPw9HRUaRziYuc\nnBxYW1tj586daNq0qUQ0vHr1CpcvXxaumiMjI1G9evV8q2ZHR0dUq1atXPO8e/cOpqamuHnzJoyN\njUWkXnQIBAK8ePGigHH9/PX79+9hYmKSz7h+/m9jY2OoqKgUO8/06dOhqKiIGTNmiFx/UFAQVqxY\ngeDgYLRp0wbAp92iFi1a4N69e1CvVg3tsrPxXcuW8Nm8WSwPAjIqHpkhloHs7Gx07twZJ0+eRBcS\ne//dou4uJ4eD/z6NN2rUCJGRkWJ5+r5y5QrGjBmDzMxMNGnYEDxyBC7OzmgVFFQpvmiWLVuGiIgI\nhIhpm7K0kERsbGy+VXNMTAzq1KmTzzhbW1tDXr50qQQCAwOho6ODWbNmiUd8MXz8+LHQVWzen6Sk\nJKiqqhYwrp+/1tHRKXdc78qVK3H79m2sXr1aRO/s0wNx//79kZycjD179ggfdHbu3IkBAwYgJycH\nurq6MDc3x+7du2FkZCSyuWVIHpkhlgHg05f32LFj8fj334Vnxb+2aoVjyso4e/YssrKyoKqqip07\nd8Ld3V0s8wcHB+PQ4MEIzsgAAOz284P39u0in0vUvH37FhYWFoiKioKFhYWk5RRKRkYGoqOjhSvm\nS5cuIT09XbidnfdHS0ur2HFu3bqFDh06ID4+HoqKiiLTRxIvX74scss4MTER6enpMDY2LnQla2pq\nChMTk3Kv+ktCaGgotm3bJrIt+jt37sDT0xMuLi747bffoKysjEePHmHYsGE4c+YMAEBVVRWDBg3C\nzz//LNL/dxnSgcwQy8jHgjlzcHXuXOTm5uK4ggKux8QgMTER3bt3x9u3bwEAdnZ2+PPPP8WyfRzs\n4wPff1eWI42N4bllC5ydnUWanUgcTJ48GRkZGVi2bJmkpZSY58+f51s1X7lyBbq6uvlWzQ0aNChw\nJNGyZUuMGzcOnp6eJZ4rMzMTjx8/LtLQJiUloVq1akWezZqamkJPT6/UK3hxcOHCBfz4448iyTYW\nEhKCwMBALFq0CAMGDEBGRgYWLlyI3377DR8/foScnBxUVFSwadMmdOvWTQTqZUgjMkMsowB79uyB\nv78/3r17B0VFRVy8eBHW1tYwNDTEu3fvoK6uLlxNLV++HE5OTiKb+3lSEs5Mnoy/jx1DposLbty9\nC4FAgMDAQPTr1w81a9YU2Vyi5MmTJ6hfvz5iY2OhqakpaTllIjc3F3fv3s23an748CEaNGiQzxHs\nrwMHcDMoCK4uLmgVFARdY2O8evWq2LPZV69ewdDQsMgtYxMTE6FnsLTz6NEjuLi4ID4+vsxj5OTk\nYNKkSdi7dy9CQ0Ph6OiIQ4cOYcyYMTA0NERkZCSqVq0KMzMzHDhwAJaWlqJ7AzKkDpkhllEoZ8+e\nRbdu3ZCeng4FBQXs27cPioqK6NWrF4BPX9okkZOTAycnJ8yePVukK9eLFy/C29sbd+/exdWrV7Fy\n5UqcPHkSfn5+CAwMRN26dUUyjygZ8L/27j4sqjLvA/iXmYFhQGFgeEcmMURAXjVTaUMtbcVJNMVo\nMVt3N3bXbH1JLN0e27aXldJsddXtMvfxZfFda9UUKlOjsLLaxLdUEkwSAUEIBAacOb/nD3WerDQj\n4JB8P9d1Lgc4zPkdgfM9933uc+6JE9GrVy/Mnj1b7VJazYULF/DJJ584Ws0ffvghkqqrsaGpCcCl\nEfav2e1wdna+7rXZgICAm2Z0b0NDA7y9vdHY2Nii3/fy8nKkpaVBr9dj7dq1qKmpwdSpU1FYWIjB\ngwdj+fLlMBgMSE9Px6JFi35w8BjdBNrtjmX62Tl48KB4d+kiY7RaGQXInFmzJDU1VQYMGCDx8fGS\nmJgofn5+YjAYxM/PT/r16yevvfZaqz30Yfz48fLkk086Pi4pKZE5c+aIv7+/3HXXXfLaa691qMnO\nCwoKJDAwUKxWq9qltBlFUeR/R492PNRi2YgRnXKuWw8PjxbNGvbBBx9It27d5Mknn5S6ujp56qmn\nxNdolKfi42Va9+7iDIibm1uL5j6mny8GMV3X/44addXEEcnJyWIymWTUqFFisVhk3bp1EhQUJD4+\nPmI2myUyMlIiIyNl5cqVP3kS9JKSEvH29v7OfKpWq1XWrFkjiYmJEhISIs8//7yUl5f/pG21lnvu\nuUdWrFihdhltquz0aVlz//1yv14ve99+W+1yVBEeHi5Hjx694fUVRZGlS5eKr6+v/Oc//5Ft27ZJ\naGiopKamysK77rrqaV2HDx9uw8qpI1J/5AN1aF3c3R2vBUBubi4uXryI8+fPQ1EU5OTk4MSJE3j0\n0UdRVVWF06dPIzw8HCtWrEBYWBgWLVqEhoaGFm27W7dumDZtGmbOnHnV5/V6PdLT05Gfn4+tW7ei\nqKgIvXr1woQJE/DRRx/94PSLbSkzMxMvvfSSqjW0Nf+QEKRv2ID+f/sb/r5kidrlqCIwMBBlZWU3\ntG5jYyN+85vfYMmSJVi7di1effVVzJw5E5MnT0Z9fT32XB4ZDQAjR45E796926hq6rDUPhOgju2b\nMzht27xZ/P39BZcyWYYOHSr9+/eXzMxMERE5ffq03HfffeLm5iYmk0lefPFFGTNmjPj5+cmzzz4r\n58+f/9Hbb2hoELPZLHv37r3uelVVVTJv3jwJDQ2V2267TVasWCENDQ0t2uefQlEUiY2Nldzc3Hbf\ndntraGiQ4ODgNpv2ryNLS0uTtWvX/uB6RUVFkpCQIKmpqTJr1iwxmUwyfvx46du3r/To0UN8fX3F\nTaeTZ/r2lQ2ctrDTYhDTj9Lc3CyZmZmiAyQFkNEajfh5ecmcOXMc6+Tl5cmtt94qBoNBhgwZIu+8\n845MnDhRvL29JTMz80dPJLFhwwaJj493TP92PTabTd544w0ZPny4+Pr6yhNPPPGdru22tmrVKhk6\ndGi7blMtr7zyigwbNkztMtrd1KlT5aWXXrruOrm5ueLn5ye//e1vxWw2S0JCgoSEhEhiYqIsWLBA\nunbtKkajUfbv399OVVNHxSCmFlk8bJjjulaqTidOTk5y5513ypEjR0TkUiAuWbJE3N3dRa/Xy3PP\nPSdFRUUyZcoU8fLykoyMDCksLLyhbSmKInfeeacsW7bsR9V44sQJmT59unh7e0tKSoq89dZb7TJ7\nUFNTkwQHB8tnn33W5ttSW3Nzs/To0UP27NmjdintKisrS2bOnPm9X7Pb7fLcc8+Jn5+f9OnTR0wm\nk3h4eMh9990n+fn5kp2dLTqdTuLj41vUS0Q3HwYxtcjG9PSrBnFptVpxcXERo9EoycnJsmvXLlEU\nRaqrq+Whhx4SZ2dn6datm3z44Ydy7tw5mTNnjvj4+EhaWtoNBdann34qAQEBLZr0/MKFC7Js2TKJ\njY2V8PBwWbhwYZtPnp6VlSUTJkxo0210FKtXr5bExMRONRfuypUrv/fnW1NTIxaLRfz9/cXFxUVc\nXV0lIyND8vfulY3p6fJ4eLi4APKHP/yh008pSf+PQUwt8s1rxxuzs8XLy0sAiKtGI09ERMgEDw+J\njYiQVatWSVNTkxw5ckSio6NFq9XK6NGjpba2Vmpra2X+/PkSFBQkycnJ8u677173YP673/1OZsyY\n0eKaFUWRvLw8SUtLE6PRKH/84x/l0KFDLX6/66murhYvLy8pKSlpk/fvSGw2m0RFRcmOHTvULqXd\n5ObmfufyQ0FBgZhMJtHpdOLi4iLTp093jObPHjfOceL6QmKiGiVTB8YgplahKIrMnDlTRl0+2Agg\nLw8eLHfddZcEBQXJ3LlzpaqqSrKzs6Vr165iMBhk6dKlInLpdqRly5ZJWFiYJCYmyvbt2783kMvK\nysTHx0eOHz/+k+s9c+aMPP300xIYGCiDBg2STZs2/eTbrb5t2rRp1+y+vNls2bJF4uPjO00rr6Cg\nQKKjo0Xk0onIxIkTxcnJSbRarUyZMkXq6+tF5NIAxszMTBnn4uL4u9iYnq5m6dQBMYipVf07NfU7\nB5zPPvtMHnzwQfHy8pI//elPcvToUcnIyBCNRiM9e/Z03I9ps9lk/fr1Eh8fL9HR0ZKdnf2dB3a8\n+OKLcu+997ZavU1NTbJ+/Xq58847JTg4WJ555hk5e/Zsq7x3cXGxeHt7d4oHXiiKIn379pWNGzeq\nXUq7KC8vF5PJJPPnzxe9Xi8AxGKxOEbqf/N3ftq0abI/P9/Rg8SR0fRtDGJqVd/ssv72AaekpEQe\nf/xx8fb2lrFjx8qWLVskLi5ONBqNpKakyLoHHpCN6ely9ssvJScnR5KSkqR79+6yZMkSxwHOarVK\nWFhYm9wedODAAcnIyBCj0Sjp6emSn5//k697pqWlyYIFC1qpwo4tNzdXevXq1aGedtYWKioqZMaU\nKZJyeXyEr9EoBw8eFEVRJCcnR+6++24JCgqSrKwsDsaiG8IgpnZXV1cnCxculNDQUBk4cKDMmDFD\nxup0jpb0+m903eXn58vIkSMlICBA5s6dKzU1NbJ161aJjIxs9a7kK86fPy8LFiyQsLAwSUhIkOXL\nlzu6Gn+sjz/+WMxmc5vV2pEoiiJJSUk37ZPFTpw4IRMnThS9Xi+jNRrH7+va+++XFStWSHR0tMTE\nxDjGRRDdKAYxqebixYuyadMm6d+/v4zv0sVxYBuj1UpmZqZUVlY61j148KCMHz9evL295YknnpBB\ngwbJokWL2rQ+u90uOTk5YrFYxGQySWZmppw8efJHv09SUtINPfzhZvDee+/JLbfcctM8b1tRFFm9\nerVERUWJTqcTjUYjAK665vuAq6sMGzZM3nzzzU41cpxaD4OYVKcoimzfskUe7dZN7ndxkT4xMaLX\n68XZ2VnGjh0rH330keMAV1RUJI888oh4eHiIwWCQ//73v+1S48mTJyUzM1N8fHzEYrHIzp07b3hg\n0rZt26RPnz6d5iA9fPhwWbx4sdpltNiFCxfk9ddfl3vuuUdcXFxEp9NJWFiY6PV6MRgM8swzz0jG\nxImS6uIi00JDZXcneIoatS0GMXUohYWFMnnyZPH09JTevXuLq6ureHh4SExMzFVdxGVlZdKnTx/R\n6/Xy4IMPttltSN9WX18v//rXvyQhIUHCwsJkwYIFP3gd0G63S69evTrNQy8++eQTCQwMbHF3vhq+\n+OILWbhwoQwdOtQRuCEhIfLYY49J9+7dRavVysiRI2XMmDFiMpnk8ccf7xS3plH7YBBTh1RVVSXP\nP/+8+Pv7S0hIiLi5uUnPnj3Fy8tLpk6dKseOHZPKykrx8fGRqVOnSkBAgIwcOVL27dvXLvUpiiL7\n9u2T9PR0MRqNkpGRIQcOHLjm+suWLROLxdIutXUEY8eOlRdeeEHtMq6pqalJ3n77bZk+fbqEh4eL\nn5+fJCQkiIeHhwwfPlx27twpycnJ4uTkJD179pR+/fqJ2WyWl19+WWpra9Uun24yTiI38TQx9LPX\n1NSEdevWISsrC2VlZbDb7bj99ttx6NAhxMTEIDQ0FMXFxdi+fTtWrlyJefPmwWw2Y/bs2fjlL3/Z\noonbf6zy8nK8+uqreOWVV9CjRw9MnjwZY8aMgbOzs2Mdq9WK7t27Y8+ePYiMjGzzmtR29OhRDB48\nGIWFhfD09FS7HADAmTNnkJOTgx07dmD37t2IjIzEgAEDUFpaindycjDRzw8RkZH4omdPvLxkCVxd\nXdG1a1cEBgYiMzMTqamp0Ol0au8G3YzUPhMguhGKosibb74pAwcOFL1eL56enjJx4kRJTEwUnU4n\naWlp8tVXX8nFixclOztboqOjJSEhQTZs2HBDk0W0hubmZtm8ebMMHjxYAgMD5S9/+ctVE1z89a9/\nlYcffrhdaukIHnroIXnqqadU277NZpP8/Hz585//LPHx8eLl5SUPPPCArF69WnJycmTs2LHi4+Mj\nTz75pPwrJeWqR7a6ubmJxWKRvXv3dppr+6QeBjH97Bw8eFCSk5NFq9WKt7e3pKakSKpOJ6nOzpIy\nfLjs2rVLbDabbNu2TQYOHChhYWGybNmydh3Je+jQIZk0aZIYjUZJS0uTvLw8qaioEKPRKGVlZe1W\nh5pOnjwp3t7ecu7cuXbbZmVlpaxZs0bS09PFZDJJbGyszJo1S9577z1pamqSrVu3yi9+8QsJCQmR\njIwMmTBhggQHB1/1RLhpoaGOh8wQtQd2TdPPVmlpKaZPn47mjRvx+uXPPdOnDzY1N6O5uRmTJk3C\nQw89hMOHD2Pu3Lk4ePAgHnvsMfz+979H165d26XGr7/+GqtXr8bixYuhd3JC/+pq+Pr4YGpuLvxD\nQtqlBjU98sgjcHd3x7x589rk/UUEBQUF2LFjB3bu3IlDhw5hyJAhsFgsSE5ORkhICKxWK/7xj3/g\n5Zdfhs1mg7OzMyorK+Hi4gJFUWA0GmGtrcWk7t0RHh6OX/79753iZ0MdB4OYfvbW3H8/xm/aBAAY\nq9WiMCoKZrMZFRUV+Pzzz5GamopHH30UGo0GWVlZ2L17NyZNmoQpU6bAx8enXWpUFAWLhw3DlN27\nAQDj3d1xKi4O/v7+CAgI+N7F398fer2+XeprK6WlpYiOjsbhw4cRFBTUKu9ZV1eHXbt2YefOndi5\ncycMBgMsFgssFguSkpKg1+tRVFSEnJwcrFixAgUFBXByckJISAgURUFNTQ1uu+02lJaWoqmpCTNm\nzMCvf/1ruLu7t0p9RD8Wg5h+9spLSpA3axYAYMCzz6Kiuhr79+/H/v37sW/fPpw6dQpOTk4wGo0Y\nPXo0UlJS8Prrr2PLli2YMGECZsyYAbPZ3OZ1bho/HuPWrgUAzI6IwL3Ll6OsrOyqpby8/KrXbm5u\n1w3qK699fX2h1WrbfB9aYubMmaivr8fSpUtb9P0igsLCQker98MPP8SAAQNgsVgwYsQIhIWF4ciR\nI3jvvfeQl5eHvXv34sKFC2hubobZbIa7uztOnTqFESNGwGg0Ijc3F0FBQcjMzERKSkqH/X+jzoNB\nTDe92tpa7N+/H9nZ2cjNzcW5c+fg4uKCuLg4aDQaHDp0CBaLBU8//TQiIiLarI4rJwxf19bif/bt\nQ+GpU9ftIhcRVFdXXzOov7mcP38eJpPpe0P624vRaGyX0eRXVFZWolevXvj444/Ro0ePG/oeq9WK\nd999Fzt37sSOHTvQ2NjoCN6kpCQUFhY6gjc/Px/e3t6IjIzEmTNncOzYMQT7+iL2zBkEBgQgdPJk\nlJ0/jxUrViApKQkzZsxAYmJiG+810Y1jEFOnU1RUhPnz52PNmjUwmUxwc3PDF198AWlqwhiDAWaz\nGXc89xzuTk5us+7K9PR0xMbGYtbllvxPZbPZcO7cueuG9ZXFarXCz8/vhlrarbX/Tz/9NIqLi7Fq\n1aprrlNSUuII3r179yImJgYWiwVDhgxBY2Mj3n//feTl5WH//v0IDQ1FUlIS7rjjDjQ0NGDBggUo\nLi6Goii47bbb0L+sDPO/+AIAMM7FBf4ZGZg2bRrCwsJaZX+IWhODmDotq9WKzZs345///Ce+/PJL\n/MrNDfMKCwEA43Q6vKHTITw8HLfffrtj6d27d6vcS/r5559j0KBBOHnyZLsNHLvCarVeN6yvfO3s\n2bPQ6XTXDOlvLn5+fnBxcbnmNmtraxEeGopnBw6E0dMTSVlZMAUG4oMPPnCEb2lpKYYPH45Bgwah\nS5cuKCgoQF5eHgoKChAbG4ukpCTE9OoF5ORAq9Fgt6sr1m7ejMbGRuj1ejg7O6OhoQF2ux33aTR4\nTVEAAP9OTcWEy2MIiDoiBjERgAMHDiB73DhHK+o+Jyds12jg6uqKMLMZ/c+fR0N9PbbZ7Yjt2/eq\ncO7evXuLunrT09MRExOD2bNnt/butAoRQV1d3XWvY19ZKioq4OHhcd1W9jtTpuCF48cBAFPNZvy7\nrg633HILBg0aBF9fX5w9exbvv/8+Tp48idtvvx133HEHoqKiYDAYcODAARw8eBBdd+/GypoaAMAo\nADu0Wvj7+yMiIgJ9+vTBkCFDMGTIENRWVjrGDSRlZXEUNHVoDGKiy8pLSrBrxgwUFxdjj6srjhQW\noqKiAiNFsPXyOqMAvG0wwN3dHTqdDvX19XByckJcXBySkpKQmJiIfv36wdfX9we3p2aruLUpioKq\nqqrvDeqSkhIcPnwYEcePY4vdDgB4rEcPfNW3Lz7++GNUVVUhNDQUXl5eaGpqQlVVFcrLy1FXV+c4\nwfH09ERQUBAGff01lnz1FQBgXVoafrV+vWr7TNRaGMRE12G327F85Ej8IScHAPAboxFv6HSoqamB\nzWaDk5MTvvkndCU4nJ2dERwcjPj4eNx9990YMWLE97acO3qruKUaGxuxefNmrFq1Cvn5+TCZTNDY\nbOh37hwUEbyl1cKqKDAYDNDpdGhubobNZkNwcDDCw8MRFxeHxMRExMXF4ZZbboFGowFw9Qh5tnTp\nZsEgJvoB1zr4W61WlJWV4dSpU46u06NHj+L06dOoqqpCc3Pzd95Lq9WiS5cu6NatGyIiItBVr8eF\njRtx99ChGLF4Mcy33tqu+9ZSIoLKykr896OPcPiFF1Df0ICSmBgcOnYMn3/+OWprawFcOjG5ss86\nnQ5WqxV2ux0RERGIi4tDdHQ0oqKi0Lt3b3Tr1s0RuESdCYOYqI2ICEpLS3HkyBF88MEH2LNnD44e\nPYrq6mrYbDYAl7q6/3N5/VEA3tBooNfr0aVLF3h6esLHxweBgYGO2468vLxgNBrh6uoKAFe1sK+8\nvta/P2YdRVFQXV2Nc+fOOZbKykpUVFQ4PtbpdEhRFKyzWgEAowFsBaDT6aDRaKDVatG7d28kJCQg\nKirKsQQHB7fr7VNEHR2DmEgFdXV1+PTTT/HW5Mn429GjAIApZjMO33orqqurUVNTgwsXLqChoQFW\nqxUi4mgtKooCJycn6PV6GAwGdOnSBW5ubnB1dYXBYLjU3SuC8MJCaDQalMTGwsXd3dGFLiJQFAX1\n9fWOpa6uzvG6sbERzc3N0Ol00Gq10Gg0jnC22Wyw2+0QETg5OV11/TzDzw+3ToC//t0AAAKnSURB\nVJ+Ofv36ISoqCgEBAQxcohvAICZS0Y1e86ytrUVxcTGKiopw8uRJHDt2DMePH0dxcTEqKipgMBjQ\ntWtXODs7Q0QwoLwc6y+3VMdqtXjbzQ0iArvdjosXL8Jms0Gr1cLJyQl2ux0GnQ73ajRw1etxIiwM\nXby9YbfbUV9fj6qqKpSWlsLV1RXh4eGIjY1FXFwcevfuDR9PTxydPx9OP1A/EV0bg5joZ85ut+PM\nmTMoKipyLHXZ2Vj45ZcALt2K9YZWC7k025pjAeAYbJYCOFq2v3J1xbnLtw5duX4bGRnZbs/lJups\nGMREN6EfamkrioLm5mY0NTWhqakJ2x9+GL/bvh0AsCk9HePWrGn3mok6KwYxEfG2ICIVMYiJiIhU\nxJv2iIiIVMQgJiIiUhGDmIiISEUMYiIiIhUxiImIiFTEICYiIlIRg5iIiEhFDGIiIiIVMYiJiIhU\nxCAmIiJSEYOYiIhIRQxiIiIiFTGIiYiIVMQgJiIiUhGDmIiISEUMYiIiIhUxiImIiFTEICYiIlIR\ng5iIiEhFDGIiIiIVMYiJiIhUxCAmIiJSEYOYiIhIRQxiIiIiFTGIiYiIVMQgJiIiUhGDmIiISEUM\nYiIiIhUxiImIiFTEICYiIlIRg5iIiEhFDGIiIiIVMYiJiIhUxCAmIiJSEYOYiIhIRQxiIiIiFTGI\niYiIVMQgJiIiUhGDmIiISEUMYiIiIhUxiImIiFTEICYiIlIRg5iIiEhFDGIiIiIVMYiJiIhUxCAm\nIiJSEYOYiIhIRQxiIiIiFTGIiYiIVMQgJiIiUhGDmIiISEUMYiIiIhUxiImIiFTEICYiIlIRg5iI\niEhFDGIiIiIVMYiJiIhUxCAmIiJSEYOYiIhIRQxiIiIiFTGIiYiIVMQgJiIiUtH/AQmz8uKnIGO4\nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1175c6a10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import networkx as nx\n",
    "n = 10\n",
    "ex = np.ones(n);\n",
    "lp1 = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr'); \n",
    "e = sp.sparse.eye(n)\n",
    "A = sp.sparse.kron(lp1, e) + sp.sparse.kron(e, lp1)\n",
    "A = csc_matrix(A)\n",
    "G = nx.Graph(A)\n",
    "nx.draw(G, pos=nx.spring_layout(G), node_size=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Fill-in minimization\n",
    "\n",
    "Reordering the rows and the columns of the sparse matrix in order to reduce the number of non-zeros in $L$ and $U$ factors is called (**fill-in** minimization)  \n",
    "\n",
    "is based on graph theory:\n",
    "\n",
    "\n",
    "- **Minimum degree ordering** - order by the degree of the vertex\n",
    "- **Cuthill–McKee algorithm** (and reverse Cuthill-McKee) -- order for a small bandwidth\n",
    "- **Nested dissection**: split the graph into two with minimal number of vertices on the separator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Take home message\n",
    "- CSR format for storage\n",
    "- Sparse matrices & graphs ordering\n",
    "- Ordering is important for LU fill-in"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Next week\n",
    "- Iterative methods and preconditioners\n",
    "- Matrix functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Questions?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>\n",
       "<style>\n",
       "    @font-face {\n",
       "        font-family: \"Computer Modern\";\n",
       "        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "    }\n",
       "    div.cell{\n",
       "        /*width:80%;*/\n",
       "        /*margin-left:auto !important;\n",
       "        margin-right:auto;*/\n",
       "    }\n",
       "    h1 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\n",
       "    h2 {\n",
       "        font-family: 'Fenix', serif;\n",
       "    }\n",
       "    h3{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "        margin-top:12px;\n",
       "        margin-bottom: 3px;\n",
       "       }\n",
       "\th4{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "       }\n",
       "    h5 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\t   \n",
       "    div.text_cell_render{\n",
       "        font-family: 'Alegreya Sans',Computer Modern, \"Helvetica Neue\", Arial, Helvetica, Geneva, sans-serif;\n",
       "        line-height: 1.2;\n",
       "        font-size: 120%;\n",
       "        /*width:70%;*/\n",
       "        /*margin-left:auto;*/\n",
       "        margin-right:auto;\n",
       "    }\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "    .text_cell_render h1 {\n",
       "        font-weight: 200;\n",
       "        font-size: 50pt;\n",
       "\t\tline-height: 110%;\n",
       "        color:#CD2305;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\t\n",
       "    .text_cell_render h5 {\n",
       "        font-weight: 300;\n",
       "        font-size: 16pt;\n",
       "        color: #CD2305;\n",
       "        font-style: italic;\n",
       "        margin-bottom: .5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\n",
       "    \n",
       "    li {\n",
       "        line-height: 110%;\n",
       "    }\n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "\n",
       "</style>\n",
       "\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"]\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"./styles/custom.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}