mscroggs/lecture7.ipynb

## lecture7.ipynb
{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "authorship_tag": "ABX9TyPQsbolA7O3qFa4n3rIBU4k",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/gist/mscroggs/7c1b4440942fa48fe4a8cab0b9cb4a49/lecture7.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Cqa9yCWDQK4X"
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "A = np.array([[-2.,4,3,4],[-2,3,5,3],[-2,5,3,-2],[-2,4,2,2]])"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "def lu(A):\n",
        "  n = A.shape[0]\n",
        "  assert A.shape[1] == n\n",
        "  L = np.zeros((n, n))\n",
        "  for row in range(n):\n",
        "    L[row,row] = A[row,row]\n",
        "    A[row,:] /= A[row,row]\n",
        "\n",
        "    for i in range(row+1, n):\n",
        "      L[i, row] = A[i, row]\n",
        "      A[i, :] -= A[i, row] * A[row, :]\n",
        "  return L, A"
      ],
      "metadata": {
        "id": "Edm481klQYLZ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "A_copy = A.copy()\n",
        "L, U = lu(A_copy)"
      ],
      "metadata": {
        "id": "nSPGk9_EV5-n"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "L @ U"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "6d0h9DT9V7MN",
        "outputId": "ffdb6f04-68b6-4d9d-dd78-172f40dd34d6"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([[-2.,  4.,  3.,  4.],\n",
              "       [-2.,  3.,  5.,  3.],\n",
              "       [-2.,  5.,  3., -2.],\n",
              "       [-2.,  4.,  2.,  2.]])"
            ]
          },
          "metadata": {},
          "execution_count": 31
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "A"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "X-ihxXmIYeDe",
        "outputId": "98253e8e-e82d-448b-9258-039f748ba876"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([[-2.,  4.,  3.,  4.],\n",
              "       [-2.,  3.,  5.,  3.],\n",
              "       [-2.,  5.,  3., -2.],\n",
              "       [-2.,  4.,  2.,  2.]])"
            ]
          },
          "metadata": {},
          "execution_count": 32
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "L, U"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "YaERSPQ2Ye8j",
        "outputId": "52776c3b-cd48-4126-aed7-4d4b86ee275a"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(array([[-2. ,  0. ,  0. ,  0. ],\n",
              "        [-2. , -1. ,  0. ,  0. ],\n",
              "        [-2. ,  1. ,  2. ,  0. ],\n",
              "        [-2. ,  0. , -1. , -5.5]]), array([[ 1. , -2. , -1.5, -2. ],\n",
              "        [-0. ,  1. , -2. ,  1. ],\n",
              "        [ 0. ,  0. ,  1. , -3.5],\n",
              "        [-0. , -0. , -0. ,  1. ]]))"
            ]
          },
          "metadata": {},
          "execution_count": 33
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "print(\"-----------------------------------------\")"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "YKdpy-yvYj0P",
        "outputId": "50f078b2-59ac-462c-bd2c-e0958143cbf0"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "-----------------------------------------\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "from scipy.sparse import coo_matrix\n",
        "\n",
        "def discretise_poisson(N):\n",
        "    \"\"\"Generate the matrix and rhs associated with the discrete Poisson operator.\"\"\"\n",
        "    \n",
        "    nelements = 5 * N**2 - 16 * N + 16\n",
        "    \n",
        "    row_ind = np.empty(nelements, dtype=np.float64)\n",
        "    col_ind = np.empty(nelements, dtype=np.float64)\n",
        "    data = np.empty(nelements, dtype=np.float64)\n",
        "    \n",
        "    f = np.empty(N * N, dtype=np.float64)\n",
        "    \n",
        "    count = 0\n",
        "    for j in range(N):\n",
        "        for i in range(N):\n",
        "            if i == 0 or i == N - 1 or j == 0 or j == N - 1:\n",
        "                row_ind[count] = col_ind[count] = j * N + i\n",
        "                data[count] =  1\n",
        "                f[j * N + i] = 0\n",
        "                count += 1\n",
        "                \n",
        "            else:\n",
        "                row_ind[count : count + 5] = j * N + i\n",
        "                col_ind[count] = j * N + i\n",
        "                col_ind[count + 1] = j * N + i + 1\n",
        "                col_ind[count + 2] = j * N + i - 1\n",
        "                col_ind[count + 3] = (j + 1) * N + i\n",
        "                col_ind[count + 4] = (j - 1) * N + i\n",
        "                                \n",
        "                data[count] = 4 * (N - 1)**2\n",
        "                data[count + 1 : count + 5] = - (N - 1)**2\n",
        "                f[j * N + i] = 1\n",
        "                \n",
        "                count += 5\n",
        "                                                \n",
        "    return coo_matrix((data, (row_ind, col_ind)), shape=(N**2, N**2)).tocsr(), f"
      ],
      "metadata": {
        "id": "nUAkUyHEbfYe"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "from scipy.sparse.linalg import gmres\n",
        "\n",
        "A, b = discretise_poisson(100)\n",
        "\n",
        "residuals = []\n",
        "def c(res):\n",
        "  global residuals\n",
        "  residuals.append(res)\n",
        "\n",
        "x, _ = gmres(A, b, callback=c, callback_type=\"pr_norm\")\n",
        "\n"
      ],
      "metadata": {
        "id": "B_MWogrCbiXR"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import matplotlib.pylab as plt\n",
        "\n",
        "plt.plot(residuals)\n",
        "plt.yscale(\"log\")"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 265
        },
        "id": "j-hc3EzpbrKc",
        "outputId": "4b54deb0-6289-4e47-ba4a-258d2918a074"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD4CAYAAAAEhuazAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3xUVfrH8c+TTg0goZcgoYhgAUR6dKUrouLuwirWNaIiYNgC6zZXXXV3DTYUUexKEVEBKYI/DEVaEIHQg7QoSBCVJjXn90cGN7KUtJk7M/m+X695OffOzZ3ncOOTM+c+c4455xARkdIlwusAREQk8JT8RURKISV/EZFSSMlfRKQUUvIXESmForwO4GyqVq3qEhMTvQ5DRCSkLF++fI9zLuFsxwR18k9MTCQjI8PrMEREQoqZbTvXMRr2EREphZT8RURKISV/EZFSSMlfRKQUUvIXESmFAlbtY2blgOeBo8Cnzrm3A/XeIiLyc8Xq+ZvZK2a228wyT9nfw8w2mFmWmQ337b4BmOScuwu4tjjvKyIixVPcYZ/XgB75d5hZJDAK6Ak0A/qbWTOgDrDDd9iJYr7vWS3dspexC7b48y1EREJasZK/c24esPeU3W2ALOfcl865o8B4oA+QTd4fgLO+r5mlmFmGmWXk5OQUKa53M3bw8LS1TF35dZF+XkQk3Pnjhm9t/tvDh7ykXxuYDPQ1sxeAqWf6YefcGOdca+dc64SEs347+Yweub45lyVWZtjElXyWtadI5xARCWcBq/Zxzh10zt3unLvnXDd7zay3mY354YcfivResVGRvHzLZSRWLUvKm8vJ/Kpo5xERCVf+SP5fAXXzbdfx7Ssw59xU51xKfHx8kYOILxvNG3dcTnyZaG57dRnbvj1Y5HOJiIQbfyT/ZUAjM2tgZjFAP2BKYU5Q3J7/STXi43j9jjacyM1lwNil5Ow/UqzziYiEi+KWeo4DFgFNzCzbzO50zh0HBgGzgHXAROfcmsKctyR6/iclVSvPK7ddRs7+IwwYu4Ste/QJQETEnHNex/A/zKw30DspKemuTZs2lcg5523MYdA7n3M81/Hg1Rfwmzb1MLMSObeISDAxs+XOudZnOyYop3coyZ7/SZ0bJzDrgc60rFeZB9/PZNC4FRw6erzEzi8iEkqCMvn7S834MrxxRxuG92zKjNU7ueH5z3QjWERKpaBM/iV1w/d0IiKMgckNee32Nuz84TC9np7P659tJTc3+Ia/RET8JSiTvz+GfU7VuXEC04d0olViFf42ZQ23vrqUb/Yd9tv7iYgEk6BM/oFSu1IZXr/9Mh65rjnLtu6l28h5TNGUECJSCgRl8vfnsM9p3oub29Zn+uBONKhajsHjVnDfO5+z9+BRv7+3iIhXgrLU86TWrVu7jIyMgL3f8RO5vDjvS56as5H4MjE8fkMLujSrHrD3FxEpCSFb6umVqMgI7rsyiQ/v60jV8jH89o0Mfv/uSvYfPuZ1aCIiJUrJ/zSa1arIlEEdGXRlEu99nk2Pp+ZrdlARCStBmfwDOeZ/JjFREfyuexPeu6c9sVER/OblJaRO+ILvD+legIiEPo35F8CPR0/w3NxNvJj+JVXKxfB43xb8oqnuBYhIcNKYfwkpExPJ77s35YP7OlCpbDR3vJbB0PErVBEkIiFLyb8QmteOZ+r9HRlyVSOmrdpJ17R0pq36mmD+9CQicjpK/oUUGxXJA10bM21wR2pXLsOgd1Zw95vL2a1vB4tICAnK5B8MN3zPpWmNiky+pz0jejYlfWMOXdLSeTdjhz4FiEhI0A3fEvBlzgH++N4qlm39js6NE3j0uubUrVLW67BEpJTSDd8AOT+hPBNS2vHQtReSsXUvXdLSGTNvMyc0U6iIBCkl/xISEWHc2j6ROanJdG6cwD+nr6ffmEVszjngdWgiIv9Dyb+E1apUhjEDWvGfX17Mhl376fn0fEbNzeLYiVyvQxMR+YmSvx+YGTe2qsOcYcl0uaAa/561gWufW8iq7O+9Dk1EBFDy96tqFeJ4/qZWvDigFXsPHuG6UQt5ZNparR0sIp4LyuQfCqWehdH9whrMTk2mf5t6vLxgC91GziN9Y47XYYlIKaZSzwBbumUvwyev4sucg9xwaW3+ck0zKpeL8TosEQkjKvUMQm0aVGH64E7c/4skpqz8mq4j05m+eqfXYYlIKaPk74G46EiGdWvClEEdqRlfhnvf/pyBmiJCRAJIyd9DzWpV5P172zO8Z1PmbthNl7R0JmqKCBEJACV/j0VFRjAwuSEzhnSiaY2K/GHSKgaMXcqOvYe8Dk1EwpiSf5A4P6E841Pa8vB1zVmx/Tu6jZzH2AVbNEWEiPiFkn8QiYgwBrStz+zUZNo1PI+Hp63lhhc+Y/2ufV6HJiJhJmDJ38zON7OxZjYpUO8ZqmpVKsPYW1vzdL9L2LH3ENc8s4AnZq7n8LETXocmImGiQMnfzF4xs91mlnnK/h5mtsHMssxs+NnO4Zz70jl3Z3GCLU3MjD6X1GZOajLXXVqbFz7dTI+n5rH4y2+9Dk1EwkBBe/6vAT3y7zCzSGAU0BNoBvQ3s2Zm1sLMpp3yqFaiUZciVcrF8J9fXszbv72cXAf9xizmT++vZt/hY16HJiIhrEDJ3zk3D9h7yu42QJavR38UGA/0cc6tds5dc8pjd0EDMrMUM8sws4ycHE2BcFKHpKrMGtqZuzo1YPzS7XRLm8ectd94HZaIhKjijPnXBnbk28727TstMzvPzEYDl5rZiDMd55wb45xr7ZxrnZCQUIzwwk+ZmEgevLoZk+/tQHyZaH77Rgb3j1vBngNHvA5NREJMwG74Oue+dc4NdM41dM49drZjw21it5J2Sd1KTL2/I6ldGzMzcydd09L5YMVX+nKYiBRYcZL/V0DdfNt1fPuKzTk31TmXEh8fXxKnC0sxUREMvqoRHw3uRGLVcgyd8AW3vbpMXw4TkQIpTvJfBjQyswZmFgP0A6aURFDq+Rdc4+oVmDSwPX/r3YxlW/fSdWQ6L6Zv5rhWDhORsyhoqec4YBHQxMyyzexO59xxYBAwC1gHTHTOrSmJoNTzL5zICOP2Dg2Yk5pMp0YJPDZjPdc9v5A1X+uPp4icnubzDzPOOWZk7uKvH67hu0NHSel8PkOuakRcdKTXoYlIgITsfP4a9ik6M6NXi5rMSe1M35Z5Xw7r+fR8fTlMRH5GPf8wtzBrDyMmr2b73kP0b1OPEb2aUjEu2uuwRMSP1POXn74cltL5fCYs207XtHQ+XrPL67BExGPq+Zciq7K/5w+TVrF+1356tajB36+9kGoV4rwOS0RKWMj2/MU/LqqT9+Ww33dvwpx1u+maNk8rh4mUUkGZ/DXs4z/RkRHcd2USM4Z0okn1Cvxh0ipuHruE7d/qy2EipYmGfUqx3FzHO0u38/iM9RzPzWVY1ybc3iGRqMig7BOISAFp2EfOKiLCuLltfWandqZjUlUenb6OG174jHU7tXKYSLhT8hdqxpfhpVta89xvLuXr73+k97ML+PcsrRwmEs6CMvlrzD/wzIxrLqrF7AeS6XNJbUbN3UyvZ+azdMupyziISDjQmL+c1ryNOfzp/dVkf/cjN11ej+E9m1JBXw4TCQka85ci69w4gY8f6MydHRswbul2umrlMJGwouQvZ1Q2Joq/XPPzlcMGvfO5Vg4TCQNK/nJOJ1cOG9a1MR+v+YYuaelMWp6tL4eJhLCgTP664Rt8YqIiuP+qRkwf0pGkhPL87t2V3PLKUq0cJhKidMNXCi031/HWkm08MWM9uQ6GdWvM7R0aEBlhXocmIuiGr/hJRIRxS7tEZqcm067heTzy0TqufW4BGVtVFioSKpT8pchqVSrD2Ftb82z/S9l78Cg3jl7Enz9Yzf7Dx7wOTUTOQclfisXM6H1xLeakJnNHhwa8sySvLFRrBogENyV/KRHlYqP4a++8stBKZaNJeXM59769nN37D3sdmoichpK/lKiTZaEn1wzo8mQ645duV1moSJAJyuSvUs/QdnLNgJlDOnFBzYoMn7ya/i8tZsueg16HJiI+KvUUv8rNdUzI2ME/p6/jyPFchnZpxF2dzidaawaI+I1KPcVzERFG/zb1+CQ1mauaVuNfMzdw7XMLWbnje69DEynVlPwlIKpVjOOFm1vx4oBW7D14hOufX8gj09Zy6Ohxr0MTKZWU/CWgul9Yg9mpyfRvU4+XF2yh28h5pG/M8ToskVJHyV8CrmJcNI9e34KJd7cjJiqCW19ZSuqEL9h78KjXoYmUGkr+4pk2DaowfXAnBv8iiSkrv6ZLWjofrPhKZaEiAaDkL56Ki44ktVsTpg3uSL0qZRk64Qtue3UZ2d9ptlARf1Lyl6DQtEZF3runPX/r3YxlW/fSbeQ8xi7YwolcfQoQ8YeAJn8zu87MXjKzCWbWLZDvLcEvMsK4vUMDZqcmc3mDKjw8bS03PL+QdTv3eR2aSNgpcPI3s1fMbLeZZZ6yv4eZbTCzLDMbfrZzOOc+cM7dBQwEfl20kCXc1a5Uhlduu4yn+11C9nc/0vvZBfxn1gYOHzvhdWgiYaMwPf/XgB75d5hZJDAK6Ak0A/qbWTMza2Fm0055VMv3o3/2/ZzIaZkZfS6pzZzUZPpcUpvn5mbR6+n5LP7yW69DEwkLhZrewcwSgWnOuea+7XbA351z3X3bIwCcc4+d4ecNeByY7Zybc4ZjUoAUgHr16rXatm1bgeOT8DV/Uw5/en81O/b+SP82dRne8wLiy0R7HZZIUArE9A61gR35trN9+87kfqALcKOZDTzdAc65Mc651s651gkJCcUMT8JFp0YJzBramZTO5zNh2Q66pqUzM3On12GJhKyA3vB1zj3jnGvlnBvonBt9puM0q6ecTtmYKP7U6wI+vK8jVcvHMvCtz0l5I4NdP2jNAJHCKm7y/wqom2+7jm9fsTjnpjrnUuLj44t7KglDLerE8+GgDgzv2ZT0jTl0TUvn7SXbyFVZqEiBFTf5LwMamVkDM4sB+gFTihuUev5yLtGREQxMbsisoZ1pXjueB9/PpN+YxWTtPuB1aCIhoTClnuOARUATM8s2szudc8eBQcAsYB0w0Tm3prhBqecvBZVYtRzv3HU5/7rxIjZ8s59eT8/n2U82cfR4rtehiQS1oFzMxcx6A72TkpLu2rRpk9fhSIjYvf8wD01dy0erdtKkegUe79uCS+tV9joskYArSLVPUCb/k7SSlxTFnLXf8JcPM9m17zC3tkvkd92bUD42yuuwRAJGK3lJqdSlWXU+fqAzA9rW5/VFW+k+ch5z1+/2OiyRoBKUyV83fKW4KsRF848+zZk0sB1lYiK5/bVlDB63gj0HjngdmkhQ0LCPhL0jx0/wwqebGTU3i7IxUYzo2ZRfta5LRIR5HZqIX2jYRwSIjYpkaJfGzBjSiSY1KjB88mr6v7SYLXsOeh2aiGeCMvlr2Ef8IalaBSaktOWJvi1Yu3Mf3Z+ax/OfZnHshMpCpfTRsI+USrv3HeZvU9YwI3MXF9SsyBN9W3BRnUpehyVSIjTsI3IG1SrG8cLNrRh9cyu+PXCE60Yt5OFpazl09LjXoYkEhJK/lGo9mtdgzrBk+repx9gFW+iaNo9PN6gsVMJfUCZ/jflLIFWMi+bR61sw8e52xEVHcNuryxg6fgXfqixUwpjG/EXyOXL8BKPmbuaFT7MoHxvFn69uxg0ta5O3DpFIaNCYv0ghxUZFktq1MR8N7kSDquUY9u5KbnllKdu/PeR1aCIlSslf5DQaV6/ApIHtebjPhazY/j3dnkpnzLzNHFdZqIQJJX+RM4iIMAa0S2R2amc6JiXwz+nr6TNqIZlf6V6UhL6gTP664SvBpGZ8GV66pRWjftOSb/Ydoc+ohTw2fR0/Hj3hdWgiRaYbviKF8MOhYzw2Yx3jl+2gXpWyPHp9czo1SvA6LJGf0Q1fkRIWXzaax/texLi72hIZYQwYu5TUCV+w9+BRr0MTKRQlf5EiaNfwPGYM6cSgK5OYsvJruqSl8/6KbIL5k7RIfkr+IkUUFx3J77o3YdrgjtSrUpYHJqgsVEKHkr9IMTWtUZH37mnPQ9deyOfbvlNZqIQEJX+REhAZYdzaPpHZqck/Kwtdna2KNQlOQZn8VeopoapWpbyy0Odvasnu/UfoM2oBj2i2UAlCKvUU8ZMffjzGEzPX886S7dSpXIZHr29BcmOVhYr/qdRTxEPxZaL5p2+20NioCG59ZSlDxmsReQkOSv4iftamQRWmD+nEkKsaMX31TrqkpfNuxg6VhYqnlPxFAiA2KpIHujZm+uBOJCWU5/eTVnHTy0vYqkXkxSNK/iIB1Kh6BSbe3Y5HrmvO6uwftIi8eEbJXyTAIiKMm9vWZ3ZqMlc2qca/Zm6g97ML+GLH916HJqWIkr+IR2rExzF6QCteHNCK7w4d5frnF/LQ1DUcOKKyUPG/gCV/M7vAzEab2SQzuydQ7ysS7LpfWIPZqcncfHl9XvtsK93S0pmz9huvw5IwV6Dkb2avmNluM8s8ZX8PM9tgZllmNvxs53DOrXPODQR+BXQoesgi4adiXDQPX9ecSQPbUT4uit++kcF9b3/O7v2HvQ5NwlRBe/6vAT3y7zCzSGAU0BNoBvQ3s2Zm1sLMpp3yqOb7mWuBj4DpJdYCkTDSqn4Vpt3fid91a8zsdd/Q5cl0Ji5TWaiUvAJ/w9fMEoFpzrnmvu12wN+dc9192yMAnHOPFeBcHznnrj7XcfqGr5Rmm3MOMOK91Szdupf2Dc/jn9e3ILFqOa/DkhDg72/41gZ25NvO9u07UzBXmNkzZvYiZ+n5m1mKmWWYWUZOTk4xwhMJbQ0TyjM+pS2PXv/fstAXPt2sslApEVGBeiPn3KfApwU4bgwwBvJ6/v6NSiS4RUQYN11en6uaVudvUzJ5YuZ6pq78mif6XkSLOvFehychrDg9/6+Auvm26/j2FZtm9RT5uRrxcbw4oDWjb27JngOaLVSKrzjJfxnQyMwamFkM0A+YUhJBOeemOudS4uPVsxHJr0fzmsxOTaZfm3q8vGAL3UbOI32jhkel8Apa6jkOWAQ0MbNsM7vTOXccGATMAtYBE51za0oiKPX8Rc7s5GyhE1LaEuObLfQBLSIvhaT5/EVC2OFjJ3h+bhYvpG+mfGwUf7mmGddfWhsz8zo08VDIzuevnr9IwcRFR5LarQnT7u9EYtVypE7MW0R+x14tIi9np56/SJg4ket4a/E2/jVzPbkOhnVrzG3tE4mKDMo+nviRev4ipUj+ReTbNzyPRz5ax3XPL2SlZguV01DPXyQMOef4aPVOHpq6lj0HjnBL2/oM696EinHRXocmARCyPX8RKR4z45qLavHJsGRuaVufNxZvo2taOh+v2eV1aBIklPxFwljFuGge6tOc9+/tQOWyMaS8uZx7316u2UIlOJO/xvxFStYldSsx9f6O/L57E+as202XJ9OZsGy7ZgstxTTmL1LKbM45wIjJq1m6JW+20MduaEH98zRbaDjRmL+I/I+GCeUZf1febKGrfLOFvjz/S07kBm9HUEpeUCZ/DfuI+NfJ2UJnp3amQ8OqPPLROm544TPW7dzndWgSIBr2ESnlnHNMWfk1/5i6lu9/PMZvOzVg6FWNKRMT6XVoUkQa9hGRczIz+lxSmzmpyfRtWZsX07+k21Ppmi00zCn5iwgAlcvF8K8bL2Z8SluiI/NmCx08bgU5+494HZr4gZK/iPxM2/PPY8aQTgy5qhEzMnfSJS2d8Uu3k6sbwmElKJO/bviKeCs2KpIHujZmxpDONKlRgeGTV9PvpcVszjngdWhSQnTDV0TOKjfX8e7yHfxz+np+PHqC+3+RxN3JDYmJCsq+o6AbviJSAiIijF9fVo85qcl0u7A6T87eSO9nF/D59u+8Dk2KQclfRAokoUIsz/2mJWNvbc2+w8fo+8Jn/H3KGg4c0SLyoUjJX0QK5aoLqjM7NW+20NcXbaVbWjqfrPvG67CkkJT8RaTQysdG8VCf5kwa2J7ycVHc+XoG9739Obv3abbQUKHkLyJF1qp+Zabd34lhXRsze903XJWWzluLt6ksNAQEZfJXqadI6IiJiuD+qxoxc0gnmteK588fZPKrFxeRtXu/16HJWajUU0RKjHOOScuzeeSjdfx49AT3XZnEPVeoLDTQVOopIgFlZvyydV0+GZZMj+Y1GDlnI9c8O5/l21QWGmyU/EWkxFUtH8sz/S/lldtac+DwcW4c/Rl//TCT/YePeR2a+Cj5i4jf/KJpdT5OTebWdom8uXgb3UbOY85alYUGAyV/EfGr8rFR/P3aC3nvnvZUjIvmt2/4ykK1iLynlPxFJCBa1qvM1Ps78rtueWWhXZ7Mmy00mItOwpmSv4gETExUBIN+kVcWekHNigyfvJpfj1msslAPKPmLSMCdn1Ce8SlteaJvCzbs2k/Pp+fz71l5s4ZKYAQ0+ZtZOTPLMLNrAvm+IhJ8zPJmC/1kWDK9L67FqLmb6f7UPBZm7fE6tFKhQMnfzF4xs91mlnnK/h5mtsHMssxseAFO9UdgYlECFZHwVLV8LGm/uoR37rqcCIObXl7C799dyfeHjnodWlgraM//NaBH/h1mFgmMAnoCzYD+ZtbMzFqY2bRTHtXMrCuwFthdgvGLSJho37AqM4d25p4rGjJ5xVd0SZvHR6t26oawnxR4egczSwSmOeea+7bbAX93znX3bY8AcM49doaffxQoR94fih+B651zuac5LgVIAahXr16rbdu2Fa5FIhLy1nz9A398bxWZX+2jywXVeeS65tSIj/M6rJDh7+kdagM78m1n+/adlnPuQefcUOAd4KXTJX7fcWOcc62dc60TEhKKEZ6IhKoLa8Xzwb0d+FOvpizIyqFrWjpvarbQEhXwah/n3GvOuWlnO0azeopIVGQEKZ0bMmtoZy6qG89ffLOFbvpGZaEloTjJ/yugbr7tOr59xeacm+qcS4mPjy+J04lICKt/XjneuvNy/n3jRWTlHKDXM/MZOXsjR46rLLQ4ipP8lwGNzKyBmcUA/YApJRGUev4ikt/J2ULnpCbTs3lNnv5kE72ens+yrXu9Di1kFbTUcxywCGhiZtlmdqdz7jgwCJgFrAMmOufWlERQ6vmLyOmcnC301dsv4/CxXH45ehF/en81+zRbaKEF5WIuZtYb6J2UlHTXpk2bvA5HRILQwSPHGTl7I68s3ELV8rH8o8+F9Ghe0+uwgkJBqn2CMvmfpJW8RORcVmV/z/D3VrN25z66NavOP/qoLFQreYlI2LuoTiU+HNSBET2bMm9TDl3S0nlz0VaVhZ5DUCZ/3fAVkcKIjozg7uS8stBL6lbiLx+u4cbRn7Fu5z6vQwtaGvYRkbDinGPy51/x6PR1/PDjMW5rn0hq18aUi43yOrSA0bCPiJQ6ZkbfVnX4v2HJ/Kp1XcYu2EK3kfP4dIOmFcsvKJO/hn1EpLgqlY3hsRta8O7AdsRFR3Dbq8sYMn4Few4c8Tq0oKBhHxEJe0eOn+D5uZt5/tMsysVG8eerm9G3ZW3MzOvQ/ELDPiIiQGxUJA90bcz0wZ1omFCe3727kgFjl7Lt24Neh+YZJX8RKTUaVa/Au3e34+HrmvPFju/p/tQ8Rqdv5tiJ004yHNaCMvlrzF9E/CUiwhjQtj5zUpPp3CiBx2es59rnFrIq+3uvQwsojfmLSKk2M3MXf/0wkz0HjnB7hwZhURaqMX8RkXPo0bwGc4Yl079NvZ/KQueWgrJQJX8RKfUqxkXz6PV5ZaFlYiK5/dVlDB4X3mWhSv4iIj6XJVbho8EdGdqlETMyd3LVk+lMzNgRlovIB2Xy1w1fEfFKbFQkQ7s0ZsaQTjSqVp4/TFrFTS8vYeue8CoL1Q1fEZEzyM11jFu2ncenr+foiVyGdGnEXZ3OJzoyKPvNP9ENXxGRYoiIMG66vD5zhiVzZZNq/GvmBq5+Zj4ZYbB8pJK/iMg5VK8Yx+gBrXj5ltYcPHKCG0cv4s8fhPbykUr+IiIF1KVZdT5+oDN3dGjAO0u20zUtnVlrdnkdVpEo+YuIFEK52Cj+2rsZ79/bgcplY7j7zeUMfHM53+w77HVohaLkLyJSBBfXrcTU+zvyhx5NmLthN13S0nl7ybaQWT4yKJO/Sj1FJBRER0Zw7xVJzBzamea14nnw/Uz6jVlM1u4DXod2Tir1FBEpAc453l2ezaMfrePHoye478ok7rmiITFRge9jq9RTRCRAzIxfta7LnNRkujevwcg5G4O6LFTJX0SkBCVUiOXZ/pfy6m2Xceho8JaFKvmLiPjBlU2r/U9Z6MzM4CkLVfIXEfGT/GWhVcrFMvCt5dz9Zga7fvC+LFTJX0TEzy6uW4kpgzrwxx5N+XRDDl3T0nlzsbdloUr+IiIBEB0ZwT1XNOTjBzpzUd14/vJBJr98cREbv9nvSTxK/iIiAVT/vHK8deflPPnLi9mcc4Crn5lP2scbOHzsREDjCFjyN7MrzGy+mY02sysC9b4iIsHGzOjbqg6fpCZzzUW1eOb/suj1zHyWfPltwGIoUPI3s1fMbLeZZZ6yv4eZbTCzLDMbfo7TOOAAEAdkFy1cEZHwcV75WEb++hJev6MNR4/n8usxixn+3ip+OOT/stACfcPXzDqTl7jfcM419+2LBDYCXclL5suA/kAk8Ngpp7gD2OOcyzWz6kCac+6mc72vvuErIqXFoaPHeWrOJl6e/yVVysXy4oCWtKpfpUjnKsg3fKMKciLn3DwzSzxldxsgyzn3pe/NxgN9nHOPAdec5XTfAbFnCToFSAGoV69eQcITEQl5ZWOi+FOvC7j24lo8MXM99c8r59f3K1DyP4PawI5829nA5Wc62MxuALoDlYDnznScc24MMAbyev7FiE9EJOQ0rx3Pm3eeMZWWmOIk/0Jxzk0GJhfkWDPrDfROSkryb1AiIqVUcap9vgLq5tuu49tXbM65qc65lPj4+JI4nYiInKI4yX8Z0MjMGphZDNAPmFISQWk+fxER/ypoqec4YBHQxMyyzexO59xxYBAwC1gHTHTOrSmJoNTzFxHxr4JW+/Q/w/7pwPQSjQiN+YuI+FtQTu+gngBA0XwAAAUJSURBVL+IiH8FZfIXERH/Csrkrxu+IiL+FdQLuJtZDrCtiD9eFdhTguF4LZzaE05tgfBqTzi1BcKrPYVpS33nXMLZDgjq5F8cZpZxrrktQkk4tSec2gLh1Z5waguEV3tKui1BOewjIiL+peQvIlIKhXPyH+N1ACUsnNoTTm2B8GpPOLUFwqs9JdqWsB3zFxGRMwvnnr+IiJyBkr+ISCkUlsm/kGsLe87M6prZXDNba2ZrzGyIb38VM5ttZpt8/63s229m9oyvfavMrKW3LfhfZhZpZivMbJpvu4GZLfHFPME3EyxmFuvbzvK9nuhl3KdjZpXMbJKZrTezdWbWLsSvzQO+37NMMxtnZnGhcn1Ot554Ua6Fmd3qO36Tmd3qRVt8cZyuPf/2/a6tMrP3zaxSvtdG+Nqzwcy659tf+JznnAurB3lrCG8GzgdigJVAM6/jOkfMNYGWvucVyFsbuRnwL2C4b/9w4Anf817ADMCAtsASr9twmjalAu8A03zbE4F+vuejgXt8z+8FRvue9wMmeB37adryOvBb3/MY8lajC8lrQ94KfFuAMvmuy22hcn2AzkBLIDPfvkJdC6AK8KXvv5V9zysHUXu6AVG+50/ka08zXz6LBRr48lxkUXOe57+MfvjHbAfMyrc9AhjhdVyFbMOHQFdgA1DTt68msMH3/EWgf77jfzouGB7kLezzCfALYJrvf749+X6hf7pG5E0J3s73PMp3nHndhnxtifclSztlf6hem5PLr1bx/XtPI2951ZC5PkDiKcmyUNcC6A+8mG//z47zuj2nvHY98Lbv+c9y2clrU9ScF47DPqdbW7i2R7EUmu9j9aXAEqC6c26n76VdQHXf82Bv41PAH4Bc3/Z5wPcubw0I+Hm8P7XF9/oPvuODRQMgB3jVN4z1spmVI0SvjXPuK+A/wHZgJ3n/3ssJ3esDhb8WQX2NTnEHeZ9eoITbE47JP2SZWXngPWCoc25f/tdc3p/0oK/LNbNrgN3OueVex1JCosj7WP6Cc+5S4CB5Qws/CZVrA+AbD+9D3h+1WkA5oIenQZWgULoW52JmDwLHgbf9cf5wTP5+W1vYn8wsmrzE/7bLW+we4Bszq+l7vSaw27c/mNvYAbjWzLYC48kb+nkaqGRmJxcPyh/vT23xvR4PfBvIgM8hG8h2zi3xbU8i749BKF4bgC7AFudcjnPuGDCZvGsWqtcHCn8tgv0aYWa3AdcAN/n+oEEJtycck7/f1hb2FzMzYCywzjmXlu+lKcDJSoRbybsXcHL/Lb5qhrbAD/k+9nrKOTfCOVfHOZdI3r/9/znnbgLmAjf6Dju1LSfbeKPv+KDpuTnndgE7zKyJb9dVwFpC8Nr4bAfamllZ3+/dyfaE5PXxKey1mAV0M7PKvk9C3Xz7goKZ9SBv2PRa59yhfC9NAfr5KrAaAI2ApRQ153l548aPN1B6kVcxsxl40Ot4ChBvR/I+qq4CvvA9epE3tvoJsAmYA1TxHW/AKF/7VgOtvW7DGdp1Bf+t9jnf94uaBbwLxPr2x/m2s3yvn+913KdpxyVAhu/6fEBehUjIXhvgIWA9kAm8SV71SEhcH2AcefcqjpH3qezOolwL8sbSs3yP24OsPVnkjeGfzAWj8x3/oK89G4Ce+fYXOudpegcRkVIoHId9RETkHJT8RURKISV/EZFSSMlfRKQUUvIXESmFlPxFREohJX8RkVLo/wFjmXzWVEnl+wAAAABJRU5ErkJggg==\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "from scipy.sparse.linalg import gmres, spilu, LinearOperator\n",
        "\n",
        "A, b = discretise_poisson(100)\n",
        "\n",
        "M = LinearOperator(\n",
        "    matvec = spilu(A, fill_factor=20, drop_rule='dynamic').solve,\n",
        "    shape=(10000, 10000), dtype=np.float64\n",
        ")\n",
        "\n",
        "residuals2 = []\n",
        "def c(res):\n",
        "  global residuals2\n",
        "  residuals2.append(res)\n",
        "\n",
        "x, _ = gmres(A, b, M=M, callback=c, callback_type=\"pr_norm\")"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "P4I_MTXBccob",
        "outputId": "5c93e5d7-eec9-49d8-c5a8-a7ca2ed8c0c7"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "/usr/local/lib/python3.8/dist-packages/scipy/sparse/linalg/dsolve/linsolve.py:407: SparseEfficiencyWarning: splu requires CSC matrix format\n",
            "  warn('splu requires CSC matrix format', SparseEfficiencyWarning)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import matplotlib.pylab as plt\n",
        "\n",
        "plt.plot(residuals)\n",
        "plt.plot(residuals2)\n",
        "\n",
        "#plt.xlim([0,50])\n",
        "plt.yscale(\"log\")"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 265
        },
        "id": "ud8Fc_pTdORk",
        "outputId": "e0fd0a47-e83a-410b-a724-774f81384b7e"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD4CAYAAAAEhuazAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXgUVb7G8e+vsyGL7CBbWASRKCoQAUVIoiIgm7vgrijivgCOjjP3jjNzR0cQFEUBFXFBFHdEEBVJAopAQJEdAgiERcAFBUQInPtHd5yQCZCkO+nqzvt5njzprlS6fkWFt09OnZxjzjlERKR88YW7ABERKXsKfxGRckjhLyJSDin8RUTKIYW/iEg5FBvuAo6mVq1arkmTJuEuQ0QkoixcuHCnc6720fbxdPg3adKErKyscJchIhJRzGzDsfZRt4+ISDmk8BcRKYcU/iIi5ZDCX0SkHFL4i4iUQ2U22sfMKgHPAvuBdOfcxLI6toiIHC6olr+ZjTez7Wa2tMD27ma2ysyyzezBwOZLgLedc7cAfYI5roiIBCfYbp8JQPf8G8wsBhgN9ACSgP5mlgQ0BDYFdjsY5HGPav76Hxk/Z31pHkJEJKIFFf7OuUzgxwKb2wPZzrl1zrn9wBtAXyAH/xvAUY9rZgPNLMvMsnbs2FGiuiZnbeLvU5czNmNtib5fRCTalUaffwP+08IHf+h3AEYBz5hZT+DDI32zc26cmW0FesfHx7crSQGPXtKa3w4c5NHpK9n9ey73dz0JMyvJS4mIRKUyu+HrnNsD3FjEfT8EPkxOTr6lJMeKi/Exql8bKsfH8vTn2fy6L5f/6ZWEz6c3ABERKJ3w3ww0yve8YWBbkZlZb6B38+bNS1xEjM947NLWVKkQywtz1vPLvgM8fulpxMZodKuISGkk4QKghZk1NbN4oB8wpTgv4Jz70Dk3sGrVqkEVYmY83LMVg7uexLuLNjPotUXs3P17UK8pIhINgh3qOQmYC7Q0sxwzG+CcywXuBGYAK4DJzrllxXzd3mY2bteuXcGUl/da3HVeCx7pcwrpq7aTNjydl75YT+7BQ0G/tohIpDLnXLhrOKLk5GQXyimds7fv5pEPlzF7zU5a1TueRy9pzRmNqoXs9UVEvMDMFjrnko+2jyc7wEPZ8s+veZ3KvHJTe8Zc05Yf9/zOxc9+wf9+sJRf9x0I6XFERLyuXLX88/t13wGGz1jFK19toG6VCvytTxLdTjlBQ0JFJOJFbMu/LFSpEMcjfU/lvds7Ub1SPINeW8Q1L85j5bZfwl2aiEip82T4l1a3T2HOaFSND+/sxCN9TmHZll+48KnZPPzeEn7QqCARiWLlttunMD/v3c+Tn63h1a82UDE+hnvOa8F1ZzUhPtaT75EiIoVSt08xVasYz9/6nMKMezvTNrE6//xoBd2fzGTmiu/x8pukiEhxeTL8y7LbpzDN61Th5Zva89INZ4LBgJezuG78fFZ//2tY6hERCTV1+xzDgYOHeGXuBp76bDV79h/k6g6J3Hf+SVSvFB/WukREjkTdPiEQF+NjwDlNSR+axlXtE3ntqw2kBv5KeH+u/kpYRCKTwr+IalSK5x8Xncr0e7rQukFVHvlwOReMzGD6kq26HyAiEUfhX0wtT6jCqwPaM/6GZOJifNw2cRGXjZnLwg0/hbs0EZEi82Sff74pnW9Zs2ZNuMs5otyDh3hrYQ5PfLKanbt/p2frejzQvSWNa1YKd2kiUo4Vpc/fk+Gfxws3fItiz++5jMtcx7jMdeQeOsQ1HRtz97ktdFNYRMJC4V/Gvv9lHyM/Xc3krE1USojlzrTmXH92EyrExYS7NBEpRzTap4zVPb4Cj116GtPv6UJy4+o8On0l5z2Rwftfb+bQIe++yYpI+aPwLwUtT6jCSze2Z+LNHah6XBz3vvkNfUd/wdy1P4S7NBERQOFfqjo1r8XUu85hxBWns3P37/R//isGTFhA9nb9pbCIhJcn+/wjZbRPcew7cJDxX6znuVlr2XvgIFee2Yh7z29BnSoVwl2aiESZ8nvDd/kU+Ok76HR3yGsK1g+7f2fUzDVMnLeR+Fgft3Y5kVu6NKVifGy4SxORKFF+b/iungHzxoa7ikLVrJzAI31P5ZP7utClRW1Gfraa1GHpvLlgIwd1U1hEykh0hr8vBg55e13eZrUrM+badrw96CzqVzuOP72zhJ6jZpO+arumixCRUhed4R8TB4dyw11FkSQ3qcF7t5/N6Kvasnf/QW54aQHXvjifZVvCM521iJQP0Rn+vlg4GBnhD2Bm9DytHp/e34W/9kpi6ZZd9Hp6DoMnL2bbrn3hLk9EolCZhb+ZNTOzF83s7VI/mC82Ylr++SXExjDgnKZkDEnjls7N+HDxFlKHz+LfH69k115vd2OJSGQpUvib2Xgz225mSwts725mq8ws28wePNprOOfWOecGBFNskfliPd/nfzRVK8bx5wtb8dn9KVyQdAJjMtbS+fHPGZe5ln0HDoa7PBGJAkVt+U8AuuffYGYxwGigB5AE9DezJDNrbWZTC3zUCWnVxxKhLf+CEmtWZFT/Nky7uzNtG1fnX9P800V88I2mixCR4BQp/J1zmcCPBTa3B7IDLfr9wBtAX+fcEudcrwIf20Nc99HFxIE7BIeiY6WtVvWOZ0K+6SLueeMbLn72C+at03QRIlIywfT5NwA25XueE9hWKDOraWZjgDZm9tBR9htoZllmlrVjx46SVeYLzKIZBa3//PKmi3ji8tPZ/uvvXDnuK255JYu1O3aHuzQRiTBl9melzrkfgEFF2G+cmW0FesfHx7cr0cF8cf7Ph3KB6JpT3+czLm3XkJ6n1ePFOet5Ln0tF4zM5OoOidx9XgtqVU4Id4kiEgGCaflvBhrle94wsC1ozrkPnXMDq1atWrIX8AXe0yL4pu+xVIiL4Y605qQPTeWq9olMnLeRzv+exWPTV7Lrt+g9bxEJjWDCfwHQwsyamlk80A+YEoqizKy3mY3btauEf+gUk9fyj/6RMbUqJ/CPi/zTRXQ7pS5jM9eSOmwW4+esZ39udNzzEJHQK+pQz0nAXKClmeWY2QDnXC5wJzADWAFMds4tC0VRwbf8A33+B8tPC/jE2pV5sl8bPrzzHJLqH8/fpy6n68gMpi3ZqukiROS/FKnP3znX/wjbpwHTQloRh03pXLIXOKzPv3w5tUFVXhvQgfTVO3hs2kpun7iItonVeLhnK9o1rhHu8kTEIzw5vYP6/INjZqS1rMO0ezrz70tbk/PTb1z63FwGvbqQdRoZJCJ4NPxD1ucfQfP7lIYYn3HlmYmkD03lvvNPInPNDrqOzOSv7y9lx6+/h7s8EQkjT4Z/0C3/P8J/f+iKimAV42O55/wWZAxNo3/7Rrw+fyOpw2bx1Gdr2PN7+X6DFCmvPBn+QYsJjO1X+B+mdpUE/nlRaz65rwud8xaSGZ7O6/M2kntQI4NEyhNPhn/w3T554V8++/yP5cTAQjLv3HYWiTUq8uf3ltDtyUw+WbZNI4NEyglPhr+6fcpGu8Y1eHvQWYy9th0OGPjqQq4YO5eFG34Kd2kiUso8Gf5BU7dPkZkZ3U45gRn3duGfF53K+p17ufS5Lxn06kLNGSQSxTwZ/ur2KXtxMT6u6diYjMDIoNlrdnDByEwefm8J23/VamIi0caT4a9un/CplOAfGZQ+NI2rOyTy5oJNpA5LZ8Snq9mtkUEiUcOT4R80dfsErXaVBP7e91Q+vT+F1Ja1GTVzDanDZvHq3O84oJFBIhEvysNf3T7BalqrEs9e3Y73bj+bZrUr89cPlnHByEzNGSQS4TwZ/qGb1VPhHyptEqvz5sCOvHh9MnExxu0TF3Hxs19qNTGRCOXJ8A++zz+woEmublSGkplxXqu6TL+nC49fehrbdu3jynFfMWDCAlZ//2u4yxORYvBk+ActNi/81edfGmJ8xhVnNmLWkFSGdmvJ/PU/0v3JTP709rds26U3XJFIEKXhX8H/WS3/UnVcvH81sYwH0rjh7Ka8+3UOqcNn8fjHK/lln7rcRLwsOsM/74ZvrmauLAs1KsXzP72T+HxwKt1OOYFn09eS8vgsXpyznt9zo381NZFIFJ3h7/P53wDU8i9TjWpU5Kl+bZh6l381sX9MXc75IzL44JvNHDqkkUEiXuLJ8A96tA/4u37U8g+LvNXEXr6pPZUT4rjnjW/oM3oOX2TvDHdpIhLgyfAPerQP+G/6HlT4h4uZkXJSbT666xxGXHE6P+05wNUvzOO68fNZvuWXcJcnUu55MvxDQi1/T/D5jEvaNmTm4BQevrAVizf9TM+nZ3P/5G/Y/PNv4S5PpNyK4vBPgAMKF6+oEBfDLV2akTk0jYGdmzH1262kDU/nX9NWsGuvRgaJlLXoDf+YBM3t40FVK8bx0IWtmDUkld6n1ef52evo/PjnjM1Yy74DGhkkUlaiN/xjE9Tt42ENqh3HE1eczrS7O9MmsTqPTl/JucPTeWdhDgc1Mkik1EV5+Guop9e1qnc8L9/Untdv7kDNygkMfmsxPUfNJn3Vdk0cJ1KKyjT8zewiM3vezN40swtK9WCx6vaJJGc3r8UHd3RiVP827Nmfyw0vLeDqF+axJCeI4b4ickRFDn8zG29m281saYHt3c1slZllm9mDR3sN59z7zrlbgEHAlSUruYhi1PKPND6f0ef0+sy8P5X/7Z3Eym2/0vuZOdw16Ws2/rA33OWJRJXYYuw7AXgGeCVvg5nFAKOBrkAOsMDMpgAxwKMFvv8m59z2wOO/BL6v9KjPP2LFx/q4sVNTLm3XkLEZa3lxzno+XrqVazo25q5zW1CjUny4SxSJeEUOf+dcppk1KbC5PZDtnFsHYGZvAH2dc48CvQq+hpkZ8Bgw3Tm3qLDjmNlAYCBAYmJiUcv7bwr/iHd8hTiGdjuZazs24cnPVvPyl9/xdlYOg1JP5KZOTTkuPibcJYpErGD7/BsAm/I9zwlsO5K7gPOBy8xsUGE7OOfGOeeSnXPJtWvXLnllCv+ocULVCjx26WnMuLcLHZrVZNiMVaQOn8Ub8zeSqyUlRUqkTG/4OudGOefaOecGOefGHGm/0M3toz7/aNKibhVeuD6ZybeeRf1qx/Hgu0vo8dRsPl3+vUYGiRRTsOG/GWiU73nDwLaghGZuH4V/tGrftAbv3nY2z13dltxDjlteyeLKsV+xaONP4S5NJGIEG/4LgBZm1tTM4oF+wJRgiwpZy//Ab6AWYVQyM3q0rscn93XhHxedyrqde7jk2S+57bWFrNuxO9zliXhecYZ6TgLmAi3NLMfMBjjncoE7gRnACmCyc25ZsEWFpOUfVwFwGusf5eJifFzbsTEZQ1O59/wWZKzeQdeRmfzl/SXs+FX3fESOpDijffofYfs0YFrIKsLf8gd6N2/evOQvEnuc//OB3/6zpq9ErUoJsdx7/klc3aExo2au4fX5G3l30WZu6dyMgV2aUSmhOKOaRaKfJ6d3CF3LH/X7lzO1qyTwj4tO5dP7upByUm2emrmGlGHpvPrVBg5oZJDIHzwZ/iHp84+r6P+saZ3LpWa1K/PcNe149/azaVarEn99fyndRmby8dKtGhkkgkfDP2SjfUDhX861TazOm7d25IXrkonxGYNeW8Slz33Jgu9+DHdpImHlyfAPibyWf67Cv7wzM85Pqsv0ezrz70tbs/nn37h8zFxufjmL7O2/hrs8kbDwZPiHptsn3w1fESA2xseVZyaSPiSNod1a8tW6H7hgZCYPvvMt3/+ie0NSvngy/ENzw1d9/lK44+JjuCOtORlDU7nurCa8syiHlGGzGD5jFb/u05KSUj54MvxD4o+Wv6YClsLVrJzA3/qcwsz7U+madALPzMomZVg6L32xnv25Ghkk0c2T4a9uHylLiTUr8nT/Nky5sxMnn1CFRz5czvkjMpiyeAuHtKSkRClPhn9ou33U8peiOa1hNSbe3IEJN55JxfgY7p70NX1Hf8GX2TvDXZpIyHky/EMi74+8DuhGnhSdmZHasg4f3d2Z4Zefzg+7f+eqF+Zxw0vzWbH1l3CXJxIy0Rv+MYHVnjS3j5RAjM+4rF1DPh+SykM9TmbRhp+4cNRsBk9ezOaf1ZUokc+T4R+SPv8/wl+jN6TkKsTFcGvKiWQ+kMYtnZvx4bdbSBuezqPTVrBrr362JHJ5MvxD0ufviwUMDmpmRwletYrx/PnCVnw+OIVep9Vj3Ox1dBk2i3GZa9l34GC4yxMpNk+Gf0iY+Vv/6vaREGpYvSIjrjiDj+7qzOmNqvGvaSs574kM3l2Uo5FBElGiN/zBP5Wzun2kFCTVP55XbmrPxJs7UKNSPPdPXkyvp+cwe82OcJcmUiTRHf4xcVrEXUpVp+a1+OCOTjzV7wx+2XeAa1+cz7UvzmP5Fo0MEm+L8vBXt4+UPp/P6HtGA2YOTuEvPVvxbc4uej49m/vf/Iacn/R3JuJNngz/kIz2gUD4q9tHykZCbAw3d25G5tA0BnZpxtQlWzn3iQz+76Pl/LRHjRDxFk+Gf0hG+0Ag/NXtI2WrasU4HurRivQhqfQ5vT4vzFlPl2GzGD0rm9/2a2SQeIMnwz9k1PKXMKpf7TiGX346H9/ThQ5NazBsxipSh89i0vyN5GpJSQmz6A5/nw+c/pNJeLU8oQovXH8mbw06i4bVK/LQu0vo/tRsPlm2TUtKSthEd/ibwl+848wmNXh70FmMvbYdh5xj4KsLuXzMXBZu+CncpUk5FP3hf0h9rOIdZka3U07gk3u78K+LW7Phx71c+tyX3PpqFmt37A53eVKORHn4x6jlL54UG+Pjqg6JZAxNZXDXk/gi27+k5J/fW8J2LSkpZaDMwt/MWpnZGDN728xuK5uDqttHvK1ifCx3ndeCjKGpXNuxMZMXbCJlWDojPlnF7t9zw12eRLEihb+ZjTez7Wa2tMD27ma2ysyyzezBo72Gc26Fc24QcAXQqeQlF4PCXyLEH0tKDk7hvFZ1GPV5NimPz+LlL7/TkpJSKora8p8AdM+/wcxigNFADyAJ6G9mSWbW2symFvioE/iePsBHwLSQncHRKPwlwjSuWYlnrmrLB3d0okXdyvzvlGV0HZnBh1pSUkKsSOHvnMsEfiywuT2Q7Zxb55zbD7wB9HXOLXHO9SrwsT3wOlOccz2Aq490LDMbaGZZZpa1Y0eQk2T51Ocvken0RtWYdEtHXrrxTI6Li+GuwJKSc9ZoSUkJjWD6/BsAm/I9zwlsK5SZpZrZKDMby1Fa/s65cc65ZOdccu3atYMoD/+0zgp/iVBmRlpgScknLj+dH/fs55oX53Hti/NYtiXIqU+k3IstqwM559KB9KLsa2a9gd7NmzcP7qAa6ilRIMZnXNquIT1Pq8drX23g6c+z6fX0HC46owGDLziJhtUrhrtEiUDBtPw3A43yPW8Y2OYdGuopUaRCXGDiuAfSuLXLiXy0ZCvnDvdPHPfzXk0cJ8UTTPgvAFqYWVMziwf6AVNCUVTIJnbTDV+JQlWPi+PBHif7J447IzBx3OOzGJOhJSWl6Io61HMSMBdoaWY5ZjbAOZcL3AnMAFYAk51zy0JRVMimdFb4SxTLmzhu2t2dadu4Oo9NX0na8HTeytrEQY0MkmMwL08slZyc7LKyskr+Aq/3g19yYNCc0BUl4lFfrt3JY9NX8m3OLlrWrcKferQkrWUdzCzcpUkZM7OFzrnko+0T3dM7+GLAw29uIqF09on+JSWfuaoN+3IPctOELPqN+4qvN2riOPlvngz/0HX7aKinlC9mRq/T6vPpfSn8ve8pZG/fzcXPfsntExeyfueecJcnHuLJ8A/pDV8N9ZRyKD7Wx3VnNSHjgTTuPq8F6at20HVEBn95fwk7ftXqduLR8NcNX5HQqJwQy/1dTyJ9aCr92jdi0vxNpAybxchPV2viuHLOk+Efupa/xvmLANSpUoF/XtSaT+/rQspJtXlq5hpSh83i1bnfcUBLSpZLngz/kFHLX+QwzWpX5rlr2vHu7WfTrHZl/vrBMrqOyOCjb7dqSclyphyEv/r8RQpqm1idNwd25MXrk4mP9XHH64u46Nkvmbv2h3CXJmXEk+Efsj5/DfUUOSIz47xWdZl+Txcev+w0tv+yj/7Pf8WNL81n5bZfwl2elDJPhn/o+vw11FPkWGJ8xhXJjZg1JJU/dT+ZrA0/0eOp2Qx5azGbf/4t3OVJKfFk+IeM+vxFiqxCXAy3pZ7I7AfSuPmcpkz5Zgtpw9N5dNoKdu09EO7yJMSiP/w1zl+kWKpVjOfhnkl8PiSFXq3rMW72Ojo//jljNXFcVPFk+IdunL+GeoqUVMPqFRlx5Rl8dFdn2iRW59HpKzl3eDpvL8zRxHFRwJPhrymdRbwjqf7xvHxTe16/uQM1Kycw5K3F9Bw1m1krt2t4aATzZPiHjIZ6ioTM2c39E8c93b8Ne/cf5MYJC+j//Fcs3vRzuEuTEoju8PfFwCG1/EVCxeczep9en8/uT+GRPqew5vvd9B39BXdMXMR3mjguopTZGr5hYT5Av5aKhFp8rI/rz27CJW0b8HzmOp6fvZ4Zy7ZxVYdE7j6vBbUqJ4S7RDmG6G75a5y/SKmqUiGO+y9oScYDqVx5ZiMmzttIyuOzePKz1ezRxHGe5snw16yeIpGlTpUK/N/Frfnkvi50Oak2T362hpRh6Zo4zsM8Gf4a7SMSmU7MP3FcrUr89YNlnD8igymLt3BIw0M9xZPhHzIKf5GwaJtYnTdv7cj4G5I5Li6Guyd9TZ/Rc8hcvUPDQz1C4S8ipcLMOPfkunx0d2dGXHE6P+05wHXj53P1C/M0PNQDFP4iUqpifMYlbRvy+ZAU/qdXEiu3/frH8FCtKxw+0R/+oGmdRTwgITaGm85pSsbQVO4+tzmzVm2n64gM/vr+Uq0rHAblJPzV+hfxirzhoXnrCr8+fyOpwzQ8tKyVafibWSUzyzKzXmV0QP9nhb+I5+RfVzj/8NDXvtqg4aFloEjhb2bjzWy7mS0tsL27ma0ys2wze7AIL/UnYHJJCi0RtfxFPC9vXeF3bjubprUq8pf3l9JtZCYfL92mkUGlqKgt/wlA9/wbzCwGGA30AJKA/maWZGatzWxqgY86ZtYVWA5sD2H9R6fwF4kY7RpXZ/KtZ/H8dcn4fMag1xZy2Zi5ZH33Y7hLi0pFmtvHOZdpZk0KbG4PZDvn1gGY2RtAX+fco8B/deuYWSpQCf8bxW9mNs25/05lMxsIDARITEws8okUSuEvElHMjK5JdUlrWZu3FuYw8tPVXDZmLhck1eWB7ifTvE7lcJcYNYKZ2K0BsCnf8xygw5F2ds49DGBmNwA7Cwv+wH7jzGwr0Ds+Pr5dEPUp/EUiVGyMj/7tE+l7Rn3Gz1nPmIx1dHsykyvPbMS957WgzvEVwl1ixCvz0T7OuQnOuanH2Cd00zuAwl8kQlWMj+XOc1uQMTSVazs2ZvKCTaQMS2fYjJX8sk/rCgcjmPDfDDTK97xhYFvQQjqxGyj8RSJczcoJ/K3PKcwcnELXpLqMnrWWLo/P4vnMdVpXuISCCf8FQAsza2pm8UA/YEooigp9y18jBkSiQeOalRjVvw1T7zqH0xpW4/+mrSBteDqTszZpXeFiKupQz0nAXKClmeWY2QDnXC5wJzADWAFMds4tC0VRavmLyNGc2qAqr9zUntdv6UCdKgk88Pa39Hgqk8+Wf6/hoUVkXv6HSk5OdllZWSV/gQUvwEeDYcgaqFwndIWJiGc455i+dBvDZqxi/c49nNmkOn/qfjLJTWqEu7SwMbOFzrnko+3jyekd1PIXkaIyMy5sXY9P7uvCPy86le9+2MtlY+Zy88sLWLXt13CX51meDH+N9hGR4oqL8XFNx8ZkDE1laLeWzFv/I92fymTIW4vZ8vNv4S7PczwZ/iGj8BcpdyrGx3JHWnMyh6Zx8zlNmfLNFtKGp/Pvj1ey6zcND83jyfBXt4+IBKt6pXge7pnEzMEpXNi6Hs+lryVl2CxenLOe33M1PNST4a9uHxEJlUY1KjLyyjOYetc5tG5QlX9MXa51hfFo+IeMwl9EAk5tUJVXB3TglZvaUzkhjrsnfU3f0V/wRfbOcJcWFp4M/9B3+5Tfd3cROVyXk2rz0V3nMOKK0/lxz36ufmEe142fz7ItQeZNhPFk+KvbR0RKky+wrvDMwSn8pWcrvs35mZ6j5nDvG1+z6ce94S6vTHgy/ENGK3mJyFFUiIvh5s7NyBiaxu2pJ/Lxsm2c+0Q6j3y4jB92R/e6wlEe/mr5i8ixVT0ujge6n0z6kDQua9eQl7/8jpRh6Tw9cw1790fnusKeDH8N9RSRcDihagUeveQ0PrmvC2efWJMnPl1NyrB0Js6LvnWFPRn+6vMXkXBqXqcK465L5p3bzqJxjYo8/N5SLhiZybQlW6Nm4jhPhn/IKPxFJAjtGtfgrUFn8cJ1ycT6jNsnLqLPM18wZ03kDw8tH+F/SH/NJyIlY2acn1SXj+/twvDL/cNDr3nRPzx0xdZfwl1eiUV3+BMY7UN0/JomIuET4zMua9eQz4f4h4cu3vQzF46azZC3FrN1V+RNHBfd4f9Ht094yxCR6JEQ6x8emjk0jVs6N2PKN1tIHZbO4x9H1rrCngz/0I320Th/ESkdVSvG8ecLWzFzcAo9Tj2BZ9PXkjosnQlfrGd/rvczx5PhH7rRPur2EZHS1ahGRZ7s519X+OQTqvC3D5dz3oh03sraRK6Hh4d6MvxDRy1/ESkbpzaoysSbOzDhxjOpdlw8Q9/+lgtGZvLx0m2eHB4a3eGvid1EpAyZGakt6zDlzk6MvbYdPp8x6LWFXD5mLgs3/BTu8g4T5eGvlr+IlD0zo9spJ/DxPZ3518Wt2fDjXi597ktue20h63fuCXd5AMSGu4BSpT/yEpEwio3xcVWHRPqeUZ8XZq9nbOZaPl3+PVd3SOTu81pQs3JC2GqL8pZ/3ump20dEwqdSQiz3nN+CjKFp9GvfiNfmbSRlWDrPfL6GPb+HZ+K4Mgt/M0s1s9lmNsbMUsvoqP5PavmLiAfUrpLAPy9qzYx7/RPHDf9kNZ0fn8W4zLXsO1C2MxEUKfzNbAKLZ2gAAAjrSURBVLyZbTezpQW2dzezVWaWbWYPHuNlHLAbqADklKzcYtINXxHxoOZ1KjPuumTeu/1sTql/PP+atpJzh6fz9sIcDpbRusJFbflPALrn32BmMcBooAeQBPQ3syQza21mUwt81AFmO+d6AH8CHgndKRyFbviKiIe1SazOqwM68PotHahVJYEhby2m56jZZKzeUerHLlL4O+cygR8LbG4PZDvn1jnn9gNvAH2dc0ucc70KfGx37o8E/gkom7sc6vMXkQhw9om1eP/2Tozq34Y9+3O5fvx8Zq74vlSPGcxonwbApnzPc4AOR9rZzC4BugHVgGeOst9AYCBAYmJiEOXBf/r8Ff4i4m0+n9Hn9Pp0O6UuH3y9hZSTapfq8cpsqKdz7l3g3SLsN87MtgK94+Pj2wV1UPX5i0iESYiN4YozG5X6cYIZ7bMZyF9hw8C2oIVubp+8F1Sfv4hIfsGE/wKghZk1NbN4oB8wJRRFhXwNX/X5i4gcpqhDPScBc4GWZpZjZgOcc7nAncAMYAUw2Tm3LBRFhazlr3H+IiKFKlKfv3Ou/xG2TwOmhbQi/C1/oHfz5s2DfCH1+YuIFMaT0zuErs9fc/uIiBTGk+GvlbxEREqXJ8M/5C1/3fAVETmMJ8M/dNTyFxEpjCfDP+RDPXXDV0TkMJ4M/5Av4K7wFxE5jCfDP2TU5y8iUihPhn/Iun3yqM9fROQwngz/0I/zV8tfRCQ/T4Z/yGicv4hIoaI8/NXnLyJSGE+Gf+iHeqrlLyKSnyfDX7N6ioiULk+Gf8johq+ISKGiPPzV8hcRKUyUh79u+IqIFCa6wx9N7yAiUhhPhr8mdhMRKV2eDP/QT+ymPn8Rkfw8Gf4hkxf+6vMXETlMdIe/xvmLiBQqusNfff4iIoWK8vBXy19EpDDRHf6+OEi6CGo2D3clIiKeEltWBzIzH/AP4Hggyzn3cqkfNK4CXFH6hxERiTRFavmb2Xgz225mSwts725mq8ws28wePMbL9AUaAgeAnJKVKyIioVDUlv8E4BnglbwNZhYDjAa64g/zBWY2BYgBHi3w/TcBLYEvnXNjzextYGZwpYuISEkVKfydc5lm1qTA5vZAtnNuHYCZvQH0dc49CvQq+BpmlgPsDzw9eKRjmdlAYCBAYmJiUcoTEZFiCuaGbwNgU77nOYFtR/Iu0M3MngYyj7STc26ccy7ZOZdcu3btIMoTEZEjKbMbvs65vcCAouxrZr2B3s2ba5SOiEhpCKblvxlolO95w8C2oIVuJS8RESlMMOG/AGhhZk3NLB7oB0wJRVEhm9VTREQKVdShnpOAuUBLM8sxswHOuVzgTmAGsAKY7JxbFoqi1PIXESld5jw4701enz9wJbCmhC9TC9gZsqLCL5rOJ5rOBaLrfKLpXCC6zqc459LYOXfUETOeDP9QMLMs51xyuOsIlWg6n2g6F4iu84mmc4HoOp9Qn0t0z+0jIiKFUviLiJRD0Rz+48JdQIhF0/lE07lAdJ1PNJ0LRNf5hPRcorbPX0REjiyaW/4iInIECn8RkXIoKsO/mOsMhJ2ZNTKzWWa23MyWmdk9ge01zOxTM1sT+Fw9sN3MbFTg/L41s7bhPYP/ZmYxZva1mU0NPG9qZvMCNb8Z+KtwzCwh8Dw78PUm4ay7MGZWzczeNrOVZrbCzM6K8GtzX+DnbKmZTTKzCpFyfQpbW6Qk18LMrg/sv8bMrg/HuQTqKOx8hgV+1r41s/fMrFq+rz0UOJ9VZtYt3/biZ55zLqo+8K8nsBZoBsQDi4GkcNd1jJrrAW0Dj6sAq4Ek4HHgwcD2B4F/Bx5fCEwHDOgIzAv3ORRyTvcDrwNTA88nA/0Cj8cAtwUe3w6MCTzuB7wZ7toLOZeXgZsDj+OBapF6bfDPvLseOC7fdbkhUq4P0AVoCyzNt61Y1wKoAawLfK4eeFzdQ+dzARAbePzvfOeTFMizBKBpIOdiSpp5Yf9hLIV/zLOAGfmePwQ8FO66inkOH+BfJGcVUC+wrR6wKvB4LNA/3/5/7OeFD/yT/M0EzgWmBv7z7cz3A/3HNcI/PchZgcexgf0s3OeQ71yqBsLSCmyP1GuTNxV7jcC/91SgWyRdH6BJgbAs1rUA+gNj820/bL9wn0+Br10MTAw8PizL8q5NSTMvGrt9irvOgKcEfq1uA8wD6jrntga+tA2oG3js9XN8EngAOBR4XhP42fnng4LD6/3jXAJf3xXY3yuaAjuAlwLdWC+YWSUi9No45zYDw4GNwFb8/94LidzrA8W/Fp6+RgXchP+3Fwjx+URj+EcsM6sMvAPc65z7Jf/XnP8t3fPjcs2sF7DdObcw3LWESCz+X8ufc861Afbg71r4Q6RcG4BAf3hf/G9q9YFKQPewFhVCkXQtjsXMHgZygYml8frRGP6lts5AaTKzOPzBP9E5925g8/dmVi/w9XrA9sB2L59jJ6CPmX0HvIG/6+cpoJqZ5S0elL/eP84l8PWqwA9lWfAx5AA5zrl5gedv438ziMRrA3A+sN45t8M5dwD/CnudiNzrA8W/Fl6/RpjZDfiXw7068IYGIT6faAz/UltnoLSYmQEvAiuccyPyfWkKkDcS4Xr89wLytl8XGM3QEdiV79fesHLOPeSca+ica4L/3/5z59zVwCzgssBuBc8l7xwvC+zvmZabc24bsMnMWgY2nQcsJwKvTcBGoKOZVQz83OWdT0Ren4DiXosZwAVmVj3wm9AFgW2eYGbd8Xeb9nH+FRDzTAH6BUZgNQVaAPMpaeaF88ZNKd5AuRD/iJm1wMPhrqcI9Z6D/1fVb4FvAh8X4u9bnYl/WuvPgBqB/Q0YHTi/JUByuM/hCOeVyn9G+zQL/KBmA28BCYHtFQLPswNfbxbuugs5jzOArMD1eR//CJGIvTbAI8BKYCnwKv7RIxFxfYBJ+O9VHMD/W9mAklwL/H3p2YGPGz12Ptn4+/DzsmBMvv0fDpzPKqBHvu3FzjxN7yAiUg5FY7ePiIgcg8JfRKQcUviLiJRDCn8RkXJI4S8iUg4p/EVEyiGFv4hIOfT/+WSqC3e0hlgAAAAASUVORK5CYII=\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "-haZiAbQdPoh"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}