{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 12: Iterative methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Syllabus\n",
    "**Week 1:** Matrices, vectors, matrix/vector norms, scalar products & unitary matrices  \n",
    "**Week 2:** TAs-week (Strassen, FFT, a bit of SVD)  \n",
    "**Week 3:** Matrix ranks, singular value decomposition, linear systems, eigenvalues  \n",
    "**Week 4:** Matrix decompositions: QR, LU, SVD + test + structured matrices start  \n",
    "**Week 5:** Iterative methods, preconditioners, matrix functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Recap of the previous lecture\n",
    "- Test\n",
    "- Short intro into sparse matrices and their storage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Today lecture\n",
    "\n",
    "Today we will talk about **iterative methods**, where they arise, how we store them, how we operate with them.\n",
    "\n",
    "\n",
    "- Concept of iterative methods\n",
    "- Richardson iteration, Chebyshev acceleration\n",
    "- Jacobi/Gauss-Seidel methods\n",
    "- Idea of Krylov methods\n",
    "- Conjugate gradient methods for symmetric positive definite matrices\n",
    "- Generalized minimal residual methods\n",
    "- The Zoo of iterative methods: BiCGStab, QMR, IDR, ...\n",
    "- The idea of preconditioners\n",
    "- Jacobi/Gauss-Seidel as preconditioner, successive overrelaxation\n",
    "- Incomplete ILU"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Direct methods\n",
    "\n",
    "Iterative methods are the main tools for solving large-scale numerical linear algebra problems. \n",
    "\n",
    "$$Ax = f,$$\n",
    "\n",
    "If the matrix $A$ is **dense**, the complexity is $\\mathcal{O}(N^3)$, and is feasible only up to $N \\sim 10^4 - 10^5$.\n",
    "\n",
    "If the matrix $A$ is **sparse**, the complexity is $\\mathcal{O}(N^{\\alpha})$, where $\\alpha = \\frac{3}{2}$ for 2D problems, and $\\alpha = \\frac{5}{3}$ for 3D problems (non-optimal).\n",
    "\n",
    "**But** the matrix-by-vector product is computed in $\\mathcal{O}(N)$ operations, **can we use it**?\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Matrix as a black box\n",
    "\n",
    "We have now an absolutely different view on a matrix: matrix is now a **linear operator**, that acts on a vector,  \n",
    "\n",
    "and this action can be computed in $\\mathcal{O}(N)$ operations.\n",
    "\n",
    "**This is the only information** we know about the matrix: the <font color='red'> matrix-by-vector product </font>\n",
    "\n",
    "Can we solve linear systems?\n",
    "\n",
    "Of course, we can multiply by the colums of the identity matrix, and recover the full matrix, but it is not what we need."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Richardson iteration\n",
    "The simplest idea is the **\"simple iteration method\"** or **Richardson iteration**.  \n",
    "\n",
    "\n",
    "  $$Ax = f,$$\n",
    "  $$\\tau  (Ax - f) = 0,$$\n",
    "   $$x + \\tau (Ax - f) = x,$$\n",
    "   $$x_{k+1} = x_k + \\tau (Ax_k - f),$$\n",
    "   \n",
    "   where $\\tau$ is the **iteration parameter**, which can be always chosen such that the method **converges**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Convergence of the Richardson method\n",
    "Let $x_*$ be the solution; introduce an error $e_k = x_{k} - x_*$, then  \n",
    "\n",
    "$$\n",
    "     e_{k+1} = (I + \\tau A) e_k,\n",
    "$$\n",
    "\n",
    "therefore if $\\Vert I + \\tau A \\Vert < 1$ the iteration converges. \n",
    "\n",
    "For symmetric positive definite case it is always possible to select $\\tau.$\n",
    "\n",
    "What about the non-symmetric case?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Optimal parameter choice\n",
    "The optimal choice for $\\tau$ for $A = A^* > 0$ is (prove it!)\n",
    "$$\n",
    "  \\tau = \\frac{2}{\\lambda_{\\min} + \\lambda_{\\max}}.\n",
    "$$\n",
    "\n",
    "where $\\lambda_{\\min}$ is the minimal eigenvalue, and $\\lambda_{\\max}$ is the maximal eigenvalue of the matrix $A$.\n",
    "\n",
    "So, to find optimal parameter, we need to know the **bounds of the spectra** of the matrix $A$,\n",
    "\n",
    "and we can compute it by using **power method**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Connection to ODEs\n",
    "\n",
    "The Richardson iteration has a deep connection to the Ordinary Differential Equations (ODE).\n",
    "\n",
    "\n",
    "Consider a time-dependent problem\n",
    "\n",
    "$$\\frac{dy}{dt} = A y - f, \\quad y(0) = y_0.$$\n",
    "\n",
    "Then $y(t) \\rightarrow A^{-1} f$ as $t \\rightarrow \\infty$, and the **Euler scheme** reads\n",
    "\n",
    "$$\\frac{(y_{k+1} - y_k)}{\\tau} = A y_k - f.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Convergence speed and condition number\n",
    "\n",
    "Even with the optimal parameter choice, the error at the next step satisfies\n",
    "\n",
    "$$e_{k+1} \\leq e_k q, \\rightarrow e_k \\leq c q^k,$$\n",
    "\n",
    "where \n",
    "\n",
    "$$\n",
    "   q = \\frac{\\lambda_{\\max} - \\lambda_{\\min}}{\\lambda_{\\max} + \\lambda_{\\min}} = \\frac{c - 1}{c+1},\n",
    "$$\n",
    "\n",
    "$$c = \\frac{\\lambda_{\\max}}{\\lambda_{\\min}} = \\mathrm{cond}(A)$$\n",
    "\n",
    "is the condition number of $A$.\n",
    "\n",
    "Let us do some demo..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((-0.0058683976325189921+0j), (-3.994131602367478+0j))"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAD9CAYAAABDaefJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXl8U1X6/99p0jRp0jRtU4pFdsUFcQMdXFEUFcQFnXFH\nGRf8iqID7v50qMuMzriNG26DjjLqLIg7uOAuLiOL4r6CgljomjRN26TJ+f1xcm/TUNoUmjYtz/v1\num22k3vy5ORzz33uc57HopRSCIKwrWAB5Dff/XSb3W3dsZOtwNLTHdhKevuPpzfav7fb3CCdtu/q\n9+4rNofeZfdEzO/AkjRT740/YkEQBCEu7Fk93QtBEASh6+iMqG8PnADslKa+CIIgCFtJKj71LKAM\nuC7hsdnAncAwYDSQH/+/M1AD3A281YX9FARBEFIgFZ/6jcC1SimWLl3K2LFjsdlsEeBCYC6bPzCc\nB/y9a7srCIIgbAYFHYv69sDaaDTKxIkTqays5KOPPiI7OzsGxJRSthdeeIFQKMT69et5++23OeCA\nA7jsssvIysqqAzzd9GEEQRC2dRR07H65CuDVV19l+fLlrF27luzsbIAPgf1ramo47rjjGDNmDLvs\nsgujRo2isbGR5uZm7Ha7A/ACDUAMyAOCQDgNHyYrvq9SYDugGCiI79MZ/58LuAB3/DFn/LYbcAA5\n8f92IJvWtlFANL6p+OdRCbdBHxAtSbct8b5lAdak1ygggrZPCGhC2yYQ3+rjjweA2vj/moTHQ2h7\nBuOPVwJ1Cf3pLtzog78XbWfDpoXxxzzo78Kwuwv9XeTEN+Nx474Vba/kCUaMFrsb/6HFxsmvNx5r\n6+zT+D4jQHN8C6Nt2oC2cR3atrXx25VoO/uBivj9SqAq/lh3hvVZ0DYsQts4H21HH9rejvjmRtu+\nIP46V/x543fhpGWsG2N0c6QaJWd8PzFa7BtB27cJbWM/LWPWGM818c0Y03W0tnMTPRM6aUHbbjv0\nmE7cjLGdn7C50HY1tMTOprqSeC3T0JEoLXYziNJii91IUTs7EvUpADfccANnnXUWubm5xuP7A6xc\nuZLCwkLmzp3LXnvthc3W6u3+hv5SUvHbPxj/H4l/iAj6g+eiDZId33LQBisA+tMiJM4U9rHFRCIR\nSzAYtDU0NNjq6uqor68nFApRU1OD3+8nEAhQU1NDIBCgoaGBhoYGwuEwjY2NNDU1EQ6HiUQiRKNR\nYjH9nWVlZVmys7PtDofD7nK58nNycrDb7Xg8HjweD7m5ubhcLvLz883N6/XicrnM53Jzc7FYWv2+\nokB1fPsl/r+WFjFqim8NtPz4jMFkRQ+8rPh/Q2iNH38uehCXoO3eDz2wu9T2SinTnvX19QSDQeI2\nz6qvr89qaGggEAhQW1trPhcMBgmHwzQ3NxONRolGo61uGxi2slgsFqvVasvOzrbZbDZsNht2u53c\n3FycTidutxu3243H4yE/P5+8vDx8Ph9erxev10txcTFOpzPR9mFaxL4u4XYgfj9Ei7A1xL+DRNs7\naBHYXFrEwR3/bwi3J27zvPhrt9rW1dXVVFVVEQgECAQC+P1+amtrqayspLa2lrq6OkKhEOFw2NLY\n2GiO7ebmZnMsGyf7WVlZWK1WS1ZWltVms1ntdnt2dnY2iZvT6SQvL4+CggKKiorMcVxQUEBhYSFu\ntxuXy0VeXp4xgUy0sXGANTZD9P3xrRqtH4aN6+N2h5YDkTHJstMixL74ZmhK4sSkJP79dCmRSISK\nigpKS0sTJ39t4Y73qwy4JpX3bs/9MhD4uampiZKSEi6//HJGjhzJkUceidOpf8f19fUMHz6cDRs2\nsNtuu/HEE0+w++67A5SjZ83RhQsXWlatWsX111+/2U5MnTqV6dOnk5+fT79+/SgoKMBut6fS/5RR\nStHQ0EBdXZ05eCsrKykvL6e6upq6ujpqa2vZsGEDGzZsIBgMUllZSUVFBXV1dSntw26343Q6cTqd\n2O12HA4HDocDu92OzWbDarVitVpRShGLxYhEIjQ2NhIKhWhqaqKpqYlAIEBjY2PKn8vlclFYWIjP\n5zN/KD6fj9LSUgoLC/F6veZ/h8NBXl4ebrcbh8OB2+3Gam1vctZ5wuEwtbW1VFVV4ff7CYVC1NfX\nU1lZid/vJxgMUlFRYYpyMBikqqqKqqoq88CY6iLn3Nxc8/Mk2zjxduKBz7B9NBolEonQ3NxMc3Mz\n4XCYUChEY2MjwWCQpqamdvftdDrx+Xz4fD6Ki4tNsXe73RQUFJjjuKCgAI/HYx4kunpsK6XMsVNX\nV0dVVRUbN25k/fr1plgbNt64cSOVlZXmAbOqqqrdz2m1WnG73eTm5ppjO3E8Z2VlkZXVokXGQTQW\ni5k2NWwciUQIh8M0NDQQDAZbHWw3R35+PsXFxebY9nq9eDweXC4XHo+H4uJiioqKzINEYWEhLpcL\nl8uF2+1OPih0KdFo1DwQhkIhGhoaqK+vx+/3mxON8vJyysvL8fv9+P1+KioqzO/A0JQ5c+a0q40D\nBw7kiy++IC8v71s6jjzcxP2SfEp1MsDixYsJBAJs2LCBefPmEY1Geeutt3jggQf4y1/+wpo1a6it\nreWpp57i3HPP5X//+x/Ar+gjoWXVqlUdGmj+/PnMnz+/1WM2m80cUA6Hg5ycHLKzszf58VosFvOH\nagwmYwAZM+qmpqaUBpLb7aZfv370798fr9fL8OHDzdsej8ecxRkDx+v1kp+fb/5Yc3Jy2n3/srKy\nNm8nE4lETDE0BkTi7NX4YRgzWeOHWlVVxWeffcaGDRuoqanp0O6A+SO12+3k5OTgcDjMWZXxnCH8\nSimUUuaP15i5GQPasHVHKKVaff6PPvqIXXfd1bRnfn6++QPNy8sjLy/PPDsxZnr5+fldfkBKpLm5\n2Tz4BwIBKioq8Pv9pp0rKipMoayqqmL16tVUVVVRV1dHJBJp970NgXS73axbt66VLZYsWUJWVpZ5\nIDIOQMa4Nr7/cDhMOBwmGAzS3Nzc7v6MmXC/fv3w+XwMHToUl8tFUVERAwYMwOfzkZeXZ56ZeL1e\nfD4fbrc7+UywSzDOxqqrq81xU11dzVFHHdXKFtXV1WzcuJHq6moqKir47rvvzINUKBTqcD+5ubm4\n3W5zkpWTk4NxRmxoh7ElaohxYIpGo4TDYZqammhsbDS3hoaGDm1uUFRUZGqEz+djhx12wOfzUVhY\nSL9+/SgvL2+3/bBhw7j99tvb1YtkEmfqyd/e/wNuevzxxzn33HNpampi3bp1DBs2jOXLl7NixQqm\nTZtmvvjhhx/m0Ucf5f333wdYARwG1Fgslg6PRtOmTeO0006jpqZmk1mcMXtqamoyxdo4tTZO/0Cf\nWmdlZWGz2UzxN07tcnJyzJmSx+MxB7DP56OkpITi4mJcLley+6jLSZ4xppOmpiZqa2vNU2m/32/O\nQuvq6szb9fX1povIGLyRSMQ8iwiHw6adjf4bM+GcnBzzzMSwtcfjwev1mrNU44dl3Ddm1d1pi+7G\nmDFXV1ebZx+G3RPddPX19cybN6+VLQ499NBNxrVh7+zsbPPAZhyIjXFtjHFDLEpLS7ttXHclqY6L\nSCRinuHV1dVRXV1t6oYx2THuGwdB44wmHA4TjUY3eX9DQ6xWq/nfOAgYZ90Oh8M8YzEOgsbE0+12\nk5+fb7rviouLcTja99yUlZW1q40JtvgGHTLe7sugfVHfCfi6vr6ePfbYA4/Hw88//8yf/vQnzj//\nfACefPJJZsyYQX5+PhUVFbz44ouMHz8edCjjeYDqyz/eziK2aEFs0YLYogWxRQtJtugSUQf4L/Db\nQCDAG2+8waGHHkp+fn6rF/j9fj766COGDx/O8OHDATYAewPrEVFvhdiiBbFFC2KLFsQWLaRL1N3A\nP4ATEx77HJiOviJ7F/qCqsFX8dd+ZexEvqQWxBYtiC1aEFu0ILZoIckW3wM7dtAkJVE3GAmMAL5D\ni7q5X2AHYBd0SNG7tI6zFFFPQGzRgtiiBbFFC2KLFtIt6luKiHoCYosWxBYtiC1aEFu0kGSLlEMa\nJfWuIAhC5pNy/G66Y51OnDNnztO33HKLsShpm2bOnDk93YWMQWzRgtiiBbFFC3PmzCFBO1PW6nS7\nX04BniopKeH444/nwQcf7LCBIAiCoEnQzrXAoA5e3i3uFyvoFXThcDryeAmCIPRdErQzZa1Ot6gr\n0HkyGhoa0rwrQRCEvkWCdqacMKhbRN3lclFfX5/mXQmCIPQtErQz5Wyo6Rb1XIDs7OwOkxwJgiAI\nrUnQzpRTTqZb1L2gMy6mmtVMEARB0CRoZ8rRL90yU7darSnlTxYEQRBaSNDOlOPU0y3qeaAromzr\nq8MEQRA6y5ZoZ7eIuiAIgtA9dItPXWbpgiAInWdLtNMQ9XSsJgVdzJVYLBY/jUjTXgRBEPoghnYS\nDw9PhXTP1D2gS09lZ2eLqAuCIHQCQzuBlJfkG6KeLrntB9DQ0IDT6UQyBQiCIKSOoZ1AY6pt0j1T\nLwIIh8PY7XZSKDQvCIIgxDG0ky2YqaeDLHQ5PKLRKDabjVAIZA2SIAhCahjaCaSsnOkU9TziF2CN\nU4hQCP73vzTuURAEoQ+R4H5JOSNiOkW9yLgRDAbjiWng7bfTuEdBEIQ+hKGdQMoZEdMp6l7jhnG0\naWiAf/4TiYIRBEFIgUybqReCdvSHw2E8Hg/hMHz5JSxalMa9CoIg9AEStROoS7VdOkU9H6CmpgaA\ngoIC+vXTT9x8cxr3KgiC0AdI1E6gJtV2aXe/BAIBAPLy8hg2DAoKYOlSeO+9NO5ZEAShl5OonUAg\nhSYWSK+oe0A7+gHcbjfZ2TBzpn7y+uvTuGdBEIReTqJ2kiEXSl3QcrSJ+4X4wx8gLw+WLIEPPkjj\n3gVBEHoxSdrpT7VdOkXdDZv4hSgoaJmtX3EFxGJp7IEgCEIvJUk7a1Ntl05Rd0JLxwoLC80nLr8c\n+vXTfvV589LYA0EQhF5KknZWp9ou7aK+ceNGAIqKzLVIeL1w11369mWXwfr1aeyFIAhCLyRJO6tS\nbZdOUc8FfbTJzs42ruCanHwyTJ4MgQDMmCELkgRBEBJJ0s7KVNulfaYeDAY3EXQAiwXmztUXTZ97\nDv71rzT2RBAEoZeRpJ3BVNulffFRRUVFK9dLIgMHwu2369sXXQQ//5zG3giCIPQikrQzLT714cDd\nwALgTFqXwCsBTgAOTXjPfIDa2tpWF0mTOfdcmDgRqqvht7+FxpRTwQuCIPRdkrSzy90vRwNfAjOB\nE4HHgAsBK3AL8CvwNPAG8A0winicenszddBumPnzYcgQ+PhjuPjiVLsuCILQd0nQzkYg5RJDqYj6\nAOA/gH3hwoXceuutxuMnAXcCV8ZiMctrr73GTz/9BLADsIr4itLKykr6GUlfNkNRETz9NDgc8PDD\ncP/9qXZfEAShb5KgnRvZwsLTm2v0RyB34cKFTJ8+nTFjxiS+fmZDQwPjx4/nqKOOYtKkSSxZssR4\nvr9SisrKynZn6gZ77w1//7u+fckl8O67qX4EQRCEvkWSdqbsT4eOZ+qFwJlKKe69915mzJjBoYce\najxXBfDiiy8yZMgQVq1axRtvvMHs2bOpqKgAsDY2NtLQ0MBtt92WUmdOP12nEYhE4JhjYOXKznwU\nQRCEvoHFYiEUClFcXAydiFGHjkV9DOBYuXIlH3zwARdccEHic1MAHnvsMS677DJGjhxJSUkJJSUl\nxOJr/6uqdF/KysqwWCyb3crKysw3ve02OPFE8PvhiCPgm28683EEQRB6Dx1p49KlS6ETKQKgY1E/\nBuCzzz4jFosxffp0zj//fH788UdisRjRaJTy8nIGDBgAwAsvvEBBQQElJSVAi6h3BqsVnnwSjjoK\nKivhyCNh7dpOv40gCEKvx+FwQCcyNALYOnh+HMC8efMYN24cM2fO5K677uLoo4/mkksuYePGjYwa\nNYprrrmG4uJivv76ax566CGjraqtrbVs/q03j90OCxbA4YfDhx/CIYfAm2/CoEFb8m6CIAi9k7io\npxz5Ah3P1H8G2H///QmHw2y33XZ4PB72228/TjrpJC677DLuvfdehg4dyuDBg3nqqafwer2go2Us\nlZU6tHLKlCkopTa7JbpfDFwuXfZuzBj48Uc46CD49tvOfDRBEITMpqysrE1NXLBgAQCXXnopdHKm\n3pGoPwYwa9YsmpubOeKIIxg2bBj3338/hYWF5Obm4nK5uOKKKzjnnHOwWq0AZcCpAHV1uqyekUu9\nsxQUwKuvwtixerXpuHHw1Vdb9FaCIAi9hiTt7FSC8o7cLwuAl0pKSo5+b9P6c8+iFyNtB5yODpB/\nFL34CAC/X+d1z8/P70yfWlFQoAtqHHssvPGGFvZXXoG99tritxQEQchokrQzuzNtO5qpK/QK0jJg\nPbpO3mvAsejol3XAx8AfgKtIEHTYtOrRluJywYsv6oumFRXaFRM/OxEEQehzJGlnR5PvVqSyorQJ\nuB69sjQfOAJ4IZU3DwaD5OTkYLN1qk9t4nTqbI5nnAH19fC730HL4lZBEIS+Q5J2OjrTNp1ZGgkG\ng0bR1C4hJwcef7wls+MVV8Cll0pJPEEQ+hZJ2mntTNt0irpqaGgwQnK6DIsFZs/WsezZ2XDHHbrY\nRrxIiCAIQq8nSTs3LUjRDmmdqYdCIXJzc9Py3qeeCi+/DIWFsHixDn385JO07EoQBKFbSdLOTrk7\n0irq4XAYu92etvcfP14L+X776VWnBx0Ezz6btt0JgiB0C0na2anwwbSKenNzc5dcJG2PgQP1atNT\nToFgEKZMgQsvhKZOrcESBEHIHJK0s7QzbXu9qIO+gPrEE9q/brfr2qcHHADr1qV914IgCF1OG6Ke\ncqx6OkXdEolEyM7uVNz8FpOVBbNmwXvvwdChsHw57LMPvP12t+xeEAShy0jSzizAm2rbtM7Uo9Go\nkTqg29hnH1i2TPvby8v1/5tukrBHQRB6D21oZ8oRJ2kVdaUUWVlp3UWbFBbqVAJXX63F/LrrdNGN\nLcgELAiC0O20oZ0pL8tPu+JaLFuUfXersdngz3/WmR4LC/X/3XbTYZCCIAiZTpJ2phyrnvaZek8z\ncSKsWAEHHqjdMRMnwsUXQ0NDT/dMEAShbdrQzsJU23a/b6QHGDwY3npLz9xtNrjnHr1YSWqgCoLQ\nS0g5Vj3top4Js3XQZfKuvho++gh23hm+/BJ+8xtdE1UuogqCkGkkaWdGRL8oq9VKNBpN4y46z957\n63DHGTMgEoHLL9cpfX/9tad7JgiCoGlDOzNiph7LRFEHyM2F++6DF16A4mJdhGO33eC//+3pngmC\nILQp6hkR/RKz2Ww0NzencRdbx+TJOnfMkUdCdTWcdBKceSbE89MLgiD0CG1opyvVtukU9WimztQT\nKS3VWR7vu08X4pg/H/bcU69MFQRB6Ana0M6Uc5hvk+6XZCwW7WNfsULXPl29Gg4+GK69VvvdBUEQ\nupNMdb80Z7r7JZmdd4YPP9RRMgB/+pMudL1mTY92SxCEbYw2tLMo1bbpFPVGu91OOBxO4y66Hrtd\nx7O/+SYMGAAffKDdMc8809M9EwRhW6EN7UylUIaC9Ip60O12EwwG07iL9DFuHHz6KRx3HPj9cMIJ\ncNFFshJVEIT004Z2OlNtm05Rr83Ly6Ouri6Nu0gvRUV6hn7HHboe6n33wb77wtdf93TPBEHoy7Sh\nnTmptk37TD0UCvWKi6Wbw2LRedo//BBGjIDPP9cpBh5/vKd7JghCX6UN7cyI6Jdar1evbPX7/Wnc\nTfdgrEQ99VSor4ezzoIzzpCYdkEQup42tDMjZup+j0dH4fRmF0wibrcumzdvnl6V+sQTeta+bFlP\n90wQhL5EG9pp3+yLk0inqDc4ndq339CHri5aLHD22Tqmfffd4bvvYOxYKCuDXhS9KQhCBtOGdmaE\n+6Xe5dIrW+vr69O4m55hp520n/0Pf4BoFK6/XpfOk5h2QRC2lja0MyMKT4eMjoVCoTTupudwOuHO\nO+GNN2C77eDdd2HUKHjgAUnnKwjCltOGdqZc7Dmt0S95eboCU1/xqW+OQw+FVavgxBMhGIQLLtBJ\nwtau7emeCYLQG2lDO1PW6nSKep3brRdB9XVRB/D5YMEC+M9/dHy7kc734YchQ+qECILQS9ga7Uzr\nTL0v+9Q3x+9+B198Acceq8Mdp0/Xs/b163u6Z4Ig9Ba2RjvTGtJYUFAAQHV1dRp3k3mUlMCzz8JT\nT+lZ+2uvwR57wL/+JbN2QRA6Zmu0MxVRzwP6AccAVwLTaUnYbkGXWdoOOAG4GjgLHShf4/F4sFgs\nfWLxUWexWOCUU/QK1AkToLJSL1w64QSoqurp3gmCkMm0oZ0pTwfbE/XdgS+AALABeB64BXgQ+Am4\nBPgBqAXWA08Dfwb+AXwLFFgsFjwezzYp6gb9+8PLL2vfusejZ/C77grPPdfTPRMEIVNpQztTjqfb\nnKhbgKeAXYPBILW1tXzxxRfcf//9rFy5EnRu378BQ+vq6qitreWTTz7h7rvv5ttvvwUYBPw/0KcR\ntbW1W/bJ+ghZWXDuuTrr48EHw8aNcPzx8Pvf6wyQgiAIySRpZ8oJtCxKtfLyWuL/twfW1tfX4/V6\naW5upqSkhO23357c3Fxee+01cnJy+OmnnxgyZAgAAwYMoH///pSUlPDiiy9isVjWA6VjxoyhuLiY\nxYsXb/2n7APEYnDPPXDVVdDYCAMH6rQDEyb0dM8EQcgkkrSzAcjtoIkCsG3myZ0APv30U/Lz8/nm\nm28oLCzEYrG0etGKFSsYNmwYn376KS6XK/n5N4HfFRUV2UXQW8jKgksu0RExU6fqvDFHHAGnnw53\n3aUvrAqCICxbtoyxY8cad1OeqW/O/RIE2H777amqqmKnnXbi0ksv3aTgxcCBA/nxxx8ZPnw4119/\nPU1NTcZTG4H/A0Jut5uysjIsFstmt7Kysk581L7Bzjvrqkq33KJXpj7xhI5r//e/JUJGELYVOtLG\nH3/80XhpyiXkNifq/wMYNGgQ77//Prfddhvl5eWcccYZADQ0NNQvWbKEMWPG8Prrr/OXv/yFjz/+\nmIsvvtho/wH6wGAuQBI2xWaDK6+Ezz6DAw+E8nIdMXPYYVKIQxAEEkvaRVJtszlR39O4sd9++zFt\n2jROO+00fvrpJwC+/fZbl+FSGT9+PNOmTeOEE04wn6clo1jAyAssbJ7hw+Htt+Ghh7T75c03dV3U\nsjIpnycI2zIJWRpTP39XrSG+XamUUs8884w69thj1fTp01Vpaam6/fbbW7340UcfNZ8vKipS8+fP\nN556IP4+b954440KUOFwWAkdU1Wl1NlnK6WdMEoNGaLUs8/2dK8EQehukrTzZ9Wiz5vblFJKJc7U\nE69yfgkwYcIEzj77bNxuNzfddBOzZs1qdUCYMmUKU6dOJScnhwcffNBwz0SBu+MvCWxL+V+6gsJC\nHQ3z9ts64+OaNTr88bjj4Icferp3giB0F0nambJPPXGmnqj4FqXU3DYOHhGl1GlKKZtS6t9tPF+v\nlDoq4X0ef+SRRxSgvv/++605aG2TRCJK3XWXUm63nrXb7UpdcYVSfn9P90wQhHSTpJ3fqC2YqbfS\nemAGMAKYBdyPXk06EngSaAZOBvYALgceAG4AdgFeTnifSp/PB0BNTU3KBxpBY7PBxRfDN9/AmWdC\nOAx//asugP3oo5KzXRD6MknaudUhjQbfoVeOzkDndfk26flVwG3ABcAc4Oek581MjcnhkELqlJbC\nY4/BRx/B/vvDhg26pN6YMfDWWz3dO0EQ0kGSdnaZqG8tDbm5ehFUX61+1J3suy+89x7Mnw/bbw8r\nV+oCHVOm6FqpgiD0HZK0c6tzv3QVgfz8fEDcL12FxQJnnKFdMjfeCC5XS5KwSy7R2SAFQej9JGln\nRlQ+Aqgtiq97F1HvWnJz4dpr4dtvtSsmGoW774Zhw+CGG0CCjQShd5OknRlRoxQgZCw+2tYzNaaL\n0lIdAvnJJ3DUUVrM58zR4n7HHbJ4SRB6K0naaU+1XbpFvd5ut5OTk0MgEEjzrrZtdt8dFi/WF073\n20+7YS69VK9WnTsXWtLyCILQG0jSzuxU26Vb1GOgHf5yobR7GDcOli6FF1/UqQZ+/RUuvBB22AHu\nvBO2oXKxgtDrSdDOjPGpx0CH5mxLxad7GosFjj4ali+HBQt09sd162D2bBg6VGeGFG+YIGQ+CdqZ\nMTN1BeB0OhMT0wjdRFYWnHiirrj0/PM6JLKiAq6+WhfnuOwynRlSEITMJEE7M0bUARH1niYrC445\nBj78EF55Raf2DQbh9tthyBC9WnXFip7upSAIyWSiqFsBHA5HYgENoYewWHSVpSVLdMWlKVN06oH5\n82H0aO2Pf+YZHR4pCELPk6CdGRP94gSwWq00NzeneVdCZxg9GhYu1JkfZ88GjwfeeQdOOEFXZbr/\nfgmHFISeJkE7N1d6dBPSLeoO0B2LyvQvIxk6VLth1q2Dv/1Nu2O+/x5mzNCpCGbPhtWre7qXgrBt\nkqCdGbP4KB/YpGC1kHnk5ek0A999p+ukjhkD1dU6DHKHHeB3v9NuG8kMKQjdx5ZoZ7pFvQB0x5RU\nU+4V2Gxw0knwv/9pv/vUqWC16tDICRNg5Eh44AGQZQeCkH62RDvTLequNL+/kCYsFu13f/xxXX1p\nzhwdBvn113DBBS0hkV980dM9FQQhkXSLuhsgGo2SldUt0ZNCGigt1UWwf/gBnnpKx7tXV2tf/G67\nwZFH6qgZuRYuCF1LgnamPF3vFvdLOBzGbk85IkfIULKz4ZRTdLz7Bx/A+eeD0wmvvqqjZoYNg5tu\nkgVNgtBVJGhnyjVKu0XUGxsbcTgcad6V0F1YLDB2rPatr12rZ+wjRujb110HgwbBqafqHDRyKUUQ\ntpwE7WxMtU26Rd0LEIlEyM5OeUGU0IsoKtJhj199pWfsU6boxUv/+hcceCDssw888ojEvAvClpCg\nnZFU23RLSKOIet8nK0tHxyxcqOPar75aC/7y5XDOOTrm/aqr9GxeEITU6KSoW6BF1NMVSO4Ccb9s\nawwaBH9DIKCcAAAgAElEQVT+sxbwf/xDz9arq+Evf9GLnSZP1qmBZT2aILTP1rhf0uX5dAA0NTWR\nk5OTpl0ImYrTCWedpWPeP/hAX2S1WOCll3SCsR131GmA5cKqILRNgnamnDwr3e6XXJDoF0FfWH3q\nKfjlF7j1Vj1jN9w0gwbpYtoff9zTvRSEzCITo19MUZeZugDQr59etPTdd3rGfvzx2g3zxBM6/v2A\nA/TqVYl5F4RW2pkxM3VHNBqlsbERl0sWlwotWK0waZJetPTjj3D55eD1wvvv6zwzO+wAt90mFZqE\nbZck7Uw5MUc6RT0HsAaDQQDcbncadyX0ZgYPhr/+VV9Yve8+Leg//aSFfvvt4aKL4Ntve7qXgtC9\nJGlnXart0inqRaCv3gIS/SJ0iNutU/5+842OjjnsMF0o+777YKed9MXVt96SBU3CtkGSdmbE4qN+\n0HK0ycvLS+OuhL5EVpYunL1kCaxapePcc3K00B96qE409s9/QiTl5RiC0PtI0s5gqu3SKeoDAOrq\n9FmDuF+ELWHUKPj737Vr5vrrobgYVq7UKYGHDtVuG/G7C32RJO3MCPdLKYj7Regaiovhj3+En3/W\nIr/rrjo88sordUjkZZfp5wShr5CJ7pd8gNr4NCo/Pz+NuxK2FRwO7Y75/HNYtAjGj4e6Op1UbNgw\nHe++alVP91IQtp4k7Qyk2q4rRT0v6f08ANXV1QD4fL4u3JWwrWOxwMSJ8PrretHSqafqx594AvbY\nQ6ciePfdnu2jIGwNSdpZmWq7jkTdChwHXAzcDLwGfAX8G9gl/poTgC/RR5Ia4MZ4u1Yzda/Xm2qf\nBKFTjBkDTz6pi3jMnKnTE7z0Ehx8sL6w+vLLEjEj9D6StNOfarv2RD0LeAl4FrgLuAo4HNgZOAn4\nH/AY8DSwS319PUopD3AtsIR4LnXjaFNQUNCJjyMInWfwYLj7bh3jft11ejHTW2/pGf3o0fDCC1I4\nW+g9JGlndart2hP1I4Ej169fz7333ss999zD+eefz5577snzzz8PulTdmdFolMsvv5yioiKOPvpo\nKisrAQ4BDgZYv349BQUFkvtF6DaKi+GGG3RumZtvhu220xEzxx4Le+4J//63iLuQ+SRp54ZU27Un\n6mMBXnjhBa6++mqWLl3KwIEDmTZtGjabzXzR888/z6pVq7jrrrs455xzOOWUU4jqnKqDASorKykp\nKdmSzyQIW4XXq3O4//AD3HGHrrX62Wc6W+Ree+lCHpL+V8hUkrTz15QbqhZI2u5TSqmZM2eqSZMm\nqaVLl6poNKqSmThxolq+fLl5f6+99lLBYNC8P378+E3aCEJP0Nio1AMPKDVggFLay67UTjsp9cQT\nSrUxtAWhx9lvv/2MmzlqU41ua1PtzdS/BzjhhBN4++23OeCAA5g0aRI1NTWtXvTtt9+y9957A7B8\n+XIGDx7cKnlXTU0NZWVlWCyWzW5lZWUpH4QEYUvJydHFsr//Hh58EIYM0SkJTj9du2WefFKyQwrd\nS0faGHdn19NFWRqXAhxyyCH4/X4aGxtxOp3ce++9ADz55JM8++yzuN1u3nzzTV5//XX+8Ic/cPvt\ntxvt6wCjU4KQMTgcMH26ThL28MM6adhnn2lx32MP+M9/xC0jZAZOpxM6EfkC7Yv6LICGhgZCIZ31\n0e/3m8v9S0tLGTNmDA899BBXX301TzzxBIsWLWLYsGEAnxFP6m5cwRWETCM7G849V+d2f+ghnXbg\nyy/h5JNh5EiYP19m7kLPkpubCzpUPHUS3DfJvpn3lVLq/PPPVxaLRXk8HrXvvvsqv9/fkRtoiVLK\nq5QKxmIxBajrrrtu651LgpBmGhuVmjtXqSFDWnzuQ4cqdf/9+jlB6C6StPNjlZo/fROfevLyjP8A\nPPDAA9TV1bFo0SKWLl2Kx+MBWIxekHQRcDvwA7AMOAuYANQCVmOGLwUyhN5ATg5ccIF2yzzyiK6h\nunq1fmz4cLj3XmhMOQOHIGw5SdrZqVHXnvvlbuBy4FeXy8UBBxyAzWaLxR8/FrgHuA+4DNgB2Ad4\nnJaDg62+vh6QDI1C7yI7G37/e/jqKx3TPmqUTh42c6Z20dxyC/g75eUUhM6RpJ2dSjLdnqjHgNvQ\n2RaHAgeh0+leAqTiaTRn6nFnvyD0KqxWOOkk+OQTWLhQR8iUl7cUy77iCi32gtDVJGlnp67spJrQ\naw3wHlDeife2GEcbcb8IvZmsLJgyBVasgFdegUMOgUAAbr1Vz9zPO0+HSQpCV5GknWkR9S3CSEgj\neV+EvoDFAkccAW++CR99pAtkR6M6v/uIEXDccVJuT+gakrQzc0TdH3c8xi+uCkKfYd99dTz7V1/B\n2WdrP/zzz+uskPvuC48/LhdVhS0nSTstnWnbLaIuBTKEvsqIETBvnq66VFYGRUWwbBmcdRYMHKiz\nRf6aetYOQQA20c6MEXUloi5sK5SUwJw5LeX29twTKivhppt0SuBTToF33hHXjJAaSdrZqVGTVlGX\nC6XCtkZuri63t2KFFvEpU7Tf/d//hnHjYLfd4K67QBZaC+2RpJ2dEtC0ul+MsJz4UldB2GawWOCg\ng3Qo5Jo1cO210L+/TkPwhz/oNMCnn67L8UludyGZJO3s1EXJtIp6IBDA4XCQnZ2dzt0IQkYzcCDc\neKN2zSxYoCNompp0VsjDD9crV+fMga+/7umeCplCknZmjKirQCAgkS+CECc7G048Uce6//ijFvKB\nA/XtG26AXXaBsWN1WuB4RJuwjZKknZkj6uFwmJycnDTuQhB6J0OH6miZ1avh1Vd1WoK8PB3//n//\np101J56oZ/YNDT3dW6G7SdLOQjoRAZNW90tTU5OIuiC0g9UKEyboBGLl5Tq+/bDDIBzW/vjf/U4L\n/DnnwJIlkgp4WyFJO7MBR6pt0ynqWfX19XKRVBBSJDcXpk7V4r12ra6rOnq0TknwyCNa/Pv1g1NP\n1fVVA4Ge7rGQLtrQzpSX5adT1C2NjY04HCkfYARBiDNgAMyapRcyffmlXsS0005QU6MF/dRTtcAf\nfzw88QTEI+CEPkIb2jkg1bZpFfVIJCKRL4Kwleyyi76Q+vXXOtf77bfrcMlwGJ57Ds44Qwv8Mcfo\n8nwbN/Z0j4WtpQ3tLEq1bVp96rFYjKystO5CELYpdtwRZs/WC5vWrdOFO/bfH0IhePFFXXt1u+10\nDpr774cNG3q6x8KW0IZ25qXaNu2Ka7F0Km2BIAgpUloKF14IS5dqgX/4YZg0SV98festmDFDC/x+\n++k4+c8/7+keC50hSTtTvjiZVlHPysoiKmXZBSHtDBigi2i/9BJUVMBjj8HRR4PdDh9+CH/8o67g\ntPvu2pXz2WeShyaTaUM7M0LUY1arVURdELqZ/Hw480ztjqms1H73c86BwkIt5nPmaHHfeWddxemD\nDyRVQabRhnb6Um2bTlFvzs7OplkCawWhx3C74dhjdebIX3/VQn/OOeDz6Yuut9yiffKDB8Oll+po\nG5nB9zxtaGdGiHrEbrcTDofTuAtBEFLFbtcuGUPglyyBiy/WqQrWrdNx8fvso0Mny8q06As9Qxva\n6U21bTpFvSknJ4dGKf8iCBmHzaZXrt51l84i+f77MHOmDo387ju4/not7qNHw513ShRNd9OGdmbE\n4qMmh8Mhoi4IGU5Wlo6Quftu+OUXnXBs2jTweHRe+Nmz9YXYY47RuWiamnq6x32fNrQzI0S9Qdwv\ngtC7sNl0auBHH9Wz8wULtE/eYtH++N/9TodJXnQRrFzZ073tu7ShnRkR/VLvdDppkBRzgtArcTh0\npsjnntMz+Dvv1GX6amrgvvtg771hr73gnnukklNX04Z2ppwZMa0zdRF1Qegb9OunKzatXKm3mTOh\noAA++URfbDUqOb35poRHdgVtaGfKJe3SKep1Ho+HcDhMkzjhBKHPsOee2v++fr1OLnbEEToPzZNP\nwvjxMGKEDpX89dee7mnvpQ3tdKfQTEF6Rb26oED79mtqatK4G0EQegKHA04+uaWS03XXwfbbww8/\n6EVNAwfCccfByy/L7L2ztKGdGZFPPWCUYwpI4mdB6NMMGaLTD6xZo1MVnHCCvrj6/PMwcaJORHbn\nnRAM9nRPewdtaKc11bZpFfX8/HwA/H5/GncjCEKmYLXqpGJPP60XNP3pTzBokJ7Jz56tZ/IXXwyr\nVvV0TzObNrTTlmrbdIp6jXEKUStVdAVhm6OkBK65Rgv6s8/CgQeC36+jZfbYA8aNgxdeAEkPtSlt\naGdGiHq9UY6pXsqyCMI2i9WqfevvvqsjZy68UBfZfucdHQM/dKie0ZeX93RPM4c2tDMjCk/Xu1w6\nCkdEXRAE0JEz996rXTO33grDh+t6rNdeqy+snnoqfPRRT/ey59ka7exKUbcBuwND4/eDeXm6WEdQ\nro4IgpCAxwOXXaaThr36qq61qpQOkRw7Fg4+GBYuhEikp3vaM7ShnSnnzmxP1LOAu4Ag8AuwLn77\nZ+A4dDD8zcBS4EsgAnwK/AgsBqxOpxNAFiAJgtAmWVkwYQI884z2vV95JXi92lVz4ok6quaWW3Re\n+G2JNrSzS0T9dOBitHiXoqtZu4CBwD+BfwBXAfsDuwD89NNPxtXao4Ar3G4dLy8hjYIgdMSgQVrA\nf/oJ/vY3XXB7/Xod8z5ggF6x+s4720a+9za0M+VI//ZEfTzAVVddxXPPPQfAxIkT+VYnWXYDv43F\nYkyYMIGxY8cyZMgQhgwZwl577WWcMozMzs4mLy9Pol8EQUgZjwcuuQS++AIWLdIhkpGIXrE6bhyM\nGQP/+Af05QSwbWhnyo6o9kTdDvo0wIiVDIVC/Pzzz+YL1q5dy7vvvssFF1zAggULqKioYPny5cZR\nZiNAcXExd9xxR6c+kCAIgsWiFy699BKsXq0vphYX63TAv/+9ntlfdZWe2fdFAoEA5S0hQSmnu21P\n1D0AP/74o1nVesKECaxbt858QSQSobm5mVdeeYX169fj8/nM+Eq0+wafz0dZWRkWi2WzW1lZWar9\nFQRhG2TwYLjxRi3gjz6qM0RWVMBf/gLDhsGUKb0zHUFH2vhRSyhQyuclHV0opbKykt12242amhrW\nrFnDQw89RGVlJTU1NQwfPpwbb7yRWCzGlVdeyaWXXmq0/Rh4HYgaviFBEIStxenUBTyWLdPVmk4/\nXcfBP/usntXvsIMW/75SqSkhp3rKIYTtifovAOPGjePAAw9k0KBBPPbYY9hsNvx+PwcffDAWi4Wr\nr76af/3rXzzyyCM8//zzRtsC9NXa2pKSki34KIIgCJvHYtHVmv75T/j5Z/jzn3Wc++rV8Mc/tsS8\nv/FG75u9J5IQ0phyCKFFqVbXkhNXLR0EvBMKhfj888/Zcccdueeee7Db7Vx11VX4/X7Tp15SUsLj\njz+OzWZjwYIFAMuBMcAXs2bN2vXvf/87dXV1W/nxBEEQNk80qkX8vvt0+gFDzIcO1T74M8/Ubpze\nwqxZs0jQzhXA6A6adJh6913gb7m5uey7774UFBSw3377sf/++wM64UxOTg5ffPEFc+bMAWDevHlG\n2/nx/wGv10swGCQqCR4EQUgjVquOeX/2WR3zPmdO69n70KH6+SefhFCop3vbMUnambKAJs7UN5db\nYGd0LPq8hMeC8de3VY3jX8AZ8U68cscddxxx6aWXUltba2YeEwRB6A6iUXj9dR0C+cwzLWGQbrdO\nD3zGGbqwhzXlxLbdxx133EGCdn4I7NdBk5SLZHwNPALsCNwJ3AKMAAqBC9GrR38EVgGX0CLoANWF\nhYUAVFRUdObzCIIgbDVWq67M9OSTeiHTfffBb36j87o//rh+bvvtdam+jz7KrIVNSdrZnGq7zuR+\n+R6YDVwN/IqOm5wLTAKGA3sAd9P6NOGX/v37A7Chr1yOFgShV1JQADNmwIcf6pwzc+boaJnycrjr\nLp1zZocddLrgVat6XuCTtDMjsjQCVPh8PkCHRgqCIGQCO+4IZWVa3P/3P72CdbvttC/+5pt1vvdd\nd9Wv+fLLnhH4JO3MiMpHADWGH13yvwiCkGlYLLDPPjrXzNq18OabcP75UFQEX38N118PI0fCzjvr\nGfyyZd0n8Ena2SUJvbqCOqPWnoQ0CoKQyVitcMgh8MAD8OuveoXq2WdDYaGe0d98sz4ADB2qS/K9\n9VZ6qzYlaWfGzNT9Rl5gmakLgtBbyM6GI4+EefP06tQlS3TFpu2206kK7rkHDj1U3z/zTJ0HvqtL\nMSdpZ8aIetDpdGK1WkXUBUHoldhscNhhLRWbPvhA530fPlznn5k/X69eLS6Go46CuXP167aWJO3M\nSbVd2kXdYrHgdrul+pEgCL2erCwdJXPLLfDddzo98G236UpN0Si88oqe0Q8cqEv3XXONjrbZklQF\nSdppT7mPnd9VpwgD5OTk0NTUlOZdCYIgdB8Wi46QufRSePttHRr56KM6Y2RuLnz6qfbD77efnsVP\nnqwvyK5Zk/o+ErTTlmqbdIt6E+jK2KHesC5XEARhCyku1hkkFy6Eqipde/Wii3S+mepqnRd+1ix9\noXWnneD//g+efhraqyGUoJ0pz9RTVv8tJAy6MvaWVMUWBEHojTgcOs/MhAlw99364uq778KLL8Li\nxTqa5ttv4cEHtUtnn33g8MP1CtexY8Eel/AE7WwrJUubpHumbgGw2+2JeYEFQRC2GSwWXUB76lT4\n97/1LP6DD+Cmm7QvPitLpyj40590ub7rr29pm6CdGTNTt4GutxeJpFxiTxAEoc+Sna1n42PHwv/7\nfzoPzTvvwGuvaZfN4YcnvtbUzuxU3z/dom4HsNlsNDennI9GEARhm8Ht1sW1J03a9LkE7cyYC6Vu\ngKysLGK9ufyIIAhCD5CgnSlrdbpF3Qm6Y6qnU54JgiD0MrZEO7tlpm6xpJw1UhAEQYizJdqZblEv\nBGSWLgiCsAVsiXamW9S9ANFolKysdO9KEAShb5GgnRmTejcPdMesmVgEUBAEIYNJ0M6Uk/ymW9SL\nAJqbm7HZ0h09KQiC0LdI0M601CjdErwA4XAYuz3lBVGCIAgCrbQz5SX56Rb1XIBIJEJ2dsoLogRB\nEARaaWfKS/K7xafe1NRETk7KOd4FQRAEWmlnyrnLuyVOvbGxEYfDkeZdCYIg9C0StLMx1TbdMlMP\nhULk5uameVeCIAh9iwTtbEjh5RboJp96MBjE7XaneVeCIAh9iwTtrEu1Tdpzv4TDYZqbm0XUBUEQ\nOkGSdqZcZSjtM/WamhoAvF5vmnclCILQd0jSznaK3rUm7TP1jRs3AlBcXJzmXQmCIPQdkrRzY6rt\n0inqFsDp9/sBKCgoSOOuBEEQ+hZJ2pkRM3UnYG1s1JE4EqcuCIKQOknamXJIo0W15Hbs6qTnA4B1\nZWVl5gOJt7dFxBYtiC1aEFu0ILZoIckWs4E7U2im0inqY4CPE5O8b+t51cUWLYgtWhBbtCC2aCHJ\nFtOAx1JoprrK/eIGxgOjaTk49Oui9xYEQdjWqUz1hR2Juhd4GPgCeAZYGL/9D3RaXQvwB/SV2deB\nZcCnwK6Ap5OdFgRBENpmQ6ov7CjJ+a3AufHbuyY8visQA34AbgJYsWIFpaWl9O/ffxRa+K9PubuC\nIAhCe3TZ4qPjAU444QSampoIBAKcfPLJxnOnATc1Nzdz8sknM3r0aA466CDefPNN4/kLO9trQRAE\noU0Cqb6wo5m6E+CHH36gqamJ8vJy3n//faMaRw7AW2+9RUNDA++99x7Dhg1j0qRJvPTSS5SWlvqM\nN+ntFzwikQjBYJCGhgbq6uqor68nFApRU1OD3+8nEAhQU1NDIBCgoaGBhoYGwuEwjY2NNDU1EQ6H\nWbRoUav3HD9+PNnZ2TgcDlwuFzk5OdjtdjweDx6Ph9zcXFwuF/n5+ebm9XpxuVzmc7m5uVtUbTyT\nUUqZ9qyvrycYDJo2r6+vp6GhgUAgQG1trflcMBg0l1RHo1Gi0Wir2waGrSwWC1arlezsbGw2Gzab\nDbvdTm5uLk6nE7fbjdvtxuPxkJ+fT15eHj6fD6/Xi9frpbi4GKfT2ettr5SiurqaqqoqAoEAgUAA\nv99PbW0tlZWV1NbWUldXRygUMsezMbabm5uJxWLm+wBkZWVhtVrJysoybZqdnd1qczqd5OXlUVBQ\nQFFRkTmOCwoKOOCAA3rSHN1KJBKhoqKC0tLSzb5GKZU4xlKeqSeKuqJ1BEwW4IrFYmzYoN05I0aM\noF+/fgSDQXPZ/6OPPsoVV1xhfiGDBw824ysNysrKuP76zXtjpk6dyvTp08nPz6dfv34UFBR0eaUk\npZQpysbgrayspLy8nOrqaurq6qitrWXDhg1s2LCBYDBIZWUlFRUV1NWllkvHbrfjdDpxOp3Y7XYc\nDgcOh6PNzxIOh6mvr6exsZFQKERTU5N5NpRsv/ZwuVwUFhbi8/nMH4rP56O0tJTCwkK8Xq/53+Fw\nkJeXh9vtxuFw4Ha7u7x2bDgcpra2lqqqKvx+P6FQiPr6eiorK/H7/QSDQa699tpWbcaMGUNVVZV5\nYEx1EpCbm2t+Hrvdjs1mw2q1YrVaW902fhhKKZRSxGIxotEokUiE5uZmmpubCYfDhEIhGhsbCQaD\nNDW1n77a6XTi8/nw+XwUFxebYu92uykoKDDHcUFBAR6PxzxIdPXYVkqZY6euro6qqio2btzI+vXr\nTbEOBoPm45WVleYBs6qqqt3PabVacbvd5ObmmmPbGM82m42srKxWBeXD4TDRaJRYLGba1LBxJBIh\nHA7T0NBAMBhsdbBN/CyJGDY1xrbX68Xj8eByufB4PBQXF1NUVGQeJAoLC3G5XLhcLtxud1oL80Sj\nUfNAGAqFaGhooL6+Hr/fb040ysvLKS8vx+/34/f7qaioML8DQ1PmzJnTrjYmkLIotDdTtxqdz8nJ\nIS8vj1gsRn19PU888QQXXqi9K+vWrWPo0KEAfPDBB4RCIYYNG2a+yW9+85sOOzF//nzmz5/fumM2\nmzmgHA4HOTk5ZGdnb/LjtVgs5g/VGEzGADJm1E1NTZsdSIm43W769etH//798Xq9DB8+3Lzt8XjM\nWZwxcLxeL/n5+eaPtaMFVnPmzDFvtxeDG4lETDE0BkTi7NX4YRgzWeOHWlVVxWeffcaGDRvMvBEd\nYfxI7XY7OTk5OBwOc1ZlPGcIvyGKxgzYmLkZA9qwdUf89a9/JT8/39xvcXExu+66q2nP/Px88wea\nl5dHXl6eeXZizPTy8/PTWsy8ubnZPPgHAgEqKirw+/2mnSsqKkyhrKqqYvXq1VRVVVFXV0ck0n6R\nGkMg3W43TqeToqIi80B04IEHkpWVZR6IjAOQMa6N7z8cDhMOhwkGgzQ3t1++0pgJ9+vXD5/Px9Ch\nQ3G5XBQVFTFgwAB8Ph95eXnmmYnX68Xn8+F2u9NyNmKcjVVXV5vjprq6mqeeeorJkycTiURoampi\n1KhRbNy4kerqaioqKvjuu+/Mg1QoFOpwP7m5uaaNjfFtnBEb2mFsiRpiHJii0SjhcJimpiYaGxvN\nraGhoUObGxQVFZka4fP52GGHHfD5fBQWFtKvXz/Ky8vbbT9u3DgOOeQQ6MRC0cQ4dWg9U7cAUaWU\nZdSoUTidTj7//HMaGxs599xzOeyww1i1ahUbN27E7XYzdOhQXn31VR5++GHjlOJrYGfoeKY+bdo0\nTjvtNGpqaqioqDBPrY0vz3BjGGJtnFobp3+gT6mN0z5D/I1Tu5ycHHOm5PF4zAHs8/koKSmhuLgY\nl8vVVnHsGLriSBP6SBlG5zUOxTc/UB3fNgAVaN+XHwii02X60adORvtw/AuyoVMT56MjidzxrT/g\nQ0cPeYDi+Gu86Pz0DiAnvtkTtlYK19TURG1trXkq7ff7zVloXV2debu+vt50ERmDNxKJEIlEaGxs\nJBwOm3Y2fuCGAOXk5JhnJoatPR4PXq/XnKUaPyzjvjGrboMouriuYe+6BFvXxW1o/PcD6+N2r4zb\nPhS3cSRu40R7xxL2Y9jeHre/YU8XUJBgcy9QEn/Miw7RNb4Td7xtTrLdAXPGXF1dbZ59GHZPdNOB\nXjVYW1trzmIjkcgm49qwd3Z2tnlgMw6Ixrg2xrghFqWlpR2N61jcNsZ4DcVv1wJVwK/xrSpu52D8\ndYH4Fo5/X4nqZoxFW9yuLrQLNyf+34Ueyz5g+7hN3ehxXRi3rfF647vJpg1Bi0QiVFVVmQfS6upq\nUzeMyY5x3zgIGmc0xhlF8pmBoSGGC8lqtZoHAeOs2+FwmGcsxkHQmHi63W7y8/NN911xcXFbxYEM\nuzcCwbKysu3b08Y5c+ZQVlbmj9snttkXtqDaE3WAecDZy5YtY+3atYwePZpJkybx7LPP4vF4TL/Z\nrbfeitPpZNasWcZs9V103PrRwDT0j+NL4FH0D6QQGA6Uor9YQ6wM8XKiB0c2+keTFe9bZ6YNipbB\nG0X/2A1xrkcLRiV64K4FauK3V6PDMlP2YbWDD5iDXoi1Hh0RlA1cDOyIjhK6Efi5C/YF2n7bA7uh\n7bsd2taGIBniZfxYHOgfoGFjY+uMrRNtbIhyAy0iYCT3V8BXaPtbgbeBb9Bj5Ij48wuBB+K3ewsW\ntG13BAajD9CD0LY3Dg6GWOXSIlYOtixNh2FvRWubG/8D6MnFevS4rkKP83Vo+68ns+w7ADgTba+f\n0AtsaoELgGPR9n0GeBo9to0x3Z+WA4MHbfc8tH3ttEx8jIOChc5rCGhbGVsUbXvjYNZIy2SiCa0Z\nxkGyDm3zL4FfgPL450u0fS4wFzg53vYR4Ha0PhyMPpg+DLwH7An8HzAQWA7cQNux6wrVGpK2gUqp\n9YkvmD17tlq7dq1qh7lKKYdS6qE2nqtUSo1oYz99cdsr/nlT4cgM6G+6tvOUUrEU7aCUUm8ppfIy\noN/p2vKVUt8qpVRTU5NqbGxszxbzlFJZGdDndG1nKqUiSZ+5USkVbMMW7yilPBnQ567c/tTG5+wM\n+3sIgk8AAAirSURBVLTxnvpPApsbhOcopb5Iem1AKbW7UupgpdQjSov42HibC5VSqq6uTl111VVq\n4sSJ6uuvvzbaVSulXJ388L1tG6iUalJKqeXLl6sDDzxQLV68WCmlVCwWU3PnzlWHHHKIevnllxPt\nuWsG9LurN4tSqlYppd59912llFLNzc1q6dKl5odevXq1OvPMM9WZZ56p/H6/8fDzGdD3dG13KKXU\nihUr1KBBg9TDDz9s2mL58uXq2GOPVZdffrlqbm42Hr45A/qcju0c4wO+9NJLaubMma3GxZo1a9RZ\nZ52lpk6dmjguXlJ6TPV037tqW6qUUkcccYSaOHGiikT08e39999XkyZNUnPmzFHRaFQppVQ4HFY3\n33yzOuyww9SHH36oEuif9J4ppQnwo90wo9CrR18B7kafYq4C3gHOBqYDH6JPb64EmDFjBu+88w6j\nR49mypQpfP3116B9lDNS2G9v5iLA/uKLLzJlyhT69+/PggULALj//vu5//772WWXXbj55pv573//\na7S5sac6m0bcQH4kEjHXN7z66qv88Y9/BPRFwKOPPpoff/yRkSNHcvjhhxMIBACOAfbvqU6nkUJg\nJsDFF1/M6aefzjnnnANAZWUlxx9/PKFQiJycHE488UTjwv5VaDdeX8IH3KeUYubMmZxxxhnY7XYm\nTJjAJ598QiwWY/Lkyfzwww+MHDmSww47zBgXk4C+FPf4M8BBBx3E4sWLWbRoET///DMnn3wysViM\nUCjE1KlTicVilJWV8fTTT7PnnntywQUX8O677xrvcc0m76pa0xVHn32UUmrjxo1q5MiR5hu/8MIL\navr06cbdhzLgKJnO7Q6llHr++efV559/rh588EF17bXXqlgspvbYYw/166+/KqWUqqioULvvvrth\nk1UZ0O+u3pxK6bOTIUOGmJ956NChSimlFi1apM466yzj86urrrpKPfbYY8bdszKg/129XaCUUm+9\n9ZbKzs5WP/30k/nZb731VnXrrbea9ydPnqyWLVtm3B2TAX3vyu0EpZR67bXXlM1mU2vWrFGxWEz9\n5je/URdccIF6+eWX1dSpU01bXHPNNerRRx817v4+A/rfVdvxSin1/fffK0DNnTtXXXfddeqBBx4w\nP/u4cePUl19+qYYNG6ZCoZBSSqlPPvlETZw40XjJ4qT3VB0tPtoShoFelDR58mTzwS+//JK9997b\nuJtycppeyv3ArGOOOQa/38+tt97KypUrWbduHV6vl/79+wPbhE1KAdasWWPGM/t8PjP09JVXXuHE\nE080X/zVV19x2mmnGXcruruz3cCBAEuXLuXggw+mubnZrBb/8ssv8/DDDwOYC/123nlno137cW+9\nj2qAXXfdlUGDBnHssceSlZXFJ598wj//+U/mzp27ybg46aSTjLt96XfyIsDw4cMpLCzk888/Z9my\nZbz00ksAZiTPxo0bGTt2LE6nE+hYN9Ih6gp0eFBtrS7W8dprr7F48WKWLFlivGZhGvabSfxq3Jgz\nZw4XXXQRbreb2tpaAoEAsViM1atXc/nll7NwoWmKZ3qmq2mlGTBD8BYvXsz777/PunXrOOWUUxg4\ncCC1tbXEYjFuv/12Bg8ezKhRo4y2b/Rct9NCNnCMUornnnuOmpoaDjzwQGKxGPPmzTN/L5FIhBkz\nZnDeeefhcrlAhwav69mudzlvA5SWlnLooYeyePFiTj75ZH755Rf+85//mLaIxWLceeedDBgwgD32\n2MNo+3rPdbvLaUZH+ngTVzvX1taSl5fHeeedx+zZs3E4HNTW1qKUYuXKlfztb39L1NIFm7yrak1X\nnFLsopS+SLr//vurHXfcUU2fPl0Fg+YF7Qcz4LQn3dupSim1bNky5XQ61Q033KDmzp2rAoGAOu+8\n89TgwYPVhAkT1C+//GLY5HOllC0D+t3Vm1XpC+rquOOOU3vvvbc677zzlNPpVKtXr1ZffPGF2nHH\nHdWIESPUnXfemXhx8MwM6HtXbwVKKRUMBpXL5VL33Xefam5uVlOmTFHnn3++WrRokRo8eLDaZZdd\n1DPPPKMS2C8D+p6OrVwppXbffXd1+umnq8bGRjV27Fh1+eWXq6+++kqNGDFCjRgxQt1xxx2J42Ja\nBvS7q7dKpZQqKipSM2bMUAsWLFCDBg1SI0eONAMpotGomjx5sho2bJj67W9/q6qqqgx7JLteUMaf\nBLqqo48ppX2pSeGPjyulsjPAkOneTldKqY8++kidccYZ6rbbblPTpk0zfahr165VsZgZ5feNUqpf\nBvQ5Xdu1KomSkhJVXl6ulFKqsbFRbdiwIfHpmRnQ53Rty5VS6sILL1SjR49Wv/3tb9XQoUPV999/\nr5RSqqamRgUCgURbTM6APqdru08pHQU0ZswYNXToUDVr1izzd9HU1JQ8Li7JgD6nY/tJKaUeeOAB\n9fHHHyullKqsrFT19fUqmSQtfUMp5W7j/fSfBLqqozlKqZsS3jeg+ubsa3ObRyllTsPb4UallD0D\n+pvOza50WF658aG/+uqrtmzxheq7s1JjO039//buWKVhKIzi+N/VDiIduzj0BXTUwd0H6OJbWHwF\n94Bb3sDBRQQXX6KD4OjgIE4S6CCEz+G0IMW0FL2Ym5wfhECnm4/2prl8Jzd0IyuKIsqyjKqqfqrF\nY0SMWzDelMcw9Idmk6eIOG7BeFMdl2uu/T4iDiJiGhHzxWd1RFxEc2vnxkTpb+2hxNcLisb2yQiY\nLs5zlCh8Bq5RyvCNLV6n2REnwB1qa/0AJqgtdhcleduUdExhB7hCbYrfPQDnqAW0RgnEPthHiesz\n9Pt4R/sz3KIU+oB+fC9OgSM0H9ygdOwnSv8uDVDjwSuaT5okn9TNVg2BMdpgpUudDNs4RL34NbrJ\nzf53ONYhntTNzDrkzzaeNjOzFkjRp55azutrOT4J5VzvpVzqnnOtc6lxk5xqv7bWq8svZmaWMS+/\nmJl1yBey6ps9dAIlYAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1125ce1d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import scipy as sp\n",
    "import scipy.sparse\n",
    "import scipy.sparse.linalg as spla\n",
    "import scipy\n",
    "from scipy.sparse import csc_matrix\n",
    "n = 40\n",
    "ex = np.ones(n);\n",
    "lp1 = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr'); \n",
    "rhs = np.ones(n)\n",
    "ev1, vec = spla.eigs(lp1, k=2, which='LR')\n",
    "ev2, vec = spla.eigs(lp1, k=2, which='SR')\n",
    "lam_max = ev1[0]\n",
    "lam_min = ev2[0]\n",
    "\n",
    "tau_opt = 2.0/(lam_max + lam_min)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "plt.close(fig)\n",
    "\n",
    "niters = 100\n",
    "x = np.zeros(n)\n",
    "res_all = []\n",
    "for i in xrange(niters):\n",
    "    rr = lp1.dot(x) - rhs\n",
    "    x = x - tau_opt * rr\n",
    "    res_all.append(np.linalg.norm(rr))\n",
    "#Convergence of an ordinary Richardson (with optimal parameter)\n",
    "plt.plot(res_all)\n",
    "lam_max, lam_min"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Condition number of and convergence speed\n",
    "Thus, for **ill-conditioned** matrices the error of the simple iteration method decays very slowly.\n",
    "\n",
    "This is another reason why **condition number** is so important:\n",
    "\n",
    "1. It gives the bound on the error in the solution\n",
    "\n",
    "2. It gives an estimate of the number of iterations for the iterative methods.\n",
    "\n",
    "Main questions for the iterative method is how to make the matrix **better conditioned** \n",
    "\n",
    "The answer is <font color='red'> use preconditioners </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Better iterative methods\n",
    "\n",
    "But before preconditioners, we need to use **better iterative methods**. \n",
    "\n",
    "There is a whole **zoo** of iterative methods, but we need to know just few of them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Attempt 1: Different time steps\n",
    "\n",
    "Suppose we **change** $\\tau$ every step, i.e. \n",
    "$$\n",
    "   x_{k+1} = x_k + \\tau_k (A x_k - f).\n",
    "$$\n",
    "\n",
    "Then, $$e_{k+1} = (I - \\tau_k A) e_k = (I - \\tau_k A) (I - \\tau_{k-1} A)  e_{k-1} = \\ldots = p(A) e_0, $$\n",
    "\n",
    "where $p(A)$ is a **matrix polynomial** (simplest matrix function)  \n",
    "\n",
    "$$\n",
    "   p(A) = (I - \\tau_k A) \\ldots (I - \\tau_0 A),\n",
    "$$\n",
    "\n",
    "and $p(0) = 1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Optimal choice of the time steps\n",
    "The error is written as \n",
    "\n",
    "$$e_{k+1} = p(A) e_0, $$\n",
    "\n",
    "where $p(0) = 1$ and $p(A)$ is a **matrix polynomial**. \n",
    "\n",
    "To get better **error reduction**, we need to minimize\n",
    "\n",
    "$$\\Vert p(A) \\Vert$$ over all possible polynomials."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "##  Polynomials least deviating from zeros\n",
    "\n",
    "Important special case: $A = A^* > 0$.\n",
    "\n",
    "Then $A = U \\Lambda U^*$, \n",
    "\n",
    "and \n",
    "\n",
    "\n",
    "\n",
    "$$\\Vert p(A) \\Vert = \\Vert U p(\\Lambda) U^* \\Vert = \\Vert p(\\Lambda) \\Vert = \\max_k |p(\\lambda_k)| \\leq \n",
    "\\max_{a \\leq \\lambda \\leq b} p(\\lambda).$$\n",
    "\n",
    "Thus, we need to find a polynomial such that $p(0) = 1$, that has the least possible deviation from $0$ on a given interval $[a, b]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "##  Polynomials least deviating from zeros (2)\n",
    "\n",
    "We can do the affine transformation of the interval $[a, b]$ to the interval $[-1, 1]$.\n",
    "\n",
    "The problem is then reduced to the problem of finding the **polynomial least deviating from zero** on an interval $[-1, 1]$ \n",
    "\n",
    "with some normalization constraint $p(c) = 1$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Exact solution: Chebyshev polynomials\n",
    "\n",
    "The exact solution to this problem is given by the famous **Chebyshev polynomials** of the form\n",
    "\n",
    "$$T_n(x) =  \\cos (n \\arccos x)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## What do you need to know about Chebyshev polynomials\n",
    "\n",
    "1. This is a polynomial! (we can express $T_n$ from $T_{n-1}$ and $T_{n-2}$).\n",
    "\n",
    "2. $|T_n(x)| \\leq 1$ on $x \\in [-1, 1]$.\n",
    "\n",
    "3. It has $(n+1)$ **alternation points**, were the the maximal absolute value is achieved (this is the sufficient and necessary condition for the **optimality**\n",
    "\n",
    "4. The **roots** are just  \n",
    "$n \\arccos x_k = \\frac{\\pi}{2} + \\pi k, \\rightarrow x_k = \\cos \\frac{\\pi(k + 0.5)}{n}$\n",
    "\n",
    "We can plot them..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x11245f5d0>]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD9CAYAAABeOxsXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXl4FEX6xz+TyZ1JMiEXgoCACCiILi6oKCggIqiw4KqA\ngC676E9FvFHXhcDqKuqqoLseuyogICrgsewCKqCoICAeXMop95H7nslMJvX7o6Ynk2FIOhjIJHk/\nz1NP99RUd1fX9Hzfrusti1JKIQiCIDRpwuo7A4IgNEos9Z2BJspJl3v4rzxJY6lVNPQHt6H/Dg2x\n/Bt6mRucyrKv63M3ljKH0Cp3BVWNwem4qHBqkN/h9CNlfvqRMj+FGM1EUsiCIAhNmNrUDNKArkAJ\nsBFwA9HAMKA18B3wKbrK0QG41pvmQ+Ag2vBcD3QGtnvjK4AW3nNEAEuBn3/lPQmCIAi1xOI3muhE\ntYNI4GngPr+4g8CzwONAql/8CmAf8IeAczwPXIc2Egabgc8CzgswFZiGNhSCIAjCqUVBzcYgBvgY\n6O/xeFizZg0tWrSgffv2vgRbt27l008/ZezYsSQlJQFQXl7OggULsNvtXHvttb60+/fv5/333+fm\nm2+mZcuWAFRUVLBkyRIKCwsZNWoUFosFtKF5uM5vWRAEQQjElDF4Gph05MgRbrzxRvbt20d+fj7z\n58/n2muvZeHChdx4441069aN8PBwPv74Y8444wz69OnD6tWr6dWrF/369WPq1KmsXbuWvn37cuaZ\nZ5KQkMCsWbPo2rUrt912G7NmzeKyyy6jbdu2zJ492zAIvdDNURVAPFAMuE5BQYQBdnRz1Rnomk6S\n95ox3m0sEAfYvHEx3n0buqksyruNRDd3+Te/KcDjDcp7P8pvH3TZWwL2Ld68hQHWgDQK3QTnAEqB\nMnTZFHpDiTe+EMj3bvP84kvR5Vnsjc8Gijj9tTEbcCa6/OOpLNNm3rgE9G9hlHsc+reI8gYj3vhs\nRZdX4LNcQWW5G1uoLOPA9EZcsP+E8Xu6gXJvcKHL1IEu4yJ02eZ797PR5VwAZHk/ZwM53rjTOUrG\ngi7DZHQZJ6LLMQVd3tHeYEOXfZI3XZz3e+N/EUPls248oyfC//6q6580fp8KKsvXjS7fMnQZF1D5\nzBrPc543GM90EVXLuYz6GYlkQZfdGehn2j8Yz3aiX4hDl6uhJZEcryv+0wEMHfFQWW4GHirLogvV\na6ep0UR/BBgyZAjt2rVj2bJlvPLKK9xxxx0cPHiQ6dOns2PHDpo3b87q1auZMGECEydOJCkpifz8\nfOLj4xk8eDCrV69m5syZLFmyhB49enDw4EFGjRrFokWL2LZtG/n5+dhsNu6++27mzJnD2LFjAb7E\n3DyI17xbt/fm3d7jYtEFGeENUeiCTgKaUylAMSaucdK43W5LcXFxuMPhCC8qKqKkpITS0lLy8vIo\nKCigsLCQvLw8CgsLcTgcOBwOXC4XTqeTsrIyXC4Xbrcbj8dDRYX+rcPCwiwRERGR0dHRkXFxcYlR\nUVFERkaSkJBAQkICsbGxxMXFkZiY6At2u524uDjfd7GxsYbRNfAAud5wyLvNp1LEyrzBQeWf1ngI\nregHNsy7NQTaEI1Y9MOfji73NPQfok7LXinlK8+SkhKKi4vxlnlYSUlJmMPhoLCwkPz8fN93xcXF\nuFwuysvL8Xg8eDyeKvsGRllZLBaL1WoNj4iICA8PDyc8PJzIyEhiY2OJiYnBZrNhs9lISEggMTGR\n+Ph4UlJSsNvt2O12UlNTiYmJ8S97F5VGoshvv9D7uZRKQXR4fwP/so+mUphjqRQVm3drCH6Ct8zj\nvWl/dVnn5uaSk5NDYWEhhYWFFBQUkJ+fT3Z2Nvn5+RQVFVFaWorL5bI4nU7fs11eXu57lo130bCw\nMKxWqyUsLMwaHh5ujYyMjIiIiMA/xMTEEB8fT1JSEsnJyb7nOCkpiWbNmmGz2YiLiyM+Pp6IiCq3\n6KLSMBvBMBYF3pCL1g+jjEu85Q6VBsx4OYukUsBTvMHQFP8XmnTv71OnuN1usrKyaNGihf9LYzBs\n3nxlAI+dII3PSNZkDGIAjB9yxowZLFy4kN69e7Nz506Sk5M5++yzAUhPTyc2NpZ3332X8ePHk5iY\nCEBqqu5S2LlzJ3379sVisZCenk5cXBwffPABo0eP9qU1zuHFsnjxYjZt2sTUqVNPmMHRo0ffblwv\nLS2NpKQkIiMja7it2qGUwuFwUFRU5Hvos7OzOXr0KLm5uRQVFZGfn8+xY8c4duwYxcXFZGdnk5WV\nRVFRkalrREZGEhMTQ0xMDJGRkURHRxMdHU1kZCTh4eFYrVasVitKKcrLyykpKcHpdFJaWkpZWRll\nZWUUFhbidDpN31dcXBzNmjUjJSWFpKQka3JycmpKSkpqixYtOjZr1gy73Y6xjY6OJj4+HpvNRnR0\nNDabDau1upfB2uNyucjPzycnJ4eCggJKS0spKSkhOzubgoICiouLycrK8ol5cXExOTk55OTk+Ayq\n2Qn1sbGxvvsJLGP/fUO0lVIopaioqMDj8eB2uykvL6e8vByXy0VpaSlOp5Pi4mLKysqqvXZMTAwp\nKSmkpKSQmpoaabfbU1NTU1NtNhtJSUm+5zgpKYmEhASfcanrZ1sp5Xt2ioqKyMnJITMzk8OHD/tE\n3ijjzMxMsrOzfYY2Jyen2vu0Wq3YbDZiY2N9z7b/8xwWFkZYWKWGuVwu3wuPUaZGGbvdblwuFw6H\ng+Li4ipG+kQkJiaSmppqPNuRdrs9LSEhIS0uLo6EhARSU1NJTk72GZdmzZoRFxdHXFwcNpst0JjU\nKR6Px2dAS0tLcTgclJSUUFBQ4HtBOXr0KEePHqWgoICCggKysrJ8v4GhKVOmTKlWG1u1asXWrVuJ\nj48fzomNgY+ajMEvwHn/+te/uPjii9m4cSM5OTlccMEF5OXl+US8pKSE++67j5deeolnn30Wu90O\nwMqVK8nKyqJz585ERkZisVjweDyMHz+eSZMmsX79es455xwAtmzZwvLly5k0aZJxbcumTZtqyj9v\nv/02b7/9dtWbCg/3PYjR0dFERUURERFx3J/eYrH4/uDGQ2g8eMYbfFlZmakH0GazkZaWRvPmzbHb\n7bRv3963n5CQ4HtrNB64q666ioyMDN/x/vu/Brfb7RNR40Hyf1s2/lDGm7PxB8/JyWHz5s0cO3aM\nvLw8U9cy/tyRkZFERUURHR1NREQEP/30U5X7WbVqlU9MjTdu4wXD+CMYZV0T8fHxvrdAm81Gamoq\n5557Lna73VcLMv7Y8fHxxMfH+2pDxptlYmLirzZk1f125eXlvpeGwsJCsrKyKCgo8JVzVlaWT2Bz\ncnL45ZdfyMnJoaioCLfbTXUYwmqz2YiJicFqtfriIiIiCAsL8xkww3AZz7Xx+7tcLlwuF8XFxZSX\nlx93DaWU755sNhsLFiwgLS2NlJQU2rZtS1xcHMnJybRs2ZKUlBTi4+N9NSG73U5KSgo2my2w5lkn\nGLW/3Nxc33OTm5tLXl4excXFlJaW+l7WMjMzyc3NJSsri507d/qMW0lJCRkZGWRmZgIEFdTY2Fhf\nGRvPt1EDN7TDCP4aYhg0j8eDy+WirKwMp9PpCw6HI2iZByM5OZnExEQSEhJISUnh7LPPJiUlhWbN\nmpGWlsbRo0erPb5du3b8/e9/N60tNfUZzADueeKJJ3jqqacoLCxk0aJFjB07lk2bNjFw4ECuvfZa\ntmzZwqOPPkr//v154403eP311+natSsOh4OXX36ZhIQEevToQffu3dm3bx9jx45l5MiRrFixggkT\nJtC3b18OHDjAzJkzadOmDcB7wI0Wi6VG63frrbcycuRI8vLyjntrNN7WysrKfCJvNAEY1VTQTQBh\nYWGEh4f7jIZRBY2KivK9mSUkJPge/JSUFNLT00lNTSUuLo7w8NrP3/P/s4SSi6iysjLy8/N9Vf6C\nggLfW29RUZFvv6SkxNeUZTz0breb999/v8q99enTx/fZePOOiory1YSMsk5ISMBut/veio0/pPHZ\neIsPBU7Vb2e8oRsCZwhYUVFRleZEQ9iNt2iHw4Hb7T7uuTbKOyIiwmcQDQNuPNfGM26ITO/evUP2\n2awr/O/vv//9r083jJck47NhPI0alFGDCSwTQ0OsVqtvaxgPo5YfHR3tqyEZxtN4YbXZbCQmJlZ5\nyYmOrr6FKSMjo1ptNFBKbQc6nehr3z3UYAxuAd7es2cPl19+OV26dOHnn3/m3//+N1dddRVZWVks\nXLiQ4cOHk5aW5jvo888/Jzc3l6FDh/qqgsXFxcybN48BAwbQtm1bX9offviBH3/8kREjRhh/9G+B\nK4GipvRANrb7a8z3BnJ/DZ0mdn91YgzCgbXARSUlJezZswe73U6rVq2CnfRb4G3gBap2aExDzy8Y\n4Rd3EHgK3bHhP0/hJeABdMeNamI/WD3mpO5pzPcGcn8NnSZ2f3ViDED3lL8KDPeL+xn4HXrI1IXo\nCWSfer9LAG5CjyZZBBzxxv8GuALYS+Xs4yjgRvQwtyXALv9MNrEfrB5zUvc05nsDub+GThO7v11U\nnfDrT62MgUELtGi7gN3oMcCnEjEGDZjGfG8g99fQaWL3Z8oY1KbX87A3CIIgCA0HU5NJZXEbQRCE\nxo2pcdS/dj2DU8nwKVOmLHr66ac5//zz6zsvp4QpU6bUdxZOGY353kDur6HTFO7PTztN6Xxt+gxO\nNzcD76SnpzN06FBee+21Gg8QBEEQNH7aeQC9zEAwfH0GodxMZAU949LlOhX+6QRBEBovftppSudD\n2Rgo0H5cHA5HfedFEAShQeGnnaam7Ye8MYiLi6OkpKS+8yIIgtCg8NNOU96BQ9kYxAJERETU6LxL\nEARBqIqfdppywRrKxsAO2gOpWS9/giAIgsZPO02NJgpcNSeUiAXt5dKM/3JBEAShEj/tNDXPIJRr\nBvGgV0BqjNPFBUEQTiW11c6QNwaCIAjCqacujEEqcB5Vhy9Z0S5TWwaktQNd8TYBebEA7YG2VJ34\nZofG6URKEAThVFNb7ayuYyENmAJ0oVKkI9Fuq1PQ7qYrgG7e7/OBuwAn8CJgLHrwKXAvep2CMd5r\nlgOPAOvRC9p39qbdCNyGdomdAnrpvpNZRUwQBKEp46edpqzCiVQ2DlgNdKzm2PNBL2R96NAh2rZt\nawfmGV8eOnTIWMbtKmCrkbndu3fToUOHcOA5I212djYWi4Xk5OTuwCbgEvS6CLjdbuLi4szciyAI\nguDFTztNuXA4UTNRf6Dj9u3bueKKK+jdu7cvXHrppaSkpOB2u1m3bh1nn302nTt35qGHHvItev7w\nww/TunVrunbtyvr16wHYvXs3559/Pl26dOGWW26htLQUgBkzZtC6dWs6duzIf//7X+P6nwDpAA6H\ng5gYU3MmBEEQBC9+2uk0k/5ENYODAG3btvWtbxwZGYlSittuu40ePXoQHh7Ogw8+yNSpU7nkkkuY\nOXMmM2bMoHfv3qxdu5YFCxbQvXt3fv/737NkyRKmTZvG6NGj6d+/PytXruSRRx7h8ccf51//+hev\nv/46/fv3Z9iwYZx55pl069Yt3sjbpk2b6qJcBEEQmhQul8tYV96cczdVFfxCRsB3at26dcpqtaod\nO3aobdu2qX79+vm+W7t2rRo3bpy6/fbb1ccff+yLv/nmm9WXX36punTpotxut1JKqWPHjqkrr7xS\nTZ8+Xb344ou+tI899pj68MMPjY8VSik1ZcoUhW7zChqmTJkSmE1BEIRGT03amJSUpEaMGKGUUodV\nVW33Dz6qG02UgV7feJoR8dxzzzF8+HA6dOjAnj17OO+88wyDwvPPP88tt9xSJX7nzp3s3buXZs2a\ncdZZZ/k6gv/+978flzY3N5fly5dz5ZVXGpcLNZfagiAIDYby8nKjmciUp8+ahun8gHZyNDkzM5OP\nP/6YJUuWANCyZUu++OILVq5cyTvvvMMFF1zAFVdcwezZs5k7dy7nnnsur776Km+88Qapqals27aN\n5cuXs2rVKoqLixk7diyHDh1i4cKFOJ1O/vGPf/Dcc8+RkJAA8C1w0UmXgiAIQhPH5XIZHcjmPH0G\n1DyCVSMeUUqp2bNnq4SEBFVaWupLPH/+fDVo0CC1cuVKX1x2dra644471J133qlyc3N98Z988om6\n5ppr1LvvvuuLczgc6tFHH1WjRo1S+/bt88/H1cZOZGSkevjhh2tZgRIEQWh6fP65UhdfrPf9tHOd\nMtFMZGYAfzuA7777jssuu6zKyJ4RI0YwYsQI/7RFycnJ8a+88op/XBaQetVVV3HVVVf5x2dHR0en\n/O1vf/OP2wn0wzuk1eVy4XK5jNqCIAiCUA3ffw8JCcdpZ5GZY80Yg0Og19SMiPB5Qs0A5gAPo4X7\nW+BZIBu4AT1xTAGzgIXoCWQPAr8FtgPPAL+ghf9u9LyGhcAbgAfoAZCXlwdAUlKSmXsRBEFo0mza\nBElJx2lnnpljzRiDGcCApKSkS72fFwLT0WNX/y9I+ve9wZ8sYFKQtCu8IRA7QGFhIQDx8eKmSBAE\noSa2bYOuXY/TzkIzx5oxBvnAZehJYG4g5+SyWSsSAIqLiwGw2Wyn4ZKCIAgNF48HNm+GSy45TjtN\ndSCbdfqjgKMnlcOTIw4qrZv0GQiCIFTPzp1QWgqtWx+nnQVmjg9VF9Y2kD4DQRAEs2zdqrft2h2n\nnflmjg9VYxADlTfUrFmzes2MIAhCqLNli96eeeZx2plr5viQNgaZmZkAJCcn12tmBEEQQp1t2/S2\nefPjtNNUP2+oGoNY0NYtIiJCRhMJgiDUwM8/621i4nHamW3m+FA1BjGge8TFEAiCIFRPRQVs3673\nY2KO085iM+cIVWOQCJCVlSVNRIIgCDVw8CCUlekmIqv1OO1s0H0GiQD5+fnSeSwIglADe/bobbt2\nehugnQ26mSgOpGYgCIJgBsMYdOmit37a6QTKzJwjVI1BDOi1kdPS0uo7L4IgCCHNvn1629G7ar2f\ndmaiJw3XSKgag1ilFNnZ2VIzEARBqIH9+/W2XTu92JifdprqL4DQNQbRDocDh8NBampqfedFEAQh\npDGMQatWEKCdpn3JmfVNBLpT92GgO7AFmIoesnQF0BrYAHinPZAC9EcvxPwplf60fwucB2z1pgft\neuJqb15WoDs7IrKysgCZcCYIglATBw/qbfPmur8AfNpZkysKC95mJLPGoD2wEi36oMV7JNoBUie/\ndLOA3ej1DqzeOBcwERiF9n5qsAxYgl7bINYv/k9AeH6+vgfxSyQIgnBiKioq+wySk2H79iraaW7J\nS8w1E1mA94DWGzdu5LbbbsPtdgOcAXQ6cuQI7777Lk6nE+BW4K9KKeuKFSv49ttvASKBV4DL8vPz\neeeddygoKAAYCLwMxK5fv55Vq1YZ1/sXfsYgMTHR7L0IgiA0OY4c0XMMUlMhOloPKwWfdpoaSQTm\njMEw4DcHDx6kT58+dO/e3bfi2fLly2nXrh2PPvooAwcOpKCgAKUUN9xwA/3792fMmDG8+uqrAPz4\n44+0b9+e+++/nwEDBrB3714AHnzwQXr27Mkdd9zBo48+alzT4nA4AKTPQBAEoRqMWsFZZ+ltdrae\nVuDVzjqtGdwHMGvWLM4//3zuvvtu3xd//etf+fTTT/nmm2/44x//yMSJE/nuu+8oLCzk559/5vvv\nv+edd97h22+/5e9//zsvv/wyGzdu5IUXXmDMmDEcPnyYlStXsm3bNjZv3szevXtZtGgRULmGQbdu\n3czeiyAIQpPj0ktBKVi/Xn8ePnw4SilDQyvMnqcmY2AFegHMmDGD+++/3/fFL7/8QlRUFJdddhlp\naWl069YNl8vF22+/zYQJE+jYsSNRUVF06NCBoqIivv/+e2666SZatGjhS7tw4ULGjBlD586diYyM\npGvXrrhcLkD3iANkZGRgsVhOGDIyMszeqyAIQqOhJm30tspE1HQeg5qMwSCA/fv3k5ubyzPPPMPg\nwYNZsmQJmZmZpKenA+B2u3nwwQe56667qsRv3LiR7du306lTJ+Li4ggLC0MpxT333MOECROqpN23\nbx+LFy9m8ODBAOzevbuWRSMIgiAYREVFQS1GjNaUsDfAnDlziImJ4Y033mDx4sVcd911rF+/njVr\n1jB9+nS++uorbr31Vnr16sXatWuZNm0a3bt3Z9OmTSxYsIDk5GQKCwuZOnUq27Zto2fPnowcOZIP\nPviAZ555hj179vD1118za9YsX/PQ7t27jZsRBEEQaklYWBhAtOkDVFUICP+nlFJffvmlioqKUt9+\n+62aOnWqOuecc5TT6VQ7duxQjz/+uNq+fbvvBB6PR82bN0+9/PLLyul0+uIPHz6sJk+erDZs2FDl\ngkuXLlVPPfWUKigoMKIOKKXUnXfeqZKTk5UgCIJwYnr1UgqUWrVKfw7QztnqeF0PDEr57/h99g92\npVSZx+NRDz30kEpLS1O333672r9/f7A8ZSqlXg4SP0cp9VWQ+BeDxC1UStmUUhW33XabatmyZbDr\nCIIgCF7atdPGwHgnD9DOxcqkMaipmSgfuDssLOz1Z555hmeeecb/u6+AeegJaVuA+YAbeBY93yAS\nWABsRs9VuBboA/yCnpxWAjwBjAOSgY+Ar42Tl5aWEhvrPxdNEARBCCTH63AiJUVvA7TTZvY8ZjoX\n/gUcAf6MFv5dwBzg30B5kPT70K4q/FHAf7zBn2xgerCLulwuIiMjTWRPEAShaVJeDkVFYLGAsbBZ\ngHaanrVrtqd5iTecNsrLywkPr43rJEEQhKbFsWPaHUV6OnjnAgdqZwuz5wpVr6ViDARBEGogsIkI\nghoDU3MNQtUYWNxut8/thSAIgnA8RV5/0Da/noEA7QwD7GbOFarGAI/Hg9VqrTmhIAhCE6W4WG+N\n/gIIqp2mRuKErDFQShmTJgRBEIQg5OXprb+n/yDamWDmXCGtthaLpb6zIAiCELLo1QAg0NN/gHbG\nY4KQNQZKmVrDWRAEocmS613h2O7XKxBEO5uZOVfIGgNBEASheg4d0tuWLatNZmquQUgbA6kdCIIg\nnBjDGJx5ZtX4AO1s0KOJlNVqxePx1Hc+BEEQQhajz8C/mSiIdjbomkGFGANBEITqOXJEb73LwgBB\njUGtRxOF0tCdivDwcMrLg7k+EgRBEJSCgwf1vn8zURDtjDNzvlCtGXikZiAIgnBicnP1pLP4+Bqb\niUwtcFMb5z8W4DygK9rb6AqgJTAaSEe7n37fm/YaoD9QiHZzvRM9vGks0BbYhHZjXQ5cCgwFPMBC\nYCPSTCQIglAtBw7obevW2mupwck2E9VkDNoBrwA9qLlH+h7gEaAC6O4XPwX4GOhLVd/ajwH70Wsc\nGDwC/A0ol2YiQRCEE5OVpbdpaVXjg2hnspnzVddMZAHeAwbgNQQ5OTksWbKEw4cPA3r40ocffsif\n//xncrT7vAuB7jk5OUybNo333nvPONf1gG316tU89NBD7Nu3D3QNoU9xcTHPP/88//znP43hUI8B\n1sjISFwul5l7EARBaHJkZuptoDEIop2mFripzhi0BroXFBSQlpZGYmIiZ511FrfddhsDBw6kvLyc\nJ598ktGjR/PDDz/Qv39/srOzKSkp4be//S0vvfQSc+bMISMjA4B58+YxaNAgNmzYwJAhQ/jpp58A\n6NevHxkZGaxatYpbb73VMAhxbdq0odjwwiQIgiBUwRhJdMYZVeNtNlugdsaYOqHfUpqB62Jeo5RS\nX331lUpLS1M//PCDys/P9yUuLS1VHTp0UE6nU7ndbrVgwQL1pz/9Sb355ptq0qRJyu12K4/Ho3r1\n6qW2bNmiunfvrvbv36/cbrdav3696tevn/r888/V73//e+V2u5VSSo0aNUp9/PHHSiml3nrrLRUb\nG1sna4QKgiA0YX5WJtZArq5m4AaIj48nMzOTcePGMW/ePF9b1Ndff82VV15JVFSUbyGF5s2bs3jx\nYsaMGeOLCwsLo7i4mGbNmtGqVatq0xrxoNu9SktLycjIwGKxnDAYNQ9BEISmRE3aOHnyZCOpqdFE\n1RmDrwC6du3KvHnzuO6663j33XeZMGECACUlJURFRQGwf/9+nn32WSZOnFgl/qWXXuKiiy4iNjbW\ntyZnfn4+9957L48//niVtB999BEOh4OLLrrIl04QBEE4OcrKyozdKDPpqxtN1A60K9SRI0cCcMEF\nF/D0008DcOmll3LfffeRmZlJTk4O8+fPJzk5meHDhzNs2DDS09O5+OKLeeaZZ6ioqGD//v3ccMMN\nHDt2jBdffJFOnToxdOhQ7rrrLjp27EjLli2ZM2eO4Xr1+5iYmAtPsgwEQRCaPH7GINJMeotSPo9G\ngTOQ/wJMW7p0KUuXLiU9PZ25c+dy0003+ZpmCgsLWbVqFYMGDaqyROXGjRuJiIjg/PPPr5KxZcuW\n0b9/f+LiKifE7dy5k2PHjnHZZZcZUa8Blvnz548fNWoUP/30E506dTJ5+4IgCI0flwuiosBqhbIy\nvTWYP38+AdpZSvWzkBVU30x0EKBbt260b9+ejRs3csstt/i3QzkTEhIYMmSIYQj+DdwMFHXv3t0w\nBB7gRuCtqKgohgwZYhiCj9HDTfd16NDB3xD8CbgDKDEMRklJiZmyEQRBaDIYE85atKhqCACCaKep\nxeSrayaaC4xs0aJF/4kTJzJx4kT/7yYCM4HfAOcAG4Dd3u8+Bq5GTz5bDpShZyZPR89D+An40Zt2\nKXAV2motB7zLO1Nq3FBpaamZ+xAEQWgy7Nypt+3bH/9dEO00tZh8dcbAjZ5wdjUwEOgA5AGv4u1c\nBr7zBn8cwIdBzrfdG/wpRxuEQIrjvSs8FxUVBflaEASh6bJ1q96ee+7x3wXRTlM+6GpyR6GAZd5w\nOimy2fSkOTEGgiAIVdmyRW+7dDn+u5PVzlD1WlosfQaCIAjB8TpwCFozOFntDFVjUJCUlARArrHi\nsyAIgkBFRfXNRCernaFqDPISEhKwWCwUGOu6CYIgCOzdq9cxaN4cUlOP/z6IdppaTD5UjUGhxWIh\nISFBjIEmISepAAAgAElEQVQgCIIfP/ygt926Bf8+iHZWVHM6n6EIVWOQB7q6I24pBEEQKjGMwQUX\nnDhNgHaaWiUsVI1BLkBycjJZxgoOgiAIAt9/r7cXVuOwJ0A7G7QxKARcycnJ5OXl1XdeBEEQQobv\nvDO7ajIGftrZoI0BQGmQRRoEQRCaLIcP65CYCGeffeJ0AdppasnIUDYGRWIMBEEQKtmwQW9/8xsI\nq0a9A7TTbebcxukCPZaGAoV2u13mGQiCIHj5yusI6NJLq08XoJ21GlpqKvFpJis1NZWioiLcblOG\nTRAEoVHz5Zd6e/nl1acL0M5a1QxCkULxTyQIgqBxOGDjRt08dMkl1acN0M7T1mfQHu3dtLlfXCRw\nOXApVZ3htfSmbROQh55AHyDGL74gMTERQEYUCYLQ5Nm4EcrL4bzzICGh+rQB2mlqNFFNXkstwP3A\nGCALOAZ08e4/D9yDdnFt8ALwDfBPINkbdwS4Hb3IzS1+ad8B5qBXNmvtjXMCo4GFQHZKSor/DQmC\nIDRZjCaiXr1qThugnXViDIYBz53gu34ADoeDDRs20KtXL6xW633Gl9u3byc8PJz27dufgV7wBrfb\nzZo1a7j44ouJiooaAYwA2L9/PwUFBXTt2jUavRDOePw8l8qIIkEQmjqffKK3ffvWnDZAO+tknsFQ\ngL/97W98/fXXAIwbN47Dhw8D8P3339OuXTuGDh3Kdddd52vbnzRpEp07d+byyy9n7ty5gBb8rl27\nMnjwYPr168ehQ4cAePnll+nQoQN9+/blqaeeMq77OhAdGxsLyGpngiA0bUpK9EiisDDo37/m9AHa\nWZ1vIh81GYMYgIKCAp94f/311z5jMG3aNKZPn86GDRsYPnw4EydOZO/evaxYsYJPP/2UHTt28Prr\nr7NhwwZmzJjB2LFjWbduHVOnTmX06NEUFRXxyiuv8L///Y9du3axceNGFi9ebFy7bWJiIkopBg0a\nZOZeBEEQGiVr1uj+ggsuAK+H6mq5+OKLAV8zkam+4ZoSJQAcOXLEFzF8+HB27dpFVlYWBw4cYPTo\n0bRv356rr76ao0ePMnv2bO6880769euHzWajZ8+eZGZmsnTpUu69917OO+88rrjiCnJzc/nwww/5\n3e9+R79+/UhMTKRfv34cO3bMuJQrOTmZjIwMLBbLCUNGRoaZ+xQEQWiwLPUuDjxgQGVcTdoIPmNg\nag3kmoyBFeDQoUN08a6vVlRUxOLFi9m7dy9nnnkmFosFpRSPP/4448aN4+DBg7Rt2xbQTUOrVq2i\nZ8+ehIeHExOjBws98cQTjBo1qkra/Px83nzzTYYOHWpce4fdbjdzD4IgCI0ao79g4MDaHef1XBpp\nJm1NHcg5ABdeeCFDhw6lqKiIo0eP0q1bNxISEvjxxx+ZPXs2n3/+OR07dmTYsGHs3r2bV199lZ07\nd7J48WLefPNNkpKSKCsr45VXXmHnzp243W5eeOEFPvzwQ/7xj38QHh7Oe++9x9NPP80ZZ5wBsA1Y\nGxlp6h4EQRAaLUeP6pXNYmLA2/pjCqvVSmFhIUCEqQNUJQQJv1dKqWPHjql58+apTZs2qQkTJqiX\nXnpJKaXUxo0b1bhx49Tq1at9J3G73WrGjBnq4YcfVrm5ub74Xbt2qfHjx6sPPvjAF1dRUaHmzp2r\n7rzzTnXgwAG/rKh2Sqn+SimVlJSk7rrrLiUIgtAUmTtXKVBq4MDaHeennQdUcH33139lUUoZriiC\n+ScKA5YBVxkRc+fOJTExkeuuuy6YbVmBd8ipH98CKcBZAfGfoCeg+bMNPW/hINAXWNGqVSv69+/P\nW2+9VbNlEwRBaGSMHQtz5sBzz8EDD5g/zk87j1F1UrA/PldENTUTVQADgeuAa4Dbb7nFN2+sHHgJ\nSAd+Av6BXqHsfOAOIA54F/if9zpj0IbigPe4Q2gDMcF7jqXoiWjGMCgFEBMTg8PhqPHGBUEQGhtK\nVfYXDAh8da4BP+001UxUkzEALc4fecP7wL1AKTAN2Bok/SbgzoC4cuBNb/BnL1CtrRNjIAhCU2Xd\nOt1n0LIleMfwmOZUGAN/VnjD6cAKEB0dTVlZ2Wm6pCAIQujw3//q7bBhYKnlQgN+2mlqJE4oey2N\nAd0jXl5eXt95EQRBOO0YTURXX119umD4aaepl/5QNgbRoG/I4zHlWkMQBKHRkJcH334LERHQp0/t\nj/fTzjqZdFafJAK+mXSCIAhNiU8+gYoKvaqZd3mCWlFb7QxlY5AE+GY4C4IgNCXef19vr7/+5I6v\nrXaGsjGIq+8MCIIg1Adud2V/wfDhp+ea/sYg1F6/bQAej4ewsFC2WYIgCHXLN99AURF06gRt2tSc\nPhh+2mlK20NZZZMAXC4X4qNIEISmxGef6a2ZtQtOhJ92nrY1kE8VSQBOp5Po6Oj6zosgCMJpw3BZ\nfTJDSg38tNNpJn0oGwM76KUyIyLMOd0TBEFo6Bw6BBs2QGysuSUuT4SfdrrNpA9lY5AIYgwEQWha\nLFumt/36aYNwsjQmYxAH0kwkCELT4n//09trrvl152lMzUTRAGVlZURFRdV3XgRBEE45bndl5/Gv\nNQZ+2mnKuVttHdUBRAEetCfSaGAY0Br4DvgUPYypA3AtunryIXp9gjDgeqAzsN0bXwG08J4jAu3G\n+mfvdWJBRhMJgtB0WLMGCgv1kNKzzvp156rtaKKajEEY8CAwAi3OdiANbWl2A22oOjlsBbAP+INf\n3EvA8+g1ETr4xW8GPgPu84t7HpiKdo/tMwZSMxAEoSnwn//o7eDBv/5cftpZJzWDKcDkwAtYrdYo\nq9V6LsDWrVv59NNPGTt2LElJSf0AysvLWbBgAXa7nWuvvRbgfoD9+/fz/vvvc/PNN9OyZcuuQNeK\nigqWLFlCYWEho0aNwmKxTEEbgmiPx4PT6SQuTiYjC4LQuFEKPvpI75+sCwqDAO0sNXNMTX0GfwQY\nPnw4nTp1olOnTqSlpXG9N6cLFy6ka9euzJ49mwEDBnDkyBEA+vXrx+jRo3n66aeZMmUKAGvXrqVj\nx468+uqrXH/99WzevBmAcePGMWTIEF577TXGjh1r+NJ4CLAWFxcDYDsZL02CIAgNiJ07YdcuaNZM\nO6f7NQRoZ5GZY2qqGWQDLSZPnozD4WD37t3ceeedPoGfPn06O3bsoHnz5qxevZoJEyYwceJEkpKS\nyM/PJz4+nsGDB7N69WpmzpzJkiVL6NGjBwcPHmTUqFEsWrSIbdu2kZ+fj81m4+6772bOnDmMHTsW\n0L3h4qROEISmgOGL6KqrIPxkenP9SExMBKjT0UR3A3Tr1o2oqChGjx5NYmIiFRUV7Ny5k+TkZM4+\n+2xsNhvp6enExsby7rvvMn78eBITEwkLCyM1NRWAnTt30rdvX+Lj40lPTycuLo4PPvjAd06r1eo7\nh0FxcTEZGRlYLJYThoyMjFoXlCAIQqhhjCK66ipz6WvSRoD4+HiAYjPnq8kY/GLsdOrUiS1btjBt\n2jQuu+wytm/f7rM+JSUl3HfffTzwwAPk5+djt9sBWLlyJVlZWXTu3JnIyEgsFgsej4fx48czadKk\nKmm3bNnC8uXLjT4GAIqKTNVuBEEQGjTl5fD553r/1/gjCqQum4l6AqxatYoHHniAK6+8kp9//pnx\n48fTu3dvJk6cyMSJE9myZQuTJ0+mW7duXHnlldx333107doVh8PB/PnzSUhIoKKigvHjx7Nv3z7G\njh3LtddeS0xMDBMmTOCbb77hwIEDLFiwgJiYGIB1QE+n01TtRhAEoUGzcSMUFMDZZ5+8l9Jg1KaZ\nCFUVAkIHpVSFx+NRW7duVVOnTlUrV65UFRUVSimlMjMz1T//+U917NixKidZtWqVWrRokfJ4PL64\noqIi9eqrr6o9e/ZUSfv999+rWbNmqbKyMiNqg1LqMaWUWrp0qQLUV199pQRBEBorf/ubUqDU+PF1\nc74A7XxcHa/tRvBRkzFAKWU2exuUUvcopTwB8VOVUvMD4g4ope5USmUGxM9USkUopaYppdS8efMU\noH7++WeTWRAEQWh49O+vjcGCBXVzvgDtvEOZMAZm+qxfB+agZw93Aw4AC9Ejja4ELkRPIPvUm34W\ncBN6pvIi4Ig3/jngCmAvlbOP3wBuBJKBJcAub9pEgPz8fABfv4IgCEJjo6QEVq8Gi0U7p6sLArSz\nwMwxZgcwOYH3vMGfVd7gTyHwryDn+M4b/CkD3g6S1g6Qm5sLQFJSkslsCoIgNCxWrQKXC3r0gJSU\nujlngHbmVpPUgncltFB1VJcCcPjwYZKSksQ3kSAIjRbDZfXAgXV3zgDtPGbmmFA1BnEA2dnZpKen\n13deBEEQThmfehvYf82qZoEEaOeR6tIahKoxiALIycmRJiJBEBothw/Djh0QH6+bieqKAO3MN3NM\nqBoDO0BeXh7Jycn1nRdBEIRTgjHr+PLLf70LCn/8tLMEk15LQ9UY2EBXdaRmIAhCY8VoIjLrgsIs\nftppaiQRhK4xiAHdI55SV93rgiAIIcYXX+jtr1n4Phh+2pln9phQNQbRSilKSkrEfbUgCI2SAwd0\nSEyELl3q7rwB2ukwe1yoGgNraalej0EWthEEoTGyfLne9u4NYXWoxAHaadrBW6gag/CSkhJAFrYR\nBKFxsmSJ3tbFEpf+BGin2+xxoWoMfDUDrxdTQRCERoPbDStW6P1Bg+r23AHaWW72uFA1BhbDukkz\nkSAIjY1Nm6C4WLusbtWqbs8doJ0N3hj4HC3J0FJBEBobxkI2l11W9+cO0M6GbwwKCvTw2ISEhHrO\niSAIQt1yqoaUwnHaaTF7XG3mvCUDDwPnAD8DzwCtgXHAGejVyV5CW6Jbgf7o2W+zgK+AdsDtwNlo\nl9cvoD2cDgOGeK+xEPjY/4aMpTUFQRAaA0rB2rV6v1evuj9/gHbWuTHoCSzHu86Al0cC0twATAZc\naMNhMA7YhzYcRsaGAVOAw0ALv7SjgfmAKigosIAYA0EQGhfbt0N2NqSnQ9u2dX/+AGOgzB5nppmo\nC/ANkLhhwwaGDx/O3r17AfB4PLz++uuMHTuWPXv2AMQDyfv27ePOO+/k2WefRSkF0EYpZVm0aBEj\nR47ku+98yxq0yM7OZtKkSTz66KO4XC6AkXBcJ4ggCEKj4D//0dt+/fSCNnVNgHaaFlAzxmAywPz5\n8xk6dCiHDh1i0aJFAEyZMoWnnnoKgOuvv55du3ZRWlpK7969+fzzz9m1axe33norSineeustJk6c\nSFlZGbfffjurV69GKcWAAQNYtGgRDoeDwYMHGwbBEh8fD0jNQBCExsWHH+rtsGGn5vwPPPAASili\nY2MBTHe6mjEGdoADBw6wZs0arrzySpKSknA6nSxYsIAdO3Ywe/ZsXnzxRf7yl78wf/58brnlFrZt\n28Zrr71Gbm4u69at4+WXX2bdunUsWrSId999l4ceeohVq1bRsWNHdu3axYsvvsi5557LwoULAYiP\njyc6OpqMjAwsFssJQ0ZGxkkUlyAIwumnoADWrQOr9dc7p6tJG5988kmoY2PwGsCkSZNwu92sWLGC\nUaNGsXbtWq644goiIiIAyMzMpEOHDixZsoSbb74ZgIqKCvLy8ggPDycpKYmWLVueMK1/PGhHSzKS\nSBCExsSKFeDxwCWXwGmStzo1Bl8ZO48//jj33XcfUVFRuN1uysv1ENbt27fzwgsv8OCDD/rilVI8\n8cQT9O3bF7vdjtutZ0UfO3aMe++9l2nTplU5x9tvv01ERAQXXXQRAJs3byYqKsr0HQuCIIQ6hguK\na645bZdshtkRRaoqBAmPK6XUJ598omJiYtSiRYvUf//7X1VcXKwuuOACdfnll6tBgwapgwcPKqWU\nWrRokWrfvr3q3r27euGFF5TH41EVFRVqwIABqmfPnuryyy9XP/30k1JKqW+++Ua1bdtW9ezZUz3y\nyCPK6XQa+fCMGDFCnX322UoQBKExUF6uVHKyUqDU1q2n7jpBtDNGBdd2IyillDIztDQW4MiRI/Tt\n25evvvqKrVu30qdPH9avX8/333/PRRddRJjX7d6wYcP47W9/i1KK1q1b+06ybNkyvv32Wy644AJf\n01LPnj3ZsGEDWVlZdOrUyUj6N+CRkpISowNEEAShwfPjj5CTA23aQOfOp+46QbQzCROurM0YgwXA\nw2PGjLGOGTPmuC97VC7c+RR6wti/W7VqdZ43bj9wC9DXYrFk/Pa3vzXSvg68CryYnJzc2295tluA\nD4FHnU4n0dHRJrInCIIQ+hhNRH37npohpQZBtLMlek5XtZgxBpuAvsCD6ElnZegF69cBbwGpwDYg\n25u+K3CRd/9b9KSHL9Gzk7sCu4GD3u/7AOej5yeso9KPhsXtdvtqEIIgCA2ZigqYNUvvjxhxaq8V\nRDtNLSRvdgbyam8Ixk8BnxWwIUi6XOCLIPGbgp20oqLC1/QkCILQkPniC/jlF+2h9FT4I/IniHbG\nmzkupNXWcirrUoIgCKeJefP0duxYPcfgVBOgnaY6X/2NQUgpb1hYGB6Pp76zIQiC8KtwOsE7l/aU\nNxFBUO2stTEIJSqsVqsYA0EQGjwffaRnHnfvDueee+qvF0Q7U8wcF6rGoDwiIsI3IU0QBKGh8u9/\n6+2tt56e6wXRzgZtDNyRkZGG0zpBEIQGyc8/w2efQWwsjBp1eq4ZRDvtZo4LVWNQFhUVhdPprO98\nCIIgnDSvvaa3I0fC6VrBN4h2mrpyyBqD6OhoMQaCIDRYiorgzTf1/v/93+m7bhDtbNDGwCHNRIIg\nNGRmzYLCQrj8cvjNb07fdYNoZ4MeTVQSExODw1GjOw1BEISQw+OBGTP0/sSJp/faQbTTlPvnUDUG\nDjEGgiA0VN5/H3bv1mscDxlyeq8dRDtNLX0ZqsagKCEhAZfLRVlZWX3nRRAEwTTl5WAswPjooxBu\n1ulPHRFEO21mjgtVY5Cb5O16z8vLq+esCIIgmGfePNi+Hdq1O31zC/wJop2m3D/X1hhEAGmcetcV\nhcaSl4WFhaf4UoIgCHWD2w3Tpun9KVOgPhwvB9FOU96QamMMRqLdVB9Dr1MwEPgjsBbYA8wDWnjD\nq8BO4DvgTrTxGAqsBPYCHwGd0S6xn0F7Pt0CTEYbnMLExEQACgoKapFFQRCE+uPf/4Y9e6BjRz23\noD4Iop2mGqrMtma1BN5USkXl5eXRrFmzM4GlAWnaog1GIP/wBn/aANcHSTsVuAFYYlR18vPzTWZR\nEASh/sjJgb/8Re8/8cTp7yswCKKdpnJitmbQBYj64osvSEtLY6HXBV9eXh5/+tOf6Nq1K59//rkv\n8WeffUaPHj0YNWoUJSUlAJSVlfHII4/QqVMn3nvvPV/a7777jr59+zJw4EAyMzNBL4Bzo7Fsm3G8\nIAhCKDNpkjYIV14Jw4fXXz6CaKepZn2zxmA7wMUXX4zdbmfq1KkA3H333eTk5DB+/HgeeOABVq9e\nzeHDhxk/fjydOnXi8ssvZ9CgQTidTp544gm+++47br31Vt58803mzJmDy+XixhtvJDU1lVGjRnH1\n1VeTk5MD0D4lRftWGjp0qNkyEARBqBfWrIE33tB9BP/856ld1rImfvOb36CUqv2LtKqEGsI3Sil1\n++23q9TUVHXkyBF17rnnqoqKCqWUUl9++aUaOXKk+utf/6pmzpzpO+kNN9yg1qxZozp27KgKCgqU\nUkodOnRIXXzxxWrevHnq3nvv9aV9+OGH1TvvvKOUUionJ0cBasqUKQq9elrQMGXKFCUIglBfuFxK\nnX++UqDUY4+dvuvWpI2DBw82klao6rVdKaVUbTqQlwOcd955ZGVlsWzZMvr37+9bUefrr7+mT58+\nfPHFFwwePBiA0tJSfvnlF6Kjo+nQoYOvlztYWqUUGzZs4JJLLgGQ9Y8FQWgQTJ8OmzbpCWZ//nN9\n56YSPzfWykz62hiDHP8P0dHRHDhwAKUU77//Pp999hl//OMfffGFhYWMHj2ahx56iMTERI4cOUJ5\neTlffvklzz//PI899pgvrdPpZMKECVx99dW0adMGvI7qBEEQQpl16yonmL3+unZVHSr4TTqrMHWA\nX62jpmaiCUopNXPmTAWoI0eOqBEjRqjk5GQ1adIk5XK5lFJKrV27VnXo0EGde+656quvvvKd/IEH\nHlDNmjVTY8eOVSUlJUoppXbu3KkuuOACddZZZ6lFixb5mpyUUncppVR8fLy677776qRKJQiCUJcU\nFCjVrp1uHrr//vrOTVUCtLNUmWgmsiiljCpETV0eo4E5u3fvZvPmzQwZMgSLxUJJSQlxcVVdX7hc\nLiwWy3FNPcHSVlRUUFZWRkxMjBE1DD0/YW/79u3p2bMn8+fPN2XYBEEQTgdKwZgxMHcuXHghrF0L\nUabcwZ0eArSzgOoXuFFQu2aiRcBP7du3Z+jQob6+Aq+4/wxch55wNjcyMtIwBHOBdPSEs5+8aQ8B\n49ATzmaGhYW5vYZgCXA28AGQCZCSkmKMLhIEQQgZZszQhiA2FubPDy1DAMdpp6mFYWozLaIU+A3Q\nH+0f+zvgF7Tfi1IqOylGA3/yfjYarT4C/gPEeDNmrNY8EXgYbZT83ew5AY/NZrMWFxfXIouCIAin\nlg8+gPvv1/tvvgmdOtVvfoJhs9nw005TIlpb30RO9Bv8e8AutKiXcHxvtZNKQ2BQ4U3rCYgvo6oh\nwHu+/PT0dI4ePVrLLAqCIJwa1q/XaxkrBX/9K9x0U33nKDgB2lndWgA+7Q5Vr6UAx9LT041ZyYIg\nCPXKli0waBA4HDBuXGgNIw0kQDtNLRkZysag0G63U1xcjMcTWJkQBEE4fWzZAv36aXcTgwbBK6/U\n7yzjmgjQTlMCGtLGID4+HgDpNxAEob7YsAH69IHMTBgwABYtqh/X1LUhQDsbvDHIbdasGQBZWVn1\nnBVBEJoi778PvXtDbi4MHgwffQQNYT5sgHaWV5vYSygbg0PNmzcH4NixY/WcFUEQmhJKwZNPwo03\ngtOp+wgWL24YhgAgQDtNNWjVk8dtU2QZnkuzs7PrOSuCIDQVsrLgD3+AJUt0v8Azz8ADD4R2H0Eg\nAdppaqWzUDYGecaKPbL0pSAIp4PFi+FPf9LNQna7nljm9aXZoAjQzjp3VHe6KTK8nBYVFdVzVgRB\naMzs2AHXX68XpcnNhf794YcfGqYhgMp1kL3aWedrIJ9uCowecakZCIJwKsjKgnvvhfPOg//8B2w2\nmDkTPvkEtAPlhkmAdjb4ZqLimJgYrFarGANBEOqUrCx44QUt/CUluj/gD3/QaxefcUZ95+7XE6Cd\npjwnhbQxsFgsgT42BEEQTprNm7UBmDtXjxIC3RT05JPQrVv95q0uCdDOSDPHhLIxcAFERUX5L9Ig\nCIJQK/bu1fMF3nkHvv++Mn7wYPjLX6Bnz3rL2inFTztN6XyoGAM72vOpvw+NMoDY2FhKS0vrJVOC\nIDQ83G49a/h//9OTxLZsqfwuMRFuuQXuuQfOOaf+8ng68NPOBlEzGAg8C3RBi/9bwP1oL3su0Osl\nlJSU1FsGBUEIXRwO2LoVfvxRb7ds0QvN+Lcsx8fDNddoD6ODBjWciWO/Fj/tjKspLdSvMRgLzAI9\n/Ck2NjbKarXeAfQGLsE7ay4yMhKXy5TTPUEQGhEejx7mmZUFBw/Cnj067N+vPxtbFWQUfceOenjo\nkCHar1CkqXfjxoWfdoZ0zSAJmKWUYtq0aUyfPp0uXbqwcOFCWrdufS7aSDwMEBERgdvtrqdsCk0R\npbQQeTxQEbCUeFiYDuHhDWtGan2gFBQWQlGRDoWF+o09L0+LfF4e5OdDQYEORrrsbDh6VHsIDSb0\n/litWvi7dYOuXaFzZ+jRA1q0OD33GMr4aacpt3r+xkBh0odFHXAbwNdff81nn33Gc889R5cuXbjp\nppv47LPPiIuL+x3wHMC6detOU5bqBqX0KAXjYc7N1X8Ah0M/6KWllX+IwkL9J8jLq0zjcIDLpc9R\nVqb33W4tTMYfw2rVXhNjYnSIitIhMVFXiWNjdUhKgoQEPZPSbq9MGxenv7fZdHxqqk4XXt+NhrVE\nKT0s8NixynIsLNRlXlJSKS4lJTqUluptbq6OLy7WW6cTystPbABOhNWqyywiQjc9xMbqMk5MrCx/\nmw2Sk6FZM/3bpKXp/dRUvW3eXB8b6obFeEvPydHllp8Phw7pt/bCQl2umZn689GjcPiw/j3KTblI\nOzF2uy6zFi2gXTsd2rSBM8+EVq10aIpv/WYIDw+nXP8AId2BPBjg9ddfZ/LkyVx11VWAnjWXlZVl\nrKucDpCRkcHUqVNPeKLevacwYEAGSUn6obHbtdjFx+s/ZnS0flgiI/WfNjxchzDvdDuldKio0KJb\nXq5F2BBnl0sLTXFx5TYvT/8pDh/Wf4rCQr09ckTHmxWTk8UQLadT56WuSErSoWVLLVTG5+RkXZax\nsdpoREfrMrbZtHExyjcmRpdxRERlvNWqhc6/nD0eXcalpbqMDZEuLtaCk59fKeilpXrf+JyTUylI\nxcX6XKeCsDCd97CwSqE28m/cgxGM5+VkiYioNMiJifo5NsrebtflHBenPxuGJjZWx/kbo8jI4Pk1\nnmujnF0uHUpK9LNbVlZpJDMz9XOcmVn5Vm887zW9pQfDZtP/xYQEHWw2fQ8pKZUvKYmJOiQk6LSG\nkUxJaXgvKKeTmrSxVatWtG/fHkxOLq6vorYDbN26lR49egCwc+dOXC4XZ511lpHG1Eyz1at1CCWi\novSfOzVV/4ENw5SQUPknNt7aExP1w2+zVb5dRkbqfUNQIyIqhclownC7K4XUqEkYVe3SUh2MN2Wj\nKu50Hi++eXm6Wl5YqPfz8nS7bEMhOhrS06uWY0qK3tpsuoyNmpBRGzLe0v1fGsLDdRkboaY3daW0\nwJaXV/0tHI6qzSGGAcvN1fFZWbq8s7N1fGam/v0OH9YhVLFY9DObmlop8GeeqQ2XUcbGM9+8uX6T\nb/0Tc24AABofSURBVNYs9BaKb0pYLBZULSx4fRkDN2jPekuXLqV58+ZkZGTw8ssvG9+/AZjq8+/T\nBy65pPKtpqBA/ymN6r/xdl9WVvnnDVZ1NdqBjbdam62yWSU+vvLtJj5eC09ycuUDb8Q1b67/HKfy\nbcZi0XmNiNDilpxcN+ctL698Iz90qOobem6uLsvS0kqjYhgT4y2zrKxqs5axbzRvWSyVebdadbka\ntQ3DQBpG0m7X92XEGeUdH18pSIaQ19ciIxZLZS3IMPQni8OhjURRkS7zrKzKNnXDsBjfFRTo9Eat\nyTBGTqcu74qKyjd4/2clPFznMy6u8iXDaNIyXkKaNdOG9Iwz9LNsvM0nJOjn2mrKqYEQKlhq2fZo\nUVVNx+lquZwCZGzfvp177rmHlJQUZsyYYbhdPQR0BoYCc/r06YPFYuHzzz+v0wz4/2EEQRAaGwHa\neSKl8+l/fdUMngN+17Fjx27Lly/3j18PXAcU4W1K8ng8RJ6CHiIxAoIgNGb8tDOkXViXABcDjwJb\ngG+B29HzCzK9aeJB35BV6qeCIAi1wk87TQ2zqM++eifwtDcEIxmgvLyccBlSIAiCUCv8tLPBr4Fs\nB3C5XKekmUgQBKEx46edplw4hLIxiAVwu91E1NeQEUEQhAaKn3aacuEQysYgHqCsrIwoGawsCIJQ\nK/y009QaAKFsDGwATqeT6KbiZlAQBKGO8NNOp5n0oWwM4gFKS0uJjY2t77wIgiA0KPy002EmfSgb\ng1iA4uJibDZbfedFEAShQeGnnUVm0oeyMYhxuVyUl5eLMRAEQagFAdppanWwUDYGsXlel5x2u72e\nsyIIgtBwCNDOfDPHhLIxiMnM1JORU1NT6zkrgiAIDYcA7cysNrGXUDUGFiCmoKAAgKSkpPrNjSAI\nQgMiQDsbdM0gBrA6nXpElMwzEARBME+AdpoaWlpfLqxroiVwMCMjwxfhv99YaMz315jvDeT+GjpN\n7P7uB144QVKf/oeqMbgI2OC/OENtVuxpKDTm+2vM9wZyfw2dJnZ/twKzT5D0pNYzCAOGA52AzcBH\n3hOdA7QBfgCyvGmjgZ7oadAbqHSh2tp7/DbgoDfO6k1rRa9nUAak1SJfgiAIwonJNpOoNsZgFjDa\n7/NWIAfo7Rf3JLADmAkkeuP2Af8H/AG4wS/tm8CHwKtAC29cITAKrysKQRAE4VdzzEwis81E7YDd\nBQUFvPbaa4wbN45k7+K7xcXF/PDDD1x66aWEhVX2R2/duhWbzUabNm18cS6XizVr1tCrV68qnkh/\n+eUXnE4nnTt3NqL+DfyxiVXl6jEndU9jvjeQ+2voNLH7Oxf46QRJfTdf02iis4E+6KYg8vLyeOSR\nR7jrrrsAWLNmDW3btmXIkCEMHz4cp9OJUoq77rqLLl26cNlll7F48WIAdu7cSadOnbj++usZMGCA\nbxzs9OnTOeecc7jiiiuYOXOmcd0/1vLeBUEQhOAUmklUXTPRk8Bj/hFnnXUWAwYMYNmyZQBMnjyZ\nV155hYsuuogPP/yQSZMmMX78eH744QdWrVrFhRdeyMCBA+nYsSMzZszg3nvvZcCAAezYsYNx48bx\n1ltvMXfuXJYtW0bPnj258cYb6dSpEwMGDPBds6FbbbfbTXFxMQ6H4//bu/boKKtr/5vJJJnJPDOZ\nBAwkkkCWErgCiheNsAjVkFAtcMXbRSmPuBAErDbLuxC8vCL4jwqI1gveKqACLi/IVVrbhoRCLYqy\ntDSAhmAErgElyeQ1z8wjmd/948z5mAkhiS2SSZvfWnvle0723mefvfc533nA5XLB4/HA6/Vi4sSJ\nUc+tWLECbW1taGtrQyAQgM/ng9/vRyAQQDAYREdHB0KhEABArVYjPj4eWq0Wer0eiYmJSEhIgMlk\ngslkQlJSEvR6Pcxms0IWiwV6vV65l5SUFJU9/COAJJxOJ1paWuDxeOB2uxWdezwetLW1wel0orW1\nVbnndruVqfsdHR3o6OiIOpaQulKpVDh8+HDU/509ezaSkpKg0+lgMBhgMBhgMplgNpthNBphs9lg\nsVhgsViQmpoKnU7X73VPEs3NzWhqaoLT6YTT6YTD4UBraysaGxvR2toKl8sFr9er2LO07fb2dsWW\nZf1Wq9WIi4uDWq2GRqNBQkIC4uPjo0in08FoNCI5ORkpKSmKHScnJ8NqtcJgMECv18NoNP5D74ES\nDAZht9uRnp5+zWdIRtpYr5ajuFYwmALgP/1+P6qrq3HzzTdDrVbDZDLhJz/5CcrLy3H8+HH4fD7M\nmjULKpUKhYWFKC0txc6dO1FSUoL8/HwAwJgxY1BXV4djx45h69at0Gg0yMrKwoYNG7B3717MnTsX\n9957LwAgPz9faTFIlJaW4plnnrmmAPPmzcPixYthNpuRlpaG5OTk674zGknFmUujb2xsRF1dHZqb\nm+FyudDa2or6+nrU19fD7XajsbERdrsdLlfXa0R1DnJbtmyBTqeDTqdDQkICtFottFotEhISoNFo\nEBcXh7i4OJBEe3s7PB4PfD4fvF4v/H4//H4/nE6nMr64N9Dr9bBarbDZbEoFs9lsSE9Ph9VqhcVi\nUf5qtVoYjUYYDAZotVoYDIbrvjd1IBBAa2srmpqa4HA44PV64fF40NjYCIfDAbfbDbvdrjhzt9uN\npqYmNDU1oaWlBU6ns9fJQ1JSkiJPZx1HHssKRRIkFScWiRMnTsDr9cLn88HtdsPv7375eJ1OB5vN\nBpvNhtTUVCVIGAwGJCcnK3acnJwMk8mkBJfrbdskFdtxuVxoampCQ0MDCgoKop57+OGH0dDQgMbG\nRiXQNjU1dStnXFwcDAYDkpKSkJCQAJ1OF2XParU6qls5EAgoCU97e7uSBLW3tyMYDCIQCKCtrQ1u\ntzsqSF8LZrMZqampim1bLBaYTCbo9Xps2rQp6tlPPvkEVqsVer0eer0eBoPhBw0mHR0dSgD1er1o\na2uDx+OBw+FQEpS6ujrU1dXB4XDA4XDAbrcrZSB9yrp167r1jRHozimoEO4qulYwuAUAKioqUFpa\nikAgwLS0NNV7772H5ORkkMSZM2eQlZUFlUqFUCiENWvWYPHixdixYwdmz54NAPjqq6/w17/+FatX\nr4bZbFb2Ml61ahUWLlyI2tpajB07FgBgt9uxd+9elJeXK0xMmDChRyl37dqFXbt2RV3TaDSKIWq1\nWiQmJiI+Pv6qSq9SqZQKLo1QGp7M4P1+f68M0GAwIC0tDYMHD4bFYsHw4cOVY5PJpGSN0uA+++wz\n/OIXv0BiYiK0Wi2effbZHmXtDYLBoOJEpSFFZsuyQsnMWVbwpqYmnD59GvX19cq6Jj1BVu6EhARF\nDpnFpaenKzqePHmy4kxlxi0zRVkRpK57gtFoVLJAg8GA1NRU5ObmwmKxKK0gWbGNRiOMRqPSGpKZ\npdls/rsD2bp165TjzuPU29vblaTB6XTCbrfD4XAoerbb7YrjbWpqwoULF9DU1ASXy4VgsPtNqaRj\nNRgM0Ol0iIuLU67Fx8dDrVYrASwUCiktnWAwqJR/IBBAIBCA2+1Ge3v32+NqNBpUVFQgLS0NNpsN\nWVlZ0Ov1SElJwZAhQ2Cz2WA0GpWWkMVigc1mg8Fg+EFaP7L119zcrNhNc3MzWlpa4Ha74fV6lWSt\noaEBzc3NsNvtqKmpgdPphNvtxubNm6N+My8v76r/k5SUpOhY2rdsgUu7lhTpQ2RA6+joQCAQgN/v\nh8/nU6itra1HnUukpKTAbDbDZDLBZrNhxIgRsNlssFqtSEtLQ11dXbfvT548WSblvZpcfK1g8EcA\ngQceeCDhgQceALr4sJyZmYmPP/4Yu3btQllZGSZNmoSpU6fi6NGjeOWVVzBhwgQcOHAAu3btQkpK\nCurr6/Haa6+hsrISer0eixcvxptvvomdO3fC4/Fg3759ePnll2G1WgExSinl008/7XFCSHFxMebM\nmYOWlparskaZrfn9fsXJyy4AmeGpVColU9FoNErQkE3QxMREJTMzmUyK4dtsNgwaNAipqanQ6/VK\noItACGKYrB8iMgcg1hX3hslx5513NgNohvjab4fo23MAcEMsO+uAaOLJ9wMQBauBWOLbDCAFYvSV\nAcDg+Ph4m9lsNpnNZlN6enpq+BkLxP4QWgCJYUqIoCjP6Pf70draqjT5HQ6HkvW6XC7l2OPxKF1Z\n0uiDwSCCwSDS09MRCASiurYAMSNSo9EgMTFRaQlJXZtMJlgsFiUrlhVSnsssvgt0QGz6LfXtitC1\nK6xD+dcB4Luw3hvDuveGdRwM6zhS35FNAan7hNLS0qQIfeoBJAMwAUjVaDQWq9U6yGq1Jod1nxa+\nZwqXU1L4vasikszQpYOTDszlcinnMoBKxyKz5mAwGNVykXqWdi0Dogzg0q6ljUsnk56e3pNdh8K6\nkfbqDR+3QtTdy2FqCuvZHX7OGaZAuLwivaK0RU1Yr3qIlQgSw3/1AMwqlcpmNpuHms3mtLAujQCs\nYd3K52XZxKMLRxgMBpUWpcvlQnNzs+I3ZJIkz2XwlC0o2YKJbIFG+hDZ1RUXF6cED9nK12q1SgtJ\nBk+ZsBoMBpjN5qgkp4tNvaTefQDcpaWlQzs/EIn8/HyUlpY6cB1mIE8CsAqALXx+BwDs378fixYt\nwunTp1FbW4vt27dj0aJFShYfCATw0ksvwel0Yvny5TCZTACAM2fOYPPmzZg1axaKiooAiCi/Y8cO\nnDp1Ck8//TQGDx4MiIrthCjYQwCWQxRqGoBHAYyBaNacg6i80slJp6eDMKp4iMqmDsv1fdIU4orR\nd4T/j3TqHghH0whh8BcBtISPLwA4ie776KYBWAhgMIBPAWyAqCgAIJ13NoC7wsflAA5+D97/VpgB\nDAUwGsBwADdBVDLpyKTTk5VMC1FxpY7j0IVz6wGROpbOvA1XnIc8lk67BULftRDBsx7CAf3QH5Zu\nAvAfEGVSB2ALgI/C9xIh9GYEkAcx6OISxIi4nj7cqSB0mwMxVycFYi7OTRC2YMEVJ5eE6CCiwfe3\nbalvIlrn8RCOVQXhvP8M4CyEk/9XAMMgAugWAH/o5f+KFWQD2AUxl6kRwAcAaiDKbBiE3adD6FoD\n4Ww1EHqWwUSF7+9DAKFnSR0QupdB0IcrSYgfwmfI4OqCsKEqAN9C2Nw3uLad2wCUQvjo7wA8A+AU\nhB9cAGAGRPn+EcBzYRk78wkwGuiGRrFnHOni2sckL3VxvbyLa1XXuDaOZH0X9zaSVHXBayySiuSL\nXcjgITmW5J4u7knsJqmOARm6IjXJV7vguZnkLSTnXUOmB2OA997Q3SQdXfC/muRMki3XkK+e5L/E\nAP890V1d8F5DcgzJy13ce5n9p86B5J+7kGEFyS+6uC4xNwb47i2NI9nYif9Wktkk93ch29ckh3b6\nDTLyIOK8O1pCspaicpwg+QaFE1tJ0hR+5g6S/01yJ8n7wtcSSC4NP/ssyUHh6zkkXyK5i+S/k1xD\nkps3b+ZNN93EixcvRjH30Ucfcdq0aVy/fj1DoZC8vPw6K/eHogUk6fV6uXbtWk6dOpUnT56Mks/l\ncrGmpobV1dV8/fXXuWHDBrrdbnl7dgzI0BUtJUm3282nn36aRUVFPHPmTJRcoVCIO3bs4OTJk/n+\n++9H3rozBvjvjvQknSRZUVHBgoICbtq0KdL2SJLnz59nTU0NDx48yFWrVvEvf/mLvHU0BmToiVaQ\n5CuvvMJBgwbxwoULUbJ98skn/PGPf8w1a9awo6NDXl4VA3z3huJIsr29nUajkevXr4+Sra6ujjU1\nNTx+/Diff/55vv3227JsfSRTYoD/niiTZIAkP//8c06cOJFlZWVSvBBJNjY28oknnuCMGTMi/Wk1\nyfiI3yEjDyLObzTdR/L3JP9A0k+Sn332GQFw6dKlCmPnzp1jRkYGp02bxieffJLz58+PrJTaGCiY\nzjSS5P+Q/CPJNyWjS5cu5V133cWVK1dy1KhRrKysVGRctmwZ1Wo1rVYrCwoKWFhYyHfffVfefjIG\nZOpMKpLfkOT8+fOZl5fHVatWceTIkVEB4a233uLIkSO5cOFCFhUVcfv27fJWeQzI0B0tIcnKykpm\nZmZyxowZXLJkCZ944gnF9r788ksCoMFg4Lhx41hYWMgFCxZI+b6OARl6otkkeerUKQLgwoULlXKr\nra1lRkYGCwsL+dRTT/FnP/tZZEDQxwDvvaFLJDl16lTGxcWxvv5KB8PQoUOZmJjIrKwsFhYWsqio\nKNJhjogB3nui50nyt7/9LTMzM/nQQw/xkUceYSSKioo4ceJErl+/nrm5uZHyLYj4HTLyIOL8RtK9\nJBXrkgiFQrzllls4ffp05drKlSsjnQjz8vIis5isPuC9OxpK0VSLgtPp5IgRIxRHcuTIEc6ZM0e5\nP2HCBL700kudXyNF9B8VA3J1pjtI0m63Mzc3V2H2d7/7XZRR3n333Tx37hxJ0TLKzs6Wty7GgAzd\n0RGSXLJkCffvFy3uUCjE3NxctrSI3qG33nqLubm5DAaDV5ea6Bbsaxl6okSSHaFQiKNHj2ZhYaHC\n/Lp167ht2zblPD8/n9XV1fL01hjgvTf0Ikm++eabBMDTp0+TJC9fvky1Wq3YZSecYnTmHKv0Ikke\nOHCAX375JV999VWuWbNGEeKrr77iPffco5xv376dq1evlqfPRvwOSbKv9zN4DIB648aNKC4uxmOP\nPRZ6+eWXoVKpcMcdd+CLL75QHjx48CBmzZoFAHA6ncqIlTB6tRDTDcQ8AOYPPvgADz30EEpKSrBy\n5Ur86U9/QmFhoTLkrqqqCrfffrvy0siRI/HLX/4SmZmZysS+MKZBrAUVa8gGgA8//BD333+/cjFS\nLjkSKTs7GwBQXV2NcePGyUd7tWZKH2IEABw6dAgzZ84EADQ2NirD/QBRZlVVVdDpdCgpKYHX65Xv\nfgHgqT7g+fvCD+CgSqXC+PHjo+pcWVmZUufkKJthw4bJ273aPSsG8B4AjB8/HgAU+VJSUpCamorh\nw4ejoKAA58+fj3znPohBI7GOrQAwffp0DBkyBBs3bsSKFSuUm+Xl5XjwwQeV807+5iqf+X0Wqvsh\n4ASAqVOn4tZbbwUA9fDhwwFcmfF59uxZtLa2QqVSobW1FYmJiXjkkUewfPlyOczwQ4iv77EEJwDc\ndtttKC4uBiCMT87QBICjR49iz549OHLkCA4fPoycnBxs27YNP/3pT9HS0oIFCxagrKxMOs7eDUy+\n8SAApWwA4Tg/+OADHDp0COXl5Rg2bBg8Hg/a29tx+fJlLFu2DLt375bv/28f8d1bKKM3HA4H4uPj\nUVxcjHXr1sHtduPw4cOYOXMmTpw4gUuXLuGNN95ASUkJfv3rXwNiLZj+4FCAsGOXde7rr79GQ0OD\nUq4mkwmLFi1CSUmJ3CzlOMSoo/4AOxC9Vs/u3bsxZ84cfP7556isrMTJkycxZcoUfPPNN4Cou/0l\n0F2WB2vXrsXjjz8OvV4Pkti9e7dSfiSxb98+1NTU4IUXXpCvvH/Vr3VqHt3oZk42w33OnfHzn/+c\n2dnZ3LNnDw8cOMB33nmHmZmZHD16NA8dOhT56Og+4Lsn0lGMooqCz+fjlClTOGLECM6fP58Ohxik\nsnbtWp49e1Z5rqOjgzk5OayoqJCXCmJApq7oVlJ8PL7nnnuYk5PDRYsW0eVykRRde+fPn+eTTz7J\njIwMTp48ObJrrzasp76WoTv6DUm+9tprzMjI4JgxY3js2DGS5IkTJ/jcc89Fle+LL77ImTNnytP3\nY4D/3tIOknz44Yc5ZMgQ7t27l++++y7379/PzMxMjho1KvLDJClGv/U1z72lkSRZVVVFANy9ezfn\nzZvH9vZ2RZiqqirqdDp52hYDPPeW5pDiG6tWq+WGDRu4detW1tbWcsmSJWxoaODYsWOZk5PD5cuX\ns62tTcq4odPvkJEHEec3mpIoPiI/EKZqkiwrK+O+ffuimLPb7fR6vZGXpsRAgVyL1CQnRshVQ4o+\n586jpEgx4iE/P5/5+fmcNGkSb7vtNvp8Pnn75hiQ51q0sxu5lMK6ePFi5Af/b3n18LZYpNslw3V1\ndfT7/Z3l45YtW5iXl8cHH3yQJpMp0mk+HwP895Z+RYoRU++8806UfI2NjfR4PJGX7uvF78USZZJk\nS0sLN23axO+++46kGCWVl5fH6dOnMyMjg8uXL5fynYoBnntLc0ny+PHjnDt3Ljdu3Mji4mLW1tYq\nhdXR0cFLl6JG9m/i1cPUSYrRILG209l4AB9DTJiIhBNiAg4A/Abie8Ml9B9oADwNYB3EpCEfgBKI\nLqDXAeDbb79FRUUFLl++jEcffVTOxn4DwMN9wnHvkABgDYDV4XM3gMchJkpNh5gIIyetAcALEDpo\nu7Fs/s3IA/BfAMaGz98BsA2iexI+nw/Hjh3D0aNH8aMf/QiTJk0CxES42yAmyPUHjIHo+um82Xhk\nnfs9xL4ktTeQr+sBFUSXyPTIi6FQCJWVlfjwww9hsVhQXFwsu5JmIfa7LyVMEN2R11qxzgcxORQQ\n3+cehdiUrDNEd28MBgNAfLgrCh+XAfg6fDwM4oPX5S7e6S8wQsw+vogr08THANgHMRM1Eu8CmAsh\nc6zDDDFb+RtcPcMxHqLs7BCOsr9BBSALItDJ/mQNRFDovNx6A4D7AXx+w7i7PsiC4FsFMev9bPj6\nzRDfPr7rI76uBzQA/g1ABoD/g1g65BVccZSAcIjLIDbb6k8YAjE7fgjEjGYdRNlthpjNbA3TBVzZ\ncbIzYjoY/DNCBZGFToIo0I8BVOCHX2phAH8fMgAUQIys+g7A2+ifAe+fDXoAhRAtuABEMlbTpxz1\nHQaCwQAGMIABDEAEg76eZzCAAQxgAAOIAfT1PIO/Bf2526Q/trz6s74l+ove+7Ou+4uOr4X+pPsf\nRNedu4kGMIABDGAA/4QY6CYawAAGMIAB4P8BUvZjS7FYViQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110d42610>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "plt.xkcd()\n",
    "%matplotlib inline\n",
    "x = np.linspace(-1, 2, 128)\n",
    "p = np.polynomial.Chebyshev((0, 0, 0, 0, 0, 0, 0, 0, 0, 1), (-1, 1)) #These are Chebyshev series, a proto of \"chebfun system\" in MATLAB\n",
    "plt.plot(x, p(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Convergence of the Chebyshev-accelerated Richardson iteration\n",
    "\n",
    "Given the roots of the polynomials $x_k$, the best parameters are determined as \n",
    "\n",
    "$$\\tau_k = \\frac{1}{x_k}.$$\n",
    "\n",
    "The convergence (we only give the result without the proof) is now given by\n",
    "\n",
    "$$\n",
    "   e_{k+1} \\leq c q^k, \\quad q = \\frac{\\sqrt{c}-1}{\\sqrt{c}+1},\n",
    "$$\n",
    "\n",
    "where $c = \\mathrm{cond}(A)$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Beyond Chebyshev\n",
    "\n",
    "We have made an important assumption about the spectra: it is contained within an interval over the real line (and we need to know the bounds)\n",
    "\n",
    "If the spectra is contained within **two intervals**, and we know the bounds, we can also put the optimization problem \n",
    "\n",
    "for the **optimal polynomial**.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Spectra of the matrix contain on multiple segments\n",
    "\n",
    "For the case of **two segments** the best polynomial is given by **Zolotarev polynomials** (expressed in terms of elliptic functions)\n",
    "\n",
    "For the case of **more than two segments** the best polynomial can be expressed in terms of **hyperelliptic functions**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "<img width=100% src='WAT.jpeg'> </img>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## How can we make it better\n",
    "\n",
    "The implementation of the Chebyshev acceleration requires the knowledge of the spectra.\n",
    "\n",
    "It only stores the **previous vector** $x_k$ and computes the new correction vector\n",
    "\n",
    "$$r_k = A x_k - f.$$\n",
    "\n",
    "It belongs to the class of **two-term** iterative methods.\n",
    "\n",
    "It appears that if we **abandon** this and **store more vectors**, then we can go without the spectra estimation (and better convergence in practice)!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Crucial point\n",
    "\n",
    "The Chebyshev method produces the approximation of the form\n",
    "\n",
    "$$x_{k+1} = p(A) f,$$\n",
    "\n",
    "i.e. it lies in the the **Krylov subspace** of the matrix  which is defined as\n",
    "\n",
    "$$\n",
    "   K(A, f) = \\mathrm{Span}(f, Af, A^2 f, \\ldots, )\n",
    "$$\n",
    "\n",
    "The most natural approach then is to find the vector in this **linear subspace** that minimizes \n",
    "\n",
    "certain **norm of the error**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## What to minimize in the Krylov subspace\n",
    "\n",
    "The ideal way is to minimize\n",
    "\n",
    "$$\\Vert x_k - x_* \\Vert,$$\n",
    "\n",
    "over the Krylov subspace $K_n(A, f)$\n",
    "\n",
    "where $x_*$ is the **true solution**, but we do not know it!\n",
    "\n",
    "Then we can minimize the **residual**: \n",
    "\n",
    "$$\\Vert f - Ax \\Vert,$$\n",
    "\n",
    "For the SPD case we can do better!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Symmetric positive definite case\n",
    "\n",
    "Consider again $A = A^* > 0$ \n",
    "\n",
    "Then we can introduce **A-norm** as \n",
    "\n",
    "$$\\Vert x \\Vert_A =  \\sqrt{(Ax, x)}.$$\n",
    "\n",
    "Let us verify that it is the norm (on the whiteboard)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Conjugate gradient\n",
    "\n",
    "The most popular iterative method today (for positive definite matrices) is the **conjugate gradient method**:\n",
    "\n",
    "It minimizes the **A-norm** of the error over the Krylov subspace\n",
    "\n",
    "$$x_{k+1} = \\arg_{x \\in K_n(A, f)} \\min \\Vert x - x_* \\Vert_A$$\n",
    "\n",
    "And we can compute the minimization without knowing the error!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Why we can compute A-norm of the error\n",
    "\n",
    "Given $x$ we can compute the A-norm of the error:\n",
    "\n",
    "$$\\Vert x - x_* \\Vert_A^2 =  (A (x - x_*), (x - x_*))$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## How to get the CG formulas\n",
    "\n",
    "Consider the one step:\n",
    "\n",
    "$x_i = x_0 + y_i, \\quad y_i = \\arg \\min \\Vert x_0 + y - z \\Vert_A, \\quad y \\in K_n$.  \n",
    "\n",
    "Using Pythagoreus theorem,\n",
    "\n",
    "$(x_n, y) = (z, y)_A$ for all $y \\in K_n$, $\\rightarrow$ $r_n = Ax_n - f \\perp K_n$\n",
    "\n",
    "After some machinery..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## CG-method: formulas\n",
    "\n",
    "We get the following **three-term** recurrence relations for the CG method:\n",
    "$$\n",
    "   \\alpha_n = (r_{n-1}, r_{n-1})/(Ap_n, p_n),\n",
    "$$\n",
    "$$\n",
    "   x_n = x_{n-1} + \\alpha_n p_n, \n",
    "$$\n",
    "$$\n",
    "   r_n = r_{n-1} - \\alpha_n A p_n,\n",
    "$$\n",
    "$$\n",
    "   \\beta_n = (r_n, r_n)/(r_{n-1},r_{n-1}),\n",
    "$$\n",
    "$$\n",
    "    p_{n+1} = r_n + \\beta_n p_n\n",
    "$$\n",
    "\n",
    "\n",
    "The vectors $p_i$ constitute an A-orthogonal basis in $K_n$.\n",
    "\n",
    "We store $x$, $p$, $r$ (three vectors)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## CG-method: history\n",
    "CG-method has a long history:\n",
    "\n",
    "1. In was proposed in by [Hestenes and Stiefel in 1952](http://nvlpubs.nist.gov/nistpubs/jres/049/jresv49n6p409_A1b.pdf).\n",
    "2. In exact arithmetics it should give **exact solution** at $n$ iterations\n",
    "3. In floating-point arithmetics it did not - people thought it is **unstable** compared to Gaussian elimination\n",
    "4. Only decades later it was realized that it is a wonderful **iterative method**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## CG-method: properties\n",
    "\n",
    "1. Convergence estimate (same convergence speed as Chebyshev)\n",
    "  $$\n",
    "     \\Vert e_k \\Vert \\leq 2 \\Big( \\frac{\\sqrt{c} - 1}{\\sqrt{c} + 1} \\Big)^k \\Vert e_0 \\Vert.\n",
    "  $$\n",
    "  \n",
    "2. CG mystery (why it is much better than Chebyshev) is **clustering of eigenvalues:**\n",
    "\n",
    "if there are $k$ \"outliers\", then in $k$ iterations they are \"removed\" (we can illustrate it on the board)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Non-symmetric matrices\n",
    "\n",
    "CG only works for SPD matrices, and (sometimes) for symmetric matrices (although we can divide by zero).\n",
    "\n",
    "It completely **does not work** for the non-symmetric matrices.\n",
    "\n",
    "For non-symmetric matrices the two main methods are\n",
    "\n",
    "1. GMRES (Generalized minimal residual method)\n",
    "2. BiCGStab (BiConjugate Gradient Stabilized)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## GMRES\n",
    "\n",
    "GMRES idea is  very simple: **we construct orthogonal basis in the Krylov subspace.** \n",
    "\n",
    "Given the basis $q_1, \\ldots, q_n$ we compute the residual $r_ n = Ax_n - f$ and orthogonalize it to the basis.\n",
    "\n",
    "A Python realization is available as ```scipy.sparse.linalg.gmres``` (although quite buggy)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Disadvantage of GMRES\n",
    "The main disadvantage of GMRES: we have to store all the vectors, so the memory costs grows with each step.  \n",
    "\n",
    "We can do <font color='red' size=50.0> restarts </font> (i.e. get a new residual and a new Krylov subspace): we find some approximate solution $x$ \n",
    "\n",
    "and now solve the linear system for the correction:\n",
    "\n",
    "$$A(x + e) = f, \\quad Ae = f - Ax,$$\n",
    "\n",
    "and generate the new **Krylov subspace**.\n",
    "\n",
    "\n",
    "**BiCGStab method** avoids that using \"short recurrences\" like in the CG method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Idea of BiCGStab\n",
    "\n",
    "Idea of BiCG method is to use the normal equations:\n",
    "\n",
    "$$A^* A x = A^* f,$$\n",
    "\n",
    "and apply the CG method to it.\n",
    "\n",
    "The condition number has squared, thus we need **stabilization**.\n",
    "\n",
    "The stabilization idea proposed by Van der Vorst et al. gives the most practically used iterative method for non-symmetric systems.\n",
    "\n",
    "Let us do some demo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x113611d50>"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAESCAYAAADqoDJEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4lOW5+PHvPTMhK5CwBYGwGFDE5RJlFRE4VkCpRcSK\nKFZt1R45xeXUI5uaWKVqbWtXtUcBlVbh59Yq0YpHWcUFRKhY0ULZERBJAgRISHL//nhnhkmYhNmS\nzGTuz3Xlmpl35n3mmTcwd5772URVMcYYY0LlauoKGGOMSSwWOIwxxoTFAocxxpiwWOAwxhgTFgsc\nxhhjwmKBwxhjTFgscBhjjAmLBQ5jjDFh8TR1BeoiIpnAE0A5sERVX2jiKhljjCG+WxxXAv9PVW8F\nvtfUlTHGGONo1MAhInNEZI+IfFbr+GgR2SAi/xKRqd7DnYHt3vtVjVlPY4wxdWvsFsdcYHTgARFx\nA3/wHu8DTBSRM4AdQJ73ZfHcMjLGmKTSqF/IqrocKK51eACwUVW3qOoxYD4wFngVGC8iTwCvN2Y9\njTHG1C0eOscDU1LgtDQGquph4IdNUyVjjDF1iYfAEfG67iJia8IbY0wEVFUiPTce+g52crwvA+/9\nHaGfvhSYCRQAVwP9aN26NSkpKfTs2ZNBgwYxaNAg8vPzuf7661m3bh1Tp05l/vz5LFy4EFVtFj8F\nBQVNXod4+bFrYdfCrkX9P9GKhxbHaqCXiHQHdgETgImhnz4TuAr4BjgDOEpp6Vbc7kNs3ryZf//7\n36SkpJCdnU1VVRWjRo1i8ODB7Ny5kxEjRjBnzhy++uorhg4dypgxY2L92YwxJm4sWbKEJUuWRF1O\nowYOEXkRGAa0FZHtwP2qOldEfgK8DbiB2ar6ReilzgKeA07xPj4bUKqqFgOHACgvL6e0tJT58+eT\nnp7O4sWL6dmzJ6tWraK8vJyhQ4eSl5fHtGnTOO2009i2bRuFhYUx+tTGGBMfhg8fzvDhw3nggQei\nKqdRA4eqBm1JqOpbwFuRlfo80BEYCSzC+UjpQC9E1tKpUycqKyvZs2cPAGVlZbhcLtasWUNaWhoA\nR44cYdSoUXTp0oXXXnuNl156yd8S6du3L1lZWXHfGhk+fHhTVyFu2LU4zq7FcXYtYkdike9qKk7n\n+JnALcBeYBTHg8ennHLKFlq2PMrXX39NVVUVFRUVVFZW1iijbdu2fPvtt6SnpwNw5plncvDgQYqL\nixkyZAiZmZmMGDECwFJaxpiE5ktVPfDAA2gUnePNIHCMB3rjtDh8KSsPUInH8yUdOnxGt25tANi3\nbx/bt2+nsrKS1NRU0tPT2bdvn788t9uNy+WiqqqK9PR0fwf7wYMH/Smtu+++mxdeeMFSWsbUQyTi\n7yQTY8G+40UkuQPHqadez7//fSrOQKxcarY6vqB9+/088shEFi9eTEVFBatWrSI1NZWsrCw+//xz\nwOkD6d+/Px999FGN8l0uZ9CZL6V12WWXsWLFCrp06cKBAwd46aWXWL16dUKltIxpDN4vpqauRtKr\n6/eQ9IFj4cKlXHnl76io8ACTcfrYA4PHWiZO7MbYsRewZs0aevfuzfnnn8/3v/99WrduzY4dOxgy\nZAhvvvkmIkJ5eTktWrTg8OHDNd4rOzubkpISS2kZEwILHPHBAkcQIqIFBQVs2VLCc8/t4HjK6m2c\n0VaO9PT/5KWXrmXMmIsAKCwspFu3bmzYsIHrrruOxx57jPfff5/U1FRatmzJV199RWVlJeXl5aSm\nppKZmcnevXv95aWkpABYSsuYOljgiA+1fw/Wx4ETOHz1z8//gTdlVQk8BCzjeKujkr59d7NmzTMn\nlFFUVMSyZcvo3bs3qsrixYspKytj5cqV5OTknJDS6tevHx9//HGNMiylZUxNFjjiQ0O1OOJh5nhM\n/O53N9OixXqcQLEMp6P8uPXrD1JUtOyE88aMGcOjjz7KTTfdRG5uLt/97nfp1asXDz/8MAsWLODA\ngQOcddZZZGdnc+WVV7J+/XoyMzPxeDxkZGQAUF1dTXV1NYcPH6ZFixa8/PLLlJaW8vnnn9OqVSuu\nvvpqpk+fzldffcXChQvZs2cPc+bMYdq0aRQVFTX8xTHGnGD+/PkMHDiQrKwscnNzGTRoEE8++SQA\nN954Iy6Xi9dfr7m+6l133YXL5eK555zvl2effRa3203Lli1p3bo155xzDq+99pr/9Vu2bMHlctGy\nZcsaPy+99BIAO3bsYPz48bRv357s7GzOPvtsf9lxramnvkc5bV4DFRT8UV2uKxR+pDBDQQN+Zmjf\nvj/ScBQUFOicOXP0nnvu0XXr1umkSZO0R48e2rt3b+3fv7+2bt1aMzMz1ePxaGZmpnbo0EFx1t5S\nQFNSUjQlJUVdLpdmZmZqdna29uvXT08//XTt3r27Xn/99bpu3TqdOnWqzp49WwsKCsKqnzHxqvb/\nzXjzy1/+UnNzc/WVV17RQ4cOqarqp59+qpMmTdLy8nK94YYb9PTTT9fx48f7zzl27Jh26tRJe/Xq\npc8995yqqs6dO1eHDh2qqqrV1dX6pz/9SdPS0nT//v2qqrp582YVEa2qqgpaj+HDh+tdd92lhw8f\n1qqqKv3000/1rbfeitnnrOv34D0e8Xdvs0lV+RQWPsEDD7wFvEHtdFVm5qccOhTZX/ixSGm53W4A\nUlNTAUtpmearrhRJUVERQ4YMITs723+spKSE999/P6R/69GeD1BaWkrnzp2ZN28e48aNC/qam266\niXbt2jFv3jw2bNhAdnY2Cxcu5IknnuDgwYPccsst/OAHP+DZZ59l9uzZLF++HIDDhw+TlZXFxx9/\nTL9+/diyZQunnnoqlZWV/pR2oJYtW/L+++9zzjnnhFT3cFmqKkSFhZPJzGyJEzTexunvKAQe4vDh\nU4Kmq0IRi5RWVVUVVVVVHD58mNTUVEtpmaQzZMgQZs6cSUlJCeB86c+cOZMhQ4Y0yvkAH3zwAeXl\n5YwdO7be16WlpTF27Fjmz58PwPPPP88PfvCDOl9fVVXF3Llzyc7O5vTTT6/xXF1/oA8aNIjJkyez\nYMECtm3bFvJnaHLRNFea+oc6mmF9+96mMNObolrqvV+gMDPsdFUowk1ptW/fvkZKy+12q8fjsZSW\naTbq+r+pqlpcXKyTJ0/WzZs36+TJk7W4uDissqM9f968edqxY8caxwYPHqzZ2dmanp6uy5Yt0xtv\nvFHvvfdeXbFihQ4ePFhLSko0NzdXjxw5ohdeeGGNVJXH49Hs7GxNSUnR9PR0XbFihb9cX6oqOzu7\nxs+GDRv8n2XatGl65plnqtvt1nPPPVdXrVoV1uepT12/ByxVFaw5vIxx4/6XY8dupfbQ3LS023j5\n5Yn+obmxFm5Ka8CAAXz44Yc1yrBRWibRnWxU1ZYtW+jRo0fU77N582a6d+8e1jlvvfUW3/ve9ygv\nLz8hfZSXl8ef//xnnn32WfLy8vjZz35Gr169uOKKKygpKeHpp59m6NChQVNVZWVl/OhHP+Lw4cP+\nTvWTpaoCffvtt9x9992888477NgRxs4S9bBUVR0KCwtPWCZ4zJiLOOusVjj9G6OAe3HSVfdy9OhE\nfv/7dxqsPuGmtP7xj3/YKC2TVEpKSnjsscfYvHkzkydPpri4OKy/douLi5k8eTKbN2/mscce86et\nQjV48GBSU1P561//GtLrJ02axK9//et601QAmZmZPPnkkyxdupSlS5eGVSdw1s376U9/yq5duygu\nrr3DdmwsWbIkNnPKommuNPUP9TSHFy5cqikplwcdXXXmmbfWeV5DinaUlsfjsVFaJiHU9X/Tl2by\npZdqPz6ZaM/3+cUvfqG5ubn68ssv64EDB/wjmnJycnTJkiX+VJWq6v79+/W9997zn1s7VXXhhRfW\nKPvuu+/W0aNHq+rxVFVlZWXQetxzzz26fv16PXbsmB44cEAnT56sp512WlifpT51/R6IMlXV5F/+\ndVYMegDPAC/V85p6L1rLluNqBQ3np23bCfWe1xgWLlyo99xzj86ZM0dnz56tkyZN0nHjxmlubq72\n7t1b+/Xrp+np6Zqenq4ul0sHDhxYI4gA6nK51OVyaUZGhmZkZOhVV12lHTt21H79+ulpp52m69at\n09mzZ+vUqVN1/vz5unDhwqb+2CZJ1PV/c+HChSd8yRcXF4f8bzPa8wP95S9/0QEDBmhGRoa2b99e\nBw4cqE8//bRWVFTojTfeqPfdd1/Q8wIDx7PPPusfjuuzY8cOTU1N1XXr1vkDR1ZWVo2fxx9/XFVV\np0yZor169dKsrCxt3769Xn755f7+j1hoqMAR930cIvKSqn6/jue0vvqfffadrF//G2oPy+3WbTNb\ntvylQeobiaKiIg4dOmRraZlmw2aOx4eEXatKROYAY4C9qnp2wPHRwG9wdv17RlUfreP8iAPHqFH3\nsmjRydeuijfRrqXl8XgQEVtLyzQZCxzxIZEDx1CcPVyf9wUOEXEDXwLfAXYCq3D2Ge8HnAc8pqq7\nvK+NOHAUFS3j+9//I0eOLDjhuVGj7uPvf38wmo/WKMIdpWXLw5t4YIEjPiRs4AAQke7AGwGBYzBQ\noKqjvY+nAajqIwHntAF+DlxMHS2SkwUOqDtddeaZ37B+/Z+i/3CNqKFSWh6Phw0bNjBkyBBycnIY\nMWIEmZmZlJWVWWvERMQCR3xoboHjKmCUqt7ifTwJGKiqU8IsVwsKCvyPfRuxB0rUdFUoYpXSqq6u\nJjc3l+LiYjp37ky7du0A2L9/v7VGTEQscMQH3+/Bt5y6T0Isqx4kcIwHRscicJys/s0hXRWKWKS0\nMjIyOHr0KGlpaaiqdbCbiFngiA/NrcUxCCgMSFVNB6rr6iCvp9yTBg4ITFfVNGxYIUuWFIbzlgkh\nFikt26zKRMMCR3xoqMDhiapWkVsN9PIGlF3ABJzO8bAVFhYGTVEF6tQpi/XrTzyellYVyVvGPV8r\nYMKECYBzjaZNm1YjpZWbm1sjpZWVlUV1dTXHjh2jf//+rFy50l9eWVkZLpeLNWvW+DvYjxw5wqhR\no+jSpQuvvfYaL730EnPmzLGUljFxrHbKKlKNMarqRWAY0BbYC9yvqnNF5FKOD8edraoPR1B2SC2O\noqJl3Hzzc+zefQq+zvGOHXfxzDM3JnQfRySCpbSOHTvGv//9bzZu3Eh+fj7bt2/n0KFDdbZGWrdu\nTWlpadA5I/n5+eTl5TFq1CgACyRJyloc8SGhU1UNJbzA8Vd27/61/1jHjv/NM89ckXSBI5AvpQWw\nZs0arrvuOmbNmsXy5cvp3r07lZWVJ+1gd7vd/g729PR0MjIyEBFatWpFZWUlAwcOJCUlxd83Yimt\n5JCogeO2226jc+fO3HvvvU1dlTp1796d2bNnc/HFF5/0tQ0VOJp8aZFofgAtKCjQxYsXa31GjvQt\nsV7zZ9Soe+s9L9nEYhkUt9utLVq0ULfbrWlpadqzZ08dNGiQDho0yJZBSSLE8Q6A3bp10/T0dM3K\nytKcnBwdM2aMbt++PehrS0tL9Y477tCuXbtqVlaW5ufn65133qn79u3zv+bFF1/UAQMG+NeXGzhw\noD7xxBP+57dv365XXnmltmvXTlu3bq1nnXWWPvvss6p68h0Cg+nevbu+++67Ib229u9h8eLFWlBQ\n0HzXqgqp8iH+4xw2rCBo4Bg2rCCk85PRwoULdf78+f5Asm7dOj3ttNO0f//+esopp+hVV12lGRkZ\n/kUZMzIyTggigIqIP9j4FmTs0KGDjhs3TidNmqSzZ8/2BxILIs1HPAeOwC/eo0eP6g9/+EO94oor\nTnhdeXm59uvXT0eOHKlffPGFqqru3btXH3roIX3zzTdVte4taK+77jqtqKhQ1fq3hz3ZIognq//J\n1PV7sMARAmtxRC+clX19AaV2S6Suzao6dOigEyZMqLGyrwWSxFbf/82FC5fqyJEzddiwAh05cqYu\nXLg0rLKjPb/2F29RUZF/RdobbrjBvyru008/rbm5uVpWVha0nJKSEs3MzNRXX3213vfLysrSdevW\nBX0uLy+vxiKIH374oW7cuFFHjBihbdu21Xbt2ul1112nJSUlNer/8MMPa58+fTQnJ0dvuukmPXr0\naNDyGypwNNWoqkZ1++0j2bRpJps2jcI3ezw9/QsGDRrW1FVLGIH9EkVFRXTq1In77rsPVaeDvUuX\nLqxcuZJOnTpx4MABqqurAaioqKBNmzbs2bPHf37tUVoej4dPPvmESy65hK5du/Lyyy9zwQUXkJeX\nx7Rp06xzvRkpKlrGHXe8zaZNxyfjbto0EyCk/sZoz/dxvjudPcIXLFjA4MGDASf3L+Kk/v/v//6P\nSy+91L9PTm2hbkHr2x52ypQpDB48mK5du/qfW758OT169KC0tNS/NNCmTZuYOXMmF110EaWlpYwf\nP57CwkIef/xxf91feOEFFi1aREZGBpdffjkPPfQQDz7YiHPSook6Tf1DiH0cqqoFBX/U9PQf12hx\n5OfPCPuvFVNT7ZTW7NmzdcKECTpo0CA95ZRT9Nprr9Xc3FxNTU2tM6XVqlUr/zLxqampmpubq9de\ne6127NjRUloJijr+0o229R+L7EG3bt00KyvLv91r586d9bPPPlNVrbGc+iWXXKLTp0+vs5xQtqBV\nrX972FD6OF577TXt27ev/3H37t31T3/6k//xm2++qfn5+UHPrf17sD4ODT1VpWrpqsYSGEh8Ka38\n/HwdPHhw0M2qcnNzTwgkmZmZltJKcHX936yrvxHqOh7a68LprwxMVVVXV+urr76qbdq00d27d9cI\nHNdcc43ecMMNdZbz5ptvqsfjCfql36VLF1269MQ/Svft26c33nijdu7cWVWDB47du3frhAkTtHPn\nztqqVSvNysrSrl271qi/r49FVXX9+vWanp4etI51/R6iDRwJv3VsqMrLg2fljh51N3JNmrcxY8Yw\nYcIEHn30UbZv306nTp2YOXMmN998M6effjr/8R//QVZWFj179uSMM87gwIEDpKenIyJ07tyZ9u3b\nU1ZWRnV1NWVlZZSWlvLJJ5+wbds2jh496k9pvfLKKzz44IPs2rXLn9KaM2eODfWNc6mplUGPjxpV\nFVLYGDky+PmRTuYVEcaNG4fb7WbFihU1nvvOd77D22+/fcI8Jp9wt6CFE7eH9aXFAs2YMQO32836\n9espLS1l3rx5/tSvz7Zt22rc79SpU8h1iIWkCRx1/YNtrrPH40Go+6+3atWKnj17UllZyZEjR8jI\nyMDjcQK97y+cI0eO4HK52LhxI99++y27du1i7969/hns7777Lo8++ihXXnmlf//1BQsW2B7sceb2\n20eSnz+zxrH8/BlMmXJJo5zv4/zR7dz+7W9/o6SkhD59+gRmM7j++uvJy8tj/PjxfPnll1RXV/Pt\nt9/y85//nLfeeovs7GwKCgqYPHkyr7zyCgcPHqS6upq1a9dSVlbmf6+pU6fy+eefU1lZycGDB3ny\nySfp1asXOTk5tG/fHpfLxaZNm/yvP3ToEJmZmbRq1YqdO3fy2GOPnVD3P/7xj+zcuZP9+/cza9Ys\nrrnmmrA+f9Siaa409Q9hpKoWLlyq+fk19x/Pz59ufRxNJHCUlm+4b35+vvbu3VsHDx6svXr10rS0\ntHpTWu3bt1fghOG+3bt31wkTJtToG7H91xtXff83Fy5cqqNG3avDhhXoqFH3RjSqKprzu3fv7p/H\n0bJlSz377LP1hRdeUFU9YcvY0tJSvfPOOzUvL88/j+OnP/2p7t+/3/+a+ragVT359rD333+/tm/f\nXnNycvSjjz7Szz//XM8//3zNysrSvn376q9+9SvNy8urUf9HHnlE+/Tpo9nZ2XrjjTfqkSNHgn7W\nun4PNPetY+sT6sxxn6KiZUyd+g779rk599wqpky5JKlnjscT3xLxqs4orYqKClatWkVqamqNlX2P\nHj1Kp06dqKio4JtvvvGf36JFCwAqKytp0aIFXbp0oV27dhQXF1NZWcmrr77qXx7eFmVseIk6c7y5\nsSVHggg3cAC88QY89RRYBiM+1bey77Zt22jVqhUHDhygrKyMqqoqysvL6devHx9//PEJZfmG+mZm\nZtKyZUsOHDhA165dbdfDRmCBIz5Y4AgiksDxq18t42c/W0Tfvh5SUyu5/faR1uqIY4GbVfkCyZVX\nXonL5eLQoUMnXSLe5XKRnp5OWVlZ0EUZbZ+RhmGBIz5Y4AjCtwPgyZZV9ykqWsbkyW+zbdvxyUP5\n+TP57W9HWfBIEKHuelhRUUHbtm3ZvXt3jfPr2mfE4/HwzTffcMEFF/j3GbGWSOQscMSH2r8H37Lq\nCbEDYEMJt8XhbCP7UJDjzWcnwGRS166HH374ISNGjODdd9+lpKSEyspKXC4X/fr144MPPqhRhtvt\nRlVJSUkhOzubiy++mPfee48ePXrY8vBRsMARH5KyxSEiY4ExQCucPTveqfV8WIFj+PBCli4tPOF4\nc90JMJkE9o34WiIffPABbdq0Ydu2bWRlZbF3716qq6vr3GckNTWViooKMjIyaiwPf+DAAUtphckC\nR3xIysDhIyLZwC9V9eZax63FYU4Q2BLxpbRmzZrF0qVLycnJqZHS8u0zkpGRUWOUlsvlwu12B01p\nbdmyhWHDhjFjxgx/Sguw1kgACxzxIaEDh4jMwWk57FXvvuPe46M5vgvgM1rHnuMi8kvgz6q6ttbx\nsIfj1l4gLT9/Br/97Wjr42jG6kpprVy5kpycHP9wXxGhoqKCysoTJ4v6FqBLS0ujZcuW5OXlsXHj\nRnr27Mmpp55KixYtrDUSwAJHfEj0wDEUOAQ87wscIuIGvgS+A+wEVuHsO94POA94DPgaeARYpKrv\nBik37FFVRUXLuO66d+jRw01urs3lSDZ1Dfdt27YtAPv27WP79u3++SAiUmMWMEB2djYHDhwgNTUV\nt9uNx+OhZ8+eHDx4kPLycoYOHervYE/WOSPBltIwTSNhAweAiHQH3ggIHIOBAlUd7X08DUBVHwk4\n53bgBzhBZa2q/qlWmWEHDoBLL4UpU+CyyyL7LKb5qGviocfj8S8Pf/jwYX9KKysrq8YS8eB0sIPT\nRwJw2WWXsWLFCrp06WJzRkxcSuTAcRUwSlVv8T6eBAxU1SlhlKkFBQX+x6EOy504ES6/HK69NpxP\nYJqz2i0RVWXRokVs3bqVzZs3B01pdejQgV27dtUop2PHjuzevdvmjJi44huG65Mww3GDBI7xwOho\nA0ck9b/tNjj7bJg8OexTTZIIJaW1c+dOcnJy2LNnD6mpqaSlpfHtt9/6y6hrzoiltExTS+QWxyCg\nMCBVNR2orquDvI4yIwoc06dDq1bOrTGh8KW0MjIyKCsrY/HixRQXF7Ny5Up69+7NsWPH/Otp1bUM\nim/OSFpaGmApLdN0EjlweHA6xy8GdgEfAxNV9Yswygxr5jg4neN33bWII0c89OljS46Y8NW3ntaO\nHTtOugwKOB3sJSUlltIyjSqhZo6LyIvAMKAtsBe4X1XnisilHB+OO1tVHw6z3BgMx7UlR0x0Ql0G\npa45Iy1atEBVLaVlGk3CtDgagk0ANPEm1Dkj4KS0zjvvPFavXl2jDN+ckdTUVETEUlom5pI+cIST\nqrIlR0xjikVKy7fkiaW0TCwkVKqqoViLwySScFNateeMuN1uRITq6mpLaZmoJH2Lw5YcMYko3JRW\n//79+eijj2qUEbgMCtgoLRM6CxwRLDny85+/w9q1boYOtSVHTNOLRUqrZcuWHDx4kNTUVFwul6W0\nTL2SPnCEOxwXYNMmGDnSuTUm3kSb0hIR/8q+GRkZltIyftbHQeQTAL/+Gs47z7k1Jp6Fm9IaOHDg\nCZtViQgiYikt45f0LY5I6l9aCl27OrfGJIpYjtJKS0tDRGqktPLz823XwyRhgSOC+h87BhkZzq0x\niSrclFZmZiZ79+6tUYbb7aa6uvqEXQ8rKysZOHAgKSkp/r4RS2k1HxY4Iqx/SgocPuzcGpPoYpHS\nCtz1sEWLFnTp0oV27dpRXFxMZWUlr776qj+lZX0jic0CR4T1b90atm1zbo1pTmKR0vJJS0vD4/GQ\nmZlJy5YtOXDgAF27drW+kQRngSPC+p9yCqxZ49wa05yFk9IKtuuh2+0mLS2NsrIym8HeTFjgiLD+\n+fmwaJFza0yyqC+l5es49+16WFFRQZs2bU7Y8dDj8SAiJyzK6PF4+Oabb7jgggv8w32tJRKfkj5w\nRDKPA+Css+DFF50NnYxJRvXterh161ZGjBjBu+++S0lJCZWVlbhcLvr378/KlStrlOPbZyQlJYXs\n7Gwuvvhi3nvvPXr06GGjtOJMs5/HISK9gTtwlmJ/W1VnB3lNxC2OAQPgD39wbo0xNQOJL6X1wQcf\n0KZNG7Zt20ZWVhZ79+6lurq6zr6R1NRUKioqaozS6tq1Kxs3bmTIkCHk5OQwYsQIMjMzKSsrsw72\nJtLsWxwi4gLmq+rVQZ6LKHAUFS3j+usX0bWrh9xc28zJmNoCU1q+vpFZs2axdOlScnJyQtpnRETw\neDxUV1eTm5tLcXExnTt3pl27dgDs37/fOtibSNwHDhGZA4wB9vp2//MeH83xTZyeCbZlrIhcDkwG\nnlbVV4M8H9FaVbaZkzHhCXW4r4hQUVFBZWXlCWWkp6dz9OhR0tPTUVV/B7tveLDNGWk8iRA4hgKH\ngOcDto1142wb+x1gJ7AKmAj0A84DHlPVXQFl/E1VxwYpO+zAYUurGxOduob7tm3bFoB9+/axfft2\nKisrg47SAmfXQ8D/Gpsz0rjiPnBA0P3GBwMFqjra+3gagKo+EnDOMOBKIA34QlV/E6TcsAOHbeZk\nTGz5hvv6WiIVFRWsWrUKj8fjH6V15MgRqqurOXbsGOeffz4ffvjhCeWkp6fjdrttzkgjiDZweGJZ\nmTB0BrYHPN4BDAx8gaouBZbG+o1TU09sQgOkpVXF+q2MSQq+lkBRURHf/e53WbNmDaNHj/aP0tq5\ncyfl5eVs3LiRPn36sHnzZjIzM0/oYD9y5Ih/CZQ9e/aQnp7O559/zplnnsnVV1/tnzOyc+dORowY\nwZw5c2xI9kq1AAAadElEQVTOSBNpqsARs2ZOYPM1lGG5t98+kk2bZp6wmdOUKaNjVSVjkpLvy3vC\nhAmAE0jGjRsHwJo1a3jmmWeYNWsWO3fu5JxzzqGyspKvvvqKzMzMGnNGAgOJx+Nh7dq1/l0PFy9e\nTM+ePVm1apV/ifi8vDymTZtmKa16+IbhxkpTpaoGAYUBqarpQHWwDvKTlBvxqKpbb32HrCw3PXrY\nZk7GNIa6Otg//PDDkOeM2K6HsZGofRwenM7xi4FdwMfARFX9IsxyI54A+MMfwpAh8KMfhXWaMSYG\nYjFnJDs7m5KSElsGJQwJMwFQRF4EhuFM5NsL3K+qc0XkUo4Px52tqg9HUHbEEwBvvhkGDXJujTFN\nJ5I5I7WXiE/xLnNdexkU2/UwuIRocTSUaALHLbdA//5w660xrpQxJirhLhF//vnns2rVqhpl2K6H\n9bPAEWH9b70Vzj8ffvzjGFfKGBMzsVgivnXr1pSWllpKK0DSB45I+zh+/GPo2xf+8z8bpm7GmNiL\ndtdDt9uNiPhHaSVbSith+jgaUjQtjttug3POcW6NMYkn3JRW//79+eijj2qU4RullZqaiogkTUor\n6VsckdZ/8mQ480z4r/+KcaWMMY0uFimtrKwsDh06RFpaGiLSrFNaFjgirP9//ReccQb85CcxrpQx\npsmFm9LKysqqsWGVy+XC5XJRVVVFRkZGjZRWcXExI0aMYMaMGf6UVqIFkqQPHJH2cUyZAqed5twa\nY5qvcFNaAwYMOGEtrcBRWh6Phw4dOvjX0iouLk6YXQ+tj4PoWhy33+5sG3vHHTGulDEmbsVylJbL\n5UrYXQ+TvsURaf3vuAN69IA774xxpYwxCSPalBY4y5+Ul5fX2PWwVatWVFZWxu0+IxY4Iqz/nXdC\nt25w110xrpQxJiGFk9I6evQonTp14tixYzWG+7pcLtxuN1VVVXG9z4gFjgjr/9//DV26OLfGGBOo\nvpTWtm3baNWqFQcOHKCsrIzq6uo6dz2E+NxnxAJHhPX/6U/hlFPg7rtjXCljTLMTmNLyBZIrr7yS\nlJQUcnJyTrrrocvlIi0tjcOHD8fFDHYLHBHW/+67ITcX/ud/YlwpY0yzF8quh4cPH66xz0igumaw\nezwevvnmmwYfpdVgOwCKyE/rOU9V9deRvmksFRYWRjQc1+WCBI6ZxpgmdLJdD7du3crWrVv9+4yk\npqb69xkZMGAA77//vr+ssrIyXC4Xa9as8Y/Sys/PZ8SIEZx66qls3ryZUaNGsWDBAsrKyqJqjcRq\nQ6c6WxwiUkjwnfoEJ3A8EPW7RymaFsfUqZCTA9OmxbhSxpikFYt9RlJTUzl27BgdO3bk4MGDuN3u\nmK+n1axTVSKSCSzB2S2wKMjzEQeOadOgdWuYPj26OhpjTDCx2GekRYsWVFdXo6qkpqYCsVkivsED\nh4ikAz8C+gDpeFshqvrDSN805MqJPAAcBL6IdeCYPh1atoQZM6KspDHGhCDU4b4iUucorTZt2rB/\n//6oO9gbI3C8DHwBXAc8AEzC+SK/PaQ3EJkDjAH2+raO9R4fzfEdAJ+pvd+4iFwCtAHSgH2xDhwz\nZ0JGhnNrjDGNqa7hvm3btgVg37597Nixwx9ETjnlFMrLy09ojQBUVlaGvUR8YwSOtap6roj8Q1XP\nEZEUYIWqDgzpDUSGAoeA5wP2HHfj7Dn+HWAnsAqYCPQDzgMeAyYDmTgtnSPAuNpRIprAce+9kJoK\n990X0enGGBMztUdpARQXF7NixQrS0tLIy8vjyy+/9Ke1zjvvPFavXl2jDI/HE3JK69FHH23wwPGx\nqg4QkeU4X+a7gY9U9dSQ30SkO/BGQOAYDBSo6mjv42kAqvpIkHNvAL5R1TeDPBdx4LjvPkhJgfvv\nj+h0Y4yJOV9LBPC3Rnr16sU111xDp06dQl5Pq0OHDuzdu7dGSss31Ldfv34sWLCgYYbjBnhaRNoA\n9wKvA1lAtH+ndwa2BzzeAQRtwajqc1G+V1A2HNcYE28C+yMmTJgAOK2RBx98sMZ6Wrm5uTXW08rM\nzPR3sKenp/tTWkeOHMHtdrN27VpUlS5duhCLAVEnDRyq+rT37lKgR9Tv6C02RuXUGIYWznwOEaiu\njlUtjDGmYQR+xxUVFdGpUyfuu+8+f1qrS5curFy5kry8PH8He3p6etBdD7du3crXX38ddZ1CSVUV\nBDz0v1hVfxbym5yYqhqEM8TWl6qaDlTX7iAPodyIU1UPPABVVfCzkD+FMcbEj2iXiG/oVFUZxwNG\nOvBd4J+RvqHXaqCXN6DsAibgdI6HLdKZ4yKWqjLGJC5fWiswpTVt2rR6U1ppaWlUVVVx7NixqN47\n7AmAIpIKLFLVYSG+/kVgGNAW2Avcr6pzReRSjg/Hna2qD4dVEaJrcTz0EBw96twaY0xzEmzOyLFj\nx9i1axdnnnkmTz31VOPOHPd2lH+sqj0jfdNYiSZwzJoFZWXw85/HuFLGGBNHAkdqjRo1iuzs7IZb\n5NBHRD4LeOgCOgBx0zNgqSpjjKlb4EitBl/k0P8Cpx/CpxLYo6rRJchiJJoWx8MPQ2kpPHLCzBFj\njGneGnJZ9TbeuwdqPdXS+6b7I33TeOBy2XBcY4yJRH2pqjU4o6kE6AoUe4/nAFuJ3ZyOJmGpKmOM\niYyrridUtbuq9gDeAb6rqm1VtS3OgoXvNFYFG4oFDmOMiUydgSPA4MB1olT1LeCChqtSeAoLCyPq\n7LElR4wxyWbJkiVhb/oUTCid44uAZcCfcdJW1wIXqeqoqN89StF0jv/617B9Ozz+eIwrZYwxcS7a\nzvFQWhwTcYbgvga86r0f0SzveGKpKmOMiUwoixx+C4S0aVMiscBhjDGRqW847m9V9Q4ReSPI06qq\n32vAejU46+MwxpjI1NfieN57+6sgzyX8V64tq26MMZGpM3Co6ife2yW+Y95JgV1U9R8NX7WGZakq\nY4yJzEk7x0VkiYi08gaNT4BnRCRuxiLZcFxjjAlNYw7HXauq54rIzUCeqhaIyGe+TZmaUjTDcZ98\nEtatg6eeinGljDEmzjXGcFy3iJwCXA0UeY81+N/qIjJcRJaLyJMiEtLeH+GVby0OY4yJRCiB42fA\n28AmVf1YRPKBfzVstQCoBg4CqcCOWBdugcMYYyJz0sChqi+p6jmqepv38SZVHR/qG4jIHBHZU2tf\nD0RktIhsEJF/icjUIKcuV9XLgGnAA6G+X6isj8MYYyITSuf46SLyroh87n18jojcG8Z7zAVG1yrT\nDfzBe7wPMFFEzhCR60XkcRHpFNB5UYLT6ogpG45rjDGRCSVV9TQwA6jwPv6MMJYcUdXlHF+S3WcA\nsFFVt3g3hZoPjFXVeap6l6ruEpFxIvIUznyS34f6fqGyVJUxxkTmpEuOABmq+pGI0wGvqioi0e4A\n2BnYHvB4BzAw8AWq+hrO+lgNwlJVxhgTmVACxzci0tP3QESuAr6O8n1j9pUdOCY5nL3HLVVljEkW\nsdpr3CeUeRz5wP8Cg3H6GzYD16nqlpDfxNm3/A3f3A8RGQQUqupo7+PpQLWqPhpW5aOYx/Hcc/Du\nu/D88yd/rTHGNCcNtue4j6puAi4WkSyc/TgO4czp2BLpmwKrgV7egLILmECES7UXFhaG1dLwsT4O\nY0yyiVXLo84WhzdQ/BjIB9YDTwFjgVk4HdshrY4rIi8Cw4C2wF7gflWdKyKXAr8B3MBsVX047MpH\n0eL485/h7393bo0xJpk0ZIvjeeAA8AEwErgROApcq6prQ30DVQ3akvBuQftWyDWNMevjMMaYyNQX\nOHqq6jkAIvIMTod4N1U90ig1C5GlqowxJjSNkar6VFX71vU4HkSTqpo/H/76V+fWGGOSSUOmqs4R\nkYMBj9MDHquqtor0TeOBpaqMMSYy9W3k5G7MijQ2S1UZY0xkQllyJK5FupGTBQ5jTLJptI2c4lk0\nfRyvvAIvvODcGmNMMmmMjZyaJevjMMaYyCR14EjgxpYxxjSZpA0ctjquMcZEJmkDh6WqjDEmMkkd\nOKzFYYwx4bPAYYwxJiwJHzgincdhfRzGmGRj8ziIbh7HW2/Bb3/rLK1ujDHJpME3cmoq4mxy/hDQ\nElitqjHdq89SVcYYE5l4TlVdAXQGKoAdsS7cAocxxkSmwQOHiMwRkT0i8lmt46NFZIOI/EtEpgY5\n9TTgfVW9G7gt1vVyuWw4rjHGRKIxWhxzgdGBB0TEDfzBe7wPMFFEzhCR60XkcRHphNPKKPGeEvOv\neGtxGGNMZBq8j0NVl4tI91qHB+DsW74FQETmA2NV9RFgnvfYq8DvRWQosCTW9bLAYYwxkWmqzvHO\nwPaAxzuAgYEv8G5Re/PJCgocWhbOFrI2HNcYkyxitWWsT1MFjph9ZUc6JtmWHDHGJIvaf1Q/8MAD\nUZXXVKOqdgJ5AY/zaICRU/WxVJUxxkSmqQLHaqCXiHQXkRbABOD1SAqyHQCNMSY0CTNzXEReBIYB\nbYG9wP2qOldELgV+A7iB2ar6cARlRzxzfMUKmDbNuTXGmGQS9zPHVXViHcffAt5q6Pevi/VxGGNM\nZOJ55nhILFVljDGhSZhUVUOKJlX14Ydw553OrTHGJJNoU1UJ3+KIlKWqjDEmMkkdOBK4sWWMMU0m\n4QOH9XEYY0xorI+D6Po4PvkEbr3VuTXGmGRifRwRsj4OY4yJTFIHjgRubBljTJNJ2sBhq+MaY0xk\nkjZwWKrKGGMik9SBw1ocxhgTPgscxhhjwpLwgSPSeRzWx2GMSTY2j4Po5nF88QWMGwcbNsS4UsYY\nE+fifln1SInIhcB1OHXso6pDYlu+tTiMMSYScRs4VHUFsEJExgIfx7p8S1UZY0xkGryPQ0TmiMge\nEfms1vHRIrJBRP4lIlPrKeJa4IXY18uG4xpjTCQao3N8LjA68ICIuIE/eI/3ASaKyBkicr2IPC4i\nnbyv6wqUqmpZrCtlqSpjjIlMY2wdu1xEutc6PADYqKpbAERkPjBWVR8B5gW87ofAnIaolwUOY4yJ\nTFP1cXQGtgc83gEMrP0iVS08WUGBQ8uGDx/O8OHDQ6qA9XEYY5LFkiVLIpq2UJdGGY7rbXG8oapn\nex+PB0ar6i3ex5OAgao6JcxyIx6Ou2ULDBsGW7dGdLoxxiSsRF1WfSeQF/A4D6fV0WgsVWWMMZFp\nqsCxGuglIt1FpAUwAXg9koJs5rgxxoQmYWaOi8iLwDCgLbAXuF9V54rIpcBvADcwW1UfjqDsiFNV\nO3bAwIGwc2dEpxtjTMKKNlWVtEuO7NwJ/fvDrl0xrpQxxsS5RO3jiJlIU1XWx2GMSTYJk6pqSNG0\nOHbvhnPPdW6NMSaZJH2LI1K25IgxxkQmqQNHAje2jDGmySR84LA+DmOMCY31cRBdH8e338Jppzm3\nxhiTTKyPI0LWx2GMMZFJ6sCRwI0tY4xpMkkbOGzJEWOMiUzSBg5LVRljTGSSOnBYi8MYY8JngcMY\nY0xYEj5w2LLqxhgTmmY/j0NEugC/A4qBr1T10SCviXgex5EjkJMDR49GV09jjEk0zXkex9nAK6r6\nI6BvrAu3VJUxxkSmwQOHiMwRkT0i8lmt46NFZIOI/EtEpgY5dSVwq4i8C/w91vWyVJUxxkSmMVoc\nc4HRgQdExA38wXu8DzBRRM4QketF5HER6QTcBNyrqhcDY2JdKRuOa4wxkfE09Buo6nIR6V7r8ABg\no6puARCR+cBYVX0EmOc99h5wv4hcC2yOdb0sVWWMMZFp8MBRh87A9oDHO4CBgS9Q1X8AVzVUBSxw\nGGNMZJoqcMTsKztwaNnw4cMZPnx4SOdZH4cxJlksWbIkomkLdWmU4bjeVNUbqnq29/EgoFBVR3sf\nTweqgw25PUm5EQ/Hdc53+jkk4kFpxhiTeBJ1OO5qoJeIdBeRFsAE4PVICop0AqCPtTqMMckiYSYA\nisiLwDCgLbAXuF9V54rIpcBvADcwW1UfjqDsqFocbjdUVDi3xhiTLKJtccTtzPFQRBs4PB5nBnlK\nSgwrZYwxcS7awNFUneNNrqhoGVVVi7j4Yg/p6ZXcfvtIxoy5qKmrZYwxcS/hA0dhYWFYo6nACRp3\n3PE2MIvly51jmzbNBLDgYYxptmI1uiopU1WjRt3LokUPBTl+H3//+4OxqJoxxsStRB1V1aTKy4M3\ntI4etV5yY4w5maQMHKmplUGPp6VVNXJNjDEm8SR84IhkHsftt48kP39mjWP5+TOYMuWSGNbMGGPi\nS8LM42hI0QzHLSpaxu9//w5Hj7pJS6tiypRLrGPcGJMUbB5HAtffGGOagnWOG2OMaVQWOIwxxoTF\nAocxxpiwWOAwxhgTloQPHNEuq26MMcnChuNio6qMMSYSzXZUlYj0EZEFIvKEiIxv6voYY4xxxG3g\nAEYDv1fVycAPmroyxhhjHA0eOERkjojsEZHPah0fLSIbRORfIjI1yKnzgGtE5Bc4uweaelg/z3F2\nLY6za3GcXYvYaYwWx1yc1oOfiLiBP3iP9wEmisgZInK9iDwuIp1U9RtV/QkwHdjXCPVMaPaf4ji7\nFsfZtTjOrkXsNPhGTqq6XES61zo8ANioqlsARGQ+MFZVH8FpaSAi3YAZQCbwi4aupzHGmNA01Q6A\nnYHtAY93AAMDX6CqW4EfN2aljDHGnFyjDMf1tjjeUNWzvY/HA6NV9Rbv40nAQFWdEma5NhbXGGMi\nEM1w3KZqcewE8gIe5+G0OsISzQc3xhgTmaYajrsa6CUi3UWkBTABeL2J6mKMMSYMjTEc90VgJXCa\niGwXkZtUtRL4CfA28E9ggap+0dB1McYYE70GDxyqOlFVO6lqqqrmqepc7/G3VPV0Ve2pqg+HW24I\n80CalWDzYUSkjYi8IyJficgiEckOeG6699psEJGRTVPrhiEieSKyWEQ+F5H1InK793jSXQ8RSROR\nj0RkrYj8U0Qe9h5PumsBzlB/EflURN7wPk7W67BFRP7hvRYfe4/F7lqoasL9AG5gI9AdSAHWAmc0\ndb0a+DMPBfoCnwUc+wVwj/f+VOAR7/0+3muS4r1GGwFXU3+GGF6LjsC53vtZwJfAGUl8PTK8tx7g\nQ+DCJL4W/w38BXjd+zhZr8NmoE2tYzG7FvG85Eh9/PNAVPUYMB8Y28R1alCquhwornX4e8Bz3vvP\nAVd4748FXlTVY+rMldmIc82aBVXdraprvfcPAV/gDPFO1utx2Hu3Bc4fVcUk4bUQkS7AZcAzgG/g\nTNJdhwC1Bw/F7FokauAINg+kcxPVpSnlquoe7/09QK73fidqjlJrttfHO9S7L/ARSXo9RMQlImtx\nPvNiVf2c5LwWjwP/A1QHHEvG6wCgwP+JyGoRucV7LGbXoqmG40bL5m/Uoqp6knktze6aiUgW8Apw\nh6oeFDn+B1YyXQ9VrQbOFZHWwNsiMqLW883+WojId4G9qvqpiAwP9ppkuA4Bhqjq1yLSHnhHRDYE\nPhnttUjUFkdM5oE0A3tEpCOAiJwC7PUer319uniPNRsikoITNOap6l+9h5P2egCoailQBJxP8l2L\nC4Dvichm4EXgP0RkHsl3HQBQ1a+9t98Ar+GknmJ2LRI1cNg8EMfrwA3e+zcAfw04fo2ItBCRHkAv\n4OMmqF+DEKdpMRv4p6r+JuCppLseItLONzpGRNKBS4BPSbJroaoz1Bm12QO4BnhPVa8nya4DgIhk\niEhL7/1MYCTwGbG8Fk3d+x/FqIFLcUbTbASmN3V9GuHzvgjsAipw+nduAtoA/wd8BSwCsgNeP8N7\nbTYAo5q6/jG+Fhfi5LHX4nxJfoqz0nLSXQ/gbGCN91r8A/gf7/GkuxYBn28Yx0dVJd11AHp4/z2s\nBdb7vh9jeS0SeutYY4wxjS9RU1XGGGOaiAUOY4wxYbHAYYwxJiwWOIwxxoTFAocxxpiwWOAwxhgT\nFgscJiGJyCHvbTcRmRjjsmfUevx+LMuPNRG5UUR+39T1MMnDAodJVL4JSD2Aa8M5UUROtkbb9Bpv\npDoknPKbQFSTsUTEvgdMWOwfjEl0jwBDvRvW3OFdKfYxEflYRNaJyK0AIjJcRJaLyN9wZtMiIn/1\nrh663reCqIg8AqR7y5vnPeZr3Yi37M+8m+RcHVD2EhF5SUS+EJE/B6uo9zWPiLPx0pcicqH3eI0W\ng4gsFJGLfO8tIr/w1vEdERkkIktFZJOIXB5QvG9zq69E5P6AsiZ53+9TEXnKFyS85f7Su6ruoJj8\nJkzyaOrp8fZjP5H8AAe9t8OANwKO3wrM9N5PBVbhbE4zHDgEdAt4bY73Nh1nLZ+cwLKDvNd4nKUa\nBOgAbMXZVGo4UIKzPLXgbJU8JEidFwOPee9fCrzjvX8j8PuA170BXOS9X413CQjgVe/7u4FzgE8D\nzt8F5ABp3s9yPs7mVq8Dbu/rngCuDyj3qqb+PdpPYv4k6rLqxvjU3qxmJHC2iFzlfdwK6AlUAh+r\n6taA194hIr7NbPI4+eJuFwIvqKoCe0VkKdAfOOAtexeA96/47kCwvpFXvbdrvK85mQpVfdt7/zPg\nqKpWicj6WucvUtVi7/u/6q1rFU4AWe1dcj4d2O19fRXO6sLGhM0Ch2mOfqKq7wQe8O7RUFbr8cXA\nIFU9KiKLcf5ar49yYqDy9S+UBxyrou7/W+VBXlNJzbRxYD2OBdyvxlnkElWtrqevRgLq9Zyqzgjy\nmqPeAGhM2KyPwyS6g0DLgMdvA5N9X6oicpqIZAQ5rxVQ7A0avamZ5z9Wx5fycmCCtx+lPXARTgul\ndjAJ1xacjZhERPKIbAvTS0Qkx7u0+lhgBfAucJW3rohIGxHpGmVdjbEWh0lYvr+W1wFV3vTQXOB3\nOCmcNd59O/YC47yvD/wL++/Af4rIP3GW5/8g4Ln/Bf4hIp+os6eDAqjqayIy2PueirOE+V4ROYMT\nRzaF8te8r9wV4mxA9E+c/dM/qaccDXJfcQLYKzib8MxT1TUAInIvsMjbKX4MmAxsC7F+xgRly6ob\nY4wJi6WqjDHGhMUChzHGmLBY4DDGGBMWCxzGGGPCYoHDGGNMWCxwGGOMCYsFDmOMMWGxwGGMMSYs\n/x8T17t/ZCr1mgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1140a5c50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#### import numpy as np\n",
    "import scipy.sparse.linalg\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "n = 50\n",
    "ex = np.ones(n);\n",
    "lp1 = -sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr'); \n",
    "rhs = np.ones(n)\n",
    "ee = sp.sparse.eye(n)\n",
    "\n",
    "#lp2 = sp.kron(lp1, ee) + sp.kron(ee, lp1)\n",
    "#rhs = np.ones(n * n)\n",
    "res_all = []\n",
    "res_all_bicg = []\n",
    "def my_print(r):\n",
    "    res_all.append(r)\n",
    "\n",
    "def my_print2(x): #For BiCGStab they have another callback, please rewrite\n",
    "    res_all_bicg.append(np.linalg.norm(lp1.dot(x) - rhs))\n",
    "#res_all_bicg = np.array(res_all_bicg)\n",
    "sol = scipy.sparse.linalg.gmres(lp1, rhs, restart=10, callback=my_print)\n",
    "plt.semilogy(res_all, marker='x',color='k', label='GMRES')\n",
    "sol2 = scipy.sparse.linalg.bicgstab(lp1, rhs, x0=np.zeros(n), callback=my_print2)\n",
    "res_all_bicg = np.array(res_all_bicg)/res_all_bicg[0]\n",
    "\n",
    "plt.xlabel('Iteration number')\n",
    "plt.ylabel('Residual')\n",
    "plt.semilogy(res_all_bicg, label='BiCGStab', marker='o')\n",
    "plt.legend(loc='best')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Battling the condition number\n",
    "\n",
    "The condition number problem is **un-avoidable** if only the matrix-by-vector product is used.\n",
    "\n",
    "Thus we need an **army of preconditioners** to solve it.\n",
    "\n",
    "There are several **general purpose** preconditioners that we can use (short list today, more details tomorrow),  \n",
    "\n",
    "but often for a particular problem a special design is needed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Preconditioner: general concept\n",
    "\n",
    "The general concept of the preconditioner is simple:\n",
    "\n",
    "Given a linear system \n",
    "\n",
    "$$A x = f,$$\n",
    "\n",
    "we want to find the matrix $P$ such that \n",
    "\n",
    "1. We can easily solve $Py = g$ for any $g$\n",
    "2. Condition number of $AP^{-1}$ is better\n",
    "\n",
    "Then we solve for (right preconditioner)\n",
    "\n",
    "$$ AP^{-1} y = f.$$ \n",
    "\n",
    "or  (left preconditioner)\n",
    "\n",
    "$$ P^{-1} A y = f,$$ \n",
    "\n",
    "The best choice of course is $P = A,$ but this does not make life easier.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Other iterative methods as preconditioners\n",
    "There are other iterative methods that we have not mentioned. \n",
    "\n",
    "1. Jacobi method\n",
    "2. Gauss-Seidel method\n",
    "3. SSOR (Successive over-relaxation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Jacobi method\n",
    "\n",
    "Jacobi method is when you express the diagonal element:\n",
    "\n",
    "$$a_{ii} x_i = -\\sum_{i \\ne j} a_{ij} x_j + f_i$$\n",
    "\n",
    "and use this to iteratively update $x_i$.\n",
    "\n",
    "From the matrix viewpoint we use \n",
    "\n",
    "$$D = \\mathrm{diag}(A).$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Gauss-Seidel (as preconditioner)\n",
    "Another well-known method is **Gauss-Seidel method**. Given $A = A^{*} > 0$ we have  \n",
    "\n",
    "$$A = L + D + L^{*},$$\n",
    "\n",
    "where $D$ is the diagonal of $A$, $L$ is lower-triangular part, and Gauss-Seidel is $(L + D)^{-1}$.  \n",
    "\n",
    "**Good news: ** $\\rho(I + (L+D)^{-1} A) < 1, $\n",
    "\n",
    "where $\\rho$ is the spectral radius."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Successive overrelaxation (as preconditioner)\n",
    "\n",
    "We can even introduce a parameter into the preconditioner, giving a **SSOR** method:\n",
    "\n",
    "$$(D + w L) x = w b - (w U + (w-1) D) x,$$\n",
    "\n",
    "$$P = (D+wL)^{-1} (w U + (w-1) D).$$\n",
    "\n",
    "Optimal selection of $w$ is **not trivial**.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Others (but very important, see it tomorrow)\n",
    "\n",
    "- Incomplete LU for sparse matrices (tomorrow)\n",
    "- Algebraic multigrid\n",
    "- Sparse approximate inverse \n",
    "- Domain decomposition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Take home message\n",
    "- Iterative methods are the core\n",
    "- Krylov subspaces are the core\n",
    "- Condition number is the problem\n",
    "- Conjugate gradient is the best for SPD\n",
    "- GMRES/BiCGStab are the best for non-symmetric\n",
    "- Preconditioners are crucial (and often do not work)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Next lecture\n",
    "- More on preconditioners for sparse matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Questions?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>\n",
       "<style>\n",
       "    @font-face {\n",
       "        font-family: \"Computer Modern\";\n",
       "        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "    }\n",
       "    div.cell{\n",
       "        /*width:80%;*/\n",
       "        /*margin-left:auto !important;\n",
       "        margin-right:auto;*/\n",
       "    }\n",
       "    h1 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\n",
       "    h2 {\n",
       "        font-family: 'Fenix', serif;\n",
       "    }\n",
       "    h3{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "        margin-top:12px;\n",
       "        margin-bottom: 3px;\n",
       "       }\n",
       "\th4{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "       }\n",
       "    h5 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\t   \n",
       "    div.text_cell_render{\n",
       "        font-family: 'Alegreya Sans',Computer Modern, \"Helvetica Neue\", Arial, Helvetica, Geneva, sans-serif;\n",
       "        line-height: 1.2;\n",
       "        font-size: 120%;\n",
       "        /*width:70%;*/\n",
       "        /*margin-left:auto;*/\n",
       "        margin-right:auto;\n",
       "    }\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "    .text_cell_render h1 {\n",
       "        font-weight: 200;\n",
       "        font-size: 50pt;\n",
       "\t\tline-height: 110%;\n",
       "        color:#CD2305;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\t\n",
       "    .text_cell_render h5 {\n",
       "        font-weight: 300;\n",
       "        font-size: 16pt;\n",
       "        color: #CD2305;\n",
       "        font-style: italic;\n",
       "        margin-bottom: .5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\n",
       "    \n",
       "    li {\n",
       "        line-height: 110%;\n",
       "    }\n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "\n",
       "</style>\n",
       "\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"]\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"./styles/custom.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}