notebooks/subspace-correction-mg.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# Robust multigrid for IGA based on subspace correction\n",
    "\n",
    "\n",
    "This notebook implements some basic 1D, 2D and 3D examples which demonstrate the\n",
    "robust multigrid method for Isogeometric Analysis described in the paper\n",
    "\n",
    "[C. Hofreither and S. Takacs. Robust Multigrid for Isogeometric Analysis Based on Stable Splittings of Spline Spaces.](http://dx.doi.org/10.1137/16m1085425)\n",
    "\n",
    "Only the two-grid method on the reference domain is implemented.\n",
    "It uses the function `twogrid`, a generic two-grid method, as well\n",
    "as the abstract `SubspaceOperator`, which is used as a\n",
    "subspace correction smoother as described in the paper.\n",
    "\n",
    "A crucial component of the method is the computation of a basis for the\n",
    "space $\\widetilde S$ of splines with vanishing odd derivatives at the\n",
    "boundary (denoted by $S_0$ in the paper). This is done using the\n",
    "function `Stilde_basis` from `pyiga.stilde`.\n",
    "\n",
    "There are also some simple tests for the Dirichlet BC case. The method\n",
    "we use here is not the one described in the paper which involves\n",
    "adaptation of the space $S_0$, but a simple algebraic approach based\n",
    "on truncating the subspace prolongation matrix: the first and last row,\n",
    "corresponding to the Dirichlet dofs, are simply dropped. While we don't\n",
    "currently have a proof for this approach, it seems to work very well in\n",
    "practice and is easy to implement.\n",
    "\n",
    "This notebook depends on the [`pyiga` toolbox](https://github.com/c-f-h/pyiga)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.sparse\n",
    "from scipy.sparse import kron as spkron\n",
    "\n",
    "from pyiga import bspline, assemble\n",
    "from pyiga.operators import KroneckerOperator, make_solver, SubspaceOperator\n",
    "from pyiga.solvers import *\n",
    "from pyiga.stilde import Stilde_basis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "%pylab is deprecated, use %matplotlib inline and import the required libraries.\n",
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "plt.style.use('bmh')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def random_rhs(A):\n",
    "    return A.dot(rand(A.shape[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# 1D Multigrid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coarse dofs: 107\n",
      "Fine dofs: 207\n"
     ]
    }
   ],
   "source": [
    "p = 7\n",
    "nspans_c = 100\n",
    "\n",
    "kv_c = bspline.make_knots(p, 0.0, 1.0, nspans_c)\n",
    "kv = kv_c.refine()\n",
    "h = 1.0 / kv.numspans\n",
    "\n",
    "M = assemble.mass(kv)\n",
    "K = assemble.stiffness(kv)\n",
    "A = 1.0*M + K\n",
    "\n",
    "P = bspline.prolongation(kv_c, kv)\n",
    "\n",
    "print('Coarse dofs: %d\\nFine dofs: %d' % (P.shape[1], A.shape[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# construct prolongator for S-tilde and its orthogonal complement\n",
    "(P_tilde, P_compl) = Stilde_basis(kv)\n",
    "\n",
    "# basis for M-orthogonal complement of S-tilde\n",
    "P_tilde_orth = solve(M.A, P_compl)  # M^-1 * P_compl\n",
    "\n",
    "A_tilde_orth = P_tilde_orth.T.dot(A.A).dot(P_tilde_orth)\n",
    "\n",
    "M_tilde = P_tilde.T.dot(M.A).dot(P_tilde)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8 iterations\n"
     ]
    }
   ],
   "source": [
    "sigma = h**-2 / 0.09\n",
    "\n",
    "smoother1 = OperatorSmoother(SubspaceOperator(\n",
    "    [P_tilde,                             P_tilde_orth],\n",
    "    [make_solver((1 + sigma) * M_tilde),  make_solver(A_tilde_orth)]\n",
    "))\n",
    "smoother2 = GaussSeidelSmoother()\n",
    "smoother3 = SequentialSmoother((smoother1, smoother2))\n",
    "\n",
    "twogrid(A, random_rhs(A), P, smoother3);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "31 iterations\n"
     ]
    }
   ],
   "source": [
    "# Dirichlet test\n",
    "sigma = h**-2 / 0.09\n",
    "\n",
    "K_dir = K[1:-1,1:-1]\n",
    "smoother = OperatorSmoother(SubspaceOperator(\n",
    "    [P_tilde[1:-1],                       P_tilde_orth[1:-1]],\n",
    "    [make_solver((1 + sigma) * M_tilde),  make_solver(A_tilde_orth)]\n",
    "))\n",
    "twogrid(K_dir, random_rhs(K_dir), P[1:-1], smoother);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "Linv = P_tilde_orth .dot (inv(A_tilde_orth)) .dot (P_tilde_orth.T) + \\\n",
    "       P_tilde .dot(inv((1 + sigma) * M_tilde)) .dot (P_tilde.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAGYCAYAAABs2slQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABSQklEQVR4nO3de3Rc5X0++mc0oxndRpexLqP7xbJkGYHBIsU2Jph4hdZZAQw557hpC46Lc2rcuHFMTkrLb9EsrzRkNQ0/zio1CaEJpGmLD10x7WppiaOAjSFOg7EdjG6WdZdH9/tImtHM7POH0VhjSVszo5nZ+zt6Pmt5YY22tl7vZ++ZL+9+33cbFEVRQERERKRDCVo3gIiIiGg5LFSIiIhIt1ioEBERkW6xUCEiIiLdYqFCREREusVChYiIiHSLhQoRERHpFgsVIiIi0i2T1g1YDZ/PB7fbDaPRCIPBoHVziIiIKAiKosDr9cJsNiMhQb3PRHSh4na78d5772ndDCIiIgrD3XffjaSkJNVtRBcqRqMRAFBbW+v/eySNjIzAZrNFfL8UPcxMJuYmDzOTSS+5eb1eXL58OajPbtGFyvztHqPRGJVCJSEhISr7pehhZjIxN3mYmUx6yy2YYRscTKuiv79f6yZQiJiZTMxNHmYmk8TcWKgQERGRbrFQUVFZWal1EyhEzEwm5iYPM5NJYm4sVFT09vZq3QQKETOTibnJw8xkkpgbCxUVMzMzWjeBQsTMZGJu8jAzmSTmxkJFxUpzu0l/mJlMzE0eZiaTxNxYqKgoKSnRugkUImYmE3OTh5nJJDE3FioqWlpatG4ChYiZycTc5GFmMknMTZcLvv3Hf/wHnnzySfh8Pvz5n/85Dhw4oHWTiIiIxPMpCi5dm8TQtBvjMx5YLSZMum78NyMpEbaU66XByMwcxmeuv5admohaexqMCbF/rp7uChWPx4OjR4/i7bffRnp6OrZs2YJHHnlEkyV/c3JyYv47aXWYmUzMTR5mJs/Z9jG8cMGLkdnWsH4+OzURh7YWYUd5ZmQbtgLd3fr5n//5H9xyyy0oLCyE1WrF5z73Obz11luatMVk0l0dRytgZjIxN3mYmSxn28dwrL4dI7PesPcx5JzDsfp2nG0fi1zDghDxQuXMmTN44IEHUFBQAIPBgDfeeGPRNsePH0d5eTmSkpJQV1eHd9991/+9a9euobCw0P91UVGRZvO+HQ6HJr+XwsfMZGJu8jAzObw+BcfP9URsfy+e64HXp0RsfyuJeEnsdDqxefNm7N+/H1/4whcWff/EiRM4cuQIjh8/jrvvvhs/+MEPsHv3bjQ0NKCkpASKsvgfv9JDiyYnJ5GQcKPmslgssFgsq//HEBER6YDXp+By3xRGZuaQmRQ4huTmcSY3jze5OjyDIedcxNoy6JzD5b4pbC6wRmyfaiJeqOzevRu7d+9e9vvPPfccHn/8cf8A2eeffx5vvfUWXnzxRTz77LMoLCwM6EHp6enBXXfdpfo7a2trMT097f96//79OHz4MPLz83H16lUAQF5eHhRFwcDAAABgw4YN6OnpwczMDJKSklBcXIwrV64AAHJzc5GQkACv14uGhgasX78efX19cDqdsFgsKCsrQ3NzMwAgOzsbZrMZ165dAwCUl5djcHAQU1NTSExMRGVlJRobGwEANpsNycnJ/n9fWVkZRkZGMDExAaPRiOrqajQ2NkJRFGRmZsJqtaK7uxvA9SllExMTGBsbg8FgQE1NDZqbm+H1epGeno6srCx0dnYCuN4LNT09jZGREQDApk2b0NLSAo/HA6vViuzsbLS3twMACgoK4HK5MDw8DADYuHEj2tra4Ha7kZqairy8PLS1tQEA8vPz4fF4MDg4CACoqqpCV1cXZmdnkZycjMLCQrS2tvqPN3DjAViVlZXo7e31H++SkhL/6POcnByYTCb//6FVVFSgv78fTqcTZrMZFRUVaGpqAgCsW7cOFosl4HgPDQ1hcnISJpMJFRUVaGho8B/vlJQU9PRc/z+J0tJSjI6OLnu809PT0dXVBQAoLi7G5OTkssfbZrOho6MDAFBYWIiZmRn/8a6pqUFrayvm5uaQlpaGnJycgOPtdrsxNDQEAKiurkZHRwdcLhdSU1Nht9v956zdbofP5ws4Z7u7u/3Hu6ioKOCcNRgM/uO9fv16OBwOTE9Pw2KxoLS0VPV4DwwMYGpqasnjnZSUtOQ5azKZUFVV5T/eWVlZSEtLCzhnx8fHMT4+joSEBGzcuBFNTU3w+XzIyMhARkaG/3jn5ubC4XBgdHR00Tm71PGenZ1d8pxNS0tDbm6u6jnb2dkJl8uFlJSUVb9H9PX1+Y/3WnuP8Hq96OzsFPcesfCc1eN7hD0/H791TKFzYBROD7C+KB/vX+nDpWEPpj3Qjd9e6UDimDHs94ilOiWWY1BC2TpEBoMBJ0+exJ49ewAAbrcbKSkpeP311/Hwww/7t/vqV7+Kixcv4vTp0/B4PKipqcE777zjH0x77tw5rFu3btH+PR4PTp8+jYqKiqj0qHR2dqK0tHTV+6HYYWYyMTd5mFn45ntH5mfezPd8XO6bwhsNQ5h0hT+OJFa++7nKVfWoeL1eXLp0Cffee++K451iOhpqaGgIXq/XX0nPy8vL8/+ficlkwve+9z3cd9998Pl8+MY3vrFkkbKQ1WqF0WiMeHudTmfE90nRxcxkYm7yMLOV3VyQWC0mXHJM4ledE5h0678YWU7OJ1OVY0WTYds3jzlRFCXgtQcffBAPPvhg0PvbtWsXZmdnceDAgYiuuWI2myO2L4oNZiYTc5OHmV23VDEy6fKgb3IOb7eNYHwVs2z06omtRTFdTyWmhUp2djaMRqO/92TewMDAol6WUNTX10elR6WioiLi+6ToYmYyMTd51kJmKw1gjYfekVDkpCbiCQ3WUYlpoWI2m1FXV4dTp04FjFE5deoUHnrooVg2JShNTU3YtGmT1s2gEDAzmZibPPGUWbzeognWAzXrsDEnde2sTDs1NeUf1Q0A7e3tuHjxImw2G0pKSnD06FE8+uijuPPOO7Ft2za89NJL6OrqwsGDB8P+ndG69UNERPEhHgawRkNOaiIObSvWpAAJVsQLlQ8++AD33Xef/+ujR48CAPbt24dXXnkFe/fuxfDwMI4dOwaHw4Ha2lq8+eabqxo9Hq1bPysN4iX9YWYyMTd59JrZWu8dCVWsx5uEI+KFys6dO1ecH33o0CEcOnQo0r864rhonDzMTCbmJo+Wma3FAayRptV4k3DwYQ0qrl27hszMTK2bQSFgZjIxN3likRl7RxazWox4aFM2brWnBb0yLXBjvMnUSD82byjXbLxJOOKiUOEYFSIimdg7srQ0cwK2l2Zgc77VX3REYkBrQ8MwNsVo6ftIiYtCJVpjVMrLyyO+T4ouZiYTc5MnlMyWm+Z7oXdt947MFyN3FKT7ez7GZj2wJUdvho3Eay0uCpVoGRoaQnFxsdbNoBAwM5mYmzxLZcZbNYtlJBnxmYos5FnNAbdktJruK/FaY6GiYnJyUusmUIiYmUzMTY75YuTi1VFYx8yc5vuJaN2qiTSJ11pcFCrRGqOy0oOSSH+YmUzMTX+C6h1pvaZtI2NEbQCrXguS5Ui81uS1eAnRGqNSVVUV8X1SdDEzmZibNrgIWiApvSKrIfFai4tCJVoaGhriZonotYKZycTcYmNhYbJWB7JqMYBVTyReayxUiIjiCHtJ9DeAlVYnLgqVaI1RsdlsEdsXxQYzk4m5hY4zbNbGrZpIk3itxUWhEq0xKikpKRHfJ0UXM5OJuam7uShZS4uhsXcksiRea3FRqERLT0+PuHt5ax0zk4m5XbeWe0nYOxIbEq81FipERBpYi70kS03zZUFCK2GhoqK0tFTrJlCImJlM8ZzbWuwlYe+Ifkm81lioqBgdHUVqaqrWzaAQMDOZ4iW3m59pE88zbVJMwI7yrDU7zVcqiddaXBQq0Zr1MzExEbF9UWwwM5mk5baWekmW6x0xjnah9pYyrZtHIZJ2rQFxUqhEa9ZPNPZJ0cXMZNJzbmuhlyScRdCaJ+Pi42PN0fO1thyeaSqqq6u1bgKFiJnJpKfc1sLqrQsLk3DHjugpMwqexNxYqKhobGxETU2N1s2gEDAzmbTMLV4Lk2gvFc9rTSaJubFQUaEoitZNoBAxM5mindtSY0omXfKnBGu5GBqvNZkk5sZCRUVmZqbWTaAQMTOZIp3bWugl0Xq6L681mSTmxkJFRXp6utZNoBAxM5lWk1u8DXaVsmQ8rzWZJObGQkVFV1eXuKWG1zpmJlMoucVbb4meeklCwWtNJom5xUWhEq11VIhIe/FSmHC1VqLwxEWhEq11VIqLiyO+T4ouZibTfG7x8vybhc+0idfVWnmtySQxt7goVKJlcnISVqtV62ZQCJiZLPOFyZXefnRMD4vrLVnLvSS81mSSmBsLFRVjY2MoKCjQuhkUAmamX9IHva6FXpJQ8FqTSWJuLFRUGAxr8w1IMmamH9LHlkgd5BorvNZkkpgbCxUV0lbvI2amB16fgn+52IeTHw+K6S1ZOCV4Ld2+WQ1eazJJzI2Fiorm5maRz0VYy5hZbC11O+f1jwYx6/Fp3bQlreUxJZHGa00mibmxUFHh9cr4v0G6gZlFl7TbObx9Ez281mSSmBsLFRUSV/Bb65hZZEkqTDjYNbZ4rckkMTcWKipsNpvWTaAQMbPIkDDOhL0l2uK1JpPE3OKiUInWyrQdHR3ilhpe65hZ6JZaZO3nV4YxPaefcSYc7Ko/vNZkkphbXBQq0VqZligeSbidw94SIpoXF4VKtBQWFmrdBAoRM1uenm/npJqAu8uzWJgIwmtNJom5sVBRMTMzg4yMDK2bQSFgZtfp+XbOUoNecwxTKMjP17ppFAJeazJJzI2FioqRkRHY7Xatm0EhWOuZ6bHXJMkIfLpCvbekoaGbhYowa/1ak0pibixUiAS7ebzJmfZx3Sy2ZrUY8fAtOfji7XbexiGisLFQUSFxqeG1bi1kNl+c/KprHPWtIxif1b7nZLVrmKyF3OINM5NJYm4sVFS0trZiw4YNWjeDQhCPmelxlk6kZ+XEY27xjpnJJDE3Fioq5ubmtG4ChSieMtPTeJNoTxeOp9zWCmYmk8TcWKioSEtL07oJFCLJmeltvElGkhG71tuwrTQj6tOFJee2VjEzmSTmxkJFRU5OjtZNoBBJy0xP4020XGRNWm7EzKSSmBsLFRXt7e3ilhpe66RkpofbOnpa/VVKbnQDM5NJYm66LFQefvhhvPPOO9i1axf+9V//VevmEK2aXm7rxPJ2DhFRJOiyUPmzP/sz/PEf/zFeffVVTdtRUFCg6e+n0OktM617ToJZbE0P9JYbrYyZySQxN10WKvfddx/eeecdrZsBt9utdRMoRHrIbOG4k/9u1mbZemmLrekhNwoNM5NJYm4hFypnzpzBd7/7XZw/fx4OhwMnT57Enj17ArY5fvw4vvvd78LhcOCWW27B888/j3vuuSdSbY6ZoaEh5Obmat0MCoEWmWm9zklKYgJ+d4MNeVYzMpISdd1zshxea/IwM5kk5hZyoeJ0OrF582bs378fX/jCFxZ9/8SJEzhy5AiOHz+Ou+++Gz/4wQ+we/duNDQ0oKSkBABQV1cHl8u16Gd//vOfi+yWorVJ69s60npNiIjCYVAURQn7hw2GRT0qd911F7Zs2YIXX3zR/1pNTQ327NmDZ599Nuh9v/POO3jhhRdUB9N6PB6cPn0aFRUVSEhI8L9usVhgsVhC+8cswev1wmg0rno/FDuxyGy+QPn/fjsQ0wGxUsabhIPXmjzMTCa95Ob1enHp0iXce++9MJnU+0wiOkbF7Xbj/PnzeOqppwJev//++/H+++9H8lcFqK2txfT0tP/r/fv34/Dhw8jPz8fVq1cBAHl5eVAUBQMDAwCADRs2oKenBzMzM0hKSkJxcTGuXLkCAMjNzUVCQgKuXr2K1NRUrF+/Hn19fXA6nbBYLCgrK0NzczMAIDs7G2azGdeuXQMAlJeXY3BwEFNTU0hMTERlZSUaGxsBADabDcnJyejt7QUAlJWVYWRkBBMTEzAajaiurkZjYyMURUFmZiasViu6u7sBACUlJZiYmMDY2BgMBgNqamrQ3NwMr9eL9PR0ZGVlobOzEwBQVFSE6elpjIyMAAA2bdqElpYWeDweWK1WZGdno729HcD1gVUulwvDw8MAgI0bN6KtrQ1utxupqanIy8tDW1sbACA/Px8ejweDg4MAgKqqKnR1dWF2dhbJyckoLCxEa2ur/3gDQH9/PwCgsrISvb29/uNdUlKClpYWANfn9ZtMJjgcDgBARUUF+vv74XQ6YTabUVFRgaamJgDAunXrYLFYAo730NAQJicnYTKZYDQa/b11NpsNKSkp6OnpAQCUlpZidHR02eOdnp6Orq4uAEBxcTEmJyf9x7uqeiP+6zeNuDQ0h18P+BDL5U5SjMCnCxJwaGcNrvX2YHa2D8muZPi8SWhuunHOGgwG//Fev349HA4HpqenYbFYUFpaqnq8BwYGMDU1teTxTkpKWvKcNZlMqKqqQkNDAwAgKysLaWlpAefs+Pg4xsfHkZCQgI0bN6KpqQk+nw8ZGRnIyMjwH2+DwYDMzEyMjo4CCDxn09PTYbPZ0NHRAQAoLCzE7OzskudsWloacnNzVc/Zzs5OuFwupKSkrPo9oq+vz3+819p7xOjoKHJzc8W9Ryw8ZyP5HnHz8b75nJ2ZmfEf75qaGrS2tmJubg5paWnIyckJON5utxtDQ0MAgOrqanR0dMDlciE1NRV2u91/ztrtdvh8voBztru723+8i4qKAs5Zg8GAtrY2/+ealu8RofSRRLRH5dq1aygsLMR7772H7du3+7f79re/jVdffdV/4a7kd3/3d/Hhhx/C6XTCZrPh5MmT+NSnPrVou2j3qDQ0NIibb77WRTqzWN/eiYfxJuHgtSYPM5NJL7lp1qMyz2AIfFNVFGXRa2reeuutkH6f1WqNSldWampqxPdJ0bXazLRa7yTZZMD/eVvemh1vwmtNHmYmk8TcIlqoZGdnw2g0+rtD5w0MDPi7+aJh165dmJ2dxYEDB3DgwIGI7ddut0dsXxQb4Wam1cBYDoi9jteaPMxMJom5RbRQMZvNqKurw6lTp/Dwww/7Xz916hQeeuihSP6qAPX19VHpUbl69aouusgoeKFmFuuBsXpatl5PeK3Jw8xkkphbyIXK1NSUfzAUcP25ARcvXoTNZkNJSQmOHj2KRx99FHfeeSe2bduGl156CV1dXTh48GBEG04UrlgvyJaSmIDfq1rHZeuJiMIQcqHywQcf4L777vN/ffToUQDAvn378Morr2Dv3r0YHh7GsWPH4HA4UFtbizfffBOlpaWRa/VNeOuH5qllFuvbO7ytEzxea/IwM5kk5hZyobJz584VpxUdOnQIhw4dCrtRoYrWrR+fL/ZLn9PqLJVZrG7vxPM6J9HGa00eZiaTxNx0+awfvRgYGEB2drbWzaAQLMwsVgUKe05Wj9eaPMxMJom5xUWhEq1bPyRTLAoUjjshIoqNuChUonXrZ8OGDRHfJ0XH/ADZs6Np+Mt//G3UBsiu9fVOooXXmjzMTCaJucVFoRIt3d3dqKio0LoZpCJWA2R5eye6eK3Jw8xkkpgbCxUVs7OzWjeBlhHt2zscGBtbvNbkYWYyScwtLgqVaI1RSU5Ojti+KDKiXaCw50QbvNbkYWYyScwtLgqVaI1RKSoqivg+KTzRLFA4MFZ7vNbkYWYyScwtLgqVaLly5Yq4pYbj0dn2Mfzvs10RH4PCgbH6wWtNHmYmk8TcWKiQbs33ovzkw76VNw4BCxQiIjlYqKjIzc3Vuglr0nyB8rPLA5hyR+42DwsU/eK1Jg8zk0libnFRqERrMK3BwA+zWIrWOBQOkNU/XmvyMDOZJOYWF4VKtAbT9vf3Y926dRHfLwWKRoHCAbKy8FqTh5nJJDG3uChUSKZoFCiWBGDv7Xb2nhARxQkWKirWr1+vdRPiVqRn8syPP3mkJgspyUkR2SfFDq81eZiZTBJzS9C6AXrmcDi0bkLc8foU/PRDB47Vt0ekSEk2GfDYFjt+9thm/NGWfAz0R3aGEMUGrzV5mJlMEnNjj4qK6elprZsQNyI9k2e5AbLMTCbmJg8zk0libnFRqERr1o/FYonYvtaqaIxDeWzL8mNQmJlMzE0eZiaTxNziolCJ1qyf0tLSiO9zLYn0OJR0ixFHdpRgR3nmstswM5mYmzzMTCaJuXGMioqWlhatmyDWmbbRiI9DOfGHt6oWKQAzk4q5ycPMZJKYW1z0qJB+eH0K/umCA/90oX/V++JKskRExEJFRU5OjtZNECVSt3pWU6AwM5mYmzzMTCaJubFQUWEy8fAE60zbKL71y45V7SMSS90zM5mYmzzMTCaJuclrcQw5HA5kZWVp3Qzd8voUXO6bwvudY3jj46FV7UttJk8omJlMzE0eZiaTxNziolCJ1vRkWt7Z9jEcP9eDIefcqvYTzEweIiJau+KiUInW9OSKioqI7zMeROI2T7QGyjIzmZibPMxMJom5cXqyioGBAa2boCten4KfnL+Gv15FkWIA8Ed35PmXvI/0bB5mJhNzk4eZySQxt7joUYmWqakprZugG5Ga0fP0Z8rw6Yro3R9lZjIxN3mYmUwSc2OhosJsNmvdBF2IxK2eWI1FYWYyMTd5mJlMEnNjoaJC4r28SDvTNopvv90R9s8bAPzhHXn4wzsif5tnKcxMJuYmDzOTSWJuHKOioqmpSesmaGq+J8WnhL+Ppz9ThsfqCmK2suxaz0wq5iYPM5NJYm7sUaEAkVobJSc1EU9sLeK0YyIiWhUWKirWrVundRNiarVroxgAPLQpG3eXZaLWnqbJ83nWWmbxgrnJw8xkkpgbCxUVSUlJWjchZiIxYDbaM3qCsZYyiyfMTR5mJpPE3DhGRUVvb6/WTYiJ1Q6YTbcY8cyucs2LFGDtZBZvmJs8zEwmibnFRY8Kl9AP32p6UmI9o4eIiNaeuChUorWEfllZWcT3qSer7UnRw62em8V7ZvGKucnDzGSSmBtv/agYGRnRuglRs5qpxwkG4H/psEgB4juzeMbc5GFmMknMjYWKiomJCa2bEBWr7Un5y/v0WaQA8ZtZvGNu8jAzmSTmFhe3fqLFZIq/w7OaMSkS1kaJx8zWAuYmDzOTSWJu8locQ1VVVVo3IaLC6UnRw9oooYi3zNYK5iYPM5NJYm689aOioaFB6yZEhNen4KcfOsIak/L0Z8pwaHsxNhdYdV+kAPGT2VrD3ORhZjJJzI09KnHubPsY/v5X3Rie9oT0cwkGfY9FISKitYGFioqsLNkf0qsZjyK1SJGe2VrF3ORhZjJJzI2Fioq0tDStmxC2cGf2SO9JkZzZWsbc5GFmMknMjWNUVHR3d2vdhLCcbR8Le40UyUUKIDeztY65ycPMZJKYG3tU4ozXp+D4r3pC/jnpPSlERBSfdNej0t3djZ07d2LTpk247bbb8Prrr2vWlpKSEs1+d7j+5WIfhqbnQv65eClSJGZGzE0iZiaTxNx0V6iYTCY8//zzaGhowC9+8Qt87Wtfg9Pp1KQt4+PjmvzecJ1pG8VPPuwL6WdyUhN18+TjSJCWGV3H3ORhZjJJzE13hUp+fj5uv/12AEBubi5sNptmzyaQFGg4g2cf22LHT/beouuVZkMlKTO6gbnJw8xkkphbyIXKmTNn8MADD6CgoAAGgwFvvPHGom2OHz+O8vJyJCUloa6uDu+++25Yjfvggw/g8/lQXFwc1s+vVkKC7uq4JYX6gMH5hwr+0ZZ8EYu4hUJKZhSIucnDzGSSmFvIg2mdTic2b96M/fv34wtf+MKi7584cQJHjhzB8ePHcffdd+MHP/gBdu/ejYaGBv+9sbq6OrhcrkU/+/Of/xwFBQUAgOHhYTz22GN4+eWXV2zT5ORkwMG3WCywWCyh/tMW2bhx46r3EW3h9KTEy3iUpUjIjBZjbvIwM5kk5mZQFCWMSayf/LDBgJMnT2LPnj3+1+666y5s2bIFL774ov+1mpoa7NmzB88++2xQ+3W5XPjsZz+LL3/5y3j00UeX3c7j8eD06dN48MEHMT097X99//79OHz4MPLz83H16lUAQF5eHhRFwcDAAABgw4YN6OnpwczMDJKSklBcXIwrV64AuH7LKSEhAa2trUhLS8P69evR19cHp9MJi8WCsrIyNDc3AwCys7NhNptx7do1AEB5eTkGBwcxNTWFxMREVFZWorGxEQBgs9mQnJyM3t5eAEBZWRlGRkYwMTEBo9GI6upqNDY2QlEUZGZmwmq1+qeSlZSUYGJiAmNjYzAYDKipqcGJ9z7GPzS6gzqm8x5an4w/3lqK9vZ2AEBBQQFcLheGh4cBXD+J29ra4Ha7kZqairy8PLS1tQG4flvO4/FgcHAQwPVnRnR1dWF2dhbJyckoLCxEa2ur/3gDQH9/PwCgsrISvb29/uNdUlKClpYWAEBOTg5MJhMcDgcAoKKiAv39/XA6nTCbzaioqEBTUxMAYN26dbBYLAHHe2hoCJOTkzCZTPD5fPD5fP7jnZKSgp6e67OgSktLMTo6uuzxTk9PR1dXFwCguLgYk5OTAce7ubkZXq8X6enpsNls6OjoAAAUFhZiZmbGf4uypqYGra2tmJubQ1paGnJycgKOt9vtxtDQEACguroaHR0dcLlcSE1Nhd1u95+zdrsdPp8v4Jzt7u72H++ioqKAc9ZgMPiP9/r16+FwODA9PQ2LxYLS0lLV4z0wMICpqaklj3dSUtKS56zJZEJVVZV/Se6srCykpaUFnLPj4+MYHx9HQkICNm7ciKamJvh8PmRkZCAjI8N/vL1eL7KzszE6OgoA2LRpE1paWuDxeJY83rOzs0ues2lpacjNzVU9Zzs7O+FyuZCSkrLq94i+vj7/8dbje8TCczYrKwudnZ0AgKKiIkxPT/vP2YXH22q1Ijs7e8X3iJGREeTl5Yl7j1h4zq7F94irV6/6P9e0fI9QFAUulwv33nvvig9KjGih4na7kZKSgtdffx0PP/ywf7uvfvWruHjxIk6fPr3iPhVFwR/8wR+guroa3/zmN1W3nS9UKioqotKj0tDQgE2bNq16P9Hg9Sl49LWPQ5rhk51iwj/+fm3c3e5ZSM+Z0fKYmzzMTCa95Ob1enHp0qWgCpWIrqMyNDQEr9frr5Tn5eXl+f/PYyXvvfceTpw4gdtuu80//uUf//Efceutty77M1arFUajMex2LycjIyPi+4yUcKYhH9pWHNdFCqDvzGh5zE0eZiaTxNyisuCbwRD4YagoyqLXlrNjxw5/132wdu3ahdnZWRw4cAAHDhwI6WfV6DXQs+1jIU1Dnl/MLZ5m9yxHr5mROuYmDzOTSWJuER3+m52dDaPRuKj3ZGBgYFEvSyTV19fj3LlzES1SAPjvReqJ2+PD//teaO2K58GzN9NjZrQy5iYPM5NJYm4RLVTMZjPq6upw6tSpgNdPnTqF7du3R/JXrUln28fwB69dxvisN6jt56chr5UihYiI4k/It36mpqb8o7YBoL29HRcvXoTNZkNJSQmOHj2KRx99FHfeeSe2bduGl156CV1dXTh48GBEG75QtG79aLV+y1LOto/hWH17SD+zlnpS5ukpMwoec5OHmckkMbeQC5UPPvgA9913n//ro0ePAgD27duHV155BXv37sXw8DCOHTsGh8OB2tpavPnmmygtLY1cq29SX18flcG0U1NTsFqtEd9vqMJ50OBjW+xrrkgB9JMZhYa5ycPMZJKYW8iFys6dO7HSjOZDhw7h0KFDYTdKL0ZHR5Gfn691M0Ke4ZOdYsIXb7dHsUX6pZfMKDTMTR5mJpPE3OStpbvGhDrDB1gb05CJiGhtiItCZdeuXdi6dWtQy+2HQutFcUK95ZOZZMIzu8rXxDTk5WidGYWHucnDzGSSmFtcFCrRmp48v5ywVkK55ZORZMRPfz++noQcDq0zo/AwN3mYmUwSc4uLQiVaPB6PZr871Fs+X727BGYT49QyMwofc5OHmckkMTd+sqlIT0/X5PeGesvnsS32Nd+TMk+rzGh1mJs8zEwmibnFRaESrTEqNpstovsLVii3fNbyDJ+laJUZrQ5zk4eZySQxt7goVKI1RmX+Ed2xFOotH87wCaRFZrR6zE0eZiaTxNziolCJF16fguPneMuHiIhoHgsVFYWFhTH9fZf7pjDk5C2f1Yh1ZhQZzE0eZiaTxNziolCJ1hiV2dnZiO5vJb/qHA96W97yWVqsM6PIYG7yMDOZJOYW8hL6ehStZ/0MDw8jLy8v4vtdytn2Mfzs48GgtuUtn+XFMjOKHOYmDzOTSWJucdGjIl0oY1N4y4eIiNYSFioqNm7cGJPf81vHZNBjU3jLR12sMqPIYm7yMDOZJObGQkVFW1tb1H/H2fYxfOuXHUFt+0htDm/5rCAWmVHkMTd5mJlMEnOLizEq0eJ2u6O6/7PtYzhW3x709ttKMqLYmvgQ7cwoOpibPMxMJom5xUWPSrRm/aSlpUV0fwuFumZKTmoiau3Ra0+8iGZmFD3MTR5mJpPE3OKiRyVas35yc3Mjvs95oayZAgBPbC3i2JQgRDMzih7mJg8zk0libnHRoxIt0byXF+yaKVaLEc/sKufYlCBJvP9KzE0iZiaTxNziokdFmlDWTPlfnynDHYXynnZJREQUCexRUZGfnx/xfYYyNiUnNRG35Vsj3oZ4Fo3MKPqYmzzMTCaJubFQUeHxeCK+z1DGpnBcSuiikRlFH3OTh5nJJDE3FioqBgeDuz0TimDHpnDNlPBEIzOKPuYmDzOTSWJuLFRiKJSxKVwzhYiIKE4KlWito1JVVRWxfYU6NoVrpoQnkplR7DA3eZiZTBJzi4tCpb6+HufOncOBAwciut/Ozs6I7YtjU2IjkplR7DA3eZiZTBJzi4tCJVpcLlfE9jUyE1yRwrEpqxPJzCh2mJs8zEwmibmxUFGRkpISsX31jgd3cnBsyupEMjOKHeYmDzOTSWJuLFRURGq++dn2Mfzkw74Vt+PYlNWTuEYAMTeJmJlMEnNjoaLi6tWrq95HKINoOTZl9SKRGcUec5OHmckkMTcWKlEW7CDax7bYOTaFiIjoJixUVOTl5a16H8EOoi3MsKz6d1FkMqPYY27yMDOZJObGQkWFoiir3octOTGi25G6SGRGscfc5GFmMknMjYWKioGBgVXvY3zWg5WGnXAQbeREIjOKPeYmDzOTSWJucVGoRGtl2tU62z6Gb/2yA74VClgOoiUiIlqaSesGREJ9fT2MRmPE97thw4awfzaY2T4JBuAv7yvjINoIWk1mpB3mJg8zk0libnHRoxItPT3BTSteSjCzfXwKkJEUF7WibqwmM9IOc5OHmckkMTcWKipmZmbC/tlgZ/sEux0FZzWZkXaYmzzMTCaJubFQUZGUlBT2z2YG2VPC2T6RtZrMSDvMTR5mJpPE3HjfQUVxcXFYP3e2fQx//6vuFbfjbJ/ICzcz0hZzk4eZySQxN/aoqLhy5UrIP3O2fQzH6tsxPO1ZcVvO9om8cDIj7TE3eZiZTBJzY49KBAX7XJ+c1EQ8sbWIs32IiIhWwEJFRW5ubkjbB/tcn69/ugR3FKaH2yxSEWpmpA/MTR5mJpPE3HjrR0VCQmiHJ9gZPGOzK98WovCEmhnpA3OTh5nJJDE3eS2Oob6+vpC253N9tBdqZqQPzE0eZiaTxNxYqERQrT0N2SnqRQhn+hAREQVPd4XK5OQkPvWpT+H222/Hrbfeih/+8IeatWX9+vUhbf+rznG4vD7VbTjTJ7pCzYz0gbnJw8xkkpib7gqVlJQUnD59GhcvXsSvf/1rPPvssxgeHtakLaF0kc1PS550eZf8frrFiGd2lXOmT5RJ7NYk5iYRM5NJYm66K1SMRiNSUlIAALOzs/B6vVCUFR4/HCVOpzOo7YKZlmw2GrCtNCMSzSIVwWZG+sLc5GFmMknMLeRC5cyZM3jggQdQUFAAg8GAN954Y9E2x48fR3l5OZKSklBXV4d33303pN8xNjaGzZs3o6ioCN/4xjeQnZ0dajMjwmKxBLVdMNOSh6Y9uNw3FYlmkYpgMyN9YW7yMDOZJOYWcqHidDqxefNmvPDCC0t+/8SJEzhy5AiefvppXLhwAffccw92796Nrq4u/zZ1dXWora1d9OfatWsAgMzMTFy6dAnt7e3453/+Z/T394f5z1udsrKyoLbjAwj1I9jMSF+YmzzMTCaJuYVcqOzevRvf+ta38Mgjjyz5/eeeew6PP/44Dhw4gJqaGjz//PMoLi7Giy++6N/m/PnzuHz58qI/BQUFAfvKy8vDbbfdhjNnzqi2aXJyEhMTE/4/Lpcr1H/Wkpqbm4PajtOS9SPYzEhfmJs8zEwmiblFdGVat9uN8+fP46mnngp4/f7778f7778f1D76+/uRnJyM9PR0TExM4MyZM3jiiSdUf6a2thbT09P+r/fv34/Dhw8jPz8fV69eBXC96FEUBQMDAwCADRs2oKenBzMzM0hKSkJxcbH/GQi5ublISEjA5OQkGhoasH79evT19cHpdMJisaCsrMwfdnZ2NoqTE5FpBsbcy7cx0wzkGKYwPu5Db28vgOuV7cjICCYmJmA0GlFdXY3GxkYoioLMzExYrVZ0d19/uGFJSQkmJiYwNjYGg8GAmpoaNDc3w+v1Ij09HVlZWejs7AQAFBUVYXp6GiMjIwCATZs2oaWlBR6PB1arFdnZ2WhvbwcAFBQUwOVy+Qcsb9y4EW1tbXC73UhNTUVeXh7a2toAAPn5+fB4PBgcHAQAVFVVoaurC7Ozs0hOTkZhYSFaW1v9x3s+TwCorKxEb2+v/3iXlJSgpaUFAJCTkwOTyQSHwwEAqKioQH9/P5xOJ8xmMyoqKtDU1AQAWLduHSwWi7/3rby8HENDQ5icnITJdP10bmhoAADYbDakpKSgp+f6+KHS0lKMjo4ue7zT09P9PX/FxcWYnJxc9njbbDZ0dHQAAAoLCzEzM+M/3jU1NWhtbcXc3BzS0tKQk5MTcLzdbjeGhoYAANXV1ejo6IDL5UJqairsdrv/nLXb7fD5fAHnbHd3t/94FxUVBZyzBoPBf7zXr18Ph8OB6elpWCwWlJaWqh7vgYEBTE1NLXm8k5KSljxnTSYTqqqq/Mc7KysLaWlpAefs+Pg4xsfHkZCQgI0bN6KpqQk+nw8ZGRnIyMjwH2+PxwOHw4HR0dFF5+xSx3t2dnbJczYtLQ25ubmq52xnZydcLhdSUlJW/R4xPzBxpfcIs9kccM4ODg5iamoKiYmJqKysRGNjo/+cTU5OFvEeMTk5ic7OTnHvEQvP2bX4HrHwc03L94hQxp4alFWMVDUYDDh58iT27NkDALh27RoKCwvx3nvvYfv27f7tvv3tb+PVV18NqpI7f/48Hn/8cSiKAkVR8MQTTyxbqHg8Hpw+fRoVFRUBq+1ZLJaI3IcbGBgIernhM22j+NYvO5b9Pmf8xEYomZF+MDd5mJlMesnN6/Xi0qVLuPfee/3/k7mcqDzrx2AIXCdEUZRFry2nrq4OFy9eDOn3Wa1WGI3GkH4mGGazOajtzraP4fu/7l3ye3wAYWwFmxnpC3OTh5nJJDG3iE5Pzs7OhtFoXDRPe2BgwN/NFw27du3C1q1b8fLLL0d0v/Ndh2rm109ZbtbPn9xVyCIlhoLJjPSHucnDzGSSmFtEe1TMZjPq6upw6tQpPPzww/7XT506hYceeiiSvypAfX19VHpUVhLM+ik/+HUv7i7L5Gq0REREYQi5UJmamvIPhgKA9vZ2XLx4ETabDSUlJTh69CgeffRR3Hnnndi2bRteeukldHV14eDBgxFteCyUl5erfj+Y9VMGnXO43DeFzQXWSDaNlrFSZqRPzE0eZiaTxNxCvvXzwQcf4I477sAdd9wBADh69CjuuOMOPPPMMwCAvXv34vnnn8exY8dw++2348yZM3jzzTdRWloa2ZYvEK1bP/Mj15fD9VP0Z6XMSJ+YmzzMTCaJuYXco7Jz584VpxUdOnQIhw4dCrtRoYrWrZ+pKfWVZLl+iv6slBnpE3OTh5nJJDE33T3rR08SE9ULjFp7GrJT1bfJSU1ErT0tks0iFStlRvrE3ORhZjJJzC0uCpVo3fqprKxU/b4xwYBDW4tUt3liaxEH0sbQSpmRPjE3eZiZTBJzi4tCpb6+HufOncOBAwciut/51SLVbCvNQG1e6qLXc1ITucibBoLJjPSHucnDzGSSmFtUFnxbK862j+H4uZ6AmT+piQn4wq25+OLtdvakEBERrVJc9KhEi81mW/Z7yy305pzz4Scf9uFXnePRbh4tQS0z0i/mJg8zk0libnFRqERrjEpycvKSrwez0NuL53rg9YX9GCUK03KZkb4xN3mYmUwSc4uLWz/Rmp7c29uLjIyMRa9zoTf9Wi4z0jfmJg8zk0libnHRoxJrXOiNiIgoNlioqCgrK1vydS70pl/LZUb6xtzkYWYyScyNhYqKkZGRJV/nQm/6tVxmpG/MTR5mJpPE3OKiUInWYNqJiYklX+dCb/q1XGakb8xNHmYmk8TcOJhWhdo+d5RnYnf1OvxX83DA6zmpiXhiaxEXetNINM4Dij7mJg8zk0libnFRqERLdXV10Ns+eocdt+Wnodaexp4UDYWSGekHc5OHmckkMbe4uPUTLWpLDXt9Cj5y3HgK5cO1OdhcYGWRojGJy0MTc5OImckkMTcWKioUZekF2862j+HR1z5Gz4TL/9r//bMmnG0fi1HLaDnLZUb6xtzkYWYyScyNhYqKzMzMRa/5l86fDlwjZcg5h2P17SxWNLZUZqR/zE0eZiaTxNziolCJ1qwfqzVwVVkuna9/N2dGMjA3eZiZTBJzi4tCpb6+HufOncOBAwciut/u7u6Ar0NZOp+0cXNmJANzk4eZySQxt7goVGKFS+cTERHFFgsVFSUlJQFfc+l8/bs5M5KBucnDzGSSmBsLFRU3r+DHpfP1T+Kqi8TcJGJmMknMjYWKirGxsYCvuXS+/t2cGcnA3ORhZjJJzI2FigqDYXHBsaM8E//HrTmLXs9JTcQzu8q5dL7GlsqM9I+5ycPMZJKYW1wsob9r1y7Mzs7iwIEDEZ35U1NTs+TrVsuNw7ZnUzbuLsvk0vk6sVxmpG/MTR5mJpPE3OKiRyVa05Obm5uXfL1n/MaKtL9bvY5L5+vIcpmRvjE3eZiZTBJzi4tCJVq8Xu+Sr/cuKFQK0i2xag4FYbnMSN+YmzzMTCaJubFQUZGenr7k6z3jswCA7NREJCfKe2R2PFsuM9I35iYPM5NJYm4sVFRkZWUtem1i1oMJ1/WKtCiDvSl6s1RmpH/MTR5mJpPE3FioqOjs7Az42utT8E7bqP/rQisLFb25OTOSgbnJw8xkkphbXMz6iYWz7WM4fq4n4Fk/77SPoq4onVOSiYiIooQ9KiqKiq4v7na2fQzH6tsXPZDQ6fbhWH07zraPadA6Wsp8ZiQLc5OHmckkMTcWKiqmp6fh9Sk4fq5HdbsXz/XA61Ni1CpSMz09rXUTKAzMTR5mJpPE3FioqBgZGcHlvqlFPSk3G3TO4XLfVIxaRWpGRka0bgKFgbnJw8xkkphbXBQqu3btwtatW/Hyyy9HfN8jM+pFSqjbERERUfDiYjBtfX09jMbIr2eyadMmzF2bDGpbW7L6U5UpNjZt2qR1EygMzE0eZiaTxNziokclWlpaWlBrT0N2qnoRkpOaiFp7WoxaRWpaWlq0bgKFgbnJw8xkkpgbCxUVHo8HxgQDDm1VHyX9xNYiPutHJzwej9ZNoDAwN3mYmUwSc2OhosJqtQIAdpRn4pld5bCaA28v5aQm4pld5VxHRUfmMyNZmJs8zEwmibnFxRiVaMnOzvb/fUd5JjrHZvDq+T4AwKNb7PiD2+3sSdGZhZmRHMxNHmYmk8Tc2KOior29PeDr0ZkbXWZ3FqWzSNGhmzMjGZibPMxMJom5sVAJwcL1VNalcJYPERFRtLFQUVFQUBDw9fD09UIlwQDYWKjo0s2ZkQzMTR5mJpPE3FioqHC5XAFfDzrdAIDMZBNMvO2jSzdnRjIwN3mYmUwSc2OhomJ4eNj/d69PwdgnY1SyU8xaNYlWsDAzkoO5ycPMZJKYGwuVII3MzGH+uYPrVlgAjoiIiCKDhYqKjRs3+v++cCBtNsen6NbCzEgO5iYPM5NJYm66LVSmp6dRWlqKr3/965q1oa2tzf/3gEKFPSq6tTAzkoO5ycPMZJKYm24Llb/+67/GXXfdpWkb3G63/+9Dzht/Z6GiXwszIzmYmzzMTCaJuemyULly5Qqamprwuc99TtN2pKam+v8+PzUZ4GBaPVuYGcnB3ORhZjJJzC3kQuXMmTN44IEHUFBQAIPBgDfeeGPRNsePH0d5eTmSkpJQV1eHd999N6Tf8fWvfx3PPvtsqE2LuLy8PADXZ/xcGZr2v56VwicP6NV8ZiQLc5OHmckkMbeQCxWn04nNmzfjhRdeWPL7J06cwJEjR/D000/jwoULuOeee7B79250dXX5t6mrq0Ntbe2iP9euXcO//du/oaqqClVVVUG3aXJyEhMTE/4/kZon3tbWhrPtY3j0xMe4cG3K//pf/NdVnG0fi8jvoMiSeP+VmJtEzEwmibkZFEVRwv5hgwEnT57Enj17/K/ddddd2LJlC1588UX/azU1NdizZ09QvSR/8Rd/gZ/+9KcwGo2YmprC3NwcnnzySTzzzDOLtvV4PDh9+jQefPBBTE/f6PHYv38/Dh8+jPz8fFy9ehXA9SpSURQMDAwAADZs2ICenh7MzMwgKSkJxcXFuHLlCgAgNzcXCQkJeOODVrzWtfx4lP3VJuyqzoPZbMa1a9cAAOXl5RgcHMTU1BQSExNRWVmJxsZGAIDNZkNycjJ6e3sBAGVlZRgZGcHExASMRiOqq6vR2NgIRVGQmZkJq9WK7u5uAEBJSQkmJiYwNjYGg8GAmpoaNDc3w+v1Ij09HVlZWejs7AQAFBUVYXp6GiMjIwCATZs2oaWlBR6PB1arFdnZ2f7nPRQUFMDlcvnn1m/cuBFtbW1wu91ITU1FXl6e/8TOz8+Hx+PB4OAgAKCqqgpdXV2YnZ1FcnIyCgsL0dra6j/eANDf3w8AqKysRG9vr/94l5SUoKWlBQCQk5MDk8kEh8MBAKioqEB/fz+cTifMZjMqKirQ1NQEAFi3bh0sFkvA8R4aGsLk5CRMJlPAI8xtNhtSUlLQ09MDACgtLcXo6Oiyxzs9Pd1fUBcXF2NycnLZ422z2dDR0QEAKCwsxMzMjP9419TUoLW1FXNzc0hLS0NOTk7A8Xa73RgaGgIAVFdXo6OjAy6XC6mpqbDb7f5z1m63w+fzBZyz3d3d/uNdVFQUcM4aDAb/8V6/fj0cDgemp6dhsVhQWlqqerwHBgYwNTW15PFOSkpa8pw1mUyoqqpCQ0MDACArKwtpaWkB5+z4+DjGx8eRkJCAjRs3oqmpCT6fDxkZGcjIyPAfb4/Hg5ycHIyOji46Z5c63rOzs0ues2lpacjNzVU9Zzs7O+FyuZCSkrLq94i+vj7/8e7r64PT6YTFYkFZWRmam5sBXH8IXDy+RwwPD8Nut4t7j1h4zq7F94jW1lZYrVbN3yMURYHL5cK9994Lk0n9LkVECxW3242UlBS8/vrrePjhh/3bffWrX8XFixdx+vTpkPb/yiuv4PLly/jbv/3bJb8/X6hUVFQgIeFG55DFYoHFYgn9H7SA16fgD//lI4zMeJfdJic1ET/ZewsfTqgjo6OjyMrK0roZFCLmJg8zk0kvuXm9Xly6dCmoQiWigy2Ghobg9XoX3QPLy8vz/59HNFitVhiNxoju83LflGqRAgCDzjlc7pvC5gJrRH83hW9hjwrJwdzkYWYyScwtKrN+DIbAHgZFURa9FowvfelLy/amLLRr1y5s3boVL7/8csi/YzkjM3MrbxTCdhQb813OJAtzk4eZySQxt4j2qGRnZ8NoNC7qPRkYGIjqSOP6+vqI96jYkoNbKyXY7YiIiCh0Ee1RMZvNqKurw6lTpwJeP3XqFLZv3x7JXxV1tfa0FZfKz0lNRK09LUYtomCEMluM9IO5ycPMZJKYW8iFytTUFC5evIiLFy8CANrb23Hx4kX/SOijR4/i5Zdfxo9+9CM0Njbia1/7Grq6unDw4MGINnyhaNz6MSYY8IUK9ULlia1FHEirMwunwZMczE0eZiaTxNxCvvXzwQcf4L777vN/ffToUQDAvn378Morr2Dv3r0YHh7GsWPH4HA4UFtbizfffBOlpaWRa/VNonHrBwBq0r14Zlc5nnu3C1PuGwNrc1IT8cTWIuwoz4z476TVmZ2d1boJFAbmJg8zk0libiEXKjt37sRKM5oPHTqEQ4cOhd0ovUhOTsaO8kz0TMziR7+5Pp/8j7bk4Q9vz2dPik4lJydr3QQKA3OTh5nJJDE3XT7rRy8KCwsBANNun/+1W/OsLFJ0bD4zkoW5ycPMZJKYW1wUKtEYowLAv4Liwts+qZbI32KiyJnPjGRhbvIwM5kk5hYXT9eL1hiVeVOuGwvkWM0sVIiIiGIlLnpUomV+7ZeAHhUWKrom8cmgxNwkYmYyScyNhUoQplwsVIiIiLQQF4VKtMaozD/Vc75HJSUxgQNpdW4+M5KFucnDzGSSmBvHqARhvkfFaomLw0VERCRGXPSoREtlZSUURYHzkx4V3vbRv8rKSq2bQGFgbvIwM5kk5sZCRUVvby9cXgVzvusL3Fk5NVn3ent7tW4ChYG5ycPMZJKYW1wUKtEaozIzMwMnB9KKMjMzo3UTKAzMTR5mJpPE3OJi0EW0xqgkJSVh0r1gDRX2qOheUlKS1k2gMDA3eZiZTBJzi4selWgpKSlhj4owJSUlWjeBwsDc5GFmMknMjYWKipaWFkwuWOwtjbN+dK+lpUXrJlAYmJs8zEwmibmxUFnBwsXe0tijQkREFFMsVFTk5OQELJ/PQkX/cnJytG4ChYG5ycPMZJKYW1wUKtGa9WMymQILFQ6m1T2TibfnJGJu8jAzmSTmJq/FS4jWrB+Hw4EpV7r/az45Wf8cDgeysrK0bgaFiLnJw8xkkphbXPSoRJNz4ZOT2aNCREQUUyxUVFRUVGBywWBaqzkuOqDiWkVFhdZNoDAwN3mYmUwSc2OhoqK/v589KsJIfDIoMTeJmJlMEnNjF4EKp9OJSdf1Wi4xwQCL0aBxi2glTqdT6yZQGJibPMxMJom5sUdFhdlsDnhyssHAQkXvzGaz1k2gMDA3eZiZTBJzi4tCJVrTk6+PUbn+rB9OTZZB4v1XYm4SMTOZJOYWF4VKfX09zp07hwMHDkR0v5cbGjA95wMAGAB4fUpE90+R19TUpHUTKAzMTR5mJpPE3OKiUImGs+1jOHZ+zv9197gLj574GGfbx7RrFBER0RrDQmUJZ9vHcKy+HRNzga8POedwrL6dxYqOrVu3TusmUBiYmzzMTCaJubFQuYnXp+D4uR7VbV4818PbQDplsVi0bgKFgbnJw8xkkpgbC5WbXO6bwpBzTnWbQeccLvdNxahFFIpr165p3QQKA3OTh5nJJDE3Fio3GZlRL1JC3Y6IiIjCx0LlJrbkxIhuR7FVXl6udRMoDMxNHmYmk8TcWKjcpNaehuxU9SIkJzURtfa0GLWIQjE0NKR1EygMzE0eZiaTxNxYqNzEmGDAoa1Fqts8sbUIxgSuUqtHk5OTWjeBwsDc5GFmMknMjYXKEnaUZ+KZXeVIumkx2pzURDyzqxw7yjM1aRetzGTi46skYm7yMDOZJOYWF4VKNJbQ31Geif9rs93/9Zfq8vGTvbewSNG5qqoqrZtAYWBu8jAzmSTmFheFSrSW0HcMDPr/vjE3hbd7BGhoaNC6CRQG5iYPM5NJYm5xUahEy5z3xt8tRh4qIiKiWOOnr4oE840V/MwmHioJbDab1k2gMDA3eZiZTBJz46evCp/hxmha9qjIkJKSonUTKAzMTR5mJpPE3Pjpq2Js0un/u4U9KiL09Kg/p4n0ibnJw8xkkpgbP31VzPlu/N1s4kBaIiKiWGOhosJkTvL/nbd+ZCgtLdW6CRQG5iYPM5NJYm789FXhdLn9f+etHxlGR0e1bgKFgbnJw8xkkpgbP31VzH4yP9loANdQEWJiYkLrJlAYmJs8zEwmibmxUFHh+WSMCntT5DAajStvRLrD3ORhZjJJzI2fwCoU4/WnKJs5PkWM6upqrZtAYWBu8jAzmSTmpstPYJPJhNtvvx233357xJfFD8X0J2NU2KMiR2Njo9ZNoDAwN3mYmUwSc9PlYxQzMzNx8eJFrZvhn55sNnJ8ihSKomjdBAoDc5OHmckkMTd2FaiY4xgVcTIzM7VuAoWBucnDzGSSmFvIn8BnzpzBAw88gIKCAhgMBrzxxhuLtjl+/DjKy8uRlJSEuro6vPvuuyH9jomJCdTV1WHHjh04ffp0qE2MCK9PgfeTwpNjVORIT0/XugkUBuYmDzOTSWJuIX8CO51ObN68GS+88MKS3z9x4gSOHDmCp59+GhcuXMA999yD3bt3o6ury79NXV0damtrF/25du0aAKCjowPnz5/H97//fTz22GOaTKdyeW4sS8seFTkWnmckB3OTh5nJJDG3kMeo7N69G7t37172+8899xwef/xx/yDY559/Hm+99RZefPFFPPvsswCA8+fPq/6OgoICAEBtbS02bdqElpYW3HnnnctuPzk5iYSEG8WExWKBxWJZdvtguLwLCxWOUSEiItJCRAfTut1unD9/Hk899VTA6/fffz/ef//9oPYxOjqKlJQUWCwW9PT0oKGhARUVFao/U1tbi+npaf/X+/fvx+HDh5Gfn4+rV68CAPLy8qAoCgYGBgAAGzZsQE9PD2ZmZpCUlITi4mJcuXIFAJCbm4vh2RuFiskAdHZ2wul0wmKxoKysDM3NzQCA7OxsmM1mf29QeXk5BgcHMTU1hcTERFRWVvpHWdtsNiQnJ6O3txcAUFZWhpGREUxMTMBoNKK6uhqNjY1QFAWZmZmwWq3o7u4GAJSUlGBiYgJjY2MwGAyoqalBc3MzvF4v0tPTkZWVhc7OTgBAUVERpqenMTIyAgD+Ys/j8cBqtSI7Oxvt7e0ArheFLpcLw8PDAICNGzeira0NbrcbqampyMvLQ1tbGwAgPz8fHo8Hg4ODAICqqip0dXVhdnYWycnJKCwsRGtrq/94A0B/fz8AoLKyEr29vf7jXVJSgpaWFgBATk4OTCYTHA4HAKCiogL9/f1wOp0wm82oqKhAU1MTAGDdunWwWCwBx3toaAiTk5MwmUwoLi5GQ0OD/3inpKT4H8JVWlqK0dHRZY93enq6//82iouLMTk5uezxttls6OjoAAAUFhZiZmbGf7xramrQ2tqKubk5pKWlIScnJ+B4u91uDA0NAbg+VbCjowMulwupqamw2+3+c9Zut8Pn8wWcs93d3f7jXVRUFHDOGgwG//Fev349HA4HpqenYbFYUFpaqnq8BwYGMDU1teTxTkpKWvKcNZlMqKqq8h/vrKwspKWlBZyz4+PjGB8fR0JCAjZu3Iimpib4fD5kZGQgIyPDf7xtNhscDod/1cyF5+xSx3t2dnbJczYtLQ25ubmq52xnZydcLhdSUlJW9R6RkJCAvr4+//Hu6+tbU+8RHo8HnZ2d4t4jFp6za/E9wuPxoKGhQfP3iFAG9RqUVQwBNhgMOHnyJPbs2QMAuHbtGgoLC/Hee+9h+/bt/u2+/e1v49VXX/VfuGref/99/Mmf/AkSEhJgMBjwzW9+07//m3k8Hpw+fRoVFRUR71HpGpvFgX+9/ubx2Q02/D/3yns+wlp07do1f48cycHc5GFmMuklN6/Xi0uXLuHee++FyaTeZxKV6ckGQ+CtEkVRFr22nO3bt+Ojjz4K6fdZrdaIr7bnXjhGhYNpxRgbG9PFRUihYW7yMDOZJOYW0U/g7OxsGI1Gf3fovIGBAX83XzTs2rULW7duxcsvvxyxfS4co2LmGBUxgi2ISV+YmzzMTCaJuUW0R8VsNqOurg6nTp3Cww8/7H/91KlTeOihhyL5qwLU19dHoUflxh0x9qjIUVNTo3UTKAzMTR5mJpPE3EL+BJ6amsLFixf9K8e2t7fj4sWL/gFGR48excsvv4wf/ehHaGxsxNe+9jV0dXXh4MGDEW14tAX2qLBQkSKYcVCkP8xNHmYmk8TcQu5R+eCDD3Dffff5vz569CgAYN++fXjllVewd+9eDA8P49ixY3A4HKitrcWbb76J0tLoDUbdtWsXZmdnceDAgYg9GyhwjIq8rrK1yuv1at0ECgNzk4eZySQxt5ALlZ07d644rejQoUM4dOhQ2I0KVTRu/bBHRSaJqy4Sc5OImckkMTd+Ai/DtWCMShILFTFsNpvWTaAwMDd5mJlMEnOLi0/gqMz6WXDrh8/6kWN+gSWShbnJw8xkkphbVNZRibWozPrx8lk/REREWuMn8DICe1Q4mFaKwsJCrZtAYWBu8jAzmSTmxkJlGW7vgnVU2KMixszMjNZNoDAwN3mYmUwSc4uLT+Coj1FhoSLG/EO/SBbmJg8zk0libhyjsoyAMSq89UNERKQJdhUsgz0qMklcHpqYm0TMTCaJufETeBkuL5/1I1Fra6vWTaAwMDd5mJlMEnPjJ/AyApbQZ4+KGHNzc1o3gcLA3ORhZjJJzC0uPoGjMpiW66iIlJaWpnUTKAzMTR5mJpPE3DiYdhnzY1QSDIApgYNppcjJydG6CRQG5iYPM5NJYm7sKliG+5Nn/bA3RZb29natm0BhYG7yMDOZJObGT+FlzN/64XN+iIiItMNP4WXMD6a1mHjbR5KCggKtm0BhYG7yMDOZJOYWF4VKdAbTXr/1wx4VWdxut9ZNoDAwN3mYmUwSc+Ng2mXc6FFhoSLJ0NAQcnNztW4GhYi5ycPMZJKYGz+Fl+D1KZjzfTKYlj0qREREmuGn8BIWPufHzDEqolRXV2vdBAoDc5OHmckkMTcWKktwc/l8sTo6OrRuAoWBucnDzGSSmBs/hZfABxLK5XK5tG4ChYG5ycPMZJKYGz+Fl+Dm8vlipaamat0ECgNzk4eZySQxN34KL2Fhj4rFyDEqktjtdq2bQGFgbvIwM5kk5hYXhUqk11FxeW6MUeGtH1muXr2qdRMoDMxNHmYmk8TcuI7KEgKenMzBtERERJrhp/ASZue8/r+PzMzB61NUtiY9kditScxNImYmk8TcWKjc5Gz7GP72TJf/65+3jODREx/jbPuYdo2ioPl8vpU3It1hbvIwM5kk5sZCZYGz7WM4Vt+OSZc34PUh5xyO1bezWBFgYGBA6yZQGJibPMxMJom5sVD5hNen4Pi5HtVtXjzXw9tAREREMcRC5ROX+6Yw5JxT3WbQOYfLfVMxahGFY8OGDVo3gcLA3ORhZjJJzI2FyidGZtSLlFC3I210d3dr3QQKA3OTh5nJJDE3FiqfsCUnRnQ70sbs7KzWTaAwMDd5mJlMEnNjofKJWnsaslPVi5Cc1ETU2tNi1CIKR3JystZNoDAwN3mYmUwSc4uLQiUSK9MaEww4tLVIdZsnthbBmMAl9fWsqEg9Q9In5iYPM5NJYm5xUajU19fj3LlzOHDgwKr2s6M8E8/sKl/Us5KTmohndpVjR3nmqvZP0XflyhWtm0BhYG7yMDOZJOYWF0voR9KO8kxsK83A5b4p/PZKB27bUIZaexp7UoiIiDTAQmUJxgQDNhdYUWjOR3a2VevmUAhyc3O1bgKFgbnJw8xkkphbXNz6iRaDgb0o0jAzmZibPMxMJom5sVBR0d/fr3UTKETMTCbmJg8zk0libixUiIiISLdYqKhYv3691k2gEDEzmZibPMxMJom5sVBR4XA4tG4ChYiZycTc5GFmMknMjYXKMlwuF/7u7/4OLpdL66ZQkJiZTMxNHmYmk9TcWKgsw+Vy4cc//rG4QNcyZiYTc5OHmckkNTcWKkRERKRbLFSIiIhIt0SvTKsoCgDA6/VGfN8+nw8pKSnw+XxR2T9FHjOTibnJw8xk0lNu879//nNcjUEJZiudmp2dxXvvvad1M4iIiCgMd999N5KSklS3EV2o+Hw+uN1uGI1GkcsCExERrUWKosDr9cJsNiMhQX0UiuhChYiIiOIbB9MSERGRbrFQISIiIt1ioUJERES6xUJlCcePH0d5eTmSkpJQV1eHd999V+sm0QLf/OY3YTAYAv7Y7Xb/9xVFwTe/+U0UFBQgOTkZO3fuxMcff6xhi9eeM2fO4IEHHkBBQQEMBgPeeOONgO8Hk5HL5cLhw4eRnZ2N1NRUPPjgg+jp6Ynhv2LtWSm3L33pS4uuva1btwZsw9xi59lnn8WnPvUpWK1W5ObmYs+ePWhubg7YJh6uNRYqNzlx4gSOHDmCp59+GhcuXMA999yD3bt3o6urS+um0QK33HILHA6H/89HH33k/97f/M3f4LnnnsMLL7yA3/zmN7Db7fjsZz+LyclJDVu8tjidTmzevBkvvPDCkt8PJqMjR47g5MmTeO2113D27FlMTU3h85//vObrP8SzlXIDgN/7vd8LuPbefPPNgO8zt9g5ffo0/vRP/xTnzp3DqVOn4PF4cP/998PpdPq3iYtrTaEAv/M7v6McPHgw4LWNGzcqTz31lEYtopv91V/9lbJ58+Ylv+fz+RS73a585zvf8b82OzurZGRkKN///vdj1EJaCIBy8uRJ/9fBZDQ2NqYkJiYqr732mn+b3t5eJSEhQfnv//7vmLV9Lbs5N0VRlH379ikPPfTQsj/D3LQ1MDCgAFBOnz6tKEr8XGvsUVnA7Xbj/PnzuP/++wNev//++/H+++9r1CpaypUrV1BQUIDy8nL8/u//Ptra2gAA7e3t6OvrC8jQYrHg3nvvZYY6EUxG58+fx9zcXMA2BQUFqK2tZY4ae+edd5Cbm4uqqip8+ctfxsDAgP97zE1b4+PjAACbzQYgfq41FioLDA0Nwev1Ii8vL+D1vLw89PX1adQqutldd92Fn/zkJ3jrrbfwwx/+EH19fdi+fTuGh4f9OTFD/Qomo76+PpjNZmRlZS27DcXe7t278U//9E/45S9/ie9973v4zW9+g8985jP+p/EyN+0oioKjR49ix44dqK2tBRA/15roZ/1Ey82r3CqKwpVvdWT37t3+v996663Ytm0b1q9fj1dffdU/sI8Z6l84GTFHbe3du9f/99raWtx5550oLS3Ff/7nf+KRRx5Z9ueYW/R95StfwW9/+1ucPXt20fekX2vsUVkgOzsbRqNxURU5MDCwqCIl/UhNTcWtt96KK1eu+Gf/MEP9CiYju90Ot9uN0dHRZbch7eXn56O0tBRXrlwBwNy0cvjwYfz7v/873n77bRQVFflfj5drjYXKAmazGXV1dTh16lTA66dOncL27ds1ahWtxOVyobGxEfn5+SgvL4fdbg/I0O124/Tp08xQJ4LJqK6uDomJiQHbOBwOXL58mTnqyPDwMLq7u5Gfnw+AucWaoij4yle+gp/97Gf45S9/ifLy8oDvx821ptkwXp167bXXlMTEROUf/uEflIaGBuXIkSNKamqq0tHRoXXT6BNPPvmk8s477yhtbW3KuXPnlM9//vOK1Wr1Z/Sd73xHycjIUH72s58pH330kfLFL35Ryc/PVyYmJjRu+doxOTmpXLhwQblw4YICQHnuueeUCxcuKJ2dnYqiBJfRwYMHlaKiIuUXv/iF8uGHHyqf+cxnlM2bNysej0erf1bcU8ttcnJSefLJJ5X3339faW9vV95++21l27ZtSmFhIXPTyBNPPKFkZGQo77zzjuJwOPx/pqen/dvEw7XGQmUJf//3f6+UlpYqZrNZ2bJli3+qF+nD3r17lfz8fCUxMVEpKChQHnnkEeXjjz/2f9/n8yl/9Vd/pdjtdsVisSif/vSnlY8++kjDFq89b7/9tgJg0Z99+/YpihJcRjMzM8pXvvIVxWazKcnJycrnP/95paurS4N/zdqhltv09LRy//33Kzk5OUpiYqJSUlKi7Nu3b1EmzC12lsoKgPLjH//Yv008XGt8ejIRERHpFseoEBERkW6xUCEiIiLdYqFCREREusVChYiIiHSLhQoRERHpFgsVIiIi0i0WKkRERKRbLFSIiIhIt1ioEBERkW6xUCEiIiLdYqFCREREusVChYiIiHTr/wd428Q6OJoTsAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "semilogy(scipy.linalg.eigh(A.A, inv(Linv), eigvals_only=True), 'o-');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh4AAAGYCAYAAAANwusSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABhI0lEQVR4nO29eXRcV5Wo/5Wm0ixL1jxPlmVHaeM43R1nwCEQg4GQAI+E0WAcFsZgSNzkddz5NeRlkU6vpslzdwcH3iIhneGFNLwkNAsvggnBjuPkvcSxTRLJkmXNckmyLEsqzaqq+/vDXWWVNVgl1a17z6n9reVFXLpVtb9zrq43Z9jHYRiGgSAIgiAIQgSIsToAQRAEQRCiB0k8BEEQBEGIGJJ4CIIgCIIQMSTxEARBEAQhYkjiIQiCIAhCxJDEQxAEQRCEiCGJhyAIgiAIEUMSD0EQBEEQIkac1QFcis/nY2pqitjYWBwOh9XhCIIgCIKwCAzDwOv1kpCQQEzM/OMatks8pqameO2116wOQxAEQRCEJXDdddeRmJg4789tl3jExsYCUFdXF/jvcDIwMEBWVlbYP9cOiJua6OwGevuJm5ro7AbW+Xm9Xt59993L/tttu8TDP70SGxtrSuIRExNjyufaAXFTE53dQG8/cVMTnd3Aer/LLZOIusWlvb29VodgGuKmJjq7gd5+4qYmOruB/f2iLvEQBEEQBME6oi7xqK6utjoE0xA3NdHZDfT2Ezc10dkN7O8XdYlHd3e31SGYhripic5uoLefuKmJzm5gf7+oSzzGx8etDsE0xE1NdHYDvf3ETU10dgP7+0Vd4rHQ3mLVETc10dkN9PYTNzXR2Q3s7xd1iUdpaanVIZiGuKmJzm6gt5+4qYnObmB/v6hLPJqamqwOwTTETU10dgO9/cRNTXR2A/v72a6AmCAIgiAI4cPrM3i3Z4T+sSmGxj1kJMaTnRJPXX4qsTGRPxMt5MTj0KFD/PCHP+To0aO4XC5eeOEFbrvttkW997XXXmPTpk3U1dVx/PjxUL86LOTk5FjyvZFA3NREZzfQ20/c1ERnNwj2O9w6yL43uugfnZ51XXZKPDuvKeb6ihURjG4JUy2jo6OsW7eORx55JKT3DQ0NsXXrVj74wQ+G+pVhJS5O30EecVMTnd1Abz9xUxOd3eCi3+HWQR54uXXOpAOgf3SaB15u5XDrYASjW0LisWXLFn7wgx/wqU99KqT3ff3rX+fzn/88GzduDPUrw4rL5bL0+81E3NREZzfQ20/c1ERnN7jg5/UZ7Huja1HXP/pGF16fYXJUF4nI4tKf//znnD59mu9///uLfo/b7WZ4eDjwZ3Jy0sQIBUEQBEEf3u0ZmXek41LOjk7zbs+IyRFdxPTxplOnTnHvvffy6quvhjS8VVdXx9jYWODv27ZtY9euXRQUFHD69GkA8vLyMAyDvr4+AFatWkVXVxfj4+MkJiZSUlLCqVOnAMjNzSUmJgav10t9fT1VVVX09PQwOjqK0+mkvLycxsZGALKzs0lISODMmTMAVFRUcPbsWUZGRoiPj6e6upqGhgYAsrKySEpKClSKKy8vZ2BggOHhYWJjY1m9ejUNDQ0YhsGKFStIS0ujs7MTuLDlaXh4mMHBQRwOB2vWrKGxsRGv10t6ejqZmZm0t7cDUFxczNjYGAMDAwCsXbuWpqYmPB4PaWlpZGdnB9wKCwuZnJzk3LlzANTW1tLS0sLU1BQpKSnk5eXR0tICQEFBAR6Ph7NnzwJQU1NDR0cHExMTJCUlUVRURHNzc6C94eIBRNXV1XR3dwfau7S0NLCaOicnh7i4uMD/s6isrKS3t5fR0VESEhKorKzk5MmTAKxcuRKn0xnU3v39/bjdbuLi4qipqQm4ZWVlkZycTFfXhUy+rKyM8+fPz9ve6enpdHR0AFBSUoLb7Z63vbOysmhrawOgqKiI8fHxQHuvWbOG5uZmpqenSU1NJScnh9bWVgAKCwuZmpqiv78fgNWrV9PW1sbk5CQpKSnk5+cH7tn8/Hx8Pl/QPRsfH099fT1JSUkUFxcH3bMOhyPQ3lVVVbhcLsbGxnA6nZSVlS3Y3n19fYyMjMzZ3omJiXPes/72rq+vByAzM5PU1NSge3ZoaIihoSFiYmKora3l5MmT+Hw+MjIyyMjICGrvkZGRQN/NvGfnau+JiYk579nU1FRyc3MXvGfb29uZnJwkOTl52c+Inp6eQHtf7hnhd1PlGTHznr3cM8Lvpsozwn/PLuYZYRgG9fX1yjwjOjs7A+29mGeE1+uloS206qX9o1OBNlzqM8IwFjdq4jAWe+Vcb3Y4Flxc6vV6ueaaa9i+fTs7duwA4P777+fFF1+cd3Gpx+Ph4MGDVFZWEhNzcUDG6XTidDqXGmqA9vZ2ysrKlv05dkTc1ERnN9DbT9zURGc3uOA3GJ/FPfubF/2eH360mnWFacv6Xq/Xy4kTJ9i0adOCAw2mjni43W7eeustjh07xre+9S0AfD4fhmEQFxfH73//e2666aY535uWlkZsbGzYYxodHQ37Z9oFcVMTnd1Abz9xUxOd3eCCX11tKdkp8Yuabsn5r621kcLUxCM9PZ133nkn6LV9+/bxxz/+kV/96ldUVFSY+fVzkpCQEPHvjBTipiY6u4HefuKmJjq7wQW/2BgHO68p5oGXWy97/TeuKY5oPY+QE4+RkZHAfB5Aa2srx48fJysri9LSUvbs2UN3dzdPPvkkMTEx1NXVBb0/NzeXxMTEWa9HisrKSku+NxKIm5ro7AZ6+4mbmujsBhf9rq9Ywfc+WMG/vNbB0IR31nU5KfF8Q4U6Hm+99Rbr169n/fr1AOzevZv169fzve99D7iwjce/UMeO+Bcr6Yi4qYnObqC3n7ipic5uEOx3fcUK7tl0cT3LhqJU/nZTGT/8aDVP3nFFxJMOWMKIx4033rjgytUnnnhiwffff//93H///aF+rSAIgiAIS8Dru/jff1GQxgdXZVkXDFF4SNzKlSutDsE0xE1NdHYDvf3ETU10doPZfp4ZxcHiLTib5VKiLvEIx5ZcuyJuaqKzG+jtJ25qorMbzPbz+C4OecTFWv/PvvURRBh/ARodETc10dkN9PYTNzXR2Q1m+017L454xMmIhyAIgiAIZhI01RIriUfEsaJ2SKQQNzXR2Q309hM3NdHZDWb7zUw8Yh2SeEQcf618HRE3NdHZDfT2Ezc10dkNZvvNnGqREQ8LcLvdVodgGuKmJjq7gd5+4qYmOrvBbL+ZIx6yxsMCQjkhVzXETU10dgO9/cRNTXR2g9l+07LGw1pqamqsDsE0xE1NdHYDvf3ETU10doPZfp4ZFcRkxMMC6uvrrQ7BNMRNTXR2A739xE1NdHaD2X7BUy3W/7NvfQSCIAiCIJiGTLVYTFaWtTXqzUTc1ERnN9DbT9zURGc3mO3nkQJi1pKcnGx1CKYhbmqisxvo7SduaqKzG8z2k10tFtPV1WV1CKYhbmqisxvo7SduaqKzG8z2k6kWQRAEQRAihuxqsZiysjKrQzANcVMTnd1Abz9xUxOd3WC2n+di3kG87GqJPOfPn7c6BNMQNzXR2Q309hM3NdHZDWb7eXwzRjxkqiXyDA8PWx2CaYibmujsBnr7iZua6OwGs/1kcanFxMbGWh2CaYibmujsBnr7iZua6OwGs/2CDomTxCPyrF692uoQTEPc1ERnN9DbT9zURGc3mO03c8QjVhKPyNPQ0GB1CKYhbmqisxvo7SduaqKzG8z28yceMQ5JPCzBMIzLX6Qo4qYmOruB3n7ipiY6u8FsP/9Uix3Wd0AUJh4rVqywOgTTEDc10dkN9PYTNzXR2Q1m+/lHPCTxsIj09HSrQzANcVMTnd1Abz9xUxOd3WC2n387bXysPf7Jt0cUEaSjo8PqEExD3NREZzfQ20/c1ERnN5jtJ1MtgiAIgiBEDJlqsZiSkhKrQzANcVMTnd1Abz9xUxOd3WC2nz/xsMMBcRCFiYfb7bY6BNMQNzXR2Q309hM3NdHZDWb7yVSLxQwODlodgmmIm5ro7AZ6+4mbmujsBrP9ZKrFYhwOezS8GYibmujsBnr7iZua6OwGwX6GYchUi9WsWbPG6hBMQ9zURGc30NtP3NREZzcI9gs+IM4e/+TbI4oI0tjYaHUIpiFuaqKzG+jtJ25qorMbBPvZ7WRaiMLEw+v1Wh2CaYibmujsBnr7iZua6OwGwX5BJ9PKVIs16FyxTtzURGc30NtP3NREZzcI9vPKiIf1ZGVlWR2CaYibmujsBnr7iZua6OwGwX7TknhYT1tbm9UhmIa4qYnObqC3n7ipic5uEOwnazwEQRAEQYgYHh3WeBw6dIhbbrmFwsJCHA4HL7744oLXP//889x8883k5OSQnp7Oxo0beemll5Ya77IpKiqy7LvNRtzURGc30NtP3NREZzcI9pv+r5NpQeERj9HRUdatW8cjjzyyqOsPHTrEzTffzP79+zl69Cgf+MAHuOWWWzh27FjIwYaD8fFxS743EoibmujsBnr7iZua6OwGwX52rOMRF+obtmzZwpYtWxZ9/d69e4P+/g//8A/8+te/5je/+Q3r168P9euXzcDAAPn5+RH/3kggbmqisxvo7SduaqKzGwT72XGqJeTEY7n4fD7cbvdlVxW73W5iZmRnTqcTp9NpdniCIAiCoA123NUS8cTjRz/6EaOjo9x+++0LXldXV8fY2Fjg79u2bWPXrl0UFBRw+vRpAPLy8jAMg76+PgBWrVpFV1cX4+PjJCYmUlJSwqlTpwDIzc0lJiYGwzCor6+nqqqKnp4eRkdHcTqdlJeXB6q9ZWdnk5CQwJkzZwCoqKjg7NmzjIyMEB8fT3V1NQ0NDcCFbUtJSUl0d3cDUF5ezsDAAMPDw8TGxrJ69WoaGhowDIMVK1aQlpZGZ2cnAKWlpQwPDzM4OIjD4WDNmjU0Njbi9XpJT08nMzOT9vZ2AIqLixkbG2NgYACAtWvX0tTUhMfjIS0tjezs7IBbYWEhk5OTnDt3DoDa2lpaWlqYmpoiJSWFvLw8WlpaACgoKMDj8XD27FkAampq6OjoYGJigqSkJIqKimhubg60N0Bvby8A1dXVdHd3B9q7tLSUpqYmAHJycoiLi8PlcgFQWVlJb28vo6OjJCQkUFlZycmTJwFYuXIlTqczqL37+/txu93ExcVRU1MTcMvKyiI5OZmuri4AysrKOH/+/LztnZ6eTkdHB3DhqGi32z1ve2dlZQVWgxcVFTE+Ph5o7zVr1tDc3Mz09DSpqank5OTQ2toKQGFhIVNTU/T39wOwevVq2tramJycJCUlhfz8/MA9m5+fj8/nC7pnExMTqa+vJykpieLi4qB71uFwBNq7qqoKl8vF2NgYTqeTsrKyBdu7r6+PkZGROds7MTFxznvW39719fUAZGZmkpqaGnTPDg0NMTQ0RExMDLW1tZw8eRKfz0dGRgYZGRlB7T0yMhLou5n37FztPTExMec9m5qaSm5u7oL3bHt7O5OTkyQnJy/7GdHT0xNo78s9I/xuqjwjZt6zl3tG+N1UeUb479nFPCMcDgf19fXKPCM6OzsD7b2YZ4RhGJw+fZqysjJOt164PwA8UxOBdjLjGWEYF5OchXAYi71yrjc7HLzwwgvcdttti7r+2Wef5c477+TXv/41H/rQh+a8xuPxcPDgQSorK00Z8Th16hSrVq1a9ufYEXFTE53dQG8/cVMTnd0g2O/19iG+f+BC0r7t6gI+9z7zppi8Xi8nTpxg06ZNxMXNP64RsRGP5557ju3bt/PLX/5y3qRjJmlpacTGxoY9junp6bB/pl0QNzXR2Q309hM3NdHZDYL9oraOx7PPPstXvvIV/vf//t987GMfi8RXzktqaqql328m4qYmOruB3n7ipiY6u0Gwn8eG22lDHvEYGRkJzOcBtLa2cvz4cbKysigtLWXPnj10d3fz5JNPAheSjq1bt/Iv//IvXHPNNYH506SkJDIyMsKksXhycnIi/p2RQtzURGc30NtP3NREZzcI9gs+JM4e22lDjuKtt95i/fr1ga2wu3fvZv369Xzve98DwOVyBRbqAPz0pz/F4/HwzW9+k4KCgsCf73znO2FSCA3/Yh8dETc10dkN9PYTNzXR2Q2C/ew41RLyiMeNN9644MrVJ554Iujvf/rTn0L9CkEQBEEQwoAdEw97jLtEkMLCQqtDMA1xUxOd3UBvP3FTE53dINhv2oYFxKIu8ZiamrI6BNMQNzXR2Q309hM3NdHZDYL9ZMTDBvgLuOiIuKmJzm6gt5+4qYnObhDsZ8fKpVGXeAiCIAhCtOCVxMN6Vq9ebXUIpiFuaqKzG+jtJ25qorMbBPt5vBfreMgaD4vw19nXEXFTE53dQG8/cVMTnd0g2C94qsUe/+TbI4oIMjk5aXUIpiFuaqKzG+jtJ25qorMbBPsFLS6VEQ9rSElJsToE0xA3NdHZDfT2Ezc10dkNgv2CttPKGg9ryM8372Q+qxE3NdHZDfT2Ezc10dkNgv1kO60NOH36tNUhmIa4qYnObqC3n7ipic5uEOwnUy2CIAiCIEQMmWqxAToPsYmbmujsBnr7iZua6OwGl061XNxOK1MtFuGb0Qm6IW5qorMb6O0nbmqisxsE+82caomPtcc/+faIIoL09fVZHYJpiJua6OwGevuJm5ro7AbBfjOnWmTEQxAEQRAEU5FdLTZg1apVVodgGuKmJjq7gd5+4qYmOrtBsJ8/8YhxQKwkHtbQ2dlpdQimIW5qorMb6O0nbmqisxsE+/mnWuwy2gFRmHhMTExYHYJpiJua6OwGevuJm5ro7AbBfv4RD0k8LCQpKcnqEExD3NREZzfQ20/c1ERnNwj282+nlcTDQoqLi60OwTTETU10dgO9/cRNTXR2g2C/wFSLTaqWQhQmHqdOnbI6BNMQNzXR2Q309hM3NdHZDYL9/FMt8TH2+efePpEIgiAIghBWZI2HDcjNzbU6BNMQNzXR2Q309hM3NdHZDYL9AomHTLVYh8Nhn8YPN+KmJjq7gd5+4qYmOrtBsJ/H659qsY9z1CUevb29VodgGuKmJjq7gd5+4qYmOrvBRT/DMJiWqRZBEARBECLBjGNaZKrFSqqqqqwOwTTETU10dgO9/cRNTXR2g4t+096Lp9TKVIuFuFwuq0MwDXFTE53dQG8/cVMTnd3gol/wAXH2+efePpFEiLGxMatDMA1xUxOd3UBvP3FTE53d4KKfZ8Zci0y1WIjT6bQ6BNMQNzXR2Q309hM3NdHZDS76Tc8Y8ZCpFgspKyuzOgTTEDc10dkN9PYTNzXR2Q0u+gVPtUjiYRlNTU1Wh2Aa4qYmOruB3n7ipiY6u8FFv5lTLfEy1SIIgiAIgplM+y7uaomVEQ/ryMnJsToE0xA3NdHZDfT2Ezc10dkNLvp5ZI2HPYiLi7M6BNMQNzXR2Q309hM3NdHZDS76Be1qkcTDOnTevy1uaqKzG+jtJ25qorMbXPSbuaslLtY+/9zbJxJBEARBEJaNzzA4ccbNm51DgdfOjU7hnZGIWEnIicehQ4e45ZZbKCwsxOFw8OKLL172PQcPHmTDhg0kJiZSWVnJT37yk6XEGhYqKyst+26zETc10dkN9PYTNzXR2e1w6yAPnvBxz/5mfvnO2cDrf2g+z5eee4/DrYPWBfdfhJx4jI6Osm7dOh555JFFXd/a2spHP/pRbrjhBo4dO8bf/d3f8e1vf5v/83/+T8jBhoO+vj5LvjcSiJua6OwGevuJm5ro6na4dZAHXm7l3Jhnzp/3j07zwMutlicfIa+w2bJlC1u2bFn09T/5yU8oLS1l7969AKxZs4a33nqLf/7nf+bTn/50qF+/bEZGRiL+nZFC3NREZzfQ20/c1ERHN6/PYN8bXYu69tE3uthYlmHZFlvT13i8/vrrbN68Oei1D3/4w7z11ltMT0/P+z63283w8HDgz+TkZFjiSUhICMvn2BFxUxOd3UBvP3FTEx3d3u0ZoX90/n9TZ3J2dJp3e6xLvkzfU9TT00NeXl7Qa3l5eXg8Hvr7+ykoKJjzfXV1dUEH+Wzbto1du3ZRUFDA6dOnA59jGEZg2GzVqlV0dXUxPj5OYmIiJSUlnDp1CoDc3FxiYmKYnJykvr6eqqoqenp6GB0dxel0Ul5eTmNjIwDZ2dkkJCRw5swZACoqKjh79iwjIyPEx8dTXV1NQ0MDAFlZWSQlJdHd3Q1AeXk5AwMDDA8PExsby+rVq2loaMAwDFasWEFaWhqdnZ0AlJaWMjw8zODgIA6HgzVr1tDY2IjX6yU9PZ3MzEza29sBKC4uZmxsjIGBAQDWrl1LU1MTHo+HtLQ0srOzA26FhYVMTk5y7tw5AGpra2lpaWFqaoqUlBTy8vJoaWkBoKCgAI/Hw9mzF+YCa2pq6OjoYGJigqSkJIqKimhubg60N0Bvby8A1dXVdHd3B9q7tLQ0UDEvJyeHuLi4wOrqyspKent7GR0dJSEhgcrKSk6ePAnAypUrcTqdQe3d39+P2+0mLi6OmpqagFtWVhbJycl0dV3I7MvKyjh//vy87Z2enk5HRwcAJSUluN3ueds7KyuLtrY2AIqKihgfHw+095o1a2hubmZ6eprU1FRycnJobW0FoLCwkKmpKfr7+wFYvXo1bW1tTE5OkpKSQn5+fuCezc/Px+fzBd2zDoeD+vp6kpKSKC4uDrpnHQ5HoL2rqqpwuVyMjY3hdDopKytbsL37+voYGRmZs70TExPnvGf97V1fXw9AZmYmqampQffs0NAQQ0NDxMTEUFtby8mTJ/H5fGRkZJCRkRHU3iMjI4G+m3nPztXeExMTc96zqamp5ObmLnjPtre3Mzk5SXJy8rKfET09PYH2vtwzwu+myjNi5j17uWeE302VZ4T/nl3MM2J6epr6+nplnhGdnZ2B9p7vGfHns15CoWdojEzP+bA+IwxjcYtXHcZir5zrzQ4HL7zwArfddtu819TU1LBt2zb27NkTeO21117j+uuvx+VykZ+fH3S9x+Ph4MGDVFZWEjPjGF+n0xmWg338D0AdETc10dkN9PYTNzXR0e3EGTf37G9e9PU//Gg16wrTwhqD1+vlxIkTbNq0acFaKaaPeOTn5wf+H4Sfvr4+4uLiWLly5bzvS0tLIzY21uzwBEEQBEF56vJTyU6JX9R0S05KPHX5qRGIam5MX+OxceNGDhw4EPTa73//e66++mri4+PN/vpZLJTsqI64qYnObqC3n7ipiY5usTEOdl5TvKhrv3FNsaVnt4SceIyMjHD8+HGOHz8OXNgue/z48cAc2Z49e9i6dWvg+h07dtDe3s7u3btpaGjg8ccf57HHHuO73/1ueAxCJDEx0ZLvjQTipiY6u4HefuKmJrq6XV+xgntvLJv35zkp8XzvgxVcX7EickHNQchTLW+99RYf+MAHAn/fvXs3AF/+8pd54okncLlcgSQELiwC2r9/P3fffTc//vGPKSws5F//9V8t2UoL0N3dTUZGhiXfbTbipiY6u4HefuKmJjq71eamXPzvnCRurMwkIzGe7P+aXrHDKbUhJx433njjgitXn3jiiVmvbdq0ibfffjvUrxIEQRAEIQTOj11c47EmN5VPXZm3wNXWEHVntZSXl1sdgmmIm5ro7AZ6+4mbmujsdn78YtXSFUn2PIU36hIP/55rHRE3NdHZDfT2Ezc10dpt/OKIR1Zy5DdwLIaoSzyGh4etDsE0xE1NdHYDvf3ETU10dhucMeKRKSMe9mChoiaqI25qorMb6O0nbmqis9vMEY/MJBnxsAU1NTVWh2Aa4qYmOruB3n7ipiY6u81c45EliYc98NeX1xFxUxOd3UBvP3FTE53dZu5qyZCpFkEQBEEQzMQ/4pGRGEecDWp2zEXUJR6ZmZlWh2Aa4qYmOruB3n7ipia6uhmGwfn/WuNh14WlEIWJR2qqdQfjmI24qYnObqC3n7ipia5uY9M+prwXCnxK4mEjOjs7rQ7BNMRNTXR2A739xE1NdHU7r8COFojCxEMQBEEQdGRgbMaOFpsWD4MoTDxKS0utDsE0xE1NdHYDvf3ETU10dZs54mHXcukQhYnH0NCQ1SGYhripic5uoLefuKmJTm5en8GJM25eOT3AOy534PWRSS9e3/wHulqJJB4aIW5qorMb6O0nbmqii9vh1kG+9Nx73LO/mYdeaec/G84FfvaLE7186bn3ONw6aF2A8xB1iUdMjL7K4qYmOruB3n7ipiY6uB1uHeSBl1vpH52e95r+0WkeeLnVdsmHwzAMW43FeDweDh48yLp164iNjbU6HEEQBEGwFV6fwZd+8R79Y/MnHTPJSYnnyTuuINbkgmJer5cTJ06wadOmBc/DUT/tC5GTJ09aHYJpiJua6OwGevuJm5qo7Ob1GfzPV9sXnXQAnB2d5t2eEROjCg37Lns1CZ/PZ3UIpiFuaqKzG+jtJ25qoqKb12fw7PEe/uPPfUx4Qo9/5qm1VhN1iUdGRobVIZiGuKmJzm6gt5+4qYlKbstNOPzY6aRaSTw0QtzURGc30NtP3NREFbfDrYP8z8MduCe9y/qcnJR46vLtUyY+6tZ4dHR0WB2CaYibmujsBnr7iZua2N3N6zN4+m0XD7zcuuykA+Ab1xSbvrA0FKJuxEMQBEEQ7Ih/WuX5d/sYmVr+OpSclHi+cU0x11esWH5wYSTqEo+SkhKrQzANcVMTnd1Abz9xUxM7uoVrWiXdGcNX3pdNycp06vJTbTXS4SfqplpGRuyzpSjciJua6OwGevuJm5rYze1Qy/mwTavcdX0ZG7JjWFeYZsukA6Iw8Th//rzVIZiGuKmJzm6gt5+4qYmd3A61nOcfXmlb9uekO2P53gcruL5iha385iLqploEQRAEwQ4cbh3kB39sW9ZnJMU5+Mxf5PG59+XbdoTjUqRkuiAIgiBEmFDLnl+KHRMOKZk+D01NTVaHYBripiY6u4HefuKmJnZwe/Z4z5KSjjRnLFuvyuf5rev44lUFcyYddvBbiKibavF4PFaHYBripiY6u4HefuKmJla6+bfMPvl2T8jv3XpV/qJGOOzed1GXeKSnp1sdgmmIm5ro7AZ6+4mbmljldrh1kB+/3sm5sdASg3RnLHddX7roehx277uoSzyysrKsDsE0xE1NdHYDvf3ETU2scDvUcj7khaQO4Avr8/jC+rmnVObD7n0XdWs82trarA7BNMRNTXR2A739xE1NIu221C2z991UztYNhSEvHrV730XdiIcgCIIgRIqljHTAhfUc76/MDH9ANiDqRjyKioqsDsE0xE1NdHYDvf3ETU0i5bbUkY7s5Dg+9778JX+v3fsu6hKPiYkJq0MwDXFTE53dQG8/cVOTSLj5i4P5llApa+fGkmXV5rB730Vd4nHu3DmrQzANcVMTnd1Abz9xUxOz3bw+g32vd4X8vpyU+EDZ8+Vg976TNR6CIAiCEEaWUhxssTU6dGBJIx779u2joqKCxMRENmzYwKuvvrrg9c888wzr1q0jOTmZgoICtm3bZllGVltba8n3RgJxUxOd3UBvP3FTEzPdDrWcD6k4WIwD/r+byuetQroU7N53IScezz33HHfddRf33Xcfx44d44YbbmDLli10dHTMef3hw4fZunUr27dv57333uOXv/wlb775Jnfeeeeyg18KLS0tlnxvJBA3NdHZDfT2Ezc1McttKYtJ/+4D5WHfvWL3vgs58Xj44YfZvn07d955J2vWrGHv3r2UlJTw6KOPznn9G2+8QXl5Od/+9repqKjg+uuv5+tf/zpvvfXWsoNfClNTU5Z8byQQNzXR2Q309hM3NQm3m9dn8PTbrpAWk/pHOszYMmv3vgsp8ZiamuLo0aNs3rw56PXNmzdz5MiROd9z7bXX0tXVxf79+zEMg97eXn71q1/xsY99bMHvcrvdDA8PB/5MTk6GEuq8pKamhuVz7Ii4qYnObqC3n7ipSTjdDrcO8sVfvBvy2StmjHT4sXvfhbS4tL+/H6/XS15eXtDreXl59PTM3ejXXnstzzzzDHfccQcTExN4PB4+8YlP8G//9m8LflddXR1jY2OBv2/bto1du3ZRUFDA6dOnA99rGAZ9fX0ArFq1iq6uLsbHx0lMTKSkpIRTp04BkJubS0xMDENDQ9TX11NVVUVPTw+jo6M4nU7Ky8tpbGwEIDs7m4SEBM6cOQNARUUFZ8+eZWRkhPj4eKqrq2loaAAulKZNSkqiu7sbgPLycgYGBhgeHiY2NpbVq1fT0NCAYRisWLGCtLQ0Ojs7ASgtLWV4eJjBwUEcDgdr1qyhsbERr9dLeno6mZmZtLe3A1BcXMzY2BgDAwMArF27lqamJjweD2lpaWRnZwfcCgsLmZycDKyjqa2tpaWlhampKVJSUsjLywsMxRUUFODxeDh79iwANTU1dHR0MDExQVJSEkVFRTQ3NwfaG6C3txeA6upquru7A+1dWloaOBUxJyeHuLg4XC4XAJWVlfT29jI6OkpCQgKVlZWcPHkSgJUrV+J0OoPau7+/H7fbTVxcHDU1NQG3rKwskpOT6eq6sGK8rKyM8+fPz9ve6enpgWnAkpIS3G73vO2dlZUVqPhXVFTE+Ph4oL3XrFlDc3Mz09PTpKamkpOTQ2trKwCFhYVMTU3R398PwOrVq2lra2NycpKUlBTy8/MD92x+fj4+ny/onp2cnKS+vp6kpCSKi4uD7lmHwxFo76qqKlwuF2NjYzidTsrKyhZs776+PkZGRuZs78TExDnvWX9719fXA5CZmUlqamrQPTs0NMTQ0BAxMTHU1tZy8uRJfD4fGRkZZGRkBLX3yMhIoO9m3rNztffExMSc92xqaiq5ubkL3rPt7e1MTk4G1pEt5xnhf5Yt5hnhd1PlGTHznr3cM8Lvpsozwn/PLuYZ4Xa7qa+vX/Yz4lDLAP92dJBQ+dL78igyBqivdy3qGdHZ2Rlo78U8I9xuN6dPn474M8IwFjfc4zAWeyVw5swZioqKOHLkCBs3bgy8/uCDD/LUU08Fgp5JfX09H/rQh7j77rv58Ic/jMvl4p577uEv//Iveeyxx2Zd7/F4OHjwIJWVlcTEXByQcTqdOJ3OxYY6L/4HoI6Im5ro7AZ6+4mbmoTDzesz+NIv3gt590p2chxPfbbO1N0rVvWd1+vlxIkTbNq0ibi4+cc1QhrxyM7OJjY2dtboRl9f36xRED8PPfQQ1113Hffccw8Af/EXf0FKSgo33HADP/jBDygoKJjzfWlpacTGxoYSniAIgiBEhKVsmYXlFwfTgZDWeCQkJLBhwwYOHDgQ9PqBAwe49tpr53zP2NhY0MgFEEgoQhhsCRvzJTo6IG5qorMb6O0nbmqyXLdQt8zCxcWkyy0Othjs3nch72rZvXs3P/vZz3j88cdpaGjg7rvvpqOjgx07dgCwZ88etm7dGrj+lltu4fnnn+fRRx+lpaWF1157jW9/+9v81V/9FYWFheEzWSQejyfi3xkpxE1NdHYDvf3ETU2W47bU81fMXEx6KXbvu5ATjzvuuIO9e/fywAMP8L73vY9Dhw6xf/9+ysrKAHC5XEE1Pb7yla/w8MMP88gjj1BXV8dnPvMZVq9ezfPPPx8+ixDwL5DSEXFTE53dQG8/cVOTpbr5T5oN5fwVfxn0SJ40a/e+W1LJ9J07d7Jz5845f/bEE0/Mem3Xrl3s2rVrKV8lCIIgCJazlJGOaCqDHgoh7WqJBP5dLevWrTNlcanH41lwta3KiJua6OwGevuJm5qE6na4dZAHXm4N6Tu2XpXPF6+yZq2FVX232F0tUXc6rX/Pu46Im5ro7AZ6+4mbmoTitpSTZrOT4/jc+/JDDSts2L3voi7xCFcFVDsibmqisxvo7SduahKK21K2zVq9ZdbufafnONoCJCcnWx2CaYibmujsBnr7iZuaLMbN6zN49nhPyCfN/t0HIrNldiHs3ndRl3jYfX/zchA3NdHZDfT2Ezc1uZzb4dZBfvx6J+fGQtuWGsktswth976LuqkWfz18HRE3NdHZDfT2Ezc1WcjtUMt5Hni5NaSkw8yTZpeC3fsu6hIPQRAEQZgLFYqD6UDUTbXMd6aMDoibmujsBnr7iZuazOXmLw4WKluvyrdd0mH3vou6EQ+blS0JK+KmJjq7gd5+4qYml7otdaTD6m2z82H3vou6xKOvr8/qEExD3NREZzfQ20/c1GSm21LKoPuxetvsfNi976JuqkUQBEEQYOkjHTkp8XzjmmLLt82qStSVTJ+eniY+Pj7sn2sHxE1NdHYDvf3ETU2mp6f5v12jIZdBBzXOX7Gq76Rk+jx0dYVW+lYlxE1NdHYDvf3ETU06OjtDLoPu3zL7xasKbJ10gP37LuqmWsbHx60OwTTETU10dgO9/cRNTX59yk3/mC+k96i0ZdbufRd1iUdiYqLVIZiGuKmJzm6gt5+4qYW/DPrvOhefdPjLoKuSdID9+y7qEo+SkhKrQzANcVMTnd1Abz9xUwfVy6CHgt37LurWeJw6dcrqEExD3NREZzfQ20/c1ECHMuihYPe+i7rEQxAEQYgepAy6/Yi6qZbc3FyrQzANcVMTnd1Abz9xszc6lUEPBbv3XdSNeMTE6Kssbmqisxvo7Sdu9kW3MuihYPe+s3d0JtDT02N1CKYhbmqisxvo7Sdu9kTHMuihYPe+i7qpFkEQBEFfpAy6/Ym6xKOqqsrqEExD3NREZzfQ20/c7IO/RseTb4f+//ZVKIMeCnbvO5lq0QhxUxOd3UBvP3GzB4dbB/niL94NOelQqQx6KNi976JuxGN0dNTqEExD3NREZzfQ20/crGepO1dA3y2zdu+7qEs8nE6n1SGYhripic5uoLefuFnLUtdzOID7FC0Othjs3ndRl3iUl5dbHYJpiJua6OwGevuJm3UsZ6Rjz42l2iYdYP++i7o1Ho2NjVaHYBripiY6u4HefuIWebw+gyePnuHBJSQd/jUduVO94Q/MRti17/xE3YiHIAiCoCaHWwf5n4c7cE96l/R+/5qO+npXmCMTQiHqEo/s7GyrQzANcVMTnd1Abz9xixzLmVq5tEaH3dzCjd39oi7xSEhIsDoE0xA3NdHZDfT2E7fIsNRFpDB3jQ47uZmB3f2ibo3HmTNnrA7BNMRNTXR2A739xM1cvD6Dp992Lan8+UI1OuzgZiZ294u6EQ9BEATB/hxuHeTHr3dybsyzpPfrWqNDB6Iu8aioqLA6BNMQNzXR2Q309hM3c1jOeo4Yx+WTDp37DezvF3VTLWfPnrU6BNMQNzXR2Q309hO38LOc9RywuJEOnfsN7O8XdSMeIyMjVodgGuKmJjq7gd5+4hZeljPSke6M5a7rSxd1uqzO/Qb294u6xCM+Pt7qEExD3NREZzfQ20/cwsdyyp9/YX0eX1i/+IPedO43sL/fkqZa9u3bR0VFBYmJiWzYsIFXX311wesnJye57777KCsrw+l0UlVVxeOPP76kgJdLdXW1Jd8bCcRNTXR2A739xC08+Ec6Qt25AhfOXNm6oTCk02V17jewv1/Iicdzzz3HXXfdxX333cexY8e44YYb2LJlCx0dHfO+5/bbb+fll1/mscceo7GxkWeffZba2tplBb5UGhoaLPneSCBuaqKzG+jtJ27LYznlz3NS4vneByuWtHNF534D+/uFPNXy8MMPs337du68804A9u7dy0svvcSjjz7KQw89NOv63/3udxw8eJCWlhaysrIA+x9gIwiCIJjLcsqfz1UUTFCHkEY8pqamOHr0KJs3bw56ffPmzRw5cmTO9/znf/4nV199Nf/0T/9EUVERNTU1fPe732V8fHzB73K73QwPDwf+TE5OhhLqvPiTHx0RNzXR2Q309hO30PEXBXvg5daQk46FioKFgs79Bvb3C2nEo7+/H6/XS15eXtDreXl59PT0zPmelpYWDh8+TGJiIi+88AL9/f3s3LmTgYGBBdd51NXVMTY2Fvj7tm3b2LVrFwUFBZw+fTrwvYZh0NfXB8CqVavo6upifHycxMRESkpKOHXqFAC5ubnExMTQ29vLwMAAVVVV9PT0MDo6itPppLy8PHCiX3Z2NgkJCYHqbxUVFZw9e5aRkRHi4+Oprq4ODGVlZWWRlJREd3c3cGE0Z2BggOHhYWJjY1m9ejUNDQ0YhsGKFStIS0ujs7MTgNLSUoaHhxkcHMThcLBmzRoaGxvxer2kp6eTmZlJe3s7AMXFxYyNjTEwMADA2rVraWpqwuPxkJaWRnZ2dsCtsLCQyclJzp07B0BtbS0tLS1MTU2RkpJCXl4eLS0tABQUFODxeALbr2pqaujo6GBiYoKkpCSKiopobm4OtDdAb++Fkx2rq6vp7u4OtHdpaSlNTU0A5OTkEBcXh8t14TCmyspKent7GR0dJSEhgcrKSk6ePAnAypUrcTqdQe3d39+P2+0mLi6OmpqagFtWVhbJycl0dXUBUFZWxvnz5+dt7/T09MA0YElJCW63e972zsrKoq2tDYCioiLGx8cD7b1mzRqam5uZnp4mNTWVnJwcWltbASgsLGRqaor+/n4AVq9eTVtbG5OTk6SkpJCfnx+4Z/Pz8/H5fEH37Pnz5xkYGCApKYni4uKge9bhcATau6qqCpfLxdjYGE6nk7KysgXbu6+vj5GRkTnbOzExcc571t/e9fX1AGRmZpKamhp0zw4NDTE0NERMTAy1tbWcPHkSn89HRkYGGRkZQe09MjJCX18fAwMDQffsXO09MTEx5z2bmppKbm7ugvdse3s7k5OTJCcnL/sZ4X+WLeYZ4b8vVXlGzLxnL/eM8LuF8xnRaWTwyOEOhqZZEl9aFUv2hIvx8cRZzwj/PbuYZ4T/nlTlGdHZ2Rlo78U8I4aGhhgdHY34M8IwFrdIx2Es9koulGEtKiriyJEjbNy4MfD6gw8+yFNPPRUIeiabN2/m1Vdfpaenh4yMDACef/55/tt/+2+Mjo6SlJQUdL3H4+HgwYNUVlYSE3NxQMbpdOJ0Ohcb6rzU19ezdu3aZX+OHRE3NdHZDfT2E7fFY3ZRsFDQud/AOj+v18uJEyfYtGkTcXHzj2uENOKRnZ1NbGzsrNGNvr6+WaMgfgoKCigqKgokHXAhMzQMg66uLlatWjXn+9LS0oiNjQ0lPEEQBMGGRKIomKAOIa3xSEhIYMOGDRw4cCDo9QMHDnDttdfO+Z7rrruOM2fOBBU0aWpqIiYmhuLi4iWEvDx0Xtgqbmqisxvo7Sdul2c5W2X9azrCnXTo3G9gf7+Qt9Pu3r2bn/3sZzz++OM0NDRw991309HRwY4dOwDYs2cPW7duDVz/+c9/npUrV7Jt2zbq6+s5dOgQ99xzD1/96ldnTbNEAv98nI6Im5ro7AZ6+4nb/Cxnq6wfs0Y6dO43sL9fyNtp77jjDs6dO8cDDzyAy+Wirq6O/fv3U1ZWBoDL5Qqq6ZGamsqBAwfYtWsXV199NStXruT222/nBz/4QfgsQmB4eNiS740E4qYmOruB3n7iNjfL2SoLoZU/Xwo69xvY329JJdN37tzJzp075/zZE088Meu12traWdMzVqHzuhFxUxOd3UBvP3GbzXIWkaY5Y/nkFTmm1+jQud/A/n4h7WqJBP5dLevWrbN94wmCIAgX8S8iXcp6DikKpj6L3dWypLNaVMbupWSXg7ipic5uoLefuF3AXxRsKYtIw1UULBR07jewv1/UnU5rswGesCJuaqKzG+jtJ24X1nP8+PVOzo15lvQ9VmyV1bnfwP5+UZd4rFixwuoQTEPc1ERnN9DbL9rd7FQULBR07jewv1/UJR5paWlWh2Aa4qYmOruB3n7R6ub1GTxzzMUzx3qX/PlWFgXTud/A/n5Rt8bDfwaCjoibmujsBnr7RaPb4dZBbn/mHZ4+1stSBvTTnbFLPs4+XOjcb2B/v6gb8RAEQRCWxnKmVhzAF9bn8YX1kVtEKtiTqEs8SktLrQ7BNMRNTXR2A739osltueet3GdC6fOlonO/gf39om6qxe4V3ZaDuKmJzm6gt1+0uC3nvJWclHjLp1YuRed+A/v7RV3iMTg4aHUIpiFuaqKzG+jtp7PbwPnznDjj5tHXO5d83srWq/J58o4rTCt9vlR07jewv1/UTbU4HPrOLYqbmujsBnr76ep2uHWQfz06xeBU85Leb+VW2cWga7/5sbuflEwXBEEQAixnAakfM46yF+yPlEyfh8bGRqtDMA1xUxOd3UBvP53cwnGMvR22yi4GnfptLuzuF3VTLV7v0o5pVgFxUxOd3UBvP13clnuMvWpbZXXpt/mwu1/UJR7p6elWh2Aa4qYmOruB3n46uIVjasVOW2UXgw79thB294u6xCMzU51fjlARNzXR2Q309lPdbbm1Oey+iHQ+VO+3y2F3v6hb49He3m51CKYhbmqisxvo7aeq23KOsZ+JikkHqNtvi8XuflE34iEIghCteH0Gzx7v4fl3+xiZ8i35c3JS4vnGNcW2q88hqEHUJR7FxcVWh2Aa4qYmOruB3n4quYVjAemta7O5rnwFdfmpSiwinQ+V+m0p2N0v6qZaxsbGrA7BNMRNTXR2A739VHE71HKeB15uXXLSARcWkO68toR1hWlKJx2gTr8tFbv7RV3iMTAwYHUIpiFuaqKzG+jtZ3e3aKrNEQp277flYne/qJtqEQRBiAaWO7UC8EWFanMI6iAl0wVBEDRDyp4LViAl0+ehqanJ6hBMQ9zURGc30NvPTm5en7Hs02Th4jH2+Z6z4QvOZtip38zA7n5RN9Xi8XisDsE0xE1NdHYDvf3s4na4dZB9b3TRPzq9rM/ZelU+n3tfPrExDurrz4QpOvthl34zC7v7RV3ikZaWZnUIpiFuaqKzG+jtZwe3cEyrpDtjuev60qC6HHZwMwud3cD+flGXeGRnZ1sdgmmIm5ro7AZ6+1np5vUZPHPMxTPHepf8GQsd7ib9pi5294u6NR6tra1Wh2Aa4qYmOruB3n5WuR1uHeT2Z97h6WO9LGd3wH03lbN1Q+Gcu1ak39TF7n5RN+IhCIKgMmZNrQhCpIi6xKOwsNDqEExD3NREZzfQ2y+SbmZPrVyK9Ju62N0v6hKPyclJq0MwDXFTE53dQG+/SLmFoxgYXJhaWWxtDuk3dbG7X9St8Th37pzVIZiGuKmJzm6gt5/Zbv7j65d7zoq/NkcoBcGk39TF7n5RN+IhCIJgd8JxfL1Op8kKehF1JdN9Ph8xMXoO9IibmujsBnr7meEWrmmV5ZY8l35TF6v8pGT6PLS0tFgdgmmIm5ro7AZ6+4XbLRzH14frNFnpN3Wxu1/UTbVMTU1ZHYJpiJua6OwGevuFyy3SO1YWg/SbutjdL+oSj5SUFKtDMA1xUxOd3UBvv3C4WbFjZTFIv6mL3f2WNNWyb98+KioqSExMZMOGDbz66quLet9rr71GXFwc73vf+5bytWEhLy/Psu82G3FTE53dQG+/5brZaWrlUqTf1MXufiEnHs899xx33XUX9913H8eOHeOGG25gy5YtdHR0LPi+oaEhtm7dygc/+MElBxsO7D73tRzETU10dgO9/ZbiFq7j69OcsWy9Kp/nvnClKRVIpd/Uxe5+IU+1PPzww2zfvp0777wTgL179/LSSy/x6KOP8tBDD837vq9//et8/vOfJzY2lhdffHHJAQuCIKiIf4vsC++dXfa0yszj6wVBNUIa8ZiamuLo0aNs3rw56PXNmzdz5MiRed/385//nNOnT/P9739/0d/ldrsZHh4O/AlXJbaCgoKwfI4dETc10dkN9PZbrJv/ULcn3+4Jy7TKF68KzwLShZB+Uxe7+4U04tHf34/X6501f5SXl0dPT8+c7zl16hT33nsvr7766oL7ei+lrq6OsbGxwN+3bdvGrl27KCgo4PTp04HvNQyDvr4+AFatWkVXVxfj4+MkJiZSUlLCqVOnAMjNzSUmJob29nZcLhdVVVX09PQwOjqK0+mkvLycxsZG4MKRwgkJCZw5cwaAiooKzp49y8jICPHx8VRXV9PQ0ABAVlYWSUlJdHd3A1BeXs7AwADDw8PExsayevVqGhoaMAyDFStWkJaWRmdnJwClpaUMDw8zODiIw+FgzZo1NDY24vV6SU9PJzMzk/b2dgCKi4sZGxtjYGAAgLVr19LU1ITH4yEtLY3s7Gza2tpwuVwUFhYyOTkZqF5XW1tLS0sLU1NTpKSkkJeXFxiKKygowOPxcPbsWQBqamro6OhgYmKCpKQkioqKaG5uDrQ3QG/vhZX31dXVdHd3B9q7tLSUpqYmAHJycoiLi8PlcgFQWVlJb28vo6OjJCQkUFlZycmTJwFYuXIlTqczqL37+/txu93ExcVRU1MTcMvKyiI5OZmuri4AysrKOH/+/LztnZ6eHpgGLCkpwe12z9veWVlZtLW1AVBUVMT4+HigvdesWUNzczPT09OkpqaSk5MTOAGysLCQqakp+vv7AVi9ejVtbW1MTk6SkpJCfn5+4J7Nz8/H5/MF3bNnzpzB5XKRlJREcXFx0D3rcDgC7V1VVYXL5WJsbAyn00lZWdmC7d3X18fIyMic7Z2YmDjnPetv7/r6egAyMzNJTU0NumeHhoYYGhoiJiaG2tpaTp48ic/nIyMjg4yMjKD2HhkZoaenB5fLFXTPztXeExMTc96zqamp5ObmLnjPtre3Mzk5SXJy8rKfEf5n2WKeEf77cqFnRPNEEo+8Pchy+UhJDN+5+UqaTzVRX39mSc+Imffs5Z4RfjdVnhH+e3Yxzwj/vwOqPCM6OzsD7b2YZ8T58+cZGBiI+DNisWXBQiogdubMGYqKijhy5AgbN24MvP7ggw/y1FNPBYL24/V6ueaaa9i+fTs7duwA4P777+fFF1/k+PHjc36Hv4BYZWVlUAEUp9OJ0+lcbKjzUl9fz9q1a5f9OXZE3NREZzfQ228ht5lbZJdbpXG5xcCWQrT2mw5Y5bfYAmIhjXhkZ2cTGxs7a3Sjr69vzlW0brebt956i2PHjvGtb30LuFBRzTAM4uLi+P3vf89NN90053elpaWZUrlUEATBbMK1RVaOrxd0JKTEIyEhgQ0bNnDgwAE++clPBl4/cOAAt95666zr09PTeeedd4Je27dvH3/84x/51a9+RUVFxRLDXjo1NTUR/85IIW5qorMb6O13qZt/AemTb8899bxYwl0MbClEU7/pht39Qt5Ou3v3bn72s5/x+OOP09DQwN13301HR0dgKmXPnj1s3br1wofHxFBXVxf0Jzc3l8TEROrq6iwpcnK5bb8qI25qorMb6O3nd/OfIvuZp/+87KQDLhQD27qh0NJdK9HQb7pid7+Qt9PecccdnDt3jgceeACXy0VdXR379++nrKwMAJfLZWvpiYkJq0MwDXFTE53dQG+/sfFxnn7bxX/8uY8Jz9JOkZ1JTko837im2BZTKzr3m85uYH+/JZVM37lzJzt37pzzZ0888cSC773//vu5//77l/K1YSEpKcmy7zYbcVMTnd1AX7/DrYP86M0pRj3LG+FIjo/hIzUr2ViWYavj63XtN9DbDezvF3VntRQVFVkdgmmIm5ro7AZ6+h1qOc8PllF1FOyxjmMhdOw3Pzq7gf39lnRWi8r495vriLipic5uoJef12fw5NEzyyp17scO6zgWQqd+uxSd3cD+flE34iEIghAKXp/Buz0jvN4xxO8azzE2vby1HLJFVoh2oi7xsPupfctB3NREZzdQ2+9w6yD73uiif3R62Z+V5ozlk1fkKHPGisr9djl0dgP7+0Vd4iEIgrAYwrGOw48c6iYIF4m6NR7+mvY6Im5qorMbqOcXznUckTzULdyo1m+hoLMb2N9PRjwEQRC4WHU0HDU5kuIcfOYv8mSUQxDmIOoSj+rqaqtDMA1xUxOd3UANv3CdrWL3LbKhoEK/LRWd3cD+flE31eI/5ldHxE1NdHYD+/sdajnPAy+3LjvpAPtvkQ0Fu/fbctDZDezvF3UjHuPj41aHYBripiY6u4E9/fxbZI+0D/Lie/3L/jwdt8jasd/Chc5uYH+/qEs8EhMTrQ7BNMRNTXR2A3v5+ddxvPDe2bCMcKTEwaf/Qs8dK3bqt3CjsxvY3y/qEo/S0lKrQzANcVMTnd3AHn7hXDjqAG5dm8115SuozU7EmRAfniBthh36zSx0dgP7+0XdGo+mpiarQzANcVMTnd3Aer/DrYPc/sw7PPl2T1hOkL3vpnJ2XlvCusI0TjefCkOE9sTqfjMTnd3A/n5RN+IhCEL0EM4iYDqu4xAEK4i6xCMnJ8fqEExD3NREZzewxs/rM3jmmItnji2/kNJCNTl07jtxUxe7+0Vd4hEXp6+yuKmJzm4QOb9wH+a2mJocOveduKmL3f3sHZ0JuFwuMjMzrQ7DFMRNTXR2A/P9wr1Txc99N5Xz/sqF49a578RNXezuF3WJhyAIehDOnSozyUmJ5xvXFMtaDkEwiahLPCorK60OwTTETU10dgNz/MJV4txPcnwMH6lZycayDOryUxddk0PnvhM3dbG7X9Rtp7X7qX3LQdzURGc3CL9fOEucJ8U52HpVPv/nS3/Bjo3FrCtMC6kQmM59J27qYne/qBvxGB0dtToE0xA3NdHZDcLjF+4S5+E6zE3nvhM3dbG7X9QlHgkJCVaHYBripiY6u8Hy/KxcOLoYdO47cVMXu/s5DMMwrA5iJh6Ph4MHD7Ju3TpiY2PD/vk+n4+YGD1nmMRNTXR2g6X5mbVwNNxFwHTuO3FTF6v8vF4vJ06cYNOmTQtu6dW35efh5MmTVodgGuKmJjq7QWh+Xp/B02+7+OSTfw5biXOANGcsW6/K57kvXBnW3So69524qYvd/aJuqkUQBPthxgjHUneqCIJgLlGXeKxcudLqEExD3NREZze4vF+4t8YuVOI83Ojcd+KmLnb3i7rEw+l0Wh2CaYibmujsBvP7+Uc5nny7JyzfE66dKqGgc9+Jm7rY3S/q1nicOXPG6hBMQ9zURGc3mO3nX8fxmaf/HLakAy7sVNm6oTCiUyo69524qYvd/aJuxEMQBGuQEueCIEAUJh4VFRVWh2Aa4qYmOrsBlJaV8/TbLi0Xjurcd+KmLnb3i7qplv7+5Vc9tCvipia6uvmnVD73H41h2xq73BLn4UbXvgNxUxm7+0XdiIfb7bY6BNMQNzXRyc1f2vz1jiF+13iOsenwjHBEcqdKKOjUd5cibupid7+oSzwWqqamOuKmJjq4mVXa3K4Jhx8d+m4+xE1d7O4XdSXTBUEIH2YtGE1zxvLJK3Jsm3AIgjAbKZk+D/X19VaHYBripiYquplV2hxg61X5/McXruSLV0WuHsdSUbHvFou4qYvd/ew9HiMIgq0wa4QDwn+AmyAI9iTqEo+srCyrQzANcVMTu7uZtWDUj93XcSyE3ftuOYibutjdb0lTLfv27aOiooLExEQ2bNjAq6++Ou+1zz//PDfffDM5OTmkp6ezceNGXnrppSUHvFySk5Mt+26zETc1saubfzrl9mfe4Z79zTz/7tmwJh3+rbHPb12nxLTKXNi178KBuKmL3f1CTjyee+457rrrLu677z6OHTvGDTfcwJYtW+jo6Jjz+kOHDnHzzTezf/9+jh49ygc+8AFuueUWjh07tuzgl0JXV5cl3xsJxE1N7OZ26fqNcO5SAT0SDj9267twIm7qYne/kKdaHn74YbZv386dd94JwN69e3nppZd49NFHeeihh2Zdv3fv3qC//8M//AO//vWv+c1vfsP69euXFrUgCGHHzPUbIDtVBEG4QEiJx9TUFEePHuXee+8Nen3z5s0cOXJkUZ/h8/lwu92XnYNyu93ExFwckHE6nWE5ca+srGzZn2FXxE1NrHQze/1GcnwMN1Wms6kq29LS5mYh96Wa6OwG9vcLKfHo7+/H6/WSl5cX9HpeXh49PYs7ZfJHP/oRo6Oj3H777QteV1dXx9jYWODv27ZtY9euXRQUFHD69OnA9xqGQV9fHwCrVq2iq6uL8fFxEhMTKSkp4dSpUwDk5uYSExNDa2srSUlJVFVV0dPTw+joKE6nk/LychobGwHIzs4mISEhcMJfRUUFZ8+eZWRkhPj4eKqrq2loaAAuLOJJSkqiu7sbgPLycgYGBhgeHiY2NpbVq1fT0NCAYRisWLGCtLQ0Ojs7ASgtLWV4eJjBwUEcDgdr1qyhsbERr9dLeno6mZmZtLe3A1BcXMzY2BgDAwMArF27lqamJjweD2lpaWRnZ3Py5EmSkpIoLCxkcnKSc+fOAVBbW0tLSwtTU1OkpKSQl5dHS0sLAAUFBXg8Hs6ePQtATU0NHR0dTExMkJSURFFREc3NzYH2Bujt7QWgurqa7u7uQHuXlpbS1NQEQE5ODnFxcbhcLgAqKyvp7e1ldHSUhIQEKisrOXnyJAArV67E6XQGtXd/fz9ut5u4uDhqamoCbllZWSQnJweGEsvKyjh//vy87Z2enh6YBiwpKcHtds/b3llZWbS1tQFQVFTE+Ph4oL3XrFlDc3Mz09PTpKamkpOTQ2trKwCFhYVMTU0FyhSvXr2atrY2JicnSUlJIT8/P3DP5ufn4/P5gu7Z06dPExcXR1JSEsXFxUH3rMPhCLR3VVUVLpeLsbExnE4nZWVlC7Z3X18fIyMjc7Z3vNPJE/+3g0MuL2OeeX4Jl4EzBm5bm8lHShPocbmIHxwntvDiPTtXe09MTMx5z6amppKbm7vgPdve3s7k5CTJycnLfkb4n2WLeUacPn2apKQkZZ4RM+/Zyz0j/L9zqjwj/FtIF/OMaGxsJDExUZlnRGdnZ6C9F/OMOHfuHCtWrFjWMyIxMXHOe/bS9s7MzCQ1NZXOzk4WWxYspAJiZ86coaioiCNHjrBx48bA6w8++CBPPfVUIOj5ePbZZ7nzzjv59a9/zYc+9KE5r/EXEKusrDRlxKO+vp61a9cu+3PsiLipSSTdzJ5OmWuHivSdmoibuljlt9gCYiGNeGRnZxMbGztrdKOvr2/WKMilPPfcc2zfvp1f/vKX8yYdM0lLSzOlcqnO1VDFTU0i4WZFwuFH+k5NxE1d7O4Xcsn0v/7rv2bDhg3s27cv8NratWu59dZb51xcChdGOr761a/y7LPPctttty34+VIyXRDChywYFQQhUphWMn337t387Gc/4/HHH6ehoYG7776bjo4OduzYAcCePXvYunVr4Ppnn32WrVu38qMf/YhrrrmGnp4eenp6GBoaWoLW8vHPu+qIuKlJuN28PoMTZ9z85I0uPv1U+EuaJ8fH8KkrcvjhR6sXVdpc+k5NxE1d7O4X8nbaO+64g3PnzvHAAw/gcrmoq6tj//79gVW0LpcrqKbHT3/6UzweD9/85jf55je/GXj9y1/+Mk888cTyDULEZmfihRVxU5NwuZl1QqyfpVYYlb5TE3FTF7v7Lalk+s6dO9m5c+ecP7s0mfjTn/60lK8wjRUrVlgdgmmIm5os183K9RuLQfpOTcRNXezuF3VntaSnp1sdgmmIm5qE6uavvdE/NsWxbjeHWodsvX5D+k5NxE1d7O4XdYlHR0eHttuoxE1NFuM2s9DXy80DDE2EfyoFLqzf+EjNSjaWZYSt4Fe0952qiJu62N0v6hIPQVAJs9dt+FH5hFhBENQi6hKPkpISq0MwDXFTk7nczF634ScSCUe09Z0uiJu62N0v6hIPt9tNWlqa1WGYgripyUw3nRIOP9HSd7ohbupid7+Q63iozuDgoNUhmIa4qcng4OCso+jNSjrSnLERP5Je977TFXFTF7v7Rd2Ih8Oh7/y1uKnBzF0pQ+MeGjqmefPNP4f9ZFg/ZiwYDQWd+u5SxE1NdHYD+/uFXDLdbKRkuqArkVoo6kcWjAqCEElMK5muOv5jrXVE3OzJpdMoZicdSXGOiE+nLITKfXc5xE1NdHYD+/tF3VSL12v+/9O0CnGznkunUXrc0/z+1DnTplFmYtcD21Tpu6UgbmqisxvY3y/qEg+7V3RbDuJmHZGeRkmOj+HDq7LIS0sgIzGe7JR4S9ZvLAa7991yEDc10dkN7O8XdYlHVlaW1SGYhrhFnkhtf/Wj4roNu/ZdOBA3NdHZDezvF3VrPNra2qwOwTTEzVz8x82/cnqAY93DPHX0jOnbX/3Ybd1GKNih78xC3NREZzewv1/UjXgIQqhEehplJiqOcAiCICxE1CUeRUVFVodgGuIWHiJ1+utC2HWh6FKQ+1JNxE1d7O4XdYnH+Pg4GRkZVodhCuK2PKwc2QDISIzlg1VZlhX6Mgu5L9VE3NTF7n5Rl3gMDAyQn59vdRimIG6Lx8ptr34SY+H9lZnkMsy6VRVaJRszkftSTcRNXezuF3WJhxDdWD2qAbOnUerr61lbaN8DnQRBEMJJ1JVMNwzD9nXsl4q4zcYO6zVg4WkUnfsN9PYTNzXR2Q2s81tsyfSoG/Fobm5m1apVVodhCuJ2AX+y8XrHEC83DzA0EfmRDf80yvrC9MsW99K530BvP3FTE53dwP5+UZd4TE9PWx2CaUSr26WjGq+3D+OesmYaZSnbX3XuN9DbT9zURGc3sL9f1CUeqampVodgGtHg5k8yBsanWZEYx7s9I7xY32/Zeg0/y9n+qnO/gd5+4qYmOruB/f2iLvHIycmxOgTT0NHNn2icGUngPw62WTqa4ScpzsGnr8zlyvxUBic8ZCUt75wUHfttJjr7iZua6OwG9veLusSjtbWVtWvXWh2GKejkZofdJ5diVlEvnfptLnT2Ezc10dkN7O8XdYmHYD/mmj755TtnLdl9MpPUhBiuLctY1AJRQRAEYXFEXeJRWFhodQimoYqbnRaDXooV1UNV6belorOfuKmJzm5gf7+oSzympqasDsE07OR26SgGwMD4tO0SjVC2vZqFnfrNDHT2Ezc10dkN7O8XdYlHf38/ubm5VodhCla62XkUYy7sdAibzvck6O0nbmqisxvY3y/qEg8hPKiUaCTHx/DhVVnkpSWQkRgv6zUEQRAsJOpKpnu9XlM+1w6E2+3Sg9TSnHG4Jy8cqPZKizUVQUPBTqMaC6HzPQl6+4mbmujsBtb5Scn0eWhra6OqqsrqMExhOW52LcwVCqruQtH5ngS9/cRNTXR2A/v7RV3iMTk5aXUIprEYt7lGMU647D1VMhdpzlhuXZsdtiJeVqLzPQl6+4mbmujsBvb3i7rEIyUlxeoQTCMpOZkTZ9yzpkYyEuPJSlZzFMNPQgzcWGXt7hOz0PmeBL39xE1NdHYD+/tFXeKRn59vdQjLYqFtqkfahxiZOm9xhOHFv07jU2sySU5KtDocU1D9nrwcOvuJm5ro7Ab294u6xOP06dO2KyU7XzIxc9RC1SmRUFho+qS+vt52/RYu7HhPhhOd/cRNTXR2A/v7RV3iYTbz7QTxT3dAcFKhezKxEKouBhUEQRCWTtQkHv6EoG0qlZPv9AaNJMz83/kShMVca/d6FlYwcxTD355LqaVh96HD5aCzG+jtJ25qorMb2N9vSYnHvn37+OEPf4jL5eKKK65g79693HDDDfNef/DgQXbv3s17771HYWEh//2//3d27Nix5KBD5XDrIPve6KJ/dDpi3xmtmDWK4fNZe2CcmejsBnr7iZua6OwG9vcLOfF47rnnuOuuu9i3bx/XXXcdP/3pT9myZQv19fWUlpbOur61tZWPfvSjfO1rX+Ppp5/mtddeY+fOneTk5PDpT386LBILcbh1kAdebjX9e6KRjMRYbqrMjEhF0L6+PrKzs8P+uXZAZzfQ20/c1ERnN7C/X8iJx8MPP8z27du58847Adi7dy8vvfQSjz76KA899NCs63/yk59QWlrK3r17AVizZg1vvfUW//zP/2x64uH1Gex7o8vU74gG/KMY6wrSAlNMsiZDEARBWAohJR5TU1McPXqUe++9N+j1zZs3c+TIkTnf8/rrr7N58+ag1z784Q/z2GOPMT09TXx8/Jzvc7vdxMTEBP7udDpxOp2hhHthkadMr4SEXQtzrVq1ytLvNxOd3UBvP3FTE53dwP5+ISUe/f39eL1e8vLygl7Py8ujp6dnzvf09PTMeb3H46G/v5+CgoI531dXV8fY2Fjg79u2bWPXrl0UFBRw+vTpwOcYhkFfXx9wobG7uroYHx8nMTGRvum0UPSiiqQYqFvp4C/LshiemIbJUTISYvjIhtV0d3UyMTRIeVISRdlFNJ5sAAj0Y29vLwDV1dV0d3cH2ru0tJSmpiYAcnJyiIuLw+VyAVBZWUlvby+jo6MkJCRQWVnJyZMnAVi5ciVOp5MzZ84AUFFRQX9/P263m7i4OGpqajh+/DgpKSlkZWWRnJxMV9eFkayysjLOnz/P8PAwsbGxrF69moaGBgzDYMWKFaSnp9PR0QFASUkJbrebwcFBHA4Ha9asobGxEa/XS3p6OllZWbS1tQFQVFTE+Pg4AwMDwIWRuubmZqanp0lNTSUnJ4fW1gtTeIWFhUxNTdHf3w/A6tWraWtrY3JykpSUFPLz8wP3bH5+Pj6fL+iefffdd3E6nSQlJVFcXMypU6cAyM3NxeFwBNq7qqoKl8vF2NgYTqeTsrKyBdu7r6+PkZGROds7MTGR7u5uAMrLyxkYGGB4eDjQ3vX19QBkZmaSmppKZ2cnAKWlpQwNDTE0NERMTAy1tbWcPHkSn89HRkYGGRkZQe09MjJCV1cXKSkprF27lqamJjwez5ztPTExwblz5wCora2lpaWFqakpUlNTyc3NpaWlBYCCggI8Hg9nz54FoKamhvb2diYnJ0lOTg7pGVFSUhLU3jExMYFnWVVVFT09PYyOjuJ0OikvL6exsRGA7OxsEhISOHXqFCkpKVRUVHD27FlGRkaIj4+nurqahoYLvzdZWVkkJSXN2d5z3bNpaWlB7T08PDzvPZuZmUl7ezsAxcXFjI2NBe7Zme2dlpZGdnZ20D07OTk5Z3unpKSQl5fHO++8Q0pKypzt3dHRwcTEBElJSRQVFdHc3Gz5M8J/zy7mGXHixAmSk5OVeUZ0dnYG2nsxz4izZ8+SlZUV8WfEYo9+C+mQuDNnzlBUVMSRI0fYuHFj4PUHH3yQp556KhD0TGpqati2bRt79uwJvPbaa69x/fXX43K5Zq2+9R8SV1lZuewRjxNn3Nyzvzmk9+jAXFMj/t03dhrFCAWd63jo7AZ6+4mbmujsBtb5mXJIXHZ2NrGxsbNGN/r6+maNavjJz8+f8/q4uDhWrlw573elpaUt+3S9uvxUslPitZxumWub6pT7PGvKi5RLKhZDUlKS1SGYhs5uoLefuKmJzm5gf7+QEo+EhAQ2bNjAgQMH+OQnPxl4/cCBA9x6661zvmfjxo385je/CXrt97//PVdfffW86zvCRWyMg53XFCuzq2WuZGKu2iHzLeycns4yvU2tori42OoQTENnN9DbT9zURGc3sL9fyLtadu/ezZe+9CWuvvpqNm7cyP/6X/+Ljo6OQF2OPXv20N3dzZNPPgnAjh07eOSRR9i9ezdf+9rXeP3113nsscd49tlnw2syD9dXrOB7H6yIeB2PhaY7Lk0qwrVL5NSpU9oOH4qbuujsJ25qorMb2N8v5MTjjjvu4Ny5czzwwAO4XC7q6urYv38/ZWVlALhcrsBCHbiwCGj//v3cfffd/PjHP6awsJB//dd/jUgNDz/XV6xgY1kG7/aMcLyplbSVeaZULl1qVU5BEARBiBZCWlwaCfyLS9etW7fsNR5z0d/fb+vCKstB3NREZzfQ20/c1ERnN7DOb7GLS2Pm/YmmOBz6jkKIm5ro7AZ6+4mbmujsBvb3i7rEw7/fWUfETU10dgO9/cRNTXR2A/v7RV3iIQiCIAiCdURd4lFVVWV1CKYhbmqisxvo7SduaqKzG9jfL+oSD3+5WB0RNzXR2Q309hM3NdHZDezvF1WJx+TkJP/2b//G5OSk1aGEHXFTE53dQG8/cVMTnd1ADb+oSzx+/vOf27pDloq4qYnObqC3n7ipic5uoIZfVCUegiAIgiBYiyQegiAIgiBEjJBLppuNv5Cq1+sN+2f7fD6Sk5Px+XymfL6ViJua6OwGevuJm5ro7AbW+vm/73IF0W1XMn1iYoLXXnvN6jAEQRAEQVgC1113HYmJifP+3HaJh8/nY2pqitjYWNuXfRUEQRAE4QKGYeD1eklISCAmZv6VHLZLPARBEARB0BdZXCoIgiAIQsSQxEMQBEEQhIghiYcgCIIgCBEjahKPffv2UVFRQWJiIhs2bODVV1+1OqSQeeihh/jLv/xL0tLSyM3N5bbbbqOxsTHomq985Ss4HI6gP9dcc41FES+e+++/f1bc+fn5gZ8bhsH9999PYWEhSUlJ3Hjjjbz33nsWRhwa5eXls/wcDgff/OY3AbX67dChQ9xyyy0UFhbicDh48cUXg36+mL6anJxk165dZGdnk5KSwic+8Qm6uroiaDE3C7lNT0/zt3/7t1x55ZWkpKRQWFjI1q1bOXPmTNBn3HjjjbP68rOf/WyETWZzuX5bzD1o136Dy/vN9fvncDj44Q9/GLjGjn23mOe+ar9zUZF4PPfcc9x1113cd999HDt2jBtuuIEtW7bQ0dFhdWghcfDgQb75zW/yxhtvcODAATweD5s3b2Z0dDTouo985CO4XK7An/3791sUcWhcccUVQXG/8847gZ/90z/9Ew8//DCPPPIIb775Jvn5+dx888243W4LI148b775ZpDbgQMHAPjMZz4TuEaVfhsdHWXdunU88sgjc/58MX1111138cILL/CLX/yCw4cPMzIywsc//nHL6yos5DY2Nsbbb7/N3//93/P222/z/PPP09TUxCc+8YlZ137ta18L6suf/vSnkQh/QS7Xb3D5e9Cu/QaX95vp5XK5ePzxx3E4HHz6058Ous5ufbeY575yv3NGFPBXf/VXxo4dO4Jeq62tNe69916LIgoPfX19BmAcPHgw8NqXv/xl49Zbb7UuqCXy/e9/31i3bt2cP/P5fEZ+fr7xj//4j4HXJiYmjIyMDOMnP/lJhCIML9/5zneMqqoqw+fzGYahbr8BxgsvvBD4+2L6anBw0IiPjzd+8YtfBK7p7u42YmJijN/97ncRi/1yXOo2F//v//0/AzDa29sDr23atMn4zne+Y25wy2Qut8vdg6r0m2Esru9uvfVW46abbgp6TYW+u/S5r+LvnPYjHlNTUxw9epTNmzcHvb5582aOHDliUVThYWhoCICsrKyg1//0pz+Rm5tLTU0NX/va1+jr67MivJA5deoUhYWFVFRU8NnPfpaWlhYAWltb6enpCepDp9PJpk2blOzDqakpnn76ab761a8G1apRtd9mspi+Onr0KNPT00HXFBYWUldXp1x/Dg0N4XA4WLFiRdDrzzzzDNnZ2VxxxRV897vfVWZkbqF7UKd+6+3t5be//S3bt2+f9TO7992lz30Vf+dsVzI93PT39+P1esnLywt6PS8vj56eHouiWj6GYbB7926uv/566urqAq9v2bKFz3zmM5SVldHa2srf//3fc9NNN3H06FGcTqeFES/MX//1X/Pkk09SU1NDb28vP/jBD7j22mt57733Av00Vx+2t7dbEe6yePHFFxkcHOQrX/lK4DVV++1SFtNXPT09JCQkkJmZOesalX4nJyYmuPfee/n85z9Penp64PUvfOELVFRUkJ+fz7vvvsuePXs4ceJEYHrNrlzuHtSl3wD+/d//nbS0ND71qU8FvW73vpvrua/i75z2iYefS6ugGoahdGXUb33rW/z5z3/m8OHDQa/fcccdgf+uq6vj6quvpqysjN/+9rezfsnsxJYtWwL/feWVV7Jx40aqqqr493//98ACN1368LHHHmPLli0UFhYGXlO13+ZjKX2lUn9OT0/z2c9+Fp/Px759+4J+9rWvfS3w33V1daxatYqrr76at99+m6uuuirSoS6apd6DKvWbn8cff5wvfOELs8p6273v5nvug1q/c9pPtWRnZxMbGzsrq+vr65uVIarCrl27+M///E9eeeUViouLF7y2oKCAsrIyTp06FaHowkNKSgpXXnklp06dCuxu0aEP29vb+cMf/sCdd9654HWq9tti+io/P5+pqSnOnz8/7zV2Znp6mttvv53W1lYOHDgQNNoxF1dddRXx8fHK9eWl96Dq/ebn1VdfpbGx8bK/g2Cvvpvvua/i75z2iUdCQgIbNmyYNVR24MABrr32WouiWhqGYfCtb32L559/nj/+8Y9UVFRc9j3nzp2js7OTgoKCCEQYPiYnJ2loaKCgoCAw9DmzD6empjh48KByffjzn/+c3NxcPvaxjy14nar9tpi+2rBhA/Hx8UHXuFwu3n33Xdv3pz/pOHXqFH/4wx9YuXLlZd/z3nvvMT09rVxfXnoPqtxvM3nsscfYsGED69atu+y1dui7yz33lfydi/hyVgv4xS9+YcTHxxuPPfaYUV9fb9x1111GSkqK0dbWZnVoIfGNb3zDyMjIMP70pz8ZLpcr8GdsbMwwDMNwu93G3/zN3xhHjhwxWltbjVdeecXYuHGjUVRUZAwPD1sc/cL8zd/8jfGnP/3JaGlpMd544w3j4x//uJGWlhboo3/8x380MjIyjOeff9545513jM997nNGQUGB7b1m4vV6jdLSUuNv//Zvg15Xrd/cbrdx7Ngx49ixYwZgPPzww8axY8cCOzsW01c7duwwiouLjT/84Q/G22+/bdx0003GunXrDI/HY5WWYRgLu01PTxuf+MQnjOLiYuP48eNBv4OTk5OGYRhGc3Oz8T/+x/8w3nzzTaO1tdX47W9/a9TW1hrr16+3tdti70G79pthXP6+NAzDGBoaMpKTk41HH3101vvt2neXe+4bhnq/c1GReBiGYfz4xz82ysrKjISEBOOqq64K2oKqCsCcf37+858bhmEYY2NjxubNm42cnBwjPj7eKC0tNb785S8bHR0d1ga+CO644w6joKDAiI+PNwoLC41PfepTxnvvvRf4uc/nM77//e8b+fn5htPpNN7//vcb77zzjoURh85LL71kAEZjY2PQ66r12yuvvDLnffjlL3/ZMIzF9dX4+LjxrW99y8jKyjKSkpKMj3/847bwXcittbV13t/BV155xTAMw+jo6DDe//73G1lZWUZCQoJRVVVlfPvb3zbOnTtnrZixsNti70G79pthXP6+NAzD+OlPf2okJSUZg4ODs95v17673HPfMNT7nZPTaQVBEARBiBjar/EQBEEQBME+SOIhCIIgCELEkMRDEARBEISIIYmHIAiCIAgRQxIPQRAEQRAihiQegiAIgiBEDEk8BEEQBEGIGJJ4CIIgCIIQMSTxEARBEAQhYkjiIQiCIAhCxJDEQxAEQRCEiCGJhyAIgiAIEeP/B7Gf0S8ivc1zAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(scipy.linalg.eigh(A.A[1:-1,1:-1], inv(Linv[1:-1,1:-1]), eigvals_only=True), 'o-');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# 2D Multigrid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coarse dofs: 841\n",
      "Fine dofs: 2916\n"
     ]
    }
   ],
   "source": [
    "p = 4\n",
    "nspans_c = 25\n",
    "\n",
    "kv_c = bspline.make_knots(p, 0.0, 1.0, nspans_c)\n",
    "kv = kv_c.refine()\n",
    "h = 1.0 / kv.numspans\n",
    "\n",
    "M = assemble.mass(kv)\n",
    "K = assemble.stiffness(kv)\n",
    "\n",
    "M2 = spkron(M, M).tocsr()\n",
    "K2 = spkron(K, M) + spkron(M, K)\n",
    "\n",
    "A2 = 1.0*M2 + K2\n",
    "\n",
    "P = bspline.prolongation(kv_c, kv)\n",
    "P2 = spkron(P, P)\n",
    "\n",
    "print('Coarse dofs: %d\\nFine dofs: %d' % (P2.shape[1], A2.shape[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9 iterations\n"
     ]
    }
   ],
   "source": [
    "# construct prolongator for S-tilde and its orthogonal complement\n",
    "(P_tilde, P_compl) = Stilde_basis(kv)\n",
    "\n",
    "# basis for 2D S-tilde\n",
    "P2_tilde = spkron(P_tilde, P_tilde).tocsr()\n",
    "\n",
    "M2_tilde = (P2_tilde.T).dot(M2).dot(P2_tilde)\n",
    "\n",
    "# compute basis for M-orthogonal complement of S-tilde\n",
    "P_tilde_orth = solve(M.A, P_compl)  # M^-1 * P_compl\n",
    "\n",
    "subsp = (\n",
    "    P2_tilde,\n",
    "    kron(P_tilde, P_tilde_orth),\n",
    "    kron(P_tilde_orth, P_tilde),\n",
    "    kron(P_tilde_orth, P_tilde_orth)\n",
    ")\n",
    "\n",
    "sigma = 1 / 0.16 * h**-2\n",
    "\n",
    "B = ((1 + sigma) * M + K).A\n",
    "B_tilde_orth = P_tilde_orth.T.dot(B).dot(P_tilde_orth)\n",
    "M_sm_tilde = P_tilde.T.dot(M.A).dot(P_tilde)\n",
    "\n",
    "B_subsp = [\n",
    "    (1 + 2*sigma) * M2_tilde,\n",
    "    spkron(M_sm_tilde, B_tilde_orth),\n",
    "    spkron(B_tilde_orth, M_sm_tilde),\n",
    "    subsp[3].T.dot(A2.A).dot(subsp[3])\n",
    "]\n",
    "\n",
    "Binv_subsp = [make_solver(B) for B in B_subsp]\n",
    "\n",
    "smoother1 = OperatorSmoother(SubspaceOperator(subsp, Binv_subsp))\n",
    "smoother2 = GaussSeidelSmoother()\n",
    "smoother3 = SequentialSmoother((smoother1, smoother2))\n",
    "\n",
    "twogrid(A2, random_rhs(A2), P2, smoother3);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "34 iterations\n"
     ]
    }
   ],
   "source": [
    "# Dirichlet test\n",
    "\n",
    "# modify subspace prolongators to omit Dirichlet dofs\n",
    "subsp = (\n",
    "    spkron(P_tilde[1:-1], P_tilde[1:-1]).tocsr(),\n",
    "    kron(P_tilde[1:-1], P_tilde_orth[1:-1]),\n",
    "    kron(P_tilde_orth[1:-1], P_tilde[1:-1]),\n",
    "    kron(P_tilde_orth[1:-1], P_tilde_orth[1:-1])\n",
    ")\n",
    "\n",
    "# keep smoothers the same as before\n",
    "\n",
    "# compute Dirichlet stiffness matrix\n",
    "M2_D = spkron(M[1:-1,1:-1], M[1:-1,1:-1]).tocsr()\n",
    "K2_D = spkron(K[1:-1,1:-1], M[1:-1,1:-1]) + spkron(M[1:-1,1:-1], K[1:-1,1:-1])\n",
    "A2_D = K2_D #+ M2_D\n",
    "\n",
    "smoother = OperatorSmoother(SubspaceOperator(subsp, Binv_subsp))\n",
    "twogrid(A2_D, random_rhs(A2_D), spkron(P[1:-1], P[1:-1]), smoother);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# 3D Multigrid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def spkron3(A, B, C):\n",
    "    return spkron(spkron(A, B, format='csr'), C, format='csr')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coarse dofs: 3375\n",
      "Fine dofs: 15625\n"
     ]
    }
   ],
   "source": [
    "p = 5\n",
    "nspans_c = 10\n",
    "\n",
    "kv_c = bspline.make_knots(p, 0.0, 1.0, nspans_c)\n",
    "kv = kv_c.refine()\n",
    "h = 1.0 / kv.numspans\n",
    "\n",
    "M = assemble.mass(kv)\n",
    "K = assemble.stiffness(kv)\n",
    "\n",
    "M3 = spkron3(M, M, M)\n",
    "K3 = spkron3(K, M, M) + spkron3(M, K, M) + spkron3(M, M, K)\n",
    "\n",
    "A3 = 1.0*M3 + K3\n",
    "\n",
    "P = bspline.prolongation(kv_c, kv)\n",
    "P3 = spkron3(P, P, P)\n",
    "\n",
    "print('Coarse dofs: %d\\nFine dofs: %d' % (P3.shape[1], A3.shape[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# construct prolongator for S-tilde and its orthogonal complement\n",
    "(P_tilde, P_compl) = Stilde_basis(kv)\n",
    "\n",
    "# compute basis for M-orthogonal complement of S-tilde\n",
    "P_tilde_orth = solve(M.A, P_compl)  # M^-1 * P_compl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "sigma = 1 / 0.19 * h**-2\n",
    "\n",
    "# 1D smoother components\n",
    "B = ((1 + 2*sigma) * M + K).A\n",
    "B_tilde_orth = P_tilde_orth.T.dot(B).dot(P_tilde_orth)\n",
    "M_tilde = scipy.sparse.csr_matrix(P_tilde.T.dot(M.A).dot(P_tilde))\n",
    "\n",
    "subsp, Binv_subsp = 2**3 * [None], 2**3 * [None]\n",
    "\n",
    "for i in range(2**3):\n",
    "    bits = [bool(i & (1<<k)) for k in range(3)]\n",
    "    num_compl = sum(bits)  # number of components which are the orthogonal complement\n",
    "\n",
    "    prolongators = [P_tilde_orth if bits[j] else P_tilde for j in range(3)]\n",
    "    P_i = subsp[i] = spkron3(*prolongators)\n",
    "\n",
    "    if num_compl == 0:\n",
    "        # inner part (S-tilde): use mass smoother\n",
    "        M_tilde_inv = make_solver(M_tilde, symmetric=True)  # 1D solver for M-tilde\n",
    "        Binv_subsp[i] = 1/(1 + 3*sigma) * KroneckerOperator(*(3*[M_tilde_inv]))\n",
    "    elif num_compl == 1:\n",
    "        # faces: use Kronecker product smoother\n",
    "        smoothers = [B_tilde_orth if bits[j] else M_tilde for j in range(3)]\n",
    "        Binv_subsp[i] = KroneckerOperator(*[make_solver(X, symmetric=True) for X in smoothers])\n",
    "    else:\n",
    "        # edges and corners (small): use direct solver\n",
    "        fullness = P_i.nnz / prod(P_i.shape)\n",
    "        if fullness > 0.4:\n",
    "            P_i = P_i.A  # use dense matrix\n",
    "            B = P_i.T.dot(A3.dot(P_i))\n",
    "        else:\n",
    "            B = (P_i.T.dot(A3.dot(P_i))).A\n",
    "        Binv_subsp[i] = make_solver(B, symmetric=True)\n",
    "\n",
    "smoother1 = OperatorSmoother(SubspaceOperator(subsp, Binv_subsp))\n",
    "smoother2 = GaussSeidelSmoother()\n",
    "smoother3 = SequentialSmoother((smoother1, smoother2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "33 iterations\n"
     ]
    }
   ],
   "source": [
    "twogrid(A3, random_rhs(A3), P3, smoother1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "29 iterations\n"
     ]
    }
   ],
   "source": [
    "# Dirichlet test\n",
    "\n",
    "# compute new subspace prolongators which omit Dirichlet dofs\n",
    "for i in range(2**3):\n",
    "    bits = [bool(i & (1<<k)) for k in range(3)]\n",
    "    prolongators = [P_tilde_orth[1:-1] if bits[j] else P_tilde[1:-1] for j in range(3)]\n",
    "    subsp[i] = spkron3(*prolongators)\n",
    "\n",
    "P3_D = spkron3(P[1:-1], P[1:-1], P[1:-1])\n",
    "\n",
    "M_D = M[1:-1, 1:-1]\n",
    "K_D = K[1:-1, 1:-1]\n",
    "\n",
    "K3_D = spkron3(K_D, M_D, M_D) + spkron3(M_D, K_D, M_D) + spkron3(M_D, M_D, K_D)\n",
    "\n",
    "smoother = OperatorSmoother(SubspaceOperator(subsp, Binv_subsp))\n",
    "twogrid(K3_D, random_rhs(K3_D), P3_D, smoother);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}