mscroggs/lecture8.ipynb

## lecture8.ipynb
{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "authorship_tag": "ABX9TyO96Vsmzh9fGjdFjzuk30IF",
      "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/04dbdf64a91a623af72e59ecc1aa4379/lecture8.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "__zOb1U8SSUV",
        "outputId": "a18e3496-1370-4266-bde1-546d92edff40"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/\n",
            "Collecting petsc4py\n",
            "  Downloading petsc4py-3.18.2.tar.gz (2.5 MB)\n",
            "\u001b[K     |████████████████████████████████| 2.5 MB 6.5 MB/s \n",
            "\u001b[?25hCollecting pyamg\n",
            "  Downloading pyamg-4.2.3-cp38-cp38-manylinux_2_12_x86_64.manylinux2010_x86_64.whl (1.7 MB)\n",
            "\u001b[K     |████████████████████████████████| 1.7 MB 43.4 MB/s \n",
            "\u001b[?25hRequirement already satisfied: numpy in /usr/local/lib/python3.8/dist-packages (from petsc4py) (1.21.6)\n",
            "Collecting petsc<3.19,>=3.18\n",
            "  Downloading petsc-3.18.2.tar.gz (18.4 MB)\n",
            "\u001b[K     |████████████████████████████████| 18.4 MB 427 kB/s \n",
            "\u001b[?25hCollecting mpi4py>=1.2.2\n",
            "  Downloading mpi4py-3.1.4.tar.gz (2.5 MB)\n",
            "\u001b[K     |████████████████████████████████| 2.5 MB 39.3 MB/s \n",
            "\u001b[?25h  Installing build dependencies ... \u001b[?25l\u001b[?25hdone\n",
            "  Getting requirements to build wheel ... \u001b[?25l\u001b[?25hdone\n",
            "    Preparing wheel metadata ... \u001b[?25l\u001b[?25hdone\n",
            "Requirement already satisfied: scipy>=0.12.0 in /usr/local/lib/python3.8/dist-packages (from pyamg) (1.7.3)\n",
            "Building wheels for collected packages: petsc4py, petsc, mpi4py\n",
            "  Building wheel for petsc4py (setup.py) ... \u001b[?25lerror\n",
            "\u001b[31m  ERROR: Failed building wheel for petsc4py\u001b[0m\n",
            "\u001b[?25h  Running setup.py clean for petsc4py\n",
            "  Building wheel for petsc (setup.py) ... \u001b[?25lerror\n",
            "\u001b[31m  ERROR: Failed building wheel for petsc\u001b[0m\n",
            "\u001b[?25h  Running setup.py clean for petsc\n",
            "  Building wheel for mpi4py (PEP 517) ... \u001b[?25l\u001b[?25hdone\n",
            "  Created wheel for mpi4py: filename=mpi4py-3.1.4-cp38-cp38-linux_x86_64.whl size=4438508 sha256=b512cced5ddd7e790159923a28c096bcd8552b3108c996aedb6d0db6f032a92f\n",
            "  Stored in directory: /root/.cache/pip/wheels/f3/35/48/0b9a7076995eea5ea64a7e4bc3f0f342f453080795276264e7\n",
            "Successfully built mpi4py\n",
            "Failed to build petsc4py petsc\n",
            "Installing collected packages: mpi4py, petsc, pyamg, petsc4py\n",
            "    Running setup.py install for petsc ... \u001b[?25l\u001b[?25hdone\n",
            "\u001b[33m  DEPRECATION: petsc was installed using the legacy 'setup.py install' method, because a wheel could not be built for it. A possible replacement is to fix the wheel build issue reported above. You can find discussion regarding this at https://github.com/pypa/pip/issues/8368.\u001b[0m\n",
            "    Running setup.py install for petsc4py ... \u001b[?25l\u001b[?25hdone\n",
            "\u001b[33m  DEPRECATION: petsc4py was installed using the legacy 'setup.py install' method, because a wheel could not be built for it. A possible replacement is to fix the wheel build issue reported above. You can find discussion regarding this at https://github.com/pypa/pip/issues/8368.\u001b[0m\n",
            "Successfully installed mpi4py-3.1.4 petsc-3.18.2 petsc4py-3.18.2 pyamg-4.2.3\n"
          ]
        }
      ],
      "source": [
        "!pip install petsc4py pyamg"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "from petsc4py import PETSc\n",
        "import numpy as np"
      ],
      "metadata": {
        "id": "8cM8BeuWSaQu"
      },
      "execution_count": 2,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "n = 1000\n",
        "\n",
        "nnz = np.full(n, 3, dtype=np.int32)\n",
        "nnz[0] = 2\n",
        "nnz[-1] = 2\n",
        "\n",
        "A = PETSc.Mat()\n",
        "A.createAIJ([n, n], nnz=nnz)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "KdEKW8wgUn1Q",
        "outputId": "cc23da30-ef57-45b1-be92-dc5270d9e9ba"
      },
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<petsc4py.PETSc.Mat at 0x7f2d75e87c70>"
            ]
          },
          "metadata": {},
          "execution_count": 6
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "for i in range(n):\n",
        "  A.setValue(i, i, 3)\n",
        "for i in range(n - 1):\n",
        "  A.setValue(i, i + 1, -1)\n",
        "  A.setValue(i + 1, i, -1)"
      ],
      "metadata": {
        "id": "BRCvxR7pU8PZ"
      },
      "execution_count": 8,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "A.assemble()"
      ],
      "metadata": {
        "id": "l8brgMSqXJvn"
      },
      "execution_count": 9,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "A[:,:]"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "9xvtnqUlXw9u",
        "outputId": "45929684-f1c2-4e10-ba0d-ac01c926025f"
      },
      "execution_count": 12,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([[ 3., -1.,  0., ...,  0.,  0.,  0.],\n",
              "       [-1.,  3., -1., ...,  0.,  0.,  0.],\n",
              "       [ 0., -1.,  3., ...,  0.,  0.,  0.],\n",
              "       ...,\n",
              "       [ 0.,  0.,  0., ...,  3., -1.,  0.],\n",
              "       [ 0.,  0.,  0., ..., -1.,  3., -1.],\n",
              "       [ 0.,  0.,  0., ...,  0., -1.,  3.]])"
            ]
          },
          "metadata": {},
          "execution_count": 12
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "A.size"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "sH1T5AJvY5JB",
        "outputId": "580af5f2-a33b-4dc8-dc70-0cd7b38166c6"
      },
      "execution_count": 13,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(1000, 1000)"
            ]
          },
          "metadata": {},
          "execution_count": 13
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "A.local_size"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "WvIrOkx6ZKAu",
        "outputId": "dffe20fa-2a4a-4fb5-d9b2-980265e00ca3"
      },
      "execution_count": 14,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(1000, 1000)"
            ]
          },
          "metadata": {},
          "execution_count": 14
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "A.isSymmetric()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "2mcsRxy8ZNh_",
        "outputId": "dd0c48b3-a669-42d8-8ec5-26c2d995d8ab"
      },
      "execution_count": 15,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "True"
            ]
          },
          "metadata": {},
          "execution_count": 15
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "A.getInfo()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "qPdfSWx2ZXUe",
        "outputId": "5b6b4963-81aa-4424-d9cd-001e49302368"
      },
      "execution_count": 16,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "{'block_size': 1.0,\n",
              " 'nz_allocated': 2998.0,\n",
              " 'nz_used': 2998.0,\n",
              " 'nz_unneeded': 0.0,\n",
              " 'memory': 0.0,\n",
              " 'assemblies': 1.0,\n",
              " 'mallocs': 0.0,\n",
              " 'fill_ratio_given': 0.0,\n",
              " 'fill_ratio_needed': 0.0,\n",
              " 'factor_mallocs': 0.0}"
            ]
          },
          "metadata": {},
          "execution_count": 16
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "5j_kWJ9lZbrK"
      },
      "execution_count": 18,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "p94SIVqNZpLa"
      },
      "execution_count": 19,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "b = A.createVecLeft()\n",
        "b.array[:] = 1.0"
      ],
      "metadata": {
        "id": "DULhwmsHZ7EZ"
      },
      "execution_count": 20,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "x = A.createVecRight()\n",
        "\n",
        "ksp = PETSc.KSP().create()\n",
        "ksp.setOperators(A)\n",
        "ksp.setType(\"bcgs\")\n",
        "ksp.setConvergenceHistory()\n",
        "ksp.getPC().setType(\"none\")\n",
        "ksp.solve(b, x)"
      ],
      "metadata": {
        "id": "O0JAivySaD6I"
      },
      "execution_count": 27,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "from matplotlib import pylab as plt\n",
        "\n",
        "res = ksp.getConvergenceHistory()\n",
        "plt.plot(res)\n",
        "plt.yscale(\"log\")"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 265
        },
        "id": "VIUl6N9raMVd",
        "outputId": "0f8f57c8-2344-470d-9ab5-26c20df6248c"
      },
      "execution_count": 28,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXhV9bn28e+TmTHIjCQhJAGU0SFOgMwoVtFaZ2tr1YrUEX1rq6fvafWc03qsbwW0TqhUbdWqVC0OLQVkEAUFrCCIZCIxASQQIBAgZPq9f7DTpimBwN7JWnvv+3NdXLKWe69166V3fjxr7bXNOYeIiES+GK8DiIhI61Dhi4hECRW+iEiUUOGLiEQJFb6ISJSI8zrAkXTt2tWlp6d7HUNEJKysXr16h3OuW+P9vi789PR0Vq1a5XUMEZGwYmZFh9uvkY6ISJTwZeGb2WQzm1VeXu51FBGRiOHLwnfOveOcm5KcnOx1FBGRiOHLwhcRkdBT4YuIRAlfFr5m+CIioefLwtcMX0Qk9HxZ+MFanl/Gcx8WeB1DRMRXIrLw567ZzK/e38Dnxbu9jiIi4hu+LPxgZ/j3f+tkenRM4t431nCwpjbE6UREwpMvCz/YGX7HpHh+9Z0h5JZW8PjCvBCnExEJT74s/FAYO6A7l52WwlNL8lm3WXf7iIhEbOED/PyigXRpl8CP31hDVU2d13FERDwV0YWf3DaeX146hK++2cuTizXaEZHoFtGFDzBxYA8uOeVEfvtBHhu27vE6joiIZ3xZ+KH+pO0DkwfRqW08985ZQ3WtRjsiEp18Wfih/qTtCe0S+O9LBrNu8x5mLdUHskQkOvmy8FvCBUN6ceGQXsxckEvOtr1exxERaXVRU/gAD14yiPZJcdw7Zy01Gu2ISJSJqsLv2j6RBy4exJri3Ty/bJPXcUREWlVUFT7A5KG9OG9gD34zP4f87RVexxERaTVRV/hmxv9cOpg28bH8ZM5aauuc15FERFpF1BU+QPcOSfxi8kBWF+3ihY8LvY4jItIqWq3wzSzDzJ43szmtdc4jufTU3ow7qTuPzPuKwh37vI4jItLigip8M5ttZqVmtq7R/klmttHM8szsPgDnXIFz7qZgzhdKZsavLh1CfGwMP/nTWuo02hGRCBfsCv8FYFLDHWYWCzwBXAAMBK4xs4FBnqdF9ExO4j8vHMinm3byh0+KvI4jItKigip859xSYGej3WcCeYEVfRXwR+CS5h7TzKaY2SozW7V9+/Zg4jXLFdkpjOrfjf/9y1cU79zf4ucTEfFKS8zwewPFDbZLgN5m1sXMngZONbP7m3qzc26Wcy7bOZfdrVu3Foj3r8yMh74zhBgz7ntzLc5ptCMikanVLto658qcc1Odc5nOuYeO9NpQPzztaHp3asP93zqJj/LKePXT4qO/QUQkDLVE4W8GUhtspwT2+dq1Z6YxPLMLv3p/A5t3H/A6johIyLVE4a8E+plZXzNLAK4G5h7LAUL9tMzmMDMevmwodc5x/5tfaLQjIhEn2NsyXwWWAwPMrMTMbnLO1QC3A/OADcDrzrn1wUdteamd2/LTSSexNGc7b6wu8TqOiEhImR9XsmY2GZiclZV1c25ubqueu67OcfWzK9iwdQ/z7x5Nz+SkVj2/iEiwzGy1cy678X5fPlrBi5FOvZgY49eXDaW6to6fvaXRjohEDl8WfmvfpdNYetd2/Pi8ASz8qpS3P/f99WYRkWbxZeF7ucKvd8OIvpyW1okH5n5J6d5Kz3KIiISKLwvfD2JjjF9fPowD1bX859vrNNoRkbDny8L3eqRTL6t7e+6Z2J9567fx7tqtnmYREQmWLwvfDyOdej8c2ZdhKcn8Yu56yioOeh1HROS4+bLw/SQuNoZHrhhGRWUNP58bFh8nEBE5LF8Wvl9GOvX69+jAneOzeG/tVv66TqMdEQlPvix8P4106t0yOpNBJ3bk/769jl37qryOIyJyzHxZ+H4UHxvDI5cPY/f+ah58R6MdEQk/KvxjMPDEjtw2Nou3P9/Cgi+3eR1HROSYqPCP0W1jszipZwf+460vKN9f7XUcEZFm82Xh++2ibUMJcTH8vyuGUbaviv9+70uv44iINJsvC9+PF20bGtw7mR+NzmTO6hIWbSz1Oo6ISLP4svDDwR3js+jfoz3/8eYX7KnUaEdE/E+Ff5wS42J55PJhbNtTyUPvb/A6jojIUanwgzAstRM3j8rg1U+LWZa7w+s4IiJHpMIP0t0T+pPRrR0//dNaKg7WeB1HRKRJvix8P9+l01hSfCyPXD6ULeUHePgvX3kdR0SkSb4sfL/fpdPY6X06c+OIvvx+RRHL88u8jiMicli+LPxw9OPzBpDepS0//dNa9ldptCMi/qPCD5E2CbE8fNlQvt65n0fmbfQ6jojIv1Hhh9BZGV24/pw+vPBxIasKd3odR0TkX6jwQ+wnk04i5YQ2/GTOWiqra72OIyLyDyr8EGuXGMfD3xlKwY59PDo/x+s4IiL/0GqFb2btzOxFM3vWzL7bWuf1wvCsrlx7VhrPfVjA37/e5XUcEREgyMI3s9lmVmpm6xrtn2RmG80sz8zuC+z+DjDHOXczcHEw5w0H919wEj07JnGvRjsi4hPBrvBfACY13GFmscATwAXAQOAaMxsIpADFgZdFfAN2SIrnocuGkldawWMLc72OIyISXOE755YCjW9HORPIc84VOOeqgD8ClwAlHCr9oM8bLkb378aV2Sk8s7SAL0r8/6lhEYlsLVG8vfnnSh4OFX1v4E3gMjN7CninqTeb2RQzW2Vmq7Zv394C8VrXzy4cSNf2Cdw7Zw1VNXVexxGRKNZqK23n3D7n3A3OuR85514+wutmAQ8CnyUkJLRWvBaT3Caeh74zhK++2ctvF+V5HUdEolhLFP5mILXBdkpgX7OF27N0jmbcST34zqm9eXJRHuu3aLQjIt5oicJfCfQzs75mlgBcDcw9lgOE09Mym+vnkwdyQrsE7n1jLdW1Gu2ISOsL9rbMV4HlwAAzKzGzm5xzNcDtwDxgA/C6c279sRw30lb4AJ3aJvA/3x7Ml1v38PTifK/jiEgUigvmzc65a5rY/z7w/vEe18wmA5OzsrKO9xC+dP6gnkwediKPfZDLeYN6MqBnB68jiUgU8eXtkZG4wq/34MWD6JgUz71z1lCj0Y6ItCJfFn4k69wugf+6ZDBrS8p59sNNXscRkSjiy8KPxIu2DV04tBcXDO7J9AU55JVWeB1HRKKELws/kkc69f7rksG0S4jl3jlrqK1zXscRkSjgy8KP9BU+QLcOiTxw8SD+/vVufveRRjsi0vJ8WfjRsMIHuHjYiUw4uQePzNvIB19twzmt9EWk5fiy8KOFmfGrSwfTMzmJG19YxWVPfcyy3B0qfhFpESp8j3XvmMT8u0fzy0sHs7W8kuue/4SrZq1gRUGZ19FEJMKYH1eTDT54dXNubvQ8S76yupbXVhbzxKI8SvceZERWF+6Z2J/T+3T2OpqIhBEzW+2cy/63/X4s/HrZ2dlu1apVXsdodZXVtfxhRRFPL8lnR0UVo/p3456J/TkltZPX0UQkDKjww9D+qhpeWl7EM0vy2bW/mvEndefuif0Z3DuyL2aLSHBU+GGs4mANL35cyKylBZQfqOb8QT2YNqE/J/fq6HU0EfGhsCr8aJ3hH82eympmL9vE8x9uYu/BGi4c0otpE/rRr4cewiYi/xRWhV9PK/zDK99fzXPLCpi9bBP7q2u5eNiJ3DW+Hxnd2nsdTUR8QIUfgXbuq2LW0gJe/LiQgzW1XHpqCneOz6JPl3ZeRxMRD6nwI9j2vQd5Zkk+v19RRE2d44rTU7h9XBYpJ7T1OpqIeECFHwVK91Ty5OJ8XvnkaxyOK7NTuX1cFr2S23gdTURakQo/imwtP8ATi/J4bWUxhnHtWWncOiaT7h2TvI4mIq1AhR+Finfu54lFebyxuoS4GON7Z/dh6phMurZP9DqaiLSgsCp83ZYZWkVl+3hsYR5v/b2ExLhYrh+ezpRRGXRul+B1NBFpAWFV+PW0wg+tgu0VzFyYy9w1W2gbH8uNI/vyw5EZJLeN9zqaiISQCl/+IXfbXmYszOW9tVvpkBjHTef25caRfemYpOIXiQQqfPk3G7buYcaCHOat30Zym3imjMrg+uHptE+M8zqaiARBhS9NWre5nOnzc1j4VSmd2yVwy6gMvndOH9omqPhFwpEKX47q8+LdTJ+fw5Kc7XRtn8DU0Zlcd3YfkuJjvY4mIsdAhS/NtqpwJ9MX5PBRXhndOyRy29gsrj4zlcQ4Fb9IOPC88M0sA/gZkOycu7w571Hhe2tFQRmP/i2HTwt30is5idvHZXHF6akkxOmbMUX8rKnCb9b/uWY228xKzWxdo/2TzGyjmeWZ2X1HOoZzrsA5d9OxxRYvnZ3RhdduOZuXf3gWvZKT+Nlb6xj3m8W8vrKY6to6r+OJyDFq1grfzEYBFcBLzrnBgX2xQA4wESgBVgLXALHAQ40OcaNzrjTwvjla4Ycf5xxLcrYzfX4Oa0rK6dOlLXeO68clp5xIXKxW/CJ+EvRIx8zSgXcbFP45wAPOufMD2/cDOOcal33j4xyx8M1sCjAFIC0t7fSioqJm5ZPW4Zzjg69KeXR+Duu37CGjWzvuGt+Pi4aeSGyMeR1PRAhypNOE3kBxg+2SwL6mAnQxs6eBU+t/OByOc26Wcy7bOZfdrVu3IOJJSzAzxp/cg3fvGMnT151OQmwMd/3xcybNWMp7a7dSV+ffmwBEol2r3WjtnCsDpjbntQ2epdOyoeS4mRmTBvfkvIE9+Mu6b5i+IIfbXvmMk3p2YNqE/pw/qAdmWvGL+EkwK/zNQGqD7ZTAvqA5595xzk1JTk4OxeGkBcXEGBcO7cW8aaOYefUpVNXUMfUPq7no8WUs3LANP9/2KxJtgin8lUA/M+trZgnA1cDcUIQys8lmNqu8vDwUh5NWEBtjXHJKb/529yh+c8Uw9lbWcNOLq/j2Ex+xeGOpil/EB5p7l86rwBigK7AN+IVz7nkz+xYwg0N35sx2zv0ylOF0l074qq6t463PNvPYB7mU7DrAaWmduGfiAEZkddGoR6SFef7Bq2Oh5+FHjqqaOuasLuG3H+SypbySM/t25p6J/Tk7o4vX0UQiVlgVfj2t8CPHwZpaXltZzG8/yKN070GGZ3bhnon9yU7v7HU0kYgTVoWvFX7kqqyu5ZVPvubJxfnsqDjIuf26cs/E/pyadoLX0UQiRlgVfj2t8CPXgapafr+ikKeXFLBzXxXjTurO3RP6MyRFd2aJBEuFL76072ANLy4vZNbSAnbvr2biwB5Mm9CPQSeq+EWOV1gVvkY60WdvZTUvfFTIsx8WsKeyhgsG92TahP4M6NnB62giYSesCr+eVvjRp/xANc8v28TsZZvYV1XDhUN6MW1Cf7K6t/c6mkjYUOFLWNm9v4pnPyzgdx8VUlldyyWn9ObO8f3o27Wd19FEfC+sCl8jHalXVnGQWR8W8NLHRVTV1nHpqb25c1w/0rq09TqaiG+FVeHX0wpf6m3fe5Cnl+TzhxVF1NY5rshO4baxWaScoOIXaUyFLxFh255KnlqczyuffI3DcdUZqdw2NoteyW28jibiGyp8iShbdh/gycV5vLayGMO49qw0bh2TSfeOSV5HE/GcCl8iUvHO/TyxKI83VpcQF2Ncd3Yfpo7OpFuHRK+jiXgmrApfF23lWBWV7ePxD/J487MSEuNi+f7wPtwyKpPO7RK8jibS6sKq8OtphS/HqmB7BY8tzOXPa7bQNj6WH4xI5+ZzM+jUVsUv0UOFL1Elr3QvMxbk8t4XW2mfEMcNI/ty08i+JLeJ9zqaSItT4UtU2vjNXmYuzOH9L76hQ1IcN5+bwQ0j0umQpOKXyKXCl6j25ZY9zFiQw9++3EZym3imjMrg+uHptE+M8zqaSMip8EWAL0rKmbEgh4VflXJC23huGZ3J98/pQ9sEFb9EjrAqfN2lIy3t8+LdTJ+fw5Kc7XRpl8DU0Zlcd3Yf2iTEeh1NJGhhVfj1tMKXlra6aBczFuTwYe4OurZP5NYxmVx7VhpJ8Sp+CV8qfJEj+HTTTqbPz2F5QRk9OiZy29gsrjojlcQ4Fb+EHxW+SDMszy9j+vwcPi3cSa/kJG4bm8WV2akkxMV4HU2k2VT4Is3knOPj/DIenZ/D6qJd9O7UhjvGZXHZ6SnEx6r4xf9U+CLHyDnH0twdTJ+fw+fFu0nt3IY7x/Xj0lN7E6fiFx9T4YscJ+ccizdu59H5OXyxuZz0Lm25c3w/LjmlN7Ex5nU8kX/TVOFrmSJyFGbG2JO6M/f2ETz7/WzaJMRxz+trmDh9CX/+fDO1df5dNIk01KqFb2bfNrNnzew1MzuvNc8tEiwzY+LAHrx3x0ievu404mNiuOuPn3P+jKW8u3YLdSp+8blmF76ZzTazUjNb12j/JDPbaGZ5ZnbfkY7hnHvbOXczMBW46vgii3grJsaYNLgXf7nrXJ649jQMuP2Vv3PBzA/567qtKn7xrWbP8M1sFFABvOScGxzYFwvkABOBEmAlcA0QCzzU6BA3OudKA+/7DfCyc+6zI51TM3wJB7V1jnfXbmHmwlwKtu9jYK+O3D2xPxNO7o6ZZvzS+kJy0dbM0oF3GxT+OcADzrnzA9v3AzjnGpd9/fsN+F9gvnNuQROvmQJMAUhLSzu9qKio2flEvFRb55i7ZjMzF+RSWLafIb2TuXtiP8YOUPFL62qpi7a9geIG2yWBfU25A5gAXG5mUw/3AufcLOBB4LOEBH1phYSP2Bjj0lNTWHDPaB65fCi7D1Rx4wur+PaTH7MkZzt+viNOokOrXrR1zj3mnDvdOTfVOff0EV73jnNuSnJycmvGEwmJuNgYrshO5YP/M4aHLxvCjr0HuX72p1z+9HKW5e5Q8Ytngi38zUBqg+2UwL6gmNlkM5tVXl4e7KFEPBMfG8NVZ6Sx6Mdj+OWlg9my+wDXPf8JVz2zguX5ZV7HkygU7Aw/jkMXbcdzqOhXAtc659aHIpwu2kokOVhTy2sri3liUR7b9hzknIwu3D2xP2f27ex1NIkwQc/wzexVYDkwwMxKzOwm51wNcDswD9gAvB6KstcKXyJRYlws3z8nnSX3juXnFw0kt7SCK59ZznXPfcLqol1ex5MooEcriHjkQFUtL39SxFOL8ynbV8Xo/t24e2J/Tknt5HU0CXN6lo6IT+2vquGl5UU8sySfXfurGXdSd6ZN6MfQFBW/HJ+wKnx9xaFEo4qDNbz4cSHPfljA7v3VTDi5O9Mm9Gdwb92tJscmrAq/nlb4Eo32VlbzwkeHin9PZQ0TB/Zg2oR+DDpRxS/No8IXCTN7Kqv53bJCnltWwN7KGiYN6sm0if04qWdHr6OJz4VV4WukI/JP5Qeqmb1sE7OXbWLvwRq+NaQnd43vz4CeHbyOJj4VVoVfTyt8kX8q31/Nc8sK+N1HheyrquFbQ3oxbXw/+vVQ8cu/UuGLRIhd+6p4blkBL3xUyP7qWiYPPZE7x/cjq3t7r6OJT4RV4WukI3J0O/dVMWtpAS8tL6SyupZLTunNHeOyyOim4o92YVX49bTCFzm6soqDgeIv4mBNLd8+tTd3jutHetd2XkcTj6jwRSLcjoqDPLMkn9+vKKK61nFpoPjTurT1Opq0MhW+SJQo3VvJ04sLePmTImrqHJeflsLt47JI7azijxZhVfia4YsEr3RPJU8uzueVT7+mrs5xRXYKt43NIuUEFX+kC6vCr6cVvkjwvimv5KnFebz6aTEOxxXZqdw2Novendp4HU1aiApfJMpt2X2AJxfn8drKQ99KevUZadw6NpNeySr+SKPCFxEANu8+wBOL8nhjVTGGcc2Zqdw6NoseHZO8jiYhosIXkX9Rsmt/oPhLiIkxvntWGj8anUl3FX/YU+GLyGEV79zP4x/k8qfPNhMXY1x3dh+mjs6kW4dEr6PJcVLhi8gRFZXt4/EP8njzsxIS4mL43tl9uGV0Jl3bq/jDTVgVvm7LFPHOph37eHxhLm9/vvnQ9/AO78MtozLp3C7B62jSTGFV+PW0whfxTv72Ch5fmMuf12yhTXws1w9PZ8q5GZyg4vc9Fb6IHJe80goeW5jLO2u30DY+lh+MSOfmczPo1FbF71cqfBEJSs62vcxcmMt7a7fSPjGOG0ekc9PIDJLbxnsdTRpR4YtISGz8Zi8zF+bw/hff0CExjhtG9uWmkX1JbqPi9wsVvoiE1Iate5i5IJe/rv+GDklx3DiiLzeq+H1BhS8iLeLLLXuYuTCHeeu30SEpjptGHir+jkkqfq94XvhmdjJwF9AVWOice+po71Hhi4SP9VvKmbkgl799uY2OSXH88NwMfjAiXcXvgaYKP6aZb55tZqVmtq7R/klmttHM8szsviMdwzm3wTk3FbgSGHEs4UXE/wadmMys72fz7h0jOSujC4/Oz+Hchxfx+MJc9lZWex1PaOYK38xGARXAS865wYF9sUAOMBEoAVYC1wCxwEONDnGjc67UzC4GfgT83jn3ytHOqxW+SPhat7mcGQtyWbBhG8lt4rn53L5cPzydDlrxt7igRzpmlg6826DwzwEecM6dH9i+H8A517jsD3es95xzFzbx96YAUwDS0tJOLyoqalY+EfGnL0rKmbkwhwUbSunUNp6bz83g+uHptE+M8zpaxApqpNOE3kBxg+2SwL6mAowxs8fM7Bng/aZe55yb5ZzLds5ld+vWLYh4IuIHQ1KSee76M5h7+whOTzuBR+ZtZOTDH/DEojwqDtZ4HS+qtNqPWOfcYmBxc17b4Fk6LRlJRFrR0JROPP+DM1hTvJuZC3N5ZN5GnvuwgJtHZXD9Oem004q/xQWzwt8MpDbYTgnsC5pz7h3n3JTk5ORQHE5EfGRYaidm/+AM3r5tBKekduLXfz204n9qcT77tOJvUcEU/kqgn5n1NbME4GpgbihCmdlkM5tVXl4eisOJiA+dktqJ391wJm/dOpxhqZ14+K9fce6vF/H0knz2V6n4W0Jz79J5FRjDoXvotwG/cM49b2bfAmZw6M6c2c65X4YynO7SEYken329i5kLclmSs50u7RKYMiqD753Th7YJGvUcK88/eHUs9Dx8kei1umgXMxfmsjRQ/LeMzuC6s1X8xyKsCr+eVvgi0Wt10U5mLMjlw9wddG2fwC2jMrnu7D60SYj1OprvhVXha4UvIvVWFR4q/mV5h4p/6uhMvnuWiv9Iwqrw62mFLyL1VhbuZMaCHD7KK6Nr+0SmBkY9SfEq/sZU+CISET7ddKj4P84vo1uHRH40OpNrz0pT8TcQVoWvkY6IHM0nBWXMWJDL8oIyundI5EdjMrnmTBU/hFnh19MKX0SOZkVBGTMW5LCiYCfdOyRy65hMro7y4lfhi0hEW55fxvQFOXy6aSc9OiZy65gsrjojNSqLX4UvIlHh4/wdzJify6eFO+nZMYlbx2ZyZXZ0FX9YFb5m+CISDOfcP1b8Kwt30bNjEreNzeTKM1JJjIv84g+rwq+nFb6IBMM5x8f5ZUyfn8Oqol30Sk7i1jGRX/wqfBGJWtFW/Cp8EYl60VL8YVX4muGLSEuK9OIPq8KvpxW+iLSkSC1+Fb6ISBMirfhV+CIiRxEpxa/CFxFppsbFH2738avwRUSOUX3xzwizD3CFVeHrLh0R8ZNwK/6wKvx6WuGLiJ+ES/Gr8EVEQsTvxa/CFxEJMb8WvwpfRKSF+K34VfgiIi3scMXvxfP4VfgiIq3kcM/jb83iV+GLiLQyr4q/qcKPabEzHj5EOzNbZWYXteZ5RUS8YGYMz+rK67ecwys/PIvUzm34+Z/XM+aRxby0vJDK6tpWzdOswjez2WZWambrGu2fZGYbzSzPzO5rxqF+Crx+PEFFRMKVX4q/WSMdMxsFVAAvOecGB/bFAjnARKAEWAlcA8QCDzU6xI3AMKALkATscM69e7TzaqQjIpGopUc9Qc/wzSwdeLdB4Z8DPOCcOz+wfX/gH6Rx2de//5dAO2AgcAC41DlXd5jXTQGmAKSlpZ1eVFTUrHwiIuGmpYq/JWb4vYHiBtslgX2H5Zz7mXNuGvAK8Ozhyj7wulnOuWznXHa3bt2CiCci4m9HGvWsKtwZ8vPFhfyIR+Gce+For2nw8LSWDyQi4rH64j8nswvL88t4akk+6V3bhfw8wRT+ZiC1wXZKYJ+IiByH+uIfntW1RY4fzEhnJdDPzPqaWQJwNTA3FKGcc+8456YkJyeH4nAiIkLzb8t8FVgODDCzEjO7yTlXA9wOzAM2AK8759a3XFQREQlGs0Y6zrlrmtj/PvB+SBOhGb6ISEto1U/aNpdGOiIioefLwjezyWY2q7y83OsoIiIRw5eFrxW+iEjo+bLwRUQk9HxZ+BrpiIiEnq+fh29m24HjfZhOV2BHCOO0tHDKG05ZIbzyhlNWCK+80ZS1j3Pu355N4+vCD4aZrTrcw4P8KpzyhlNWCK+84ZQVwiuvsvp0pCMiIqGnwhcRiRKRXPizvA5wjMIpbzhlhfDKG05ZIbzyRn3WiJ3hi4jIv4rkFb6IiDSgwhcRiRIRWfhmNsnMNppZnpnd53WeIzGz2WZWambrvM5yNGaWamaLzOxLM1tvZnd5nakpZpZkZp+a2ZpA1ge9znQ0ZhZrZn83s3e9znI0ZlZoZl+Y2edmtsrrPEdjZp3MbI6ZfWVmGwLfye07ZjYg8O+0/tceM5sWsuNH2gzfzGKBHGAih75ndyVwjXPuS0+DNcHMRgEVwEv1XxDvV2bWC+jlnPvMzDoAq4Fv+/HfrZkZ0M45V2Fm8cAy4C7n3AqPozXJzO4BsoGOzrmLvM5zJGZWCGQ758Lig0xm9iLwoXPuucAXNrV1zu32OteRBLpsM3CWc+54P4D6LyJxhX8mkOecK3DOVQF/BC7xOFOTnHNLgdB/W3ELcM5tdc59Fvj9Xg598U2TX1zvJXdIRWAzPvDLt6sbM0sBLgSe8zpLpDGzZGAU8DyAc67K72UfMB7ID1XZQ2QWfm+guMF2CT4tpXBmZunAqcAn3iZpWmBE8jlQCsx3zvk2KzAD+AlQ53WQZnLA31NbfeoAAAHUSURBVMxstZlN8TrMUfQFtgO/C4zMnjOz0H9DeOhdDbwaygNGYuFLCzOz9sCfgGnOuT1e52mKc67WOXcKkAKcaWa+HJmZ2UVAqXNutddZjsFI59xpwAXAbYHRpF/FAacBTznnTgX2AX6/tpcAXAy8EcrjRmLhbwZSG2ynBPZJCATm4X8CXnbOvel1nuYI/PF9ETDJ6yxNGAFcHJiL/xEYZ2Z/8DbSkTnnNgf+Wgq8xaFRql+VACUN/oQ3h0M/APzsAuAz59y2UB40Egt/JdDPzPoGfkpeDcz1OFNECFwIfR7Y4Jx71Os8R2Jm3cysU+D3bTh0Ef8rb1MdnnPufudcinMunUP/vX7gnLvO41hNMrN2gYv2BEYj5wG+vcvMOfcNUGxmAwK7xgO+u9GgkWsI8TgHmvkl5uHEOVdjZrcD84BYYLZzbr3HsZpkZq8CY4CuZlYC/MI597y3qZo0Avge8EVgNg7wH4Evs/ebXsCLgTsdYoDXnXO+v90xTPQA3jr085844BXn3F+9jXRUdwAvBxaBBcANHudpUuCH6ETglpAfO9JuyxQRkcOLxJGOiIgchgpfRCRKqPBFRKKECl9EJEqo8EVEooQKX0QkSqjwRUSixP8HAaNqMDRyN4sAAAAASUVORK5CYII=\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "x = A.createVecRight()\n",
        "\n",
        "ksp = PETSc.KSP().create()\n",
        "ksp.setOperators(A)\n",
        "ksp.setType(\"preonly\")  # Direct solver\n",
        "ksp.setConvergenceHistory()\n",
        "ksp.getPC().setType(\"lu\")\n",
        "ksp.solve(b, x)\n",
        "\n",
        "res2 = ksp.getConvergenceHistory()\n",
        "plt.plot(res)\n",
        "plt.plot(res2)\n",
        "plt.yscale(\"log\")"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 265
        },
        "id": "dRiZbuapavea",
        "outputId": "b8514f06-161a-4d5b-e144-2b3f2e4c0f1e"
      },
      "execution_count": 34,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD4CAYAAAAKA1qZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de3Sc9X3n8fdXd8u6GMuSLWkkyzdsyxgwHommlABNyJqDF7K9BXO2zTa0XtqQJmnTLuylOTltD+k2TTctNMQFh5IGE+LQxgeckAsxkMQJlg0BW7ZBcQBLFpYvwfIFX2R9949nNBoJSZY1o3lmNJ/XOTrj+Wnmma99QB89v6u5OyIikpvywi5ARETCoxAQEclhCgERkRymEBARyWEKARGRHFYQdgEXY9asWd7U1BR2GSIiWWX79u2H3b16pO9lVQg0NTXR1tYWdhkiIlnFzN4Y7XvqDhIRyWEKARGRHKYQEBHJYQoBEZEcphAQEclhCgERkRymEBARyWE5EQInu1/l+Qc+xqG39oddiohIRsmJENjz831c+9YjfOafH+b+H3Rw+tz5sEsSEckIORECK6++Ds8r5OaZXfzd03t5398/y6afHUAH6ohIrsuJEKBwGjZnOTdV7ufRP7yaimmF/MmGF/mtB7by0v63w65ORCQ0uRECAJEWOLCDX22awZMf+zX+9jeX88aRU3zw/h/xya+9RPexd8KuUEQk7UINATObb2YPmdnGSf+whlY4dwp62snPMz7U0siWP7+eP75+AU+90s0Nn9vC57/7KqfO9k16KSIimSLlIWBm682sx8x2DmtfZWZ7zazDzO4GcPd97n5HqmsYUSQaPHa+EG8qKy7gL1Yt4ft/eh3vXzqbf/z+a9zwuS18Y3sn/f0aLxCRqW8y7gQeBlYlNphZPnA/cBPQDKwxs+ZJ+OzRzZgL06uh891bUTfMLOW+269i453vYU5FCX/29Z/xwX/+EdteP5rWEkVE0i3lIeDuzwHDf3q2Ah2x3/zPAo8Bt47nema21szazKzt0KFDEy/MDCKt0Llt1JdEm2by7398Df/woSvo6T3Dbz+wlY9+dQf7j56a+OeKiGSwdI0J1AOJK7U6gXozqzKzB4AVZnbPSG9093XuHnX3aHX1iAfjjF8kCkc64NTov+Hn5Rn/ZUWEZz51HZ94/yKe2dPD+z7/LH/77T0cP30uuc8XEckwoQ4Mu/sRd7/T3Re4+72T/oGRluBxhC6h4UqLCvjE+y/lmU9dx+rltXxxy8+54XNb2PDCm5zXeIGITBHpCoEuoCHheSTWll71V4HljdklNFxt5TQ+/6Er+Y+PXsPcqunc88Qr3PyPz/PjjsOTWKiISHqkKwS2AYvMbJ6ZFQG3AZvS9NmDiqbD7GVDZgiN15UNM9h453u47/YVHD/dx+0P/pQ/fKSNXxw+OQmFioikx2RMEd0AbAUWm1mnmd3h7n3AXcDTwG7gcXfflerPHpdIC3Ruh/6L3z/IzFh9eR3f/7Pr+PP/tJgfdxzmA//wLH/1ZDvHTmm8QESyj2XT/jnRaNTb2i7cnz+mlzbAf9wJf/wTqFma1KV6jp/m8995la+17WfGtEI+eeOl3N7aSEF+7izEFpHMZ2bb3T060vdy76fVwODw/ovvEhqupryEz/7m5Tz5sV9jyZwK/vKbu1j1hefZsrcn6WuLiKRD7oVA1QKYdslFDQ5fyLK6Sh79w6tZ97sr6Tvfz3/78jY+vP4FXjt4PGWfISIyGXIvBMxi4wJJdiu967LGB5bN4TufvI7/ffNSdrz5S1Z94Xn+8ps7OXrybEo/S0QkVXIvBCAIgUN74PSxlF+6qCCPP7h2Pls+dT23tzby1Z++yXV/9wMefH4fZ/v6U/55IiLJyNEQiAIOXdsn7SOqyor5qw9exrc+fi0rGi/hr5/azQf+4Vm+s+stHWYjIhkjN0OgfiVgKe8SGsmls8t55COtfPn3WyjIz2PtV7Zz+7/8lPYDvZP+2SIiF5KbIVBSCdVLUjo4fCE3LK7hWx+/ls/csozdb/Vy8z89z93feJlDx8+krQYRkeFyMwQg6BLq3AZp7JopzM/jw7/axLOfuoGPXDOPjds7ueFzW/jnLR2cPnfxi9dERJKVwyHQAu/8Eo78PO0fXVlayP9Z3cx3PvlefmV+Ff/323t5/+ef5amXuzVeICJplbsh0NAaPKaxS2i4+dVlPPjhKF/9g6spKy7go4/u4He+tJWXO98OrSYRyS0FYRcQmlmLobgi2EzuyjWhlnLNwlk89SfX8njbfv7+O3u55b4f8Z75VSyPVLKsroJldRXMm1VGfp6FWqeITD25GwJ5ecHW0iHeCSTKzzPWtDay+vJa1j23jy17D/Hwj17n7PlgbcG0wnyW1JazrK6C5togHBbPKaekMD/kykUkm+XeBnKJnvkbeP5zcE9nsM10hjl3vp+OnhO0H+hl14Fedh04Rnt3L8dP9wFBcCysLguCoa6CZXWVNNdVUDmtMOTKRSSTjLWBXO7eCUAwOOz90LUD5l0bdjXvUpifx9LaCpbWVvCbK4M2d2f/0XfYdeBYPBh+9PPDPPHi4Bk9kUumxbqRKuOPsyuKMVN3kogMleMhEAvGzm0ZGQIjMTMaq0pprCrlpuW18fZDx8/Q3t0bD4f2A708vetg/PtV04uG3C0sq6tgXtV08jTOIJLTcjsESmdC1cK0rByebNXlxVxXXs11l1bH206c6WNP92BX0q4DvTz0w32cOx90AZYW5bO0tiI++LysrpJFs8soLtA4g0iuCDUEzGw+8L+ASnf/rVCKiLRAx/eCRWNTrLukrLiAaNNMok0z421n+/p5red4/G6h/UAvT+zo4pGtbwBQkGcsrClL6EoKxhvKSzTOIDIVTTgEzGw9sBrocffLEtpXAV8A8oEH3f2zo13D3fcBd5jZxonWkbRIFH62Ad5+Ay5pCq2MdCkqyIv9gK+Mt/X3O28ePTXkjuHZVw/xjR2d8dfMrSqNzUwaHGuoqSgJ468gIimUzJ3Aw8B9wCMDDWaWD9wP3Ah0AtvMbBNBINw77P0fcffwj+CKDCwaa8uJEBhJXp7RNGs6TbOmc/Plg+MMPb2ngzuGhLGGza+8Ff/+rLLiIV1JzXUVzJ1ZqnEGkSwy4RBw9+fMrGlYcyvQEfsNHzN7DLjV3e8luGvIPDXNUFgaDA4vD6dHKlPVVJRQU1HCDUtq4m29p8+xOz5lNQiIHz23j77+YJyhrLiApbXl8VBorq3g0tnlFBXk7uJ0kUyW6jGBemB/wvNO4OrRXmxmVcDfACvM7J5YWAx/zVpgLUBjY2NqqwXIL4C6q1Jy5nAuqCgp5Or5VVw9vyredqbvPK8dPJEwbbWXx9v2c+pssCleYb6xqKY8YYyhkqW15RpnEMkAoQ4Mu/sR4M4LvGYdsA6CxWKTUkhDC/z4n+DcO1A4bVI+YiorLsjnsvpKLqsfOs7w+pGT8VDYdeAYz+zp4evbB8cZmqpKB+8YYgFRU65xBpF0SnUIdAENCc8jsbbMFmmB/j7ofhkaR71xkYuQl2fMry5jfnUZ//mKOiBY6NZz/Exwx9AVhMMrXcd46pXu+PuqywfHGQa2x2jUOIPIpEl1CGwDFpnZPIIf/rcBt6f4M1Iv0hI8dr6gEJhEZsbsihJmV5Tw60tmx9uPvXOO3QnrGdoP9PLD1w4PGWdorh28W2iuq2BRjcYZRFIhmSmiG4DrgVlm1gl82t0fMrO7gKcJZgStd/ddKal0MpXVwIy5GbOZXK6pnFbIr8yv4lcSxhlOnxs6ztDePXScoSg/j0Wzy4bMTFpaW0FZcW6vfxS5WMnMDhpx/2V33wxsnnBFYYm0wJtbw65CYkoK81keqWR5ZHCc4fyQcYbgjuH7u3t4vC0YZzCDpqrpg3cMsTUN1eXFYf01RDKefm0aEGmBnRvhWBdU1oddjYwgP89YUF3GguoybkkYZzjYe2bIhnovd77NUy8PjjPUxMcZBvdNapxZqg31RFAIDIqPC2xTCGQRM2NOZQlzKkt439Kh4wztCXcM7d29PPfaYc7HxhnKiwtYqnEGEYVA3JzlkF8chMCyD4ZdjSSpcloh71lQxXsWDB1nePXg8SHnM3xt237eOTd0nKF5YFO9+kqWzNF6BpnaFAIDCoqg7koNDk9hJYX5XB6ZweWRGfG2xHGG9jHWMwzZhrtW+ybJ1KEQSBRpgRf+BfrOBqEgU95o4wzD1zPs7Bp536TmhL2TtG+SZCOFQKJIC2y9Dw6+AvUrw65GQjLaeobe0+fi22/HT3XrGFzPMD3hfIaBOwedzyCZTiGQaGBweP82hYC8S0XJu9czDOyb1J6wDffG7Z2c3BqMMySez5A4CF2hcQbJEAqBRJX1UF4XGxcYc0sjEWD4vknBjin9/c4bR08NOepz+PkMDTOnsSy2LcbAXYPOgZYwKASGa2jR4LAkJS/PmDdrOvNmTWf15XXx9p7e0+zq7k3oUjrGt3cNjjMMnAMdH4SurWDerOnka5xBJpFCYLhIC7R/E070BNtJiKRI/HyGxYP/XR0/fY49bx1nV9fg9hjrf/iL+DnQ0wrzWVJbPmRDvcVzyikp1DiDpIZCYLjERWNLbg63FpnyyksKaWmaScsI50C3Jxzc880XD/BvP3kTGJjRND1hU73gruGS6ZrRJhdPITBc7RWQV6gQkNAkngP927G2/n6n85fvBCugY11KP9l3lP946UD8fXWVJbHupMr4grfIJdM0ziBjUggMVzgtWD3c2RZ2JSJxeXlGY1UpjVWl3LR88BzoIyfOxEOhPbYd9zN7eojNWqW8pCC+kd7A7KSFNWUU5mt7DAkoBEYSaYEXvwLn+4LjJ0UyVFVZMdcuqubaRdXxtnfOnmfPW4Oh0H6gl0dfeIPT5/qBYHuMS+cE22M012p7jFynn3AjaWiFF74EPe1Qe3nY1YhclGlF+axovIQVjZfE2873O784fCIeCu3dvXwvYRtuGLY9Rmy8oaZc01anOoXASCLR4LFzm0JApoT8PGNhTTkLa8q59cpgl9yBbbjbu4PtMQbuHIZuj1EUWwU92J3UVKVpq1OJQmAkM+bC9OogBFruCLsakUmRuA338O0x9nQfH7IN90M/3Ddk2urS2vJgEFrTVrNeqCFgZkuBjwOzgO+7+xfDrCfODCKtWjQmOamipJDWeTNpnTd02mpHz4n47KRdB8aettpcG9w5zNS01YyXzBnD64HVQI+7X5bQvgr4AsEZww+6+2dHu4a77wbuNLM84BEgM0IAgi6hvU/BqaNQOvPCrxeZwooK8uKrmQe4J0xbja1p+Okvhk5bra0sSQiGoFspcsk07baaQZK5E3gYuI/ghzcAZpYP3A/cCHQC28xsE0Eg3Dvs/R9x9x4zuwX4I+ArSdSSevFFY21w6QfCrUUkA5kZDTNLaZhZyqrLBqetHj15lt1Dpq0eY8urh0Y81W0gILTbaniSOWj+OTNrGtbcCnS4+z4AM3sMuNXd7yW4axjpOpuATWb2FPDo8O+b2VpgLUBjY+NEy714dSvA8oIuIYWAyLjNnF7ENQtncc3CWfG2xFPdBtY1fL1tPyfPDt1tNTEYltVWUlmqaauTLdVjAvXA/oTnncDVo73YzK4HfgMoBjaP9Bp3XwesA4hGo56qQi+ouAxmL9O4gEgKjHSq28Buq0EwBF1KP3ztME/s6Iq/pn7GtKHBUFdB/Qytgk6lUAeG3X0LsCXMGsYUaYFXNkJ/P+RphaVIKiXutnrz5YPdSYeOnwm6k7oHd1v93u6DeOxXwIqSgiGDz821wSroogL9PzoRqQ6BLgY2VQ9EYm3ZKdICbevh8F6oWRp2NSI5obq8mOryat576eAq6FNn+9j71vEhW2RseOFN3jkXdCcV5eexaHbZkEHopTq8Z1xSHQLbgEVmNo/gh/9twO0p/oz0ibQGj53bFAIiISotKhhlFfTJIcHwg709fH374CroxpmlQ4Khua6C2soSdSclSGaK6AbgemCWmXUCn3b3h8zsLuBpghlB6919V0oqDUPVAiiZAftfgKt+L+xqRCRBfmwweWFNGbdcMcrhPd297D7Qy9Ptb8W7ky4pLRwSCs21lcyvnp6zm+olMztozSjtmxllkDfrmAVdQtpRVCRrjHR4z4kzfex9q3fI7KRHtr7Bmb7YpnoFeSyeXT4YDHUVObOpnraNuJCGVuj4Hpw+BiWVYVcjIhNQVlzAyrkzWTl3cOFn3/n+d3UnfXf3Qb7WNjjBcW5VaXy31YFwmFMxtbqTFAIXEokCDl07YMENYVcjIilSkJ/HotnlLJo9dFO9nuNnhtwxtHf38q2dg5vqDXQnLZ0zGAwLqrP3jAaFwIXUrwQsGBxWCIhMaWbG7IoSZleUcMOSsbuTvvKThO6kYWc0NNdVsrQ2O7qTFAIXUlIJ1Uu0aEwkh423O+n7w85oyIbZSQqB8YhEYc+T4B4MFotIzhtXd9IIs5NmlBa+a5whzO4khcB4DBw3eeTnMGth2NWISIYarTvp5Jk+9gxb7Hah7qQlteVpWeymEBiPhoRFYwoBEblI04sLWDn3ElbOHVzs1ne+n9ePnAyO/IyFw1jdSb+6oIpoU+q3tVcIjMesxVBcEYTAlSMujxARuSgF+XkjHvl56PiZERe79RxvVAiEJi8P6q+CzhfCrkREpjAzG3Gx28kzffF9klItOye2hiHSAgd3wdmTYVciIjlmenEBs8qKJ+XaCoHxirSC98OBF8OuREQkZRQC4xWJBo/71SUkIlOHQmC8SmfCzAXaTE5EphSFwMVoaA1mCHn6TrkUEZlMCoGLEYnCyR54+42wKxERSQmFwMWItASP6hISkSlCIXAxapZBYak2kxORKSPUEDCz683seTN7wMyuD7OWcckvgLqrNENIRKaMCYeAma03sx4z2zmsfZWZ7TWzDjO7+wKXceAEUAJ0XuC1mSEShbdehnPvhF2JiEjSkrkTeBhYldhgZvnA/cBNQDOwxsyazWy5mT057KsGeN7dbwL+B/CZJGpJn4ZW6O+D7pfDrkREJGnJHDT/nJk1DWtuBTrcfR+AmT0G3Oru9wKrx7jcL4ER10Sb2VpgLUBjY+NEy02d+tiisc5t0Hh1uLWIiCQp1WMC9cD+hOedsbYRmdlvmNmXgK8A9430Gndf5+5Rd49WV1entNgJKZ8NMxq1mZyITAmh7iLq7k8AT4RZw4REWuHNrWFXISKStFTfCXQBDQnPI7G2qSXSAr1dcGzq/dVEJLekOgS2AYvMbJ6ZFQG3AZtS/Bnhiy8a03oBEcluyUwR3QBsBRabWaeZ3eHufcBdwNPAbuBxd9+VmlIzyJzlkF+sEBCRrJfM7KARz1l0983A5glXlA0KiqDuSm0fISJZT9tGTFSkJThgpu9s2JWIiEyYQmCiIlE4fwYOvhJ2JSIiE6YQmKhIa/CoLiERyWIKgYmqrIfyOm0mJyJZTSGQjEhUM4REJKspBJLR0BqcMnaiJ+xKREQmRCGQDC0aE5EspxBIRu0VkFegEBCRrKUQSEbhNJhzuWYIiUjWUggkK9ICXTvgfF/YlYiIXDSFQLIiLXDuJPS0h12JiMhFUwgkq0GDwyKSvRQCyZoxF6ZXa1xARLKSQiBZZkGXkI6bFJEspBBIhUgLHOmAU0fDrkRE5KIoBFJhYNFY1/Zw6xARuUihhoCZXWtmD5jZg2b24zBrSUrdCrA8bSYnIlknmeMl15tZj5ntHNa+ysz2mlmHmd091jXc/Xl3vxN4EvjXidYSuuIymL1MM4REJOskcyfwMLAqscHM8oH7gZuAZmCNmTWb2XIze3LYV03CW28HHk2ilvBFWoLuoP7+sCsRERm3CYeAuz8HDB8JbQU63H2fu58FHgNudfdX3H31sK8eADNrBI65+/GJ1pIRIi1wphcO7w27EhGRcUv1mEA9sD/heWesbSx3AF8e7ZtmttbM2sys7dChQykocZJoR1ERyUKhzw5y90+7+6iDwu6+zt2j7h6trq5OZ2kXp2ohlMxQCIhIVkl1CHQBDQnPI7G2qW9g0dh+hYCIZI9Uh8A2YJGZzTOzIuA2YFOKPyNzRVrg0B44fSzsSkRExiWZKaIbgK3AYjPrNLM73L0PuAt4GtgNPO7uu1JTahZoaAE82FpaRCQLFEz0je6+ZpT2zcDmCVeUzepXAhZsJrfghrCrERG5oNAHhqeUkkqoXqzN5EQkaygEUi3SEswQcg+7EhGRC1IIpFqkBd75JRzdF3YlIiIXpBBItYFFY9pMTkSygEIg1aqXQHGFFo2JSFZQCKRaXh7UX6UQEJGsoBCYDJEWOLgLzp4MuxIRkTEpBCZDpBX8PBx4MexKRETGpBCYDJFo8KguIRHJcAqByVA6E2Yu0GZyIpLxFAKTpaFVi8ZEJOMpBCZLJAone+DtN8OuRERkVAqByaKTxkQkCygEJkvNMigsVQiISEZTCEyW/AKo06IxEclsCoHJFIlC98tw7nTYlYiIjEghMJkiLdB/Drp/FnYlIiIjCjUEzKzZzB43sy+a2W+FWcuk0OCwiGS4ZM4YXm9mPWa2c1j7KjPba2YdZnb3BS5zE/BP7v5HwO9NtJaMVT4bZjQqBEQkY034jGHgYeA+4JGBBjPLB+4HbgQ6gW1mtgnIB+4d9v6PAF8BPm1mtwBVSdSSuSIt8OZPwq5CRGREyRw0/5yZNQ1rbgU63H0fgJk9Btzq7vcCq0e51Edj4fHESN80s7XAWoDGxsaJlhueSCvs/AYc64LK+rCrEREZItVjAvXA/oTnnbG2EZlZk5mtI7ib+LuRXuPu69w96u7R6urqlBabFgPjAl1t4dYhIjKCZLqDkuburxP7LX/KmrMc8ouD4yabbw27GhGRIVJ9J9AFNCQ8j8TacldBEdRdCZ26ExCRzJPqENgGLDKzeWZWBNwGbErxZ2SfSAt0vwR9Z8OuRERkiGSmiG4AtgKLzazTzO5w9z7gLuBpYDfwuLvvSk2pWSwShb7TcPCVsCsRERkimdlBa0Zp3wxsnnBFU1GkNXjsbIP6leHWIiKSQNtGpENlPZTXadGYiGQchUC6RKLBDCERkQyiEEiXhlZ4+w040RN2JSIicQqBdIlvJqepoiKSORQC6VJ7BeQVQKe6hEQkcygE0qVwWrB6WHcCIpJBFALpFGmFrh1wvi/sSkREAIVAekVa4NxJOLQ77EpERACFQHpFosGjpoqKSIZQCKTTJU0wvVrjAiKSMRQC6WQWdAlp5bCIZAiFQLpFonDkNTh1NOxKREQUAmk3sJlc1/Zw6xARQSGQfnUrwPLUJSQiGUEhkG7FZVCzTDOERCQjKATC0NASdAf194ddiYjkOIVAGCItcKYXDr8adiUikuPSFgJmNt/MHjKzjWO15YT4jqLqEhKRcI0rBMxsvZn1mNnOYe2rzGyvmXWY2d1jXcPd97n7HRdqywlVC6FkhgaHRSR04z1j+GHgPuCRgQYzywfuB24EOoFtZrYJyAfuHfb+j7i7TlMZEF80ppXDIhKucYWAuz9nZk3DmluBDnffB2BmjwG3uvu9wOpUFWhma4G1AI2Njam6bPgiLbDle3D6GJRUhl2NiOSoZMYE6oH9Cc87Y20jMrMqM3sAWGFm94zWNpy7r3P3qLtHq6urkyg3wzS0AB5sLS0iEpLxdgclzd2PAHdeqC1n1K8ELOgSWnBD2NWISI5K5k6gC2hIeB6Jtcl4lFRC9WLNEBKRUCUTAtuARWY2z8yKgNuATakpK0dEosEMIfewKxGRHDXeKaIbgK3AYjPrNLM73L0PuAt4GtgNPO7uuyav1Cko0grv/BKO7gu7EhHJUeOdHbRmlPbNwOaUVpRL4ovGtkHVgnBrEZGcpG0jwlS9GIrKtZmciIRGIRCmvHyIrNTKYREJjUIgbJEWOLgLzp4MuxIRyUEKgbBFWsDPw4EXw65ERHKQQiBsiYPDIiJpphAIW+lMmLlAm8mJSCgUApkg0hLMENKiMRFJM4VAJmhogZM98PabYVciIjlGIZAJNC4gIiFRCGSCmmVQWKoQEJG0UwhkgvwCqLtKISAiaacQyBSRKHS/DOdOh12JiOQQhUCmiLRA/zno/lnYlYhIDlEIZAoNDotICBQCmaJ8NsxoVAiISFopBDJJpEUrh0UkrdIWAmY238weMrONCW1LzewBM9toZn+UrloyVqQFejuh90DYlYhIjhjv8ZLrzazHzHYOa19lZnvNrMPM7h7rGu6+z93vGNa2293vBH4HuOZii59yIq3Bo7qERCRNxnsn8DCwKrHBzPKB+4GbgGZgjZk1m9lyM3ty2FfNaBc2s1uAp9AxlTBnOeQXKwREJG3Ge8bwc2bWNKy5Fehw930AZvYYcKu73wusHm8B7r4J2GRmTwGPjvd9U1JBEdReAfsVAiKSHsmMCdQD+xOed8baRmRmVWb2ALDCzO6JtV1vZv9oZl9ilDsBM1trZm1m1nbo0KEkys0SDa3Q/RL0nQ27EhHJAeO6E0gFdz8C3DmsbQuw5QLvWwesA4hGo1N/r+VIFLbeBwd3Qv1VYVcjIlNcMncCXUBDwvNIrE2SoUVjIpJGyYTANmCRmc0zsyLgNmBTasrKYZURKK9TCIhIWox3iugGYCuw2Mw6zewOd+8D7gKeBnYDj7v7rskrNYdcuQaql4RdhYjkAPMsOtIwGo16W5tW1IqIXAwz2+7u0ZG+p20jRERymEJARCSHKQRERHKYQkBEJIcpBEREcphCQEQkhykERERymEJARCSHZdViMTM7BLwxwbfPAg6nsJzJlk31ZlOtkF31ZlOtkF31ZlOtkFy9c929eqRvZFUIJMPM2kZbMZeJsqnebKoVsqvebKoVsqvebKoVJq9edQeJiOQwhYCISA7LpRBYF3YBFymb6s2mWiG76s2mWiG76s2mWmGS6s2ZMQEREXm3XLoTEBGRYRQCIiI5LCdCwMxWmdleM+sws7vDrmcsZrbezHrMbGfYtVyImTWY2Q/MrN3MdpnZx8OuaSxmVmJmL5jZz2L1fibsmi7EzPLN7EUzezLsWi7EzF43s1fM7CUzy+jTn8xshpltNLM9ZrbbzN4Tdk0jMbPFsX/Pga9eM195jk0AAALYSURBVPtESj9jqo8JmFk+8CpwI9BJcDbyGndvD7WwUZjZe4ETwCPuflnY9YzFzGqBWnffYWblwHbggxn8b2vAdHc/YWaFwA+Bj7v7T0IubVRm9qdAFKhw99Vh1zMWM3sdiLp7xi/AMrN/BZ539wdjZ6SXuvvbYdc1ltjPsi7ganef6KLZd8mFO4FWoMPd97n7WeAx4NaQaxqVuz8HHA27jvFw92533xH783GCs6brw61qdB44EXtaGPvK2N+CzCwC3Aw8GHYtU4mZVQLvBR4CcPezmR4AMe8Dfp7KAIDcCIF6YH/C804y+AdVtjKzJmAF8NNwKxlbrHvlJaAH+K67Z3K9/w/4C6A/7ELGyYHvmNl2M1sbdjFjmAccAr4c62p70Mymh13UONwGbEj1RXMhBGSSmVkZ8A3gE+7eG3Y9Y3H38+5+JRABWs0sI7vczGw10OPu28Ou5SL8mrtfBdwEfDTWtZmJCoCrgC+6+wrgJJDpY4VFwC3A11N97VwIgS6gIeF5JNYmKRDrW/8G8FV3fyLsesYrdvv/A2BV2LWM4hrgllg/+2PAr5vZv4Vb0tjcvSv22AP8O0FXbCbqBDoT7gI3EoRCJrsJ2OHuB1N94VwIgW3AIjObF0vT24BNIdc0JcQGWh8Cdrv758Ou50LMrNrMZsT+PI1gssCecKsambvf4+4Rd28i+G/2GXf/ryGXNSozmx6bHECsa+UDQEbOcHP3t4D9ZrY41vQ+ICMnMyRYwyR0BUFwWzSluXufmd0FPA3kA+vdfVfIZY3KzDYA1wOzzKwT+LS7PxRuVaO6Bvhd4JVYPzvA/3T3zSHWNJZa4F9jsyzygMfdPeOnXmaJ2cC/B78XUAA86u7fDrekMX0M+GrsF8N9wO+HXM+oYqF6I/DfJ+X6U32KqIiIjC4XuoNERGQUCgERkRymEBARyWEKARGRHKYQEBHJYQoBEZEcphAQEclh/x+8Q7GJNrF5DwAAAABJRU5ErkJggg==\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "print(\"---------------------------------------------------\")"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "tsvZHsbYbE_h",
        "outputId": "fe5baf90-b7bb-40c8-9165-891a27e10b2b"
      },
      "execution_count": 35,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "---------------------------------------------------\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import matplotlib.pylab as plt\n",
        "import numpy as np\n",
        "from scipy.sparse import coo_matrix\n",
        "from scipy.sparse.linalg import gmres, LinearOperator\n",
        "import pyamg\n",
        "\n",
        "def discretise_poisson(N):\n",
        "    \"\"\"Generate the matrix and rhs associated with the discrete Poisson operator.\"\"\"\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\n"
      ],
      "metadata": {
        "id": "TpqA5Y5edrmx"
      },
      "execution_count": 36,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "A, b = discretise_poisson(200)"
      ],
      "metadata": {
        "id": "m_t5sxE8d0Wt"
      },
      "execution_count": 37,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "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",
        "plt.plot(residuals)\n",
        "plt.yscale(\"log\")\n",
        "plt.show()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 265
        },
        "id": "0TdPKeLAd3b4",
        "outputId": "e78eec73-c8e8-4ab8-d367-3cf6123a4c63"
      },
      "execution_count": 40,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3RVVdrH8e+TRqQFkaYUqVIjKk1ACSoloIiiDiDjOKOCgA0z7zjgODM6jqK+IwqKInZUUF5ABUS6JnQpjiShBpQmkCBKb4H9/pGLxkgJuTX3/j5r3cU9O/ee85y9Vh5O9nnO3uacQ0REwl9UsAMQEZHAUMIXEYkQSvgiIhFCCV9EJEIo4YuIRIiYYAdwJhUqVHA1a9YMdhgiIsXK8uXLdznnKhZsD+mEX7NmTZYtWxbsMEREihUz23Sqdg3piIhECCV8EZEIoYQvIhIhlPBFRCKEEr6ISIQIWJWOmZUCXgGOAl865z4I1LFFRMTLK3wze8vMss0so0B7spmtNbMsMxvsae4BTHDO9QVu9Oa4IiJy7rwd0nkHSM7fYGbRwEigC9AI6G1mjYBqwBbPx457edwzSl2Xw+tpG/15CBGRYserhO+cSwN2F2huCWQ55zY6544CHwLdga3kJf0zHtfM+pnZMjNblpOTU6S4Zq/ayXMz1vD9T4eK9H0RkXDkj5u2VfnlSh7yEn1VYBJwi5m9Ckw53Zedc6Odc82dc80rVvzNk8GF0q9dbQBemru+SN8XEQlHAavScc4dcM79yTk34Gw3bM2sm5mN3rNnT5GOVb18Sfq0upiPlm7hmy0/FWkfIiLhxh8JfxtQPd92NU9boTnnpjjn+iUkJBQ5iIeuq0eVsvHc+fZXZH5ftP84RETCiT8S/lKgnpnVMrM4oBcw+Vx24O0VPsD5peIY1+9KzouNps8bS1j1/d4i70tEJBx4W5Y5DlgE1DezrWZ2t3MuF7gfmAGsBsY75zLPZb++uMIHuPiCUnzoSfo9X1vE3DU7vdqfiEhxZs65YMfwG2bWDehWt27dvuvXe3/jddtPh+g3ZhmZ3+/l3qTa/LljfeJi9JCxiIQnM1vunGtesD0ks56vrvBPqlruPCYOaEPvljV4LXUjPV5dQFb2fp/sW0SkuAjJhO8P8bHRDO2RyKjfN2Pbj4e44aV5vLd4E6H4F46IiD+EZML3xU3b00luUoXpg9rRomZ5/v5JBne8+RVbfzzo8+OIiISakBzDP6l58+bOX0scOucY+9Vmnv5sNWbG365vSK8W1TEzvxxPRCRQitUYfiCYGX1aXcz0Qe24tFoCQyal0+eNJWz64UCwQxMR8YuQTPj+HNIpqHr5krx/dyuG9kgkfeseOr+YxutpGzl+InT/8hERKYqIHdI5lR17DvPYJ+nMXp1N02oJPHvrpTSoUjZgxxcR8QUN6RRClYR4Xv9Dc17qfTlbfzzEDSPmM2zWOo7k+nU2ZxGRgFDCL8DM6Nb0ImalJNGt6UWMmLOeG0bMZ8XmH4MdmoiIV0Iy4QdyDP90ypeK44Wel/H2H1tw4Egut7y6kH9NWcXBo7lBi0lExBsawy+EfYeP8dz0tby3eBPVzj+PZ3pcylX1KgQ7LBGRU9IYvhfKxMfy5E1NGH9va+Kio/j9m0t4ZMI37Dl4LNihiYgUmhL+OWhZqzzTHrqage3rMHHFNjq8kMr0jO3BDktEpFCU8M9RfGw0jyQ34NP72lKxdAn6v7+CAe8vJ3vf4WCHJiJyRiGZ8EPhpu3ZNKmawKf3t+WR5PrMWZNNh+dT+WjpZk3GJiIhSzdtfWBjzn4GT0rnq29306bOBTx9cyI1K5QKdlgiEqF009aPalcszYd9r+Tpm3+ZnuG11A3kHj8R7NBERH6mhO8jUVHG7a1qMCsliXaXVGTo52u46ZUFZGwL3WEpEYksSvg+ViUhntF3NOOVPlewY88Ruo9cwDOfr+HwMU3PICLBpYTvB2ZG18QLmZ3SjluuqMqo1A0kv5jGog0/BDs0EYlgSvh+VK5kHM/d2pQP7mnFCQe9X1/M4Ikr2XNID2yJSOCFZMIvDmWZ56Jt3QrMGNSOe9vVZvyyLXQclsr0jB3BDktEIozKMgMsfeseHpm4ktXb95LcuAr/6t6YSmXjgx2WiIQRlWWGiMRqCUz2PLA1d2021w1L5cOv9MCWiPifEn4QxEZHMbB9XWYMakejC8syeFI6vV9fzHe7tJ6uiPiPEn4Q1apQinF9r2Roj0Qyv99L5xfTePVLPbAlIv6hhB9kUVFG75Y1mJ2SRPv6FXl2+hq6j9QDWyLie0r4IaJy2Xheu6M5o35/Bdn78h7YGvr5ag4d1QNbIuIbSvghJrnJhcx+OInbmlXjtdSNJA9PY2HWrmCHJSJhIGAJ38xqm9mbZjYhUMcsrhJKxvLMLZcy9p5WANz+xhL+OmGlVtgSEa8UKuGb2Vtmlm1mGQXak81srZllmdngM+3DObfROXe3N8FGmjYnH9hKqs2EFVvp8EIqn6drhS0RKZrCXuG/AyTnbzCzaGAk0AVoBPQ2s0ZmlmhmUwu8Kvk06ggSHxvNkC4N+fS+tlQqU4IBH6zg3veWsXOvVtgSkXNTqITvnEsDdhdobglkea7cjwIfAt2dc+nOuRsKvLJ9HHfEaVI1gU/va8vgLg34cm0OHZ5PZeySzZw4oQe2RKRwvBnDrwpsybe91dN2SmZ2gZmNAi43syFn+Fw/M1tmZstycnK8CC/8xERH0T+pDtMHtaNx1bI8+nHeA1sbc/YHOzQRKQYCdtPWOfeDc66/c66Oc27oGT43GngCWBEXFxeo8IqVkw9sPdMjkVXb95I8fB6vfJnFMT2wJSJn4E3C3wZUz7ddzdPmNefcFOdcv4SEBF/sLiyZGb1a1mBOShLX1q/Ec9PX0v3lBaRv1QNbInJq3iT8pUA9M6tlZnFAL2CyL4IKt+mR/alS2XhG3dGMUb+/gl37j9B95HyenqYHtkTktwpbljkOWATUN7OtZna3cy4XuB+YAawGxjvnMn0RlK7wz11ykwuZlZJEzxbVGZ22kc4vprFAD2yJSD6aDz8MLdywi0cnpfPdDwe5rVk1Hru+EQklY4MdlogESLGaD19DOt5pU6cC0we1o39SHSZ9vY3rhqUyLX275twXiXC6wg9zGdv2MHjSSjK27aVjo8o82b0JVRK0wpZIONMVfoRqUjWBTwa2ZUiXBqSty6HjsFQ+WLJJD2yJRCBd4UeQ73Yd4NGP01m44Qda1irP0B6J1KlYOthhiYiPFasrfPGPmhVK8cE9rXjulktZs30vXYbPY+QXemBLJFIo4UcYM+N3LaozOyWJDg0r8b8z1tLtpfl8s+WnYIcmIn4WkglfY/j+V6lsPK/0acZrdzTjx4NHufmVBTw5dRUHjuQGOzQR8RON4Qt7Dh3j2elrGLtkM1XLncdTNzehfX3NaC1SXGkMX04r4bxYnr45kf/r35r42Cj++PZSBn34NbsPHA12aCLiQyGZ8DWkExwtapZn2kNX8+B19fgsfTsdhqXyydfb9MCWSJjQkI6c0tod+/jrxJX8d8tPtLukIk/d1ITq5UsGOywRKQQN6cg5qV+lDBMHtOHxbo1Y9t1uOr2Qxpvzv+W4HtgSKbaU8OW0oqOMP7atxayUJFrVLs+TU1fR49WFrN6+N9ihiUgRKOHLWVUtdx5v/7EFw3tdxpbdB+n20nz+M2Mth49pzn2R4iQkE75u2oYeM6P7ZVWZk5JE98uq8vIXWXQdPo8lG38IdmgiUki6aStFMm99DkMmpbP1x0P0aVWDv3ZpQNl4zbkvEgp001Z86up6FZn5cDvuuaoW477aTKdhaczM3BHssETkDJTwpchKxsXw2A2N+HhgW8qVjKXfe8sZ+MFysvcdDnZoInIKSvjitabVyzHlgav4S+f6zF6dTYfnU/lo6WY9sCUSYpTwxSdio6O475q6TH/oahpeWJa/Tkyn9+uL+XbXgWCHJiIeSvjiU7UrlmZc3ysZ2iORzO/30vnFNF75UnPui4SCkEz4Ksss3qKijN4tazAnJYnrGlTiuelrufHlBZpzXyTIVJYpfjcjcwf/+DSDnH1HuKttLVI6XULJuJhghyUStlSWKUHTuXEVZqUk0btlDd6Y/y2dXkgjdV1OsMMSiThK+BIQZeNjeermRMbf25q4mCjufOsrUj76r+bcFwkgJXwJqJa1yjPtwat58Nq6TFn5vebcFwkgJXwJuPjYaFI61WfqA1dTo3xJBn30X+58eylbdh8MdmgiYU0JX4Im/5z7yz1z7r+etpFclXCK+IUSvgRV/jn329S5gKemrebmVxaS+b1KckV8TQlfQsJF5c7jjTub8/Ltl7N9zyFufHkBz05fozn3RXwooAnfzG4ys9fN7CMz6xTIY0voMzNuuPQiZqckccsVVXn1yw0kv5jGwqxdwQ5NJCwUOuGb2Vtmlm1mGQXak81srZllmdngM+3DOfeJc64v0B/oWbSQJdyVKxnHc7c2Zew9rXDA7W8s4ZEJ3/DTQZVwinjjXK7w3wGS8zeYWTQwEugCNAJ6m1kjM0s0s6kFXpXyffUxz/dETqtN3QrMGNSO/kl1mLhiGx2GpTLlm+9VwilSROc0tYKZ1QSmOueaeLZbA4875zp7tocAOOeGnub7BjwDzHLOzT7NZ/oB/QBq1KjRbNOmTYWOT8JX5vd7GDIpnZVb93Btg0o8eVMTqpY7L9hhiYQkf02tUBXYkm97q6ftdB4AOgC3mln/U33AOTfaOdfcOde8YsWKXoYn4aLxRQlMGtCGx65vyKINP9BpWCrvLPiW4yd0tS9SWAG9aeucG+Gca+ac6++cG3W6z2m2TDmVmOgo7rm6NjMfbkezmuV5fMoqbh21kLU79gU7NJFiwduEvw2onm+7mqfNK865Kc65fgkJCd7uSsJQ9fIlefdPLXihZ1O+23WA60fM4/mZa1XCKXIW3ib8pUA9M6tlZnFAL2Cyt0HpCl/Oxsy4+fJqzPlze25sehEvzc2i64h5LNn4Q7BDEwlZ51KWOQ5YBNQ3s61mdrdzLhe4H5gBrAbGO+cyvQ1KV/hSWOVLxTGs52WMuaslR3NP0HP0Yh79OJ29h48FOzSRkBOSC6CYWTegW926dfuuX78+2OFIMXHwaC7Pz1zH2wu+pWKZEjxxYxOSm1QJdlgiAXe6Kp2QTPgnacUrKYpvtvzE4EnprN6+l+TGVXiie2Mql40PdlgiAaMVryRiNK1ejsn3t+WR5Pp8sTabDs+n8v7iTZxQCadEuJBM+LppK96KjY5iYPu6zBjUjsRqCTz2SQa/e20R63eqhFMil4Z0JOw555iwfCtPTVvNgSO5DGxfl4HX1KFETHSwQxPxCw3pSMQyM25rXp3ZKUl0TbyQ4XPWc/2I+SzftDvYoYkEVEgmfA3piD9UKF2C4b0u5+0/teDQ0ePcOmoRf/8kg30q4ZQIoSEdiUgHjuTyn5lreWfhd1QuE8+TNzWhY6PKwQ5LxCc0pCOST6kSMfyzW2MmDWhDuZKx9B2zjIEfLCd77+FghybiN0r4EtEur3E+Ux64ir90rs/s1dlcNyyVsUs2q4RTwlJIJnyN4UsgxUZHcd81eSWcjS8qy6Mfp9Pr9cVsyNkf7NBEfEpj+CL5OOcYv2wLT322msPHTnD/tXXpn1SHuJiQvDYSOSWN4YsUgpnRs0UNZv85iY6NKzNs1jpueGkeyzf9GOzQRLymhC9yCpXKxDPy9it4887m7Ducy62jFvLPTzPYfyQ32KGJFJkSvsgZXNewMrNSkrizdU3GLN5Ex2GpzF61M9hhiRRJSCZ83bSVUFK6RAyP39iYiQPaUDY+lnvGLOO+D1aQvU8lnFK86KatyDk4mnuC0WkbGDE3i/iYKB7t2pCeLapjZsEOTeRnumkr4gNxMVHcf209pj90NQ0vLMvgSen0Gr2YjSrhlGJACV+kCGpXLM24vlfyTI/EvIVWhs/j5bnrOZp7ItihiZyWEr5IEUVFGb1aeko4G1bmPzPX0e2l+Xy9WSWcEpqU8EW8VKlMPCP7XMEbf2jO3sPH6PHqQh6fnKkSTgk5SvgiPtKhUWVmPtyOP1x5Me8u+o5Ow1KZs1olnBI6QjLhqyxTiqsy8bE80b0JE/q3oXR8DHe/u4z7xq4gZ9+RYIcmorJMEX85mnuCUakbeHluFufFRfO3rg25rXk1lXCK36ksUyTA4mKiePC6ekx76GrqVy7DIxNXcvvrS/h214FghyYRSglfxM/qVirNh/2u5OmbE8n4fg/JL6bxypdZHDuuEk4JLCV8kQCIijJub1WD2SlJXFO/Es9NX8uNLy/gmy0/BTs0iSBK+CIBVLlsPKPuaMZrdzRj94Ej3PzKAv41ZRUHVMIpAaCELxIEnRtXYVZKEre3qsFbC76l0wtpfLE2O9hhSZhTwhcJkrLxsfz7pkQm9G/NeXHR/OntpTw47mt27VcJp/hHwBK+mTU0s1FmNsHMBgTquCKhrnnN8nz24FUM6lCP6Rk76DAslf9btoVQLpmW4qlQCd/M3jKzbDPLKNCebGZrzSzLzAafaR/OudXOuf7A74C2RQ9ZJPyUiIlmUIdLmPbQVdStWJq/TFjJ799cwncq4RQfKuwV/jtAcv4GM4sGRgJdgEZAbzNrZGaJZja1wKuS5zs3Ap8B03x2BiJhpG6lMoy/tzX/vqkJK7fsobNKOMWHCv2krZnVBKY655p4tlsDjzvnOnu2hwA454YWYl+fOeeuP9vn9KStRLIdew7z+ORMpmfuoOGFZXmmRyJNq5cLdlhSDPjjSduqwJZ821s9bacLoL2ZjTCz1zjDFb6Z9TOzZWa2LCcnx4vwRIq3Kgkq4RTfignUgZxzXwJfFuJzo4HRkHeF79+oREJf58ZVaF3nAp6bvoa3FnzLjMwd/PumJlzToFKwQ5Nixpsr/G1A9Xzb1TxtXtNsmSK/9psSzneW8oBKOOUceZPwlwL1zKyWmcUBvYDJvgjKOTfFOdcvISHBF7sTCRsnSzgf7nAJMzJ2cN3zqYxXCacUUmHLMscBi4D6ZrbVzO52zuUC9wMzgNXAeOdcpi+C0hW+yOmViInmoQ71mPbQVVxSuTSPTFhJnzdUwilnp/nwRYqxEycc45Zu5plpazh6/AQPdahH36trExuth+gjWbGaD19X+CKFExVl9Gl1MbP//MssnN1ems9/NQunnIKu8EXCyIzMHfzj0wxy9h3hzjY1+Z9O9SlVImDFeBIiitUVvogUzclZOPu0uph3Fn6XNwvnGs3CKXlCMuFrSEek6MrGx/LkTU2Y0L81JfOVcGohddGQjkgYO5J7nFFfbmTkF56F1K9vyG3NtJB6uNOQjkgE+qWE07OQ+gQtpB7JQjLha0hHxLdOtZD6yC80C2ek0ZCOSITJ3nuYf07O5POMHTSoUoZnbrmUyzQLZ1jRkI6IAFCpbDyv/r4Zo+9oxk8Hj3HzKwt4Ykom+zULZ9hTwheJUJ0aV2FWSjvuuNJTwjkslblrdgY7LPGjkEz4GsMXCYwy8bH8q3teCWepEjHc9c4y7h+7QiWcYUpj+CICwNHcE4xK3cDLcz0lnF0bcltzlXAWRxrDF5EziouJ4sHr8pVwTlQJZ7hRwheRXzlZwjm0R14JZ2eVcIaNkEz4GsMXCa6oKKN3yxrMSUmiQ8NK/O+MvFk4v978Y7BDEy9oDF9EzmrWqp38/ZMMdu47zJ2ta/I/netTWrNwhiyN4YtIkXVsVPnnEs53F+WVcM5ZrRLO4kYJX0QK5ZcSzjaUjo/h7ndVwlncKOGLyDlpdvH5TH3gav7c8RJmZu7kuue/5KOlm7WQejGghC8i5ywuJooHrqvH54OupsGFZfnrxHR6v75YJZwhTglfRIqsTsXSfNg3r4Qz8/u9P5dwHs1VCWcoCsmEr7JMkeJDJZzFh8oyRcSnZq3ayT8+zWDHXpVwBovKMkUkIDo2qszMh9vxB5VwhhwlfBHxuTLxsTxRoITzvrEryN53ONihRTQlfBHxm/wlnLMyd9Lh+VSVcAaREr6I+NWpSjh7jV7Mxpz9wQ4t4ijhi0hAnCzhfKZHIqu37yV5+DxenrteJZwBpIQvIgETFWX0almD2X9OomPDyvxn5jqVcAaQEr6IBFylMvGM7HMFb/yhOXsPH6PHqwt5fLIWUve3gCZ8MytlZsvM7IZAHldEQlOHRpWZlZLEna1r8u6i7+g4LJXZq1TC6S+FSvhm9paZZZtZRoH2ZDNba2ZZZja4ELv6KzC+KIGKSHgqXSKGx29szMQBbSgbH8s9Y5Zx3wcq4fSHQj1pa2btgP3AGOdcE09bNLAO6AhsBZYCvYFoYGiBXdwFNAUuAOKBXc65qWc7rp60FYksR3NPMDptAyPmZhEfE8WjXRvSs0V1LaR+jk73pG2hp1Yws5rA1HwJvzXwuHOus2d7CIBzrmCyP/n9p4BSQCPgEHCzc+43t+fNrB/QD6BGjRrNNm3aVKj4RCR8bMzZz5BJ6Sz5djetapVnaI9EalcsHeywig1/TK1QFdiSb3urp+2UnHN/c84NAsYCr58q2Xs+N9o519w517xixYpehCcixVXtiqUZpxJOnwt4lY5z7p2zDedotkwR+VUJZ6NfSjhXqISzyLxJ+NuA6vm2q3navOacm+Kc65eQkOCL3YlIMVapTDwjb/+lhPOWVxfyz08zVMJZBN4k/KVAPTOrZWZxQC9gsi+C0hW+iBSUv4RzzOJNKuEsgsKWZY4DFgH1zWyrmd3tnMsF7gdmAKuB8c65TF8EpSt8ETkVlXB6JyQXQDGzbkC3unXr9l2/fn2wwxGREHTs+AlGp21k+Jz1xMdEMaRrQ3o2r05UlEo4vS7LDAbV4YvI2WzM2c+jH6ezeONuWnpKOOtEeAmnVrwSkbB0soTz2VsSWbN9L11UwnlaIZnwddNWRM6FmdGzxW9LODUL569pSEdEwk6kL6SuIR0RiRgnF1LPPwvnLJVwhmbC15COiHirTHzsr0o4+45ZxoD3l5O9N3JLODWkIyJhL38JZ4mYKIZ0aUivFuFbwqkhHRGJWLHRUdx3TV1mDGpHk4sSePTjvIXUs7IjayF1JXwRiRi1KpRibN9WPHfLpazduY+uw+fx0pzIKeEMyYSvMXwR8Rcz43ctqjM7JYlOjSvz/Ky8Es7lm8K/hFNj+CIS0Wav2snfPSWcd1x5MX/pXJ8y8bHBDssrGsMXETmF/LNwvrd4Ex2HpTEjc0eww/ILJXwRiXgnZ+H8eGBbypWM5d73lnPve8vYGWYlnEr4IiIel1Uvx5QHruKR5Pp8uTaHDs+n8v7iTZw4EbpD3+ciJBO+btqKSLDERkcxsH1eCWditQQe+ySDnqMXkZW9L9iheU03bUVETsM5x8QV2/j3Z6s4eOQ4911TlwHt6xAXE5LXyj/TTVsRkXNkZtzarBqzU5Lo3KQKL8xex/Uj5rF80+5gh1YkSvgiImdRoXQJXup9OW/9sTkHjx7nllcX8dgn6ew9fCzYoZ0TJXwRkUK6tkHeLJx3ta3F2CWb6fB8KtMzik8JpxK+iMg5KFUihn90a8THA9tyQekS9H9/Of3GLGPHntAv4VTCFxEpgqbVyzH5/rYM7tKA1HU5dByWynshXsIZkglfZZkiUhzERkfRP6kOMx9uR9Pq5fj7Jxnc9toi1u8MzRJOlWWKiPiAc46Pv97Gk1NXceDIcQZeU4cB7etQIiY64LGoLFNExI/MjB5XVGNWShLJTarw4uz1XD9iPku/C50STiV8EREfqlC6BCN6X87bf2rBoaPHuW3UIh79OJ09h4JfwqmELyLiB9fUr8SslHb0vboWH361mQ7DUpmWvp1gDqMr4YuI+EnJuBj+dn0jPr3vKiqVKcHAD1bQd8xytu85FJR4lPBFRPwssVoCn97Xlke7NmB+Vg4dh6Xx7sLvOB7gEk4lfBGRAIiJjqJfuzrMHJTE5TXK8c/Jmdw2aiHrAljCqYQvIhJANS4oyZi7WvJCz6Z8u+sA14+Yx7BZ6ziSe9zvxw5Ywjez9mY2z8xGmVn7QB1XRCTUmBk3X543C2fXxAsZMWc9XYbPY8nGH/x63EIlfDN7y8yyzSyjQHuyma01sywzG3yW3ThgPxAPbC1auCIi4eOC0iUY3uty3r2rJUdzT9Bz9GIGT1zJnoP+KeEs1JO2ZtaOvGQ9xjnXxNMWDawDOpKXwJcCvYFoYGiBXdwF7HLOnTCzysAw51yfsx1XT9qKSKQ4eDSX4bPX88b8bzm/ZByv3XEFzS4uX6R9ne5J25jCfNk5l2ZmNQs0twSynHMbPQf4EOjunBsK3HCG3f0IlDhDoP2AfgA1atQoTHgiIsVeybgYhnRtSLemF/Hs9DVcfEEpnx+jUAn/NKoCW/JtbwVane7DZtYD6AyUA14+3eecc6OB0ZB3he9FfCIixU6Tqgm8d/dpU6lXvEn458Q5NwmYVJjPmlk3oFvdunX9G5SISATxpkpnG1A933Y1T5vXnHNTnHP9EhISfLE7ERHBu4S/FKhnZrXMLA7oBUz2RVCaD19ExPcKW5Y5DlgE1DezrWZ2t3MuF7gfmAGsBsY75zJ9EZSu8EVEfK+wVTq9T9M+DZjm04jQGL6IiD+E5NQKusIXEfG9kEz4IiLieyGZ8HXTVkTE90J6EXMzywE2FfHrFYBdPgynOFNf5FE//EJ9kSdc++Fi51zFgo0hnfC9YWbLTjWXRCRSX+RRP/xCfZEn0vohJId0RETE95TwRUQiRDgn/NHBDiCEqC/yqB9+ob7IE1H9ELZj+CIi8mvhfIUvIiL5KOGLiESIsEz457jWbrFzqjWGzay8mc0ys/Wef8/3tJuZjfD0xUozuyLfd+70fH69md0ZjHPxhplVN7MvzGyVmWWa2UOe9kjsi3gz+8rMvvH0xROe9lpmtsRzzh95ZrbFzEp4trM8P6+Zb19DPO1rzaxzcM7IO2YWbWZfm9lUz3ZE9sNvOOfC6kXemrobgNpAHPAN0CjYcfn4HNsBVwAZ+dqeAwZ73g8GnvW87wp8DhhwJbDE014e2Oj593zP+/ODfW7n2A8XAld43pchb43lRhHaFwaU9ryPBZZ4zlbPR2sAAALLSURBVHE80MvTPgoY4Hk/EBjled8L+MjzvpHnd6YEUMvzuxQd7PMrQn+kAGOBqZ7tiOyHgq9wvML/ea1d59xR4EOge5Bj8innXBqwu0Bzd+Bdz/t3gZvytY9xeRYD5czsQvKWm5zlnNvtnPsRmAUk+z9633HObXfOrfC830feNN1Vicy+cM65/Z7NWM/LAdcCEzztBfviZB9NAK4zM/O0f+icO+Kc+xbIIu93qtgws2rA9cAbnm0jAvvhVMIx4Z9qrd2qQYolkCo757Z73u8AKnven64/wqqfPH+KX07elW1E9oVnGOO/QDZ5/2ltAH5yeWtXwK/P6+dz9vx8D3AB4dEXLwKPACc82xcQmf3wG+GY8COey/ubNGLqbc2sNDARGOSc25v/Z5HUF8654865y8hbbrQl0CDIIQWcmd0AZDvnlgc7llAUjgnfb2vthridnuEJPP9me9pP1x9h0U9mFktesv/AOTfJ0xyRfXGSc+4n4AugNXnDVicXOsp/Xj+fs+fnCcAPFP++aAvcaGbfkTecey0wnMjrh1MKx4Tvt7V2Q9xk4GR1yZ3Ap/na/+CpULkS2OMZ7pgBdDKz8z1VLJ08bcWGZ6z1TWC1c25Yvh9FYl9UNLNynvfnAR3Ju6fxBXCr52MF++JkH90KzPX8NTQZ6OWpXqkF1AO+CsxZeM85N8Q5V805V5O83/25zrk+RFg/nFaw7xr740VeNcY68sYw/xbsePxwfuOA7cAx8sYW7yZv3HEOsB6YDZT3fNaAkZ6+SAea59vPXeTdjMoC/hTs8ypCP1xF3nDNSuC/nlfXCO2LS4GvPX2RAfzD016bvESVBfwfUMLTHu/ZzvL8vHa+ff3N00drgS7BPjcv+qQ9v1TpRGw/5H9pagURkQgRjkM6IiJyCkr4IiIRQglfRCRCKOGLiEQIJXwRkQihhC8iEiGU8EVEIsT/A/XpLuJDLP5WAAAAAElFTkSuQmCC\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "ml = pyamg.smoothed_aggregation_solver(A)\n",
        "res1 = []\n",
        "x = ml.solve(b, tol=1e-12, residuals=res1)"
      ],
      "metadata": {
        "id": "T52zo_Ntd4Po"
      },
      "execution_count": 46,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "plt.plot(residuals)\n",
        "plt.plot(res1)\n",
        "plt.yscale(\"log\")\n",
        "plt.show()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 267
        },
        "id": "S4Ynv7zaeZms",
        "outputId": "b12ccb58-6f7d-4ca5-f3c2-d14008e863c3"
      },
      "execution_count": 48,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD6CAYAAACiefy7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAf10lEQVR4nO3deXDc5Z3n8fdXrcuy7pZ86bAOG2z5AGOBLRMC4RhMwuCdkM1CspPswpSHTKhNZmpri1RqN5Oa2nVmJvfAJOMlDoGd4DCETbzghHA7sWVjm3AZYyPJlzBgS8bGHMbXd//4/dRqyRaW3ZJa+vXnVeVS96NW99O/gk8//fye5/szd0dERKIvK90dEBGRkaHAFxHJEAp8EZEMocAXEckQCnwRkQyhwBcRyRAKfBGRDJGd7g58lIqKCq+rq0t3N0RExpTNmzd3uXtl//YRC3wzmwl8BagAnnD3H53pb+rq6ti0adOw901EJErMbNfp2lOa0jGzFWa2z8xe7te+2My2mVmbmd0B4O5b3f024LPApam8roiInL1U5/DvARYnN5hZDLgLuA5oAm42s6bwdzcAjwCrU3xdERE5SykFvruvAQ70a74EaHP3Dnc/CqwEloSPX+Xu1wGfH+g5zWypmW0ys0379+9PpXsiIpJkOObwq4A9Sfc7gQVmdgXwaSCPjxjhu/tyYDlAc3OzKruJiAyRETtp6+5PA0+P1OuJiEhfw7EO/3WgJul+ddgmIiJpNByBvxGYbmb1ZpYL3ASsGobXERGRs5Dqssz7gVbgfDPrNLNb3f04cDvwKLAVeMDdt6Te1bPw9k5oe3xEX1JEZLRLaQ7f3W8eoH016Vx6+cN54Cfhbw+lrQsiIqNNNGvp+Ml090BEZNSJZuCLiMgpFPgiIhlCgS8ikiEU+CIiGSKagd+0JN09EBEZdaIZ+OWNkJWT7l6IiIwq0Qx8My3NFBHpJ6KBnwWo0KaISLLoBr5G+CIifUQz8LHgh2uULyLSI5qBb+HbUuCLiCRENPB7Rvia1hER6RHtwNeJWxGRhIgGfs+Ujkb4IiI9ohn4aEpHRKS/aAa+TtqKiJwi4oGvEb6ISI+IBr6mdERE+oto4GuELyLSXzQDHy3LFBHpL5qBr5O2IiKnGLHAN7MGM/uJmT04/C+mwBcR6S+lwDezFWa2z8xe7te+2My2mVmbmd0B4O4d7n5rKq93Fh0LfmoOX0QkIdUR/j3A4uQGM4sBdwHXAU3AzWbWlOLrnB2VVhAROUVKge/ua4AD/ZovAdrCEf1RYCUw6IvMmtlSM9tkZpv2799/bh3TKh0RkVMMxxx+FbAn6X4nUGVmcTP7MTDPzL420B+7+3J3b3b35srKynPsgqZ0RET6yx6pF3L3buC2EXkxnbQVETnFcIzwXwdqku5Xh20jR1M6IiKnGI7A3whMN7N6M8sFbgJWDcPrDEyrdERETpHqssz7gVbgfDPrNLNb3f04cDvwKLAVeMDdt6Te1bPpmEb4IiL9pTSH7+43D9C+GlidynOnpCfwtSxTRCQhmqUVEqt0FPgiIj2iGfhapSMicoqIBr5O2oqI9BfRwNdJWxGR/iIa+Brhi4j0F9HA1yodEZH+oh34GuGLiCREMvBP9gzsFfgiIgmRDPxfPf8GAIc/OJrmnoiIjB6RDPzKknEA3P7zzdzXupNDHxxLb4dEREaBSAb+ZdMnAhAvyOG//3oLC/7X43x15R95+MW9HHpf4S8imWnE6uGPqHBZ5nf+/Rz+k53HLzbu4eEX3+BXz+8llmXMqynlivMrueL8CTRNLiYry87whCIiY1+kA9+AudWlzK0u5Zs3zOL5PQd5ett+ntm+n2//bjvf/t12KgrzuPy8Sj42Pc6ixgomFuent+8iIsMkooF/6rLM7FgWzXXlNNeV81+vPZ/9hz9kzfb9PL19P0+8+ha/fK4TgMbK8Vw6rYJFjRW0NMQpKchJxzsQERly0Qz8QVzTtrIojxvnV3Pj/GpOnnReeeMd1rV3sbatm3/b1Mm9rbswg9lTSlg0LRj9X1xXRkFuRA+ZiEReNNPrLKtlZmUZs6tKmF1VwtKPN3L0+Ele6DzI2rYu1rV3s+IPO/iXZzrIiRnzaspYNC3OpdMquKC6lNzsSJ73FpEIinjgn9vGq9zsLC6uK+fiunK+ejW8f/Q4m3a+zdr2Llrbu/nBE6/x/cdfoyA3xsV15SxqDL4BNE0pJqYTwCIySkU88E8MydMV5Gbz8fMq+fh5lQAcev8Y63d0sy78BrDsN68CUJyfzcKGePABMK2C6RMKMdMHgIiMDhEN/OG94lVJQQ7XzprEtbMmAbDv8BFa27tpbe9mXXs3v3vlLQAqCnNpaawIvwHEqS0v0AeAiKRNRAN/ZKtlTijKZ8mFVSy5sAqAPQfep7Uj+ABY29bF/3thLwBVpeNoCcO/pTHO5HBHsIjISIh24KepeFpNeQE15QV8trkGd6ej6z3WtXfT2t7FE1vf4sHNwRLQhorxLAw/ABY2xKkozEtLf0UkM4xY4JtZA/B1oMTdPzO8LzZ6rmlrZjRWFtJYWcifL5zKyZPO1jffSUwBrXp+Lz/fsBuAGZOKwm8AFVxSX07JOO0BEJGhM6jAN7MVwPXAPnefndS+GPgBEAPudvdvDfQc7t4B3GpmD6bW5cEYvVe8ysoyZk0pYdaUEv7isgaOnzjJS68fYl17N+s7urn/2d38dO1OsgxmV5UkPgC0B0BEUjXYBLkHuBO4t6fBzGLAXcA1QCew0cxWEYT/sn5/f4u770u5t4M1zCdth1J2LIt5tWXMqy3jy5+YxofHT/D87oPhFFDvHoDsLOPCmtJw/r+CebWl5OfE0t19ERlDBhX47r7GzOr6NV8CtIUjd8xsJbDE3ZcRfBtInzF8xau87BgLGuIsaIjz19f07gFo7QhWAN35VBs/fLKNvOws5k8tS5wAnltdSk5Mm8BEZGCpzBFUAXuS7ncCCwZ6sJnFgf8JzDOzr4UfDKd73FJgKUBtbe259WwMB35//fcAvHPkGM92HEh8AHz7d9vDx2kTmIh8tBGbFHb3buC2QTxuObAcoLm5+dzmZGz0zuGnqjg/h6ubJnJ1U1Dz/8B7R9kQhv+69i6W/WZ/+LhgE1jPOYDzJmoTmEimSyXwXwdqku5Xh23pF6ER/pmUj8/lujmTuW7OZAD2vXMksQeg/yawBT27gBsrqItrE5hIpkkl8DcC082sniDobwI+NyS9StUIb7waTSYUn34T2Pr2bta2d/HIi8H1fieX5NPS8w1gWgVVpdoEJhJ1g12WeT9wBVBhZp3AN9z9J2Z2O/AowcqcFe6+Zdh6ejYyaIR/Jv03ge3o2QTW0c3T2/fz0B+DL2VT4wWJD4CWxjgTinQhGJGoGewqnZsHaF8NrB7SHg2JsbMscySZGQ2VhTRUFvIfw01g2/cdZl1bMP3zyEtvsHJjcB5+2oTCRA2gBfVxysbnprn3IpKqaO7kGUU7bUezrCxjxqRiZkwq5paP1XPipLNl76HEHoAHN/deCGbmpOJEHaCL68spztcuYJGxJqKBH91VOsMplmWJawDfdnlwIZgXO3s3gd23fhc/+cMOsgzmVJfSEp4EbtYuYJExIZr/l2oOf0jkZvdeB/i/XDWdI8dO8Nzut1kfrgC6+/cd/PiZdnJiwS7gloY4CxvjXFRbpl3AIqNQRANfI/zhkJ8TY1FjcIH3v6F3F3DPSeCeXcC52VnMr+27C1iXghRJv4gGfuYuyxxJp9sFvHHHgcQU0Hce2w6PBbuAm8NdwC0NcWZXlWgXsEgaRDvwNcIfUcX5OVw1cyJXzQx2Ab/93lE27OjdBPat8FKQRXnZLGgoDy8HWcGMSUVk6QNAZNhFO/BPDs01beXclI3PZfHsySyeHe4CPnyE9R0HwmsBdPH41qCAallBTlIZiDiNlSoDITIcoh34mtIZVSYU5XPDBVO44YIpAOw9+EEQ/mEpiN+8/CYAlUV5vZvAGuJMVRkIkSERzcAfxRdAkV5TSsdx4/xqbpxfjbuz58AHtHZ0Jc4BrAqvBTylJD+8FGQFLY1xlYEQOUfRDHxtvBpzzIzaeAG18Vr+w8W1uDvt+98LR/9dPL1tPw89d5oyEA1xJhSrDITIYCjwZVQyM6ZNKGTahN5rAW9763BiCii5DERj5fjE6H9hQ5xylYEQOa2IBr6mdKImK8uYObmYmZN7y0C8svcd1rV30drRzS+f6+S+9buA4GLwPR8Auhi8SC8FvoxJsSxjTnUJc6pL+MvLGzl24iQvdh6iNfwA+NcNu1ixdkfvxeDDKaCL68oZnxfN/+xFziSa/+VrlU7GyYkF1/idP7WM26+czofHT/DH3QfDJaDdrFi7g39ZE1wM/oKa3jpAF01VGQjJHNEOfI3wM1ZedoyFDcGc/l9fAx8cPcGmXQcSm8B+9Ew7dz7VRm4si3m1pYkpoAtrVAZCoivaga+NVxIalxvjsumVXDY9KANx+MixsA5QMAX0/Se2873HYVxOjOa6ssQKoDlVJWTH9AEg0RDtwNeUjgygKD+HT8yYwCdmTADg4PtH2bDjQGIK6B9+uw2AwrxsLqkP6gAtbIjTNLlYZSBkzIpm4GvjlZyl0oJcrp01iWtnTQKg690PWd8RTP+sb+/myVeDMhAl43JY2FAengSu4LyJKgMhY0c0A1/r8CVFFYV5XD93CtfPDcpAvHnoCK0dXYlzAI9ueQuA+PjcRB2glsY4DRXj9QEgo5YCX2QQJpXk82fzqvmzedUA7DnwPq0d3YmLwTzy0hsATCzurQO0qLGCmvKCdHZbpI+IB76mdGR41JQXUFNewGeba3B3dna/n9gF/Ie2Ln71fFAHqKp0XKIKaEtjnMklqgMk6RPRwNccvowcM6O+Yjz1FeP53IKgDlDbvndp7ehmXVs3j299iwc3dwJQFy+gJVwC2tIQp7IoL829l0wyYoFvZjOBrwAVwBPu/qNhfLHgpwJf0sDMmD6xiOkTi/hCSx0nTzqvvnmYde1drO/o5uEX9nL/s7sBmDahMHElsAWqAyTDbFCBb2YrgOuBfe4+O6l9MfADIAbc7e7fGug53H0rcJuZZQH3AsMX+BBO62gOX9IvK8tomlJM05Ri/uKyBo6fOMmWve8krgPw4OZO7m3trQPUM/+vOkAy1AY7wr8HuJMgqAEwsxhwF3AN0AlsNLNVBOG/rN/f3+Lu+8zsBuBLwH0p9vvMLEsjfBmVsmNZXFBTygU1pdyWqAN0MHEO4OcbdvPTtTvJMpg1pSSxAujiunIKVQdIUmA+yJUsZlYHPNwzwjezFuBv3f3a8P7XANy9f9if7rkecfdPDfC7pcBSgNra2vm7du0aVP9O8XeVsPCv4Jpvntvfi6RJnzpAHd08v/sgR0+cJJZlzK3uLQTXPLWccbmqAySnMrPN7t7cvz2V4UIVsCfpfiew4CM6cAXwaSAPWD3Q49x9ObAcoLm5OYU5GUNTOjIW9akDRFAH6LndYRmI9m6Wr+ngn59uJydmzKspY2F4DmBebakKwclHGrHvh+7+NPD0SL1eMKWjwJexb1xujEunVXDptAoA3vvwOBt3HkicA7jzydf44ROvkZcdVAzt+QYwt1qF4KSvVAL/daAm6X512DY6aA5fImp8XjZXnD+BK84P6gAd+uAYG3f0fgB857Ht8FhvIbieSqCzpxSrEFyGSyXwNwLTzayeIOhvAj43JL0aChrhS4YoGZfD1U0TubppIgBvv3eUDTu6E+cA/v63rwJQlJfNxfXliW8AKgSXeQa7LPN+4Aqgwsw6gW+4+0/M7HbgUYKVOSvcfcuw9fRsaYQvGapsfC6LZ09m8ezJAOw/HBSC6ykFkVwIbkF9eWIZqArBRd+gAt/dbx6gfTUfcQI2rQwFvghQWZTHn14whT+94NRCcK0d3fzulb6F4HpOAjdWqhBc1ER3Ua82Xomc1kcVgmvt6C0EN6EoL1ECoqUxTm15gT4AxrhoB75G+CJn1L8Q3K7u91kXhv/atm5+nVQILrkUdFWpCsGNNdEOfF3iUOSsmBl1FeOpSyoE177/3cR1AJ589S1++VxQCG5qvCAx+m9piDOhOD/NvZcziXbga0pHJCVmxrQJRUybUMSfh4Xgtr11ODH//8hLb7ByY7D/sqFyfFgIroKFDeXEC1UJdLSJbuBjmtIRGWJZWcbMycXMnFzMLR+r58RJ55W97yROAv/f517n/6wPKoGeP7EoMf2zsD5OSYEKwaVbdANf6/BFhl0sy5hTXcKc6hKWfjwoBPfS64dobe9mfUc3Kzfu5p51OzGDpsnFtDTEWTQtKARXlK8PgJGmwBeRIZMTy+Ki2jIuqi3jy5+YxofHT/DCnkPhFFAX967fxd1/2EEsy5hdFRSCW9QYp7mujILc6MbRaBHdI6xVOiJpl5cd45L6ci6pL+crTOfIsRM8t+vtRBmIu3/fwY+fCQrBXVBdmjgBfNHUMhWCGwYRDnzN4YuMNvk5MRZNq2BRUiG4TbveTpwEvuupNv7pyTZys7O4qLaUloagDtCFNSoENxQU+CKSNuPzsrn8vEouP68SgMNHjrFx5wHWtQUfAN9/Yjvfexzyc7JonhqUgVjYEGdudQk5KgR31iIc+FqWKTLWFOXncOWMiVw5IygEd/D9o2zYcSBxEvgfH90GwPjcWJ9CcLOmlBBTIbgzinbga4QvMqaVFuRy7axJXDtrEgDd736Y+ABo7ehm2W/CSqD52WEhuApaGuLMmFSkSqCnEe3A105bkUiJF+bxyTmT+eScoBLovneOBHWAwpPAj28NKoGWFeSwoL63DMT0CaoEClEPfE3piETahOJ8llxYxZILqwDYe/CDxOi/tb2b3255E4CKwjwWNpQnVgHVV2RmJdDoBr522opknCml47hxfjU3zk+qBBp+AKxr7+LhF4NKoJOK81nYUJ64GlhNeUE6uz1iohv42nglkvESlUAvDiqB7uh6LzH6/0NbF79KqgTaM/pfNC3O5JJoVgJV4ItIRjAzGioLaags5PMLpuLuvLYvqAQazP+/xYObg0qgdfGCxBLQlsY4E4qiUQk0woGvKR0RGZiZcd7EIs6bWMQXFwWVQF9983D4DSCY/rn/2aAS6LQJhYkloAsb4pSPz01z789NhANfyzJFZPCysoymKcU0TSnm1rAS6Ja9hxLXAvjlc53ct34XADMmFSWmgBaMoUqgCnwRkdOIZRlzq0uZW13KX14eVAJ9sfNQYgnozzfs5qdrg0qgs6YUByeAG+JcXF9OYd7ojNbR2auhoCkdERlCObEs5k8tY/7U3kqgz+8+mDgJfM/anSxf0xGUjK4qoaUxrAQ6tZxxuaOjEFyEA18jfBEZPnnZMRY0xFnQEOerV8ORYyfYnFQI7n+v6eBHTweVQC+sKaWlIc7CxjgX1aavEuiIBb6ZXQH8HbAFWOnuTw/vC2rjlYiMnPycGJdOq+DS01UCbe/izqfa+GFYCXR+bVliF/AF1SNXCXRQgW9mK4DrgX3uPjupfTHwAyAG3O3u3/qIp3HgXSAf6DznHg+WRvgikkb9K4G+c+QYG5PqAH3v8e189zEYlxOjua4ssQR0TtXwVQId7Aj/HuBO4N6eBjOLAXcB1xAE+EYzW0UQ/sv6/f0twO/d/Rkzmwh8F/h8al0/A63DF5FRpDg/h6tmTuSqmb2VQNd3HEicBE6uBNpcV84d181g5uTiIe3DoALf3deYWV2/5kuANnfvADCzlcASd19G8G1gIG8DA17O3syWAksBamtrB9O9gZ5JI3wRGbVKC3JZPHsSi2cHlUC73v2QDR0HEheELxiGE72pzOFXAXuS7ncCCwZ6sJl9GrgWKCX4tnBa7r4cWA7Q3Nx87kN0MzipwBeRsaGiMI9PzZ3Mp+ZOHrbXGLGTtu7+EPDQSL1eMKVzbMReTkRktEvlzMDrQE3S/eqwbXTQSVsRkT5SCfyNwHQzqzezXOAmYNXQdGsIKPBFRPoYVOCb2f1AK3C+mXWa2a3ufhy4HXgU2Ao84O5bhq+rZ0k7bUVE+hjsKp2bB2hfDawe0h4NFY3wRUT6GJntXemgnbYiIn1EO/A1whcRSYh44GuELyLSI7qBDxrhi4gkiW7ga0pHRKSPiAe+pnRERHpEPPA1whcR6aHAFxHJEBEP/BPp7oWIyKgR4cA3zeGLiCSJcOBrp62ISLJoB75G+CIiCREOfFXLFBFJFuHA1yodEZFk0Q18XcRcRKSP6Aa+5vBFRPqIeOBrhC8i0iPiga+NVyIiPSIe+JrSERHpocAXEckQEQ58rdIREUmWPVIvZGaXAZ8PX7PJ3RcN7wuqtIKISLJBjfDNbIWZ7TOzl/u1LzazbWbWZmZ3fNRzuPvv3f024GHgZ+fe5UHSCF9EpI/BjvDvAe4E7u1pMLMYcBdwDdAJbDSzVUAMWNbv729x933h7c8Bt6bQ58HRskwRkT4GFfjuvsbM6vo1XwK0uXsHgJmtBJa4+zLg+tM9j5nVAofc/fA593iwFPgiIn2kctK2CtiTdL8zbPsotwI//agHmNlSM9tkZpv279+fQvc0pSMikmxEV+m4+zfcfd0ZHrPc3ZvdvbmysvLcX0zLMkVE+kgl8F8HapLuV4dto4N22oqI9JFK4G8EpptZvZnlAjcBq4amW0PAwremUb6ICDD4ZZn3A63A+WbWaWa3uvtx4HbgUWAr8IC7bxm+rp4lBb6ISB+DXaVz8wDtq4HVQ9qjoWIW/PSTRHlDsYjIYEU3CfsEvoiIRDjwe96apnRERCATAl8jfBERQIEvIpIxFPgiIhkiuoGPTtqKiCSLbuBrHb6ISB8ZEPga4YuIQEYEvkb4IiIQ6cDXHL6ISLIIB76mdEREkinwRUQyhAJfRCRDKPBFRDJEBgS+rnolIgIZEfhalikiApEOfC3LFBFJFuHAj+5bExE5F9FNRZ20FRHpI8KBrykdEZFkEQ58jfBFRJJFN/BVD19EpI8RC3wzazKzB8zsR2b2meF/QY3wRUSSDSrwzWyFme0zs5f7tS82s21m1mZmd5zhaa4D/sndvwR84Rz7O3hahy8i0kf2IB93D3AncG9Pg5nFgLuAa4BOYKOZrQJiwLJ+f38LcB/wDTO7AYin1u1B0E5bEZE+BhX47r7GzOr6NV8CtLl7B4CZrQSWuPsy4PoBnurL4QfFQwO9lpktBZYC1NbWDqZ7AzyRRvgiIslSmcOvAvYk3e8M207LzOrMbDnBt4R/HOhx7r7c3ZvdvbmysvLce6fAFxHpY7BTOilz952EI/cRoZO2IiJ9pDLCfx2oSbpfHbaNDtp4JSLSRyqBvxGYbmb1ZpYL3ASsGppuDQGN8EVE+hjsssz7gVbgfDPrNLNb3f04cDvwKLAVeMDdtwxfV8+SRvgiIn0MdpXOzQO0rwZWD2mPhopG+CIifUS3tEKiPLJW6YiIQCYE/kltvBIRgUgHfiz4qZ22IiJAlAM/Kwz8k5rDFxGBKAe+aumIiPQR3cBPjPAV+CIiEOXA1xy+iEgf0Q38caXQtASKpqS7JyIio8KIFU8bcSXV8Nl7z/w4EZEMEd0RvoiI9KHAFxHJEAp8EZEMocAXEckQCnwRkQyhwBcRyRAKfBGRDKHAFxHJEOY+ei8QYmb7gV3n+OcVQNcQdmcs07EI6Dj00rEIRPU4THX3yv6NozrwU2Fmm9y9Od39GA10LAI6Dr10LAKZdhw0pSMikiEU+CIiGSLKgb883R0YRXQsAjoOvXQsAhl1HCI7hy8iIn1FeYQvIiJJFPgiIhkikoFvZovNbJuZtZnZHenuz1AzsxVmts/MXk5qKzezx8zstfBnWdhuZvbD8Fi8aGYXJf3NF8PHv2ZmX0zHe0mFmdWY2VNm9oqZbTGzr4TtmXgs8s3sWTN7ITwW3wzb681sQ/ief2FmuWF7Xni/Lfx9XdJzfS1s32Zm16bnHaXGzGJm9kczezi8n5HH4RTuHql/QAxoBxqAXOAFoCnd/Rri9/hx4CLg5aS2fwDuCG/fAfx9ePuTwG8AAxYCG8L2cqAj/FkW3i5L93s7y+MwGbgovF0EbAeaMvRYGFAY3s4BNoTv8QHgprD9x8CXwtt/Bfw4vH0T8IvwdlP4/0weUB/+vxRL9/s7h+PxN8DPgYfD+xl5HPr/i+II/xKgzd073P0osBJYkuY+DSl3XwMc6Ne8BPhZePtnwL9Lar/XA+uBUjObDFwLPObuB9z9beAxYPHw937ouPsb7v5cePswsBWoIjOPhbv7u+HdnPCfA1cCD4bt/Y9FzzF6ELjKzCxsX+nuH7r7DqCN4P+pMcPMqoFPAXeH940MPA6nE8XArwL2JN3vDNuibqK7vxHefhOYGN4e6HhE6jiFX8XnEYxsM/JYhNMYzwP7CD602oGD7n48fEjy+0q85/D3h4A40TgW3wf+G3AyvB8nM4/DKaIY+BnPg++kGbPe1swKgV8CX3X3d5J/l0nHwt1PuPuFQDXBaHRGmrs04szsemCfu29Od19GoygG/utATdL96rAt6t4KpycIf+4L2wc6HpE4TmaWQxD2/+ruD4XNGXkserj7QeApoIVg2io7/FXy+0q85/D3JUA3Y/9YXArcYGY7CaZzrwR+QOYdh9OKYuBvBKaHZ+VzCU7ErEpzn0bCKqBndckXgV8ntX8hXKGyEDgUTnc8CvyJmZWFq1j+JGwbM8K51p8AW939u0m/ysRjUWlmpeHtccA1BOc0ngI+Ez6s/7HoOUafAZ4Mvw2tAm4KV6/UA9OBZ0fmXaTO3b/m7tXuXkfw//6T7v55Muw4DCjdZ42H4x/BaoztBHOYX093f4bh/d0PvAEcI5hbvJVg3vEJ4DXgcaA8fKwBd4XH4iWgOel5biE4GdUG/Od0v69zOA4fI5iueRF4Pvz3yQw9FnOBP4bH4mXgf4TtDQRB1Qb8G5AXtueH99vC3zckPdfXw2O0Dbgu3e8thWNyBb2rdDL2OCT/U2kFEZEMEcUpHREROQ0FvohIhlDgi4hkCAW+iEiGUOCLiGQIBb6ISIZQ4IuIZIj/D4A6/jgHAOtdAAAAAElFTkSuQmCC\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "from scipy.sparse.linalg import LinearOperator\n",
        "\n",
        "ml = pyamg.smoothed_aggregation_solver(A)\n",
        "\n",
        "M = LinearOperator(\n",
        "    matvec=lambda b: ml.solve(b, tol=100),\n",
        "    shape=A.shape, 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": {
        "id": "_SF45lSnetMd"
      },
      "execution_count": 59,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "plt.plot(residuals)\n",
        "plt.plot(residuals2)\n",
        "plt.yscale(\"log\")\n",
        "plt.show()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 265
        },
        "id": "pv2EHUKXfcPr",
        "outputId": "7078b7e5-1ce7-422a-eb09-da4a9b244ade"
      },
      "execution_count": 60,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAdoUlEQVR4nO3deXSc1Znn8e+jzfIqWZatKluSLWNjW65yAIt9M4uXkpnQneGkIcwkJ5Djk3RyJjM90z2ks9BJp5t0ciYdepJOQid0QqYDoRMS3FiyMVsgxAEbCFTJ+4KxTJXkBe/YsqQ7f9RbslBZtqytSm/9PufouOqqllsX/NP18756H3POISIi/peX6QmIiMjwUOCLiOQIBb6ISI5Q4IuI5AgFvohIjijI9ATOpby83M2YMSPT0xARGVFee+21/c65yT3HszrwZ8yYwYYNGzI9DRGREcXMdp9tXCUdEZEcocAXEckRCnwRkRyhwBcRyREKfBGRHKHAFxHJEQp8EZEc4cvAX7fjAD95eVempyEiklWy+hev+uvfN+zhiTf20ungnutqMj0dEZGs4MvA/8Z/XsDxtna+9tRG/rDzAP9zyRzmBMZneloiIhnly5JOUUEe3797IX+5dA7rdhxg2YMv8umfvcYLW1rp6FSHLxHJTZbNLQ7r6urcQK+lc+hEG//y0k4efXUPB4+3MbWkmDvqqvjTS6dRUz52kGYqIpI9zOw151xd2rjfAz/lVHsHz2xs5Rcb9vDStn04Bx+qLOHDl0zjtgVBKiYUD8r7iIhkWs4Hfnfxw+/z1JtxnnxzL7G9RzCDuukTqQ8HiYSCBEoU/iIycinwe7Fj3zFWvRWnIRpnc+IoZrCwOhn+9WGFv4iMPAr8PtjeeoyG6JnwB1g4PRX+AYIlo4dtLiIi/aXAv0A79h2j4a04qxT+IjLCKPAHYOe+5M5/VTTBpvgRAC6rLu0q+0wtVfiLSPZQ4A+SnfuO0RhLsOqtOBu98L+0upTl4SCRcJBpCn8RyTAF/hDYtf94V82/6d1k+F9SlQr/AJUTx2R4hiKSixT4Q+zt/cdZ1SP8P1RVyvJwgPpwUOEvIsNGgT+M3t5/nIZYMvxje73wryzpqvlXlSn8RWToKPAzZPeB4zREEzRE40T3HgYU/iIytBT4WeCdAye6dv5vNSfDf4EX/ssV/iIySBT4WeZs4R+edib8qycp/EWkfzIe+GY2E/giUOKcu6Mvz/Fz4He35+CJrrN93vTCPzRtQlf4T5+kq3qKSN8NKPDN7GHgNqDVORfqNr4MeBDIB37knPtGH17rlwr83u05eILGWPKXvN7ccwiA+VPPhP8MXdJZRM5joIF/A3AMeCQV+GaWD2wFFgPNwHrgLpLh/0CPl7jHOdfqPU+B30fN752gMZpgVTTOH73wrw1OYPmC5AFfXc9fRM5mwCUdM5sBPNUt8K8G/sY5t9S7/wUA51zPsO/5Ogr8fth76H0ao8lr+7zxTjL85wUncJvCX0R6GIrAvwNY5pz7lHf/vwJXOuc+18vzJwF/R/JfBD/q7QeDma0AVgBUV1cv3L17d5/ml0t6C//UL3nNnDwuwzMUkUzKeOD3h3b45/fuofe7Dvi+7oX/3MD4rqt6zpqi5u0iuaa3wC8YwGvuBaq63a/0xmQYTS0dzaeun8mnrp9J/PD7NEYTNMbi/OMzW/n22q3MnjKu65e8Lq4Yh5llesoikiED2eEXkDxoewvJoF8PfMw51zRYk9MOv/9ajpxkdSz5G76vvn0Q5+CiyWO72jjOC45X+Iv41EDP0nkUWASUAy3A/c65H5tZPfAdkmfmPOyc+7vBnLQCf3C0Hj3JmqYWGt6K88quA3Q6qCkfSySUrPnPnzpB4S/iIxn/xav+UOAPvv3HTvF0UwsN0Tjrdh6go9NRXTaGSDjA8nCQ8LQShb/ICKfAlzQHj7fxdFOChliC32/fT3uno3LiaK/sE+CSqlKFv8gIpMCXczp0oo2nN7bQGI3zu+37Od3hmFpSTMQ74HtpVSl5eQp/kZFAgS99dvjEaZ7ZlCz7vLRtP20dnQQmFBPxzvNfWD1R4S+SxRT40i9HTp7m2U0tNEQT/HbrPtraO5kyflTXAd+6GWXkK/xFsooCXwbs6MnTPLe5lYZonBe27ONUeyeTx49i2fwAkXCAK2smKfxFsoACXwbV8VPtPLe5lcZYnOc2t3LydCfl44pYMj9AfSjIVTPLKMjPy/Q0RXKSAl+GzIm2dp7fvK8r/E+0dTBxTCFLapM7/2tnlVOo8BcZNgp8GRbvt3Xw263J8H92UyvHTrVTMrqQxbUV1HvhP6ogP9PTFPE1Bb4Mu5OnO/jdtv00ROOs3dTC0ZPtjB9VwK21FURCAW64eDLFhQp/kcE2FBdPEzmn4sJ8bq2t4NbaCtraO3l5x34ao3Ge3tjCr9/Yy9iifG6eV0F9KMCiOVMYXaTwFxlK2uHLsDvd0ckfdh6gIZpgTVOCg8fbGF2Yz01zJxMJBbl57hTGjtJeRKS/VNKRrNTe0cmrbx/0LuucYP+xU4wqyOPGiydTHw5y87wpTCguzPQ0RUYUBb5kvY5Ox2u736MhGmd1LEHiyEmK8vO4fnY5kXCQxfMqKBmj8Bc5HwW+jCidnY439hyiMRqnMZZg76H3Kcgzrp1VTn04wOLaAGVjizI9TZGspMCXEcs5x1vNh2mIJVs57jn4Pvl5xtUzJxEJB1g6P0D5uFGZnqZI1lDgiy8452h69wiNsTgN0QS79h8nz+CKmjLqw0GWzQ8wZUJxpqcpklEKfPEd5xxbWo7SEE22ctzeegwzqJs+kUgoyLJQgKmlozM9TZFhp8AX39vmhX9jLM7mxFEALq0upd4L/6qyMRmeocjwUOBLTtmx71hXE/emd48AsKCyhEgoSH04wPRJYzM8Q5Gho8CXnLX7wHEaYwkao3HebD4MQG1wAvXhAJFwkIsmj8vwDEUGlwJfBGh+70TXzv/1dw4BMKdifFc3r4srxmd4hiIDp8AX6SF++H1WxxI0RhOs330Q52DWlHHUh5I7/7mB8WriLiOSAl/kHFqPnGRNU4KGaIJXdh2g00FN+diuVo7zp05Q+MuIocAX6aP9x07xdFMLjbE4v99xgI5OR1XZaOpDQSLhIB+qLFH4S1ZT4Iv0w3vH21i7sYWGWJyXt+/ndIdjakkxkXDybJ9LqyaSpz6+kmUU+CIDdPjEaZ7ZlNz5v7h1P20dnVRMGEUkFCQSClA3o0xN3CUrKPBFBtHRk6d5bnMrDdE4L2zZx6n2TsrHjWJZqIL6UJAratTEXTJHgS8yRI6fauf5La00RhM8t7mV9093UDa2iKXzK4iEglx90SQ1cZdhpcAXGQbJJu6tNEQTPLupheNtHZSMLmRJbQX14SDXzJqkJu4y5BT4IsPs5OkOXtqW7OO7dmMLR0+1M764gMXzKoiEg1w/u1xN3GVIqIm5yDArLsxncW0Fi2srONXewe+3H6DBa+L+hNfE/ZZ5FdSHA9x4sZq4y9Abth2+mc0DPg+UA886575/vudohy9+dLqjk3U7DtAYi7OmqaWrifvNc6cQCQe4aY6auMvADKikY2YPA7cBrc65ULfxZcCDQD7wI+fcN/rwWnnAI865/3K+xyrwxe/aOzp5dddBGmJxVsdaupq4L5rjNXGfO4XxauIuF2iggX8DcIxkUIe8sXxgK7AYaAbWA3eRDP8HerzEPc65VjP7MPAZ4GfOuZ+f730V+JJLOjodG94+mLyyZyxOy5FTFOXnccPF5URCQW6traBktMJfzm/AB23NbAbwVLfAvxr4G+fcUu/+FwCccz3D/myvtco5t7yX760AVgBUV1cv3L17d5/m9wFPfhb2vg5/vu7CnyuSBZJN3N9LNnSJxnn38EkK870m7qEgi2srmKgm7tKLoThoOw3Y0+1+M3DlOSawCPgIMApo6O1xzrmHgIcgucPv18za2+D0iX49VSQb5OUZC6eXsXB6GV9aPo83mw/TGI3TEIvzV796i/xfG9dcNIlIKMiS+RVq4i59MmxHhpxzLwAvDMubmUEWn24qciHMjEuqSrmkqpT7InNpevcIDdE4DdE4f/3rKF/6TZQrayZRHw6wVE3c5RwGEvh7gapu9yu9sSyg65mIP5kZoWklhKaV8JdL57A5cZTGaJxV0ThffrKJr6xs4vLpZUTCAZaFAgRL1MRdzhhI4K8HZptZDcmgvxP42KDMalBohy/+ZmbMC05gXnACf7FkzgeauH/1Pzby1f/YyGXVpdSHk03cKyeqiXuu61Pgm9mjwCKg3Myagfudcz82s88Ba0iemfOwc65pyGZ6IcyU95JzZleM5/MV4/n8rbM/0MT966s28fVVm/hQZUnyss6hINWTFP65yJ+XVvj1Z+Dt38H/iA7+pERGmLM1cZ8/dQL14eRlnWeqibvv5OClFbL3B5nIcJo+aSyfvvEiPn3jRew5eMJr5RjnW2u28K01W5gbGE8klGzoMltN3H3Nn4Gvs3REzqqqbAyfun4mn7p+5geauH/n2a384zNbk03cvW5ecyrUxN1v/Bn4GNrhi5xbsGQ0n7y2hk9eW0PrkZOs9nb+331uG//07DZmlo8lEg4QCamJu1/4M/D1/6XIBZkyoZiPXz2Dj189g31HT/H0xuTO/we/3cn3nt9BddkYIuEA9aEgC9TEfcTyZ+CDSjoi/TR5/CjuvnI6d185nYPH21i7MUFDNMGPX9rFD3+7k2mlo1kWCqiJ+wjk08BXSUdkMJSNLeLPLq/mzy6v5vCJ06zd1EJjNM7P1u3mx7/bRcWEUSybHyASDnK5mrhnPX8Gvg7aigy6kjGF3LGwkjsWVn6giftj6/fw03W7KR9XxNL5AerDQa5UE/es5M/AVxFfZEiNLy7k9kumcfsl0zh+qp0XtuyjIRbnidf38m+vvMPEMYUsqQ0QCQe45qJyigoU/tnAp4EPKumIDI+xowpYviDI8gVBr4n7Phpjyev7/GLDHiYUF7C4Nlnzv252uZq4Z5A/A18lHZGMGF2Uz7JQ8sJtJ0938Ltt+2mMJVi7McGvXm9m3KgCbpk3hUgoyKI5k9XEfZj5M/BV0hHJuOLCfG6treDW2gra2sP8fsd+GqMJ1mxM8OQf32VMUT43zZ1CfSjITXMnM6bIp3GURXy8wtrhi2SLooI8Fs2ZwqI5U/h6R4hXdh70mrgnWPVWnOLCPBZdnGzirj6+Q8efga+SjkjWKszP47rZ5Vw3u5yv3R5i/dsHaYzGaYwlWN2UUB/fIeTPwNd5+CIjQn6ecdXMSVw1cxL3/6f5vP5Oso/v6licZza1qo/vIPNn4OvXvkVGnLw8o25GGXUzyvjybWfv43v1zElEvFaO6uN74fwZ+KCSjsgIdq4+vl/8dYwv/ybGFTVl1IeDLJ0foEJ9fPvEp4Gvko6IX/TWx7cxluArTzZx/8om6qZPJBJKtnKcWqo+vr3xZ+CrpCPiS2fr49votXL82lMb+dpTG7mkqpR677LOVWVq5didPwMfVNIRyQGzK8Yzu2I8/+2W2ezcdyzZyjEW5+8bNvP3DZsJTyvpuqZ/TfnYTE8343wa+CrpiOSamZPH8dmbZvHZm2ax5+AJGmNxGqIJvrl6C99cvYV5wQnUh5LX95k1JTdbOfoz8FXSEclpVWVjWHHDRay44SL2Hkq1cozz7We28n/WbmX2lHFEcrCVoz8DH7TBFxEAppWO5t7rarj3uhpajpzs+u3eVCvHmvKxRELJsk9omr9bOfo08FXSEZF0FWdp5bg6luCHL+7kn1/YQeXE0dSHg0RCAS6pKvVd+Psz8H32H0lEBl/3Vo7vHW9j7cYWGmNx/vXlXTz04k6mlhSzLBQkEg6wsNofrRz9Gfigs3REpM8mji3io5dX8dHLqzj8/mme3dRCQzTB/3tlNw+/vIsp40exzCv7XFEzcls5+jTwR+Z/DBHJvJLRhXzksko+ctmZVo6rYwke37CHR7xWjkvmB4iEAlw1cxKFI6iVo08DH1TDF5GB6t7K8USb18oxGuc3b+zl56+8Q+mYQpbUVhAJB7l2BLRy9Gfg6/LIIjLIxhQVUB8OUh8OcvJ0By9u3Zf8Ra9ogsc3NDO+uIDFtRVEQkGun12eld28/Bn4IiJDqLgwnyXzAyyZH+BUewcvb99PQzTB000Jnnh9L2OL8rllXgX14QA3XjyF0UXZEf4+Dnzt8EVk6I0qyOfmuRXcPLeCtj8Ns27nAVbH4qxpamHlm+8yujCfm+dOYVko2c1r7KjMxa4/A18lHRHJgKKCPG68eDI3XjyZv729k1d3HaQhFmd1rIVV0TijvO/Xh4PcPG8KE4a5leOwBb6ZLQL+FmgCHnPOvTCE7zZ0Ly0i0gcF+XlcM6uca2aV89UPh3ht93s0ROOsjiV4emMLRV6rx0gowOLaCkrHDH03rz4Fvpk9DNwGtDrnQt3GlwEPAvnAj5xz3zjHyzjgGFAMNPd7xn2mHb6IZIf8POOKmjKuqCnjK7fV8saeQ13X9H9ucysFecY1s8qp98J/0hB18zLXh9KHmd1AMqwfSQW+meUDW4HFJAN8PXAXyfB/oMdL3APsd851mlkF8G3n3N3ne9+6ujq3YcOGC/g4nrVfgVd+CF9qufDniogME+cc0b2HaYgmL+u8+8AJ8gyumjmJLy2vpXbqhH69rpm95pyr6znepx2+c+5FM5vRY/gKYLtzbqf3Bo8BtzvnHiD5r4HevAf0+uPLzFYAKwCqq6v7Mr3eJt3/54qIDAMzY0FlKQsqS/nfy+awMX6ERi/8J4we/Ir7QF5xGrCn2/1m4MreHmxmHwGWAqXAd3t7nHPuIeAhSO7w+zc1XTxNREYWM2P+1BLmTy3hfy2dMyTvMWwHbZ1zTwBPDMub6eJpIiJpBvJ7wHuBqm73K72x7KCSjojIBwwk8NcDs82sxsyKgDuBlYMzrYFSSUdEpKc+Bb6ZPQqsA+aYWbOZ3eucawc+B6wBNgGPO+eahm6qF0AlHRGRNH09S+euXsYbgIZBndFgUUlHROQDsvtanv2mHb6ISE8+DXxQDV9E5IP8Gfi6eJqISBp/Br5KOiIiaXwa+KCSjojIB/kz8HVapohIGn8Gvko6IiJpfBr4Hh24FRHp4s/AT5V0FPgiIl38Gfgq6YiIpPFp4Kdohy8ikuLPwFdJR0QkjT8DXyUdEZE0Pg38FO3wRURS/Bn4qQ2+SjoiIl38Gfgq6YiIpPFp4Kdohy8ikuLPwNdZOiIiafwZ+CrpiIik8Wngp2iHLyKS4s/AV0lHRCSNPwNfJR0RkTQ+DfwU7fBFRFL8Gfgq6YiIpPFn4KukIyKSxqeBn6IdvohIij8DXyUdEZE0/gx8lXRERNL4NPBTtMMXEUnxZ+CrpCMiksafgX/mgvgZnYWISDYpGK43MrPrgbu996x1zl0zhG+W/FM7fBGRLn3a4ZvZw2bWamaxHuPLzGyLmW03s/vO9RrOuZecc58GngJ+2v8p94UO2oqI9NTXHf5PgO8Cj6QGzCwf+B6wGGgG1pvZSiAfeKDH8+9xzrV6tz8G3DuAOZ+fdvgiImn6FPjOuRfNbEaP4SuA7c65nQBm9hhwu3PuAeC2s72OmVUDh51zR3t7LzNbAawAqK6u7sv0zvYqqZn38/kiIv4zkIO204A93e43e2Pnci/wr+d6gHPuIedcnXOubvLkyf2bmXb4IiJphu2gLYBz7v7hfD/t8EVEzhjIDn8vUNXtfqU3lnna4YuIpBlI4K8HZptZjZkVAXcCKwdnWgOlGr6ISE99PS3zUWAdMMfMms3sXudcO/A5YA2wCXjcOdc0dFO9ANrhi4ik6etZOnf1Mt4ANAzqjAaFdvgiIj3589IK2uGLiKTxaeCnPpYCX0QkxZ+BnyrpuM7MTkNEJIv4M/BV0hERSePPwNdBWxGRNP4MfO3wRUTS+DPwRUQkjU8DXyUdEZGe/Bn4KumIiKTxZ+Cr45WISBp/Br52+CIiafwZ+Krhi4ik8Wfgm0o6IiI9+TPwU1TSERHp4u/AV0lHRKSLPwNfB21FRNL4M/B10FZEJI0/A187fBGRNP4MfP3ilYhIGp8Gfop2+CIiKf4MfJV0RETS+DPwddBWRCSNPwNfO3wRkTT+DHzt8EVE0vgz8HUtHRGRNP4M/BSVdEREuvg08FXSERHpyZ+Br4O2IiJp/Bn4+k1bEZE0Pg38FO3wRURShi3wzazWzB43s++b2R1D/GbJP5X3IiJd+hT4ZvawmbWaWazH+DIz22Jm283svvO8TAT4v865zwAf7+d8+0gHbUVEeiro4+N+AnwXeCQ1YGb5wPeAxUAzsN7MVgL5wAM9nn8P8DPgfjP7MDBpYNM+Dx20FRFJ06fAd869aGYzegxfAWx3zu0EMLPHgNudcw8At/XyUp/1flA80b/p9pUO2oqI9NTXHf7ZTAP2dLvfDFzZ24O9Hxh/DYwFvnWOx60AVgBUV1cPYHqgko6IyBkDCfwL4px7Gy/Iz/O4h4CHAOrq6vqX2F0lfAW+iEjKQM7S2QtUdbtf6Y1lAR20FRHpaSCBvx6YbWY1ZlYE3AmsHJxpDZAuniYikqavp2U+CqwD5phZs5nd65xrBz4HrAE2AY8755qGbqr9oJKOiEiXvp6lc1cv4w1Aw6DOaFCopCMi0pM/L62g8/BFRNL4M/B1Hr6ISBqfBn6KdvgiIin+DHyVdERE0vgz8HXQVkQkjT8DX+fhi4ik8Wfgp6ikIyLSxaeBr5KOiEhP/gz8seVQ+ycwdnKmZyIikjWG7WqZw2ryHPjoTzM9CxGRrOLPHb6IiKRR4IuI5AgFvohIjlDgi4jkCAW+iEiOUOCLiOQIBb6ISI5Q4IuI5AhzWXy9GTPbB+zu59PLgf2DOJ2RTGuRpHU4Q2uR5Nd1mO6cS7vUQFYH/kCY2QbnXF2m55ENtBZJWocztBZJubYOKumIiOQIBb6ISI7wc+A/lOkJZBGtRZLW4QytRVJOrYNva/giIvJBft7hi4hINwp8EZEc4cvAN7NlZrbFzLab2X2Zns9gM7OHzazVzGLdxsrMbK2ZbfP+nOiNm5n9k7cWb5nZZd2e8wnv8dvM7BOZ+CwDYWZVZva8mW00syYz+7w3notrUWxmr5rZm95afNUbrzGzV7zP/AszK/LGR3n3t3vfn9Httb7gjW8xs6WZ+UQDY2b5ZvaGmT3l3c/JdUjjnPPVF5AP7ABmAkXAm0Btpuc1yJ/xBuAyINZt7JvAfd7t+4B/8G7XA40kG/1eBbzijZcBO70/J3q3J2b6s13gOgSBy7zb44GtQG2OroUB47zbhcAr3md8HLjTG/8B8Bnv9p8DP/Bu3wn8wrtd6/2dGQXUeH+X8jP9+fqxHn8B/Bx4yrufk+vQ88uPO/wrgO3OuZ3OuTbgMeD2DM9pUDnnXgQO9hi+HUj1dfwp8Cfdxh9xSX8ASs0sCCwF1jrnDjrn3gPWAsuGfvaDxzkXd8697t0+CmwCppGba+Gcc8e8u4XelwNuBn7pjfdci9Qa/RK4xczMG3/MOXfKObcL2E7y79SIYWaVwHLgR959IwfX4Wz8GPjTgD3d7jd7Y35X4ZyLe7cTQIV3u7f18NU6ef8Uv5TkzjYn18IrY/wRaCX5Q2sHcMg51+49pPvn6vrM3vcPA5Pwx1p8B/groNO7P4ncXIc0fgz8nOeS/ybNmfNtzWwc8CvgvzvnjnT/Xi6thXOuwzl3CVBJcjc6N8NTGnZmdhvQ6px7LdNzyUZ+DPy9QFW3+5XemN+1eOUJvD9bvfHe1sMX62RmhSTD/t+cc094wzm5FinOuUPA88DVJMtWBd63un+urs/sfb8EOMDIX4trgQ+b2dsky7k3Aw+Se+twVn4M/PXAbO+ofBHJAzErMzyn4bASSJ1d8gngyW7jH/fOULkKOOyVO9YAS8xsoncWyxJvbMTwaq0/BjY5577d7Vu5uBaTzazUuz0aWEzymMbzwB3ew3quRWqN7gCe8/41tBK40zt7pQaYDbw6PJ9i4JxzX3DOVTrnZpD8u/+cc+5ucmwdepXpo8ZD8UXybIytJGuYX8z0fIbg8z0KxIHTJGuL95KsOz4LbAOeAcq8xxrwPW8tokBdt9e5h+TBqO3AJzP9ufqxDteRLNe8BfzR+6rP0bVYALzhrUUM+Io3PpNkUG0H/h0Y5Y0Xe/e3e9+f2e21vuit0RYgkunPNoA1WcSZs3Rydh26f+nSCiIiOcKPJR0RETkLBb6ISI5Q4IuI5AgFvohIjlDgi4jkCAW+iEiOUOCLiOSI/w9w+maFKBihVQAAAABJRU5ErkJggg==\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "lc2f6bvIfySU"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}