Jacobi.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chebyshev-PSOR Jacobi method for solving linear equations\n",
    "\n",
    "Compairions of Jacobi method, Constant SOR, and Chebyshev PSOR: \n",
    "\n",
    "The linear equation $$P x = 0$$ is assumed without loss of generality where\n",
    "$$P = I + M^T M \\in \\mathbb{R}^{512 \\times 512}$$ is a diagonally dominant matrix.\n",
    "Each element of $M$ follows ${\\cal N}(0, 0.03^2)$. The minimum and maximum \n",
    "eigenvalues of $B = I-A$ are given by \n",
    "$\\lambda_{max}(B) = 1.922$ and $\\lambda_{min}(B) = 0.6766$, respectively.\n",
    "\n",
    "The 3-curves represent the errors in Euclidean norm from the fixed-point, i.e., $||x^{(k)} - x^*||$ where\n",
    "$x^*  = 0$.\n",
    "The plain Jacobi method, Jacobi method with the optimal constant SOR factor, \n",
    "and Jacobi method accelerated by the Chebyshev PSOR are compared. \n",
    "\n",
    "The optimal constant SOR factor \n",
    "is given by $$\\omega_{SOR} = 2/(\\lambda_{min}(B) + \\lambda_{max}(B)) = 0.7697$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Required packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<torch._C.Generator at 0x7fbd9f38c990>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import torch\n",
    "import torch.nn as nn\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "torch.manual_seed(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Global constants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 512 # dimension of matrix \n",
    "max_itr = 25 # maximum number of itaration\n",
    "T = 8 # period of Chebychev step"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generation of symmetric diagonally dominant random matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = torch.normal(torch.zeros(n,n), 0.03)\n",
    "M = M.t() @ M\n",
    "P = M + torch.eye(n)\n",
    "D = torch.diag_embed(torch.diagonal(P))\n",
    "R = P - D\n",
    "Dinv = torch.inverse(D)\n",
    "B = -Dinv @ R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluation of eigen value of $I-B$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lmax = 1.921623706817627\n",
      "lmin = 0.676575243473053\n",
      "condition number =  2.8402217275249186\n"
     ]
    }
   ],
   "source": [
    "lmax = torch.max(torch.eig(torch.eye(n)-B)[0].sum(dim=1)).item()\n",
    "lmin = torch.min(torch.eig(torch.eye(n)-B)[0].sum(dim=1)).item()\n",
    "print('lmax =', lmax)\n",
    "print('lmin =', lmin)\n",
    "print('condition number = ', lmax/lmin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## q = 0 is assumed without loss of generality"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = torch.zeros(n,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Jacobi method\n",
    "\n",
    "Fixed-point iteraiton of Jacobi method is given by \n",
    "$$\n",
    "\tx^{(k+1)} = D^{-1} \\left(q - (P - D) x^{(k)}  \\right).\n",
    "$$\n",
    "where $D$ is the diagonal matrix whose diagonal elements are identical to the diagonal elements of $P$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Jacobi(nn.Module):\n",
    "    def __init__(self, P, b):\n",
    "        super(Jacobi, self).__init__()\n",
    "        D = torch.diag_embed(torch.diagonal(P))\n",
    "        self.R = P - D\n",
    "        self.Dinv = torch.inverse(D)\n",
    "        self.b = b\n",
    "\n",
    "    def forward(self, num_itr):\n",
    "        n = self.R.size()[0]\n",
    "        x = torch.normal(torch.zeros(n,1), std=1.0)\n",
    "        for i in range(num_itr):\n",
    "            x = self.Dinv @ (self.b - self.R @ x)\n",
    "        return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SOR Jacobi method\n",
    "\n",
    "A constant SOR factor $\\omega$ is assumed.\n",
    "The SOR iteration is given by\n",
    "$$\n",
    "x^{(k+1)} = x^{(k)} + \\omega\\left( f( x^{(k)} ) - x^{(k)}  \\right).\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "class SOR_Jacobi(nn.Module):\n",
    "    def __init__(self, P, b, omega):\n",
    "        super(SOR_Jacobi, self).__init__()\n",
    "        D = torch.diag_embed(torch.diagonal(P))\n",
    "        self.R = P - D\n",
    "        self.Dinv = torch.inverse(D)\n",
    "        self.b = b\n",
    "        self.omega = omega\n",
    "    def forward(self, num_itr):\n",
    "        n = self.R.size()[0]\n",
    "        x = torch.normal(torch.zeros(n,1), std=1.0)\n",
    "        for i in range(num_itr):\n",
    "            tmp = self.Dinv @ (self.b - self.R @ x)\n",
    "            x = x + self.omega*(tmp - x)\n",
    "        return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Chebychev-PSOR Jacobi method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Chebychev root\n",
    "def c_root(k, a, b, T):\n",
    "    return  (a + b)/2.0 + ((b - a)/2.0) * math.cos(math.pi * (2*k+1)/(2*T))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Chebychev_SOR_Jacobi(nn.Module):\n",
    "    def __init__(self, P, b):\n",
    "        super(Chebychev_SOR_Jacobi, self).__init__()\n",
    "        D = torch.diag_embed(torch.diagonal(P))\n",
    "        self.R = P - D\n",
    "        self.Dinv = torch.inverse(D)\n",
    "        self.b = b\n",
    "    def forward(self, num_itr, T):\n",
    "        n = self.R.size()[0]\n",
    "        x = torch.normal(torch.zeros(n,1), std=1.0)\n",
    "        for i in range(num_itr):\n",
    "            tmp = self.Dinv @ (self.b - self.R @ x)\n",
    "            gamma = 1.0/c_root(i % T, lmin, lmax, T) # Chebychev step\n",
    "            x = x + gamma * (tmp - x)\n",
    "        return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_pure_jacobi(P, q, color):\n",
    "    norm_list = []\n",
    "    jacobi = Jacobi(P, q)\n",
    "    for i in range(max_itr):\n",
    "        r = jacobi(i)\n",
    "        norm_list.append(torch.norm(r))\n",
    "    plt.plot(range(max_itr), norm_list, '.-', label='Jacobi', color=color, linestyle = \"solid\", linewidth=1, marker=\"s\") "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_sor_jacobi(P, b, color):\n",
    "    norm_list = []\n",
    "    omega_opt = 2.0/(lmax+lmin) # optimal constant SOR factor\n",
    "    print('SOR omega = ', omega_opt)\n",
    "    sor_jacobi = SOR_Jacobi(P, q, omega_opt)\n",
    "    for i in range(max_itr):\n",
    "        r = sor_jacobi(i)\n",
    "        norm_list.append(torch.norm(r))\n",
    "    plt.plot(range(max_itr), norm_list, '.-', label='Constant SOR (optimal)', color=color, linestyle = \"solid\", linewidth=1, marker=\"^\") "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def plot_chebychev_sor_jacobi(T, P, b, label, color):\n",
    "    chebychev_sor_jacobi = Chebychev_SOR_Jacobi(P, q)\n",
    "    c_norm_list = []\n",
    "    for i in range(max_itr):\n",
    "        r = chebychev_sor_jacobi(i, T)\n",
    "        c_norm_list.append(torch.norm(r))\n",
    "    plt.plot(range(max_itr), c_norm_list, '.-', label=label, color=color, linestyle = \"solid\", linewidth=1, marker=\"o\") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot of error curves (iteration step v.s. error in Euclidan norm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SOR omega =  0.7697639935449304\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEGCAYAAABGnrPVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOydd3iUxdbAf5NADAm9BSlJAAkCAULooJAgqCiKKIoUC6ioeEGsKNgFvCpeAYUrVbCBDXvDaLiCop8gKJ2oQEB6aAkQ0s73x+xuNsnuZgnZlp3f88zz7s5b5swm2ZOZ05SIYDAYDAZDeRHiawEMBoPBULEwisVgMBgM5YpRLAaDwWAoV4xiMRgMBkO5YhSLwWAwGMqVSr4WwNvUrVtXYmNjy3TvyZMniYyMLF+BAoRgnjsE9/yDee4Q3PO3n/vatWsPi0g9d+4LOsUSGxvLmjVrynTvihUrSEpKKl+BAoRgnjsE9/yDee4Q3PO3n7tSape795mtMIPBYDCUK0axGAwGg6FcMYrFYDAYDOWKUSwGg8FgKFeMYjEYDAZDuWIUSyk0aNCAYUqxU4XSK7kPO1Uow5SiQYMGvhbNYDAY/BKjWEqhz4EDzANiKSAEIZYC5ln6DQaDwVASo1hKYSpQPDQq0tJvMBgMhpIYxVIK0WfZbzAYDMGOUSylkO7kIxIg98Zh8N13UFDgXaEMBoPBjwl4xaKUqqGUWqiU2u2J50/kJk4SUaTvJBGMYQoPv38Bu4fehTRtCo8/Dn/+WXjR229DbCyEhOjj2297QjyDwWDwO/xKsSilpimldiilRCkVb9cfp5RarZTabjm2sJ4TkeMiMgrY5gmZ3q/0AHcwl53EUIBiJzHcwVwWhF5F4xeGcOHJr+l+8N/8tjIN6dkTLr4Y7rhDt127QEQfR482ysVgMAQFfqVYgI+BXkDxZGevAbNEJA6YBczxlkC5uW15R4YTKzv5IfV7YmUn78hw8vLa8sADbTh+vBmdbruQrj/OoLG8wc8XD+bU/Plw+nTRB506xe6bb/aW2AaDweAz/EqxiMgqESmypaWUqg8kAkssXUuARKWUW+mbPU2lSopXX+3A4cN16NS9Dj3+fTPhTq5tXFAAX3wBmZklT5qtM4PBUEFQIuJrGUqglNoJDBCRjUqpjsAbItLG7vxmYISI/GZ5PwsYCHwOPC8iO4o9bzQwGiAqKqrj0qVLyyRXVlYWVatWdXnNP/8U0HPE9cRypMS5owAdOlB9yxaymjfnaGIixzp04LwDB4ibPp3QM2ds1+afdx7bHnyQg337lknW8saduVdkgnn+wTx3CO752889OTl5rYh0cutGEfG7BuwE4i2vOwKbip3fDCSW5dkdO3aUspKamurWdUNBsggT0RYWEcv7oSD5+fkip06JpKSITJwo0q2biFJFrrW1mJgyy1reuDv3ikowzz+Y5y4S3PO3nzuwRtz8nvWrrTAn7AYaKaVCASzHhpZ+v2QJDbiD14oZ/BeyhBc4//w23Dl+PF/l5HDmiSdg9WoKnKwaC3btArtVjMFgMAQCfq9YROQgsB4YaukaCqwTkUO+k6o0HmcJQ2nKTkIpoCk7WcJg4ErOnNnItm238+STrxEVFcWQIUNId/KUHICoKBgyBJYsgWPHil5g7DIGg8EP8SvFopSaqZTaAzQGUpRSmyyn7gLGKqW2A2Mt7/2WSpUuhhIm/POoVCmfP/4IpX37zvz55ycMHLiXDh2uYyJwstjVJ4FRANu3w6WXwjvvQHS0fj17No9Wq8bJESOKuDSfHDGCu2vU8MIMDQaDwTl+VfNeRMYB4xz0bwW6el+ispGb29bJGd0/YwZMmgQzZkTw4os3cIQc4BGmcopojpFOTSZylCVA+rXX0qpVK1olJ9Pm1ltJOHCA+j/9xOSsLEKLPT0SmHDihOcmZjAYDG7gV4olmKhfH6ZMgYcfhpo1N7GEdSyhNhAKnAKaAQeYPHkyW7ZsYcuWLXz11Vds3bqVI0eO4MBhGXCdw6xBgwYccJCVOSoqiv3795/znAwGgwGMYvE5eufq30BTYCRasYQD7wEDSUpKIikpqcg9mZmZpFevTqyD5wnwbatW7L3yShpedhlt27YlKioKpZRDpQI47TeKyGAwlAW/srEEK3XrtgVuAipbekKAHsAuhjnIc1mtWjWndpn7gDrNmzPwtddodN11PHfBBTStV48+ffoA2vNhB5BvOQ7FOWeriMpKgwYNUEqVaKaYmsEQmJgVix9w/fV/sGAB5OQU9oWFVWLEiOq0bw/33QdZWTBqFNx6KzRuXJiGYCohRCOko5hIAUuAmZ9/Dnl51Pz6a6bPm8fLP/zAP+HhvAZMpLC+TCwwz/K6ffv21K5dmzp16lC7dm1q167tlblD2RSYWU0ZDP6LWbH4AatXF1UqoN//9huMGwe//w7vvgt79kC7dnDFFVCjxiiWUJmmvEIouTRlJkvQX6wAVKoEAwagPvmEkI0badKzJ0/hvGjZ4sWLefzxx7nxxhvp2LEjNUrxLps+fTqrVq3i5Mmi6yZvrT7Ksq1nVkUGg3cwKxY/YN061+eVgs6ddfvPf+CDD+DkyQVs2rSAY8cgPx+qVLmHv/++B4ffk40awaRJhD72mMPnRwOxCQkl+idOnOhUprS0NJYsWcLGjRtp2rQpnTp1onPnzi6/8Ddu3Eh6enqJ5or4+HjbCsq+nS3e2tYzGAxGsQQcERFw8826DRsG772n+7Ozdd/nn0NYmON7/wkJoYmDomQnlKJmVhacRT6kWbNmAZCTk8PGjRtZs2YNa9ascXnPDTfcQHR0tK1ddtllREdHl3BOsGfJkiUcOXKkSDt69KjLcS644ALq1q1LvXr1bM1gMHgPo1gClH374KOP9GoFdIzkd99Bw4aF5WCaNSt6T5M33tB1YU6dKuwMD6dmQgI0bw733w/33GNTMFFRUU7tGFbCwsJITEwkMTGR0aNHM2/evBLXW9m8efNZz7NtW8cxQc8995zTe7766isOHTrEoUOHOHz4MIcOuU7ScPToUWrVqnXWshkMBscYG0uA8uyzJSsiV6qkA/NzcqBLF+jfHz75BPLyLBcMHw5z55LXSOcwy2scA/PnayNPaiqsX68VzPPPQ1YW+/fvtyWVS01Ntb0ub+O4vaJyp780WrRoQY8ePRg4cCC33XYbjzzyiMvrY2Ji6NGjB08//TS//PIL+RZtbW+XSU5ONnYZg8FNjGIJUJwZ/LdsgZdegt279VbZCy/oNGJPPw3//AMMH864q3dSOaSAe6/eqZUNQOvWOh9ZcQWzcCHExtK7Tx+P5SOzV2D2zZUCK09ldPDgQZ555hkyMzO5/fbbiYqK4sYbbzR2GYOhrLibBrmiNG+kzfc3fv9dZMwYkVq1RC6/XCTMktG/ShWRffuc3LRpk07pXzyVf0SEyFtvOR0rKipK0HGaRVpUVJRnJucmZyPX7t27ZcGCBQ6vtzZf4u3P2N3fe3/92Z8rgfp3Xx5U5LT5hnOkXTuYNQvS03XF5Nxc3Z+TA04cxfQKZt++kv2nTulEZ04oy+rDG5yNXI0bN2bUqFEunzdy5EgWLFjAtm3brDWCAO+4NfvrSqqs8UjGDbziYRRLEJGZCb/8opceoA3/Cxbo3bD16x3c4MwVuBQX4WCga9eupKamctlll1G/fn2uueYapk2b5vEv/QIHXn32ZGRkOOz39Rd4TvF9WwsmHunsCJT5G6+wIMKRwT8sDP78E666CmJitFPYdddZXJajo3Va/uIoBY8+Co8/rv2fg5C77rqLu+7S1Rv27NnDjz/+yKpVq1zek5+fT2ho0ZzUpWUQOH36NGvWrGHVqlX8+OOPrF692uUYzZo1IzQ0lBYtWnDBBRfQokULWrRo4bE8cbm5ufz888+kpKS4lCsyMpJq1aoRFRVFVFQU9evXL5M9zF9Xa94iUOZvFEsQ4czgn5MDO3bAp5/Cq69qr+M77oDxD06h9oRi7skREdoj4McfoU0bfcOVV3p3Il7CHXdr0FtnQ4YMYciQIbz66qtOn1etWjWaN29Oy5Ytbc3VF0X37t35448/aNOmDT179uTWW29l3rx5NGzY0OkYx44d4/Dhw6SlpfHnn3+SlpbGZ5995nKe69ato0WLFkXqupcW6JqSkkJKSgorV66kRYsW9O3b1+UYZ86c4ciRIxw4cICDBw9y4MCBUr8MW7ZsaVNEDRo0KLOXoEn/432U/f5wMNCpUycpLZDPGStWrHAZzFdR2LQJZs/WTmJPtniboRsmUfd0OkerRlPntSmFnmTffgtjxmgjzowZOolZBcXdn71Syum5zMxMtm/fzrZt22xt6dKlTq///vvv6dKlC5GRRRPxuBrD2d+zq3vatWtHWloaNWvWpEWLFsTFxTF//nyn1zdr1oy+ffvSt29f+vTpQ506dTwi15YtW9i/f79NCe3fv99l/NKwYcOKBOBGR0cTExPjMj2RM7n8VRmV5TM+F+x/75VSa0Wkk1s3umvlrygtGL3Cysrx4yKTJ4sopR3CKlcW2bKl2EWnT4s8+aRInToi//mPSG6u9hqLidE3xsS49CILFNz92XOWnmRne71I2byvShsnPz9f0tPT5bvvvpP//ve/fiPX2Vz/xhtvyOTJk2X06NFy+eWXS6tWrSQyMtLlPRkZGVJQUHDOcnnaIy4rK0umTJniUi5H8zhXyuoV5vMvem83o1jOjrvvLnRPDgnRr2++WWTlSpEiv8fbtolccolIkyYi4eFyNi7KgYCnXG7L8gVeFryh8MqCpz+vgoICl/dUr15dqlevLu3atZOBAwfK+PHjZfr06X7zeZ05c0ZeeeUVadCggQwZMsTlOF27dpXly5eXq4Ipq2IxNhaDU/btg9dfL7TLFBRASIiOk7z9dv169Gi46SaoExent8bq19eJy+yxuihbt9AqMP66Z++uvcjbnO3ndbbzcLV1BHD8+HGOHj3Kjh072LlzJzt27CAtLc3lPW3atCEiIqJIK2/y8/N55513ePLJJ2nZsiVffvklHTp0YMWKFU7nP378eMaOHUv9+vV59tln6d27d7nL5S4BrViUUjWAl4F+ItLE1/JUNBx5kRUUwOHDOsJ/5UqYOxeeekrb70ePVvTKyMDhn7JxUXaIt77wPf0F7i08obhr1apFrVq1SExMtPVZk6w64r333uPUqVNF2rJly5xe/9BDD9GrVy8uuuiiIjnpnNlxatasSePGjalWrRqvv/56EQVR2vwHDx7MO++8w6hRo2jatCnr1q3jyJEjJa7zuK3I3aVNeTVgGrp4oQDxdv1xwGpgu+XY4iyemeLutWYrzH0SEoruaFlbQkLR6zIyRKZPF2ndWmRPpRjHN8XE+GIK5Uaw/eztCfS5+9r288wzz0jfvn2latWq0r59exk3bpx88MEHLu/59NNPz2lLKycnR+bNm3fOW3SBFHn/MdALKB4g8RowS0TigFnAHOsJpVRzpVRKsfaQ90QOTtatK9QMqakrbK+L14+pXRvuvRc2boSTk6aQHVp0a0BCK8GUKV6U3GAoxNe56B5//HG+/fZbMjIyeO2112jYsCELFy50ec9VV11V6jaeKypXrsztt99e5vvPFa9vhYnIKii696mUqg8kAv0sXUuAV5VS9UTkkIj8Bbh2lDf4HKUg7qnhLFwFl3w3iSaks4fGhOQX8PsbGbTqXjKVv8Hgj9grHXdczd0tMdGtWze6devGhAkTzklx+Ds+i2NRSu0EBojIRqVUR+ANEWljd34zMEJEfivlObOAgcDnwPMissPBNaOB0QBRUVEdXcUOuCIrK6tIEFkw4e7cMzLCGDasKzk5hRHmF1T6m19CuzM6dB5/xV1M//776dXrEOHhrtOT+BPmZx+ccwfPzT85OdnpudTUVL8Yw37uycnJ/h/HAuzEYmMBOgKbip3fDCSW97jGxlI23J27vXuytYWFibw46EcpqFdPvnlpg/TvrzMtjx4t8vPPhW7Le/eK9OrlIuOyDzE/++DFU/PHCy7d5zpGINlYHLEbaKSUCgWwHBta+g0BhLO0MW/v6IF6+WUufeUqvlx0kD/+0LnJhg+H+HhdQ2biRFi1SnujGQwVnfIucOerMRxRqmJRSrVUSr2plPpJKfV/1laeQojIQWA9MNTSNRRYJyKua8oa/A57g799W7cOrUVGjIBBg2hcN5uJEyEtDf77X511edEi7c48b55OjGkwVGS8UWLCV2Us3FmxvA/8ATwGPGTXyoRSaqZSag/QGEhRSm2ynLoLGKuU2g6Mtbw3VDSefhoaNtRZLkVQCnr1grp1LRmV0aWUW7eG66+HDz/UNWQMBkPg4I5iyRORF0XkexH5n7WVdUARGScijUWkkog0EIvBXkS2ikhXEYmzHLeVdQyDHxMSAosXw9atMHUqUDLCXwRCQ6FbN72aadhQL3Q+/7zkNtu+fdC7N/hpwLvBEJS4o1i+Vkpd7nFJDMFDRITO0T9nDnzwgdMI/7//hpQUrYO6d4fnn4fzz9fpZFJS9Mrm2WeNXcZg8DfcUSwpwPtKqeNKqYNKqUNKqYOeFsxQwTn/fPjkE7j7bo5/t8ahwf+nn/TrqChdgGzlSl3pslUrXWesQQOdUqagQK94zKrFYPAP3FEsc4GRQAegM9DJcjQYzo0OHWDePN4+eQ0y8xUkJhZRIfr41tslIvwBmjSBBx6AX3+Fyy4r7D99GgYMgG1mA9Vg8DnuRN4fEZEPPC6JITi55hpYulTnhLEG6+7apdMmg9OMyPv2wbJlkJ9f2Pf773DxxdqNecQIuPFGvdoxGAzexZ0Vy8dKqbuUUrWVUhHW5nHJDMHD6tWFSsWKNdW+ExzZZUJC4LrrdFqytWuhZUu4/HJ46y3IytLXGGO/weB53FEsk4HZwGEgE8iyHA2G8mG3kzhYF6n2nQVi/vwzXHopvPEG/PMP3HKLXhA1bqwXP7ffboz9BoOncalYlFIhQAcRCbG0UOvRS/IZgoHo6LPrp5RATAuRkTB0qHZTTkuDNm3gq6/0SmfuXNiwoZznYTAYgFIUi4gUAPO9JIshWJkyRbsg2xMRUa6p9uvVgz17oHJl/b6gADp10tUvf/213IYxGAy4txW2RSkV62E5DMHM8OF6CRETo3Pvh4ZqY345ljJ2Vma5aVO44QYdjPn22yW31wwGw9njjmKpB/yhlPpSKfWetXlaMEOQMXw47Nypv/GXLYOPPirXb3lnQZgZGTov2cSJWvHExOisM/bG/X374N57E4zB32BwE3cUy1J07q53gS/smsHgGa66Cpo3h+nTy+2Rzoz9P/2kF0hXX62j+VNStFJp3Vq7LP/yi1ZKGzbUMAZ/g8FNSo1jEZHF3hDEYLChFMyYAV27aut7kybn/EhHwZaOaNNG5yebOlWvYK6/XttmRBQLF8Ljj+uIf4PB4Bx30ubXVUottaZyUUq9o5Sq5w3hDEFM8+bwr3/pMHsfUKsW3H8/XHGFXtEAZGdDjx4643J2tk/EMhgCAne2wuYA24EEdF36NEufweBZJkyANWvg2299Mvy+fToRc15eYd+ePfCf/+iMy7fdBt9/XzT63/5eE4hpCFbcUSzNReQJEflHRPaIyJNAM08LZjBQpQrMnKlXLmfOeH14RwZ/pSAhQcfAtGqlF1TR0fDgg4WxNdZ7TSCmIVhxR7GEKKXqW99YXvtLSWNDRWfAAJ2b5eWXvT60K4N/o0aFyuTbbyE8HK69VttoHnkEFi40WZcNwYs7CmIasE4pNVcpNQdYC7zgWbEMBjtmzIBp01ymePEE9tH9qakrHEb3g/YgmzxZ14+ZP19H+lsXWLm52n3ZYAgmSlUsIvIGcBm6PPFG4DIRecvTghkMNpo2hXHjtDXdj1FKi/rXX4V9eXm6ntnjj+uYGYMhGHBrS0tENorIqyLyiohs9rRQBkMJHn5YV/n65htfS+ISR3aZSpXg/ffhggtg5Ejtj2AwVGTccTfuoZRapZTa628VJJVSrZVSc5RS85VSi5RSytcyGTxEeLjeEhs71ieGfHdxZJfJzdV+CGlpcOGFOrV/t246nb91KsaLzFCRcGfFsgCdNv8iyqGCpFJqmlJqh1JKlFLxdv1xSqnVSqntlmOL0p4lIptF5E4RuR2ItDRDReXKK7VB46WXfC2JU1xlXa5bV3tQ//23TiHzxhvao2zSJL0gM15khoqCO4rltIi8IyJ/i8guazuHMT8GegHFn/EaMEtE4oBZ2MXKKKWaK6VSirWHLOeSlVJvo+vFnDoHuQyBwPTpOpBk17n8CvoWawqZ5cvhhx/gwAG9eiko0Mb/PXt8LaHBcG64o1i+VEr1L68BRWSViBSp7GRxYU4Elli6lgCJ1gh/EflLRPoWay9azqWKyHAgDx3EaajIxMbC+PFw332+lqRcaNkSwsJ0A71tFhenVzZbt/pWNoOhrCgpXhK2+AVKHQLqoKtGngEUICJS3+WNpQ2s1E5ggIhsVEp1BN4QkTZ25zcDI0TkNxfPSAIGW2SqDPxLREqkxFVKjQZGA0RFRXVcunRpmWTOysqiatWqZbo30PGnuaucHDrfdht//utfHOna1Stjemr+GRlhDBvWlZycwtp5lSsXcMUVe1m5sh7nn5/N5ZfvJzn5IJGRDkL8vYA//ex9QTDP337uycnJa0Wkk1s3iojLBsQ4aqXd58ZzdwLxltcdgU3Fzm8GEs91nOKtY8eOUlZSU1PLfG+g43dz/+orkfr1RZo0EVFKJCZG5K23PDacp+Z/990iYWFFLTJhYSJjxojk5op89pnIoEEiNWqI3HyzSGqqSH6+vnfvXpFevUT27fOIaDb87mfvZYJ5/vZzB9aIm9+z7mQ39sZm9m6gkVIqVETylVKhQENLv8FQkowM3ayJunbtgtGj9etyLBDmaVxF91eqpBMPDBgABw/qQmRjx8KpU9ptedu2QoP/rFm+kd9gcIRfpGYRkYPAemCopWsosE5EDvlOKoNfM2lSyeyPp07p/gDClReZPfXra7PSH3/Au+9q12V7g/8///hGfoPBEV5XLEqpmUqpPUBjIEUptcly6i5grFJqO7qw2F3els0QQDhL77JrF0yZor+Bi9sP335bG/9DQvTx7bc9LWW5oxR06gSRkUUN/hdeqKd94IBv5TMYwI1CX+WNiIwDxjno3wp4xxJrCHyiox27HNevr/eNrrlGr2ise0kHDsA99+hVDQTs1hnoYMrXXy/cQhPRymXzZq1gLrsM7r4bevXSishg8DZurViUUpcopUYrpcZYm6cFMxhcMmUKREQU7YuI0DEuM2bohF1ffaWL2E+dCqNGFSoVKwG4dQaO08aIQM2asGMH9OypFUt8PLz6Khw/XnidifA3eAN3UrosAmZQGHlvjb43GHzH8OEwd65WHErp49y5hasPpXSU/sMPw8qVzp/j5YzJ5YErg3/NmtrAv2mTNuivXKl3/e68U6daM3ViDN7Ana2wHkAbEcn1tDAGw1kxfLj721jOts4AkpLghht0QZUAKGhf3LDvCKX0tJKS9OpkwQKdEWffPr26WbhQZ1wOgOkaAhB3tsKMy68h8HG2dbZwoXa3+uknbaBISoLZswv3iiwG/959+gSswb9BA73jd9VVOp0MQHY2tGundwm3bPGtfIaKhzsrlu3Ad0qpj4Fsa6eIzPaYVAZDeWNd2UyapLe/oqO1srH2Dxyov22/+UbnuJ84Ec4/X9tqcnNREPAG/8WLdX0YK8ePw59/Qr9+UK0aDBqkF20dOxYa/fftg3vvTeCbb8zqxuA+7iiWcOAvoK1dn+s8MAaDP1La1ll4uFYwViUTHa3dreyxGvwDTLE4MviDTuefnq5rxHz0kZ7W6dNayQwaBEuXwoYNNUwQpuGscCfyfqQ3BDEY/IrwcDh82PG5Xbv0dtl110FUlHflKiOuDP4hIdCli27WrbGPPtJFOzdsAFAsWACPPaYXcQZDabjrbtxSKXWDUupma/O0YAaDz4mOdtxfrx78+KNOTZycDP/9b9HIRD8MxHQ3wt/qTDdpElx0EVSurPvPnIG2bfVUT5zwvvyGwMIdd+NxwDJ0vZThluNQlzcZDBUBZwb/l1/WykIbILT/bsuW0KePTuJ1xx16VSNSaJfxA+VyNliDMO13AjMz4YsvtGf3HXfAr7+WTG5gMIB7K5bRQBcgXUQus7w+6lGpDAZ/wC5WRhzFylSpoiP8rUpm3Dh47z1tpLAnAAMxndlkYmL0VlmzZjBkiDb0z5mjlY4VE4RpcEexZIvISSBEKaVEZCPQ3MNyGQz+wfDhsHMn//v+e9i507nR3qpkiisVK+npAfXvvSubTIMG8Oij2qPsuee0I110tA7CXLvWBGEa3FMsp5RSlYHfgeeVUmOBiFLuMRiCE2d2mdBQXRpyyhTY7f+hYfY2mdTUFQ5tMiEhOi/ZsmU6T1l0tHaoe+01vdpZuNCsWoIVdxTLGCAMeACoDfQGbvKkUAZDwOLMLrNokd4y27MHEhL0N/LSpXqF44fG/rPl/PP1bt+AAUWDMJOS4IcfAmqxZigH3HE33mh5eRK43bPiGAwBTmmBmF266ESZH3+s/6W//XbtcmWNXKxgQZh//62nCDoP6C23GJflYMDpikUpda/l+IKj5j0RDYYAw2KXoaDAsV2mShUYOhSWL9dZI+2/iSEgjf3g2OCvFPTtC2+8oZMYtG4NV18Nn3xS1OPMGPwrFq62wqzpW046aQaD4VzZu9dx/65dJRWOn+PM4L96NXTrBvPmafPSoEHw4ovQpIlOPr1tmzH4VzScboWJyBzL8WnviWMwBBnOsi6fd57+937iRL3isUYq+jHuZF2uWlWH+owcCVu36liZiy+GjIxCg7/Juhz4uNoKc7gFZrbCDIZyxJmxf8ECHTPzxhs6+HLevJLLgQDnwgvh+ed14ssQyzdRdrbeOvvtN9/KZjg3XG2FWbe8GgBDgMqWdgNQw/OiGQxBgKuCZUlJ8P33Wrl88AG0aKFzlGVnVwhPMnBs8N+2TbstJybq6dpXwDQEBk4Vi4g8bdkGqwskish9InIf0BFo5C0BDYYKT2nG/osu0lGI770HX36p3apGjgz4tDHg2OAfEqJrxzz/PM4psMQAACAASURBVPzvf1rX3nKLtsFY3ZaNsd+/cSeOJVpEMqxvLK9jPSbRWaCUilVKbVFKvaaUesbX8hgMHqVrV/j8c5152Vk6/wDDlcG/Xz94911IS4P27XV+slatYNo0bXoyxn7/xR3FskUpNV8p1d3S5gJbyzqgUmqaUmqHUkqUUvF2/XFKqdVKqe2WYws3H5kJnAf8XVaZDIaAwj6Tsj3p6d6VoxxwJ+tyvXpw//06un/BAvi//9PxpgUFMH++jjk1+BfuKJbbgGPAq8As4Dgw6hzG/BjoBRR3hXkNmCUicZZx5lhPKKWaK6VSirWHgF0i0sUizwClVOw5yGUwBAbO0sY0qtg71EpBz55Qty6Ehem+3Fzt2/DII9o2Y/APlPgo14JSaicwQEQ2KqXqo0sg1xGRfKVUKJABtBCRQ24+bzbwXxHZ4ODcaHSWZqKiojouXbq0TDJnZWVRtWrVMt0b6ATz3MG/5l8/JYWW06YReuaMra8gNJS8iAg2/vvfnGjdulzH86e5Z2SEMWxYV3JyQm19lSsX0L//XlaurE/Dhqfp338/SUkHiYzML5cx/Wn+3sZ+7snJyWtFpJNbN4qIwwZcbzmOcdSc3eduA3YC8ZbXHYFNxc5vRjsNuHpGErAQmAf8x51xO3bsKGUlNTW1zPcGOsE8dxE/nP9bb4nExIgopY9vvSXy6aci9eqJzJlTrkP509zvvlskLKzoxllYmMiYMSI5OfojGDRIpEYNkZtuEvn+e5H8fH3v3r0ivXqJ7Nt3dmP60/y9jf3cgTXi5ve7q1xh8cD7QGdH+sgtreVhRGQFsMLHYhgM3mf4cMe5xFat0qHtv/4Kr7yiDf0VCFfp/CtX1t5kV10FBw9qJ7l774WsLO1V9uefhQb/WbN8I3+w4Mrd+EnLcaSDdi42FkfsBhpZtsCwHBta+g0Gg7vExcHPP8OxY9oft4JZtt0tsVy/Ptx3H/z+O7z/vvbifustbfCfN08rGYPncBV5f4WrVp5CiMhBYD2FJY+HAuvETfuKwWCwo1o1HfNy7bU6m/L//udriXyGUrrKZZUqhQb/vDxo00Yv+L7+OuBSsgUErrbCHnJxToAvyzKgUmomcC06oj9FKZUhIm2Au4DFSqkn0KWPby7L8w0GA/obdcIEHb4+ZIgu+Vi3rvN0/hWYfft0TjLrFpqIDsJs0waeeELHmg4bBjffrONlit97770JfPONyV92NrhKQpnsiQFFZBwwzkH/VqCrJ8Y0GIKWfv20YSIpCf75B/ItnlIBXPflbHEU3V9QoD+O//s/nQzzzTd1Ov8aNbSCGTYMGjbU927YUMPYZc6SUuNYlFI3KaVq2b2vrZSq2L+JBkNFomlT/W96fjH32wCN1j9bXBn8QSfDnDIFduyAmTNhyxa9mklK0gGYIorXXzfpY84GdwIkHxSRo9Y3InIEeNBzIhkMhnLHmRE/AKP1zxZ3Df4hIVqZLFigVzOVKxfaX3Jy4MknvS56wOKOYnFEaOmXGAwGv8FZtL6z/iDn+PGiSS/z83XS6Wee0Qs9g2vcUSz7lVLXWt8opa4DDnpOJIPBUO44qvtSqZLuN5TAkV2mcmWdo6x5c3jpJThp6ug6xR3FMh54Tin1p1LqT2AKcI9nxTIYDOVK8bovTZro9PsVLM6lvHBkl8nN1cb9r7/WoULNm8MLL+gATENRSlUsIrIFaA1cDQwE2oiISfdmMAQa9nVf0tP1t+fs2fDOO76WzO+wt8ukpq4oYpdp314HXaak6PfNm8Nzz0FmpqkTY8Udr7DWQEvLWwFaWvoMBkMg06gRfPEFjB8PK1b4WpqAIz4elizRH92mTVrBDBxo6sSA6wBJK1/YvQ4HotAp75t6RCKDweA94uNh6VK44QZITdV+toazolUrnS5m5UpITtYLwtdeg8aNdWxqs2a+ltD7uLMV1tSunQ/0Az7wvGgGg8Er9OkD//kPXHEF7N3ra2kCliVLINTiLxsSAosX6/oxLVrA2LG6+Kcjg39F3D47a3djEfkO6OkBWQwGg68YMQLuvFMrlxMnfC1NwFE8bUxenjZj/fYbfPCBXr385z86Lcwll8CLL8Iff2i7zbPPVrztM7dsLHYtXik1EqjhBdkMBoM3efRR6NoVBg/WLlAGt3HknpyfD5Mna2P/hAnw/fd6QTh+vPahGDRIO+bNnavvrUjR/e6sWL6wax8BNwC3eFIog8HgA5TSCbHCwnQeMR9Vlw1ESksbY6VaNV0vZtYs+OsvvXqxcvo09O2r7wn0j75U472IGCO9wRAsVKqkjflJSXrlsnYtvYMsG3JZKJ4exh327YNly4qmcNu2TRclE9G7kzfdpL3NAg1X9Vjso+3bFTs32pNCGQwGH1K1KowaBR99BLt2oUQKsyG//bavpaswONo+CwnRCamXLoUjR6B7d+0A8Npr+r0Vfzf4u9oKe8zu9aJi5+4qf1EMBoPf8MILJfdjgiQbsrdwtn22ejV06qQzLf/zjzZ9pabqJNXXXQcffwxPPeXfBn9XikU5ee3ovcFgqEg4y3ocBNmQvYU7WZcrV4YBA+Ddd/WisX9/+Pe/Cw3+8+frfn/DlWIRJ68dvTcYDBUJkw3Z76hZE26/XRcFrVxZ9+XmQlycNn19+KH/JMZ0pVjqKaXGKKXG2L+2vK/rJfkMBoMvMNmQ/RJrvIzVG9xaZrldO5gzR7svDxqkK2IePVr0Pm/aZFwplhSgs6XZv+4MfOd50QwGg8+wy4YsSukIv6pVoV49X0sW1Dgrs5yeDsuXF8bHfPihTmR92WVa4Tz6qHdtMq5q3o/0jghlRynVBRhleTsASBaRNB+KZDBUHIYPh+HD+d+KFSQlJekIvxEjdDh5gwa+li4oKS1epnZtuPlm3bKydIr/t96CTz7R519/HR5/3PM/vrJWkCwzSqlpSqkdSilRSsXb9ccppVYrpbZbji1Ke5aI/J+I3AU8DvxmlIrB4EH69NFuyDffXPLfZoNXcLfMMugF5uDB0LChjnkFHTPjjVWL1xUL8DHQC50h2Z7XgFkiEgfMAuZYTyilmiulUoq1h+zuvQuY62nBDYag56mntNvxiy/6WhKDGxTPYZaT453UMUp8lDtAKbUTGCAiG5VS9YHtQB0RyVdKhQIZQAsROVTKcyqjbT5JIuLw3yhLQOdogKioqI5Lly4tk8xZWVlUrVq1TPcGOsE8dwju+Ref+3kHDtDx7rvZOHkyJ1pX/NJMgfyzf/nlFnz55fnk5RWuISpVKuDKK/cxfnzpGzz2c09OTl4rIp3cGlhEfNKAnUC85XVHYFOx85uBRDeeMwK4191xO3bsKGUlNTW1zPcGOsE8d5Hgnr/DuS9bJhIbK3L0qNfl8TaB/LNPSHC0cab73cF+7sAacfN7ttRcYUqplugo/ObYGftFpItbmsvDiMhbvpbBYAg6Bg3StXlHj9bRe8rETPsjZclhVh64U0FyKfA+8DqQX8q1ZWU30EgpFSqFW2ENLf0Gg8EfeeklnWZ//ny44w5fS2PwI9xRLCEiMtWTQojIQaXUemAo8JbluE5Ksa8YDAYfEh6usyX26gU9epiyxgYb7niFrS6e3fhcUErNVErtARoDKUqpTZZTdwFjlVLbgbGYRJcGg//TqhU8/7wu7n7qlK+lMfgJ7iiWrsCvSqk/lFL/Z21lHVBExolIYxGpJCINRKSNpX+riHQVkTjLcVtZxzAYDF5k5EidU+S++3wticFPcGcrbLzHpTAYDIGLUrpgSGIivP8+XH+9ryUy+JhSVywi8j9HzRvCGQyGAKF6dViyBG67DRo10pkRY2NNYbAgxemKRSn1vIhMUEq9j4M0+SJyg0clMxgMgcX27ZCdDZmZ+r216iSYksZBhqutsFWW4+feEMRgMAQ4kyYV5nO3Yq06aRRLUOEqu/FnluNi74ljMBgCFmfVJXftgu++g6QkCA31qkgG3+CLJJQGg6Ei4qy6ZK1a8PDDukDII4/A5s1Fz7/9trbHGLtMhcEoFoPBUD44qjoZEQGvvAJr1+riICLQrx906gQzZ8J//6vtMLt26XNWu4xRLgFNqYpFKVXDG4IYDIYAx67qJErp49y5hfaV+HgdTJmeDs89B2vWwD33lAystNplDAGLS8WilFLAD16SxWAwBDrDh+v6uAUF+ujIaB8aqlctb7zh/Dnp6SUdAQwBg0vFYkmV/LdSqpaX5DEYDMGEM7tMpUpQrx4MHKi30rZu1VtlYGwyAYA7kfdZwDql1JeW1wCIyMMek8pgMAQHU6Zom4r9dlhEhN5C69dPe5OlpOiKlSLQtCn88kthSUQTK+OXuKNY/rQ0g8FgKF+symDSJL39FR2tlY21f+hQ3UQgLQ26dStUKlZMrIzfUapiEZGnvSGIwWAIUoYPL10pKAVxcXDsmOPzu3bBrFk6T1n9+uUvo+GscMcrLEIp9ZxS6ldLZuMpSqmI0u4zGAyGcseZTaZePfjpJ618Lr9cOwacOFF43thlvIo7cSyvoKs5jgfus7x+1ZNCGQwGg0Ocxcq8/LJWFv/8A7feCh98AE2a6BXM+PEmVsbLuKNYOovILSLyo4j8CIwCOnlYLoPBYChJabEykZFw443w6aewYwdceqlO6W9iZbyKO4pFKaUi7d5HAMpD8hgMBoNr3ImVAahdG+64o6Sx34qz3GaGc8YdxfIWujzxRKXUo8BPgIvIJoPBYPAjnNllAEaNgl9/9Z4sQYI7hb6eBx4GagN1gQki8qKnBTMYDIZywZldZtYsbey/4Qadu2z+fDh5svAai8G/d58+xuB/lrh0N1ZKhQI/i0hn4GvviOQ+Sqk44N/AfuCIiDzmY5EMBoO/UVqszEMPwfLlOiHmhAm6v3FjePppOHVK7/ubQMyzorSULvnAYaVUeHkNqJSappTaoZQSpVS8XX+cUmq1Umq75djCjcddCSwWkTFAVaVUx/KS02AwVCBc2WVCQ6F/f23wX7dOl1l+9FFj8D8H3LGxpAE/KKUeVkqNsbZzGPNjoBewq1j/a8AsEYkDZgFzrCeUUs2VUinF2kNoW09fpdRLQCMg9hzkMhgMwU50NEyeXJiXrDjG4O8WSpx9gNYLlHrdQbeIyKhzGlipncAAEdmolKoPbAfqiEi+ZQsuA2ghIofcfN7rwAsissXBudHAaICoqKiOS5cuLZPMWVlZVK1atUz3BjrBPHcI7vkH49y73Xgj4QcOlOjPqVGDn5Yt04GWQYD9zz45OXmtiLgXaiIiThsQCjzt6pqyNmAnEG953RHYVOz8ZiCxlGc0AhZY2n3ujNuxY0cpK6mpqWW+N9AJ5rmLBPf8g3Lub70lEhEhotcuup13nkhMjEiHDiKffy5SUOBrKT2O/c8eWCNufr+7NN6LXj30Oksl5zVE5B/gtnN9TkFBAYcPH+bYsWPk5+c7va5GjRps2VJiQRQUBPPcwffzDw8Pp3HjxlSuXNlnMgQVdgZ/SU9HWQ3+w4bBxx9rI//kybr16aODNQ023Mlu/IVS6kG0PcM+bf4p57ecNbuBRkqpUCncCmto6fc4e/bsQSlFbGwslStXRjn5JcnMzKRatWreEMnvCOa5g2/nLyJkZGSwZ88emjZt6hMZghJLcsz/rVhBUlJSYf+gQXD11fDee3D33dCokVYwO3c69zwLMtxRLC/YHQUddS/obbJyQUQOKqXWA0PRAZlDgXXipn3lXDl58iQtW7YkJEj2TQ2BhVKKOnXqcOiQV/4cDO4QGqrT+V9/vU54OXAgHD2qvc4g6N2T3QmQDLFrodZjWQdUSs1USu0BGgMpSqlNllN3AWOVUtuBsZb3XsMoFYM/42wVbfAxlSrp6P3IyEKlYiWI3ZOdrliUUtEi4tC3TimVKCK/lWVAERkHjHPQvxXoWpZnGgwGg0/Z7WTXPkjdk139m/6x9YVS6v+KnZvvGXEChwYNGqCUKtEaNGhQLs+PjY1l48aN5fIsKzt37qRu3boOz+3du5fk5ORyHc9gCBqc5SOrXdu7cvgJrhSL/dq7uCtK0K/LDzjwcXfV7+80bNiQ1NRUX4thMAQmjvKRhYdrW8z48ZCb6xu5fIQrxSJOXjt6b/AQL730Ep07d6ZDhw50796d9evX286tXr2aiy66iPbt29O+fXuWL18OwK+//kr37t1p164d3bt359di2VsffPBBunTpQtu2bVm5ciXgejVjMBhKwVGdmPnzYetWSEuDvn0hQP/pLBPOAlzQAYqtgNbFXrcGNrsbKONvzVGA5ObNm0v0oZWnR5o7xMTEyIYNG+TgwYO2vm+//Va6du0qIiIZGRkSFRUlP/74o4iI5OXlyZEjR+TMmTPSpEkT+fbbb0VEJCUlRZo0aSJnzpyRHTt2CCCLFy8WEZEVK1ZIo0aNJDs7W3bs2CF16tRxKs+JEyfckrui4g/zd/R76g2CMkDSjnOef36+yBNPiDRuLLJ6dbnI5C08ESAZAXxp997+dYVfsejPsSj2sQyuvHQc3VtW1q5dy9SpUzly5AghISFs374d0KuV1q1b06NHDwBCQ0OpVasWGzZsICwsjL59+wJwySWXEBYWxrZt26hWrRphYWGMGDECgN69e1OlShW2bdtG9erVy01mg8FgR0iIzpTcqZOOf5k8udAVuYLiVLGISKwX5TA4oKCggMGDB/PDDz+QmJjI3r17adSoEeBceYmIQ6XnTBE6u95gMJQzV10Fq1bpAMtff4VXXtF2mAqICd4oI1FRUWfVX1by8vJo0qQJALNnz7b19+jRg82bN7N69WoA8vPzOXr0KBdeeCFnzpyxGeJTU1PJzc0lLi4OgJycHN555x0AVq5cSXZ2Ni1btixXmQ0GgxPi4uDnn+HYMejdG2bO1EXEQkIqVDExdyLvDQ7Yv3+/R5+fl5dHREQEzzzzDJ07dyY6Opr+/fvbzteuXZtly5Zx//33c/LkSUJCQpg2bRp9+/blww8/ZNy4cZw8eZLIyEg++OADwsLCAKhTpw5paWl07dqVU6dOsWTJEts5g8HgBapV0+lghg6Fe+8t7K9A0fqlps2vaHTq1EnWrFlTpG/Lli20atWq1Hu9lS9q3759XHjhhezfv58qVap4fDx3MLnCfD9/d39Py5sVxXNlBRkem39srFYmxYmJ0XnH/AD7uSul3E6bb7bC/IyZM2eSnJzMtGnT/EapGAwGD+AsKr8CROubrTA/Y9y4cYwbVyLjjcFgqGhERztesTiL4g8gzIrFYDAYfIGjaH3QtpcAxygWg8Fg8AWOovUff1xH7Fu8PQMVsxVmMBgMvsJSTKwI3bvDNdfAt99Cu3a+kescMSsWg8Fg8Cf699fBk/376zxjAYhRLH5Kbm4uTzzxBHFxcbRp04ZWrVrxwAMPkOuBLKnTp0/n4MGD5/ycRYsW2VLOOCI1NZWuXbuSkJBAq1at6NOnDwWW4kgFBQVMnTqVuLg4WrVqRevWrZk3b57t3p07d1KpUiUSEhJo164dnTp14qeffnI6VkFBARdddBF79uw5pzk99dRT5OTk2N4/8cQTvPvuu+f0zOIopcjK0lW/k5KS2LFjR7k+3xCA3HCDTgPTr5/zWi/+jLtJxSpKczcJpSMcJSLcu1ekVy+RffvceoTbDB8+XK699lrbmDk5OTJnzhzJzMws34GkMOGlK9xJwti7d2/57LPPHJ7Lzc2VWrVqye+//27r++2336SgoEBERJ555hnp1auXZGRkiIhIenq6tG7dWt58800RkRJJMmfPni0XXnihU1mWLl0qo0ePLlXm0gAkMzPTo0korWOIiHz00Udy8803O7zOJKH0DT6d/7RpIi1bitglo/UmZU1C6fMvem+38lYsd98tEhIiMmaMW49wi+3bt0tERIQcOXLE4fm8vDx54IEHpE2bNtKmTRt54IEHJC8vT0REbrnlFrnzzjslOTlZLrjgArnppptsX95z5syRCy+8UNq3by9t27aVLVu2yOTJk6Vy5crSsmVLad++vWzatElSUlKkW7dukpCQIPHx8bJkyRLb3Hv37i0PPvig9OzZU5o2bSoTJkwQEZGFCxdKZGSkNG3aVNq3b2/LrmzlyJEjEhYWJvv37y8xn+zsbKlSpYps2bKlSP+XX34pzZo1E5GSimXz5s0SERHh9DPs06ePrFixwvY+LS1N+vTpI23btpUOHTrIV199ZTsHyJNPPindu3eXuLg4+eCDD0REZMyYMQJI27ZtpW3btnL06FG55ZZb5JVXXhERkSeffFKGDBki/fv3l+bNm8sNN9wgv/32myQnJ0uzZs3kwQcftI0xbdo06dSpkyQkJEi3bt1k3bp1Rca3KpacnBypV6+ew981o1h8g8/n/9hjIh06iBw75vWhjWLxgWLZu1ckPFx/ilWqlN+q5d1335X27ds7PT979my55JJL5MyZM3LmzBnp06ePzJ49W0S0YunZs6ecPn1azpw5I61bt5bly5eLiEj16tUlPT1dRPSX+cmTJ0Wk5IrlyJEjNkW1f/9+adSokezatUtEtGK54YYbJD8/X44dOyZ16tSR7du32845W7GIiNxzzz1SrVo1GTBggDz33HM2WTZs2CA1atQocf2RI0cEkBMnTpRQLFOmTJFRo0Y5HCcnJ0fCw8Pl1KlTtr4uXbrI/PnzRURk06ZNUqdOHVtJAkCefvppERHZunWr1K5dWw4cOGA7Z79iKa5YLrjgAjl27Jjk5eVJu3bt5NJLL5Xs7GzJysqSevXq2T4bZ+UP7MewkpycXETxWTGKxTf4fP4FBSL/+pfIRReJWP5mvUVZFUvA2FiUUjWUUguVUrtd9ZXfeCVb9erVirxv2BCys/X1p0/D+ec7vq94Kw39M3ROSkoKt956K2FhYYSFhTFy5EhSUlJs56+55hrCw8MJCwsjMTGRv/76C4A+ffowcuRIXnnlFf755x8iHPnQA4cOHWLw4MHEx8dz2WWXceTIEdLsjIjXX389ISEh1KhRg1atWtmeXxqvvvoq69evZ+DAgfz666/Ex8eTlpbmdL7WrMvW47Fjx0hISKBRo0bMmDGDxx57zOF9hw8fJiwszJa5IDMzk/Xr1zNy5EgAWrduTUJCAj///LPtnttuuw2Ali1bkpiYWOScKy677DJq1KhBaGgo7dq1o1+/fpx33nlERkbSsmVL22ezdu1aevXqRXx8PPfff3+Rgm3FadCgwTnbhgwVCKVgxgxo2hSuuw4WL/b7xJUeVSxKqWlKqR1KKVFKxdv1xymlViultluOLUp7logcF5FRwDZXfeWFXocUbSdOZNpe791bMuN1lSqwb5/je+1baSQmJpKWlsbRo0edyFYy1b39+3A7wUJDQ8nLywNg2bJlTJ06lZMnT5KcnMxXX33l8Pl33303SUlJbNiwgfXr19O4cWPOnDlT6vPdoVmzZtx+++18+OGH9OjRg88++4wWLVqQk5PD1q1bi1y7evVqmjZtStWqVQGoWbMm69evJz09nSFDhnDjjTc6HKNKlSpkWzU+zhV1eZQSKP5ZOPpscnJyGDx4MNOnT2fjxo18/fXXRT7P4mRnZ5t0PoaihITAwoVw6BDcdpuO2BcpTFzpZ8rF0yuWj4FeQPG8Ba8Bs0QkDpgFzLGeUEo1V0qlFGsPeVjOs+bZZ8Hi0GQjP1/3nystWrTg6quv5s477yQzM9Py7HxmzJhBVlYW/fr1Y9GiReTm5pKbm8vixYtthb2ckZeXx99//02XLl145JFHuPTSS1m3bh0A1atX5/jx47Zrjx07RmxsLEopvv32W/7880+35C7+HHuysrJYvny57Uv+2LFj7Nixg6ZNmxIeHs6ECRO46667bMp09+7dPPTQQzz11FMlnhUaGsoLL7zAvn37+OSTT0qcr1mzJlFRUey0JPKrXr06CQkJLF68GICtW7fy+++/07VrV9s9r7/+OgBpaWmsX7/edq5atWpO5+Qu2dnZTssfOGLLli20b9/+nMY0VEAqVYKDB/UXjT2nTsGkSb6RyQkeDZAUkVVQ9D9DpVR9IBHoZ+laAryqlKonIodE5C/A9bfkWaKUGg2MBl0vZcWKFUXO16hRw/YF7or8/HzbdatWRZCTE1rkfE4OrFyZT2bmqXOW+dVXX+Xf//43HTp0ICwsjIKCAi699FJycnIYOnQomzdvtn35XHLJJdx4441kZmaSm5tLdna2TU7r+2PHjnHTTTdx/PhxQkJCaNSoEZMmTSIzM5PRo0dzyy23EBERwYIFC3jiiSe4//77mTJlCvHx8cTHx1NQUEBmZib5+fmcOnWqiMKzvh8xYgSPPfYYL7zwAs8++yzJycm2+WRmZjJjxgzuuecewsPDycvLY/DgwfTt25fMzEzGjx9PXl4enTp1olIl/Wt59913M2jQIDIzM8nKykJEivycJk2axJNPPkmfPn1KfH5XXnkln3zyCaNGjQJgzpw5jB8/npdeeolKlSoxZ84cwsPDizyvW7duZGRkMH36dKpUqUJmZib/+te/SEpKIjw8nC+//LLI53vmzBlycnJKfNbFPxulFJMmTaJTp040btyYfv362T4T+89HREhPT6egoICYmJgSv5PZ2dklfne9QVZWlk/G9Rf8af699+zB0Vpa0tP5nwdkLPPc3TXGnEsDdgLxltcdgU3Fzm8GEt14zixgD3rF09RZn6tW3l5hwUKgzf3vv/+Wrl272jziXEEx47kjvDX/Rx55xOZkUBxjvPcNfjX/mBjHO+wxMR4ZrsIb7wFE5B4RaSwid4nIDmd9BkPTpk154IEH2Ldvn69FOSsaNmxoczIwGErgKHFllSq634/wRa6w3UAjpVSoiOQrpUKBhpZ+g6HcuP766926TtzxqPASY8eO9bUIBn/Gmlds0iRdtyUyEhIS/K7ipNdXLCJyEFgPWHNDDwXWicghb8tiMBgMAcfw4brCZEEB7N8Phw8Hl1eYUmqmUmoP0BhIUUptzqAeogAADZNJREFUspy6CxirlNoOjLW8NxgMBsPZEBkJS5bA+PHgZjyZN/C0V9g4oEQ5RBHZCnQteYfBYDAYzoqEBHjsMRg2DFatgsqVfS1RYBnvDQaDweCAceOgbl144glfSwIYxWIwGAyBj1Lw+uvwxhvw3Xe+lsYolnPi7bc9lrPHVT2WRYsWMXjw4LN+5lNPPcWDDz5YbjLGxsaycePGcnueI6699lp+/fVXpkyZQkJCAgkJCVStWpVmzZrZ3u92s17FwYMHufLKK2nXrh0XXnghY8aMsaWjefjhh1myZIknp2IweJb69WHRIrjlFp36xYcYxVJW3n5b5+jxUM6ekSNHsmnTJtauXcumTZv4448/aNmypcscUxWNX375hZMnT9K5c2cmTZrE+vXrWb9+PZ06dWLmzJm299ZUKaUxdepUWrVqxR9//MGGDRtYu3Yty5YtA7CljykonqfHYAgk+vXTtpZRo9xLTOghjGIpK5Mm6Rw99pRTzp60tDQ++ugj5s+fT7Vq1QCoXLkyo0ePtiVkPHHiBEOGDKFNmzb07NmT/fv32+5/4YUX6NKlC4mJiVx11VVFzqWnp3PFFVfQpk0brrvuOo4fP87p06c5//zziwQTjhs3jqlTpwI6GeSll15K+/btad++PcuXL7dd995779G9e3diY2N59dVXbf3btm2jf//+dO7cmfbt29tycT377LPcd999tusyMjKoU6cOJ0+eLPE5zJ07l2HDhp3TZ2mPUorMzEwKCgps6VgaNWoEQL169WjWrBnf+cE2gsFwTkyerN2Q7f4evY67IfoVpbmd0qX0JMVlb6VQWj2W119/XWrWrGmrZ3L77bfLxIkTRUTkzTfflDvuuEPy8/NFRNduGTZsmIjo+iENGjSwFdsaOXKkPPDAAyIiMmHCBHnqqadERGy1RA4cOCAZGRkSFRVlK9yVl5dnK0AWExNju3/Hjh0SGRkpmZmZkpubK4mJibbCXSdOnJC4uDjZsmWL7Nq1Sxo0aCC5ubkiIjJz5kwZOXKkw3k2a9asRPEvEcd1X8aOHSvt27d32P78808REcnIyJCkpCSpX7++REZG2oqUWXn66adL9Fnxh5Q2JqWLbwjI+aelidStK2JXsbUslDWliy8i7wMDB8vIzMxM2wqC2Fi9/VWcmBgdvHROQ5e+hO3Zs6dtC6hbt258++23AHz66aesWbOGxMREQGc1rlGjhu2+AQMGEBUVBegaJNZI73vuuYeLLrqISZMm8eabb3LppZdSv359vvjiC1q3bm3L9hsaGkqtWrVsz7Omro+NjaVWrVrs2bOHgoICtmzZUiSt/ZkzZ9iyZQuDBg2idevWfPnll1x99dUsWrSI6dOnO5zjnj17bLKWxsyZM0u95v3336ddu3Z89913ZGZm0r9/fz744AObvapBgwb88MMPbo1nMPg1F1wAL70EN94Ia9aUTAPjYYxiKStTpmibiv12WEREueTssa/HYv8lbo+zmigiwmOPPWbL6usKkcK6I02aNKFz58588sknzJ49mzlz5tiucYUjOZRS1K1b12kxq1tvvZXFixfTrFkzjh8/zsUXX+zwuuJ1VVwxbtw4p0rhww8/pHnz5rzyyissXLjQVqRs4MCBpKam2hSLqYNiqFDcdBMsXw733Qdz5pR+fTlibCxlZfhwmDtXr1CU0se5c8slZ09p9VhccfXVVzN79mxbXZMzZ87w+++/285/8cUXHLJ4jCxatKhIavuxY8cyfvx4KlWqRPfu3QHo0aMHmzdv5pdffrHJ4awAmZWWLVsSERHBm2++aevbunUrJ06cAOC6667jhx9+YNq0adx6661On9O2bVu2bXOvhpu9Mb94a968OaATU3799dcA5OTkkJKSQny8rf6cqYNiqFgoBbNnw0cfQb16Xq04aRTLuWCfs2fnznJNBLd48WJatGhBx44diY+Pp23btuzevZvzzjvP5X033XQTw4cPp3fv3rRr146OHTvy448/2s5fcskljBo1ijZt2nDkyBEef/xx27nevXsTHh7OmDFjbH21a9dm2bJlTJw40fa8tWvXupShUqVKfPbZZyxdupR27drRpk0bxowZQ05ODgAREREMHDiQN998k5tvvtnpc6699lq++eYbl2OdDdOnT2flypW0bduWhIQE4uLiuOOOOwC9Mvvuu+8YOHBguY1nMPiczz6DzEydT8wD3qtOcdcYU1GaqcfinL///lvOP/98OXnyZIlzvpj78ePHpW3btnLq1CmPj/X111/LiBEjnJ73h5+9Md77hoCe/znWbwmKeiwGz/HEE09w8cUX89JLLxHhZUOfM6pXr85LL73Ejh2eL7Nz4sQJnn/+eY+PYzB4lfT0s+svJ4zx3gDAM888wzPPPONrMUpgLePradyt3WIwBBTR0Y69V6OjPTqsWbFYMBHXBn9G3HBBNxhK4KjiZDl5r7rCKBYgMjKSf/75h5ycHPMHbPA7RISMjIwirt0Gg1t40HvVFWYrDGjcuDGHDx9m165dtngQR2RnZwftH3cwzx18P//w8HAaN27ss/ENAczw4V4vXWwUCxASEkL9+vWpX7++y+tWrFhBhw4dvCSVfxHMcwczf4PhbDBbYQaDwWAoV4xiMfx/e/cWakUVx3H8+9NIqMhTWdlFjEANEoS0HiyjIIhCxaIowewEEV196UJRdO8p6Z6gmNjFLphdFHvIMqWrXTUzVAgPRRllSSF0gfr3MOvotNmns/dp3INn/T4wnDVrZvZa/7OW5+/M3nvGzKxSTixmZlYpJxYzM6uUcvt4raQfgSbfGGrJCGBHhd3Zl+QcO+Qdf86xQ97xl2MfHRGHt3JQdonl/5D0cURMqrsfdcg5dsg7/pxjh7zjH2jsvhRmZmaVcmIxM7NKObG0Z0HdHahRzrFD3vHnHDvkHf+AYvd7LGZmVimfsZiZWaWcWMzMrFJOLC2QNFbS+5K2pp9j6u5TJ0nqkbRZ0vq0nF13n/YWSXMlbZMUksaX6rOYA/8R/6CfA5IOk/SapC2SPpf0kqTD07ZBPf79xN7+2Lf6DOOcF2A1MCuVZwGr6+5Th+PvAcbX3Y8OxXoaMKox5lzmwH/EP+jnAHAocEZp/X7giRzGv5/Y2x57n7H0Q9IRwEnAc6nqOeCk3mxug0tEvBMR35TrcpoDzeLPRUT8HBFrSlUfAKNzGP++Yh/o6zmx9G8U8G1E/AWQfn6X6nOyJJ0iz5PUVXdnOsxzoJDNHJA0BLgKWE5m498Qe6+2xt6JxVoxJSImACcDAh6ruT/WebnNgUeBXQz+OJtpjL3tsXdi6d83wDGShgKkn0en+iz0XhqJiD+AecCp9fao4zwHMpoDkuYCY4CLIuJvMhr/JrEPaOydWPoRET8A64GZqWom8FlE/FhfrzpH0oGShqeygIspfh/Z8BzIZw5Iug+YCMxIf0izGf9msQ907P3N+xZIOgF4EjgE2AnMjogt9faqMyQdDywDhqblS2BORGyvtWN7iaRHgPOBkRS3C/8pIk7MZQ40ix+YRgZzQNKJwBfAVuC3VL0tIs4b7OPfV+zA9Qxg7J1YzMysUr4UZmZmlXJiMTOzSjmxmJlZpZxYzMysUk4sZmZWKScWy1a6a+v4VO6WNHYvtNEl6aaGuoWSplTdVpO2Z0g6ZW+3Y9bIicWs0A20nVgkDUlfHOtLF/CvxBIRl0fE2+22NQAzACcW6zgnFsuepMuAScAj6XkTZ6X6myR9KOlTSSskjUz1d0p6RtIrwAagKz3H5CNJGyS9Kan3zrCPp+3rJb2Xjl8jaWoqHynp5XSDv42SZpf61SPp7vT8jx5J1/bR/8mpj+slbZI0Mz0zYzpwc6qfnfa9VNI6SZ9IWi1pXKrvlrRK0rIUw2pJx1T/27Ys1P0cAC9e6looPWcCWANMLW2bBSwAhqT1q4AlqXwn8DUworR/uXw58HwqHwfsaGh3d1vAC8A9qXwUsL3Upx5gbul1dgEHNYnjVeCSVBbQlcqLgWtL+00BVgLD0vo5wLup3E3xjetxaf0O4MW6x8jLvrnsN9CEZDbITac4i/k0XenaD/iltP21iNhRWj9H0jXAQWnfVp1FcdsMImK7pJXAmRS31wB4Pm3rkbQTOBbY3PAabwG3pLOkVRGxro+2pgETgHUpJlHcoqTXO7HnNiULgY1txGG2mxOLWXMC7o2IRX1s37V7x+IP+oPAyRGxTdJk4Nk22mq8r1J5/fdS+S+a/JuNiIckraBIUo9Kej0ibmvSjoBFEXF7C31Sk36ZtcTvsZgVfgWGl9aXA1dLOgRA0jBJE/o49mDgT+D79JCkKxte9wBJff0n7g3gitTGSOBcijOQlkkaGxFfRcR84GH2vGHfGNMKYLakY9NxQyVNLG0/VXue5d7dbj/MevmMxaywAJgr6Qbgxoh4WtIIYG26bDSE4lkUGxoPjIiNkpYCmyjee1kLnJ62/SxpCbBR0s6ImNxw+BxgvqTPKc4Sbo6ITW32fY6kMymS2x/Adan+aWCxpAuBByLiKUm3AsvTM0X2B5YCn6T91wJ3pTvd/gRc0mY/zADf3djMKD4VRvGBggvq7ovt+3wpzMzMKuUzFjMzq5TPWMzMrFJOLGZmViknFjMzq5QTi5mZVcqJxczMKvUPL71JOBL0tMgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.rcParams[\"font.size\"] = 11\n",
    "plt.grid()\n",
    "plt.yscale('log')\n",
    "plt.xlabel('Iteration step')\n",
    "plt.ylabel('Error in Euclidian norm')\n",
    "plot_pure_jacobi(P, q, 'black') \n",
    "plot_sor_jacobi(P, q, 'blue') # maybe this is optimal for original SOR\n",
    "plot_chebychev_sor_jacobi(8, P, q, 'Chebychev (T=8)', 'red')\n",
    "plt.legend()\n",
    "plt.savefig('jacobi.pdf')"
   ]
  }
 ],
 "metadata": {
  "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}