lectures/lecture13/lecture-13.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 13:  Iterative methods and preconditioners "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Previous part\n",
    "\n",
    "- Arnoldi orthogonalization of Krylov subspaces\n",
    "- Lanczos for the symmetric case\n",
    "- Energy functional and conjugate gradient method\n",
    "- Convergence analysis\n",
    "- Non-symmetric case: idea of GMRES"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Now we will cover the following topics\n",
    "\n",
    "- More iterative methods: MINRES, BiCG and BiCGStab \n",
    "- The concept of preconditioners"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## What methods to use\n",
    "\n",
    "- If a matrix is symmetric (Hermitian) positive definite, use CG method.\n",
    "- If a matrix is symmetric but indefinite, we can use MINRES method (GMRES applied to a symmetric system)\n",
    "- If a matrix is non-symmetric and not very big, use GMRES\n",
    "- If a matrix is non-symmetric and we can store limited amount of vectors, use either: GMRES with restarts, or BiCGStab (the latter of the product with $A^{\\top}$ is also available)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### More detailed flowchart from [this book](http://www.netlib.org/templates/templates.pdf)\n",
    "<img src=\"flowchart_iterative_methods.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## MINRES\n",
    "\n",
    "The MINRES method is GMRES applied to a symmetric system. We minimize\n",
    "\n",
    "$$\\Vert A Q_j x_j - f \\Vert_2 = \\Vert Q_j \\widehat{x}_j + h_{j, j-1} q_j \\widehat{x}_j - f \\Vert_2 = \\Vert Q_{j+1} \\widehat{H}_{j+1}  \\widehat{x}_j - f \\Vert_2 \\rightarrow \\min$$\n",
    "\n",
    "which is equivalent to a linear least squares with an **almost tridiagonal** matrix\n",
    "\n",
    "$$\\Vert \\widehat{H}_{j+1} x_{j} - \\gamma e_0 \\Vert_2 \\rightarrow \\min.$$\n",
    "\n",
    "- In a similar fashion, we can derive short-term recurrences.\n",
    "\n",
    "- A careful implementation of MINRES requires at most $5$ vectors to be stored."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Difference between MINRES and CG\n",
    "\n",
    "- MINRES minimizes $\\Vert Ax_k - f \\Vert_2$ over the Krylov subspace\n",
    "- CG minimize $(Ax, x) - 2(f, x)$ over the Krylov subspace\n",
    "- MINRES works for indefinite (i.e., non-positive definite) problems.\n",
    "\n",
    "CG stores less vectors ($3$ instead of $5$). \n",
    "\n",
    "Now, let us talk about non-symmetric systems."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Non-symmetric systems\n",
    "\n",
    "- The main disadvantage of GMRES: we have to store all the vectors, so the memory cost grows with each step.  \n",
    "\n",
    "- We can do  **restarts** (i.e. get a new residual and a new Krylov subspace): we find some approximate solution $x$ 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** from the residual vector. This spoils the convergence, as we will see from the demo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0xa1495f128>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAakAAAEcCAYAAAB9K0udAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3W2MY9l5H/j/U9Uvekn33OL0xBoHKfXcklc71kIekdVIZpHNTqdJOasPK8Qmu+0vXmSjJp0FDAjwhlTFkISRI7dY8sICEsBglQADAYSoi4xjZBMhcrHV2mSVCVDFmoFheazYRY0G61iGZoqcGUga9XTXsx/OPaxL8vL97d6q/w8gqury8txz61bx4TnnueeIqoKIiCiMlhZdASIiol4YpIiIKLQYpIiIKLQYpIiIKLQYpIiIKLQYpIiIKLQYpIiIKLTOLboCYSEiLoA0gAMAcQBbqtpcbK2IiM424c28hojUVDXhfe8A2FbVzIKrRUR0pp2qlpSIpAE0VbUa8FwcwDqAOgAXQN3u5z13ZPdV1aaIJIc55pUrV/Tq1atTqD0R0dlRq9VeU9UnBu13aoKUF1S2AXS1fryuvKKqpnzbyiJSV1UbtDq79o5EJK6qB/2Oe/XqVezv709+AkREZ4iIfG+Y/SKfOCEiroiUYALNUY/dcgBKHdtKAIre97EZVY+IiCYQ+SClqnVVzanqVp/dbEKE3763HTDBzel4noGLiGjBIh+kBvGSILpaWTZzz+sKrCMgKA3q6iMiotk69UEKXvDpk07udgYjL3B1JV8QEdF8nZrEiT46u/F6uS0ieZzcJ3W7144ikgWQBYDV1dWJK0hERMHOQktqKKp6oKqbqlr1vva8kVdVt1R1XVXXn3hiYAZlu81N4P799m3375vtRETU5iy0pACYsalQzCBx7Rpw8ya+/dnP4t++9Rb+10uX8KHnnwd2dhZdMyKi0DkLQarufY3Bdy+Ul1Dhf34+rl/Htz/7WfyNX/s1/G0Afx3At//5P8eHrl+fazWIiKLg1Hf3ea2nOoJTzJvezbxz9Vf//t/jVQDXAewA+LdvvcUuPyKiAKc+SHmqMFMi+cWxoAy+p37mZ/A/AHgE4H8H8D8fHAA3bwLnzjFQERH5nLYgFUNwNl8B3dMl5bzt83X/Pp76V/8Kf/J3/y6WAFwA8LcrFfzxz/4scOcOAxURkU/kg5SIOCJSFJEyTIAqikjJm2wWQKvLr+Dtl/ZSzYuL6OrD3h6ws4P/++//ffw/AJYBKIC1//gf8Z+fe+4kUOVyDFZEdOZxqY4Jra+v6zgTzH77X/wL/PVf+zU8APA3YALVAwD3rl7FR3/wA5w7dw74zGeA73wHWFsD8vkp15yIaHG85ZE6h2G6RL4lFUn37+NDzz+PP0uncQHAXwIQmK6//+WVV6A//CG+trKCh5/5DHD3LltWRHRmMUgtwt4esLGB//Gb38SfpdM4B+D/gwlUAnNfgA1Wf7q0hEe/8RvtwYoBi4jOCHb3TWjc7j5sbpqgc+cO/vNzz+FnKxW8B8B5mEAFmC5A65V3vxtPPniA8+fPY/niRXYFElGkDdvdxyA1obGDFHASqJ5/Hg+Pj/Ev3/9+/PK3v413eU8HBauHAHavXsXP/8VfYBkAPvEJ4NEjE6yuXQO++lUGLiIKPQapOZkoSAEmUB0eAh/8YCtY/eGVK0i+8grOe7uIb3d7td4BgKUlnD8+xiMRfO+JJ/DU0RGWLl4EPvc508qyGLSIKGQYpOZk4iBlDRGsjmFS1q2gK/dIBPBdU11awo+ffhqXn3oK+Omfbt+ZwYuIFoRBak6mFqSsjmD16MEDPHz4EP/t3Dlcffvt1m62dfVjAO/2vg+6kg9gsgYBM8OFpUtLeO3jH8eTnbO4M3AR0RwwSM3J1IOUZYMV0Baw/uzd74bbaLRaV4AJRMvoni2435W1wYuBi4gWgUFqTmYWpPwCugL/39VV/E8vv4wlVTwQwbFqK+EC6O4a9Ot1xbsClwiWlpexdPu2Sc7wY+AiogkwSM3JXIKUZYPVL/2Sudfq8BBYXga+/GUcq+Lh8TGWVAHVngGqV/AaOnABUBEcOg7e+847+P5734sfr6zgv/vgB/G+n/opswMDGBENwCA1J3MNUkH83YKAaW195jPAT34CADg+PoYcHwMwGYHnA4oYtdXlH+cCTABr7ecLYMNaXl6G4zh4z7vffbLxAx9gogfRKcYgNScLD1KdOoMWAPzFXwCvvw7Uaq3MPxu8xglcwOCW1zQ8AgARiAiOAbzkjZW974c/HOr133/ve3HkD3wAfvhTP4Unf+d38Oyzz06plkQ0DgapOQldkOolKHh5XYXDBq5HGG4p51n8Rdm8xnf13aubv5V3DKAGM1X+uy5cwNLysgmC/heI+clu+6u/9te6At333/te/MHP/MyINentfe97Hz7ykY/gxRdfxPe///2hX/Mrv/IrDLYUWQxScxKZIBVkQOA6Pj6GHh9jCb27Cq1hA9g0jPsXO+gceukV6CZxCOC/BWz74ghliAh+7ud+Do899hiWlpYgXqvTPjq3df5sy5i1cYLwqOUzYEcPg9ScRDpIBRmixYXjY/PwPFpexlJn9t+UvAPTqhknuPQyyV/8uIEuyAPv6xJOgt8SgB8AWJvSMYDpBMSwswH78uXLXcHaPt/rYZ9fWlrq+t7+PMisA/GizeKDAIPUiEQkDrPEvAPgGoDCMIsinrogFaTXONfhIeA4bWNdnQFsGD3Ht7w3kfMjltdp0PjaJCb973kEU7d3YFqiDzHdgOw/zrRbg/2chcB41iwtLeHWrVv40Y9+hEajARHBU089hWw2O1bwYpAagYg4AG6q6pb3cxJASVUHfqA9E0Gqn34BbBIf+ID5+vWvt03zhIcPRyvn/HnAl2k47l/7LANdL7P6z5xma7CfabYUGfTC6+LFi7h///7IgWrYIDWvYYSxecvAN1W1GvCcbf3UAbgA6kH7DcEFUACw5f28D8AVEcdbep56mWVK+OZmdxr6KAHQtvLOeX/mvnT8kZw/j+URUupHsrRkWp7vfz/wve+1PTWr0aJ+2ZfTDIwXYILTEoDjpSU8e3yMhyI4N8YH46cRfK/e//H22/jRj3/ctT+D2vw8ePAA3/zmN2c2JhjqIOW1aLYBZAKecwEUVTXl21YWkfow3XR+qnogIinfpnWYwMgAtUiTBsDNTeDDH27fNmorLyDQjdqd2ZNt5S0vmwC1vNw9s8c02YDYx7QDo219LnvHPT9Bz03bm9W5c/jvGw0AwYE1KKjtvv/9uNGx37QyNU/zmNTR0RG+9a1vmUSqgOt34cIFPPfcczM7fiiDlBeACjA9BEc9dssBKHVsKwEoIiCoDdIR2HIAbo9aBoXMNFp50wh0QRwHeOklIJE4GdN79Kire3JqZlXuIH1aihPxnUuvwNoZ1D72yiumPjYpQgT4xCfwv3HKr4FeeOEFfPOb38Tjjz/eFojnkVkZ+jEpETkEkOvsxvO2p/zBxRtbaqiqeD9n0b/7ezeg3CyAI1WtDFO/Mz8mRePZ3DSLVNrprQAT/H7wA+CJJyYPgn5TSG4Zi7+l+OjRdFqKQ7QGR65fQODiXJWzd2oSJ4KClA1GAFY6u+RERAGsjdrl5702CQCjjGsxSFHozSq5pR/bUnzmGRMgHz0ywWWSFt04rx01qPUKXIkEcOUKp+qaolOTONFDDAD6jBm5MMkUQ/OSMI5U9cD7OT1sa4oo1BbxJrq5CfzWb5mWou0unaSlOM7YYK+g1rEwaBu7v7/88+eB//JfzPdBra6/9beAN988mW9ybc20kr/6VQaxKYhqSyoOoGa79Tr2V5huwKFbQ94YWOetI/VeKehel2AWAFZXVxPfm2ZfOxF1G7U12KuLs1fgGqcb0ZblL/PcOWB93RwbMC2wpu+zNCdObjnt3X1TDVKTYHcfUQiNMnPKNAJXUOusX7mdXYnNgE6hUx7QTnt3HwAzNsU0cSLqEvRGvrkJ/MN/2L7NBq7ObsR+419BwSvow36v1/fqSuz0p38aHNB+7/fMz6c8iFlRDVJ2vCkGoBWkvIQK//NERMakgQvoHbzs9gsXzKwow7bABiWC9Ato/iB27hxw8aJZSy6XM8+fkoAVySClqk0RqaN7+rEYzE24DFJENNiwgcu/Jlu/VtfSkglUb7+NvsYZAwsKaLYODx8CTz0FfOpTJsieooAVySDlqcLMDHHg2xb3thMRjafXm3nQjd0f/CBw967J7rPJGgDwrncBDx4EB6JRU+n7ZSP6vfyy+frOOycBS/UkC9EGLCsigSsKiRMNALc708G9rr1yx7RIuzBJFnNrSTFxgogADJeBOOqN1cMEtEFBLKQ3LUc6u88LQBsw9zulYcaYqjAzRFR8+8UB3AKw5+17MK+sPotBioiGNkoq/TABrTOIDdvq8qfPX/CmHD53znQRPnw4l0AV6SAVJQxSRDQTgwKaP4gtLZlWkO3yA4YPWCLA5csmQN25A2xszCVQMUjNCYMUES2EP4h98IPA88+fjIN1BiyrV8LG8jLw2GPAb/zG3AIVg9ScMEgR0cL1Clg2HX6YeROXlkzrbE6BikFqThikiChUJplt4/x506La2DCBamcHuH59JtVkkJoTBikiCr1egev3fs/MVP/SS6bV9PAhcOmSuc9qZ8fst7c3k9bUsEFqaepHJiKicMnngVKp/XH1KvD5zwN//ufm63veY/Z96y3g8cfN9zdvmqy/zc2FVZ1BiojoLMrnTcvJdu199rMngeo73wE+9rGT5xYYqBikiIjOqs5A9e/+HfDcc+a5t982aekLDlQMUkREZ5kNVHYM6o//GPibf9N8/8MfAp/73EmgunZt7tVjkCIiOutsYsTNmyYg/fjHJh0dMPMSfv7zJ0Fszq0pBikiIjJZfLbFtLFhUtTtrBUXL5p9FpBIMVEKuog8A7M8RhdV/cbYBUcIU9CJ6NTY3DRByAaqT38a+NGPzHMXLpjt9v4pYKL09JneJyUivwhgx/4YsIuq6vLIBUcQgxQRnSr+QLWzY9LV7941z73nPSa5AjCtqglu9p31fVJFAN8FcBNAIuAx8MBERBRCnYkU9+4BH/qQ+f5HPwK+9rWTAAXMvOtv3CDlAsir6r9W1ReDHtOsJBERzVFnIsVf/dXJgo+//dvAr/7qyfMzzvgbN0hV0WMsioiITgF/IsXODvClL5lECgD4wheAf/AP5tKaGjdI5QDkROTXReQZEbna+ZheFedPREqLrgMR0UJ1dvvdvAkUi+b7hw/NLOsvvjjz1tS4iRM3AJQBPNZrn6gmTohIEmYF4KCEkC5MnCCiU29z8yQQfexjZjaK8+fNJLVf+9pYyRPDJk6cG7lkwwun+CKAgHWPp0dE0gCaQcvCe8vHr8MsL+8CqE+yfLy3bP0RgOa4ZRARnTr5PHD/vmk1/c7vAP/4H5slPs6NG0KGN+4R4gCyqvrlaVamk9eq2QaQCXjOBVBU1ZRvW1lE6qpaH/OQSVWtiAzViCIiOjv29kzX33/9r+bnmzeBr38d+OpXZ7bmFDD+mNQBTItjJkTE9caF3D7HyQHoHDsq4aSVN+oxkzAJIURE1Mlm/H3qU+brL/8y8G/+DfD7v29aWTMybkvqNoAdEWkAqKnqm1OsE7yWUA4ARKTQY7c0uoPUPoBd+4OIZAGs9TnUrqpWbTefqrKbj4iol7094ItfBG7fNtMmXb9uWld7ezNrTY0bpO4BcOC1PAK6x1RVZ9ZZ6QWVrlaWqjZFBCLiqmpdVbeGLDIJICYidhDP8QJcdYKuQyKi0yWfN7OkAyZIASY4zbC7b9xAsg1gkevOxwATlHo878IkUwxFVSv+n0WkNEKAIyI6O2yjxAapGRs3SH1+2l18I3JmUajXQst63+cBVIJaUl4rKwsAq6urs6gKEVE4LXmpDHMKUuMmTtRE5NenWpMQUNWmqm6qqnhfA1tjqrqlquuquv7EE0/Mu5pERIsTkSD1DQC/NM2KjMNr+RAR0bzYIDXBMk8jHW6cF6lqDoCIyH8QkfdPuU7DsC2ctvkDfUGLyQ5ERLMw55bUWGNSImKz6uIA6vPO7vOy+OroHpuKwcxOwSBFRDQLEUmcWHR2H2DS39dhbiy24uANuUREsxOFlpSqFgBARP4ezJRFNuW7PINl42MIzuYrwExy608Vz3kPIiKahSgEKQAQkd+FCQhNmAB1DUBWRMqqOlFShTe2tAET/BwARRFJwcwQUQFaXX4FESkC2PP2LbKrj4hohuacODHumNRtmACVUtV7vu1pmOmS9lT1/xq3Ut5Nur2mQ/Lvd4D27j4iIpqliKSg52CWj7/n3+i1cj4F4FcnrRgREYVQRIJUHL1bMIcwXW9ERHTazDm7b9wgVYWZhTzILbALjojodIpI4sSnAOx790dt4WRl3ByAXwSQ6v1SIiKKrCgkTqjqgYjchAlQNuXb3tGbm0EaOhERhUFEWlI2SaLirWj7FExral9V35hW5YiIKGSiEqQsVeUMD0REZ0VEpkUCAIjI5V7PLXi9KSIimoUotKRE5CMAdtA/1Xx5rBoREVF4RSFxAmaC2ccBfBHmvigiIjoLotCSgrmZN6uqX55mZYiIKOQiMuNEFcDRwL2IiOh0iUjiRAFmIlkFcI9JEkREZ0REuvvuwSyhUQGAea/MS0RECxKhxIlFr8xLRETzFoXuPrsy72njrYfVYhdYJCIin6Wl0CdOnDoikgdagakKszIwERF1Egl3S2qevNZNM2j6JRGJA1jHySzs9XGmabLL1avqCtBaGTgxUcWJiE6rObakQh2kvMlrtwFkAp5zARRVNeXbVhaRuqrWRzzUOoC6DYgw94FVxiiHiOj0O+vdfSLiikgJpnXU636sHIBSx7YSgOIYh3RhAlPVa4ltAdgdoxwiotNvaSn02X0z5bVgcgAgIr2SNNLoDlL78AUXEckCWOtzqF0vKNVhugqb3vGbXqB02ZoiIurA7r7+vDGkrlaWF1xgg4uqbg1ZZFAgak5aTyKiU+msd/cNIQa0EhyC9JudvYvXWmp6wc8GwXHGtoiITj9m9w3kzKDMDIANETmE6SLsStawvG7ELACsrq7OoCpERCEWpu4+ETnG6LNLRG5aJK/VNNRNyl434hYArK+vc+YNIjpbQpY48UWEdAokEXH6dPkREdEshKklFdIpkOxYUQy+BAc7poTgRAgiIpoGJk7057We6ugem4rBzE7BIEVENCtRSZwQkV9Aj0w6Vf3tScoeQhVmpogD37a4t52IiGYlTN19QUTkKZgbZ1e8TQrAv6iUAphWkIohOJuvAKAML4HBk/MeREQ0KyFLnAhSAvBdAElVfVFEjgD8PQCvANjxHmOzE77CtNIcAEURScHMEFEBWjfuFkSkCGDP27fIrj4iohkLe0sKQBJegPJ+rgOtwFEE8AUAXx63Ut6Y08CEDVU9QHt3HxERzVoEEieaAB7z/VyHGR8CTMsnPkmliIgoxOaYODFukKoC+CXfzzsACl4iRRGc946I6PSKQEuqADMmBaC1mu0rACowiQ63J64ZERGFU9gTJ1T1uwA+1bEtJSJPec8REdFpFYGWVCAGKCKiMyDs2X0icnnQPqr65jhlExFRyIU9SMEkRgzqkFwes2wiIgqzCEyL9KsIDlLrMEkT2bFrRERE4Rb2llSfZdm3ReTPYRYNJCKi02iO2X2zmAW9CiA9g3KJiCgMoprd50mix8zoRER0CoS9u2/AkvICoDZ2jYiIKNwikDjRb0n5Q0w4CzoREYVY2FtSIV1SnoiI5iHiiRNERHSahaklNczsEkGiNuOEiMTRvtxIhQsoEhEFCFOQwnCzSwSJ2owTSVXdtD+ISAlcip6IqFvIEieCZpdwYNaNKqB97agPAPgnAPJTqR0AEUkDaKpqNeA52/qpw6S914P2G1JORLa8VYGJiKiXMLWkgmaXEJHfBZBV1a4l4kXkNQCPT6NyIpIEsA0gE/CcC6CoqinftrKI1MfspisC+K6I2KQQJocQEQVZWgIePpzPocZ83S0ARz2eO8CEM06IiOt1t7l9jpMDUOrYVoIJNiPzgvEdr9wczOKNRETUKUwtqT6yAH4/YHsGE77Bey2hHAD4Wjad0ugOUvsAdu0PIpJF/3kEd233oIjkvTGpTe91uwNeS0R0NkUgSH0BwBdE5D/ABIrvwgSmHEyQGqs1MywRcRDQylLVpohARFxVrfeZCLezvCRMC9CWsyUiayISV9WDPi8lIjp7QpY40UVVN0UEMEvIfxQmsUJgkigKqvrFqdUwWMyrR68kBxcmmWJYRzBzDrYlXTBAEREFiEBLCr6usTiAp2Ay616cWs36c6ZZmKoeeONgdh0sB8DdXvt7+2UBYHV1dZpVISIKvznOODHJmBSAVmsj8i0OVa2MsO8WgC0AWF9fn8+VIiIKizAt1SEifygiDzu2HYvIoz6PueQmemNTREQ0TyHr7qui/YZdoP8s6PNgx5ti8NXNF7Q4nRER0ayEKUj5pwrybVvoja5eFl8d3WNTMZjZKRikiIhmZY7ZfWPdzCsiV0XkmY5t/6eI3BWRfzSdqg1UxcmEsFYcHRl6REQ0ZWEak+qhbfJVEflDmHuj1gBsichvTaFuVgzB2XwFdE+XlAOnMyIimq0IrCeVgjezg4h8BOYeo4KqrgPYQMBce6MQEUdEiiJShjeZrYiUvMlmAbTukSp4+6VFJA8zlx+7+oiIZilMY1J92GCQhEmisLM7HMLcTDs2G4CG2O9UpL8TEUVKBLr76gBuegsi3gJw4FvkcNTZHoiIKErCPi0SzHRIOzCtHYFpTVk5AOUJ63Uq/eQnP8HR0RHeeustPHr0aNHVITqTlpeXcenSJcRiMVy8eHHR1YmmsHf3qWpFRNbgZdOp6hu+p4swAYx8fvKTn+DVV1/FysoKrl69ivPnz8Ob/5CI5kRV8c477+DNN9/Eq6++itXVVQaqcUQgcQIw41DrAO563X4QkcdgAhenCupwdHSElZUVXLlyBRcuXGCAIloAEcGFCxdw5coVrKys4Oio13J11FfYx6RE5AZMgkTKe9hZyd+ASUP/wrQqeFq89dZbuHz58qKrQUSey5cv46233lp0NaIp7EEKpkuv4qWcdzYJtmCSKcjn0aNHOH/+/KKrQUSe8+fPc2x4XGGfcQKmS6/XUhYNTHkpjdOCXXxE4cH/xwlEoCV1gO4piawcODUREdHpFYH1pL4AYEfMRxEF8Jg3l18RJh09MaX6ERFR2EQkBf0mzPiTwLSs7PLxN1X1pelVkYiIQiUC3X1Q1YqqxmBaTTcBJFQ1pqr/emq1o1Pp4OAAmUwGKysrEBGsra0hk8ng4OCgbR8RgYigWu3fe2z329w8WVXG/3r/I5FItO1nVavVwP39j2bzZFm1ZrOJTCaDtbW11jnkcjnU65xsZRq2trZav9tEIoFKpXvh7Hq9jlQq1XcfmpGwt6T8VPVFAC/6t4nIM2xNUZDNzU0UCgW4rouNjQ24rou9vT1UKhXcuHEDjUaj6zXlchnJZDKgNAx8YyoWi0inzbzEzWYT+/v7KBaLKJVKqNVqcJz2HJ9sNotUKhVYlt23Xq9jbW0Nrusil8vBcRwcHh62gl+pVOr/S6C+tra2UCgUsLGxAcdxUCqVkMlkUCqVkM1mAZhrmUgk4LouyuUy9vb2kMlkUC6XW9ebZmiO2X1Q1YEPAJdhsvleB/BnAP5Rn/12ADwaptzT8EgkEjqMP/mTPxlqv9Nsd3dXAWgymez5vFWr1Vr7Oo7Ts8xkMqnJZFIBaLFY7Hp9qVTqek2j0VDHcTSdTnfVLWj/oGO6rhtYbq1WG/j6Rcjn83p4eBiJch3HaftbUFWNx+Ntfwf5fF4dx9FGo9Halk6n+/6tBOH/5Zg++UnVy5cnKgLAvg7xHjtsd98BzE2792DGnrZE5Nf9O4jIHZj08yTM3H5EbWyrY3d3N/D5oNZSJpNBs9kM7PKz2zOZ0VaGcRwHyWRyYDdiL/v7+4F1dRwH8Xh8rDJnqdlsYnNzs607NazlNptNuK7b9fvd2NhAs9lsHatSqSCbzba1hDv3oRkK05iUiPwT79unVPWmqn4AwEsANkXkkoh8QkReh0k9/5SacakvzrDOFEGVSgX1eh3FYnGk18ViMcTj8cAutJ0dM0Vkr67Afg4ODhCLxUZ+na3T/v7+WK+l/hzHwb179wK3A6artdlstrpc/VzXrBDEazMHYQpSMLNHfEHbJ5G1s583YTL8tmGCWCSCk4h0fZQXEVdE8iKS9L5G9obkF154AXfu3MELL7yw6Kq07O3tAQBu3rw58mtv3boVOPZkx6qGDTa25ZVIJFCv11EojLeIc7FYxMHBAVKp1NitsV4KhQJWVlYAmJbn2tpaK2GjUqkgkUi0EgX8x65UKq1EAxFBKpVCvV5vJagAplUqIq2fAbTt43/doDp9/OMf71vuJDrHCW09AfOBxM63Z4NS5+v8CS40IyFLnIgD6Pxosud93YJpPb2BCBCRJMx6V0EfvcuqmvD224cJvBOtMDysT37yk3jppenkmbzxxhv4oz/6IxwfH2NpaQkf/vCH8dhjj01c7jPPPIMvfelLY7/evskEvQENkk6nUSgUUK1WW60mG3AGJSnkcjnkcrmu7cVisTUIP8z+2Wy2dax0Oo1SqYRcLtcKFOl0GsViseuNc1yJRKJVH8dxUKlUkMlkkM/nsbGxgbt37yKVSuHw8BAAWs/dunUL9Xodu7u7aDab2N7eRi6XQyqVQqlU6mp12nPa3t6G4zgoFApIJBKBCSz+OmWzWezv7/csd9RMx2Kx2LertFgsIplMwnGcgeW+/vrrQx+XxhTC9aTaPpqo6hvelCKlaQQob1n4pqp2fSwVkTjM7BZ1eAsqBu03DPs6EWl7Z/OOceTbr+kFtMh54403cOz98RwfH+ONN96YSpCa1CSfbl3Xheu6bW+GtqtvUMvMvnED5g15a2sLjUajZ7D0799ZB79sNotsNotqtYpyuYydnR2sra1hd3d3rO5Hv2aziVgs1jZ2d/v2beTz+VZ3aTqdRiKRQKlUwrVr1wCglQ0Xj8fbMtxs3WOxWNf8V6qbAAAYVUlEQVR5dHa/bm9vt9K5/WUE1alfudPMcNzc3MTR0VHXWCZbTAsUgRknpsYLBoGtFhFxARRVNeXbVhaRuqpO84YUFx2BGMCRiMTVLFE/U5O0UDq98MILuHHjBh48eIALFy7gK1/5Cp599tmplT+ueDyOarWKZrPZFiASiUTbQHetVgv8RJ1Op9vub7JdfY7j9H2zWltba5VXLBaxtbWFnZ2dwFZU5/7DSCaTSCaTKBaLSCQSyGQyga2QUfnf5A8ODlpJCp33eMViMWxsbAAAbty4gVu3biGZTI6dwGFfF9RaWURqfaVSQaFQQK1WawVC+/fTa5mNxx9/fG71O7NCNiYFAL/oJUi0HjDTId3s3O57vi9vDKgEEyB6LeqSA9D5n1GCmX5pmsYbQQ+hZ599Fvfu3cNv/uZv4t69e6EIUABan/ZtC8gql8uo1WoDEypsF5wNdNVqNbBbrh/HcZDP58ceixpUdqFQmFp2mb9lYgNGrVbD4eFh26NcLrdlTNruul73egWxXYmJRKLvuNK0ujKHdXBwgNu3b3d9cLFjkJ0fTibpUqYRhWxMCgB6JUT0SjVXAF/uV6DXEsoBgIj0etdIoztI7QNotftFJAuzhlUvu0N0Dx6he+b2yAauZ599NjTByUqn03Bdt2ssyL7xDRpnsF1+5XK59eY0zk2bGxsbrZtFR800HGRWb+L+cnsdI5lMolarodlsts5vc3MT+Xy+b9k2UaJQKGB93cwZbceeJjHpmFS9XseNGzdQLpe7WoWO48B13dZ4nGU/HEza3UpDsN19qmZ8aoaGCVILmSzWy67ramV540UQEVdV66q6NYXD1REQlObR1XeWlEolpFIppFKpnvdK9ZPL5dqSF8bhOA42NjbaZjQYpx7FYrHrtXbbtO+VisfjrZZa0LiMvx62tXj37t1WRmWv7Md6vd4aU/PPyjGsflmVk3QNNptNpFIpbG9v9ww46XQaW1tbbccplUqtDzM0YzYwhSFIqZn2aBHsar+9/mtcmOAyMVU98K8t442FcbmRKbNjN4VCAWtraygUCojFYtjb2xsqldtm+W1tbaFcLo9dj3w+jzt37uD27dtd5fQLnjbdfWdnB1tbW0in07h27Rpef/311n1g/vJsC2VjY2PiwLW9vY1MJoNMJoNcLodms9l6U04kEigWi8jlcnBdF/V6HQcHB63uUBvE7Bt6qVTC7u4uXNdtBT+73yity17lTspOd3R0dIStrfbPoMlksjWl1ubmJlKpFHK5XOtvaJK/CxrBkjdSdHx88v2sDDMtxawfMEvRJzu2xU31AvfXzv2HPE4cQN57fdFfhu+5pPfV6VNOFqbbcX91dXXg9B+qnH7Fr1artaawAaCu62o6nQ6cFqlcLre91nVd9f4uWhqNxkjTIqmqFotFBdCaxshOi9TvYevSaDS0WCy2pupxHEeTyWTXlEiD6hAkn893nZ+/PDsFlOM4ms1mtdFoaKPR0Hw+3/rduK6r+Xw+sNzO58rlcus6xONxLZVKXb/LfnXqVe647LXs9fD/Lg8PD1u/D9d1u/5WhsH/yzH9s39mOvsePBi7CAw5LZLonNII+xGRQwA59Y0deWnhNVXtakuKiAJI6Zip6NO0vr6uw9zh/vLLL+Ppp5+eQ40oTNbW1gLHVSgc+H85pjt3gH/6T4G33wYuXhyrCBGpqWqvxXNbZtxOm1yUZ34gAsAARaePv7tv1oea+RHGZ8eb2kZnfUGLC/dQqDWbzalnEBKFgh3DP8tBSk3CRB3BqeFNne7NvERT5zgO1zai04ktqZYqzJRIfnEw846IaHFskJpDTkNYglQM3S0mwMy23jldUs7bTkREizDHltTC5u7zxpY2YO53cgAURSQFM0NEBWjduFsQkSLMzOt2Lj929RERLcpZCFLemNPAFpGaWR848wMRUVhwTIqIiEKL2X1ERBRaZzBxgoiIooLdfUREFFoMUkREFFoMUkREFFpMnCAiotBi4gSdZgcHB8hkMlhZWYGIYG1tDZlMprX8t93HW4F54IKIdr/Nzc3A1/sfiUSibT+rWq0G7u9/+FetbTabyGQyWFtba53DqEumE0XWWbiZl86mzc1NFAqF1uqqrutib28PlUoFN27cQKPR6HpNuVzuuYx4pVLpe7xisdi2NPr+/j6KxSJKpRJqtVrXEvDZbBapVCqwLLtvvV7H2toaXNdFLpeD4zg4PDxsBb9Jlk4nioQ5BqmFr8ob9Ucikei7+qTFFUBPVr9NJpM9n7fsqrbJZFIdx+lZZjKZbK3OOuzKvI1GQx3H0XQ63VW3YVbRTSaT6rpuYLmdq/OeJfl8Xg8PDyNTrir/L8f2la+YlXm/852xi8CQK/Oyu4/mxrY6dnd3A58Pai1lMhk0m83ALj+7PZPpnIO4P8dxkEwmB3Yj9rK/vx9YV8dxzuwCh81mE5ubm21dtmEulybExAkay+YmcP9++7b79832BatUKqjX6yMvAhiLxRCPxwO70HZ2dgAEB7dBDg4OEIvFBu/Yo077+/tjvZboVGAKOo3l2jXg5s2TQHX/vvn52rXF1gvA3t4eAODmzZsjv/bWrVuBY092rGrYYGNbXolEAvV6HYXCeCu+FItFHBwcIJVKjd0a66VQKGBlZaVVvk32CGpJdO6ztbU10jEA07pdW1trJYVUKhUkEolWmf7zq1QqrUQREUEqlUK9Xm8lwQCm5SsirZ8BtO3jf92gOn384x/vWy4t0Byz+5g4EQaf/CTw0kvTKeunfxr4+Z8HnnwS+Mu/BJ5+Gnj+efOYxDPPAF/60tgvt29KnYkKw0in0ygUCqhWq61Wkw04g5IUcrkccrlc1/ZisYhsNjv0/tlstnWsdDqNUqmEXC7XehNPp9MoFotwXXfk8+tkMweLxSIymQwKhQIymQwODw9b+1QqFWQyGWSzWezu7uLg4AC5XA61Wm3oxI1EIgHgpBvWlpnP57GxsYG7d+8ilUq1jmufu3XrFur1OnZ3d9FsNrG9vY1cLodUKoVSqdTVsrX12d7ehuM4KBQKSCQSgUky/jpls1ns7+/3LJcWiNl9NLaVFROgXn0VWF01P4eAP317VK7rwnXdtjcq29U3qGVm31QB82a5tbWFRqPRM1j69++sg182m0U2m0W1WkW5XMbOzg7W1tawu7s7lTdTf1aizRxsNputet++fbstcCaTScTjcaRSKWQymYF1aDabiMVibeODt2/fRj6fb3XJptNpJBIJlEolXPNa4xsbG62xN1s/4OT3E4vFun5XnV2829vbSCQSqFQqbWUE1alfubRAzO6b7QNmYcXObXEAWQB5AGUA7jBlhS677xvfUL1yRfXTnzZfv/GN+Rx3gHw+rwC00Wi0bY/H4wqg9bDZcTY7r1wut73esll9qiarDkNk99n9gjL4RsnuC9JoNNR13b6ZiMMI+j2VSqW2bbauQZmEjuNoNpsd6hj+jDn7+wp6JJPJ1u8uHo9rsVjsOvbh4WHb9Rqk83oF1WmcckfF7L4x/cEfmOy+g4OxiwCz+7qJSFJEsgCSHdsdAOuquqWqmwBKAIJT0MLMjkHt7ACf+5z56h+jWiD7Sdy2gKxyuYxarTYwocJ2wVWr1VZXX1C3XD+O4yCfz489FjWo7EKhgGazOZVMtH7dorbrNKhl4bru0DcU+19vX1Or1XB4eNj2KJfLbVmZtruu1/1kQWxXYiKR6DuuxNZSRJy17D4RSYtIYP+EiMRFJGsDTK/9hqGqVVUNGl120b5K8D4A1wte0bG3ZwLT9evm5+vXzc9e0sIipdNpuK7bFYxc10U8Hh/45mS7/Mrlcts40Kg2NjYAYCaBal5vsPY4QcGoXq+PVQ//a+zv2j5swEwmk6jVamg0GigWi6hWq4Gzd3RKpVIoFApIpVLY3t7GvXv3Rq4fhcxZmhbJCzrbPZ5zARS9Fo4NMDlv+9SoWaLe/7FwHUBTzRL30ZHPnwQo6/p1sz0ESqUS6vX6SJ/A/Wyiwu7u7lgBCjAtlI2NjdYYz7j1CHptsVicy71S6+vrALpntqhUKq2ki1HF4/FWa7BT57naFmk8Hm9lbfbKsKzX66hWq61ElWE+kPiNe5sAzdhZSJzwAk0BQA3AUY/dcjBdb34lAEUAo/8n9qGq/o+lOQC3p1k+mU/ixWIRhUIBa2trKBQKiMVi2NvbGyqV22b5bW1toVwuj12PfD6PO3fu4Pbt213l9LrR2NY/FothZ2cHW1tbSKfTuHbtGl5//fXWfWD+8mya+8bGxlQDl+M4rezCo6OjVrZdoVBAOp0eO3Fje3sbmUwGmUymFYhLpRJc10UikUCxWEQul2t1KdqMQlsn4CRwlkol7O7utlpiNvg5jjPSvXK9yqUFO2uJEwAOASR7bHc7tjmm2q2fszBBq9cjqFztU5csgPSwdQ9d4kQE1Go1TafT6jiOAlDXdTWdTgdOi9Q5YO66blsCherwiRN+xWKxLfnAJiP0e9i6NBoNLRaLGo/H1XEcdRxHk8lkVzLBoDoE6UwQUe1OnLB2d3dbiSeu67ad/6jH8NfZTjNlkzAajYY2Gg3N5/Ot37/ruprP5wPL7XyuXC63rnU8Hm+dT1DiRL/6Bh1zUvy/HNPXv24SJ771rbGLwJCJE6Jz6FMcREQOAeRUterb5gBoAFjRjm43EVEAa9re+hnleKqqErA9CZixq2HLWl9f12FmH3j55Zfx9NNPj1RPir61tTWUy+UzO11S2PH/cky7u8BHPwr8p/8E/J2/M1YRIlJT1fVB+y18TKqPGAB0BiifqY5LiUgcwJENUCIy3qAHUQcGKDp1OOMEANOtN1VeIEp63xdh7peqeuNj9wA4Iq0GVh1A4DoQXhp7FgBWV1enXU06JZrN5shzFRJFwllInFgENVl8BwA2O7bXAQw9NYOaLMMtwHT3TbOOdHo4jjN2FiJRqHGC2RORu1eJiOi0Y5ACYLrbAG9syvIFLa7TTUS0CGdtxokgXsJEHd1jUzGYG20ZpIiIFuEszTgxQBVm9ge/uLediIgW4Qx298UQnM1XQPfMEjm0z7MXGWG4J42IDP4/TuAsZPd5Y0sbMPc7OQCKIpKCSQuvAKbLT0QKXrr4nrdvMYpdfcvLy3jnnXdw4cKFRVeFiAC88847WF5eXnQ1omdzE7h82Xxvg9T9+2Yi6xnME7qwIOWNOQ1sEfnSxiPt0qVLePPNN3HlypVFV4WIALz55pu4dOnSoqsRPdeuAb/wC+b74+P2JYJmICzdfadeLBZDo9HAa6+9hgcPHrCrgWgBVBUPHjzAa6+9hkajwVnWx3H9OmBvUv/KV04CVOcKDFNypm7mXaSLFy9idXUVR0dHeOWVV/Do0aNFV4noTFpeXsalS5ewurqKixcvLro60fTRjwJPP22C06c/PbMABTBIzdXFixfx5JNP4sknn1x0VYiIxvfd7wI/+IEJUL/7uyZIzShQsbuPiIiG5x+D+tznzNebN832GWCQIiKi4e3ttY9BXb9ufvZWaZ62UKwnFWXDridFREQnTsN6UkREdMYxSBERUWgxSBERUWgxSBERUWgxSBERUWgxu29CIvIDAN8b8+VXALw2xeosEs8lnHgu4cRzAd6vqk8M2olBaoFEZH+YFMwo4LmEE88lnHguw2N3HxERhRaDFBERhRaD1GJtLboCU8RzCSeeSzjxXIbEMSkiIgottqSIiCi0GKSIiCi0uOjhAohIHMA6gDoAF0BdVauLrVV/IuIAyAKoADgCEAOQA7Drr3tYz01E0gCaQXUZps5hOq9e5xK1a+SdhwtgzftaUtVKxz6RuDaDziVK10ZEkgBSAF6HOZ+aqm517DO/66KqfMzx4V2s3Y5tZQDuous2RL3V92gASEfh3AAkvfomx7keYTqvIc4lEtcIQBpA3PezA+AQQDZq12aEcwn9tfH+vpId22oA8ou6Luzum78cgFLHthKA4gLqMqoUgBUAa6q6oh2fehGycxMRV0RKMP8wRz12G6bOCz+vIc8FiM41clX1wP6gqk2vDv66ReLaYLhzAaJxbXIB26od2+d7XeYVoflofZo4RMenCZhPXrroug2ot4uAT+9ROTevbkGtj4F1Dtt59TmXSFwj73g1AE5A/dXWLwrXZoRzicq1KQModmwrwnT5DV3PaZ4LW1Jz5PVLd30SVvPJCyLiLqJe0xDFcxumzlE8r17Cci7e8VzvESgq12aYcxlGGM7FO15GVQsdm9MA7g5bz2mfCxMn5isGnFysAC7MIGNYud4AMWDO5UhPuiyieG7D1BlD7BOm84rENVLVlYDNSZiEkLp9I4vCtRl0Lv76ROHa+IlIFsCBqm56m+b+P8MgNV/OoiswgSMA0PaMpbKI2G1RPLdh6hyl84r6NcoBuON9H/Vr4z8XIGLXxgumKcC0rnxPzf26sLuPhqKqTe1IQ0V0Ej7OhChfI+8T+5HvE3tkBZ1L1K6NqlZUNQegICI1L518IRikFsDrsz0N6jBdGK3zieK5DVPnKJ6XJ/TXyOvay6lqKuC5SF2bfucSIPTXxuuyKwG4598+z+vCIDVfth825t/ou5hhGttoIyL5gM12YNTfxxylcxumzpE5rwhfoyKAGx3bonptgs4lytcGMCnojneT79yvC4PUHHmfSuro7rONoXuQNTS8T4fFgKwc+0dYj+K5DVPnqJxXVK+Rd+9XoXOQPYrXpte5ROXaeJl5jT5de84irguD1PxVYaYK8Yt720PJ+6PKBfxxJWEyf+w/ZeTODcPVOfTnFcVr5I3dFP11FpGk7808Mtem37lE6No4aG8JWfZ62BuW53tdJr35i4+Rb5Zz0D1dyC7CPy1S2l9HnNzE2DkdTCjPDQHT0Axb57CdV59zicw18uqa9d647CMJM+ddpK7NkOcSiWsDII/uG5N34bvBd97XhetJLYDXnL4FYA/mU8qBhmAS1kF8k2g+DvNH2PbJ0dsnNOfm9YFvePVIw3xCrML88/hTgQfWedHnNcK5hP4aeefS6PF0XVXXfPuG+tqMeC6hvzZeHbIwE8sOmmB2LteFQYqIiEKLY1JERBRaDFJERBRaDFJERBRaDFJERBRaDFJERBRaDFJERBRaDFJERBRaDFJERBRaDFJ0Knnzpql393xoiUheRHYXXY9FEpGiiHBWAQrEIEW0WNdg5nkjogAMUnSmeJ/aO5dMWNhxVTWjqjLv+hBFBYMUnRneZKB5mFmqT/1xiU4DBikiIgotBik6E0SkjJMlFcpeUkWjY5+0iNS852rectn2uaLdX0RKInIoIo63mmnZW9FURWTX36036LhBSQMiEvfKsfXIdjxftCuoduw3dEtt2DJ61C/r7e/0KK/td+j9juwxDr0lK4LqZOvS8PZLdjw/8vUZ9vdB4cUgRWfFbQAp7/sczDo5Cfuk98ZZhlmjKQOzXlNbwPH2q8GsOFpSs6JqrqP8GMxidkMdt5NXj5p3/BSAuwBK3tLkfo5X35JXruv9PIpplBFU3p2O8uwj4+3X6xjbXl0K3s+t3/8E14eibl4rPvLBxzwfMBlzCiDr2+Z624JWtW3At/qot61mtwEoeq/dHXDceOcxBhy3aP4N2+pR6nEuyY66pDvLQceqqn3qOVQZnfXztmV77Wfr2LEt79uW9rbFBxzD/s5K07w+fETvwZYUnXleF5cDIO91JanXxWWXAvfLdRXgo6oH3rcjZxB63VcOTGvCX2YVQBMnLRHLv8rp4ajHm2IZfvsB5fmPYVeijfUrRM2KtVUA69O8PhQ95xZdAaIQsAElARMM/I78P2jHct9AqyvqllfOJOnt9rVdx/C2tZWtU+jOmkYZQ5Q37jHqMF13E10fijYGKSJfUBj1Tc6bLcKF6W6yrYha71cMVQ8XwEHHcy7aWylngQvzOxn7+lD0sbuPzpKjoI1eF10TJtC06Zch5g3aJwEUVHXLKyfoTTTwuAFsEGrrsvJaajYpIQx6Jn5Mi/d7TwKojnt96HRgS4rODFVtiggA5OxXVbWZd7dhUsRttpsDEyzq6DHOoap1EWkCKHrlNXGSmTbscTv3y8Fk88VgMvtsK63ijU3N2yFg5hiEGSNKwiROTJ33uy/g5JybMJmCwBjXh04HtqTorNmEeaMtwtelpqoVmBaCA2AXJh26joCg0+E2TBJA2SuzV2sn8LidVHULJvXcpm/nYFpqnUkT87IDU197bo9jdkHhLszvfhem9ZmwY1wTXB+KOFHl5MNERBRObEkREVFoMUgREVFoMUgREVFoMUgREVFoMUgREVFoMUgREVFoMUgREVFoMUgREVFoMUgREVFoMUgREVFo/f+/q9kMB+mYyQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import scipy.sparse.linalg\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rc(\"text\", usetex=True)\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "\n",
    "n = 300\n",
    "ex = np.ones(n);\n",
    "A = -sp.sparse.spdiags(np.vstack((ex,  -(2 + 1./n)*ex, (1 + 1./n) * ex)), [-1, 0, 1], n, n, 'csr'); \n",
    "rhs = np.random.randn(n)\n",
    "\n",
    "res_gmres_rst = []\n",
    "res_gmres = []\n",
    "def gmres_rst_cl(r):\n",
    "    res_gmres_rst.append(np.linalg.norm(r))\n",
    "    \n",
    "def gmres_rst(r):\n",
    "    res_gmres.append(np.linalg.norm(r))\n",
    "\n",
    "sol = scipy.sparse.linalg.gmres(A, rhs, restart=20, callback=gmres_rst_cl)\n",
    "sol = scipy.sparse.linalg.gmres(A, rhs, restart=n, callback=gmres_rst)\n",
    "\n",
    "lim = 300\n",
    "plt.semilogy(res_gmres_rst[:lim], marker='.',color='k', label='GMRES, restart=20')\n",
    "plt.semilogy(res_gmres[:lim], marker='x',color='r', label='GMRES, no restart')\n",
    "plt.xlabel('Iteration number', fontsize=20)\n",
    "plt.ylabel('Residual norm', fontsize=20)\n",
    "plt.xticks(fontsize=20)\n",
    "plt.yticks(fontsize=20)\n",
    "plt.legend(fontsize=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### How to avoid such spoiling of convergence?\n",
    "\n",
    "- BiConjugate Gradient method (named **BiCG**, proposed by Fletcher, [original paper](https://link.springer.com/chapter/10.1007%2FBFb0080116)) avoids that using \"short recurrences\" like in the CG method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Idea of biconjugate gradient\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. improves the stability (later in the lecture)\n",
    "\n",
    "Let us do some demo for a simple non-symmetric matrix to demonstrate instability of BiCG method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0xa14f7b2e8>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAakAAAEcCAYAAAB9K0udAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsvXl8Y1d5//8+2r3L8oxntWdGnjV7bE8mgRBCRqYU8m1Y7IRCS79tiU37LS1QsJkC/cK39EU0KV3oAva0/EpZ2oxNKFsJ2CEEyDq2s88kmVizb56xLe+LlvP74y6WLHmTbVn2nPfrdV+S7j333CNZ1uc+z3nO8wgpJQqFQqFQZCKW5R6AQqFQKBTToURKoVAoFBmLEimFQqFQZCxKpBQKhUKRsSiRUigUCkXGokRKoVAoFBmLEimFQqFQZCxKpBQKhUKRsSiRUigUCkXGYlvuAax01qxZI7du3brcw1AoFIoVRUdHxxUp5drZ2imRWiBbt26lvb19uYehUCgUKwohxKm5tFPuPoVCoVBkLEqkFAqFQpGxKJFSKBQKRcaiREqhUCgUGYsSKYVCoVBkLEqkFAqFQpGxqBB0HSFEOVAJuIG9QIOUMrC8o1IoFIqrGyVSgBDCDVRKKZv01z6gFShb1oEpFApFCnRdHuJCcIzbd6xZ7qEsmFUlUkKIaiAopWxLcsywlAKAFwjEtPMCDUCT/rod8Aoh3FLK4NKPXKFQKBaP/V9+HICTD7xrmUeycFaNSOnWzyGgJskxL+CXUlbF7GsWQgSklAEpZacQoirmlEo0sVMCpVAoFMvIig+cEEJ4hRCNaNZQ7zTN6oDGKfsaAb/xYsr8Ux1w/2KOU6FQKNLBpYGx5R7CorLiRUq3hOqM+aRpqAY6p+xr1/fHIYSoBR6SUrYs4jAVCsVVyOhEhA99/Vne6B5M2zVfODPpAApHomm77lKx4kVqNvSgiAQry3Dl6a5Ao60Pba5KCZRCscyEI1E+8s0OfnX88nIPJWXO9I3wy9cv03GqL23XfPFsv/l8LKxEaiXggUlRSoIXzMCKXiOYQg/CUCgUy8T/vHyRR165yD/+/I3lHkrKjExEAM2iShcne4Zjrh9O23WXilUTODED7tka6NbUo4BbCGHsDgBJLSrdJVgLUFpaujijVCgUcRz6pTZNXLY2Z5lHkjqGSIyG0mfRDI5NCtPYhLKkVgX6vFahlFLEbNOukZJSNkkpK6WUlWvXzlqzS6FQzJPBsRAvndPcVkPj6bNCFpuRccOSSp9FMzQ+ea3R0Mr97AyuGpHS56YUCsUKoGdownw+OBZaxpEsjBFdJNIpFoNjIZw2S9qvu1RcDSJlhJZ7YnfGiJZKfaRQZBhXhsbN50NjK3deZdR096VTpMIU5zv16yuRynj0gIkAiXNTHrQFu0qkFIoM44puSW0uzIqbY/n2M6d4KSZ6LdOZDJxI39zQ0FiY4jyXdt3QyhV4g1UvUjptaFkkYinX9ysUigyjZ1izpLatyTHdfZcHx/nM917m9//9yHIObV4YIjWWJksqGpUMTYQpzjMsKRU4kWl4SB7N10BiuqQ6fb9CocgwrgxqltSWomwG9UCAn796CYA818oJSjbcbYa7743uIar+9nEef31p1n4NTYSRkkmRWgVzUivnrz0N+tzSAbT1Tm7Ar+fhazUW5Uopg0KIBiGEHziit/UrV59CkZn0DI9TkGXHk+1gaDxMNCppPdoNwIYC1zKPbu4MG3NSExEmwlE+9G/PcL5/jMde7eatOxc/MtiYv1urRCpz0OecZrWIpJSdJKZGUigUGUjP0ARrch3kumxIqQUDPPHGFQD6R1dOtJ9hSY2EIlweGud8v5ZXLzYwZC680T1EjtPKhoKshGPHLw0yHo5y3aYCc/7OnJNaBYt5V5u7T6FQrAIuD41TlOskz2UH4LkzfYyGIggBwZGVI1LmnNREhOGY9Utn+kbn1c+ffKeTv/rR0YT9fcMTVP3dL7n7H38NwNC49tmszVdzUgqFQrFk9AyNszbXac4/PRXoAeD27WtWlCU1EjMnZYjU+nwXZ3tH5tXPxYExApeHE/Y/+LPX4l4P6JZUQZYdh82yKtx9SqSWiVcvDvD1X59ASrncQ1EoMo6e4QmKch2mJfV0Vw9Zdis3lxYyNB4mtEKyexsh4JpIaYKxe0MePcMTcZaVQSgSTfhNiEQl/aMhTveOJByLTb47FoqY7r58l40suzVtUYVLyYqfk1qp/NezZ/j3J0/yZNcVrt1YwK71eQyNh3nlXD++a9axoziP4jwnFouYvTOFYhUxEY4SHAlRlOMk16n9RL1wtp+KLYUU5TgAGBgNUZTrnLGfsVCESFSS41y+n7lYd5+RrmjX+jx+8dplzgVH2bkuz2wbiUr2f/lx7ttbwv9523Zz/8BoCCm1vnqHJ8z33T8S4kzvKNuLc3mje4iL/WNm4ESu006W3boqFvMqkVom/vLua8jPsvOdZ07z81e7ieo3SDaL4BtPnQLAYbNQUpjFunwX3YPj7Nvm4U/u2k6+y76s/3jzQUrJ+f4xRifClK3NRQjB4FiIXx2/oi06zHcyMBbmfHCUGzYXsHNdHmtm+fFRrD6klHzjyZO8/dr15pzTtrU55MeEm1+/qQB3tmZZfeyh5ynIsvPJt+9i65ocsw8AIQR9wxO896tPIoAf/entZDuW5//FzN0XipjJZveszwfgTO+IKVKf+d5LeHIcnO4d4bnT8WU9+kYmU0T9x1OnePP2NdyyzcMr57VFzb4963ije4gL/WPmmrI8l41sh9VMy7SSWRm/dKsQi0XwiaqdfKJqJ2OhCMcvDTEWjnDNhnyOnOzlTN8oZ3tHCFwZ5tLAGCWFWbR0nOXbz5zGIuA9N2/m6UAPnhwH24tz+Y1r13HX7nX0jUzw1V90MR6O8tzpPk73jnDHjrVcuzGf23es4ebSwrS9x3PBUT793Rf51XEtKqvUk81NJW6eeOMKPcMTSc8RAt51/QYa3rGbPJcNd7YjbeNVLB+vnB/g8z88yhd/fIzP/9a1AFyzId+0pAD27ykmrN/NGd+p8XCUQx/S1un/zr89gyfHyT/+9s18quVFzvWNMhGJ8jc/fZ2//F/XzHj9y4Pj/Ol/PsfgeIgf/sntxFRDWBAjursvHJWm+O5arwnTGX1e6nxwlG8/c9o8Z+rcU19MoMg/PHqcR16+yE8/fgevnB8AwLenmK893sXFgVGGxsNYBGQ7rLh0S+rZE728fmmQ37l1y6K8p3SjRCoDcNmtXL+5wHx9567ipO1OXBnme51nefn8AN/tPEvFlkJynDZ+dfwy33vuHPkuG1GpuUtcdgvXby7gphI3//PSBR555SJfbn2deys382e+nQRHJjjVM8KZ3hEu9I9xzYZ8fuPa9ThsFpw2y4LdjM8Eevijb3cyHopQ/45dFGTZ+ekrl2g/2cvNpYV8+C3b2FDg4vLgOFkOK+vzXbx0rp+nA71848mT/OjFCwD81o0bufuGDbxlx1qyHFaklLx2aZBHj3Wz0e1i57o8Xjk3wEQkyu3b15h31Vcr4UiUH790gWyHDafNgs0qKC8txGW3JrQ9clKrA7p3qyfhWLp5qksLjAhHJf6fvIrLbmHbmpy4if/bvEVmZnQAi4AXz2pl4o5dGOCJN7Q+fHuKaTt2iY/7dvLqxQEeefnCrCJ16FcBMzhjeCISJ44LIdbdZoSdl3qysVkElwa111MLIp7uHSEUiWK3aiEDwZH4G7rXLg1yLjjKy+f72Vjg4pqNmmWmWVJhcp02hBBkOay0HbtE2zFtEfRv31KKdQVOHyiRWkFsW5PDJ96+i0hU8tK5fm7YVIDFIghHovzq+BV+dvQi0Sj84Vu2xfm6H3jfDQyOhfja413882NdHG4/G9evy25hLBSl/rsvArDJncUf3L6N3711Cw7b3GNrxkIRnjnRy2OvdvOtp09R6snm0Eduo2xtLgAf3Jd4J7elaFJU7txVzJ27ivngvlIe7jzHwFiIbz19ih+8cB6HzUK+y46UclorDOBNZUWUlxZyz00b2RHzGUxlaDyMw2ox31+sqyiW8XCEK0MTdJ7qwyIEd+0uJhyNkuey82TXFZ4/E+Qd167Hq7/H5WRwLMTvff1ZOk/H1/fcXJjFp39zN4XZDl48289L54KcC47xwpkgDquF+nfsYt+2orgbpXTzZNcVvGtzWJPj5NmTvdxU4sZqEeQ4NHHdWpSNzWqJs6xrKkp4qP0MlwbG+K9nT+OwWXDZLHzi8As4rBY+eGsph9vP8JOXL9I9OEae006Ww8q3nj7F469f5p8/UG7+/WOtl96hiUUTqZGJiOZ2m4hwZWjctHLW5jnpHtBEqvN0Hy67hbtv2MhoKMKPX7zA6d4R8/+mL0nI/c9f7eaZQC/lW9xkO2zku2zanNR42Aw2yZpyY3Khf5TNhdmL8r7SiRKpFYjVIripZDL7k81q4W27i3nb7uQWGECey86nfmM3vj3rePXiIO4sO6VF2ZR4sslz2ug41ceRk32EI1F+/cYV/upHR/n2M6f44zu389rFAV440487205RrgOnzcqWomyePdHLkZO92CwW1uU7ee3SIGOhKE6bhd+8fgNffPd1FGTZ5/3+SjzZ/JlvBwCf+o1ddJ7q4xevX2ZwLEwkGuWGzW5+49r19AyPc+zCANdsKCDLbuX7z5/j4efO8dXHu/ja412864YN7NtWxJaibEKRKN0D45zvH+XpQA9HTvaRbbeyZ0M++Vk2nj+j3aHfVFJAiSebiXCU45eGeOV8P8Mxd8N2q8AiBNdvKqBdvwP++7bjHPpQJbvX53EuqE1k57vm/75TpeNUH78+foVHXrnI65cG+dt7b6Qw24HVIhgNRXjgJ6/yJ995LubzzaKkMJs/vrOM1qOX+OKPj2G3Cpo+VMnbprHiU2E8HOHRY91UbCmkd3iC5vaznOkb4S/vvoYSTzZPdfXwme+9xDuv38CzJ3p5T/kmbtjs5tmTvaZ1IITgRx+9nS1F2o9r7Pfp3Tdv4qH2Mzz+2mVaOs7yrus38P69JTT9MsBNJW7W5DrZrbvWar72FGOhCN//P7fzTz9/g4sDYzz401f5zLs0C+tM7wgWAVEJvSMTlBYt/MdcSsloKMLmwixGeke5MjRBjkOzcorznHQPagt7O08HuXGzm7+puZHnTvfx4xcvELg8bIqUYUndV1nCljXZ/NezZ/jXXwW4ODBmZq3YUJDFhf4xzvSOUOKJX/C7Z0M+xy4McLp3RInUSkavzluNlpWiHGiaoeT8iuXm0sKk81KVWz1U6m6fP7lrO4+91s0Xf3yMTza/gNUiqCgt5LVLgwyfDjM8HmE0FMGdbadqzzrGw1H6RiZ4/95S3rprLbd5i5K6l1LBZbfypu1reNP2NQnH1uY52a1PQgN8dP8OPrp/B73DE3zl0eN877lzfP/58wnn7VqXR+0dXoIjE5y4MszZvlHKS93YbRa6uod44o0erBbBng15vPvmTVy/qYBSTzYneoZ57eIgQ+Nhjl8a4uO+nbz75o3c/x/t/N7XnzX7t1sFVdes4zeuXQ9od9OvXRyk83QfD7z3BkqLsvn5q90UZtu5ffsa03rrHZ6g8ZddfLfjHDdsLuAz79pj/lAZ9I+G+LvW1xkaDzMejnKub4TO00GEgB3FufzD+2/i7hs2xp2zf3cxz57sJRyRXL+pgMKcSWvkY76dnLgyzCcOP0/tf7TzmXfu4e4bN1KU4zDHFY5E6Z8mmu5C/yjN7Wc5en6AmsrN3Oot4uHOszx3OsjloXFz7gjAabMgJXxRHOWjd+3g977+LFkOK//02Btk2a28t3wzu9bl8fVfn2B/zA3XdZsmLbzYQArD2vrij48yEorwR3eWsXNdHvu8RWYb4/txqkeb//lUywvkumwwAP/57Bn+4p17ADjTN8KNJW6eOx2kd3h+2SCmYywURUrw5Dg50zvKlaFxsp3a/8XaPBdn+0YYC0V45Vw/99/hBcC7Rvt7n7gyBKwDtO+F1SJ44H3XI4R2k/TAT14F4A5DpNwuXjgTpHtwnIZ37AbghTPaz9eHbtvCgYdf0ubApi3lmrkokZqkWUpZASCEaAcOkZiU9qpACMFdu9dx585iXjrXjyfHQYln8g4sHIkyOBYmx2mblzswXXhyHHz+t67lL+++hgsDY5y6MozTbqE4z0VxvhOnbWYBjeiT81P998mEEuCh2tv4n5cvMDoRodSTzTMnevnec+f4n5cumm1sFkGuy8Z7v/oE2Q4bvbrL0renmMJsB693D/HqBW1u7W27innudB/v+ecnaPjN3dyxYy3/8OhxTvdogTTBkQnW5jlx2jT316d/czcf2Fc6rfVms1p4U1nysTtsFnatz+M7H76VP/p2B5//4VE+/8Oj5LtsvGXHWj6wr5RvPnWKnx69yB071vLbt5RQlOvkudN9CARf+skxolL7zB955SIOm4WJcJRcp42h8TB/un8HeU4b2U4r77xuA9959jQP/vQ1Ok4FcWfb+enH7qB7cJxNhVmmi+2Rj90x7d/GZp38vmU5rOzfXczPjl7iPTdvinNxG2wocJHnsjE4Fua6Tfk8e6IXixDYLIKh8TBXhiawCO1G4sbNmkjFFlycD1JKng70cqvXgxDCzNtnhM1fGRw3o3KL87XP8KVz/YSjknL9xrEg205xntMMigDN3VeYbTdvGj502xaafhlgba7TTJP0v27YyC9e09ZM3blLE66iXAeD42HuuWkjn/3vlzndO8LxS4P86vgV/uD2bZzpHeFnRy/xB2/eumiBIkvBqhIpIUQ1Wo2ohBIcQohytHIdAbQEswGjnX6s12irJ6T1pWfUmYvFIrixJDGpvM1qibsbz1QsFsEmdxab3In5zmZivpPLhTmOuPm2t1+7noZ37OaN7iGcdgvZDisFWXZ6h7XIy76RCX5n3xY6T/fx923HcWc72Lkulw/u28IH9pWwvTiPM70jfOyh5/nM914GtHnDm0rc3LKtkN9/87ZFD3YoyLbzrT/cx3NngrxwJshrFwf52dGL/PglLYDlndevp/NUkI98Kz79pW9PMZ+7+xo2urP47+fO8XSgl9++pYSbSws52zcSN+cI8Ie3b6N7YIwfvHCeg/feSGGOY97fpS/X3MiOdZrF0fShSsZCEbMS7VSEEOxZn8/x7kHev7eUz/639nneuWstv3jtMid7hs0ABcOF3jeSmkg99lo3f/Dv7Xzu7mv4rRs3mgtpTZEanmB3niakxXlOeoYnePaE9rNzc+nk/9mt3iKe7OpBSokQguDIRNxcXLbDxr///l4sMcLy7ps38dXHuxgaC5suzm99eB+XBsbIdtjY5M7idO8o33z6FP/x1CnuvmEDX/jhK7Qd6+b6TQVsKcrm4c5zfOStXlOwJsLRjLgJXTUipYtKUutHd+X5pZRVMfuahRABPRO6F5jq2usVQpTriWkVinnhsFnMeRWDbIeNv37P9ebrN21fQ91by8wfyVhKPNm0fOQ2ng708qvjl7n7ho0J/S02FougYkshFVu0u/ovhK7lJy9foH8kxO+9aSuRqOTx17VFqJvcWRw52cfHfDtM125NZQk1lSVmf1MFCjT37RfuuY4v3HNdyuN8X8XmhD5n4sA7dzM0HibbMdnu9u1rNJG6MoxTP3/3hjwcVsuMgTkzYVhgX3u8i7/+sebSBNio3yRNhKPmGIwEsD975SJbi7Lj1ga+eXsRP3jhPMe7h9i5Lo++kQkKs+Ot5Bs2x988Wi2Cr//eXoYnwqbIbC7MNuegSj3ZnO4doW9Y+8nvPN1nhvMfbj/D9uJc/I+8StU169henMsPXzjPR//zOR7/1J1J/47pZMWLlC5ADUAHMdbQFOqAxin7GgE/mqgtfwyu4qokmUAZCCG4rayI28qKpm2zlLjsVt5z86Qg2KyC/XvWma9jn2cyxhxsf0yU3L5tRVgtgpM9w+ZC35LCbApz7PSlKFKGC/eyHlr+/z1xAoBbtk3+vBguTaOUxgtn+3nvzZvi+nmz7lb+9fErmkgNh+YUyDFTmxJPNj95+QLZuiB3nOrjkh5d+KMXz/PO6zYAcLp3mO3Fufzns9q6rRNXhk2RGgtF+MaTJ/ngrVsWLfpxLiy/LbdApJQBKWWdlLJphmZGQEQs7fp+0MQtWXl5hUKxSijItpvisGVNNiWFWZy8MsLRCwNsKHCR47ThyXGaYmNwqmdYm3uLzpxn0/jRNxgYC7OhwBW3ds+ck8qbtJzeuiu+rtTmwmy2F+fyPy9d4HxwlNe7B7lmw8Ks6Os3FRAcCZmlQo6c7OPklWE2ubMYC0V57DWtVpcRYGJ8BtGYXIHNHWf50k9e5eAjry5oLPNlxYvUbOhFEb1MsbKMyD3dEguQRJSUq0+hWF3sKM6lIMtOvsvO1jU5BK4M83RXj2mtFuU4Etx9H/zXZ2h8PEDgytCMfV8aGDOfG+67Gze72ZDvMi0Pw2orzp8UqXdevyGhr/sqS2g/1Yf/kVeREqqnuDjny63eyZ+3zYVZPH8myGgowjuu0yJQjbVYhkgZ2TEGRieT4La0nwHgm0+f4o3uwQWNZz6sepFCF58Zwsm9U8VIF66E4AuFQrGy+eC+Ldz/lm0AXLexgGMXBugZnjCjHwtzHHHuvnAkylm99tOZ3plrQF0cGGPfNg9/d9+NfPZdWmj7DSXagvvr9TD6XCMEPdfJrnV5/PV7rkvq8n1fxWYcVgvff/48t3o9cdG1qbBtTY5pvf3Z/h3m/jt3rY0LOjnVoy1qNoJHjLIob3QP8sLZfu6t3IyU8PqlmQV7MbkaRCoxPC059wsh6vUAjGrg/iUck0KhWAbedcMG/kQPaPj9N281LZzpLKkn9HRNoKUrMpBS8trFQS72T1pPF/vH2OjO4j03b+Ztu4q5ffsac62ckc3DECSb1cJPP35H0iwsoIX0/+MHbubTv7mbB957w4LftzG/CdpnYFC2NteMlLQIONWrrd0aD2ulUAyRMsLbjfFOdYkuJSs+cGKxmFJefkYrSghRC9QClJaWLvHIFArFUlCU6+Rzd+/hV8evmMsU1he4GBwL0z04hlUI/uLhl9hQ4OLiwFicSP381W7+8BvtCAGP/Nkd7FyXS/fgmOnGK8i2860P7zPbby/WhOD0PIodGgK3WPzRnWVUbikk22Gj5SO30dx+lvX5Lnaty+flcwPcVOLm5XMD/OzoJfOcAV2kng70sG1NDnv0ubFUg0tS4WqwpABzbmpRkFI2SSkrpZSVa9eunf2EaRjoP7NYQ1IoFClw395S/ukD5eZrYyHsT166yB9/u5MrQ+M0/m4FO4pz4wTGWN8kpVZ4sHd4glBEsj7flfQ6t+jr2q7ftHz5EXevz+d3b9sKaBlm/NU3YLEIrtuUj0Vo6/smIlH+9D8nU2j1j4YIR6I8E+jlVm8RDpuFPKeN3hTXkqXC1WBJBfRHDzFroWJEK5BwRhr4u+++j5/2v84PPvgkDuf0iVAVCkX62LUuj40FLv7vD14B4O/uu5EbNrsp9WSbpTUAnjsT5KYSN/2jIZ7q6qEoV1tsO51IbV2Tw9MH9pvRhZnEb99SSsWWQvZsyOeGTQWc6Rvh9UtDPPZaNwNjIY5eGGBwPGy6C6fO2y01q96S0gMmAiQPMQ/qi3nTzr4t+zlnhf969M+X4/IKhSIJQgju0dctffZde8x1YiX6YlgpJeFIlJfO9nNTiZtbvR4efbWbjz/0Ap4cR1yewamsL3BlZKkMl93KDZvd2K0W3rR9DfftLeVzd19DUY6D/tEQJ65owRR79EwWhTkOepNkZl8qVr1I6bShpUSKpZxljOB7U+UfcxtZfO3SkwROPLpcw1AoFFP4uG8nrR+/gw+/xWvu21qUw8hEhH/5RRdff+IEo6EIN5YUcKuezNa3p5inDty14Ci8TKIgy07/aJgLenDIBn3ezpOd+oLnVFht7j4PyaP5GoBmIHbBb52+LRt/uf8r/M7PPkzdY3/GV6N/y/ayty/ncBQKBVpKq6m1yN598yZaj17iwZ++BmiZ7vdu9bCxIIv8LC2b/UzZQ1Yi+Vl2jl3QIhjznDYzErIwx5HWEPQVL1L63NIBtAW7bsAvhKgCWqWULWAmjG0QQviBI3pb/3K5+gw2b76Vxtsf4CO//jQf/OUn+PM39vOeu76E3b567sYUitVAQZadb/7hLZzsGcEqBAXZdrO21WLW4Mok8l12BkZDXOgfZYN7cq7Nk+1IOQlvKqx4kdLnnBrm0C42xDxj2LXzbv4zfzOfeeTD/NXFn9P4zVv4wNq9VN/xBQoKVHi7QpEpCCHYtmZ5k62mk4IsO4PjYc4FR1lfMFlJoDDHwchEhLFQZNHqxs3E6rJPVyjr19/Eod99mn/Zcz9eazZ/39vO/offySe/dQePPulnYjx9KUgUCoUCNHcfaNklNsRELXr00iPpsqZWvCW1WrBYbbzllj/lLbf8Ka+9/iOan/sXfjZymp8e/xY5r32TvdZ89q29gX073s1279sRFnV/oVAolg7DnTkRjrK+YFKkCvXaVn3DIbPo4lKiRCoD2bXzbj67824aQiM889y/8WjgRzwzcp5fdD8B3U/g+dUn2WcvYt+6Cm7ZXUNJyW3LPWSFQrHK2Fw4KUAbChItqXSlRlIilcHY7dncfstHuf2WjwJw/nw7zxw7zDMXj/Ds+GV+cr4VzreyKQK3uNZxg3s71625nrINldg9ZZBdBBlcFlqhUGQu+2LqYMUuQs7S56GMysNLjRKpFcTGjZW8Z2Ml7wFkNMqJU4/xzGvf45nLz/Ho+EW+d/kSXH4CxyuSknCIkghstmaz2bWGktxNbC3cwca112DzeMFdClmFy/2WFApFhiKEoOUjt3Hg4ZfMBLkAxkxDRM5cX2uxUCK1QhEWC95t+/Fu289vo2VlPnPlKC+f/TXHLr/IqcEznB27wjPhIUajF2HgIgx0YDspKQ2F2RoKsTVqYWtWMdsKd7B1fTnuzbdA8TXgWtoy5QqFYmVQudVD6yfeGrfPontn5EoQKSHETUxTwVZK+fOF9K2YH0IIStdeS+naa3lnzH4pJT1jPZwZOM3JnmOcvPwyJ4NdnBy5wC/BNO9OAAAgAElEQVQn+gkzAIMdMNiB+9WvsTUUZosli615JWxzl7EtfxslhTuw562HnDXa5rh6wnAVCkU8RmqnSDQ910tJpIQQ7wMOGy+TNJHA0gfQK2ZFCMGarDWsyVrDzevK446Fo2HOD53nZP8JTnS/wMnulzg1cIInx3r4/vgpuHQKLoFNSkpCYbaEQmwOh9koLWyy5bLJ4Wazs4ic7CLI9mjuw6xCyNKfZ3vAVQD2LLBng82lbSoyUaFYsRiWVKa7+/zACbRFtMuatUGROjaLjdL8UkrzS7mjJN6kH5oY4lTPqwQuv0Sg93UCAyc4PXKJZyaCjEqjpHQQZJC8oS48/RGKwmE8kQieSISCaBSnlLiiEoeUOGM3iwOn1YHD6sRhdeC0unDYnDhtLpzWLBz2LBy2LKz2LE3gDHGzu8CWBTbn5H67/tqWpR93Te63OsBiBYsdLDaw2lUgiUKxQAxLKhrNbJHyAtVSyocXczCKzCHXkcu1Gyq5dkN8Xl4pJcHxIOeGznF26CznBs/RPdJN71gvvaM9nBy9QsdYD/0Tg0SZy5c4DAzpGxAFxrXNJiUOCU4pccio/ijJjkqyZZSc6OTz7Kgkx3iMRsmWkmz9MS8a1bZIlBwEFlO0bNpjstczHbPYdPGLeRRWzUI0nxvHrPrr2HbTvBZT+7RMaWebco3YdgsYixJuxTywGpZUhotUG9PMRa1UhBDlaJnS3cBeoGG5c/tlIkIICl2FFLoKuW7NdTO2DUVDTEQmGAuPaY+R+MeJyATjkXHzcTwyTigaMp/HHQ+PMxEeZTw8wkhomJHQCFdCw4yERxkJjzIcGSMkZw+JFUCusJNvsZMrrORgJQuBCwsuBFmAS0KWFGRJSZYM4ZITZEUlrlAUl4zgjEZxRScfXdEozmgEVySCPRpGyAhEoxANg4xoj9HwbENbPsRU4VtMwUyD0AqL1lYYm3XyusISv5n7Fnhe3LnTnbc6xd94W5nu7qsDDgshCoBHiSkmaCClPLmAcaUVPUltpZSySX/tA1qBsmUd2ArHbrFjt9jJsacn0CIUCTESHmEkNMJwaJiRsPY4FBpiaGKIgYkBBicGGQoNMTgxyMD4AKPhUQbDo1wKjzAWHmM0PMpYRHuMSn1mWOjbjFNpVsCKwInL5sJp1R5dVu250+bEZXFq+61OXFY7TosDl9WB02LHabGTpT+6hE3fZ8UlbNpmseLEhlNYcAkrTiHIwopNgiAK0UiMKEaSvzYF0zg29bzwFHGNbRee0naqCBvXmJhlLNOdN6XdqkBMEa6pwmhNIoLJXlsnRS+pGFonj00nstOKrtHvTMIs4sbhnojwMdsp9rz+BGz537B255J+igtx95UBB2dosyiBE0KIarTihAm1n2Ksn4A+pkCydnPAiza/ZpTyaAe8Qgi3nsBWsQKwW+0UWAsocC68RLeUklA0xGh41NzGwmOMR8YZi4wxHtYezX0zHdP3jUfGGRrvTXo8nKKlZREWXfhcmhBaXbhsLuwWOzaLTduEbfK5vpnHzWOuhON2i12fM9TmDl1Wl/naac4nOuP2Ga8tYgHBMVImEdAZhNbYopGY1/rzaDTJvkj8eXH7jEeZuG9O58VeMzLN2GYba7LzZOIx87gx1lDi8YTPYGofMsn7mnrN+PHmyigfs0l4HSh/c8aKlF9/fBDoWqSxJKBbNIeAmiTHjHIbVTH7moUQgfm66aSUnXp5D4NKNGFUAnWVIoTAYXXgsDoWRfRmIxKNmII1Fh6LE7ZURDEUDRGOhglHw4xERsznsfvD0TBhmbhPzmkucWbsFnuCgBkWZrLXxvNsWza5jlxy7Dnk2nPJtmeTa88l167ty3HkkG3PXpgIKhbElcFx9v51K1+85xp+Z+fWJb9eqiJVDtRKKf91MQdjoAtQA9AB9E7TrA5onLKvEU1AE0RtNqYIWx1w/3z7UChSxWqxkm3JJjsDaolFohFNvCKhhDnD2HlCY25xLsfGI+OmmI5HxhmcGDTbGcdGI6Nztihz7DnkOfLId+RT4CxIeCxwFuBxevBkeVifvZ51OeuUsC0SWnSfICItaVlOkqpIdTK9eCwYXTDqAIQQ09WKqiZRpNrR5pLQz61l5nml1qnuQf2ch4yCiQrF1YbVYsWKFafVSS65ab12KBIy5xGHQ8MJz2PnGAcnBhmYGKB/vJ9TA6cYGB+gf6Kf8ch4Qr9Ztiy25m9lR+EOrim6hmuLrmVn4c6MuClYaayU6L770QIn+oAOKeXAIo5pVvRABy9ThFKvwIsQwiulDBiBEPPo10fq81oKhWKB2K123FY3bpc75T7GwmMEx4P0jfXRO9bL+eHzBIIBTvSf4MnzT/KDrh+YbddkraEkr4SSvBI25202n5fklVDoLESs0gi9hWAYpNEMj+57FC1Uuw1I9oeUUsqlzAvo0S8y3ZyRl3kuMtaDMHr1Cr4IIaqVNaVQrDxcNhfrbetZn7M+4ZiUku6Rbo72HOWN4BucGTzDmcEzPHPhmTjxAsiz57G9cDtl7jI2525mc95mNuduZlPuJgqcBVetgBmWVKaL1CFYhNnV1En9NisJ+hzYo4A75osXAJKKlO4SrAUoLVUl3hWKlYIQgnU561iXs463lb4t7th4ZJxzg+dM4To5cJLjfcdpO9VGcDz+fjjHnsOm3E1syNnA+pz15rYuex3F2cWszVq7al2JKyJ3H/DX6XbxLSX6HNic61bobsQmgMrKyuUUa4VCsUg4rU68bi9etzfh2HBomLODZ7VMK/rjuaFzXBy+yPOXn6d/vD/hnFx7Lm6nG7fTTYGrgAJHAWuz1lKcXUxxTjHFWcUUZxfjcXlWlKBZVogl1SGE+JqU8suLOpp5otYxKRSKdJBjz2GXZxe7PLuSHh8JjXBx5CKXRy7TPdJN90g3l0cv0zfWR/94P/1j/ZzsP8mV0StJAztcVhdul5tCZ6EmbC43HpcHt1Pfpx8rdBWyLmcd+Y7lK6czaUlltkj9HHg/sFwiZcw3eYjJdqEHVMQeVygUiiUn256Nt8CLtyDRCotFSsnAxIApZN0j3fSN99E3pm3B8SB9432cvXKW4FiQwdBg0n4KnYVsyd+SsG0t2IrT6kx6zmJhMdIiZbJISSnrhBDtQohHgDop5alFHtds1w8KIQIkzk150BbhKpFSKBQZhxDCXMe1o3DHrO1DkRD9E/30jvUSHAvSO97LxaGLnBw4yamBUzx1/im+3/V9s73D4qB8XTlv3vhmbtt4GzsLdy56gIcWQZ3hRQ+FEEbodzkQWIboPtAiCyvR1mwZlOv7FQqFYsVjt9rNenDTMRIa4fTgaU4OnOTFyy/y1Pmn+HLHl6EDCpwF7CzcaW7eAq8Z2GG32lMel1WIjE8wm87oPg/Jo/kagGYm8+2BtgC4Lh2DUigUikwg257Nbs9udnt2846t7wDg4vBFnjr/FC9cfoHjfcd5+PjDjIZH487zuDxsyt3EbRtvo3pHNRtyN8z5mhaLSFt0n1iIySaEuAstBZGxLql5McrG63NLB/R+q/W+29AyRLTEtCsH7gOO6G07070Qt7KyUra3t6fzkgqFQjEvojLK2cGznBw4qQV3jHZzeeQyXcEunr/8PDZh48G3PshdpXfNqb89n3uE371tC3/xzj0pj0kI0SGlrJytXcouOSHEV9GsliCaiOwFaoUQzVLK96faL5iLdKdLhxTbrpN4d59CoVAopmARFrMK91TOD53nU49/ik/84hN84U1f4J7t98zan9Ui0hY4kVJ2QCHE/WgCVSWl9EgpK6WUHjSr5l4hxJ8v5iAVCoVCsTRszN3IobcfYu/6vXz2ic/y8PHZC64Lkb7ovlRT2NYB9VLKR2N36q64TwMfWejAFAqFQpEesu3Z/PP+f+bNm97M55/8fEKKqKlYLSJt0X2pilQ507vZutDmhxQKhUKxQnBYHfz9nX/PLRtu4XNPfI7vv/H9adumM7ovVZFqQwtoSMZ9qHkihUKhWHG4bC7+8a5/ZN/6fXz2ic/y/576f0kzZKQzui/VwIlPA+36+qgmJsu31wHvA6qmP1WhUCgUmUqWLYt/2v9PfKXzK3zj6Dc4P3Ser9z1FRxWh9nGKgTRTJ6T0qPq7kVLjdQB9KFZT3XARxYjDF2hUCgUy4PD6uCTez/J52/7PE+cf4KvvfC1uOMWQcYv5jWCJFr0QoHb0KypdillYjpghUKhUKw43rfzfXR2d/L1l7/OWza/hZuLbwY0d19GW1KxSCnbpJSHpJSPKoFSKBSK1UXDLQ1syt3Exx/7OEd7jgJadF+6SnUsSKSEEPnTbYs1QIVCoVAsH/mOfL5y11ewCisf+PEHeOz0Y3p0X3qun+pi3puFEMfR5qKm21YsQojG5R6DQqFQZApl7jIevudh9nj20PCrBqK2S2lz9y0kwWwR8CDauqhVgz7HVotKVKtQKBQmBc4CvnLXV7jnv+9hKOdhItH6tFw3VZEqB2qllP+6mINJhhCiGq1GVELiWD3BbCWTIfCBhSSY1RPb9hJTSFGhUCgUGmuz11J3Yx1/0/439Ee7gIolv2aqItWG9mO+pOhWzSG0TOtTj3kBv5SyKmZfsxAisICihz4pZctiFwlTKBSK1ULlei1x+Xia4uRSDZxoAPxCiPcsRZCEEMKrzwt5mV4M64Cpc0eNgD/Fa/pQBRMVCoViRiy6bERlJC3XS9WSehStEGELkKw88YIq8+qWUJ3e93QlO6pJFKl2oNV4IYSoBcpmuFSrlLLNcPPpJUIUCoVCMQ0WYYhU5gdOpKsybwK6qCRYWVLKoBACIYRXShmQUjYl7yEBH+ARQhgFuNy6wLUtwHWoUCgUqw7DKInK9CTvS0mkpJSzFiRcYjz6OKazfIxKwXMittovaCHo8xA4hUKhuGqYdPelR6QWnHFimXAvRadCCLcQol5/Xq8HZyRrVyuEaBdCtF++fHkphqJQKBQZicWiyUZEiVT6kVIGpZQHpZRCf0xqjUkpm/RqxJVr165N9zAVCoVi2TAsKalEanb0uSmFQqFQpAkzcAIlUjNhWDie2J0xoqWCHRQKhWIJmAycWAEJZpcLPWAiQOLclActO4USKYVCoVgCrMIKqMCJudCGlhIplnLUglyFQqFYMibXSSmRMvCQPJqvgcR0SXX6foVCoVAsAQLN3ZeuwImUs0IsJfrc0gG09U5utBRMVWgZIlrAXLjbIITwA0f0tn7l6lMoFIqlI92W1KwiJYSIMv/sEgtNixRkDhaRlLIT6Ez1OgqFQqGYH2bgRJqi++YiJA+yjCmQFAqFQpE5GIETMlNy92VACiSFQqFQZAjpTjC7EgInFAqFQpEhTLr7MrtUBwBCiPeiBSwkIKX8m4X0rVAoFIrMYzItUoa4+5IhhNiGVrupUN8lgdiiUhJQIqVQKBSrjJWyTqoROAFUSCktQD/aQloP2mLausUZnkKhUCgyCUOkZJri6VIVKR9QL6V8Tn8dADN03I8SKYVCoViVTIpUZltSQaAg5nWAyRRFbjSrSqFQKBSrDCNwItNLdbQB7495fRho0AMp/GgiplAoFIpVRrrrSaUa3ddAjEtPStkihKgDWtAE6sOLMLa0I4Sojn09tay8QqFQXO1M1pPK4Og+KeUJ4NNT9lUJIbbpx1Ycetn4gC64buBRNNFVKBQKhY7m7hNpm5Na1ASzSyFQunUTlFImlOAQQpSjzYUF0NZrBZK1m8M13MABKWUhmAEgFQsaeBLGx8fp7e1lcHCQSCQ9C+EUCkU8VquVvLw8PB4PTqdzuYezIhGIjF8nlT9bGynlQCp9T7mODzhEYkkOhBBG1vOqmH3NQohACpnQK4GAIYhogR8ti5lRfXx8nNOnT1NYWMjWrVux2+3mBKRCoUgPUkpCoRADAwOcPn2a0tJSJVQpILBkvCUVZPaks9YU+zYEqAHoAHqnaVaHtl4rlka0wI0EUZsFL3rBRL0ESLt+7bJ59jMtvb29FBYWsmbNmsXqUqFQzBMhBA6Hw/w/7O3tZcOGDcs8qpWHQKQtd1+qIvURkotUJXA/UJvyiADdgqkDEEJMl+C2mkSRagdajRdCiFpmFppW3T0YQHMVBvXrB4UQXiGEd7GsqcHBQbZu3boYXSkUikUgPz+fkydPKpFKBSEQQhKNSiyWpfUIpRo40TTNoUNCiDdYRAskGfockpcpVpYuLhjiMsM4p5JMiBY1jD4SiWC32xezS4VCsQDsdruaG04RLQw9SkRKLCytSC1FFvQ2NCtnKfGAGeCQjKRJb6dDt5aCuvgZIpjK3NaMqDkohSJzUP+PqSOwADItLr+lKB/vY54ikQLuJeizBjgghOhCswSnndfS3Yi1AKWlpUswFIVCochchBAgJNE0xE6kGt03U0l5gRZ0sKLQraY5FXjU3YhNAJWVlapqsUKhuKowLKlIBltSM5WU70JLk7TkCCHcM7j8FAqFQrEECDRLKhJdepFKaU5KStkgpfz0NNshKWX/Yg90CsZckSd2pzGnRPJACMUy0tnZiR7UErdVVFRw8ODBhPYHDx6kqqoqSU9aXzU1NRQWFiKEoLCwkKqqKtrakq/jntq+rKyMmpoaOjs7F/U9KhRXC0Loc1KZKlLLjW49BUicm/KgZadQIpWh+P1+urq66OrqorW1FZ/Px5e+9CXKysoIBieN4iNHjiQVnaamJioqKujs7OTAgQM0Nzdz4MABent7aWycuiJBE7up7aurq+ns7GT//v1L+l4VitVKOt19SCln3ID8VLbZ+p3rhuY+9CXZ3wjUTtlXDTQv1rXnslVUVMi5cPTo0Tm1W610dHRIQDY2NiYc6+rqkoCsra2dsY/W1lYJSJ/Pl/R4X1/fvNq3trbOcfSK1crV/n+ZKrd+661y59//b3lpYDTlPoB2OYff2LnMSc0lu0QyUs44MQUPyaP5GoBm9AAGnTpUwcUVh9frpby8nPb29hnb1dXV4fV6aW1tTXrc7XYntHe73dO29/l8qQ1YobjKsWDRUsxmSHRfsuwSbrT0Qw3EL3rdDnwKqF/IoIyEr2ih7G7AL4SoQssQ0QLmwt0GIYQfOKK39Uvl6ltxBAIBOjs7qa+f/mvT0tJCIBBI6tJbjPYKhWLuGHNSGRHdJ5NkbRBCfBXN1favSY5dAYoWMiipzTnNGg4upewE1Oz3CiUQCNDW1kZDQwM+nw+/328ea2ho4ODBg2am5SNHjgBw7733zqnv+bZXKBRzx4KxTioDRGoa7iMmR94UOoGvoVlCijnwhR++wtHzC04av6RcszGf//u/rl1wP3V1ddTVxXtkGxsbqa2dOd1jIKAZyFNdeovVXqFQzB3TkspgkQIt48LDSfbXMCU0XKEwqK+v57777gO0DNSdnZ00NDTQ3NxMc3PztKISG/k3F+bbXqFQzB2LEBi5+5aaVEXqAeABIcQjaFF2J9CEqQ5NpPwznKuYwmJYKCuFsrIyysvLzdc+n4/a2loqKirYv38/HR3Jk5WUl5fT1tZGMBick3U0XXsjHN2go6MjbjwKhWJ2BBYQZlT1kpLqYt6DaOXjbwG+i5YGqQ2oAhqklMrVp5gzbrcbv99PZ2fntAts9+7dC8Dhw3NLZjJd++bmZjo6OuLmvxQKxfywmO6+NFwr1ROllAellB60GlL3AhVSSo+U8sFFG53iqsGwdnp7k9e4rK6uxuv1zllcpmtvhLt7vUudA1mhWL2kc05qwRknpJSdUsrvSimfW4wBKa5OGhq0YM6Z1i41NjYSCASoqUmeoH6qFWa0ny69kkKhSA0zui8T5qSEED8D7pJS2mL2zZQFHUDGtlcoDILBoBl5FwwGaW9vx+/3EwgEpl10a+Dz+WhubqampoaKigrq6urweDzmue3t7fT19cW19/v9NDQ0UFZWRkNDAx6PZ9qUSwqFYm5YhAVBNGOi+9pIrFI7UxZ0hWJaGhoaTKsJJt1vra2tc3LBVVdX09XVZfZjBEb4fL6kQRf19fVmfkCjvXFNv9+vgiYUihQQQgucSEd0n0hHdMZqprKyUs6Wzgfg2LFj7NmzJw0jUigUc0X9X6bG3d+toetilG/dfYiKLamtOBJCdEgpK2drl2rRw62AW0r5fMy+TwJ7gZ9JKf8tlX6XEyFEOVoQCGipmFpUiiWFQqFIxIIFiKQlui/VeaNGtFIZfwTmvNV+4DmgSQhRJqX8i8UZYtrw6aH1AAghGlHJahUKhSIBbTFvZkf3VaGnRRJC3Az40NZHVaKlQ0oefpUCQohqIUTSkC8hRLkQolYI4TMeF3CpupiiiQqFQqGYBouwgIhmRnTfDBiuMB9aEIWRiLYLLSP5gtFF5xBJRE8IYWQ9r4rZ1yyECKTopvMDJ4QQxqz+rAluFQqF4mrEsgJy9wWAe4UQAbRks51SSiNDqpcFlm/XBagBLZNF8tWdmituah2GRjSxmbclJ6Vs0i0pw8WXLKpRoVAornosaYzuS1WkPg0cRhMSgWZNGdShFSNMGd0SqgOIsWymUk2iSLUTk51dCFELlM1wqVYpZZvetl6fkzqon9c6y7kKhUJxVWJYUumIDk9JpKSULUKIMqAcaJNS9scc9qMJ2JKhWzxeplhZeiFEhBBeKWUgWS2safrzEVOXSreqyoQQ5XrNKoVCoVDoaIt5Mzx3H9o8VCXwkBAiH0AIUYAmXEstrx4wiyMmY75zYr1o445DCZRCoVAkYrNYcDkE2Q7r0l8rlZOEEPuBn6GFnN+MJhoDUsp+3cJ6APjjRRtlIosahSel7BRCeHU3n9H/Q4t5DYVCoVgt5Djt7FyXy5u3r1nya6U6J+VHW+x6n57HL5YmtLmipRSpRUdK2TLXtrqY1QKUlpYu2ZgUCoUiE7EKK1GZBl8fqbv7ypne0uhjkS2d6ViudU1SyiYpZaWUsnLt2rXLMQSFQqFYNoQQyDSlb01VpDqZTCE0lTq08O2lxAhxj0saFSNaKp2RQqFQLBEWLGmzpBZSPv6wEHpuDCgQQtyE5gb0ARWLNL6k6FF8ARItNg8QVDn3FAqFYumwiAwXKT0E/V60+SeBZlkJtMWv98Ymnl1C2tCsudgIvHKW3opTKBSKqxohRMbPSSGlbNHLx1cQXz7+u4s2Og0Pyee4GkjMLFGHSmeU8XR2dlJTU0NhYSFCCMrKyqipqYmrrNvZ2Ym+5m3WAoVGu4MHzfzAcefHbhUVFXHtDNra2pK2j92CwckVD8FgkJqaGsrKysz3UFdXZxZ0VCyMpqYm87OtqKigpSUxrsmoujxTG8XSkM7AiQVXz9XLxseVjhdC3LQQa0qfWzqAtt7JDfiFEFVoGSJa9OsGhRANQgg/cERv61euvszm4MGDNDQ04PV6OXDgAF6vlyNHjtDS0sL+/fvjKusaNDc3T1tWfrYfJr/fT3V1NRBfCbixsZGOjg7c7vj7n9ra2mnLzRttA4EAZWVleL1e6urqcLvddHV1meLX2Dg1EYpiPjQ1NdHQ0MCBAwdwu900NjZSU1NDY2MjtbXaKpFgMEhFRQVer5fm5maOHDlCTU0Nzc3N5t9bsXSkM3ACKeWsG5CPFs3XAxwH/nCGdoeByFz6XQ1bRUWFnAtHjx6dU7vVTGtrqwSkz+eb9rhBR0eH2dbtdk/bp8/nkz6fTwLS7/cnnN/Y2JhwTl9fn3S73bK6ujphbMnaJ7um1+tN2m9HR8es5y8H9fX1squra0X063a7474LUkpZXl4e9z2or6+Xbrdb9vX1mfuqq6tn/K4kQ/1fpkbDLxvkO1resaA+gHY5h9/Yubr7OtHKczyKNvfUJIT489gGQogvoYWf+9By+ykUcRhWR2tra9LjyaylmpoagsFgUpefsb+mZn75hI1y87O5Eaejvb096VjdbndGlqMPBoMcPHgwzp2aqf0Gg0G8Xm/C53vgwAGCwaB5rZaWFmpra+Ms4altFEuHBUvmhKALIT6lP90mpbxXSrkdeB4tEWueEOLDQogetPmgT0ttXurBJRyzYgXS0tJCIBDA7/fP6zyPx0N5eXlSF9rhw1qKyOlcgTPR2dmJx5Na2WuPx0N7e3tK5ypmxu128+ijjybdD5qrNRgMmi7XWLxeLRua+tssPZkWOHEf8ICMTyJrZD8PokX4HUITMSVOiqQcOXIEgHvvvXfe5953331J556Muaq5io1heVVUVBAIBGhoSC3Gxu/309nZSVVVVcrW2HQ0NDRQWFgIaJZnWVmZGbDR0tJCRUWFGSgQe+2WlhYz0EAIQVVVFYFAwAxQAc0qFUKYr4G4NrHnzTame+65Z8Z+F8LUeUJjnKDdkPT2anmlDVGael5sgItiaci0wIlytBIYsRzRH5vQrKd+FKnzk0/DxZeWexQzs/56+M0HUj7d+JFJ9gM0G9XV1TQ0NNDW1mZaTYbgzBakUFdXR11dXcJ+v99vTsLPpX1tba15rerqahobG6mrqzOForq6Gr/fn/DDmSoVFRXmeNxuNy0tLdTU1FBfX8+BAwd46KGHqKqqoqurC8A8dt999xEIBGhtbSUYDHLo0CHq6uqoqqqisbExweo03tOhQ4dwu900NDRQUVGRNIAldky1tbW0t7dP2+98Ix39fv+MrlK/34/P58Ptds/ab09Pz5yvq0gNi7CkpUwHzD26L+7WRGqJZAEalUAp5sJC7m69Xi9erzfux9Bw9c1mmRk/3KD9IDc1NdHX1zetWMa2nzqGWGpra6mtraWtrY3m5mYOHz5MWVkZra2tKbkfYwkGg3g8nri5u/vvv5/6+nrTXVpdXU1FRQWNjY3s3bsXwIyGKy8vj4twM8bu8XgS3sdU9+uhQ4fMcO7YPpKNaaZ+FzPC8eDBg/T29ibMZSqLafkQQhCRkbRca8Eh6IpFYAEWykqhvLyctrY2gsFgnEBUVFTETXR3dHQkvaOurq6OW99kuPrcbveMP1ZlZWVmf36/n6amJg4fPpzUiprafi74fD58PlY9wjgAABgBSURBVB9+v5+KigpqamqSWiHzJfZHvrOz0wxSmLrGy+PxcODAAQD279/Pfffdh8/nSzmAwzgvmbWyHKH1LS0tNDQ00NHRYQqh8f0x3H5TKSoqStv4rlYyKnBC5316gIS5oaVDunfq/pjjCoWJcbdvWEAGzc3NdHR0zBpQYbjgDKFra2tL6pabCbfbTX19fcpzUbP13dDQsGjRZbGWiSEYHR0ddHV1xW3Nzc1xEZOGu266tV7JMFyJFRUVM84rLZYrc650dnZy//33J9y4GHOQU29OFuJSVsyPTEyLNF1AxHSh5hL41/kPR7Faqa6uxuv1JswFGT98s80zGC6/5uZm88cplUWbBw4cMBeLzjfScDaW6kc8tt/pruHz+ejo6CAYDJrv7+DBg9TX18/YtxEo0dDQQGWlljPamHtaCAudkwoEAuzfv5/m5uYEq9DtduP1es35OAPj5mCh7lbF7GSaSC1psljF1UNjYyNVVVVUVVVNu1ZqJurq6uKCF1LB7XZz4MCBuIwGqYzD7/cnnGvsW+y1UuXl5aallmxeJnYchrX40EMPmRGV00U/BgIBc04tNivHXJkpqnIhrsFgMEhVVRWHDh2aVnCqq6tpamqKu05jY6N5M6NYWjIqcEJqaY8UigVjzN00NDRQVlZGQ0MDHo+HI0eOzCmU24jya2pqorm5OeVx1NfX86UvfYn7778/oZ+ZxNMIdz98+DBNTU1UV1ezd+9eenp6zHVgsf0ZFsqBAwcWLFyHDh2ipqaGmpoa6urqCAaD5o9yRUUFfr+furo6vF4vgUCAzs5O0x1qiJjxg97Y2Ehrayter9cUP6PdfKzL6fpdKEa6o97eXpqamuKO+Xw+M6XWwYMHqaqqoq6uzvwOLeR7oZg76QycWPa0QsuxoeUAnLrPC9SjZcyoB9xz6UulRZo/HR0dZgobQHq9XlldXZ00LVJzc3PcuV6vV2pf20n6+vrmlRZJSin9fr8EzDRGRlqkmTZjLH19fdLv95upetxut/T5fAkpkWYbQzLq6+sT3l9sf0YKKLfbLWtra2VfX5/s6+uT9fX15mfj9XplfX190n6nHmtubjb/DuXl5bKxsTHhs5xpTNP1myrG33K6Lfaz7OrqMj8Pr9eb8F2ZC+r/MjUefPZBufdbexfUB3NMiyRkmky2TEAI4UMTo0YppZhyrENKWaE/dwOHpJSz5tuprKyUc1nhfuzYMfbs2ZPawBUrlrKysqTzKorMQP1fpsbfdvwt3zn2Hdp/J/XsHvpv7nTFc00yIgRdCFGNVqwwwecjhChHqxsVQBOYQLJ2c8E4TwgR5zDXr9Eb0y6oC5pCsWCUQClWGyuhMu+ioYvBIRJrQyGEMMpvVMXsaxZCBOTiluTwMmXBMtArhCiXUqpslYqUCAaDix5BqFBkAukMnEi56OFCEUJ4dYvGS4wVM4U6YGqYUCNamfrFJLVMowrFDLjdblXbSLEquSoyTuiWUB2AEGK61ZXVJIpUO2CGEAkhaoEypqd1Du7BXhKr/yrhUigUiiRYhJZxQkqJniJvyVh2d9906MELCVaWPl+EEMIrpQxIKZuS9zAvAiQRJeXqUygUikQsQnPCSSSCpRWpZXP3zQEPaKI0zfFFW7E3VYz0ubDFrcGgUCgUqwSLLh3pCJ7IWEuKRPfbgtGj+Hz6cz/xrsD7hRD1aFWIy4H7Z+inFqgFKC0tXexhKhQKRUZjWlJpCJ7IZJFadHSLqRM4OMMxmMWK0l2MTaCtk1rkYSoUCkVGY8xDRWQEO/YlvVYmu/sAc25KoVAoFBmCYUmlw92XySJlrIOKC2iIEa3FXCelUCgUijliFVaAtNSUyliR0gMmAiQPDQ8u8mJehUKhUMwRI6LvarekQJsbmprbqRwVeadQKBTLxtXo7vOQPJqvgcR0SXX6foVCoVAsA0bgxKoOQdfnlg6grXdyA34hRBVaWHgLmAt3G/Rw8SN6W79y9SkUCsXycVVYUlLKoJSyQUpZI6UUUsoyKWWdIVAx7Tr1di1SyoOpZkBXZA6dnZ3U1NRQWFiIEIKysjJqamrM8t9GGz2zyKwFEY12Bw9OriyIPT92q6ioiGtn0NbWlrR97BZbtTYYDFJTU0NZWZn5HuZbMl2hWKmkM3DiqlonpVh+Dh48SENDg1ld1ev1cuTIEVpaWti/fz99fX0J5zQ3N09bRrylpSXpfgO/3x9XGr29vR2/309jYyMdHR0JJeBra2upqqpK1pXZNhAIUFZWhtfrpa6uDrfbTVdXlyl+CymdrlCsBNLp7lv2KrkrfVOVeeeOUf3W5/NNe9zAqGrr8/mk2+2etk+fz2dWZ51rZd6+vj7pdrtldXV1wtjmUkXX5/NJr9ebtN+p1XmvJurr62VXV9eK6VdK9X+ZKi2vtcjr/v06eWHoQsp9MMfKvJkSOKG4CjCsjtbW1qTHk1lLNTU1BIPBpC4/Y39NzawFlONwu934fL5Z3YjT0d7ennSsbrf7qi1wGAwGOXjwYJzLNpP7VSyMq2JOSnF10dLSQiAQmHcRQI/HQ3l5eVIX2uHDh4Hk4jYbnZ2deDypVWPxeDz/f3tnsNs4cq3hv+I7aNyFG2zlztYDUJsMEGQhdd8XGCmLu5kNNX4CS0AeQISfoMF5A8mbu8miLb2B1ECA7NKWN1lMNhYmg+wmYwneGoOTRVXRFEVJlKwWi9b/AYTtIlV1Do/Fw6o6dQo3N7tvm01I2TnkcB+dFDkInz59AgB89913W3/2/Pw8c+7JzlXldTa251Wv1zGdThGGu61kiKIIt7e3aDabO/fGVhGGId68eRPXb4M9snoS6Wv6/Xy71tg2AN27rVarcVDIcDhEvV6P60zqNxwO40ARpRSazSam02kcBAPonq9SKv4bwMI1yc9tkunbb79dWy8pjjhwgglmj4PobxH+cf+PosVYy+8qv0P4v7svT7MPpXSgQh6CIEAYhhiPx3GvyTqcTUEKnU4HnU5nqTyKIrTb7dzXt9vtuK0gCNDr9dDpdOKHeBAEiKIIvv/8HWRs5GAURWi1WgjDEK1WC3d3d/E1w+EQrVYL7XYbo9EIt7e36HQ6mEwmuQM36vU6gKdhWFtnt9vF5eUlPnz4gGazGbdrz52fn2M6nWI0GmE+n+Pq6gqdTgfNZhO9Xm+pZ2vlubq6gud5CMMQ9Xo9M0gmKVO73cbNzc3KeklxxD0pvOB1UuS4SIZvb4vv+/B9f+FBZYf6NvXM7EMV0A/Lfr+P2Wy20lkmr0/LkKTdbqPdbmM8HmMwGOD6+hrVahWj0WgvD9NkVKKNHJzP57HcFxcXC46z0WigVquh2Wyi1WptlGE+n6NSqSzMD15cXKDb7cZDskEQoF6vo9fr4d27dwCAy8vLeO7Nygc83Z9KpbJ0r9JDvFdXV6jX6xgOhwt1ZMm0rl5SHNxP6sh4Tg+lLNRqNYzH44UHLYCloazJZJIZfBAEwcL6JjvU53neWgdYrVbj+qIoQr/fx/X1dWYvKn19HhqNBhqNBqIoQr1eR6vVyuwhbEvSyVSr1YVz9j6me3z2fqwL2U+S7HHd3t7GQQrpdWSVSgWXl5cAgG+++Qbn5+exU9wF+7msNWUM3y8HDJwgLw77Jm57QJbBYIDJZLIxoMI+kO0DejweZw7LrcPzPHS73Z3nojbVHYYh5vP5XiLR1g2L2od7Vs/C9/3cC4qTn7efmUwmuLu7WzgGg8FCVKYdrlu1niwLO5RYr9fXziuxt1QO6KQ+M0qppRhopVRNKdVWSnWVUgOzhTzZE0EQwPf9JWfk+z5qtdrGh5Md8hsMBgvzQNtiewSfw1Ed6gFr28lyRtPpdCc5kp+x99oe1mE2Gg1MJhPMZjNEUYTxeJyZvSNNs9lEGIZoNpu4urrCx48ft5aPuAWd1GdCKdUwW783UuUegLci0heR7wH0AGQv5iE70+v1MJ1Ot3oDT2IDFUaj0U4OCtA9lMvLy3iOZ1c5sj4bRdFB1kq9fas3BkgPjQ2HwzjoYltqtVrcG0yT1tX2SGu1Why1uSrCcjqdYjwex4EqeV5Ikuy6TIB8XmzgxNGkRVJKBdB7RC3F8yqlatDbdUyhE8xOs67Lg/2cUio98O1DZ1a3Mbw3AHyllCd6XyuyB+zcTRiGqFarCMMQlUoFnz59yhXKbaP8+v0+BoPBznJ0u128f/8eFxcXS/WsWmhs5a9UKri+vka/30cQBHj37h1++eWXeB1Ysj4b5n55eblXx+V5XhxdeH9/H0fbhWGIIAh2Dty4urpCq9VCq9WKHXGv14Pv+6jX64iiCJ1OJx5StBGFVibgyXH2ej2MRqO4J2adn+d5W62VW1UvKZZDBk4UnlYIulczA9DIOOdDZ0VPlg0A+M9sU7LaSsuUpy6mRdqeyWQiQRCI53kCQHzflyAIMtMiDQaDhc/6vi/GfjGz2WyrtEgiIlEUCYA4jZFNi7TusLLMZjOJokhqtZp4niee50mj0VhKibRJhiy63e6Sfr1eTwDIbDZbKB+NRlKr1eJ7mNR/2zaSMts0U57nSbvdltlsJrPZTLrdbnz/fd+XbrebWW/63GAwiG1dq9VifZLyrpNpVb37gN/L3fjLT3+R3///7+XvP/995zqQMy2SkgMsxsrCzPmEACbmZ0dSPSS7RYckMqMrpRrm2u3HNJ7qEBFRa84PAHyQVEb2LN6+fSt5sg/88MMP+Prrr7cTlJSearWKwWBwtOmSXIffy93467/+ij99/BP+/H9/xh++/MNOdSilJiKS3tR2icKG+0TvCdUBAKXUqlnsAHp+KMkNEvNFZo6pitWM0s5vHaa+XA6KkDzQQZGXxiEDJ5yYk8rCBDP4AO6T5aI3QoRSyheRqYjkywWTr80GnjHnRUiS+Xy+da5CQsrAIQMnXI7uqwDaKa04v9d4XxOgcS9PwRW7hY8RYvA8b+coREJchj0pzfZJ3jZgHFHD/B7BDAWa+bGPADz7hgAdTZg55GeGBNsAcHZ2tm8xCSHEaZgW6TMhIrcAbgF8nyqfAsidXtkMMfYBHTixTxkJIcR1uJg3gZmbIoQQ4gh0Uhqb82VhyXnCaeVLUEYIIWSvWCd1iCVMzjopEzAxxfLcVAU6O0XpnFRRa9IIIcvw+7g7h9xPylknZRhDp0RKUjPlpeLk5ASPj49Fi0EIMTw+PuLk5KRoMUrJIQMnXHFSFWRH84UA0pklOqa8VJyenuLh4aFoMQghhoeHB5yenhYtRin5zW+OILrPzC1dQq938gBESqkmdFj4EIgX7oY2PZK5NirjUF+lUsFPP/0EAHj9+jW++OILJMLdCSEHQETw+PiIh4cHzGYzLiHZkaMIQTdzTht7RImw8VLz6tUrnJ2d4f7+Hj/++CN+/fXXokUi5Cg5OTnB6ekpzs7O8OrVq6LFKSWvX73GH7/6I7787y8/e1uFJZh9KeRNMEsIIeSJvAlmXZmTIoQQQpagkyKEEOIsdFKEEEKchU6KEEKIs9BJEUIIcRY6KUIIIc5CJ0UIIcRZ6KQIIYQ4CxfzPhOl1M8A/rnjx/8HwL/3KE6RUBc3oS5uQl2Ar0RkY8oKOqkCUUrd5FlxXQaoi5tQFzehLvnhcB8hhBBnoZMihBDiLHRSxdIvWoA9Ql3chLq4CXXJCeekCCGEOAt7UoQQQpyFTooQQoizFLYz7zGjlKoBeAtgCsAHMBWRcbFSrUcp5QFoAxgCuAdQAdABMErK7qpuSqkAwDxLljwyu6TXKl3KZiOjhw+gan72RGSYuqYUttmkS5lso5RqAGgC+AVan4mI9FPXHM4uIsLjgIcx1ihVNgDgFy1bDrklccwABGXQDUDDyNvYxR4u6ZVDl1LYCEAAoJb42wNwB6BdNttsoYvztjH/X41U2QRAtyi7cLjv8HQA9FJlPQBRAbJsSxPAGwBVEXkjqbdeOKabUspXSvWgvzD3Ky7LI3PheuXUBSiPjXwRubV/iMjcyJCUrRS2QT5dgHLYppNRNk6VH9Yuh/LQPOK3iTuk3iag37ykaNk2yO0j4+29LLoZ2bJ6Hxtldk2vNbqUwkamvQkAL0N+sfKVwTZb6FIW2wwARKmyCHrIL7ec+9SFPakDYsall96ERb95QSnlFyHXPiijbnlkLqNeq3BFF9Oeb45MymKbPLrkwQVdTHstEQlTxQGAD3nl3LcuDJw4LBXgyVgZ+NCTjK7imwliQOtyL09DFmXULY/MyHGNS3qVwkYi8iajuAEdEDK1D7Iy2GaTLkl5ymCbJEqpNoBbEfneFB38O0MndVi8ogV4BvcAIIsRSwOllC0ro255ZC6TXmW3UQfAe/N72W2T1AUomW2MM20CuneVOHVwu3C4j+RCROaSCkNFeQI+joIy28i8sd8n3thLS5YuZbONiAxFpAMgVEpNTDh5IdBJFYAZs30JTKGHMGJ9yqhbHpnLqJfBeRuZob2OiDQzzpXKNut0ycB525ghux6Aj8nyQ9qFTuqw2HHYSrIwYUyX5jYWUEp1M4rtxGhyjLlMuuWRuTR6ldhGEYBvUmVltU2WLmW2DaBD0D2zyPfgdqGTOiDmrWSK5THbCpYnWZ3BvB1GGVE59p9wWkbd8shcFr3KaiOz9itMT7KX0TardCmLbUxk3mzN0J5XhF3opA7PGDpVSJKaKXcS80/VyfjnakBH/tgvZel0Qz6ZnderjDYyczdRUmalVCPxMC+NbdbpUiLbeFjsCVmsPeyC5cPa5bmLv3hsvVjOw3K6kBHcT4sUJGXE0yLGdDoYJ3VDRhqavDK7ptcaXUpjIyNr2zy47NGAznlXKtvk1KUUtgHQxfLC5BESC3wPbRfuJ1UApjt9DuAT9FvKrTiQhHUTiSSav4X+J1x4czTXOKObGQO/NHIE0G+IY+gvTzIUeKPMReu1hS7O28joMltxeioi1cS1TttmS12ct42RoQ2dWHZTgtmD2IVOihBCiLNwTooQQoiz0EkRQghxFjopQgghzkInRQghxFnopAghhDgLnRQhhBBnoZMihBDiLHRShBBCnIVOirxITN40MavnnUUp1VVKjYqWo0iUUpFSilkFSCZ0UoQUyzvoPG+EkAzopMhRYd7a01smFNauiLRERB1aHkLKAp0UORpMMtAudJbqF98uIS8BOilCCCHOQidFjgKl1ABPWyoMTFDFLHVNoJSamHMTs122PRfZ65VSPaXUnVLKM7uZDsyOpqKUGiWH9Ta1mxU0oJSqmXqsHO3U+cjuoJq6LndPLW8dK+Rrm+u9FfUt3ENzj2wbd2bLiiyZrCwzc10jdX5r++S9H8Rd6KTIsXABoGl+70Dvk1O3J82DcwC9R1MLer+mBYdjrptA7zjaE72jaidVfwV6M7tc7aYxckxM+00AHwD0zNbkSTwjb8/U65u/t2EfdWTV9z5Vnz1a5rpVbVwZWULzd3z/n2EfUnYOteMjDx6HPKAj5gRAO1Hmm7KsXW1nSOw+asomtgxAZD472tBuLd3GhnYj/TVckKO3QpdGSpYgXQ9Su6qukTNXHWn5TFl71XVWxlRZN1EWmLLahjbsPevt0z48ynewJ0WOHjPE5QHomqEkMUNcdivwJJ2lChKIyK35desIQjN85UH3JpJ1jgHM8dQTsSR3Ob3btr091pHkJqO+ZBt2J9rKukpE71g7BvB2n/Yh5eO/ihaAEAewDqUO7QyS3Cf/kNR230A8FHVu6nlOeLv97FIbpmyhbtnDcNY+6shR365tTKGH7p5lH1Ju6KQISTiFbR9yJluEDz3cZHsRk9WfyCWHD+A2dc7HYi/lGPCh78nO9iHlh8N95Ji4zyo0Q3RzaEezwLoIMTNp3wAQikjf1JP1EM1sNwPrhBaGrExPzQYluMDKwI99Ye57A8B4V/uQlwF7UuRoEJG5UgoAOvaniNjIuwvoEHEb7eZBO4spVsxziMhUKTUHEJn65niKTMvbbvq6DnQ0XwU6ss/20oZmburQ3AE6xyD0HFEDOnBi75h7H+JJ5zl0pCCwg33Iy4A9KXJsfA/9oI2QGFITkSF0D8EDMIIOh54iw+mkuIAOAhiYOlf1djLbTSMifejQcxu+3YHuqaWDJg7FNbS8Vrff4vM5hQ/Q934E3fus2zmuZ9iHlBwlwuTDhBBC3IQ9KUIIIc5CJ0UIIcRZ6KQIIYQ4C50UIYQQZ6GTIoQQ4ix0UoQQQpyFTooQQoiz0EkRQghxFjopQgghzkInRQghxFn+A1m8NpchYcQKAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "res_all_bicg = []\n",
    "def bicg_cl(x):\n",
    "    res_all_bicg.append(np.linalg.norm(A.dot(x) - rhs))\n",
    "    \n",
    "sol = scipy.sparse.linalg.bicg(A, rhs, x0=np.zeros(n), callback=bicg_cl)\n",
    "plt.semilogy(res_all_bicg, label='BiCG')\n",
    "plt.semilogy(res_gmres_rst[:n], label='GMRES, restart=20')\n",
    "plt.semilogy(res_gmres, label='GMRES, no restart')\n",
    "plt.xlabel('Iteration number', fontsize=20)\n",
    "plt.ylabel('Residual norm', fontsize=20)\n",
    "plt.xticks(fontsize=20)\n",
    "plt.yticks(fontsize=20)\n",
    "plt.legend(fontsize=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "##  BiConjugate Gradients\n",
    "\n",
    "There are two options:\n",
    "\n",
    "1. Use $\\mathcal{K}(A^* A, A^* f)$ to generate the subspace. That leads to square of condition number\n",
    "2. Instead, use two Krylov subspaces $\\mathcal{K}(A)$ and $\\mathcal{K}(A^*)$ to generate two basises that are **biorthogonal** (so-called biorthogonal Lanczos).\n",
    "\n",
    "The goal is to compute the Petrov-Galerkin projection \n",
    "\n",
    "$$W^* A V \\widehat{x} = W^* f$$\n",
    "\n",
    "with columns $W$ from the Krylov subspace of $A^*$, $V$ from $A$ (cf. with CG case). \n",
    "\n",
    "That may lead to instabilities if we try to recompute the solutions in the efficient way. It is related to the **pivoting** (which we did not use in CG), and it is not naturally implemented here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Notes about BiCG\n",
    "\n",
    "A practical implementation of BiCG uses **two-sided Lanczos process**: generating Krylov subspace for $A$ and $A^{\\top}$\n",
    "\n",
    "In partucular\n",
    "\n",
    "1. $\\alpha_j = \\frac{(r_j, \\hat{r}_j)}{(Ap_j, \\hat{p}_j)}$\n",
    "2. $x_{j+1} = x_j + \\alpha_j p_j $\n",
    "3. $r_{j+1} = r_j - \\alpha_j Ap_j$\n",
    "4. $\\hat{r}_{j+1} = \\hat{r}_j - \\alpha_j A^{\\top}\\hat{p}_j$\n",
    "5. $\\beta_j = \\frac{(r_{j+1}, \\hat{r}_{j+1})}{(r_j, \\hat{r}_j)}$\n",
    "6. $p_{j+1} = r_{j+1} + \\beta_j p_j$\n",
    "7. $\\hat{p}_{j+1} = \\hat{r}_{j+1} - \\beta_j \\hat{p}_j$\n",
    "\n",
    "Now we move to the stable version of the BiCG method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## BiCGStab\n",
    "\n",
    "- BiCGStab is frequently used, and represent a **stabilized** version of BiCG. It has faster and smoother convergence than original BiCG method.\n",
    "\n",
    "- The formulas can be found, for example, [here](https://en.wikipedia.org/wiki/Biconjugate_gradient_stabilized_method)\n",
    "\n",
    "- It is a combination of BiCG step followed by GMRES(1) step in order to smooth the convergence.\n",
    "\n",
    "- For more details, please consult the book [\"Iterative Krylov Methods for Large Linear Systems\"](http://www.lmn.pub.ro/~daniel/ElectromagneticModelingDoctoral/Books/Numerical%20Methods/VanDerVorst2003%20Iterative%20Krylov%20Methods%20for%20Large%20Linear%20Systems.pdf) by H. Van-der Vorst.\n",
    "\n",
    "A short demo to compare \"Stabilized\" vs \"Non-stabilized\" versions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dimension of the linear system = 300\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaIAAAEcCAYAAACFy7BqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJztnXt8XHWZ/z/PTNOElraTtEi5tWWiCxVBmrRgEKFpJotowSJJi7sW/emSUERlF0noIrqCUibQ34+fcnFSL/xwdaGpgOyiSG5VkCi51Ct4S/CCLGtpkrZq6SV5fn885zvnzMyZycxkrsnzfr3Oa2bO+Z7v+Z45M+c5z+X7PMTMUBRFUZR84cn3ABRFUZTZjQoiRVEUJa+oIFIURVHyigoiRVEUJa+oIFIURVHyigoiRVEUJa+oIFIURVHyigoiRVEUJa+oIFIURVHyypx8D6AYWLJkCa9YsSLfw1AURSkqBgcHX2PmE6Zqp4IoCVasWIGBgYF8D0NRFKWoIKLfJ9NOTXOKoihKXlFBpCiKouQVFUSKoihKXlFBpCiKouQVFUSKoihKXlFBpCjK7KGtDejtjVzX2yvrlbyhgkhRlNnDmjXAxo22MOrtlc9r1uR3XLMcnUekKMrsobYW2LlThE9NDdDXJ59ra/M9slmNakSKoswuamuBdeuA//xP4IorVAgVACqIFEWZXfT2At/5jrzv6Ij1GSk5RwWRoiizB+MT+tCH5PMtt0T6jJS8oIJIUZTZQ3+/+IROPlk+n3mmfO7vz++4ZjkarKAoyuyhpUVee3rk9dAhYP169RPlGdWIFEWZ2bjNHXrxRXl9/fXcj0eJQQWRoigzG7e5Q08+Ke9VEBUEappTFGXm0dYGDA8DV11lzx3asAE4/nhg/36ZQ9TbK6Y5Je+oRqQoysxjzRrg4YdlnlBPDzAxAfz1r8Arr4hgmj9f2qlGVBCoIFIUZeZRWws8/jhw7BhQVwfU14swAoBnngH++Ed5rxpRQaCCSFGUmUltLXD++fbnU0+V1y1bgF/8Qt6rRlQQzEpBRESd+R6DoihZprcXePZZ+/Mrr8hrRQVw0knyXgVRQTCrghWIKADADyCQ77EoipJFenslOMHJ5KS8Pv88wCzv1TRXEBSdICKiBgDjzNzlsq0KwGoAIxCBM+JsZ94TUShHw1UUJR/090vE3K9+BQwMSKCC4Te/AQ4elPeqERUE0xJERHQugAq3bczcM52+4xwvAGAHgEaXbX4AQWaud6zrIKIRZh7J9FgURSlgTAaF+nrgjW8EfvITe1t5OfDzn8t71YgKgrR8RER0JRFNABgE0OWyZNQHQ0R+S4vxAxiN06wZQLSmEwIQzORYFEVJgXxXRB0fl7xy8+bZ637/ezuCTjWigiDdYIUggJcAbARQ7bKszsjoLJh5hJmbmbk9QbMGAENR6was9Yqi5IN8V0QdGwN8PmDuXHvdSy/Z702qHyWvpGua8wNoYOZHMzmYdCEiH1y0JWYeJyIQkV/Nc4qSB0xWgyuvBD78YeDBB3NbEXVsTExxS5aIdjRnjswtMpSV5WYcSkLS1Yi6EMc3lCcqABE8cbb7czgWRVGc1NaKL+buu2UOT66E0OSkCJ/ycuD002WdiZYzqCAqCNIVRM0AmonoRiI6l4hWRC+ZG2JS+JJpRERVRNRivQ9awQ/x2jYR0QARDezduzdT41SU2UdXl/hiqqqABx6IX4Qu0/6kgwdFGJWXA4sXy7q/+zt7e1kZsHdv7vxVSlzSFUR+AJUA2iABC8MuS8HBzEPM3MbMxMytbiHgjrbtzLyamVefcMIJuRymoswcensljBoAVq4Us1y8iqiZ8Cc5hdnYmLz29wN/+IO8/+Uv7balpbI+V/4qJS7TCVYAgLsAXBtnyTmWr0hRZhf5jkxLRH8/cP/98v6vf7V9Rm4VUc22Sy8F3vpWEUKp+JPa2sQHZISZEUSPPAIcOSLv3/Meu/3+/ZKNW4vi5Z10gxWqADQx85cyOZhpYAIRKgCE/UQOwaSBCsrMxWgS5qZtNImdO/M9MpnP86tfyfu//EVea2vj3/xra4FFi4Cf/hS49dbkhIQp+XDGGcC2bZJf7l3vAjzWc/bllwPf/ra8/973gE2bRDi96U3An/40vfNTMkK6GtEQ4s/nyTlWkMIIYn1FFZAsDCqIlJmL0SSuuEJuwKlqEtnGCCBndoN49PYCr70GLFyY2J/kxJR8uO02CYbo7haf1N/+Jtsvu0w0JUC0re5uEXIvv2yn/VHySrqC6BoAQSKqJaKFmRzQNOhC7PylKmu9osxsamuBpUuB73wnt5FpyWDS6UwliIwmN2cOcPRoYn+SE1Py4fBh4PbbYyepfuxj4g/avBn4j/8Atm4VoXXFFdI2GWGnZJV0BVE3JFihC8AYEU1ELcem2H86VMA9Sq4Vsal/mq31ijKz6e0V8xRR8ppErjAakXmNR38/8LWviT/n0CExscXzJ0VTWyvBEIYzzwRKSuT9oUPApz4FvOUtEkK+bZt8P6b9D3+Y+jkpGSVdH9EOADxlqwxh+Xq2QqL1fBBtrB5AJzPvAsKTV1uJKAig32obVLOcMuMxmsSSJcCrrwJf+UphmeeSNc21tET6bPbtS+xPctLba+eTmzNHCt/NnSvzhubMET9VyMoAtmqVCDczh+j661M7n0zQ1iYmRee59fbKuEyevFlEuoLoc8x8IKMjSYDlA5pSs2HmIcSm+VGUmU1/vzjf16+Xz2eeaWsSxSSIAGDU4Xretw847bSp9zElH5jl3K+5BvjEJyS/XDAoQujRRyWM3Ai22lrg3ntl/0OHgAULUj+v6ZDLAJMiEHrpmuYGiejGjI5EUZT0aGmRG43JJP3aa3LTyedNxhlS7jTNBV1yEDvb7ttnr7/vvuSOZUo+lJcDb3+7pPC5+27gH/9R3odCkSY+c7zjjpPPxk+Uy3B3E2By2WXyPpsabL7z/SVBuhpRD4CrAGzP4FgURUmXV1+13xdCJhDnE78JVgCAc85J3NbM/QHsKqpT0dIi2bS/9CXZx00AO0185nhNTfJ5927gxhtzH+5eWyta4u7dyYeqp3uc//f/JKLy4x8HvvzlwjHbWqQliJi52Up/8xSAZmb+fYbHpShKKjgF0Wuv5W8chtpaMRf+/d9H5nNb7ZKY32gH69fLBFPDiScmf7y9eyUUOxnhZY5nKrjecAPwzW/m/sZs5jYBEmCSrD8sHbxe0fyCwewKvTRJSxARkTHkVgEYIaLoJszMRVf9NeOYiXZOzjjDnuBnqKwsGFutUqT8z//Y7wtBIwIkKODYschoub/8BXBLmVVbK8EFf/6zvc5pppuK//5veU1Wi6qtFbPY179um8dyhckAcfvt9rqtW2Uc//mf2RnLD34gr1demX2hlwZFETVXtKxZA3zuc3YRrokJeSoxYaVer8xv+NSnbIdiLhyI0QLS5AJ7+GF7nQrHwmCqa/XjHwPnngucfbbdZmhI9svH9XM6xr/1rdjt8QIWurslU7Zh3jxbs3NzrDv/L07/WHe31B+a6n/U2ytzrgAJZPjgB3N3Y16zRrS/Sy8VTQyQkPLbbstOgElvL/B//o+8v+QS4CMfKayoSgBg5rQXAOsAPADgu9bruun0V6hLdXU1p0QwyNzUxNzTI8vChcweD7PE9UQul1/OPG8e89veJu22b5f9e3qkj2AwtWMnGo9Ztm9nnj+fec4cWTweZiL7c0kJ85YtkftkaixKapjfz/z5spSUMJeW2p/NsmqVXMclS2R7T0/+xrtkifzGFi2K/b3fd5/7PscdF9mOSM6pqUn6iz4f53GWLGHeuFH2u+UW9/Zu+3Z2yve5adPU+2Sanh7msjL7fJ96KnvHCgaZr71WjnPXXfbxc/B/BjDAyciSZBq57iiCZxKS6mfAep0A8HC6fRbqkrIgsm4eR48/nr/+wQ/y4fJydyFkhNMb3yivRgisXCnb5s9nfve7mdevjxQGbj8is975eft25ksvldeFC+WHX1Iir/EEIyBtzGtJiYxr/vzsCMl84nxgMDQ1MZ95prwaenqYzztPFnMdzENGLgR1Tw/z8cfLzdn52ykrkwcG86Dj9cpy2mn2fvm4Rj09zHPniiCJ/m1dc01ku2BQlo9/PPa/QSTCLPoamX2bmpgXL2a+4gp7v2QEivO/snIl84YN+fmuTjzRHvfevdk91q23ynE++cnsHieKZAURSdvUIKJrAIQA1DNzt2N9A4CdAG5i5hkTUbd69WoeGBhIaZ9f3Hsvln/0o5gHiZFnABGeNI8nMs9VeXlkxBAgdmQiSXfi8Ygp741vBH7zG3H6Ok0ZPh8wMACcfLKYNYikDPIFF4ip5vTTp1cWec4cseF/4APAV78qJsULLxR/xJ132ip+b6+YjVI17aXrT0vVzOhsf8YZwGc+I2Yd41v405/EhFpWBqxYIfnKXn5Zzt3rlW2mwieRXYJ6zhzgscfcTR3OY5q5LOZ7cjsn535r1gDveIdMVt2/P3J7RYX8Zj70IYmEMpx5pmS8zqX5xWmSe+21WD+Q+b3ffjvwyU9GzpuprQU6OuSzgUhu0Z/8JLBunWy76SZJ4/Pyy/LdXXWVRJz9+tfyezx8WNL4vOUtU//2zHg//3n5jb3wQm7m1pjjHjggwRJlZWKuv+oqST8EZGccH/+4nOtHPyqvOYKIBpnZJUIlimSkVfQC0YA+EWdbC4DfpNNvoS4pa0TMfMcdd/CTjie8wwC/7vHwpFPj8HojnwKdT7xu640WM2dOYg3LLGecEV/rcdsn3vEXLoz8PGeOmBLnzZOncq9XPq9fL58XLbK1p2Qw2lu0Ccp8V+az0RCNZrJ9u7w6zYxer226MvuYsQSDokGUlcn6RYtkzNHn67Zu/nzmyy5z//4WLbKfsKOfrKPPraxMzLHOMcT7rowJ6f3vd/+9AKIpuY01WjNw0/4yqdmasfb0MN97b/zf25w5zFdfHTu+u+/msKna+bs8/nhpt2WLfPb55Ls0mpKz71NPld/w9u3Jj3fTJvmddHbmxjxnTK3z5smY3/EOeS0rs7XsbIzjgx+U41x9dWb7nQJk0zQHMcm5+oMAXAlgIp1+C3VJRxD9/Atf4L8BfMRaXgf4CYDvA/iwx8P7zzpLfnwrV3LYpOA0SUwlQNJd3ExyZWVy8061L7f95s2z7fbJCCNzo16yRG5G0b4CIyCNWbG0NNK0eP75ttCKXkpLbfPi295mf9eXXx5fmMf7zhKZMs1N0Xne550nJr716+1zmzs3sk8jhNxuPMZ8ZG7Qzht5vGu1eTPzOefI+/e9L/Y7NubZ7dtlSeehIRE9PXL9410P53LrrbKPMWt+7GNy7RcujNzf67U/m2tQUyNjdvudbNmS/I3cmDwB5oqK7Ashp1l77ly5lkYgzZsn55UtYWjMl+95j/uYnGTQTJltQfQ0gPvjbNsJoD+dfgt1SddH9Ppxx/E6Il4L8Ji1fBHgtQC3EPGXzzqLj/h89lOguTnG+yO73aTdFiPUzGL+bEbAGc3BaBFlZbI4gxfcxpHoZuxcTNBFImFkntBNu09+MrFw8HjsG3my4zDCyPn5ggvc282ZE/sAEO84buvNvitX2sIPkO/1bW9z12acNx7nn98pnC++OPIYW7aI8Hjb22zhDDCffnrkk7Z5wnb2tX27vd1cX+e26WpHv/99ct+V+Y2sX29rN295i3xH7363rDPOdbOccELktTrzTPs7P+00DgvieD7UeHzgAxwhuLN5Yzbazn/9l2itzt+n+W0bAe3GdMa2bp30v3ZtZH/m+nd2yjoTyJQhYZhtQVRlaUX3AzgXwELr9QFIwMKMip6bTtRcKBRir9fLay0hdBPAsJabAP4XIv6b18s/WbaMRz76UTuSxgQJmD9y9Gu8xdzwTDvzZzXCaeVKW1CYp9FoR3xTk23CMlqE6S9eoEP0jdzcgEtLpS/n0+B550XehKJvOpnQCJPpI3q8zigmwI5Gc9u3rExu6uaG4rbU1cWui/7u3v/+yAgw8z0tWcJ8222R45w/P1J4b9ki26qr7WtgzuGmm2R8l19umyfjmfGMtuQ0Mcb7XZsbnzHzOYM3Tj45tn/z+4n33RvtDJDAA+dN8Mwzp76GixdLX5s3p65N9PTI/uZ7cDONZdpU1tNjC5/S0shrsmBB4mNNZ2yrV8sxVq2K7W/zZvv3mqxpM0myKoikfzTAjpSbsATTJIBr0u2zUJd0THNOnnvuOd6wYQN7PJ6wEHIKo7XW+xYifmr5ct53wQVyszbmE/OUPZWPKPqp+6yzOCzUFi2SG5fTZ5II87Q0b548pZonpXhRd8bnEX2jMU/s5hzMzXL+fDGrRfcTvZ95dQoNc5zofb1e2Wb+6M59orVEs+/KlZFa05w5zG94Q6RWY9qa9c4oQuMXixcZGX1MIlsDjhYGS5bId21u7E6NtKrKvnEYYWTMoMyyn/GJEYmGYY6zZUukJhRtapwzxxZC8SIyncc0vyMTqed8YIpe5s61BWX0Uloq358RBoDcFJnlWEbAOK9j9AOGczzxzJxumO+yu1v8Tu95j73vk0/KNb3hhsybynp67HO49VaxBJhzMf6wZITR9ddHXn/ndrf/9pveJMfw+yPXP/JI5Pdpvv8MkXVBFO4ACEAK5dUBWDTd/gpxma4gMiQSSM7F4/HwRRddxN88/3z++Re+EDkPyDjr3/1ueWI0y/r1tl/CfDah21OFfsfDaQowgskZBGC0IKcZav78qU2IZWWR/hKzlJfLzcloUsap72baMje+5csjb1CXX26bruLtY8Zs/AnO79KEbW/fLue3bJn9UOA0YZnFrIv2bbgJoZISGR8Rc3NzZPtFi2QcRst6z3si940WFPHMNOedJ/2b5dprIwUtkX0M5039yisjNbNLL4007Rmz31vfmvjaRn/XS5bYAsXtOjiv4Zo1dgCJucn29Mg60zZaEF12Waw2mczv2/n9nXcecyBg7/vUU3b/iUxl6dDUJOdw7rkybmMaNMvRo1OfQ329tP2Hf0heQ3rDG2SfxYsj13/ta/axly/PuODNmSCaDUumBJHhueee42uvvZY3bNjAJSUlSQmlDRs28LXXXsvPPfdcRsdixnPHHXdM3Xf0xFhzYzzvPNtnYZzzp57KYaGTzI3LPJU7b/jRc6C2bJGbtxEkxt9jJt+uX28LGOc+zjFccEFyEWvJ2OOjI+IWLYrVdsx3EH0eTU3iIHdGRBoB4dRegNSe9nt6bI05WlM5/vjIY7iZIo2Pxnxv69fLeS1bFntebuZP5yRWY1rbvNkOZDBajJvQ/u53IwMqnILRaGDmYcMpzNI1JZlrvGpVpGnu/POl7xNOiDQVTtd/ZM7H45GJt+az8zvYt2/qPsx1W7xYzn3xYts0GW/OlXkomzOHeXLSPpeaGvvYRkPNoDBKVhClNY/IkKhMOOewXlG2SWceUbL09fWhra0NTzzxBCad84ri4PF4cOGFF6KioiIjxx8dHcUPfvADTExMwOv14u1vfzve/OY34+qrr0ZNTU3yHZn8Wdu2Sd6sbduA5cuBwUF7/o0bRJKO//bbZU4PkaQ8OnZM5lGYOThXXWWncwFkHskf/wgEAnZbQOZg3HWX1LEx+wwPy/wgIpmTs2bN1HN4Ujnv6HlJR47I3K/Fi2Wez9y58t186lP2vKhHH5Xv6Y47ZB7T3/4m68vKZH/zW6ivl3YbN8qr81zdMLV5Dh6U24thxQrJ83b++dLG9P/tbwPve5+dq86kn/J4ZMwmLZXZdvSoPc7XX4+8tmYuT2mpjPdrX5Px/+hH8t1v2iTfUW8v8OSTcgyzr8cjNYEee0w+m5IN5jf13vfKZ68XeOghOY+eHqCqSvpN5xqauUyrVwNPPQXMny/f2ZEj8j3PmSNJQrdtszNzx6shlMxcrbY2mQd45ZXAF74gBfl6eyUxrJmbNjIic/4Sjffcc4GuLkkq+5GPyG/4Rz+S5K2XXy7pg267Dfjwh4E9e2Sf116T3+GRI5Jm6Uc/sn8nb3qTXP+//U3mIDr/f9Mk2/OIVgH4DWz/UMySTr+5WiCVXfOmEbmRrNkuV4vX6+VQKJTaSURHaDnNVW6aAlHk034mIrfyiVNjdJpLojMvRH9PS5a4BxKcfbZ7VF08nKaZxka7H683Uts05l2nc94tHY/bmIz2Zp6u462L9nstXBh5PPN7MGYmQEy68cLYnRgN6dZbp//03tMjQQJu51leHnv9nnpKftPvf3/isPvoY5hr99OfSt87d9rbjdkMYN6zJ35/5v1550nbX/1KvgvzH/P5ZLvxrwHyvX/mM/LeZHJ45RXp75JL5PPFF8v1q6vLaOg2c/IaUbo3cpPS506IfyhmSaffbC8Qf1aTyN/CEkQGY7a76KKL8i6UiIg3bNgQYbKb0owXHVnljMAzIeLGyW0mr2ZyLkux4CaMnFFnRgil4oA3372JQvN6bfOaiYpzfs9O09eSJeLIdrshO82HZm6W6dtMBZg/X8x3xqfnnKDp/E04g2CM+c3pd3SmVnIjG1Ft//IvtiCMPndj/jLH6Ouzt7n5j6Ya39NPy77f/769T2Wl3Wdv79T9GWH/xS/KNuO3+7d/s9s65yWaQBqz3403Shufzz7u0qXyO8kw2RZEkwD+KZ194/TXACAQZ1uVJTyMEHFtl+LxOJX2uRRETpy+pHwKJo/Hw29/+9sjxjBnzpywxpSUcHLzLc32hKqJ5vhEa4qp+CFMpKNTA3HTNhNpsE4BZDJnLFxoazdm0q7zGprjTxU0YLY7b7Lm9zCVUMn0PJ/oEOZov9cll0SO6brrZP2b3xx/rMaPYwIpnG0eekj2//WvI31U5ni33x57Lj09op2eckrkA8u118o2E9hy/fXyOToziJlnZfxEl15q+5MACZgoKZGHnwyTrCBKN9fc0wC+yMyPprxzbF8BAB0AGpm5K2qbH0CImesd6zoAtDLzyDSOycwcU0QpHtn0EaVCX18fHnroIbzqLIKWAZYuXYpVq1Zhz549eOGFF/DMM88gmd+Fx+PBJz7xCXz+85/H4cOHUVZWhu7u7gjfUl9fH3bv3o21a9em5nOaLRjfmvEtfe5zsv6WW+wSIanY65053wyJcpc5j08E/MM/SL0awymniH9t+3apL2T8dJnIhZbqWDON8bls3SrnPzEhPqiSkshaSLfeKgX75swB/vVfxQ/2trcBjY1yfaJrCB06JL4Ws+9tt9nb7rpLzu3AAckNuXGj5If86U9l+9y54pP6l3+JHOfmzXIdbroJ2LFD8kx+5zvAO98p+R8fekj2XbZMfJZvfjPwi19Enq/x6wHS3+OPi4/ohReAf/on8Rt1dmY0N2GufERXAFiYZh9+SOLUJgDDcNF0AAQBNEStCwDoSOeYjj44lfb50ojyhZmEiyS1Jef7O+64I9zPs88+y6Wlpez1evm4445LOeLve9/7Hn/uc5/LSqRgQZHtPHDJHt85T8iEwM+fb8/TyeR8mkLAbYKu8T+ZeTfHH29rjM5Jz+Xl8bXVL35R2pSW2pqH6f/GG+U77u6215eW2nOpNmyINcU5NVWnOW3XLmnjzD4O2NMfTjtNcumZ9SbK1JSeMfv9+c/i81q6NOO/N2TZNDcKMc/FC1Y4lmJ/8QTRMAB/1DqfU5BYgiyYYHHrl1MZ32wTRMzpB0+YUPMNGzbwaaedFrFt5cqVSYehd3Z2Miw/VTpCTEkD543ZeSNMdf5ZseI0YzpriH3oQ/LZ+JCMies733HvwzlH7n3vs0Ply8okIKC83M7z19PD/IlP2O1vu03WVVSICc20MyHlN91kt33wQTlmdAaP0lLxO519tj1up8n1zjsjw+mPHmX+yEdi5xhlgGQFUbqmuaB1c4kLM9+cQn/DAJrZYZojIh+AMQDlzDwe1Z4BVHKa5rliNc3lg2hz4NKlS7Fw4UJs374dE/FCspPA4/HgvPPOw+LFi+HxeODxeOAsOd/X14f/cZS/XrlyJc4444yIPpYuXZp6mLmixCO66us3vwncd5+YLEtLJcz5Zz+z22/cKCHUxpw4PCwVc0tLgWeekTaXXCLbmSVE/fBhWT9/vpj0TD+XXirh7u94h1SZPfNMCeWuqgLuvlvKYTz3nIRr33ef7HfffRKSf++9MpXBOa6dOyND4h9+GPj3f5cQ7YYGKbtxxRUStn7okJgc77pLzMOU9K1xSpI1zaVVKpyZW9PZL0UqrGONx9nuB5C2n0hJjpqaGtcb/YEDB/DFL34x7X4nJyfxwx/+MOn2L774Il50qafU3t6ecF6V8X/t27dP/VRKYoxfyvhILr4Yh77xDRw3NgZ+/XWQEUKmttKjjwLXXSc39Te+UebMzZ0LLFoEnHiizM357nfF58RsCyHA9gE55yF961siwJYsER8SAPz+9/Jqao8NDdl9/PWvIjBffjnyPP7rv+SYPp8IytpaWa66CnjPe6SeGSA+o1NPlfcLF4ov8vXXZV5fjklLEOUIX6Y7JKIqiI/JaHWdHBUgoSTH1Vdfja985Ss4cuRIwnYl1gTJo2YiZIaZnJzE97///aTaejwenH/++XjDG94QoX3FQzWu2U33pz6Ft46N4XUA5Wbl4sWiVRw6JDfuBx4QwfPDH8qNfXISeOUVubEbnL99r1eE1ec/L4Jq50578u7JJ4sAOnBAAiOOHZOik+vX2xN/9+yx+3r4YRFwixfLhFXDxReLprNlS2TgQW0tcM45djHHfftkX8Ae74EDKoiyDTMPARgC0DZVWyJqgvifsGzZsiyPrPioqanB7t278dBDDwFAOOrOGdFnbuQAwua90dFRPPvss0llkcg0k5OT6OvrS2mf9vZ2VFVV4YQTTsDcuXNjBJgz4jDZaEYVcEVAby/O274dnwLwaQDHAHgB0L59wPvfLzd6c/Pfv982Zxmt58AB0UrOOAP4+c9lndcrJrlPf1qi9B55RLSUNWtEM3JWaDaZFhYvlqwjN94on03UGyD93n23CLrKSjENLlwo42IG3vCG2CjEE08EfvlLeb9vH7B0qbxftMg+lxNPzMx3mAIFL4iIyJfAPJc1mLkdQDsgPqJcH78YiGe2i9fWkGwYeqKbfK4E2uTkJLLhHwy6eFhaAAAgAElEQVSFQli1ahV8Pl/YR2b8ZEQUfu/xeADAVYNTs2MW6e/HLy6/HJ+2Uvv865w52DZnDrxHj4qvZd48oLzcFh7MkULChIH/9rf2upISEULbtsnrr35lC4mtW21h46SzE/i//9d92yWXAGvXitA7elRCsr/2NQkLB4A//xn4t3+z0xMBInh275b3r70mZdWBSI0oDxSyIDL+nwoAYUFkBTE4tytFRioCLBFTCTRnHr1Cg5kx5LT3TxOTJ3CqHIQqvJKkpQVLP/IRPALgD2Vl+FxZGbyPPiqmsa1bRfM5dAiorhbfUEUFMDoq+y5dKjf0yUkRTlu2iA/plltkTpGZGxYK2cc7dixSsBm+9z3gs5+V94sWicZi3vt8MpeIWQTK//pfwK9/LfOBAJn79Y//aJvn2trErDg6KkEJ+/bJ53e9S/LdAbYgyuV8LhSwIGLmcSIaQayvqALAeLoRc8rMIRmBZibULl68OCXzWT5NiOkwMTGRtK8MEA3LbQKyYvNSQwOuvf9+fMbrRcmjj8oNvb9fNJpbbpHko7/9rSQafeIJe8eDB6VtV5cIoRUrJDhh1Srx67hNUF6zRvxA8+bZCXB9Punf+JjMRFlmOcZLLwHG1FxTI+a9d73LFkQTE2L6cx7j9tvl/c6dEuzw7W/LultvlfUHDkQmc80R08q+nbFBuIRvW+tDAAYtM5lZ1wBgEzM35mp8szl8ezaTjAkxVR9RoQm48847D/fcc0+MMNKMGMB3v/tdvPOd70RZWRkOHTpkb3CGeZusFIcOiamstjb1rBjRN35nNu4FC0Sz+utfpc8PfECESyAg0Xtf+Qpw9dWSWeHUU0UQmazoxx8vodvOgIXPflaEjsmU/s//DFx2mUTTHTwoQvW555LPKD4FWQ3fzgIVcI+Sa4Wk/2l3rGu2FkXJKpkyIUaTiVRNmTI7Pv/887jwwgvDIfBGsN5www14/fXX4fV6cd9996GpqWlaxylGDluBB687fT9AZJh3W5sIhquusm/ciTQfN/r7I2/8mzYBX/+6BEAcPChlIV56SYSeOc5b3yrReeYBecECeTW+xNZWEVbRZSpMgIUJqnj5ZZlP1NgoQu2JJ0RQZTDNTzLkTRBZvp6tkPlAPgBBIqqHhFTvAsLmuVYr1LrfahtUs5xSzGTSR5aK2TGe8EoUAn/s2DFce+21GB4exoEDByImNs/0yL/Dznk/8XATNGbeTrI4++jtlblHJvAAECEEAPfcY/d76qlionvveyUUfOFCEX4lJcAnPiFh5YGAHR5u9jvttMhjd3RIuPZZZ4kWVVkp+6Z6DtNkStMcEU1iiiwKLjAzF4q2NW3UNKfMFIw29uUvf3nac7viFWl0C4goRlPf17/+dbz//e8HIMI6mbln08KZhHXbNin+99RT9hyib3xDAiH6++X93r0S8HDZZcD998t+mzbJOreifWZddbUIOzMxt6xMTHVmjtTll4tQe/xx2W8aiW4zaZq7C6kLIkVRChCjjV199dUpVQZ2I5nJxF6vF2effTZ+9rOfYXJyEh6PJ/0qwDnGqRH99a9/xfHHH5/dAxoTnXnds0eCCQwPPiiZFd77XhEer75qByvs3CnmNhOcUFsbqw2ZfvfuBZ5+WoTQ5s0SRv6Tn0hGh+efl8rKRJL252c/k0i/LAcvFESwQqGjGpEyU3H6q6IDKbxeLwARONm4TxiNqry8PELbyKXZL5Gm9sADD+C6664DALz00ktYsWJF1scTxqkdmTlERBKMsG6d5J6bnBRT3F/+Itu3b48sHxGP5mbReD7+cTHDGaFnSlwAkuNuaEgiA19+Oe3ghWQ1IhVESaCCSJktGMEEIJwVI9oPlYvIP5MU94QTTggLxGimK7D6+vqwbt06vP766ygrK0NPT09EX/fccw/++Z//GQDQ39+P1aunLquTMZy1ms4+W7IobN4s/pxHH5Vghe5u27xWVyfh4lMRbbJzft67V0x70WzeLPOU0phTlJOoOSJ6LySAIAZmvns6fSuKknvcAincbvTxIv8yFc2XbFJcZ9LbVAXT7t27w+a3w4cPY/fu3RH7OvMomsjBnJkSzU2/t1dMcLfcIr6fnTvF/HbZZbJ9chL4u78T01pv79RaS3SEntOE19IC/O//bc9DAiRi79//XVIJZZNkakVELwBOB7APUpPI1CWadH5Op99CXWZjPSJFSRdTOj4UCoVL3V977bXhz9kse+/xePiiiy5KqubVc889xyUlJQyAvV5vTPvbbrstot+c18VyllKP/rx2rdQSOussu4bSVGXWp2L7dimV/sY3RtY32rIl7b6R5cJ4TwMYALDK+jwK4FxIGPbTAP4pnX4LdVFBpCiZ5bnnngsLKeeSSSHl8Xh4w4YNCYXHjTfeGC7oGM0tt9wS0Z/X642oQJx1nAUKDT09zOvXy627vj5SCLlVi00WUyhv+3bmz31OCvMBzH7/tAojJiuI0jXNBSCVT01O8hFLuxq35vzcCeBLafatKMoMJ9FcqmQm/Cbjp5qcnMTjjz+OJ554Iq75bsmSJeG20cd+yczfgaREmjt3LtauXZvqqaZPPJ9MT48dmOAMakh2Aq0b/f1Sx8j4jaxku/jjH21fVTbnFSUjraIXiAZ0hePzTlhaEIAroaY5RVGyjFOrSkWT8nq9HAqFmJl569atDIBXrFjBn/3sZ+P2sXbt2sIoVx9PS8pUGXen+Q9g3rx5WiY/ZFkj6gJwFYDHrM87AWwjolEAQTiyZSuKomSDaK3KaDMvvPBCQm1pYmICW7ZsASCVhgHgd7/7HT75yU+6ti8tLcU555xTGHOeMpHJIRHOYAaPB1i2LHY+UhZIVxC1wpHvjZl3EVEzgF0QIfRPGRiboihK0jgFU19fX8IJu5OTk7j++utRVVU1Zb8LFizAwYMHMz7egsQp6DweCVfIQbqftAQRM78E4OaodfVEdLq1TVEUJW/U1NTgsccei5mw+8wzzxh3Ao4ePYpf/epXMft6vd5w+Pkpp5wyuwSREzNHKQdkNB+cCiFFUQqJaPNde3s7rrvuurCgGR8fD1fC9Xq9+NCHPoRNmzah1tIAFi1ahPnz56sgyjJpCSIiWjhVG2bOT81ZRVGUODQ1NWHPnj0IhUJhzejUU0/Fli1bwml+mBkejweTk5MoLS3Nu0aUt4SxhS6IIH6gqXIDueflUBRFySNXX301HnzwwXCdoTPOOANbt24NbyciHH/88Thw4ADmzp2LBQsWYO/evXkZa19fH9auXYsjR47guOOOy21F3SIQRNfCXRCtBnANgNlXRUtRlKKgpqYGTz75JOrq6gAAfn9sljIjiPKtEe3evTtcruPIkSMxaYiySqELInaU7o5iBxH9FkBl+kNSFEXJLuvWrcOiRYuwf/9+LFwY62mYP38+AORdEK1duxYejwcTExO5n1CbQ0HkyUKfXQAastCvoihKxjjllFMAwFUQmdpD+RZENTU1uPjiiwEADz/88Iz1EWVDEAUQJyN3viGiKiIKEFEDEYWIqCDHqShK9jn55JMByDyhaIwgMj6iI0eORGTjziXz5s0DALzpTW/K7YEL3TQ3RflwAjCY9oiySzeA01ly4lUA6ABQnecxKYqSBxJpRNGmOQA4ePAgFi9enLsBWhw6dAgAsH///tweuNAFERKXDx+GpPxJGiJqADDOzDGVnYioChIEMQLRtEbc2iXJ6cxs0g+NptmHoigzgFRMc0D+BNHf/vY3ACqIYmDm1kwNgIgCAHYAaHTZ5gcQZOZ6x7oOIhph5pFUj+UQQoCkKMrYeSiKUlyY8O2XX345Zlu0aQ5A3vxEs0EQZcNHlBRE5CeiEETLiaedNAMIRa0LQRKrTue4LQA6p6FZKYpSxPT19eHee+8FANx8883o6+uL2B7PNJcPZoMgmlIjSiaLghtTZVawNJpm6xjxNJMGxAqiAQCdjvE1IXG4eITAsY7bRkRNRNTp1LYURZkd7N69O5zm59ixYzHzc+KZ5vKB+oiEZLIouDGtzApE5IOLtmQFGoCI/Mw8kmBOU3R/fgANzNxmrdoJIGT6mc5YFUUpLtauXYu5c+fiyJEjrvNz3ATRXXfdhccffzyisF4uyJtGRFRQgsgti4IPYh5rRWTtoTcCuAlAmmUCI6gAYvw6TvywKsMmiR/A4qjP4yqEFGX2UVNTg+7u7rg53Jw+ot/+9rcAgO7ubnR3d+OrX/0qent7cyaM8mqa43R0kNSZUhC5aRxE9ACAJmaOKQdORK8h8oafLr4M9BGGmbuIyGeZ8gCgHkBdvPZWuyYAWLZsWSaHoihKAZCoXLnTR/TLX/4yYlsuU+1MTk6GgypMEb+cUQTBCpsQP8BgCAWaWYGZdzFzu7U0MvNQgrbtzLyamVefcMIJuRymoih5xmmaq62txdy5c8PbcplqxwghYGb7iKYTNRcvsWkjLLNaJrB8RYqiKDnDaZqrqanB7t27sXTpUpx00kl5McsBKojcuBPA3xPRU0R0BRGdS0TriOgRiIBKKoBgCozvJkKoOQST+nYURckKxjTX3d2Nvr4+1NTUoK6uLiyYcoUKogRYkWc3AzgPwDchKX26IH6XVmbemmD3ZI8xDhE20RpRBTTIQFGULPLSS1Js+vHHH0ddXR36+vqwfPlyvPzyyzh27FjOxmFCtz0eD373u9/FzHfKKoUuiAARRsxcAcnVthFANTNXMPNdGRudCLfVUeuqrPWKoihZwUTKTU5OhoMTli9fjomJCbzyyis5G4fRiCYnJzE6OhoWijmhGASRgZn3MPM3mXnPNLqpgHuUXCtiU/9oah5FUbLK5ZdfjrKyMni93nBwwvLlywEAH/vYx3DFFVdgy5YtWRcKTtMcYEfs5YQCm9CaFSxfz1bIfB4fgCAR1UMyIewCwpNXW4koCKDfahtUs5yiKNmkpqYGPT09EfOMHn74YQDAt771rXC79vZ2XHDBBTjrrLPwgQ98IMZ/1NfXF3euUjIYQeT1ejExMYGSkpLcFccrJEFERE8DWMfMcxzrEpWBAAB2to/TYBxJaDZWiHXcMGtFUZRsED3PyJjrnExOTuLZZ5/Fs88+ix07duDCCy9ERYXEV42OjuIHP/gBmBmlpaXo7u5OWRgZH9HmzZvx4IMP4v7775+1pcK7EJk9AUhcBkJRFGXGUVdXh9tvvz1ugbzJyUl8//vfd9126NAh3HDDDbjnnntSEiRGI1q3bh0efPDB3JahKCRB5MjN5lynPhpFUWYVZj7RQw89hBdeeAHPPvssJlO4UT///PO48MIL8fa3vx2LFy8GM+PIkSMoKSnBySefHJHDzpj0TKJVU53VrWRF1igkQeQGEa0A4GPmHzvWfQLAGgBPM/OXMzI6RVGUAsJpruvr60tZKE1OTuKZZ55x3dbe3o4LL7wQAMImPa9XckcvX74cXq8Xf/rTnzJ0JklQ6IIIUpphBMAWIOxHqgOwB0A7EVUy879mZoiKoiiFh5tQevXVV2PaPfnkkzh69OiU/bmZ9thKOrpgwQKcdNJJKoiiqIeVT46IVgEIAGhh5rutonPXAFBBpCjKrCBRAtW+vj60tbXhiSeeSMmU5+S4447DqaeemltBVGBlIOJhQqgDkMAFk9ZnGBJmrSiKMuupqanBY489FldrGh0djWvamzt3Lo4dO4bnn38eZWVl+MlPfhJOOZR1ikAjGgGwkYhGIJm4hxwVWVOtE6QoijLjmUprcgqpE088EQ8++CAOHz4MQAr5TUxMYGJiAhdddBHuu+8+NDXFyzudIYpAEN0MqXDaCoAgWpGhGUDHNMelKIoya4gWUn19fWhvt3NHHz16NOwvOnbsGLZs2YLh4WH4fD4sXrwY+/btS3vSbFw8HiBHefXSEkTMvIuIKmHlfWNmZ1rYIERIKYqiKGkQncbH45FsbBMTEwAksKGtzZ5ZQ0QoKytLetJsUhkfikAjAsQvtBrANUS0kZkPENEiiHB6JCOjUxRFmYWsXbsWZWVlOHz4MDweD+677z4AwJYtW1x9ScyMQ4cO4aGHHppSED3++OO48sorMTk5CY/Hg1WrVuHUU0/FSSedFDGXqeAFERHVAXgaEq69CpK09AAz77c0pTsBXJexUSqKoswiampq0N3d7aq1xBNGgMxFeuGFF8JphgxLly7FW97yFvz85z/HI488Et5/cnISg4ODGBwcBAB89atftQv/FboggpjfdjHzJivvnJN2yDwjFUSKoihp4hbcYAIUrr/+ete5SYnSDCWDye5dLIKoCsAdcbaNwb2kg6IoijJNmpqacPbZZ2P37t1YvHgx9uzZgx07doT9R9PBlLwAUBSCaAjiH3rUZVsztHCdoihK1ojWllatWoXrrrsuJWFUUlKCd7/73QCkJHpFRQX+4z/+o3h8RBAf0E4iIkjQwiIiOhdisgtAqrYqiqIoOcBoSfHSDC1duhSrVq3Cnj178Oqrr2Lp0qURgQmXXXYZ/vjHP0aaAgtdEFnh2xsh/iCCaEgEKRex0ZkMVVEURck+iSbMTsXy5cvx7LPPRq4s9FLhRLSQmXcxcwVE+9kIoJqZK5j5mxkdYQYhoiARMRGNEVEnEWkqIkVRZj3Lly/H+Pg4Dhw4YK8sdEEEYJCIbgQAZt7DzN9k5j0ZHFe2GGZmYuZyZq7XkuOKoigiiADg1ltvRV9fn6wsAkHUA+CqTA2CiBqIKBBnWxURNRFRwLxm6riKoigKwprQvffei7q6OhFGhS6ImLkZABHRU0S0fDoDsATLjjjb/ACCzNzOzF3M3A6geRomNZ8RepaZTsPMFUWZ9YyMiHFocnIyPJeo4AUREY1Csmz/PYARIpqIWqbMlEdEfiIKWf2MxmnWDJkc6yQEic5Lh3bLt9UFSUPUnWY/iqIoM4b169ejrKwMXq/XnktUBPWIdkDCttPG8s80AwARtcZp1oBYQTQAoNN8IKImAJUJDtVpCR4w87jj+EOW2c/nXK8oijLbuOCCC9DT0xOZUuiBBwCe1m0+adIN344nODKGZTaL0ZaYeZyIQER+Zh6xzHXJ9FcFYAczR8xxUiGkKIriEv5d6Ka5HFEBJBQUqfqJRuDQrizf1K70hqYoijLDKfQJrTkio4EEliY1YpnyADHnXROvvdWuCQCWLVuWyaEoiqIUPiqIsoPxFSXZth2SOQKrV6/OjaFUURSlUFDTnI2GWCuKouQBFUQAxKcDWL4ig0MwaVYERVGUbKGCKBykMIJYX1EFgHFNz6MoipJFikEQEdEKItpGRN8looXWukVE9ID5nAG6IHWPnFRB6x0piqJkl0IXRERUB2AYQD2k/pAJtd4PiUa7M8UuK+AeJdcKoDFqXbO1XlEURckWhS6IICl2djHzakgdIiftADZN1QER+ax8bx0QIRQkohARNZg2lnmu1WrXQEQtkNxzapZTFEXJJkUQvl0F4I4428aQxBwgI2SSaDcEKbynKIqi5Ioi0IiGEOu7MTRDfTiKoijFTRFoRHcC2ElEBEl+uoiIzoWY7AKQqq2KoihKsVLogoiZdxHRRog/iCAaEgEYB7CRmX+cuSEqiqIoOacIykCAmXcB2EVEqyAJSEeKpFy4oiiKMhWFrhER0UJmPgAAlvBRAaQoijKT8HhyVo8o3WCFcSJ6nog+nMHJq4qiKEqhUOgaEYCbAWyEVGptJ6IuAF9k5scyNrIi5fDhwxgdHcXBgwcxMTGR7+EoRYjX68WCBQtQUVGB0tLSfA9Hma14LD2FWfxFWSTdYIU2AG1EdDok88FGAN8kIoYUm3t4Ngqlw4cP4w9/+APKy8uxYsUKlJSUgLJ8AZWZBTPj6NGjOHDgAP7whz9g2bJlKoyU/GAE0eQk4PVm91DT2ZmZX2LmNivDQjlEIFVillY+HR0dRXl5OZYsWYK5c+eqEFJShogwd+5cLFmyBOXl5RgdHc33kJTZilMQZftQmeiEiN4LYKe1VGGWBi8cPHgQCxeqy0zJDAsXLsTBgwfzPQxltlIMgoiI3mtl3p6AaECVEN9RpaUhzTomJiZQUlKS72EoM4SSkhL1Myr5I4eCKN3w7UlIRoX9AO4CEGLmlzI5sGJFzXFKptDfkpJXCl0QQaLlvqgTWBVFUWYohS6ImLk50wNRFEVRCohC8hER0dNEdCxq3SQRTSRYjsXrT1EURSkCCkkQQUo6PBq17q4plrszOEalSBgaGgIRxSzV1dVoa2uLad/W1ob6+vq4fTU2NqK8vBxEhPLyctTX16Ory73CSHT7yspKNDY2YmgospTV+Pg4GhsbUVlZGW7X3NyMkZHM1lpsbW1VH49S3BSSac6avBq9rihLdRORW9G+VrdzVNInGAyioUEK7Y6MjKCzsxPbtm1DKBTC4OAgfD65BP39/a6Cpb29Hc3NzfD7/di6dSv8fj9GRkbwyCOPIBQKIRAIRLRva2tDa2trRPv+/n7s2rULdXV1GBsbC4+lsrISfr8fzc3N8Pl8GB4eDgvJUCiUza9FUYqLHAoiMPOsWCACKBC1riWZfaurqzkZXnjhhaTazVQGBwcZAIdCoZhtw8PDDICbmpoS9tHZ2ckAOBAIuG4fGxtLqX1nZ2f4fSAQYL/f79rn4OBgzPqWlhYeHh5OON54tLS0sPy9psds/00peeT++5kB5ldfTbsLAAOcxD02rXlERLSPiFbE2fY0EcUrIx6vvwYiCsTZVkVETUQUMK+pj1hg5vDjNxE1YJZmgMgHfr8fVVVVGBgYSNjOaEKdnZ2u24025Wzv8/nitndqTwMDAzHalOmzqqoqYt34+Dja2tpiTHuKMmswpuUC8RG5UZ5gWwiSfy4pLMGyI842P4AgM7czcxcztwNottanBDOPO/r1Aahg5sw6BpS4jIyMYGhoyFUQGHbt2oWRkRG0tiZn+TXtg8FgUu0rKiqmFISKolg4k55m+1DJNiSiFVY2hfdaq+rMZ+cCoAlSKG+q/vxEFLLaxkuo1QwRbE5CkJLk02ErJB1R0dLX14dt27ahr68v30NJyMjICNrb21FdXY1AIBAhNKId+v39/QCAjRs3JtV3qu2DwSCGhoYSBj0ACAc9mPcmWMJ5Ts7AiPr6+rjBDuZ45eXlqKysTHhcRSkoCilYwcG1AFqs94w4WozFlM5/SxtpBgAiivcI3IBYQTQAIGyHIaImSHqheHQ6TXIWAc5xwMUNN9yAH/84MxXU9+/fj5/+9KeYnJyEx+PBOeecg0WLFk2733PPPRf33HPPtPtpbm5Gc3PkVLNQKISmpqaE+5mbebT5LVPtGxoaEAqF0NzcHBYIDQ0NCAaD8PvtZ6cdO3agubkZ9fX1rsERJqhhx44d8Pl8aG1tRXV1dTgowsk111yDrVu3orGxEcFgEPX19RgeHo44nqIUJIUoiJj5ZiLaBoAgGkw9gOjHwFFm3p+JgVnmsxhtiZnHrbBgPzOPWOa6VPoNRPdZbOzfvx+T1o9jcnIS+/fvz4ggyhQtLS3YtGkTAMlIPjQ0hNbWVnR0dKCjoyOu4BgfH3ddH49U2wNAU1MTmpqa0NXVhY6ODuzcuROVlZXo7OwMCxyfzxcWFBUVFTFCI9oUuGPHDlRXV2PXrl3haEHD4OBg+H0gEEBlZSWCwaBG6CmFTyEKIgAwQoaIdgHoZ6tceJaosI4Z727jR6wgTIYqAKnfwaZJJjQNQ19fH+rq6nDkyBHMnTsXX//611FTU5Ox/qdLZWVlhPM/EAigqakJ1dXVqKuri7g5O6mqqkJXVxfGx8eT0nLita+uro4IMhgcHIwJRggEAmFTYXV1NRobG101mmQwfU81F8nv9yMQCKifSikOCmxCawzMvDHLQgiIne+TSaYUYFaE3gARDezduzeLQ0mdmpoadHd34/bbb0d3d3dBCaF4+Hy+sI8mXiTamjVrAAA7dybnvovXvqOjA4ODg0kFMRjT2vj4eEoRcrt27UJjYyOqq6sj/EdTYeZEKUrBU6gaUTRWcIKrsZuZCzK7Aic5edUy+bUDwOrVq7MfNpIiNTU1RSGAnBitJV6xt4aGBvj9fgSDwSn9SYnaG1Nasjf8VP01JjihtbUVq1dLxZPq6uqk9h0ZGVH/kFIcFLpGRESnE9E+yDycNkgUW5tjmW5Um/NY2dSMlBxiwrIThXCHQqFwVJob0VqLaR8vVZCT5uZmV79SMBiMmUtUUVHh2sfIyAi6urrCwq+qqippwTI+Po6urq6E568oBUMRaEQhAC9Bos/2ENEogHUAfge7Uut0MY+zFXD4dByCSe0bBcr4+HhYGxkfH8fAwACCwWA43U8iAoEAOjo6wmav5uZmVFRUhPcdGBiI8OUYP09raysqKyvR2tqKioqKmPRB4+Pj2LlzJ9rb29HQ0IA1a9Zg37594blIHR0dEeMw2psJKgiFQujs7ITf7w+b80y7RCZAEy1nNCifz4etW7em8G0qSp4o9BQ/ACYBrHN8HgBwrvW+DhLIkEp/w4hKv+NYXxW1zg9gLJ1xp7toip/kMCl+ohe/388NDQ0x6XISpcEZHh7mhoYG9vl8DIB9Pp9rH85jO9ubYzpT/IyNjXEwGOSqqir2+Xzs8/k4EAi4pvdxjs/v93NLS0t4fUdHR/g4VVVVHAqFGAAHg8GYfTs6Otjv94fTEKWaMmi2/6aUPNLRISl+fvaztLtAkil+iNOYNWtpQB9m5seszzsBPM3MXyKiKwHsZGZvCv0NA2jmqPk+1oTXQXaEaFupeTYxc9LZG6bL6tWrOZlIpxdffBErV67MwYiU2YL+ppS88eijwJVXAj/5CXDOOWl1QUSDzLx6qnbppvjpAnCV4/NOAK1W8EIQqYdHV8A9Sq4VsemCmq31iqIoSrYo9GAFiCB4yXxg5l0Q/9AuiFC5ZqoOiMhHREEi6oAIoSARhSyNx/Q7DhFwQSsxagsk95z6hxRFUbJJoQcrMPNLAG6OWldPRKdb25LpYxxJaDbMPARAU1KxKsUAABDmSURBVCAriqLkkiLQiFxJVggpiqIoBU6hC6JM1yNSFEVRCozZVI9IURRFKUByWI8oaR+RpQE5M0fWEZFblsik6hEpiqIoBUyBBitktB6RoiiKUsAUoiDiHNcjUhRFUfJIIQoiIOf1iBRFUZR8UehRc8y8EUAFEW0jou8S0UIAIKJFRPSA+awoiqIUKYUuiIioDpKQtB5AAHY11f0AKgHcmakBKoqiKHmg0AURJJ/cLiuZHUVtawewaVqjUhRFUfJLEQiiKgCPxNk2huyW+VYURVGyTREIoiEA8VJ7N0OycyuzjKGhIRBRzFJdXY22ttiI/ra2triVVYeGhtDY2Ijy8nIQEcrLy1FfXx9R7C5R+8rKSjQ2NsZUdFUUJUmKQBDdCcCEczOARUR0LhF9F8CV0DINs5pgMIjh4WEMDw+js7MTgUAA27ZtQ2VlZUSp7ugqqob29nZUV1djaGgIW7duRUdHB7Zu3YrR0dFwxVQnbW1tMe0bGhowNDSEurq6rJ6rosxYCr1Cq1VMrwEyn2gSwIT1OgrgynT7LNRFK7Qmh6nQGgqFYrYNDw8zAG5qakrYR2dnZ7iaqRtjY2MptXdWaC1GZvtvSskjQ0NSofXxx9PuAklWaE2rDIQlwHYB2EVEqyApfUaYec90hKIyc/H7/aiqqsJUlW6bm5vh9/vR2dnput3n88W09/l8cdsHAoH0Bqwos50iMM2FYeY9zPxNpxAionOn268Sh7Y2oLc3cl1vr6wvYEZGRjA0NJRQMOzatQsjIyNobU3OsmvaB4PBTA1TURRDoQkiIlpIRI9Y5R9+Q0QfTtBuJ4DBjI4yQxBRFRG1WNVeg0RUfNF9a9YAGzfawqi3Vz6vWZPfccVhZGQk7PMJBAIRQqO1tRVEdvR/f38/AGDjxo1J9Z1qe0VRUiCHZSCSNc0NQSatdkFCt9uJyMfM200DK3ChBcB+RFVvLQQsodPBzJXW5yHIfKjmnAzghhuAH/84M32dfDJwySXASScB//3fwMqVwGc+I8t0OPdc4J57pj285uZmNDdHfq2hUAhNTU0J9xsZkdSF0ea3TLVXFCUFCinXHBHdZL09ne1cc4MA2ojITF4NQia23szMd6U6CCJqADDOzDEhVERUBQkVH4Hti0onPDwAR5JWZh4hoibkShBlkvJyEUJ/+AOwbJl8LiBaWlqwaZPMaR4dHcXQ0BBaW1vR0dGBjo6OuILDGVGXDKm2VxQlBQqsHtEmAHdyZFbtVgBPAxiHCKA2ANs4jczbRBSAlJSIKaZHRH4AQWaud6zrIKIRZo7O/D0V47BSEUUfI42+UicDmkYYY4679VbggQeAT38aqK3NXP/TpLKyElVVdumqQCCApqYmVFdXo66uDoOD7pbbqqoqdHV1YXx8PCktJ157E8ptGBwcjBiPoihJUGA+oioA0aFO/dZrO4ByZr45VSFERH4iCkG0nNE4zZohFV+dhCAaWEoYLcr4hSwBCBRbFggjhHbuBG67TV6dPqMCxefzIRgMYmhoKO4k0zWWn2vnzp1J9RmvfUdHBwYHBzWIQVGmQyGZ5iwibCDMvN9yMofS0YKsPkZgmcWIKF6YVANiBdEAgHCsrmVeq0xwqE4jhJi52gpUGIFtpsu+NpRJ+vtF+BgNqLZWPvf3F5RW5IbRWkZH3Z87Ghoa4Pf7EQwGp/QnJWrv90uBYONDUhQlDQpQEOUcS3OJ0ZaYedxKHeNn5hFmbk+lX2v+kzH7jTBzcTkaWlpi19XWFrwQAhAOy04Uwh0KhVBfX4/GxkZ0dHTEbB8aGoows5n29fX1cecSKYqSBgUoiK4komjNhwFsJCLXnHPM/KVpjcwuLRFPUPiRojZDRGPMbDz7zdBURFlhfHw8rI2Mj49jYGAAwWAQIyMjUwqLQCCAjo4ONDY2orq6Gs3NzaioqAjvOzAwgLGxsYj2wWAQra2tqKysRGtrKyoqKuKmD1IUJUkKKcUPJHVPqstEMmkdHMcYBhCIWlclw3Ntz9HtkzxOE8Tc1wSgIYm2AwAGli1bllQ6i9mejsWk+Ile/H4/NzQ08PDwcET7lpYWtq5xDMPDw9zQ0MA+n48BsM/nc+3DeWxne3NMTfGjKGnyxz9Kip8dO9LuAkmm+CGeIjTPSuGTjoBLOt0PEQ0DaGZHWLYVtj3IzNH1jkBEDKCe0wvjTpnVq1fzVKlpAODFF1/EypUrczAiZbagvyklb7zyCnDKKUAoBCThs3WDiAZZ6tYlZErTXCoCJRtYE2eLy4+jKIpS7BRY+Ha+MP6fiLk/jrQ8GhKlKIqSLVQQhYMURhA7z6cCkoVBBZGiKEq2UEEUpguxlWCroBVgFUVRssssFEQVcM9w0IrY1D8adq0oipJtCnAeUcaxfD1bIfOBfACCRFQPyYSwCwhPXm0loiAkrZDJPadmOUVRlGwyGwSR5QOaUrNh5iFIGQpFURQlV+SwHlGhmOZmDFPNy1KUZNHfkpJXclgGQgVRBvF6vTh69Gi+h6HMEI4ePQqv15vvYSizkbY24Jln5L3RiHp7ZX0WUEGUQRYsWIADBw7kexjKDOHAgQNYsGBBvoehzEbWrAE2b5b3k5N2+Rmr9EqmUUGUQSoqKjA2NobXXnsNR44cUdOKkjLMjCNHjuC1117D2NgYKipiajkqSvaprQW+8Q15//TTdg20LGX5L9gyEMVIaWkpli1bhtHRUfzud7/DxMREvoekFCFerxcLFizAsmXLUFpamu/hKLOVdeuAlSuB3bulGnQWS82oIMowpaWlOOmkk3DSSSfleyiKoijp8/3vA3v3ihB64IGs1j1T05yiKIoSifEJ7dwJ3HabvG7cKOuzgAoiRVEUJZL+/kifUG2tfO7vz8rhpqxHpCRfj0hRFEWxSbYekWpEiqIoSl5RQaQoiqLkFRVEiqIoSl5RQaQoiqLkFRVEiqIoSl7RqLkkIKK9AH6f5u5LALyWweHkEz2XwkTPpTDRcwGWM/MJUzVSQZRliGggmfDFYkDPpTDRcylM9FySR01ziqIoSl5RQaQoiqLkFRVE2ac93wPIIHouhYmeS2Gi55Ik6iNSFEVR8opqRIqiKEpeUUGkKIqi5BUtjJcliKgKwGoAIwD8AEaYuSu/o0oMEfkANAHYBWAUQAWAZgCdzrEX6rkRUQOAcbexJDPmQjqveOdSbNfIOg8/gErrNcTMu6LaFMW1mepciunaEFEAQD2AfZDzGWTm9qg2ubsuzKxLhhfrgnRGresA4M/32JIYNzuWMQANxXBuAALWeAPpXI9COq8kzqUorhGABgBVjs8+AMMAmort2qRwLgV/bazfVyBq3SCAlnxdFzXNZYdmAKGodSEAwTyMJVXqAZQDqGTmco56ekWBnRsR+YkoBPlTjMZplsyY835eSZ4LUDzXyM/MQ+YDM49bY3COrSiuDZI7F6A4rk2zy7quqPW5vS65ksKzaYE8Kfmj1vnk687/+BKM2w+Xp/BiOTdrbG5axJRjLrTzSnAuRXGNrOMNAvC5jJ/N+Irh2qRwLsVybToABKPWBSHmuaTHmclzUY0ow1h24pgnWpYnKBCRPx/jygTFeG7JjLkYzysehXIu1vH81uJKsVybZM4lGQrhXKzjNTJza9TqBgCPJDvOTJ+LBitkngrAviAu+CGOvULFbzllATmXUbbNC8V4bsmMGUm0KaTzKoprxMzlLqsDkCCMEXOzKoZrM9W5OMdTDNfGCRE1ARhi5jZrVc7/MyqIMo8v3wOYBqMAwJGRQB1EZNYV47klM+ZiOq9iv0bNALZZ74v92jjPBSiya2MJzHpAtCTHppxfFzXNKWGYeZyjQjhRPEEWs4JivkbWk/eo48m7aHE7l2K7Nsy8i5mbAbQS0aAVip0XVBBlCcuGOhMYgZgbwudTjOeWzJiL8bwsCv4aWWa4Zmaud9lWVNcm0bm4UPDXxjKvhQB0O9fn8rqoIMo8xi5a4VzpuGCF5GuIgIhaXFYbZ6TT5ltM55bMmIvmvIr4GgUB1EWtK9Zr43YuxXxtAAnf9lkTXXN+XVQQZRjr6WIEsTbUCsQ6NgsG6ykv6BLtYn5oI8V4bsmMuVjOq1ivkTU3qjXasV2M1ybeuRTLtbEi3sYSmOF8+bguKoiyQxck7YWTKmt9QWL9cJpdfkABSESN+eMV3bkhuTEX/HkV4zWyfClB55iJKOC4YRfNtUl0LkV0bXyI1GgM5nqYSbu5vS7TnRyli+uEMR9iU190ovBT/DQ4xwh7Il90apOCPDe4pFRJdsyFdl4JzqVorpE11ibr5mSWACRHW1FdmyTPpSiuDYAWxE7O7YRjkmuur4vWI8oSluq7CUA/5GljiAsgMehUOBI7Lob80CKeAK02BXNulk16qzWOBsiTXhfkD+IMo51yzPk+rxTOpeCvkXUuY3E2jzBzpaNtQV+bFM+l4K+NNYYmSLLTqZKe5uS6qCBSFEVR8or6iBRFUZS8ooJIURRFySsqiBRFUZS8ooJIURRFySsqiBRFUZS8ooJIURRFySsqiBRFUZS8ooJIURRFySsqiJSixcrzxdYs8YKFiFqIqDPf48gnRBQkIp09r7iigkhRss8aSF4yRVFcUEGkzDisp+/odPx5Oy4zNzIz5Xo8ilIsqCBSZhRWgsoWSHbkGX9cRZkJqCBSFEVR8ooKImXGQEQdsNP1d1iBDGNRbRqIaNDaNmiVRjbbgqY9EYWIaJiIfFZVyw6rsiUTUafTBDfVcd0c9URUZfVjxtEUtT1oKmlGtUta40q2jzjja7La++L0F/EdWt+ROcawVQ7BbUxmLGNWu0DU9pSvT7Lfh1K4qCBSZhLXAKi33jdD6qxUm43WzbEDUuOnEVLvJ0KoWO0GIZUnQyyVNZuj+q+AFDxL6rjRWOMYtI5fD+ARACGrDLUTnzXekNWv3/qcCpnow62/bVH9maXRahfvGDussbRan8Pf/zSuj1Ls5KoqoC66ZHqBRKIxgCbHOr+1zq266RgcVSitdYNmHYCgtW/nFMetij7GFMcNyl8tYhyhOOcSiBpLQ3Q/iKqumWCcSfURPT5rXVO8dmaMUetaHOsarHVVUxzDfGehTF4fXYpvUY1ImRVY5igfgBbL7MOWOcqUfXbSHNOBA2Yest6mHJlnmZp8EK3A2WcXgHHYGoXBWe1yONXjZbAPJwMu/TmPYSqSViTqhKVyaReA1Zm8PkrxMSffA1CUHGGERjXkhu9k1PmBo0o7A2Gz0Sarn+mEhpt9Y45hrYvomzNgespEH0n0l+4xRiBmtmldH6W4UUGkzBbCN69Ub2RWVgQ/xDRktIHB+HskNQ4/gKGobX5EahuzAT/kO0n7+ijFj5rmlJnGqNtKy5w2DhEmESSKvLIc5QEArczcbvXjdqN0Pa4LRtBEmJcsjcsEAhQCcYMtMoX1vQcAdKV7fZSZgWpEyoyCmceJCACazSszm4i2ayDh1SaKzAcRCCOI43dg5hEiGgcQtPobhx3xlexxo9s1Q6LkKiARc0bb2mX5inLNMCA58SA+mwAkWCHjWN99K+xzHodE4AFpXB9lZqAakTITaYPcTINwmL+YeRfkSd8HoBMSSjwCF8ESxTUQx3uH1Wc8rcX1uNEwczskbNuEPjdDNK7oQIVcsRMyXnNui5G9G/8jkO++E6JFVhuf0zSuj1LkELMmxFUURVHyh2pEiqIoSl5RQaQoiqLkFRVEiqIoSl5RQaQoiqLkFRVEiqIoSl5RQaQoiqLkFRVEiqIoSl5RQaQoiqLkFRVEiqIoSl5RQaQoiqLklf8P2Iu4UcatykgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "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",
    "n = 300\n",
    "\n",
    "ex = np.ones(n);\n",
    "A = -sp.sparse.spdiags(np.vstack((ex,  -(2 + 1./n)*ex, (1 + 1./n) * ex)), [-1, 0, 1], n, n, 'csr') \n",
    "rhs = np.random.randn(n)\n",
    "\n",
    "# ee = sp.sparse.eye(n)\n",
    "# A = -sp.sparse.spdiags(np.vstack((ex,  -(2 + 1./n)*ex, (1 + 1./n) * ex)), [-1, 0, 1], n, n, 'csr')\n",
    "# A = sp.sparse.kron(A, ee) + sp.sparse.kron(ee, A)\n",
    "# rhs = np.ones(n * n)\n",
    "\n",
    "print(\"Dimension of the linear system = {}\".format(A.shape[0]))\n",
    "\n",
    "res_all_bicg = []\n",
    "res_all_bicgstab = []\n",
    "def bicg_cl(x):\n",
    "    res_all_bicg.append(np.linalg.norm(A.dot(x) - rhs))\n",
    "\n",
    "def bicgstab_cl(x):\n",
    "    res_all_bicgstab.append(np.linalg.norm(A.dot(x) - rhs))\n",
    "\n",
    "sol2 = scipy.sparse.linalg.bicg(A, rhs, x0=np.zeros(A.shape[0]), callback=bicg_cl)\n",
    "sol2 = scipy.sparse.linalg.bicgstab(A, rhs, x0=np.zeros(A.shape[0]), callback=bicgstab_cl)\n",
    "res_all_bicg = np.array(res_all_bicg)/res_all_bicg[0]\n",
    "res_all_bicgstab = np.array(res_all_bicgstab)/res_all_bicgstab[0]\n",
    "\n",
    "lim = 350\n",
    "plt.semilogy(res_all_bicgstab[:lim], marker='.',color='k', label='BiCGStab')\n",
    "plt.semilogy(res_all_bicg[:lim], marker='x',color='r', label='BiCG')\n",
    "\n",
    "plt.xlabel('Iteration number', fontsize=20)\n",
    "plt.ylabel('Retative residual norm', fontsize=20)\n",
    "plt.legend(loc='best', fontsize=20)\n",
    "plt.xticks(fontsize=20)\n",
    "_ = plt.yticks(fontsize=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## \"Nonlinear GMRES\" or Anderson acceleration\n",
    "\n",
    "We can apply the GMRES-like idea to speed up the convergence of a given fixed-point iteration\n",
    "\n",
    "$$x_{k+1} = \\Phi(x_k).$$\n",
    "\n",
    "This was actually older than the GMRES, and known as an Direct Inversion in Iterated Subspaces in Quantum Chemistry, or **Anderson Acceleration**.\n",
    "\n",
    "Idea: **use history** for the update\n",
    "\n",
    "$$x_{k+1} = \\Phi(x_k) + \\sum_{s=1}^m \\alpha_s (x_{k - s} - \\Phi(x_{k - s})), $$\n",
    "\n",
    "and the parameters $\\alpha_s$ are selected to minimize the norm of the residual\n",
    "\n",
    "$$ \\min_{\\alpha} \\left \\| \\sum_{s=1}^m \\alpha_s (x_{k - s} - \\Phi(x_{k - s})) \\right\\|_2, \\quad \\sum_{s=1}^m \\alpha_s = 1$$\n",
    "\n",
    "More details see in the original [paper](https://users.wpi.edu/~walker/Papers/Walker-Ni,SINUM,V49,1715-1735.pdf)"
   ]
  },
  {
   "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, 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_R$ and/or $P_L$ such that \n",
    "\n",
    "1. Condition number of $AP_R^{-1}$ (right preconditioner) or $P^{-1}_LA$ (right preconditioner) or $P^{-1}_L A P_R^{-1}$ is better than for $A$\n",
    "2. We can easily solve $P_Ly = g$ or $P_Ry = g$ for any $g$ (otherwise we could choose e.g. $P_L = A$)\n",
    "\n",
    "Then we solve for (right preconditioner)\n",
    "\n",
    "$$ AP_R^{-1} y = f \\quad \\Rightarrow \\quad P_R x = y$$ \n",
    "\n",
    "or  (left preconditioner)\n",
    "\n",
    "$$ P_L^{-1} A x = P_L^{-1}f,$$ \n",
    "or even both\n",
    "$$ P_L^{-1} A P_R^{-1} y = P_L^{-1}f \\quad \\Rightarrow \\quad P_R x = y.$$ \n",
    "\n",
    "The best choice is of course $P = A,$ but this does not make life easier.\n",
    "\n",
    "One of the ideas is to use other iterative methods (beside Krylov) as preconditioners."
   ]
  },
  {
   "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\n",
    "3. SOR($\\omega$) (Successive over-relaxation) and its symmetric modification SSOR($\\omega$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Jacobi method (as preconditioner)\n",
    "\n",
    "Consider again the matrix with non-zero diagonal. To get the **Jacobi method** 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",
    "$$ x_i^{(k+1)} = -\\frac{1}{a_{ii}}\\left( \\sum_{i \\ne j} a_{ij} x_j^{(k)} + f_i \\right),$$\n",
    "\n",
    "or in the matrix form\n",
    "\n",
    "$$\n",
    "x^{(k+1)} = D^{-1}\\left((D-A)x^{(k)} + f\\right) \n",
    "$$\n",
    "\n",
    "where $D = \\mathrm{diag}(A)$ and finally\n",
    "\n",
    "$$\n",
    "x^{(k+1)} = x^{(k)} - D^{-1}(Ax^{(k)} - f).\n",
    "$$\n",
    "\n",
    "So, Jacobi method is nothing, but simple Richardson iteration with $\\tau=1$ and left preconditioner $P = D$ - diagonal of a matrix. Therefore we will refer to $P = \\mathrm{diag}(A)$ as **Jacobi preconditioner**. Note that it can be used for any other method like Chebyshev or Krylov-type methods."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Properties of the Jacobi preconditioner\n",
    "\n",
    "Jacobi preconditioner:\n",
    "\n",
    "1. Very easy to compute and apply\n",
    "2. Works well for diagonally dominant matrices (remember the Gershgorin circle theorem!)\n",
    "3. Useless if all diagonal entries are the same (proportional to the identity matrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Gauss-Seidel (as preconditioner)\n",
    "Another well-known method is **Gauss-Seidel method**. \n",
    "\n",
    "Its canonical form is very similar to the Jacobi method, with a small difference. When we update $x_i$ as\n",
    "\n",
    "$$x_i^{(k+1)} := -\\frac{1}{a_{ii}}\\left( \\sum_{j =1}^{i-1} a_{ij} x_j^{(k+1)} +\\sum_{j = i+1}^n a_{ij} x_j^{(k)} - f_i \\right)$$\n",
    "\n",
    "we **use it in the later updates**. In the Jacobi method we use the full vector from the previous iteration.\n",
    "\n",
    "Its matrix form is more complicated."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Gauss-Seidel: matrix version\n",
    "\n",
    "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 with zero on the diagonal.\n",
    "\n",
    "One iteration of the GS method reads\n",
    "\n",
    "$$\n",
    "x^{(k+1)} = x^{(k)} - (L + D)^{-1}(Ax^{(k)} - f).\n",
    "$$\n",
    "\n",
    "and we refer to the preconditioner $P = L+D$ as **Gauss-Seidel preconditioner**.\n",
    "\n",
    "**Good news:** $\\rho(I - (L+D)^{-1} A) < 1, $ where $\\rho$ is the spectral radius,\n",
    "\n",
    "i.e. for a positive definite matrix GS-method always converges."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Gauss-Seidel and coordinate descent\n",
    "\n",
    "GS-method can be viewed as a coordinate descent method, applied to the energy functional\n",
    "\n",
    "$$F(x) = (Ax, x) - 2(f, x)$$\n",
    "\n",
    "with the iteration \n",
    "\n",
    "$$x_i := \\arg \\min_z F(x_1, \\ldots, x_{i-1}, z, x_{i+1}, \\ldots, x_d).$$\n",
    "\n",
    "Moreover, the order in which we eliminate variables, is really important!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Side note: Nonlinear Gauss-Seidel (a.k.a coordinate descent)\n",
    "If $F$ is given, and we optimize one coordinate at a time, we have\n",
    "\n",
    "\n",
    "$$x_i := \\arg \\min_z F(x_1, \\ldots, x_{i-1}, z, x_{i+1}, \\ldots, x_d).$$\n",
    "\n",
    "**Note** the convergence result for block coordinate descent for the case of a general functional $F$: \n",
    "\n",
    "it converges locally with the speed of the GS-method applied to the **Hessian** $$H = \\nabla^2 F$$ of the functional.\n",
    "\n",
    "Thus, if $F$ is twice differentiable and $x_*$ is the local minimum, then $H > 0$ can Gauss-Seidel converges."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Successive overrelaxation (as preconditioner)\n",
    "\n",
    "We can even introduce a parameter $\\omega$ into the GS-method preconditioner, giving a **successive over-relaxation** (**SOR($\\omega$)**) method:\n",
    "\n",
    "$$\n",
    "x^{(k+1)} = x^{(k)} - \\omega (D + \\omega L)^{-1}(Ax^{(k)} - f).\n",
    "$$\n",
    "\n",
    "$$P = \\frac{1}{\\omega}(D+\\omega L).$$\n",
    "\n",
    "- Converges for $0<\\omega < 2$. \n",
    "- Optimal selection of $\\omega$ is **not trivial**. If the Jacobi method converges, then \n",
    "$$\\omega^* = \\frac{2}{1 + \\sqrt{1 - \\rho_J^2}},$$ \n",
    "where $\\rho_J$ is spectral radius of Jacobi iterations\n",
    "- Note that $\\omega = 1$ gives us a Gauss-Seidel preconditioner.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Preconditioners for sparse matrices\n",
    "\n",
    "- If $A$ is sparse, one iteration of Jacobi, GS and SOR method is cheap (what complexity?). \n",
    "\n",
    "- For GS, we need to solve linear system with a sparse triangular matrix $L$, which costs $\\mathcal{O}(nnz)$.\n",
    "\n",
    "- For sparse matrices, however, there are more complicated algorithms, based on the idea of **approximate** LU-decomposition.\n",
    "\n",
    "- Remember the motivation for CG: possibility of the early stopping, how to do approximate LU-decomposition for a sparse matrix? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Remember the Gaussian elimination\n",
    "\n",
    "- Decompose the matrix $A$ in the form \n",
    "\n",
    "$$A = P_1 L U P^{\\top}_2, $$\n",
    "\n",
    "where $P_1$ and $P_2$ are certain **permutation** matrices (which do the pivoting).\n",
    "\n",
    "- The most natural idea is to use **sparse** $L$ and $U$. \n",
    "\n",
    "- It is not possible without **fill-in** growth for example for matrices, coming from 2D/3D Partial Differential equations (PDEs).\n",
    "\n",
    "- What to do?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Incomplete LU\n",
    "\n",
    "- Suppose you want to eliminate a variable $x_1$, and the equations have the form\n",
    "\n",
    "$$5 x_1 + x_4 + x_{10} = 1, \\quad 3 x_1 + x_4 + x_8 = 0, \\ldots,$$\n",
    "\n",
    "and in all other equations $x_1$ are not present. \n",
    "\n",
    "- After the elimination, only $x_{10}$ will enter additionally to the second equation (new fill-in).\n",
    "\n",
    "$$x_4 + x_8 + 3(1 - x_4 - x_{10})/5 = 0$$\n",
    "\n",
    "- In the Incomplete $LU$ case (actually, ILU(0)) we just throw away the **new fill-in**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Incomplete-LU: formal definition\n",
    "\n",
    "We run the usual LU-decomposition cycle, but avoid inserting non-zeros **other** than the initial non-zero pattern. \n",
    "\n",
    "```python\n",
    "    L = np.zeros((n, n))\n",
    "    U = np.zeros((n, n))\n",
    "    for k in range(n): #Eliminate one row   \n",
    "        L[k, k] = 1\n",
    "        for i in range(k+1, n):\n",
    "            L[i, k] = a[i, k] / a[k, k]\n",
    "            for j in range(k+1, n):\n",
    "                a[i, j] = a[i, j] - L[i, k] * a[k, j]  #New fill-ins appear here\n",
    "        for j in range(k, n):\n",
    "            U[k, j] = a[k, j]\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## ILU(k)\n",
    "\n",
    "- Yousef Saad (who is the author of GMRES) also had a [seminal paper](https://www-users.cs.umn.edu/~saad/PDF/umsi-92-38.pdf) on the **Incomplete LU** decomposition\n",
    "\n",
    "- A good book on the topic is [Iterative methods for sparse linear systems](https://www-users.cs.umn.edu/~saad/IterMethBook_2ndEd.pdf) by Y. Saad, 2003\n",
    "\n",
    "- And he proposed **ILU(k)** method, which has a nice interpretation in terms of graphs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## ILU(k): idea\n",
    "\n",
    "- The idea of ILU(k) is very instructive and is based on the connection between sparse matrices and graphs.\n",
    "\n",
    "- Suppose you have an $n \\times n$ matrix $A$ and a corresponding adjacency graph.\n",
    "\n",
    "- Then we eliminate one variable (vertex) and get a smaller system of size $(n-1) \\times (n-1)$.\n",
    "\n",
    "- New edges (=fill-in) appears between high-order neighbors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## LU & graphs\n",
    "\n",
    "- The new edge can appear only between vertices that had common neighbour: it means, that they are second-order neigbours.  \n",
    "- This is also a sparsity pattern of the matrix $A^2.$\n",
    "- The **ILU(k)** idea is to leave only the elements in $L$ and $U$ that are $k$-order neighbours in the original graph.\n",
    "\n",
    "- The ILU(2) is very efficient, but for some reason completely abandoned (i.e. there is no implementation in MATLAB and SciPy).\n",
    "\n",
    "- There is an original [Sparsekit](https://people.sc.fsu.edu/~jburkardt/f_src/sparsekit/sparsekit.html) software by Saad, which works quite well."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## ILU Thresholded (ILUT)\n",
    "A much more popular approach is based on the so-called **thresholded LU**.\n",
    "\n",
    "You do the standard Gaussian elimination with fill-ins, but either:\n",
    "\n",
    "- Throw away elements that are smaller than threshold, and/or control the amount of non-zeros you are allowed to store.\n",
    "\n",
    "- The smaller is the threshold, the better is the preconditioner, but more memory it takes.\n",
    "\n",
    "It is denoted ILUT($\\tau$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Symmetric positive definite case\n",
    "\n",
    "- In the SPD case, instead of incomplete LU you should use Incomplete Cholesky, which is twice faster and consumes twice less memory.\n",
    "\n",
    "- Both **ILUT** and **Ichol** are implemented in SciPy and you can try (nothing quite fancy, but it works)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Second-order LU preconditioners\n",
    "\n",
    "- There is a more efficient (but much less popular due to the limit of open-source implementations) **second-order** LU factorization [proposed by I. Kaporin](http://www.researchgate.net/profile/I_Kaporin/publication/242940993_High_quality_preconditioning_of_a_general_symmetric_positive_definite_matrix_based_on_its_UTU__UTR__RTU-decomposition/links/53f72ad90cf2888a74976f54.pdf)\n",
    "\n",
    "- The idea (for symmetric matrices) is to approximate the matrix in the form\n",
    "\n",
    "$$A \\approx U_2 U^{\\top}_2 + U^{\\top}_2 R_2 + R^{\\top}_2 U_2,$$\n",
    "\n",
    "which is just the expansion of the $UU^{\\top}$ with respect to the perturbation of $U$. \n",
    "\n",
    "- $U_1$ and $U_2$ are upper-triangular and sparse, whereare $R_2$ is small with respect to the drop tolerance parameter.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Summary of this part\n",
    "\n",
    "- Jacobi, Gauss-Seidel, SSOR (as preconditioners)\n",
    "- Incomplete LU, three flavours: ILU(k), ILUT, ILU2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Next lecture\n",
    "\n",
    "- Toeplitz and circulant matrices, FFT\n",
    "- Matrix functions and their applications"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Questions?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "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": 41,
     "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": {
  "anaconda-cloud": {},
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  },
  "nav_menu": {},
  "toc": {
   "navigate_menu": true,
   "number_sections": false,
   "sideBar": true,
   "threshold": 6,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}