{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Sparse Matrix Factorizations and Fill-In"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import scipy.linalg as la\n",
        "\n",
        "import matplotlib.pyplot as pt\n",
        "\n",
        "import random"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Here's a helper routine to make a random **symmetric** sparse matrix:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def make_random_sparse_matrix(n, row_fill):\n",
        "    nentries = (n*row_fill) // 2 # because of symmetry\n",
        "    data = np.random.randn(nentries)\n",
        "    rows = np.random.randint(0, n-1, nentries)\n",
        "    cols = np.random.randint(0, n-1, nentries)\n",
        "    \n",
        "    import scipy.sparse as sps\n",
        "    \n",
        "    coo = sps.coo_matrix((data, (rows, cols)), shape=(n,n))\n",
        "    \n",
        "    # NOTE: Cuthill-McKee applies only to symmetric matrices!\n",
        "    return (100*np.eye(n) + np.array(coo.todense() +  coo.todense().T))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next, we will take a look at that matrix from a \"birds eye view\". Every entry with absolute value greater that $10^{-10}$ will show up as a 'dot':"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "794 non-zeros\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "<matplotlib.lines.Line2D at 0x7f32dd4099b0>"
            ]
          },
          "execution_count": 9,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQQAAAD7CAYAAACMu+pyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF3NJREFUeJztnV2obdV1x///IHlIAtaKXkFtEoxYWwImtFYChR0KfvXh\nhoJi+9CYIgixL32qPt37mDwklBAkkIjcgGLNQxKFEq3Y8+BLYrHWFI25ebjVq9wToRponkwdfdh7\nn7Pvcu2951prfowx1vjBgX3W2WetMcYc8z/nWmt+UEQQBEEAAB9pbUAQBHoIQQiC4IgQhCAIjghB\nCILgiBCEIAiOCEEIguCIJoJA8naSvyD5S5L/2MKG3JA8R/I/Sf4HyZ+tjl1G8lmSr5N8huSlre1M\nheQjJA9JvrJxbKs/JL9F8izJl0ne1MbqdLb4d4rkeZIvrX5u3/jbQyv/XiN5axury1NdEEh+BMC3\nAdwG4I8B/DXJP6xtRwE+ALAQkc+JyM2rYw8CeE5EbgDwPICHmlk3nEexLKNNev0heQeA60TkegD3\nA/hOTUNH0ucfAHxTRD6/+vkJAJC8EcDdAG4EcAeAh0mynqn1aNFDuBnAWRH5bxF5H8ATAE42sCM3\nxIfjeRLAmdXnMwC+VNWiCYjICwDe7Rzu+nNy4/j3V//3UwCXkjxRw86xbPEPWJZjl5MAnhCR34nI\nOQBnscxjd7QQhKsBvLnx+/nVMesIgGdIvkjyvtWxEyJyCAAicgHAFc2sy8OVHX+uXB3vlulbsFum\nD6xue763cUvkyb+dtBCEPgX2MH76CyLyJwDuxDKp/hw+/ErBS5k+jOWtz00ALgD4xuq4F//20kIQ\nzgP4g43frwHwdgM7srJqMSEi7wD4EZZdysN115nkVQB+3c7CLGzz5zyAaze+Z7JMReQdOZ7c810c\n3xa48C+FFoLwIoDPkPwkyY8CuAfAUw3syAbJj5H8xOrzxwHcCuDnWPp17+prXwbw4yYGjoe4uHXc\n9OdeHPvzFIC/BQCStwB4b31roZyL/FuJ3Jq/AvBfq89PAbiH5EdJfhrAZwD8rJqVFbmk9gVF5P9I\n/j2AZ7EUpEdE5LXadmTmBIAfkhQsY/qYiDxL8t8BPEny7wC8AeCulkYOgeTjABYALif5BoBTAL4G\n4Addf0TkX0jeSfJXAH4L4CttrE5ni39fXL0y/QDAOSzfmEBEXiX5JIBXAbwP4KsbPQlX0KlfQRCM\nIEYqBkFwRDFB8DgaMQi8U+SWYTUa8ZcA/gLLp7EvArhHRH6R/WJBEGSjVA/B62jEIHBNKUHwOhox\nCFxT6rXj3pFdq1d0QRA0QER6J2eV6iEkjUY8deoUls8wBCL+ftb+dX88+LvNNy8/nv3bRSlBGDQa\ncY+N7pibv4EditwyyIjRiGRUlODDRF7UpdjQZVkuLnHDru8sFouN7/sr/E3/PPtWklYxK+1fNx+0\n5Eezocskpe/aWgJTmtp+ziWuwX5IQio/VBzNuqcwlCkLWrVYDKt25QwxCFJQJwjAcfIOqahTEr50\nZfG5+p5eNMVbky0pqBSENWN7C9ro3iuW8MlDnHIxVuBLxNBaz0y1IAB+RGGNSJkksZZ4ORmaH9u+\nnyuGlvNVvSAA/kQhyMvQilxaPFPOrzWfTQgCkC4KWgMdlCdn2ZfOo1RRqp3PZgQBSAvinLvOcydn\n2WvJo9p2mBIEIHoAwfwo9SC6D3OCkPt5wpBzhRiVj8G+83sqg22+dI+XehDdhzlBAPKKwpBA1+6+\naUz+1g/ktHTlc7DNl5Y+mhQEIN48BEEJzAoCMG5EoyU8tYa76Cs/y2Vq6W1HF5WCMDQI0VsoQ62Y\n9gnfEDHUVvaW33aoFIQxQfAuCnOYgDWW3HZ6zqN9qBSEbewrKG2isGnLVLusVM5W5Cr39TTx7vla\n5lXNa5sShNSBSVpEYdPeqNBLSpXNmPj22bI+T/d8Lcuv5rVNCUIqUfn0oqlsNNmiBZeCAOjpJXhk\nV2zHxj3KSwduBUHbfaAndrWsY18F951Ta3nlfF5Rm33XdCsIwIdFQUMXUWuS5yZHrFuXV+l1E1r4\nt++argUB0Dd4qXWSt6J1/MdcX0tZxVuGTHR7B62Tcs60rlytrz+FeMuQib5XRyEK4/EUuzGL7Xjy\nfxuuBaGPEIXxWG5lu4xZbMeT/9uYnSAAIQrBhxkyqtRz7sxSEIB5qH2QzpBRpZ5zZ7aCAPhWek/k\nnBMS7MakIORcLUlDgmmwQTOe54RoK3uVgpAyqzEXGpZ395bkQ9BWIYC6Nmkre5WC0GIj1JoiVAuN\nlQ3QN3q0i0abaqFSEFqgbUTjLlJt1JrYJe2KyVXTCEHosKu3oCVpNFb0KbHJGdexsdn3f63KPtZU\nVMA2URiabFoEpAZTFijRKHBdWtk45bpj8i8EYQtdURize472RG8tWNrjY50x8Q1B2MGmKNTcPacW\nsfFM0CUEYQ/eRKAlEUv9hCAkEC1bMBSrOROCgP1vFXItxzbm7YXVxGqFlnhZ7Q2ZFITchd4tvG1v\nGKYOqBmzuafVxGpFxGsaJgWh1g7E2+bDa2mFSuDZtzFo3XeylA0mBaEFqcuxaUiWKbQYoKM5ZkP3\nnWy5H2YOQhASSV2OzXuXtYR/nmKm3ZdZL8NeGi3Tp3Nh1Rerdk+lxErSIQgTKSEKrd46aG/dtqHJ\nbutTp0MQMpBbFOKtQ3k0bTqriRCETFhPhLkR5dXPJVP+meQ5AL8B8AGA90XkZpKXAfhnAJ8EcA7A\n3SLym4l2moDMm2i5z1cLq3YH03sIHwBYiMjnROTm1bEHATwnIjcAeB7AQxOvoZKUwUtTsVapLE1n\nDvqZKgjsOcdJAGdWn88A+NLEa6hk16jDGg+WND5ZDyGwz1RBEADPkHyR5H2rYydE5BAAROQCgCsm\nXmMQGrbqrjGisaTwaBSboA6TniEA+IKIXCB5BYBnSb6OpUgkcfr06aPPi8UCi8Viojm6tupeV9pS\nLae18wLxfKEFBwcHODg4SPouJVPpkDwF4H8B3Iflc4VDklcB+DcRubHn+7Lv2l6Sx4sfnihVJhbK\nmiREpLcfOPqWgeTHSH5i9fnjAG4F8HMATwG4d/W1LwP4cfo5L/59yjp9mvA2otEDFntXNRjdQyD5\naQA/xPIW4RIAj4nI10j+PoAnAVwL4A0Ad4nIez3/v7eHoIkcym+h9Qj8s6uHkO2WYSjWBCEIvFDk\nlqEWU7va2rrqmuzRZEugAxWCMKcJO9ueJ7SonLlio0lYStuiwdeSNqgQhM3EzPFgUTt9opCyjFsf\nGhK0RRlt87vWalotmTLOxdx6CBoCXoNdC7cOefg4l3h1mavfm+Rc13ONOkGoTcuWuKv029ZytETL\nHkvta3u83uwFoUVL3He74GUYckqcNK9FMMS22sJd43qzF4QW9N0Dehq8tO+hqeYeUG7bUspUU7mH\nIDRk18KtuZKkReXbJnhzJMXvlr2qLiEIGRhbWLvWVLBegXL5oKn1bElfLEvEJgQhA2MTf8xOTlbQ\nNOvUK7HIqiP2bfQSLaMdPJWVWUEg9RVEjifUm68eh5xPWyzmhKdejFlBEKlbECkVLrc91t88aLW9\npl3WxiqYFYQ+PM6JSB2mWmqxjynUjlmqvTXtqr0P5FTfXAmCp65bl9TegoZVn3PaoHmg0FRabKy7\nD1eC0IqaO/626CkMJacNGvwphcayDEHoUGIDzZzUeK6g9d7fIrtiqVHsQhA6WBhMU1oUNE1ntoC1\nSr+LEIQC1EgCa4m2D8v+WLa9SwhCQTR17Vu1wC2uO+WamnsqMf3ZOKVfOVkYp9ByclXt/83BmNuP\nnDkQgtCInGP9WwyaysGcei2pTF0Fyf3AJM2FV5rU4dkpopAjjrnLovVgLy/knGauXhC8Fd4QusOz\nU7qTpRYf9TAlexvWG52c5aJeEIJjUhfbKD0k1noF6uJV6MYQguAQb+MUvAmQZkIQVnhLul2iYM3X\naMHrEYKwonXSlVjfwcJrydpEPHZzSWsDgiU1tydvLX4tmbPvKUQPYQZEq3gxGlfbSmGWIxUtFlQr\nNhN73yvJiOsxtVfbysUsN2qxWFCt2EzsfXFrLQqWli2zMoIylmEPBtHtGk/ZNXgqNZd4m3qtVhun\nDLU7lmEPBrGtazymt6DxlqNlb9JrTzYEYQclKoGWitUVBY3LeQX1CUHYQYlKoKlibYpCy9uJ3Hjw\nIRdDYxGCMHO6PYUcKwG3rpAtRbe1793rD41FCMKM2ewdpG4hl2s3431Y2+BkTdf32qtauZ/+HJRj\nM3lSX1/WorYdNUeK1rZhCCEIE2ndRcyJJ1+CcYQgjCDnCjWt2bxVaD14KWhPCMIIUkWgVOUqsV3b\n2F2nLeDNn5K4EQSNhZ6z91CzV+LpFSRgvxdXEzeC4L3QWy1n7kUUgjTcCEJQhhCF4ViOVwiCM1KT\nceg265aTvDa5NuhRuR08yUdIHpJ8ZePYZSSfJfk6yWdIXrrxt2+RPEvyZZI3lTJ8LF4TuzsEeR9j\nZtZ5jV1NhqxerXU7+EcB3NY59iCA50TkBgDPA3gIAEjeAeA6EbkewP0AvpPR1ix4fdYQG8zaQ2M8\n9wqCiLwA4N3O4ZMAzqw+n1n9vj7+/dX//RTApSRP5DH1GA0tlQYbcpG7G6slNlrsqEEuX8c+Q7hS\nRA4BQEQuALhydfxqAG9ufO+t1bGstFrAYqgNrRjq+9DhtSW7uiXGWMyBXL7mfqjYV5xbTfW0mYgm\naoxTKFV2cy43DYxdhv2Q5AkROSR5FYBfr46fB3DtxveuAfD2tpOcOnUap08vPy8WCywWi5HmBLVZ\ni0JUYP0cHBzg4OAg6buUhBIl+SkAT4vIZ1e/fx3A/4jI10k+COD3RORBkncCeEBE/pLkLQD+SURu\n2XJOEZFZJpU3n735s42pfmqJE0mISG8fb68gkHwcwALA5QAOAZwC8CMAP8CyN/AGgLtE5L3V978N\n4HYAvwXwFRF5act5JUWMWqCl4CwRMWtLSvyPJ7BNEIRSaBaEKcy5YpTwvVU8PZfjLkFwO1Kx1Sun\nMUnk5fVYiYeNrSrlXAdiuRWEVq3KGFJnF2qYTp3yytFLRfLaQ9iFW0FoQekNQqws8zXHiuQF14Jg\ntaWyavcmHnwojcZ9P1wLgsWWqvswy2rF8nTrMJVSW86VOKdrQbBIt0BzTaVtwdDdobxiqWGajSBY\nSMaScxDGnD8Hmw9MS7ySzOGThdyoxWwEYWwy1kyWWmsl1qRvM5gc51ufM4dPcXtzjEtB0DRjbu6J\n1t0MZuotRMs3Lbvs9VLOLgVB0+Cg1ku2a6Nvg1kr7LJ3yhZumnApCGNonZytr1+TOXTRrS5lH4Kw\nBUsFadHWOYgCoEvoU+IdgrAFTQW5D6u2WrLbAynxDkHAPFoqrUTsdeFCEKYm1ZgHQpHIeZjLrYMV\nXAhCX1doSpJZmMDjqRJ5FYXUhqXE5jpjMSkIKYFpXWFL09K/Eolp9an8LlIaqiGDq2qUuUlB8F7Z\ntVMy/l57C2u0565JQejiOYFqoSmGmkRBix21aC4IOQK+TXXnVphT0NZyDRWFWiNNvedUc0Eoee+o\nLclT2RYL7clYYj3F1HOWKOu+a1vNqVSaC8Ia74EewrZYaI+RxgU/rF67FWoEoQXaW9wpePLNky/a\nmbUgeG4BPPmm6SGjd2YtCIEd+kQhRCI/JgUhVyJEQtUlxxDz7opJGrGcVyYFIVcieBvtp51cy50B\nuuOnVahSMCkIHpibGJVYT3HoOTWLiBZCEGaI1Rasb4n6IZV8iN9zFQ8TgqB5OvJcE2coJRdHLfGw\nceqiq1YxIQhDpiPXLiSrra0FhoxSbPGw0WPZmxCEIXgspLmSexPaYD+mBMFjFw3QfUtkCY8xqu2T\nKUGo0QrUmFjUt0jGPqIF3I+HEY1jciPlPKmoFoSWexGmHs95jaHkeO1mvQIB019HaqL1GBvVghCt\n4m6Gxqfv+xY2jN3Hth2zNdqqHdWCEOjDkkjX7C14EZ8QhOBDeEluYLootFygpQWzEQRPSV4aL8m9\nZooo1IpFDdFKYTaC4C3Jg+HDlku+KZrKlPzMmduzEYTAHzkeqta6thVCEJwTt0oXs46HlrhosWON\nCUHQFrSaTB3FaP29fG7W8dDSwmuxY40JQcgZtFKVY9d5ve8zaYm1GIRI9mNCEIayr8Uswa7ztq6w\nra+viXUsrC+wUsoWl4Iw5wqgKWk30WhXyohGrWs4lrLFpSCkYHV3pH3UngyTiuZu+i7bNIlADfYK\nAslHSB6SfGXj2CmS50m+tPq5feNvD5E8S/I1krdONbDWnn37jnunht+aY6tZsGqS0kN4FMBtPce/\nKSKfX/38BABI3gjgbgA3ArgDwMPktDB7mHyTwhi7rfqqlRCFBEEQkRcAvNvzp77QnQTwhIj8TkTO\nATgL4OZJFg5gyOukzYLXkARjWs+5iGVN5i4KU54hPEDyZZLfI3np6tjVAN7c+M5bq2NVGLvkluau\nbE7m4udU5hynsYLwMIDrROQmABcAfGN1vE9bs4U3lhrTjdX9Hkqdt3Q8Spz/kjH/JCLvbPz6XQBP\nrz6fB3Dtxt+uAfD2tvOcPn366PNiscBisdhz3RTbLv69exuhaZSaN0qMAtx3vtLLu2uZdDTl/AcH\nBzg4OEj6LiXhrCQ/BeBpEfns6verROTC6vM/APhTEfkbkn8E4DEAf4blrcK/Arheei5Csu+wCUJU\nypMzxvvOtevvNcu61rVIQkR6+xd7ewgkHwewAHA5yTcAnALwRZI3AfgAwDkA9wOAiLxK8kkArwJ4\nH8BXa9f6HEFt1SoFx9ScmZhS1jUqq4a8SuohFLlwox5CtO75Gfp2R0v8h9qiyfYp7OohzG6kooYC\n9fbgs9SGKjVGTw79vqayK2GLa0EoNQNxakFoECULaIyTJlEoER/XglDqOYDGRB2CloSuTS6/NYlC\nblwLghVig9pptFgZ2asohCAowFsFrU2r+HkstxCEIJiAt15CCIJS9iWat0S0irdbhxAEpUwZTKOZ\n2pWnxvU8iUIIwkiGTJ/Wniw17astZLWu52WDWROCoLHCDZk+rb0132efpiTXZEsfJXoLNX02IQit\nK5z2JCxN7viWXpa+NblFoabPJgShNRaS0BJ98dQmujlGo1pcv0O1IGgLVmn6/E09VuLaJeleT5vo\n5rAnRRS0+a1aELQFqzR9/qYeK3HtksylbK35qVoQAtvMrYe3DUtxCEEIimGtdSyFpXEKIQjBZKwk\ne0usiEIIQpBErQ10LVSasVgQBVWCUGpBk2A6tUf8eUX7iEZVgpC6oElqMLUGfQxrX1r6NOTanmJf\nghy9hRKvpFUJQiqprUjK96wk7tqXli1oqbUTgXrlULO8U8Yg5B61OTU/TApCTjR1Ua2Ik2Vqlnfq\nxkKayl2VIGgKTAs0iVNt5jwfRZMoqBKEOVeIoCzac0uLfaoEIShLvMUZjqZnDmO/O4QQBOMMSYw5\nbE+Xu6LUfuaQKtql7DIpCNGaHeOlIufCejz6RIGst42cSUGwXuhjsTIGIZhGd/CSSL2cNyUI25LS\n8pJVQ7AyBqEVWsttLC3ePpgShG0BShnFmOteO9CLlnLLOQKxtiiYEoQUtiWFlmTZh7dWbo7kzrWa\nomBOEKxU7CGMfXrsaU6HBRtrsO3hYS1RMCcIHhkrcjnndLTGgo01aP1qOAQhCAxRupcQgtCA6B4H\nYyl96xCCMIIca/YH7dEmzKn2lBSFEIQRjFmspUuLeQXaKkBrdglzjTEvU/amKCUKIQgTKfFAcMyb\nhpTkmHPPZGjlqfH6euq5SizHFoJgnNorKdUaLZobz2KYs7cQghAMwvrAL6/kEoUQBOXU2DDUw16R\nFikxVXuWi6zOidR1+UpfQ8M5vVEq7lNEIQQhCJwxRWhCEAJXxK3KkrFxCEHYQiSWTeJWZcnYW4cQ\nhC3MJbFC+PwyRhRCEGbOXIRvrgwVhaaCcHBw0PLyxfHsX9e3MT0Nzb0TT2U3ZESja0EonXD7zr/L\nP82VIYWub2N6Gpp7J54EYU1Kb8H1LUPphJtyfs2VIfDLvrxzLQhBEAyD0qipIhltZBA0QkR6bx6a\nCUIQBPqIW4YgCI4IQQiC4IgQhCAIjghBCILgiBCEIAiO+H80xGvObK8Y2wAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7f32dfff2cf8>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "prec = 1e-10\n",
        "\n",
        "np.random.seed(15)\n",
        "random.seed(15)\n",
        "\n",
        "A = make_random_sparse_matrix(200, 3)\n",
        "print(\"%d non-zeros\" % len(np.where(np.abs(A)>prec)[0]))\n",
        "pt.figure()\n",
        "pt.spy(A, marker=\",\", precision=prec)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next, let's apply the same visualization to the inverse:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "7148 non-zeros\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "<matplotlib.lines.Line2D at 0x7f32dd371f28>"
            ]
          },
          "execution_count": 11,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQQAAAD7CAYAAACMu+pyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztfU+obslx36+N8MI2TBQhjUBSbCMLZRIMY5NMhCFwRUD/\nshgTkFCyiGUjEFjZZJWZ1ZvZyQubYIww2EKMQUKRF7YlCJIilLvQxlZQFDmRLE8WE2kk5lkQyRCv\nZE9ncb/zXk+9+vOrPn2+c777+geX+33nVFdVd1dXV1X3fa/UWjExMTEBAD+2twITExPHwXQIExMT\n9zAdwsTExD1MhzAxMXEP0yFMTEzcw3QIExMT97CLQyilvKuU8hellL8spfyHPXQYjVLKC6WU/1FK\n+e+llD87PXt1KeULpZRvlVI+X0p5ZG89WZRSPlZKuVtK+XrzzOxPKeW3SynPl1K+Vkp5fB+teRj9\nu1NKebGU8tXTz7uad0+f+vfNUso79tF6e5zdIZRSfgzA7wB4J4B/DOBfl1L+4bn12AAvA7iqtf5C\nrfWJ07OnAHyx1vpWAF8C8PRu2uXxcdzMUQu1P6WUdwN4c631LQA+BOB3z6loJ7T+AcBv1Vp/8fTz\nOQAopTwG4H0AHgPwbgAfLaWU86l6PuwRITwB4Pla6/+ptf4IwKcAPLmDHqNR8OB4PgngudPn5wD8\n8lk1WoFa65cB/EA8lv15snn+B6d2fwrgkVLKo+fQsxdG/4CbeZR4EsCnaq1/W2t9AcDzuLHjW4c9\nHMIbAHyn+f7i6dmlowL4fCnlK6WUD56ePVprvQsAtdaXALx2N+3G4HWiP687PZdz+l1c7px++JT2\n/H6TEt2m/rnYwyFoHvg23J/+pVrrPwHwHtwY1T/H7egXg9sypx/FTerzOICXAPzm6flt6V+IPRzC\niwD+QfP9jQC+t4MeQ3HaMVFr/T6AP8ZNSHl3CZ1LKa8H8Ff7aTgEVn9eBPCmhu4i57TW+v16/497\nfg/304Jb0T8GeziErwD4uVLKT5dSfhzA+wF8Zgc9hqGU8hOllJ86ff5JAO8A8Oe46dcHTmS/AuBP\ndlGwHwWv3B3b/nwA9/vzGQD/FgBKKW8D8MMltTg4XtG/k5Nb8K8A/M/T588AeH8p5cdLKT8L4OcA\n/NnZtDwjXnVugbXWvyul/DsAX8CNQ/pYrfWb59ZjMB4F8EellIqbMf1ErfULpZT/BuDTpZRfA/Bt\nAO/dU8kMSimfBHAF4DWllG8DuAPgIwD+UPan1vqfSynvKaX8bwB/A+BX99Gah9G/t5+OTF8G8AJu\nTkxQa/1GKeXTAL4B4EcAfr2JJG4Vyi3t18TERAfmTcWJiYl72Mwh3MbbiBMTtx2bpAyn24h/CeBf\n4KYa+xUA76+1/sVwYRMTE8OwVYRwW28jTkzcamzlEG7rbcSJiVuNrY4dw5tdpyO6iYmJHVBrVf84\na6sIgbqNeOfOHdz4iYpa9R/v3RY/lrwePe7cuaPy0Xhl+G9F29u3rAyWvkf33v4u87K0z8wdqwfT\nNkPTu248bOUQ6NuItd78WAj0Hw5L3ig9rP5m+G9F24usDJa+R/fe/i7zEtniGt2ZthmaLdbNJg6h\n1vp3AJbbiP8LN3866t5GPNdfl28hp+WZ4W/RZp9nabbCIruUcXrs1Z9nn913LIFY/hb6bXYPodb6\nuVrrW2utb6m1fkSjubq6auj9DrJGFi0mxnNmDbrl2X5+9tkrV6+Wtn2fiVIk3xERAdP3du6kbG0n\n7TVetj+jFs99uitzfkZA8vPsw5Idve/SK8optkIppWqyS8kZgUeb4dWDaBF78qN3LSK6xZlGY6Hx\n2mq8l8/asxHw+I6QY/HQxtGjzdhntGlZ/czOIVBQjaLi4RzCzbv8BGciAMmTmdDeSR+FLY18S5nn\nkDHaIbDtzzX3o1GK7RAO+bcMvaEQO4leOyZ0l3RWitHm1IwuGr21I1l0mgxPP/m+/TxioVp1BWY3\n3ALWeGjfe3Sw+hm1YWgzNpeR32J3h+ANfCYcyobfvc5D4xOdHMjFHIWdrEPMVrwtHWUNJLtIvfnT\n6grMnEY1JUnb6sKcFFjR3tLee6/Jlfw9R27pH/XZS5E0GktHD7s7BK8A1eZI7TvPm1sOgPGmLb38\nbKEnitH6lMklZZ+8RdDK0sZG1iFayHaZYufSzos8vJ0xyuGtZ4wj0WyodRStY5A0kS5Sd8+BaOPt\n2bs1D57tZO3zUDUEz7AvPa87ql4PC442/qPrHjnZF1ZD0NB6bgkmB8uEu1qe2ZOzsfoweSf7zANL\nz+a0ss2ovDvDr6Vjwm1Wlyhakp8jWk0fK2qKIoRI/hocKkKYmJjYHhcfIcg8tGdXzj5fSyvps7vz\nUbHleDFteudxbXS0db+99r2RVw8uwiEsYO8iLJ+jgWVCvshAPJlegdLTQ3vnhcWs8TJ98PTw6Bla\nT05vKrjAK4ou7yMdWjoL2XGI3kdpZY+deAjfX2rKwE7gOVHKWH0yfYxkj9CN5TF6HCZ8ZMf7olKG\no4fMR8Icq4nROJxD8M6e27BZXnRp37NgQ3BJa4WNmcs0kbyFnxfqZu5MRGf6XojfRipMejHiYlMP\nLZv/Mzl6xNN7xtaOsmmoZyejcDiHANgDYN0ys26GRdAuCVk8rFt3Gt25kHU+Wj/lYvduXMrPFo0F\nyylFMrNysm3kfLY6yd9ZWVl7lG2k3bXYIkI8pENYIAeC3QGi59G7iJ7ZDbfClmlCpmgbgZ2jPcDc\ngMwWQCP+DDJFyh6bZ7C7Q9A6qXlCJmxiQ1opg7kebUUhXojofWYqxhqyaUlPuGnNQ1S01BxlFOLL\nCzqZcZI8PHnsIrYiFq1v0ZVtZkOSYyFlWfpb0YJ2jTmDiz1l8HlvewWa5ZGl7cU5ZIySLw3/UvSO\n6Ef0ZWu7uu88L+zfQ5iYmNgOF3XsyIAJhXorzmtk9tC29FkZbHo0On8fWZuR7dbonrELhi6ijeYt\nqglocrYY2wx2dwgjjbcdXK1qzg68Bi1PtXJF9rabzPu0P3Cxfrc1k0g+e0vP+9zq5VW7I1kWvLqR\nJitjN5p+8rn22br9aM2bJ8uSL4+WvRqRV7vS5GsI3x89ZejNzdrBtY6UtpId8QSOcZMvm8+Poslg\nqzqD5dg8W9m75jEKF5MyaN5L7sxs9dY6suw5F+6Ftet41eReGZndXxuL3nsELE1v5OdVzbVIypIl\n7Ubj20ZccrfObiSeLhqN/PF0lzxGphK7Rwg9Xpdt43n4Nd5eRh+j++DtXlvvUlm9evhkdAHWnwr1\nyLXsZGSEt9YGNT0injftDhwhMOfVWVg7h3WXIIsod7R0ab+zRi7bZndgJqpq0ZuetbIYWvl5rS4Z\nRxWNiecM2gjP4p95l7FxJtqIEKZ8e0cI+8i+HbngxEQPLqaGsKA3N4pycpkranRWm8zOtgW8WkAP\nD4vXqH5qfDM7ZIY+ysEz+krZbD7fg8heo7Zr+6m+fxgjhImJhxkXFyFIsDt/D9+ROe+RsJWua6IS\npm0UFRylXyN4ZWoC57K9i3AI2oWdCFla7busNnttWH4evadzNnTNjNmI1KOV1crsOaqTYO4+tPI9\nntZceJetGLnee+9ZpHOUvrC6sHjVuubbwCr6WYahnR6wF2i0QV6TyWjt2aMh7Vl2LDy6bNWc6Qv7\nzuqH1CHqV+YI0hs3b04iG5DvrTFi5i26+MScZEUyM7h1NYTeBe21W+skGH6jZSw8gfPIymCU/Ow4\nMnJZ3bJOYyTWb1oHryFkK80eem+TaanBIt+KJDIyvDZeuBeFvT2httVXhqcGNn3K1AKYMDw71+z9\ngag/zIL07tdEvBkdGV49OIRDsAZPhn2ZXCojU4Z/Vg7c4wTa51FY54WY2kKILsl4siw6Nl+Xba05\n0ELrRW9tPnsN2VuozBhoNiD5L+8yY5S5xCTHZaGRfbNSL8krotFwCIeQBZOj9/Ky3klDYBc1o0t2\nYcuF1OuoesLO7OJioiTLaWTlswXgNU6HeddT37Get+PBRINrneuhHALbqWhgMru6DNPanUDz1BFf\naRhyQrWwkDE0K0LQdNUQ7RitnMyuFvH0nLeMxqQeGfktr2jD8MZi4cNEGtqc9Mxt+10bD2Z+Nd5M\nvx9od4lFRSaHY3gAcVFqhKzRaI0kk0fvUQRr+Y6aN8ZhZeWM7P8RbabF/CfUJiYm7uHwpwxHwKjq\nbq+s0TIWnlvwXYsj6rRglG5H7qOHQzgEK++SOVV0ypApqEgZbVgr+UW5LlvzsBao1y+vb9b7BW3V\n2qtdRHqxz1o5kTOy5rgHmbmwnrO2w861Zk8abVRzYOsr0RphcQiHYBU/FoNmi4Rt8YWpylqLKAPm\n9IEp7mhHTlEVWhYtIx17+rbAcz5tPcMbi3ZutMKoxlti1JxZz6xx0moTFh85h5kit3Xk2W5azFj1\nHjvuXkMYWWiyfks6+blX3lY4R4Fv+Q48OGbLs9GyMjyzenhza8m2+mzxYpy1p1f0LjNG1vhEPG/a\nHbiG0HPEZT23dh1NlhXSeWG6FX1kU4IoTLTerQ2tvR0t2iE1feTztn+ZIy8rDLd2Sw3asaPnVLw+\ne06ipbXmuneONL7e2HiRh9X3yOHsHiE8+Hy745/RO++WUQbLezTdKIyM/Lbi34tRstfy6Wl/+AhB\ng7WjejssU0xZO/gav54doSdiYCIPVp7FN9KjJ79ndWR5MnPfo6fFU3sX5fC9fc5GxAxtVofDOAQv\nP2NSAXaXXBNyS17sc/nM6p82Bl7KoPVbMwhvfLy6gfadGb9srh2lgNozzRZYh9/TD+15lIp4PNbK\nYZGOII6WMrySZnxYuFWoeeTwFjhGiHsEZJxVT8qZHaMM/bhUZd5UnJiYOMFzCKv+xaRSygsA/hrA\nywB+VGt9opTyagD/CcBPA3gBwPtqrX+9Rs6Dcu9/PlJhaXSRyKp0e20tmhERDBsqb1kYPlfbXuzR\n93Ze1spfW0N4GcBVrfUXaq1PnJ49BeCLtda3AvgSgKcZRj25aXRMw8IqxLB85YRY79lcsf2eyYvb\n39b7NfB4sHm5hlF1nZG8RvBZy4NJTyTd2nle6xCKwuNJAM+dPj8H4JcZRkxHsvcLGGgelSlYWe+0\nhZFdrNbiyharWGSMKHIKWX4LfaawGPHqbdsrd03Emil0e8XpUVHJWodQAXy+lPKVUsoHT88erbXe\nBYBa60sAXssyk4PjHTF5PG5k688lNLq1pxFR5ML0i7mj0OoaRTAWDfssAsNf6sHOJSPXG4tM9Bn9\n1uDZGhMVenPpzW12/JhxWPuvLv9SrfWlUsprAXyhlPIt3DgJCs8888y9z1dXV6j16hXvo2MwDVYe\nxea32QiD1YWNNli+7W+2T+y7nl2OiX68HLdXj56qv1ZvYX6zdRgmUmJtYI3dLO+vr69xfX0NAGiW\nnIphpwyllDsA/h+AD+KmrnC3lPJ6AP+11vqYQn+IY8dM4cZqFxV12EKfZqheoZHRndFtQU/O2iNz\noWEcGttHdk6z/c3aSrYNOw4sb07HDW4qllJ+opTyU6fPPwngHQD+HMBnAHzgRPYrAP6E5/nK77Lz\nmZTBes54WSaEl55b2xXkwvB00XYYS2425/d2rOU9U5/wZLdjFsmURjvCGUj5Gl2rW+QwejeOdq4Z\nx6Q5AytlyDq5rtSv/58xKz8L4I9wkyK8CsAnaq0fKaX8fQCfBvAmAN8G8N5a6w+V9oe4hxBNjBYF\nyHeSX2Tg2R1tq2Hq4Z1pk41g1mIUr96ocQ0vhj7rIG0+82LSxMTECYf/J9S80GZtca+n0qrRMukK\ny1PjM7qf2TZtmDuCz1r08lnTjrGVLP9eexvFM4vdHcLozvUYZ2sMnhFo+RnrzKTBeYsx6kNLF8ln\nF72spvciKoqxjjWrR89iisbWK8D26BM5HKvP1rxtsZHu7hCigpYs0mg0Ef+Wt1dYk/TtM624F+mj\nFR6touPyIxe7VmhiF0srTys8McXOFlE/IycXFfQyMq1xstpofbPmJ9NPS56Wx7P1FEuu9X3EHC7Y\n3SEs8HZhuVi9gT1ncS5TXIvaScOWxtlW0jOLyxtXyTfile1vVCCLvkfz7EU61nNm5/ZoPd5am4iH\njJh6bDtj1xHtYRyCtXMDcYTQLhTLGD3DtEJ7NuS3DEIuRm231L57pxeto2DCb23Baycn2gKzFo+3\nk1m7tbXLeeOqtZPjLmUxizhahB4PK+oagXauNF2YjcWiZbG7Q7AWYPs+Cpm03ZL1ml6YqC0c+S7S\nbcHSD8thRMZspT6RE9HGz4oamL5oUYmVjniOUkuXZPsoApDjl5lzCRl9abK8DUfy9sbd00M+l1Fi\nyycTIbHY3SF4O2M2B2tp2NAtyi0ZudYEacbB9MGb8JbGC1U9vp6cyDExRqgtrh7dLKegtYvG19t4\nrAUn3zEpmsZ3AWNHTKTMRAa9juHi7iGs7bDGL2O0W8ttF0Jm5xup2yg5PTpFbc7Vzww8nXrslU27\nsrrcpzn4PQTA3q3kb3a3lp+X71mvrdHK3caTaelhPbd2RRmyymdMZBTxZXbjaLFa0UE0F1Ya40Vd\nFk1PtGT1w3rn6Sw/R87C4iejE0uPDMI5vrQIYWJiYh0uIkJY0OP1zsHL4s9EBFZbjyebO7N97N1R\n1iAjL6vbOeaWfd6rSxSBsHJHjvPhHIJEtHA8WDn6KD2sSnp3Qae8snjlpSPWiYWHKHRl9Ms6O5Z3\ndsx6ahPa5yx/r+gXyWVgFS7ZZz1rosWhHYKXa2faZ+iiGgRDZ/GXtRAtD/Z2iOwJisVPG9dIfhZW\nPmzpxegYyWKiq6xNMbzks6y9Mn3trZVkcTiHIM9dW48pd0hZ5GONof0s6ZmjroVO08mjlWfs8qhT\n7g7eKUTbJtKzjTg0feQzaxFHWMafOWKVcqwjv+gsX+PrPbeOD61iZraYJ/mwc2TxZ6LONUehD7w/\nSlFRm3y77fqwfC2kEWsLrleXll82JfH6p+kctYl0smgiPVmcY649x7PnsWmGFzM/99tdQFFxlPFE\n37f0f9Zu0Bvaecd2TBvvHbtQWvpo/JjIIIOeCEXTJ4KVn6+tuUSy1ub7vbQeDuEQvPxNe+eFhV4I\nz0yGlu9r9MziyuSlVpiqpU9eDYDVQfLWagje+EUyZFE0qhloPJlw25PJyvJoMvanvWPSC2uMLB5s\npNjjuA6TMozn3x8unkN2FO4x4SKwTx/ZtGQUz5HYwi681A/IzeMWY/sgn/lPqE1MTJxwETWEBd5J\ngBcKLiGSRWvJsXhGsjUdPBmart537bmkZ/sraS39PbpIhie75S2fZXlbbbUxsXh7NNa49qRpUdiu\nyfLmyLIZT49s+rT2P2oZCqmsLFJ5Oa32zgvBowIYUyBjjZgJhNbknmuQLTSOlNMjzyoAAvZ4SFvw\neEheTJusbXk82JqFRbugN704VITAVunZ3V/yZnXQePQYoqdLtBAljceHcQyjquZM9NDufJ6jaz9n\njTcqvEb0UZ8tO4wWpxXtRRjlJNux0NZNxP9QDgHQDcva/b0dof3d8s1MEDOxrQEyjsgL92V/14ak\nDG3PuGSeyzFki6YRPGcSpUYMT7aN1EFzSMwcZeeB4W+tG5fXUYqKaw1kT7Th3pgqcG6HytJIugxf\nj1aGvdIh91T3z2EX3ngD5zuhWUOb0+GWnTKMXBi3EVv3nXEKvTxv+7wdoX8XdcrAoKdIp2FEUe6I\nYOsKa/j3vGN4RjWkS5+zvZ1BhIt0CKMQFe2Ydz2wctyRcrw6yx6Las2itlKbLJ9MzYXhxby/NAd2\nCIfgFc96Jt17HhWOeos6PQ6Frdgz7Rj5Xl/XGK7FNzMmvfMwol2m76NpR9o4K9PlfYk1hImJiX5c\nbA1hi6OYUTuRfL7mSK83upA/Ed1aZI/uNLnWzjwiQmH6mhmPniim5c2mKOw49OicxeEcArtgo/w4\nexHJ4sNOkFfIYy9bRYvJ4skW+UakHxG88WxrAGur7ZqOHj9pF9ZdgXbctXP8SF57D4E9CWMuuFnt\no3bZDfBwDgF4cKHJ53Igl2fMbqQNvmdczPl05ISiHNu6lSiNpTdS8oy/fZbZrSxZntwWzOWrSB9t\nN/bmy3K4meiOaZeNbNdGwEPvSTwsNYTIQNfuWD26yB0T0Bdvr16a49rrrJ/pa4ZHD5+W3opW2LFh\n7YmhW3Q6F27dxaSHFYyxjjLotfS92MNZnROZ/m01FhddVGRy3JGFGCvk08LTLLQwuJUT9ZXJNb17\n9FrNQ8rO1EGyYNKR3nw64su0G11P0dqzKY1FO0ovC4dzCNqgRAPlXcJh8npZoNPqBzLPZwzXqhXI\nPkl5UX7I1kFYPSMea9Bj7L0yFr49/WjnxXP4TI0gWzvovSC3xZweziF4C2CBtfAZOskzM3kyB5Wf\npSFFRbVogrWoJIoYPMegOdI2z9VoLH20d5KuDXm9IrCmvzU/jNNnC5Ya5DhYY6a9W3j3Oju5Wchx\n0Tak9r1sY+no4XAOQYIxAG9yIwOOQmwrjNQmJ1pQnq5RP63oQqOz+FoOx3NqUo7VPyua8vh6BUJr\nfhiHFTmV5Z2lW+Rgo42G3WS0MWqfy43H04GJhBnMouLExEOGiy0qLtgr/2UwovjVkxv2yI52u73G\ncC2y47fGnpjx2mr+RrfXcGiHYIVLGp1XvdW+j5rwnrN02c4Lk7Mha48uy+feglzE25IXtd9Klyid\nae1JaxvZJVskzNS8JNhaRdZZHsoh9Hp2awK8/JPJRddmNFaxbI+d2MqjrRObtcicwkQ8GFlWUS3L\n33KMazaBtdD6xjqDTA0GOJhDsHbKyMu1xaSeIxw50NGJRc8OwBT+mCPGto8Zw4xovYhF6snAi4SW\nz9EpTEZWe1qS5SVPDrS5ZiJQi7cFpjDobUxRFNbjuGZRccLFiEhp4li49UXFSy2GXQLOkeJk6iQj\naip79GcNr3Pa90U4hN4z2NHY2/GMPG3J5tpWAW4t74W/bMvUILT2WXmsjDX8s7JkenjOCC10CKWU\nj5VS7pZSvt48e3Up5QullG+VUj5fSnmkeffbpZTnSylfK6U8nlWIzVs1w7RqC16btWAWSnZX03hq\nxa62duLVISSNpGNPUywDtU4r1u7mI3d+WZhjx0x7tqYPbP3J0kvjbdXDGH0kmAjh4wDeKZ49BeCL\ntda3AvgSgKdvhJV3A3hzrfUtAD4E4HcJ/q8AUxW1Cn/a90whK6tb+4y9xec983hqxa12HCJ+Wf20\n5wxdNB7Z4mamSs6OqadfJDtrQz36Zk5XJH1mfDWEDqHW+mUAPxCPnwTw3Onzc6fvy/M/OLX7UwCP\nlFIejWRkwEQG0XOPV9vOiy60zxl+1ntLJhuB9ERA2TDW6iM7Dt5zJupj9GRoM8ikYkyk4tF60YzH\nK9KT6UNvDeF1tda7AFBrfQnA607P3wDgOw3dd0/PhoHddZbQNcurfe553p5dMXP8Jb/Ld1b431ME\nHJGjLvp5ea833pbOa+oDa/tljXEkP1P/8Oik3XjpB2NjTB9G/+/PWtdMNZhFq4Fpkw27RtNmMHoM\nttBTM86elIDlH9GdA2tlReMzKsUciV6HcLeU8mit9W4p5fUA/ur0/EUAb2ro3gjgexaTO3eewTPP\n3Hy+urrC1dVVpzoT50YbDUwcG9fX17i+vqZoqYtJpZSfAfDZWuvPn77/BoD/W2v9jVLKUwD+Xq31\nqVLKewB8uNb6L0spbwPwH2utbzN41lqrW+m9ocsbXVvplm0jftb7XuOX7az+ZvhvRdvLK6sPwDmU\nnn6ybeQ8aOkBO3eePtEzpp1GI3XN8Fn1byqWUj4J4ArAawDcBXAHwB8D+EPcRAPfBvDeWusPT/S/\nA+BdAP4GwK/WWr9q8E3fVLQmvf0+YgF4x3OenK318vRdw3uEg2H7qtGNGidv/NfinNGQHA9gnCO5\nGe8L+EdWGQ8/wsP2Go3VTttpos/L97Zd1IcRhjHKWUb8NV2jdpLee67RLFibyliFOzYCyTg3ufFY\nNtLqYcnynj3Y7gKuLrNFJZYPcwwTwWvb6msVirTPrZFr4aoHy1gkjXVUN8oZRLLlmESyM4vAassU\n61p4fbD09xZsRrZFK8dW9s07ZeiRr+EwDoFZFO3ALD8RbcYwpS7WjmV9tnTyUpD2mbaQZc5t8Wl1\nyEZAzNhHO+Siw0KTccQ9NYgW1lwz9uEhsh1rHNr3ng5av9s2PdHGmg0QOFDKMDExcR4cOmWQno3d\nrZjn3o6b9ahWW2Y3tKIGSxcv+mB5sXpY8jT6qK+ZcWGiKY3G6z87p8x7Zqwj3pnx8n57OnryI/3U\n90eOELQwdVQ+zCATsmX5td9lbu6lBxYfja59Z72P2ntyI9pI5hZzyKZI5zwxAPyiIJPmeTpn7fTQ\nEYJE5MFG5MMLnedhtUWalSc9vTfpy3eZD0d5Y2QAURHMMr6IV4SlLzJqsHhl5i161mNDlp49kP20\nxk0rGnrOwOqnVVDt6cuhHALrudcsTo2PZqjegmQNMFrUckI1BxKFj14l2jMwD727Z4YvG9rK/kQG\n7xU8vU3AGkdrrj395QK3oDmhHuc+woEvOJRDiKqyCzxPGXlj2c4Kv7XvlkfOOigZDUg+msdnIyNG\nN+tUY3mWdaySr/XOixAsnWR/ot3Xcvztj6ebHEfrFICp+LMLstXLc3JW31knwuDQNYR+3utz11E1\nAyCXD/bKiQx0hDyLT2YBZPPdtYgipMy4rbWrkf315sIa3/upx8FrCL07kkXf4zGjfL83r/R2Yym/\np6+tDItnb74uMWInYuoyaxClEG10FqWDVnSgwRpjxvmw6VP7nYmErSjMwyEcwtrdnOG1pt1anrIf\nkZOI8kPLaVnyIz6tM+pxTF46psn3CmSevqwu2QXWPoucZ8bxZcE6GyZd6cUhHIIFdndtsSYHXmSy\nxSpWF9abZ42tXcSeDMsh9dQqNFhj7jlBtqbAyrfasu+YOWLtKpOStDRsnYh91+Nkb2UNYWJiwsbh\nawgSPWFPNtSNwnKWJgNGzpY4l5wMRuu05Tw9DDikQ4gKNFaRhckfNRls0c1yEGzOvwYjeETFR0a2\nNzbW/ET8A62mAAATHUlEQVQ6ZcYsqoUwPBidIt4W1hSg1/AbRXNIhyARZRZRtdiDl/t6ea332Tqx\niByPVeSz2lrGIhewR6vx9t5bC9g76pI6ye/emEW6avUTj0c0FpJOzp3mHGVfWtrI0WrFV4+fbN/S\nRbiYUwYWmQgA6Ks2t88to414s0Ui61irNQLLGKyjM4/OKpB5RcmWxjPwlkdL1/Yja8CWnBZaYTYq\nRjMFXGu8mEJsS8tEn9bmIjcIpkC5NjqZRcWJw4FxUBP9uLii4oKodiDp2OdRG0Zeryf2Qk42P+6N\nkrSQU9NhVC5uyYr08iJBZoyy+keRD/PMep/J7S1ZPTr0YneHEOW2Xo7Y0lnvvHaWTCbkt+4rWPLa\nMFLLT2XI3cqRPDO7Z9QH706C/J5JsazURZOh6RjVKjQdsm003azUiK09tPLkXFmOKpojCWYDsRC9\n390hRKcDlrFEiHaLNSGptYC1gtDyXFsAXo1C0zez22TA5PaeEUqjH7WjMX3WahUMmPqCdNo9fJk5\njnhpDrrXfsNa0awhTEw8XLjYGoKFUfnTVnnYEXCb+tZbP5nI4yIdwqjA4jYHKLepb22ofZv6dURc\npEPowcO4s1iV/T3GYq1MLf/Onvb0nEZleV06LsIhSKP2CjfWZ/amnzyKi3iPOpLTPls0Ee3yzqrs\nM5dl5Hv5w9C3vNmjRO/EQrvclDmi9C55ZeeRuSQUvWMKsN6Yb2GPhysqZo/VfBnbhphrwlhPN5av\nXCiMkd6WkDsav7X9tHhoz3vlreXVHrXmTkEuqKg40mDZ+wnMe4v/muMfa7dg+UZ3JiTP9vm5MTqS\nYu6ByLbskah2pNzK1XRhdJbQjq177GnUvRTggA7Bw1bn8ED/ws6Eh167LRbsuaOBNXqPyO8ZWjbq\nsuYrE2llUqs1fEbicA5Byz8zOa41UVbePGJSGP16w3lvDJjx0XZSy/lYNZOeRenVAjRaq8bTIltD\nkG2jXV9+tmik/uzzSK+IlzeWTHsGuzsE2UktDI52z6WNFoa3n722Uid2B2fCe00v2U9rwTI35Rj5\njKHLMWRv6Xnhc8vDe2fpZbWRz9gQPnovnZPmRDJRQk84by1ixnFF8sPN62hFxYmJiW1xa24qevmd\nLBpFqYBHZ4W72XCaiVBYHoxOLA9Wh6juId9p6R4Dbw7YtlJ+Lyx+zPjI5z01BMu2rXesTix2dwjs\nwgX8MFI7c88EIFY+Gw0wG5K14WDPArDSiwhRnu610cbZMkI5B1rtgqkrRLI0ukwe7tFoqaDXH0+n\nZUy8vrLpZjRv2hrpdQy7OwTZYekRGcNgdodoocrfWi6pvdMMyZItjYyJTCL9o/eec9MMNlqsXr4v\n23p1FzkmEtEuq+3ibJHR2v01fbOLbeGr6ZFZ2B6YmklkYxYusoag7Xb9enC8pNF4O+5WQzqy3718\nLVrpJLPRGdt2hK5bYvRYsnxbJxQXge0awkU6hEvCHkZ5ZKxxHGvl3QY5I3BriooWrDC3J49iCkIs\n3y139Iwesl3Ec21hyuPPHp316qClcL3plvWOTaUi3qysnhQyK3vBoRwCa7ASzBl7Jj/zeEbvWPnM\nJFvGoBUY1xi95DnCuNeCTeOYQl22EBvpEYfkMa0211rhco1t92xGh3IIzCJkqrI92NrYNSNl2rD9\ntXJRr2CWxd4h8Zo5WuMUsrqwc6t9txxahtcazBrCxMRDhoutIXhHTqPajto1MnmoDPPlMRjDlx0D\nje+a/LwXsrYTjclW+nhjLfXwxo39zKRy7PycI2U7tENYoOWJ8n1EYyGT4zF0zHstzPfO4jXI+oEl\nR7svYY1f1MfsUaKWC8tTBk3niK/G06LJgg3d5V2K6LmGZSxkOie/Z2zPK6oz43IRDiHCmkIYW5Dx\n7hyszSelHj0LMOIb1RC8MYwu5Fi6RzWhzP0PqUv0jEHWNhhdIoeX1YGZF+tZT93q0A6BLSRKtDvn\nFgUZyTdySJYzke3ZyEG+i4wmMozMaUWEbBQRtRtdZmId0YJMGrfwt2RFJyKRvY4oDEc4tENgYXlu\nbQDPXcf0vHg2pOvZcbIG7bVZs2it8Dqrg/zeW0uR6VSka6SvJ7tt3zMfDO9ROIRD8ELRNYU2bwdq\neTNFod4ClvaMDePlc2nATA6uGX0mpLX4Wjpa7+Ri6BlPyUPb7b10yHNoIxdXNoXM7vzM5tbbn9Ah\nlFI+Vkq5W0r5evPsTinlxVLKV08/72rePV1Keb6U8s1Syjv61LoBuyNZqYU30UyI3hPGe3r20HvP\nMobE5PQMXRbe2LLtmOdZGq0NE8Vk09CMLpmx3yIVBrgI4eMA3qk8/61a6y+efj4HAKWUxwC8D8Bj\nAN4N4KOl8L6KCWWz4S67G8qKvRUpZJAtFGZlRDm4jHBkZHGOEFYbA5Z/Zrdko7lo97Z4ZMaktaGe\nvrJR6ogUUyJ0CLXWLwP4gSZPefYkgE/VWv+21voCgOcBPBHLeOVvWRRkik/LwHm7pjfoWmRgLWLW\nWUieXvgc7fa9ITYbdcgxz8j0+GnfRxXHNJ0jRI7TSm16ozF2A/DsTfJiU5yeVGhNDeHDpZSvlVJ+\nv5TyyOnZGwB8p6H57ulZCtoAtL+tNl4RSL4f+bnHGLXPXj/ZZ55sJv3oDXettlEqwqQx3nu5+HpS\nHs02vDR0q7A+Mz6sHtmx6HUIHwXw5lrr4wBeAvCbp+eaP6JUskJ1+Vvbka2w0StieeFlJjXx5ERh\noNUvq48eDwbZkH5NiDoivPX6p43d2pRB4xeNiXzG6sHI8vhZY7M2jXiV/1pHrfX7zdffA/DZ0+cX\nAbypefdGAN+z+DzzzDPNtyvUevUATRsORoUmzeMu7dswOArpMp53kSH5tjQyDNdC8mi3iFKnNr2y\n3kV9ZSMyLcXR+qfJyYb3UTqmhfIM/yg8l++tsfV49O7g1ljJz9rYaPyur69xfX0NALhzB3j2WVsf\n6o+bSik/A+CztdafP31/fa31pdPnfw/gn9Za/00p5R8B+ASAf4abVOG/AHiL9ldMI/+4yZoobfH1\nGCYrr5fOo/eesfTWO5k3t0Y2aGpc2d5cMHpoObg2LpF9aLYheXp9sfgytBFfT/dW5wy8P24KI4RS\nyicBXAF4TSnl2wDuAHh7KeVxAC8DeAHAh26Uq98opXwawDcA/AjAr2dWPTPQCxij93bdHjntTjti\nsTPvR+bdzHh4u7KnZ2SgTL7c8o/0aHWNoqdMZOk5iczOL9tkaBmHJvvMOBNGn93//DkypJ6dPDu4\nWblrPXnPrhPpKmWM4s3wyeyIWVpgG8fLRiBsSkIttsAWM7Zp6cHIv6g/f/YKJVpxRdKxzsDjAdgh\nbCRTYtnttL5EPCP5rWzGGXhja73PwJMlw/Dlc69sOY5sP9pdXpvj9scbU4tvRBchmiOW1tIx0ukQ\nDiEKczOhWl9O1ce3Z+e1Fm8mpI54MXw0vozBeDKiefR2NtYBazKjVDBaUJpj1Rxbz1x735k2I/kw\n2D1lOApGhdVr5ANjdfCinXP3dYTM3hC7pQe2TZ8uAReVMgB2imCFQ9mwUWsb7SxWCJzZUb2dqt2d\nNLmZUHKBLM5FoTqTpkSw2jD5dQRrp/TSB0nfGwVpvDKhvEUnx9+bgxHzE+o0I4SJiYcLFxchTExc\nCrbYpffE4RzCqAHOhG2jwfBkCl1M+9tikL2nDFvJYenWpg5eihPxG5nS3Xt/G1OGLQpYe/PZW8Yo\nyKPJSxqb6Cjy6HNw/yTpQlMG1vPK7+y5uwftJpjFR55fe3w0HmwhKbObMHI0GSyvSA5TXIzuhGR1\n0/SwjjOZCMyyKw3R/QtPz16wxecFjMM6tENYwBiaNiDRLTOLX0urfZa0S/XauoRktbXO7707D2su\n82i4v2usdxja3YBsldw7EfHaRAtYzgXjKLQTmkgXOddr78VEJydSzurI+Agpw8hQrDd0O2c46Bnk\nJd5DONccnYP3JaSK6/t48JTB2/W0d2uKdhYdEw5a7a2w3tod5Y4sJ5jZjTK7RrTrZuGlaEwaxN4m\nZPoo07XongFjQ5rc3kiBkeNFwNqYMqlWLw7hEFp4DsAKizRjY6/CenTeZ2txybB+oW8NhA1ZLTms\nc4j6IheTx6ttZzkU6yKUJdcDswvKeffSuui7xqul9ULyrDP3wNiT5Gl97sEhUoaJiYnz4fApwyVg\ni/CMkZNNfY6EI+p0bjDREMtnJJ2F6RBInCuYYU4cmHZHwBF1OjdG/f3Euexgd4fAFF5Gy9sba/PM\noyGj+1a05+C/Zrcf9UdVix5bYXeHYBVN2vcSa04Z2GJjhmcWUWFqjXzWaHsWg8U/sytlC4XZFKpH\nF3YsIltl2jMFQKvA3PLZCocqKnqVZabq7NGx7UeilbmH/In7ONr4S33W2H62b4cuKkY7wkITHWd5\nZ8ctb+/ITbazjgBlWyYa0W5IMv2RNJpO7DGeNbbRMaGlt8bLoovmJxpzi44dK+aGqnwf9Z8ZT0t/\naY+WfUeyvLZW3zx0/b8MI8EU0SLvZ91d0Hi3g+jReXKjCCTy/lqo2urljYF1Pm5Bo2v1sZyW9Z0x\nOHYOM/PqvcvcGZCbi8XbswuLZ6RzpJOml7z0FdkkI9/r++4RggdrcLI1h5ZXb9jo5Y9RyJaVqTkK\n6733LKKJ5Ghg0zYtQpDvPf5M1GXZB4PoMlGWx8j2vTUmoL8v99rPGsJ28KKG0TL26F8Ge+uXtZ8o\nyhutzxrbz+rm1RAO5RAmJia2x6GLihJsQUVrExXgIjlRG4uGKbJpz6KCWeZZBK9IFumVkWHp6KUQ\nnr69xT/vu6ajNafa84xdRvozemrv2XnKzuUhHIKXL8qw26KTv9m8PhNKyrbLZ7aQJAtaVsFIo7H0\nYYyTycdb2b01D0/v6LtWD2DqNi2NtUhkIU0bb022fMcW7iQv712rs1bws+pFUQ3EmofIQRzCIUSV\n7sgzajsPWxHXdGGNKgNpgNYCYfi2C4EttFqLixmnRS9rR5aOoKVbU1D13muGrjkIrT0z/paNebYa\n6Sj5Wjzl5sM4MKl3pmD5Cv2OUkPIGE+PoW2JducfoZfFh+GfiXgyURLTRysCykDb/faa6y1kH8HO\nL6KGEFVR29/e+56dm5Hn5cCW52YjCm3H1WgY42ANyNuVpF7e7mTR5qreD+qVCXmtyMVqE81FuyN7\ntsfYowQb6i+0bb+YyNgDQ3MYhwCs7xTjVDK6RKlMD1+tTWR4Uj5TN8hC4x+FvBGflr7HOY4y8gwd\n4M+75rx64KUkTN2AecdsNBKHcghR3tObF1mInAubzzPvotyTWUiSdlQ4GRkOs0trubsVWUT8JT8P\nWj1F6infMXNm6ZgZc23TkM52Tc7f6hRFIywO5RAAffA0A9PaeLuLNjnWhGkDrOW2Gb3lb4ueST9a\nHlG4LGm1dz3OyJNlGSozZhH/aP7bflq0THqRiU48W/MiTc2Bana2ODRr/Lw5zOIwRcWJiYnz4CKK\nihMTEzFG1o007O4QmHD3XHp48kbp4YXuGR5b0Paip2C7Bd/eNks7NvXq1WNEkZQpQLOyNOzuELzi\nkVdZtei13+37KO9jijtrFpnV3wz/zAlHVPeIxqwXmXyXac/Serm21U4WHz07YXXzalYW2LFhCt7R\n8bCF3R2C9My9RssWHi35rK6e7pLGKyJaz7S2VhvveaQ74J8KaG16ojnL8Vp0siDoORZpG5Z+3umD\npPHmkj1SZZ0rawfMRqd9ZnR4QKcjFRW1iYoWOKv+wicbdWTg8dYWG3tvIlpETL9G9LsnYmsX2shx\n9+wi83xrsPp49mm9055zdnUhRUUtnPYGMxMJsMd53jPPE7Pyre8aveyjttuy5+Os84l4MP32dmCt\nDTOuHm+po8ZrhDPsofP01PTTjmotRxHdQ9B0C9OWI0UIExMT/WCjoIs9dmRyJu85kwOztQUtX9U+\ne3oxuWEkx+pD1A+GxtKbkeP1w6PxZGbmy9KXocvQsGPA0Gft22vfRhIefRiBzQhhYuJ2Ia4jXGiE\nMDExkceaffYQDiETgm0Nq8DVG9ZpbbKhuMWnV4/ROIIOC3p16S0ajgabkq3h4+EQDoH5iy4WayfM\nOr7JnhJIXVp+Gq8er87kjT3I8sucYGztzJiTJ01udNkrOt0YBeaIvffyEjPeh3AIFjIeccGoCWON\n3NMtOnbrRa9ziop8y+8tjF5bUKPkRIafGS/LwZ6j3GXZUu+tQ+to09VhFhUnJm43pJM/bFHx+vp6\nE75rQ1H2GC7mc32W/PqcOfwiS85dTw7eEwGOACNP698eta6eI0j5LPO3L7s7hJ48hwmVNLpMqK8F\nL9rdf89Q7ty57gr32PN39vzZ4tmj19IfuWCYG6CSLsqJe+8csPbh4e1vv36gzTLOmuzeGlD0zBub\nKM3TUqVobC7iP3uN2vTwXpOtZIo7EQ/2XfRszdXlNXqNoB2hy+jaxJ07vE698tbwsvrL2IHnFA5d\nVJyYmDgvdi0q7iJ4YmLCLCru5hAmJiaOh5kyTExM3MN0CBMTE/cwHcLExMQ9TIcwMTFxD9MhTExM\n3MP/B7z4nHdwwEh8AAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7f32dd4602e8>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "Ainv = la.inv(A)\n",
        "print(\"%d non-zeros\" % len(np.where(np.abs(Ainv) > prec)[0]))\n",
        "pt.spy(Ainv, marker=\",\", precision=prec)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "And the Cholesky factorization:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "1819 non-zeros\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "<matplotlib.lines.Line2D at 0x7f32dd3511d0>"
            ]
          },
          "execution_count": 12,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQQAAAD7CAYAAACMu+pyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGWpJREFUeJztnU/oZtV5x79PkCySgLWiI6hNghFrS8CE1kqg8IaC/7qY\nUFBsF40pghC76arOamaZLBJKCBJIRCagWLNIolCiFftbuEks1pqiMZPFVEeZX4RqoFmZztPF+74/\n79w5/+8595x77vcDw/ze++f8ufee7z3nOc9zrqgqCCEEAD5SuwCEkHagIBBCjqAgEEKOoCAQQo6g\nIBBCjqAgEEKOqCIIInKniPxCRH4pIv9Yowy5EZGzIvKfIvIfIvKz3bYrROQ5EXlDRJ4VkctrlzMU\nEXlURA5F5NXBNmt9RORbInJGRF4RkVvqlDocS/1Oisg5EXl59+/Owb4Tu/q9LiK31yl1eWYXBBH5\nCIBvA7gDwB8D+GsR+cO5y1GACwA2qvo5Vb11t+1hAM+r6k0AXgBwolrp4nkM23s0xFgfEbkLwA2q\neiOABwF8Z86CJmKqHwB8U1U/v/v3EwAQkZsB3AvgZgB3AXhERGS+os5HjR7CrQDOqOp/q+oHAJ4E\ncLxCOXIjuPR6Hgdwevf3aQBfmrVEE1DVFwG8N9o8rs/xwfbv7877KYDLReTYHOVMxVI/YHsfxxwH\n8KSq/k5VzwI4g+1z3B01BOFaAG8Nfp/bbVs6CuBZEXlJRB7YbTumqocAoKrnAVxVrXR5uHpUn6t3\n28f39G0s954+tBv2fG8wJOqpfk5qCIJJgXvwn/6Cqv4JgLuxfaj+HH3UK4Re7ukj2A59bgFwHsA3\ndtt7qZ+XGoJwDsAfDH5fB+CdCuXIyu6NCVV9F8CPsO1SHu67ziJyDYBf1ythFmz1OQfg+sFxi7yn\nqvqufhjc8118OCzoon4h1BCElwB8RkQ+KSIfBXAfgKcrlCMbIvIxEfnE7u+PA7gdwM+xrdf9u8O+\nDODHVQqYjuDit+OwPvfjw/o8DeBvAUBEbgPw/n5o0TgX1W8ncnv+CsB/7f5+GsB9IvJREfk0gM8A\n+NlspZyRy+bOUFX/T0T+HsBz2ArSo6r6+tzlyMwxAD8UEcX2mj6uqs+JyL8DeEpE/g7AmwDuqVnI\nGETkCQAbAFeKyJsATgL4GoAfjOujqv8iIneLyK8A/BbAV+qUOhxL/b64mzK9AOAstjMmUNXXROQp\nAK8B+ADAVwc9ia6QTutFCEmAnoqEkCOKCUKP3oiE9E6RIcPOG/GXAP4CW2vsSwDuU9VfZM+MEJKN\nUj2EXr0RCemaUoLQqzciIV1TatrR69m1m6IjhFRAVY3BWaV6CEHeiCdPnsTWhqFQ7e/fvn7jfz3U\n11a3Xv71XD8XpQQhyhvRU8buWFt9yXIoMmTQBG9EETYUcil8LualmOuybheXuMl1zGazGRzf380f\n1q/nupWk1jUrXb/x89DK81HNdVlE1JR3KxemNHPXcy3XlfgREejMRsVk9j2FWKYsaFVjMay5GyfF\ngITQnCAAHz68MQ11ygNfurH0ufpeu7R0vVsqSwhNCsKe1N5Ca4zHiiXq1MN1ykWqwJe4hkvrmTUt\nCEA/orBHtcxDsrQHLyexz4ft+FzXcMnPa/OCAPQnCiQvsQ25tHiGpN/q87wIQQDCRaHVC03Kk/Pe\nl36OQkVp7ud5MYIAhF3ENXed107Oe9/KczR3ORYlCAB7AGR9lDJEm1icIOS2J8SkRTEqfw186fd0\nD2x1GW8vZYg2sThBAPKKQsyFnrv71uLDX9sg10pXPge2utSs4yIFAeDMAyElWKwgAGkejUuip7eh\nC9P9W/I9XdJsx5hFC8Keqb2FJT98JZnrupiEL0YMW7t/S57t6EIQgGmisIQ38RoCsFLJXc7WBGZO\nuhEEoD27wrAsU8u1lMZZi1xTc/sw8XFatZ6rOaccgc4EAWhLFIaNmA16S6l7kzI1ZyrLPo1xWrXu\n35xTjkCHggCw8bVMS/empbK0QpeCALTTS+gR17VNve68X23QrSC0NA7sDdebNXUq2JRmq/crV7lq\n1M+XZ7eCAFwqCi10EVt9yHOT41rXvl+l102oUT9fnl0LAtCe81Lth7wWta9/Sv6t3CvOMhSgpdmH\nNWAK0JmaxhRyNu65hwxzln01ggBQFKYy16K3OdOwkbLYzv43hwwdQVFIp5UudA5SFtvpqf42VicI\nAEWBXEqMV2nPz84qBQFYh9qTcGK8Snt+dlYrCEDfSt8TOWNCiJtVC0Ls0KHUw7ikh7zGEmq9vpFt\ngUs1n4dVCwIQt7x7qQez1Qd+jsYZOz05h0fjnOtATF0LIjerFwQgTBRabbQuWg25zu09mrucS7zX\nuaAg7GjNo9FFDYeWnNT2L8h5Xm9QEEa4egutPDQtNvRWlrBLvTa+82oukDInFAQDNlHIsQBHr0xx\nTW5R4MbUXCAllZTnj4JgYSwKKUtZtf6g1xas1q/P0km5vhQEB0NRmHspqzngh2fIGAqCh95EoCa8\nlu1DQQiAbzYSy1KfGQpCAPuhw1QX2pTZi5IPVgtBPC06FeVIY6m9IQpCIEu9wVOI8QrMNe3YgkjN\nlW+I6/Lc9RWt9KSLiNbKeyol3Zhr03PdUjBdjxau0ZQyiAhU1Sg17CEkkOq8tIQhw9wOOjmGYuP0\npmJyrQ51t55ryFLMrZw9hHRaeFMQEsNWbNhDKEJvKy8ttS5LLfdUSqwkTUGYSAlRqDWEWGpvp6Vy\nzylOJepNQchAblEI+TISmUaphrv0+0NByMTSH4S1wftl5rIpJ4vIWQC/AXABwAeqequIXAHgnwF8\nEsBZAPeq6m8mlnMR5DYyLtVoudRyk+k9hAsANqr6OVW9dbftYQDPq+pNAF4AcGJiHovBN3SInXYs\nbbTM6W053BcjBilRpL70aqaxdAPnVEEQQxrHAZze/X0awJcm5rEoXI3Y1lBSbAa13Gtz2zdyR5HW\nXpJt6T2jqYKgAJ4VkZdE5IHdtmOqeggAqnoewFUT81gcJmeW3G+Okr2Hud5yuXsGOXsHts+4zc2U\neqWcN8mGAOALqnpeRK4C8JyIvIGtSARx6tSpo783mw02m83E4rTFvtGWWkthiatAD+0LrfUMhum0\n8hm3HL2Vg4MDHBwcBJ2TzVNRRE4C+F8AD2BrVzgUkWsA/Juq3mw4fvGeiqHQyDYvIde71D0Zp9vi\nvS8SyyAiHxORT+z+/jiA2wH8HMDTAO7fHfZlAD9OzaMXSnXvOWS4OK0Yo6bPNpM6ZAjtWYSmP/eQ\nIbmHICKfBvBDbIcIlwF4XFW/JiK/D+ApANcDeBPAPar6vuH81fQQ9tjeFq63SMk3TErarjoA6QvR\nDu0uvuChKftDy5Qa4Vgq/5y4eggMbiJkZTD8uSFyz7sP053jnJ5pZWahJhSEmSnpV1CLlspeuiwt\n1LVkGSgIFTAZGU1j1hBCx+4pi6GEUmPkF+v8FUKIQbCFUe6Uzw76zqEgVMIkCkMLd+iDl/u4pdBb\nfVLI7WkKUBCqMlZ6k1PMUmIZ9raR2DRjl1DzxVKUimWYss2VT0wacwxXKAgNkBL/MGauh8g1bRrq\nkTks0/gc3/kl14rYN1DXUGFc9illsF0vX8wLbQgrYCgKJZbGynVOSDq+t17J7v5UG8IUQTIRci9j\n73fJ60dBaIhh7EOut0CNsXbo23QNTPGaHDLX7AYFoTGGomAiZYw+N6au95zutzVJ8WUIOcblHZoT\neioSsjLoqbhAYv0Qch2Xk9ZmSEoxR1mmzmiEnkNBaJRQO0Ko5XnOzlhqoFMMLXUu5yhLLicp+iEs\nmJKrIgH53zBj20Ho+DlXHUPTmuLr4Dsu1oYQ67eQ4xjn+bQhLIPSIbGt5DmFXOWdq961ri9tCB0Q\n2lvI2aNIfVhzliEmraXFZqTEn5SGgrAgQkQhNUgqJ6WHOi0TM2TwOaLVEAwKwsKIbWwpb7scD1qt\nt3Vtf4dY1+XYc1LzCoWCsEBKv4FbCmeOpUbjcpV9STYYgIKwWJb2oPlYcn2WXPYxFIQFk3vKKuXY\nKWnlzmdKcNgwHVMauXtkoelPqVfStDKnHZdN6NTV0qYQfbRUn5xliUlrGAgXFybPacduifVobIkp\nb91SY/4Ucl7bqUbJqR6rFIQOcImCq7uZ0gWN8c7zlSm3SE31egw5f7h/br+QOYYzHDJ0RqmudEpX\ndu58W84jhXL3kkOG1VBqSjJntzQ131K0KAZAnXJREDqkNz+FtXo91oCC0CkhdoVaxObf6hu8RygI\nHTNXTIHJyJXL8JV77j0kfDtn+iUMjzn8Laxp06hIyLqgUXHl1B4itEbORVnmZI4yUxBWQOzQIdb2\n4PJDKOW6nJrnHtsqxilz/TnrmOLnYfONSPE94ZBhRdScb58z76l51bpOsfmmlpNDBgKg3KfAQoyK\nruNj0vbt8zWSEKNiyPmhdYgxKsbFI5SBPYSVMtfbiLQHewjkEsZ2hVDvQ9tb2TeetTHuTZSeNgxJ\n21WX2DF6Sk/JZxMoZZcB2ENYPeM3fw89gR7qkAvTtWAPgVgZ9xRyLOxZe0qvphjUrvs4/9hrQUEg\nUZ+iz/U1Yx8pDWvKkCGXd+WUVa9ziMnUa3/Z9CKQHpi6AnBuSnymLPbcHNciZcGTmrCHQI6o3d1N\nZanlnkqJelMQyBEtfWAlphw5Pk4T451ZahgQe2zoWpoxUBDIReQWhdQpspQPtNhsCCmNcuwsFGN4\njcln+Nt07W1TuqUW1uW0I7Eyx/RdjjyGjYkfew3Jl9OOJIE5hhA5DXe9fey1BhQE4qQlu8JSWPL1\noiAQL7nCp1OYo3HlziM1uGrKsbnwCoKIPCoihyLy6mDbFSLynIi8ISLPisjlg33fEpEzIvKKiNxS\nquBkXlJ6ClNjC/b52tJ25etiXC5bnEYORyWX56DPABo6e5IrFgQI6yE8BuCO0baHATyvqjcBeAHA\niW1mcheAG1T1RgAPAvhOQPpkIYSOeYdj+lLjZFe6vjxt5Rpv8/0OIbacKfUabg+puwuvIKjqiwDe\nG20+DuD07u/Tu9/77d/fnfdTAJeLyDFfHmQ5uN5S46m/2tGOU9yYc21Lic70RTva0vIdE0KqDeFq\nVT0EAFU9D+Dq3fZrAbw1OO7t3TbSCbahw/CtG2P1358X+iD73qAxb8uYPFK3xZTXdE5ofXLNtOQ2\nKppua4OTK2QKrgY8tRG2OBW3JlKDmw5F5JiqHorINQB+vdt+DsD1g+OuA/COLZFTp04d/b3ZbLDZ\nbBKLQ+YmpyMQKcvBwQEODg6Cjg3yVBSRTwF4RlU/u/v9dQD/o6pfF5GHAfyeqj4sIncDeEhV/1JE\nbgPwT6p6myVNeip2QqwwxPjhp7gw+6b9fDMXqft9+dnyLnX97OfbPRW9giAiTwDYALgSwCGAkwB+\nBOAH2PYG3gRwj6q+vzv+2wDuBPBbAF9R1Zct6VIQOiLkwc9xbsz+nD2YEr0hm8Ck1CHkuA9dvCcI\nQikoCP2RGlPge2uXEpfUMpU4byox+TKWgcxCql1haKR0OfKY9rv2TV0wJWUKNDYUOyY6M2d0qA0K\nAsmKrSH5GnKKR6KPqf4EPq/C0Dxt+32+A6ECF+v74SwfhwyErAsOGcjs5IoDyB3gkyO9mGFLSe/C\n2J5YyH4KAinCfugQMpXncnKKNSj6RCTVNjBOY3icKyApJj3XsaZ62YKyQvMyQUEgxfAZ2Pbj9FRr\nvul3SIMMmbYMedOHeleGBCW5cPkv+M6LPYeCQIpimgKMcfBx9R7GPRBXw3N1z21CNU4jpFwpDMuW\nM9oxNpoSoCCQmRg2ypiurqshDkUgdBbD1kjG6fp6CKFv99AAL1u+U8qQIloUBDIbMVGNrjR820Mj\nDk1jf5twxDbWVFxv9RT/jlgoCGRWTD2FVKZa2U29j9DzQx2oTH4HrjKFkHvm5aK06YdAalDLxXdu\nWqwn/RBIc+QYPgxxvYltv109DNM5KUOGHGIQ468Qcx2M57OHQMi6YA+BNEsOb8ThG7S0J2LscaV6\nQaWgIJCq+IxxMQ5CY6Y0oFTX6xhRCjEyxqY1dSaEgkCqY7InuKb/fGmFHutKw+SbYMPlG+FzNHK5\nJPvyNZ2T4tg0hIJAmiD0LR8zrThXdz0mTDpmGjKltzG1h0CjImmOFqfqeoJGRbIopk5Jppwb8rYO\nSX/8Vo8pS4yjVWjPKfZaUBBIk9jG1zZCA4Rc+bm2x5THZgNIiY+wRTnabByhHpQ2KAikWaY27NA3\nrm1bakCTLe3Q5eV89oCUUOhQaEMgZGXQhkAWTamZA3IpFATSPDkjJIkbCgJZBKHLmpFpUBDIYvB5\nD+4NgTHRgcNzXedMMUbazrEdaxsihaQ7OS6ERkWyRFzTi3M7NuXMz5aWacg0nH6MWbAVoFGRdIbN\nL8D1hg99I8dGMboacExaLjFwrR8ZI0ZcZJV0S8qqwqFpTu16l+yhlEybgkAWjW+J9eFxpijG/T7T\n8aH5jomNkvRFKs45/KEgkMWTc/lyk1vwfvt42BHbgF0GQt9Qx3VMaB4hadCoSMjKoKciWQW2SMMY\n/wVTT8B2TkxcgqlMtre/zcgZO+2YYgehIJBuCOlwhnTNbWmlRF8O/x7OIpimD23l80Us2oSJH2oh\nq8fm5rxvzKaG55vjH6YdslbDuDG6loMzRWeOzzXZQHwCFjstu+cy925Clsd4TUPfuoyuNQxdMxCp\n4dBjURqnNxYjU34m0Qgtu3PoQaMi6Zm5vBbn9o6cAo2KZLWELsc2NsrZ/nfl40o3FJehMiatGHvH\nEA4ZSPf4wqfHxj7T2Lzk2z/UkclkXIwZLoTAHgJZBaEuySmzCD1BGwIhK4M2BEJ2hDgCxZ4/9bhU\nO0OJXgoFgawKUzCUb95/j8vWMP7bt/rzML+Q4YwrhiLEy9GX/h4KAlkdLkNhaAMyOQsNt4cYCk0u\n0jamLkkfCmcZyCoJsc6nWPxjrP2xMwM2d2dbmVN8I9hDIKvF1lW3BRSZXKFN+8ZpuX4P03IFQpnK\n4rMl2KZQXbCHQFZP6KIptgZpOy4kKGlsQxiea+vFuFyax2nH4u0hiMijInIoIq8Otp0UkXMi8vLu\n352DfSdE5IyIvC4it8cXiZD5MXW7bY1qygpHIXETrpgF098+o2UMIUOGxwDcYdj+TVX9/O7fT7YF\nk5sB3AvgZgB3AXhEpFcXDtIb4ze1rVseYwy0MUzDlZ5t+DL87TNaxuAVBFV9EcB7pvwM244DeFJV\nf6eqZwGcAXBrXJEIqYfrrW8yGo4ZiokrbVdIdGwZXTaEWKYYFR8SkVdE5Hsicvlu27UA3hoc8/Zu\nGyGLYdjwbQ3bJg7j7SHDh9D0QoYYIftcpArCIwBuUNVbAJwH8I3ddpNW0T+ZLI5xN374v8u6Pz7O\nZXj0DRlcZfDlb9vvI2mWQVXfHfz8LoBndn+fA3D9YN91AN6xpXPq1KmjvzebDTabTUpxCMnO2MJv\ns/j7zg/dP8VfwTdEOTg4wMHBgbO8R2mFBBiJyKcAPKOqn939vkZVz+/+/gcAf6qqfyMifwTgcQB/\nhu1Q4V8B3GiKYmJwE1kCQ+Oib/ZhfI7pWNs0oklwXILhMyi6jnEFN3l7CCLyBIANgCtF5E0AJwF8\nUURuAXABwFkAD24z1tdE5CkArwH4AMBX2erJkvHFDdj2h8wADPMwiYgpXmJ4rK8s4/xDejgMfyYk\nEFcjNzVWX9felJ6r0Y6FwNQLsR1z8T6GPxMyGZur83Cfy73ZlqbvGFPjtjV6U5li8qIgEBJBjB9C\nqI+Byzbhc4/2lTO7YxIh5GJcb1rT2N+Xjm82wiQsMb0C1+8xFARCErA1epPbsw3bON/kJj083pVu\nzLHG82lUJGRd0KhISCF6C92jIBAyAdfQIEQsTEOEkOOGx8aIEmcZCCmMTRRCIxdDgqHGx43dqofb\nTX+Pz7FBQSAkA8NpvtAAo5DgJVtaJk/FcdqmeAn2EAiZkRA/hdBzXL4EoX4I43PYQyBkZlxTkqbh\nwXD60XbM+Hjf2gg270YfFARCCuBr3LbhQOwxvvxjoSAQUgiXJ6KtS28LebY5K41/T3VM4jLshBQk\nZuGT0GNc7simmYcY2EMgpDChLsypx4TOaoRAQSCkML7xfGpPwbRvajQABYGQGQg18sV4N4b2GmKg\nIBAyEyZnonHDHR7jC7E2GR1jZiJMUBAImRHTuon73+NjfIucuM5JhYJASAWGU5JjN+MQ12XftKMN\nui4T0iipqxu5HJ5i8xxDQSCkIjbnI9N+GyExCqFQEAipjCkqcU+snWCqMNBTkZAGCF001bQv50qE\n7CEQ0giuRVvH+13bXWnSqEjIQjA5L5mWXHfZDFzhzyFQEAhpiPH0o8kr0fVFJpcQhHgvUhAIaYzx\nlGIuW4HrAy97aFQkpFFMrs6hn2izLbjigz0EQhomZG3EVP8FExQEQhon1aMx9ruOAIcMhCwCVxSk\nadbBZIQMGTJQEAhZCDbjYugSaiHCwCEDIQti6hJpPigIhCyIkAVUfPtcUBAIWRghH5i1DQ/oukxI\nh7g+6GLbFgKNioQsFN8iKSm+COwhELJwbEFRY/ixV0JWgulbkqZjfFAQCOkEU6TkkBC7Am0IhHSE\nq6dAxyRCVsiUMGkKAiEdQsckQsgRppkHrodAyIpxLe9ug4JASMeEfnV6T1VBODg4qJl9cXquX891\nA/qqX8x6CBSEgvRcv57rBvRZv5DeAocMhKwIui4TQoIRzflhuJiMRepkTAiBqhoHD9UEgRDSHhwy\nEEKOoCAQQo6gIBBCjqAgEEKOoCAQQo74fwCCU6dFyKJgAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7f32dd3c4668>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "L = la.cholesky(A)\n",
        "print(\"%d non-zeros\" % len(np.where(np.abs(L) > prec)[0]))\n",
        "pt.spy(L, marker=\",\", precision=prec)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Cholesky is often less bad, but in principle affected the same way."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Reducing the fill-in"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Define the *degree* of a row as the number of non-zeros in it."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "2\n",
            "4\n",
            "3\n"
          ]
        }
      ],
      "source": [
        "def degree(mat, row):\n",
        "    return len(np.where(mat[row])[0])\n",
        "\n",
        "print(degree(A, 3))\n",
        "print(degree(A, 4))\n",
        "print(degree(A, 5))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Then find an ordering so that all the low degrees come first.\n",
        "\n",
        "The [Cuthill-McKee algorithm](https://en.wikipedia.org/wiki/Cuthill%E2%80%93McKee_algorithm) is a greedy algorithm to find such an ordering:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 29,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def argmin2(iterable, return_value=False):\n",
        "    it = iter(iterable)\n",
        "    try:\n",
        "        current_argmin, current_min = next(it)\n",
        "    except StopIteration:\n",
        "        raise ValueError(\"argmin of empty iterable\")\n",
        "\n",
        "    for arg, item in it:\n",
        "        if item < current_min:\n",
        "            current_argmin = arg\n",
        "            current_min = item\n",
        "\n",
        "    if return_value:\n",
        "        return current_argmin, current_min\n",
        "    else:\n",
        "        return current_argmin\n",
        "\n",
        "def argmin(iterable):\n",
        "    return argmin2(enumerate(iterable))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 30,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def cuthill_mckee(mat):\n",
        "    \"\"\"Return a Cuthill-McKee ordering for the given matrix.\n",
        "\n",
        "    See (for example)\n",
        "    Y. Saad, Iterative Methods for Sparse Linear System,\n",
        "    2nd edition, p. 76.\n",
        "    \"\"\"\n",
        "\n",
        "    # this list is called \"old_numbers\" because it maps a\n",
        "    # \"new number to its \"old number\"\n",
        "    old_numbers = []\n",
        "    visited_nodes = set()\n",
        "    levelset = []\n",
        "\n",
        "    all_nodes = set(range(len(mat)))\n",
        "\n",
        "    while len(old_numbers) < len(mat):\n",
        "        if not levelset:\n",
        "            unvisited = list(all_nodes - visited_nodes)\n",
        "\n",
        "            if not unvisited:\n",
        "                break\n",
        "\n",
        "            start_node = unvisited[\n",
        "                    argmin(degree(mat, node) for node in unvisited)]\n",
        "            visited_nodes.add(start_node)\n",
        "            old_numbers.append(start_node)\n",
        "            levelset = [start_node]\n",
        "\n",
        "        next_levelset = set()\n",
        "        levelset.sort(key=lambda row: degree(mat, row))\n",
        "        #print(levelset)\n",
        "        \n",
        "        for node in levelset:\n",
        "            row = mat[node]\n",
        "            neighbors, = np.where(row)\n",
        "            \n",
        "            for neighbor in neighbors:\n",
        "                if neighbor in visited_nodes:\n",
        "                    continue\n",
        "\n",
        "                visited_nodes.add(neighbor)\n",
        "                next_levelset.add(neighbor)\n",
        "                old_numbers.append(neighbor)\n",
        "\n",
        "        levelset = list(next_levelset)\n",
        "\n",
        "    return np.array(old_numbers, dtype=np.intp)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 26,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "cmk = cuthill_mckee(A)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Someone (empirically) observed that the *reverse* of the Cuthill-McKee ordering often does better than forward Cuthill-McKee.\n",
        "\n",
        "So construct a permutation matrix corresponding to that:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 31,
      "metadata": {
        "collapsed": true
      },
      "outputs": [],
      "source": [
        "P = np.eye(len(A))[cmk[::-1]]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "And then reorder both rows and columns according to that--a similarity transform:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 32,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<matplotlib.lines.Line2D at 0x7f32dd2976a0>"
            ]
          },
          "execution_count": 32,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQQAAAD7CAYAAACMu+pyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFSpJREFUeJztnU3IbdV5x///EDJIAtaKXkFtEoxYWwImtFYChRMKfnVw\nQ0GxHTSmCELsJKN6R/cOzSChhCCBROQGFGsGSRSKWrHvwElisdYUjbkZ3OpV7hshGmhGpj4dnH3e\nu+92n3P2x/p41tr/H7y8593vPns9a+/1/Nez1l4fNDMIIQQAfCS3AUIIP0gQhBBHSBCEEEdIEIQQ\nR0gQhBBHSBCEEEdkEQSSt5H8BclfkvynHDaEhuRZkv9F8j9J/qw5dinJZ0m+TvIZkpfktnMoJB8m\neUjyldaxrfkh+W2SZ0i+TPLGPFYPZ0v+TpI8R/Kl5ue21v9ONPl7jeQteayOT3JBIPkRAN8BcCuA\nPwXwtyT/OLUdEfgAwMrMPm9mNzXHHgDwnJldD+B5ACeyWTeeR7B+Rm1680PydgDXmtl1AO4D8N2U\nhk6kL38A8C0z+0Lz8zQAkLwBwF0AbgBwO4CHSDKdqenIESHcBOCMmf2Pmb0P4HEAxzPYERriw/fz\nOIDTzefTAL6c1KIZmNkLAN7tHO7m53jr+A+a7/0UwCUkj6Wwcypb8gesn2OX4wAeN7Pfm9lZAGew\nLsfVkUMQrgLwZuvvc82x0jEAz5B8keS9zbFjZnYIAGZ2HsDl2awLwxWd/FzRHO8+07dQ7jO9v2n2\nfL/VJKopfzvJIQh9ClzD+OkvmtmfAbgD60L1l6gjX0Oo5Zk+hHXT50YA5wF8szleS/72kkMQzgH4\no9bfVwN4O4MdQWlqTJjZOwB+jHVIebgJnUleCeDX+SwMwrb8nANwTeu8Ip+pmb1jFyb3fA8XmgVV\n5G8IOQThRQCfJfkpkh8DcDeAJzPYEQySHyf5yebzJwDcAuDnWOfrnua0rwD4SRYDp0NcXDu283MP\nLuTnSQB/DwAkbwbw3qZp4ZyL8teI3Ia/AfDfzecnAdxN8mMkPwPgswB+lszKhHw0dYJm9n8k/xHA\ns1gL0sNm9lpqOwJzDMCPSBrW9/RRM3uW5H8AeILkPwB4A8CdOY0cA8nHAKwAXEbyDQAnATwI4Ifd\n/JjZv5K8g+SvAPwOwFfzWD2cLfn7UvPK9AMAZ7F+YwIze5XkEwBeBfA+gK+1IomqYKX5EkJMQCMV\nhRBHRBOEGkcjClE7UZoMzWjEXwL4K6x7Y18EcLeZ/SJ4YkKIYMSKEGodjShE1cQShFpHIwpRNbFe\nO+4d2dW8ohNCZMDMeidnxYoQBo1GPHnyJNZ9GAaz4T9jz8/1s8lfjT815632/O0iliCMGo24x8bZ\n5wshhhGlyWAjRiOSawff/O77nxAiDdHGIZjZ02Z2vZldZ2YP9p2zWq2OHH4jChdfY3ca3peoWK1W\nuU2IRs15A+rP3zayDV1edyra3qhg4/SKHoQIA0lY4k7FQfQ5dDdSMNt+3ly8RxhCpMbNXIauCHSP\nxWCXqEgsxBJxIwh9ItDXrxCCIddUc0QsETeCsI0xojD0PDm7EP24E4Qh/QpjvjsUNRGEcCgI2xgq\nClMdW1GDEAUJAjDMaeXYQkynKEEA8oT2ak6IpeBCEMY43K6mQyzHHRt1SEBEqbgQhL6RifvO7ztv\n6PiFnOMbhPCMC0HYMGY48q5IYd81PDisogjhEVeCMHUadJ9zzXG4FM7aNzdDiNy4EoSpzJkpuavp\nEZJ9Du8hahGiCEHoOtM2Jx7bOdn+HRt1TIoSKEIQuk2Dbc4Va+7DXKbY5DVi8Hh/RTiKEIQNQ95G\n5JwQtQ2vzr2PVM0p4YeiBKHNroIZo9Au0RGWmOelU6wgtAn9lkGMQ/e6HqoQhDkzJIUQF6hCELYR\nUhRizrQsnalNi6XeL89ULQhAuOHMIWZaygGEd6oXhA37ooUUHWjqpLsYNev8sRhBAFQAp6CJYMvC\nvSCEnupcoyjUlh+RD/eCsGtUYptQayp4gfSxaKxq8GURazv45NS2Yax3+0SduIwQUtXeu9KZsvR7\naLu9RzGiPlwKwpTaceoEoqmLrPSdp1pdlI5LQZjCVGf03J8ggbmA12dUG8UIwtwCsW8MggrcBTze\nC4ljGooRhLkFYug6ixqinBbdS18UIwi7CFmohkQLNdZWfbtvp6DGe1kyrgUh9nv41IusjCF1+ikd\nM/e9FdtxLQixC+m+RVZyFlzVnCIHrgUhN6lFYddK0DUhsfOLBGEPOdrTchiRi+IEYembvcayxVMe\nRT6KE4QhtWfo4cRTmw5D9pOYYksMxrx2FfVSnCAMIcZw4imi0E07lC0xhGbDUBtjzuEQ+ahSEOYw\nZFWl3A4QS2imkiL93Pd8KUgQOgxdOzFmAS2h8Kfemi636C0FCcJEYorCkMVaS5tqLYcug0ULQoha\nK0dtXvPOVCVERzWzaEEI4QRtUZhTmOd2WI65jmen8yJMS2XRghCKjSjMKcwh34akSqsPz2Ij9iNB\nCIRqtjW6D2Uza5FVkmcB/BbABwDeN7ObSF4K4F8AfArAWQB3mdlvZ9o5w8bd4XWIAry5TqjrlcLS\n8rsE5kYIHwBYmdnnzeym5tgDAJ4zs+sBPA/gxMw0ZpFi2/j2HIQck6FyITGoj7mCwJ5rHAdwuvl8\nGsCXZ6bx4UQdt1NTisKcjsWpeL73Yj5zBcEAPEPyRZL3NseOmdkhAJjZeQCXz0zjw4lmqJnGbgQz\n9jshKW3jFomMH+Zu1PJFMztP8nIAz5J8HWuRGMSpU6eOPq9WK6xWq5nmxGOKI3juV/Bklxc7auXg\n4AAHBweDzqUFehokTwL4XwD3Yt2vcEjySgD/bmY39JxvodK++Lr+ClhOmzzeD5EXkjCz3rhscpOB\n5MdJfrL5/AkAtwD4OYAnAdzTnPYVAD+ZmsY4e9a/PRb+WP0KCrVFaCZHCCQ/A+BHWDcRPgrgUTN7\nkOQfAngCwDUA3gBwp5m91/P9KBHCdnvniUWImtZbbe3Nni7e7SuVXRFCsCbDWFILghBiTZQmg5jG\nUhZSHYvuhw8kCIkZ8+Yht5PkTl+kp2hBKLHAtoc576MtGrVMs/aQlthO0YJQYiGaOsxZMxRFCooW\nhNLJPaKxa8dcYqzkJNJShSDkLoQxV17KnbdddG0zi7eytEhDFYKQO5wOvfJS6GvHYopt++6n5wFm\nS6AKQYhJ6o41zzVjio1mxtxvz/eqVCQIe0gdys8VhZhOolq7fiQIe9jmBF6nGC/JaZeU11RUJQhj\nasehbdlczE2/9JWbSl89ulSqEoQpuwltK1QeOirnFPjSBxXlXj16FzULUVWCMIbNK7IchaovzW1v\nGGoufMIfixWE3PS9w+/Dy+AlsQwkCJkYG5nUFC2UnI/a12iQIAzASwGuRRS6DlVankqzdwwShAFs\nC9tzzUBMna7XnaFzjNfI1e+UCgnCCIZ2BqawI2W6Xmv0mh0zF3OXYV80OQvkUtMORQ15iIEihMCk\nHOrspaYuEd27fiQIgUk51LnG/oRU7HoeteRxCosUhKkP3GNByd2fAPSn7/Feif0sUhCm1tZe2525\nX0f23Zccbw5CXcPrc07BIgUhNKGdccr1xoxo9Fx7L9kZPSBBCEDo13Jzpz/vS78Gp/MsaiUjQYhA\nCIebO9OxJocpbWm5kpEgOGVuga9JFIbeC63HOB8JQgJyLU9ekyiINGikYgI0qjANS8prLBQhLICQ\nUUKOiGNXmoqAwiJBWABTmw5eOvNU86dDgrAQpojCmO+opq4DCUIAPDnDLlumLMc2tHZeUi3u6XmH\nRoIQAE8TZYauVlxLoW7nN9WCKTWLnwQhMqEKbGhyi0KMtGt21FRIEBLircDmFAVv90KskSBkxEPE\nkDtS8IKX5fFyI0HoYehiqrVMs5UoiA0ShB66jppjw9fU1JSXUCzxnkgQxBGKEi5mifdDgjCQJRSO\n1E2HUDtwx7JZEcKC2RSqHLtBeyKlKOzryPN4z2uvGKoUhDlLkKUohN4mCHWZs8FsDR2tQ0Z71kqV\ngrBv5GBulc81QSjmXIb2d1LQjSRCPtPanX4XVQrCLmrYm2/OvoRTvpNbQPsIbZfHPOZgcYJQA6kF\nzbMoxLiWx7ymYq8gkHyY5CHJV1rHLiX5LMnXST5D8pLW/75N8gzJl0neGMvwWvFaGLui4NXOEJQe\nQc5hSITwCIBbO8ceAPCcmV0P4HkAJwCA5O0ArjWz6wDcB+C7AW1dBJ4LY9s2b3Z6s6dU9gqCmb0A\n4N3O4eMATjefTzd/b47/oPneTwFcQvJYGFOH4aHm8vQuf+g5IdMLjYdn2sWjTSGY2odwhZkdAoCZ\nnQdwRXP8KgBvts57qzmWDA81RcpXl0PXPwhFjv4Er30YNRK6U7HvsTlw0TDsG7yUkhzC1xYhjwu/\n9KFRjOOYugz7IcljZnZI8koAv26OnwNwTeu8qwG8ve0ip06dOvq8Wq2wWq0mmpOGfYOXyHkF13sh\n6/YhlGCzAA4ODnBwcDDoXNqAJ0ry0wCeMrPPNX9/A8BvzOwbJB8A8Adm9gDJOwDcb2Z/TfJmAP9s\nZjdvuaYNSXsqMXbxSekAY9LK7Zi509+FZ9tyQRJm1hs77RUEko8BWAG4DMAhgJMAfgzgh1hHA28A\nuNPM3mvO/w6A2wD8DsBXzeylLdeNKgjiw8R0Dq+O59WunMwShFhIEC5QS6GtIR815GEfuwRBIxUD\n4aGjMTQp5j4IX0gQAjGnVvFaI5U29yHE+glen0UqJAgiOLmcSpvKzEeCUAEew3SPNnUpwcbUSBAq\nwGONp/6EMpEgiGhIFMpDgjARFfRhbFuOLfWEKy9b23tHgjARFaZxdKOFWIubiHlIEBxRe9ShJoR/\nJAhbyLEY676argZnyiUKsaOIGp4NIEHYisfFWEPbk2shF0UKfqleEDytYeCNlDse9+2X6WFJ/KGU\nYudcpq6HUAwpN2CpgZT3qaRnss/WkvKyi+ojhDap1x8Uu8kZIexKd45NpZefRQlC6vUHYxSO0gtc\nm00/jbdOxlpq+yksShBSs62NPscBaiysqUQhRRqlPx8JwkymrBmwa03GVHakutZQ5mww20ZTm+ch\nQZiJlxF3XuwIkXZNzaLSkCA4IeZY/9KYIwopX6XWiAQhIWPC2Vy1tBfnCRkpTLmXXu5DaiQICSmh\nDevJxtzNhyWKQvUDk0TZeF+OrTYUIVRKrK3ba97sdYkRQRcJAsYXhBIKTqyt25e22WsJzzokEgSM\nd5ilhpMbcuQ/hSj05WsJU9LbSBBEMaQQhaVvTiNBqJDS+wx2EWpE477rbxjSF1NTxChBcMZmrsOc\n2Xil9xkMYapd274zxNmXsOqSBMEZm7kOIWfjdUVmSljskSmisC0vY/LowXFjIUFYAF2R2VX45wjH\nLmKG+N7WviwZDUwSF9ENkbvONtUZYjrRxs6aHTUVEgSxk1KcrBQ7d+EhD2oyVMyuzsmYYfZSBhGF\nTs9D34QEoWJ2dU7GDuGnMmdFqdT9CR5q9NBIECqh3RFY0vLmXebuh+FpObYS3+aoD6EStNz8BdqD\nl2Ldj1rvsyIEMYnu60mPEYnXQVWekSCISXTHNXjdZi6UKCxFWCQIIglj+za8Db+utYnQRYIgkrBv\nSHaKqc2pXxPuGxHqEQmCcEGKGjh0Gl6WzQ+JBEEUS8jZjiHTyHnduUgQRHBSFfYptezYpsOU9RFK\nRoJQIN4HH3mfShxqI5gal1fTwKQC8dr+bNN2Bo87LaeYIVnCc+qiCEFEYchCL7mZuhxbiJo/19uW\nfUgQFsqYZkfuQrqNXIOXPIvcXPYKAsmHSR6SfKV17CTJcyRfan5ua/3vBMkzJF8jeUssw8U8xtTg\nuR1gm3B5G7y0jyHXz32vh0QIjwC4tef4t8zsC83P0wBA8gYAdwG4AcDtAB4ivdYvIjSxOjrbwhWz\nNGnuwwBBMLMXALzb86++W3ccwONm9nszOwvgDICbZlkodhLzbcPYV2xT+gy8TRFeuijM6UO4n+TL\nJL9P8pLm2FUA3myd81ZzTEQiZsddiiXIc4fIfSwpr12mCsJDAK41sxsBnAfwzeZ4n7YWcBuER3KO\ntUiZrqeIZNI4BDN7p/Xn9wA81Xw+B+Ca1v+uBvD2tuucOnXq6PNqtcJqtZpijpjJ0DEDqckZvte0\nkvPBwQEODg4GnUsbkGOSnwbwlJl9rvn7SjM733z+OoA/N7O/I/knAB4F8BdYNxX+DcB11pMIyb7D\nQvTSdc6pzjr2eylEIbXwkISZ9Urt3giB5GMAVgAuI/kGgJMAvkTyRgAfADgL4D4AMLNXST4B4FUA\n7wP4mrxebGOMI3TPm+rUY0tjzOXYNtf05CGDIoQoCStCEDPIEc6HTjNXk2RXhKCRigWRqpPNUyfX\nNuYu9T41zSnb3JW0a7QEwTHdhUxDFaCUu0fHZopz75rSPOS7tXQ29iFBcEzfBq0hCmKujVRi0JeX\nmIOdNmKwLw2vb272IUEQo/DWCdZm44T7ZhLOEbS2MHsSxlBIEEQ17BOq0JvZeFznYS4SBFEMHmvk\noTZ5tL0PCYIoBo9hukeb5iBBEMmZsyakx1B8lyjs69fwhtZUFMnpDkHuHiuRFBvMpkARgsiK57cW\nU+hGC3Pzlro5IkEQYgK7mjwl9ytIEER1dLep73PeEJu/7nvtWKIoSBBEdXS3qe9z3hgzF9u/N2mU\nJgoSBCECsG3Q01zhSd2/IkEQYgRTXpeGiBJSRRoSBCFG0G1+DBGIkpoOGocgRMOU14VjzithjIIi\nBCEaYu9HOSdSSBVlKEIQiyH2dnBD8D6iUYIgFkOIZddCTp32KApqMggxgBhNiaHNgJTzPSQIQvQw\nZ0bmGLy9gZAgCNFD7A7GblohlncLgQRBCAfsWt4tZT+DOhVFsZS038EQQm1XNwcJgiiW0PtU5BSS\n9vLuOe1Qk0EsHg+LtPQt7x5iz4mxSBCEcMLmrUZ78FJq1GQQIjBdR56zw3XqJoQEQYjAjN0abtf3\nUouCBEGIyOxy+iGOnlIUJAhCRKbvVeLYNxua7ShEhczZX1IRghALZMiW87HEQYIghDN2RRGxmw0S\nBCEKojtGIXSkoIFJQogjJAhCFEh7uHXIZoQEQYjCCflKUoIgRAWEEgUJghCVEEIUJAhCVMRcUZAg\nCFEZc15FShCEqJCpUYIEQYgKmdp0kCAIUSlTREGCIETFjBWFrIJwcHCQM/no1Jy/mvMG1JW/MSMa\nJQgRqTl/NecNqDN/Q6IFNRmEWBD7XklKEIQQR9Ay7VBBstANt4QoHzPrbTxkEwQhhD/UZBBCHCFB\nEEIcIUEQQhwhQRBCHCFBEEIc8f9WiLZtL4A3sQAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7f32dd2dfda0>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "A_reordered = P @ A @ P.T\n",
        "\n",
        "pt.spy(A_reordered, marker=\",\", precision=prec)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next, let's try Cholesky again:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 33,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "1188 non-zeros\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "<matplotlib.lines.Line2D at 0x7f32dd26bcf8>"
            ]
          },
          "execution_count": 33,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQQAAAD7CAYAAACMu+pyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFTFJREFUeJztnU+oJcd1xr/PCC9sw0QR0ggkxTayUCbBMDbJRBgC1wT0\nL4sxAQkli1gOAoGVTVaZWc0s7YVNMEYYbCHGIKHIC9sSBGsilLvQxlZQFDlIlseLyWgk5lkQyRCv\n5MzJ4vV909NT3V3dXVVdXf394PLu7a7uqrrv1tenqk6doplBCCEA4CNzF0AIkQ8SBCHEARIEIcQB\nEgQhxAESBCHEARIEIcQBswgCyXtJ/oLkL0n+4xxlCA3J8yT/k+R/kPxZdex6kmdJvkXyBZKH5i6n\nLySfILlH8vXasdb6kPwWyXMkXyN5dJ5S+9NSv1MkL5J8tXrdWzt3sqrfmyTvnqfU8UkuCCQ/AuDb\nAO4B8McA/prkH6YuRwQuA9iY2efM7Fh17ASAF83sTgAvATg5W+mG8yT2/0d1nPUheR+A283sDgCP\nAvhOyoKOxFU/APimmX2+ev0EAEgeAfAggCMA7gPwOEmmK2o65rAQjgE4Z2b/bWYfAngGwPEZyhEa\n4trv8ziAM9X7MwC+lLREEzCzlwG83zjcrM/x2vHvV9f9FMAhkodTlHMsLfUD9v+PTY4DeMbMfmdm\n5wGcw/7vuDjmEIRbALxd+3yxOrZ0DMALJF8h+Uh17LCZ7QGAmV0CcONspQvDTY363FQdb/5P38Fy\n/6ePVd2e79W6RCXVr5M5BMGlwCX4T3/BzP4EwP3Y/1H9Ocqolw+l/E8fx37X5yiASwC+UR0vpX69\nzCEIFwH8Qe3zrQDenaEcQamemDCz9wD8CPsm5d7OdCZ5M4Bfz1fCILTV5yKA22rpFvk/NbP37Mri\nnu/iSregiPr5MIcgvALgMyQ/SfKjAB4C8NwM5QgGyY+R/ET1/uMA7gbwc+zX6+Eq2ZcB/HiWAo6H\nuPrpWK/Pw7hSn+cA/C0AkLwLwAe7rkXmXFW/SuR2/BWA/6rePwfgIZIfJflpAJ8B8LNkpUzIdakz\nNLP/I/n3AM5iX5CeMLM3U5cjMIcB/JCkYf87fcrMzpL8dwDPkvw7ABcAPDBnIYdA8mkAGwA3kLwA\n4BSArwH4QbM+ZvYvJO8n+SsAvwXwlXlK7U9L/b5YTZleBnAe+zMmMLM3SD4L4A0AHwL4as2SKAoW\nWi8hxAjkqSiEOCCaIJTojShE6UTpMlTeiL8E8BfYH419BcBDZvaL4JkJIYIRy0Io1RtRiKKJJQil\neiMKUTSxph17PbuqKTohxAyYmXNxViwLwcsb8dSpU9gfwzCY+b+Gpp/rtatfia+S61Z6/bqIJQiD\nvBF7yjg5vRDCjyhdBhvhjUhe29Bdx4QQ8Yjmumz7wSXu7Eqz2Wxq6a8VgD4xyF0w6vUrjZLrBpRf\nvzZmc10maa68m418F5dG1oMQYSAJSzyoOJqdpVD/7Gr4IcSgzCBYQownO0EArjT22A22S1QkFmKN\nZCkIO3bWgqtxTmmwPteqOyLWSNaCAAxrmHqqCzGN7AUBuHZcYXfMlc73fk0kJkIsRBAAtyi4GNuw\n1UUQYkGCAPg1WjVsIcazKEEA5jHt1Z0Qa2FxgtDVdYjVcIdaHRIQsVQWJwhA/3hCX4Oc079BiJxZ\npCAA3TMPfQ0yhwYrK0LkyGIFAej2aIztuDQV19oMIeZm0YKww9dPoc4uvashxrAg+hp8DlaLEEUI\nAnC1KPiMIfh2L8YQQmRkMYg5KEYQAHdMhbZ0sctRZ0zjztVikFCVTVGCAPh5NI75UU9pCLk27j5S\ndadEPhQnCEAcC2GNDWGNdV47RQoCkP8sQw55hmLJZRdXU6wg+C6Gars2NXoaixwoVhCAaaLQJOZK\ny6UzVszW+n3lTNGCAPiHYwvhJ+Dr+yBErhQvCDv6rIUUJru6BVcT0oITYViNIAD6AY5BC8HWxaoE\nARgvCqmFZEh+scom8VwfqxMEYJibc/2aKcQUoVhP2bY9MUS5RNvKLXdirmXoyi9WeiFCsEoLYceU\nICuudKFNbJnsIjWrFoSu8YQxId1D7gQlMRBzsGpBANLNPIzpMqjbcAUJZBpWLwhAOlFYyo86VTmH\n5CNxTIMEoWLIBrNLWiORM/ILyQ8JQgOfH+kcIdZiU88/ZXQniWReSBAcDHlyhZplcEVZSikSUxqm\nwsOVgwShBV9R8J1lGJO/np4iNRKEDlL3cYd6Ty4VCV2+SBB6SPnjTe09KUQTCYIHOT2xtZBJxESC\n4IFv16GZpu/z2LLEwDXtGntgUyKUHxIET8ZMR/Z9HksModnRHCRtK3OI2RV1jfJDgjCAIc5LKcrR\n9nkqQ+s3Jf8x07siHhKEEcSefZj7xx8ifmSTkIvIRDwkCCOJKQo+wVqXttRaDXoZSBAm0BSFkI2q\nq9HHaFwpG2yIOBQiDhKEiYSyFFzb2Xc1Up/BPt+8UjOmXiINEoQA+O463XePUGVJlZeLucVGTEOC\nEAg92fbR97BsJgVZJXkewG8AXAbwoZkdI3k9gH8G8EkA5wE8aGa/mVjORdC0Evo+910L+Kf3yWNo\nOfruNdYq2l3nun6qpSWmMdVCuAxgY2afM7Nj1bETAF40szsBvATg5MQ8FkNzPGGIv4ArbV96l3k+\ntcswZPXm2IbbtWZDYjAvUwWBjnscB3Cmen8GwJcm5rEoYg0ytuU19tqxaIygbKYKggF4geQrJB+p\njh02sz0AMLNLAG6cmMfiaPNoHDLdNmVH5SmbzvZNo7q6NWNZy3LvJTF1o5YvmNklkjcCOEvyLeyL\nhBenT58+eL/ZbLDZbCYWJy+a/eQU022+XYY24UjRZWher25CXLbbLbbbrVdaWqD/BslTAP4XwCPY\nH1fYI3kzgH8zsyOO9BYq79yZc6AsdN4a9Fs+JGFmTrtsdJeB5MdIfqJ6/3EAdwP4OYDnADxcJfsy\ngB+PzaMUxowr+Lgnz2FqSwzKZrSFQPLTAH6I/S7CdQCeMrOvkfx9AM8CuA3ABQAPmNkHjutXYyHs\nGDLtGCOPZpq2aT/Az9lqapm7yhHi/sJNl4UQrMswlDUKghA5EKXLIMahkXU3+j7yQIKQmCHrHuZu\nJHPnL9IjQZiBMXs+aABRpECCMBNDZx5KX6Eo8ckDCcKM5BqjcSypt58T4ZEgZECbtdDVwIYe70sz\nZrCzec3Q7edClUOEQ4KQCXVR8GlgQ4/3pRnjRjzmmq7VoK57i7RIEDIiROSlmKTYaGaqhSGmIUHI\njKnLp2M2khD7L+Toji2uIEHIkBxiM4bGt3sxpssiwiFByJS+Efu287vjsQbnXPESxizc6ju3hOjR\nJSJByJSYVsKU7dNCL6XuyitXllTWoUgQMqZrPMF33wbfBhy78ffdd8oshwiHBCFzhjovTZ12HMLU\niEddrtlq7PMgQVgIoYK35kBKayQ0OU8Lh0CCsCBiikJKsWkLyrIUllbeIUgQFsZOFFIOyHW5GIfI\nf0qE6bGMvXaoe/bSkCAskL6wZn3MHT696/olPH2XUMaxKISaECtDIdQKJZSDT8xyxLxnri7eS0aC\nsGD6/BTqf0ORaoZg6h6VU65ds1hIEBZO7L0kmwOYvoOauVsma270XUgQCiCEKHTFV2ieax4LuYtz\niLByQ+NVDjlXOhKEQhji0egz9uDDmIVIvtOVOS9uKtm6kCAUho+1MOXp6Ip4NKS/H8qymNNjsGQL\nQoJQIDE9Gtu2XBvCUMugbSBTFkJ4JAiFkrLBDH1iDl2l6Rs/sh6LMiayEMQiiSEKpTwdS6lHaCQI\nhRP6aZbz01GxFKYjQVgBIZ+GczxZQ82KiH4kCCtgbNchl7gFepqnQ4KwEsaIwpBr9KQuAwnCivBx\nXhoSysx3F6Z6Wl+X5ymOSXPEhygFLX9eKSWGAptSpxK/jza0/Flcg88GszF3WfLJZ0wXR0zjurkL\nIObDtZekj0OQ73nfa6dsXCvCIgth5ZQUzXkKIVyyS0CCIHpFoXluZ+6Hio3oGmjMNWJSMxZEaWhQ\nUYiVoUFF4UXJT74xrPH7kCCIA5pdh6XMMsQKIbdGA1aCIK6iLgp9TkkhGkzzPnM5G/mSSzliIUEQ\n11D3aOyL6tzEt8HU4x7U8xk6BRmqgfqI0i5dyaIgQRCt+Lotuxr0kECnbdujDXGbTtlIS+5KSBBE\nJ67GFrO/PcTCmOqq3PZ5iKiVhgRB9NIUhZhh2FMFOXGFlt+RKhRbjvQKAsknSO6RfL127HqSZ0m+\nRfIFkodq575F8hzJ10gejVVwkZahMxApSBlIdi34WAhPArincewEgBfN7E4ALwE4CQAk7wNwu5nd\nAeBRAN8JWFYxM0PWOaQgt/KUQK8gmNnLAN5vHD4O4Ez1/kz1eXf8+9V1PwVwiOThMEUVOZBLCLWu\n+ApTy5jzJjGxGTuGcJOZ7QGAmV0CcFN1/BYAb9fSvVMdE4VQ7zqEaoB1fPdgaA4q1t/7WAsaVHQT\nelDR9RXKmCsMnxH+ub0Mp9y/a1Cx9AHHsfEQ9kgeNrM9kjcD+HV1/CKA22rpbgXwbttNTp8+ffB+\ns9lgs9mMLI5ITZ8o+Hg5Nv8209Tv0+eY1JwlGBrLwXc8Yup05xxst1tst1uvtF6rHUl+CsDzZvbZ\n6vPXAfyPmX2d5AkAv2dmJ0jeD+AxM/tLkncB+Cczu6vlnlrtWAhtDbqrEfeda3YHgHZvwq7uQ1te\nLpfp5jW+1y6NrtWOvYJA8mkAGwA3ANgDcArAjwD8APvWwAUAD5jZB1X6bwO4F8BvAXzFzF5tua8E\noSBSNxIfK2OOci2BSYIQCwlCeZTQv16DgCgegkjClIaUUwSjmMu+c0eCIIIydqHRFFfnodOEfasZ\nYwaXzR0JggjOXA3GN9+SG/RUJAgiCiG9BV3dguaxoYFVuq4vuUvQhwRBRGFqjIKmX4DLb2DK3g5j\n/RBKR4IgolF6dKESkSCIqDQH/EIFVt2lbetO+Nw3VHCVktBWbiIJTeehMS7Pvmljm/wldylkIYhk\n+HQhhloIU+8hrkaCIJLiO64QatQ/tPdkmyNUKSIkQRDJaYZed8UiqAtHX2PuczTqS+Nzvn4/dRmE\nCExbQx0SdGVIIJSpjThGMJgc0eImIVaGFjeJbPH1GOyLijTGU7Hts28MxyHnl4IEQcxKl8egz/Gx\n1/VFRQKGiUwpxq4EQcxO28zDmGCpY9PWRcA17uAbwHXpwiBBEFng8mj0MeNdodtc17lCpnWVoyvP\nruuX3mXQoKLIjilRi9YQ8WgqGlQUi6LZhfBdwtx2bcinuRyThJgBH8eknZOQSxT60nR1IfoWTIXY\njyJXJAgiW3zdnH26CH2zGa4BxDYPyq77L727otWOImtyD8dWGrIQRPa09dm7ZhSa6Ybk47PuYeld\ngzYkCCJ72sYJ6uMDuzEB1yxD2wYufbgWYPXNYixdKDTtKBZDjlOKOZapD007iiJIEaNx6P1Lixsp\nQRCLwidG45TgKi7Px773S7MQupAgiEVSfzL7uDh3CUgzTf1Y33RkV75DycHSkCCIxeJ6Mrcda/oJ\ntA08Nhc3+VobpVgJEgSxaJoNvM/NuctC6Pvc5T2peAhCZMIQN2efezXvWf9cuqeiph2FWBmadhSr\nYKhH4phrp16T8n5jkCCIYmiGd6/jMv+b17rSuu5T94wsDQmCKArfCEpNfIRk6BM8hyf+ULTaURRH\nfcahbXrRdzv4tvcxyMHikIUgisXVhahPD7ativRxdIpBDhaFBEEUzZhpSF/vxBKRIIjiGRp5qS1W\n4+7c2OXUS0CCIFaBr2NR2+fmfYD+BU5LjJsgQRCrwXfQznczlimDgEMFJBUSBLEqfNYy+E5T+ubl\ney4Hq0GCIFZFX9dhl8b1HmhfCdl3n6UgQRCro60xu2ImdFkULu/GKY5MOSDHJLFKmpZCX5/eN57C\nFEemHCwKWQhitbT5GPQ1zBBP/rY85rYqJAhi9bTtB9kVht1njUSRg4oknyC5R/L12rFTJC+SfLV6\n3Vs7d5LkOZJvkrw7VsGFCIkrIEqb+e8KxuJj7vsEZp272+BjITwJ4B7H8W+a2eer108AgOQRAA8C\nOALgPgCPkznonhD99C2Lrh/zWS49d+MeQ68gmNnLAN53nHJ9DccBPGNmvzOz8wDOATg2qYRCJMQn\nVmKzuzAkAnTfubkfn1PGEB4j+RrJ75E8VB27BcDbtTTvVMeEWAx9swuuhVE+swx9XQzfuI8xGSsI\njwO43cyOArgE4BvVcZe+LdBwEmuna9l0W/opebnynINRfghm9l7t43cBPF+9vwjgttq5WwG823af\n06dPH7zfbDbYbDZjiiNEcFyDjDu6BhWnhlaLYSFst1tst1uvtF5Rl0l+CsDzZvbZ6vPNZnapev8P\nAP7UzP6G5B8BeArAn2G/q/CvAO5whVdW1GWxBOoNvOm27OOv0BQSl3D07fkQmq6oy70WAsmnAWwA\n3EDyAoBTAL5I8iiAywDOA3gUAMzsDZLPAngDwIcAvqpWL5aMb6Otn/f9xe/SpgzT1of2ZRDCk5CR\nll3WRqpIztqXQYgA1D0UfQYdXV6OzfvV/+aABEGIAcQeMJx7lkGCIMRAxjZsnynLua0FCYIQI/Bp\nuEPiNO6Y20JQPAQhRuIbL8E3fuLc1gEgC0GISfg+0ed+8vsiQRBiAr57PgDLEAUJghATca12bE5P\nuhyZXHEa5xYNCYIQAcih/x8CDSoKEZAuj8MxA42pvBd3yEIQIjDNmIu+wVPUZRCiUNoiOtdxnZu7\n6yFBECISYxr33BaCxhCEiEjbuEGuzkqyEISIjO+6hq4xhFSWgywEISLjCqySgzXgQhaCEAnockzq\nCuCaekxBgiBEIpqWQt8Gss1zKcRBgiBEQlzTkV1h2PvCv4dGgiDEDLRtGzf3no8SBCFmwiewasrw\n7IAEQYhZadsXci4kCELMTDMU+5xIEITIgFz8EyQIQmSCa68HzTIIsVLarIOU3QgJghAZ0ebm7LOc\nOgQSBCEyo9n4UwZO0eImITIl5QawO2QhCJExGkMQQlxFSlGQIAixAFKJgsYQhFgIWssghLiK2DMO\nEgQhFkRsK0GCIMTCaG4EExIJghALJJalIEEQYqHEcGeWIAixcEJOSUoQhCiAUKIgQRCiEEKIggRB\niIKYKgoSBCEKY8oMhARBiAIZayVIEIQokLFdBwmCEIUyRhQkCEIUzFBRmFUQttvtnNlHp+T6lVw3\noKz6DfFolCBEpOT6lVw3oMz6+VgL6jIIsSL6piQlCEKIA2gpYzzXMybnyVgIATNzdh5mEwQhRH6o\nyyCEOECCIIQ4QIIghDhAgiCEOECCIIQ44P8BrraVtgFHGlkAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7f32dd3c4b70>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "L = la.cholesky(A_reordered)\n",
        "print(\"%d non-zeros\" % len(np.where(np.abs(L) > prec)[0]))\n",
        "pt.spy(L, marker=\",\", precision=prec)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": []
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.5.1+"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}