Instantly share code, notes, and snippets.

# tobydriscoll/TB_Lecture_37_julia.ipynb

Created February 3, 2017 19:03
Show Gist options
• Save tobydriscoll/2511c290c58f98c6989672082897d47e to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
 { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Lecture 37: Conjugate gradients\n", "\n", "When $A$ is hermitian and positive definite, a famous method is available: conjugate gradients. There is a lot to say about CG overall, but to us it is a Krylov subspace method that minimizes a quantity other than the residual over $\\mathcal{K}_n$. This quantity is $\\|\\epsilon_n\\|_A$, where $\\epsilon_n=A^{-1}b-x_n$ is the error and the norm is defined by $\\|u\\|_A^2=u^*Au$. \n", "\n", "We won't use the iteration formulas here, just MATLAB's implementation of pcg. " ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAABcSAAAXEgFnn9JSAAAA\nB3RJTUUH4AsHFQ0SbG3WFAAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJ\nbmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAwNy1Ob3YtMjAxNiAxNjoxMzoxOHLK3lwAACAA\nSURBVHic7d17WFTVwj/wBcNVBhWMbAAFQgdBsUJ537ykIsZgXjp4KTpigpgGlpfS6heWlzJL0/fV\n85jlgym8pw6ZetJAC1FR8QIYlub1gIwKcpFQIC4yzMzvj9XZZzc3GZiZfZnv569h7b1nryXjfNlr\nrb22g1arJQAAAFxz5LoCAAAAhCCQAACAJxBIAADACwgkAADgBQQSAADwAgIJAAB4AYEEAAC84MR1\nBUB4qqqqCgoKSkpKSkpKtFptSEhIRETErFmzXFxcuK4aWMW5c+euXbvW3Nw8ePDgUaNGmd65Cx8P\nnUMGDRo0YsSIF154QSKRWLopwG9aAHPs2LHDw8ND/4MUEBDwj3/8g+vagYXdvXt33LhxzG/5xRdf\nNL1/Fz4exg4JCQnJycmxQpuAvxBI0FkPHjx48cUXTfxx4+DggEwSmTlz5rB/xSYCqQsfj4ce4uzs\nnJeXZ/1WAl84aLF0EHTOxo0bly1bRl/37dt31qxZI0eOdHR0PHjw4M6dO9VqNSHEzc2tuLh4yJAh\nnNYULKZPnz719fWEkNTU1FWrVrm4uPTq1cvgnl34eBg85N69e999911ubq5KpSKE9OrV68qVKzKZ\nzNotBV7gOhFBGGpqanr27Ek/M48//nhZWRl764EDB5ju/rfeesvgOyiVyvLyco1G08nTlZaWqtXq\nbla7o6Pj2rVrjY2Nndm5oqLi2rVr+oWVlZUWP1cnG3jjxo3S0lKVSmWp8+oz8Xtpbm5mviiOHz9u\n4k268PEwfci6deuYU2/durVrTQPBQSBBpzB/yRJCDh8+rL9DXFwc3RoSEsIuv3r16sSJE5k/qz09\nPZ999tmLFy+y95k+fbqfn5+fn9933333448/Dhw4kO7s7u6elJT0+++/093i4+PpbvPmzWMfvnbt\nWlqelJTEFJ49e/aZZ57p0aMHIcTR0TE8PPzDDz9kf+2yT1pSUhIZGUkIkclkdKtGo/nggw/8/f1p\nTfz9/Xfv3r1jxw56yPz589kVMOtcJhrIuHLlikKhYP7RJBLJ888/X1xcrLPbQ89rwkN/L35+fr6+\nvswv3cfHx8/P729/+5vBd+vCx8P0IQ8ePAgNDXVycnJycnruuec60yIQAQQSdEpUVJTBvGE0NTXd\n+TemMCMjw83NTf+63NnZedu2bcxuzzzzDC1/4YUX9OdiTZs2je62bds25tvzwYMHzOGDBw+m5bt2\n7aIlH3/8scEJWjNmzGhpadE56Zo1azw9PelrJpBeeOEFnWOdnJzGjx9PX0+fPp05u1nnMt1A6quv\nvqIZo1+BQ4cOmXVeYzrze9HfSgj56KOPDL5hFz4eDz0E7BACCTqlb9++9Otj1qxZnTzk5s2bzOyp\n8ePH79q1KzMzU6FQ0BJXV9erV6/SPZnva0KIp6dnamrqokWLpFIpLXF0dLx165ZWq62trXVy+uNG\nhR9++IEeW1paSktcXFzu37+v1WpPnTrl4OBAT5GWlvbPf/5z6dKlPj4+dLc1a9bonJSpZEBAQGRk\npFar/f7775n6zJw5c8+ePdu2bQsNDWUKmUAy91ymG6jVau/cucOUjx07NjMzc8uWLcOGDaMlUqm0\noaGh8+ftzu8lOzt7z549TLU//fTT7Ozs0tJSS308unAIiB4CCR7u3r17zBfTO++808mjEhIS6CEj\nR45sb2+nhR0dHc8++ywtnzp1Ki1kvq/79OlDQ0Wr1R47dow5KXNlEBMTQ0tSUlJoycaNG2nJlClT\naMlTTz1FS7766iumMvv376eFPj4+9AKCHRJTpky5e/cuszOziamhVqutq6vr3bu3TiCZe66HNnDe\nvHm05Kmnnmpra6OFVVVV7u7utPyf//xn58/bzd/L77//ztTwzJkzxt6wCx8PE4fExcUF6KmoqOjM\n24LQYaUGeDg3Nzemd6itra2TR505c4a+SE1NdXZ2pq8lEslrr72mswPjueeeY0Y1IiMjmZPW1NTQ\nF/Hx8fTFgQMHtFotIYT5FqYTiBsaGs6fP08IcXR07Nev39l/69Onj5eXFyHk7t27R44cYZ/0kUce\n+eabbx555BGm5Oeff6YvFixYwBT26dNn+vTp7AO7cK6HNvD48eP0xcKFC11dXenrxx577NixY9nZ\n2dnZ2YMHD+7Cedm68HsxrQsfDxOHVFdX39TT0dFhVpVAoLBSAzycm5tbUFAQ7RxTKpUG9ykvLy8p\nKaGvJ06cSEvoj+zOLkJIWFgYfXH37t27d+8yHU2EEPbsXg8PD39//5s3bxJC2tvbaWFcXNyrr77a\n3t5eWVn5008/BQUFnTp1itZw6tSphJBff/2V7qnRaMaMGWOwqlVVVewfJ06cyFx/0Fo1NTXR1wEB\nAew9g4OD2T924VymG9je3n7jxg269cknn2Qf+N///d/Ma9pks87LaGlp6cLvxbQufDx69Ohh7JCB\nAwe2tLQQQtra2q5du9bJOoA4IJCgUwYPHky/Po4cOdLY2MhM2GUsW7Zs3759hJCePXvW19c3Nzdr\nNBq6SWfsnf0jvdeEwVwTUI6OulfwvXv3VigUdIxn//79AwYMoDe4TJw4kU5MePDgAXPsgAEDDLaF\nqRjl5+fH/rFnz54ODn/cn9fQ0MDexJ4G3bVzmW5gc3MzbQ4hhLl20deF8zI6Ojq68Ht5KHM/HiYO\nycjIoC8OHDjw/PPPm1UNEDoEEnTK5MmTaedYU1PTli1bVqxYwd5669YtputszJgxEomkZ8+e/fr1\nu337NiGktLT0iSeeYHZmpiH07t2bPbG4k+Lj49mBRAuZG/6ZGXdarfann35iJgiYoPPV7+rq+thj\nj9ErjLy8vJEjRzKbDh8+zN6zC+cyzcvLSyaT0VNfv3596NChzKaMjIwLFy4QQiZMmBAREdHl81rp\n92Lux+OhhxBCsrKyzKoDiADGkKBT5s6dy3wPvv/+++vXr2e69e/fv//6668zf9ozAy3M3LAdO3aw\n3yo9PV1nB7NMnTqV9rBdvHgxJyeHENKjR4/JkyfTrX379qU3D2m12gMHDjBHaTSaxYsXT5s2bdq0\nadXV1Q89BX2xcePGgoICQoharV65cmVRURF7N4ucS8fw4cPpiy+//JIpvHfv3sKFCzdt2rRp0yat\nVtvN81rj99KFj4eJQzo6OtavX/+Pf/zD3GqA4HE2ncJ8RUVF8+fPj46Onjlz5ubNm5k5SGAbBQUF\ndKox1aNHj3Hjxo0fP56Ze0YImTRpErP/r7/+ylx8JCQkHDly5NixY8nJybTE0dGxqKiI7slMQnvv\nvffYZwwKCqLl27dvZ5fPmDGD/RmeOXMme+vf//53Wt6zZ88FCxbk5eXt2bOHyZioqCjTJ9VqtWVl\nZcz8ckKIv78/nSfN3ELEzLLr5rn0G/jzzz8zp05MTDxx4sSuXbtGjBhBSx577LHW1tbOn9egzv9e\nOjnLjjL342HskJiYGP3rM6VSafrsIA6CCaQTJ06EhIQkJCT83//93/r168PDw+fOndvJm9LBUo4c\nOaIz4sIWHh5eW1vL3n/Dhg3640DUypUrmd3MDST2/TGEkG+//Vannvq3tVKPPvro5cuXTZ+UaWmf\nPn3Yx06YMGHt2rX0NfvG2O6cy2ADP/74Y/bXNMPZ2bmgoMCs8xrTyd+LWYGkNf/jYfoQuVyelJRE\nXyOQ7IRgAmny5MlxcXHM2l8HDhyQy+Xs/59gG7/99ltiYqLO37De3t6bNm0yuORaQUHBk08+yXzD\nOjg4DBky5OjRo+x9zA2klpYWZuBEKpUavOdm586d7G86R0fHuLg45lZcEydlVFVV/f3vf1+yZMlb\nb721b9++Bw8erFmzhh6SkJBgkXMZa+Dx48eHDh3KjqXo6Ohz586Z20YTOvN7MTeQtOZ/PAwe4ujo\nOHXq1Pv37//4448IJLsijNW+6+rqRo0a9c477zB/MalUqoiIiISEhLfffpvbutmt2tra8+fPOzg4\nhIWFMWu+GdPc3Hz58mWNRjN48ODuD/53Xl1d3aVLlzw8PB5//HFvb+9OHnX69OnKykpCSGBgIF3j\njhozZszJkycJIe+++y5ztdTNc5nQ1NR08eJFZ2fngQMHsju+LHhe6/1ezPp4sA9xdXUdNmwYs5gT\n2BVhBNL58+fj4+M///xzZv0rQkhsbKxcLt+yZQuHFQNR+uSTT9555x1CiFQqLSkpoWuhFhcXjxw5\nsqOjw8HB4ZdffgkPD+e6mgBiI4xZdvRGOZ3HSvbo0UPnvhAAi5g1axa9g/X3338PDQ0dNmxYZGTk\nqFGj6DSwefPmIY0ArEEYgUQv43Qu5mifI0c1AjHz9/c/ceLE6NGjCSFqtbqkpOTcuXMqlapHjx7b\ntm3bvn071xUEECeb3hjb0NAQFxc3d+5cZnlHtry8PLqcsFQqlcvlycnJzMIt9L6T1tZW9v6tra14\njiRYyYABA06ePHn58uWSkpKqqiqpVBoQEDBu3DiDD4YAAIuwaSClp6dXVlYa7Gdbu3ZtZmamq6tr\nWFhYY2Pj7t279+/fv3XrVvpXKg0e9ppXarW6oqJi1KhRtqo72KOwsDBmhTcAsDard9mp1WqlUnn8\n+PGlS5ca6+soLCzMzMzs37//oUOHsrKyDh48uHnzZpVKlZaWRlcC9vX1lclk7IVb8vPzVSoVc1u7\nWXYVG153EgAAOGT1K6QbN24wy7oYs3fvXkLIsmXLmJsqFApFbGxsTk7OqVOnoqOjCSEvv/zyJ598\nsnbt2pkzZ16/fn3jxo39+vVjT7rrjKjPSvLL7hNCVueWH0uJCPQ28NBMAADghNUDSSaTMTOzL1y4\nwKyXxVZYWCiRSMaOHcsujI6OzsnJOXv2LA2kOXPmNDQ07Ny5MzMzkxASHh6+bt06g49hZoSEhOiU\n3B71FvEKIoQo69vkad/3rDjb53p2NxoHACA8vH2uh9UDSSqVMk9HZq8PxlCpVLW1tf7+/jrpQp89\nU1FRQX+USCRLly5dtGhRZWWll5dXJ++b0/l3T8q6wvTXqXr0+U0+6fXXXlupCDKzTdwLCQnh7Ueq\nO8TaLiLepom1XUS8TdP/S50/uJ/23djYqNFomMdoMmiJzgNpJBJJ//79u3wX98qYIJ1uulW55UFr\nTyvrO/sUVAAAsBLuA4k+Ckz/4omWmPugMNMCvd2OpUSsivnTJZGyvi1qW8nqH8steCIAADAX94FE\nn9alHzy0ROehlt0X6O22UhGkn0mrcsuRSQAAHOI+kAze9MqU0K0Wt1IRVJ420mD3nSBiSZRd20S8\n7SLibZpY20VE3TTe4j6QpFKpp6fn7du3mWdKUvQ2WOutxWCs+25VbnnUZ+cxqgQAYGPcBxIhJDQ0\ntL29/dKlS+zCkpISQohV75M32H1HCMkvu4dRJQAAG+NFIMXExBBCNmzYwJRUVVV9/fXXEolk/Pjx\n1j477b4zeKmECXgAADbDi0CKj48fMGBAUVHRtGnTduzY8dFHH82YMaOlpWXu3Lk6j560EnqptDM+\nVGdUCRPwAABsxqaLqxrj7OyckZGxZs2avLw82nHn4eHx5ptvJicn27IaiZGyccFeGcVVq3L/k0D0\nUim/7L5+XAEAgAXx64mxKpVKqVR6eHjIZDIHB4duvluXb7SmF0Y6nXWB3m4rY4ISI/HACwAQMD6v\nQMGLLjuGs7PzwIEDfX19u59G3WFsAl5S1hV03wEAWAm/Aok/6KjSsVTdFcEx0wEAwEoQSKaMC+6N\npYYAAGwDgfQQWGoIAMA2EEidYmKpIXTfAQBYBAKps7BSOACAVSGQzIDuOwAA60EgmQ3ddwAA1oBA\n6goT3XfIJACAruHF0kFCRLvvCCE66wxFbSvZGR82Lrg3d1UDsJGQkBCuq2CneLvUQjchkLplpSJo\nTqSMfWGkrG9LyrqcOFxG4wpA3MT6zchnIv47AF123UW778YFezElmOYAANAFCCQLCPR22xkfqjOk\nhEwCADALAskyAr3d5kTK9DMp6rPzXFUJAEBYEEgWY/Aupfyye5gODgDQGQgkC6NPnmWXYDUHAIDO\nQCBZXmKkTOfOWUxzAAB4KASSVdCpd/qrOWBICUBwVq1aNXnyZP3XXL2PiCGQrMXgag4YUgIQnNbW\n1qamJv3XXL2PiCGQrMjYYqwYUgIA0IdAsjr9xViV9W27zlUhkwAA2BBItqA/pIRpDgAio9Vqua6C\n4CGQbMTgkBIeWgEgXM8///zSpUtbWlpSUlJCQkL27dvHdY0ED4FkOyaGlJBJAIJz/fr1mzdvJiYm\nZmRk9OrVy9fXl+saCR4CydZWKoKOpUawSzDNAUCgDh06dPPmzWvXrhUVFY0YMYLr6ggeAokD44J7\nY5oDgAi0tbVt3769X79+XFdEJBBI3DD20IqkrCsc1goAzCKTyZ544gmuayEeeEAfZ+hDKzKKq9jP\nnN1VXJVfdk9/lQcAAckvu59RXMV1Lcw2Nrh3YqTMrEP8/PysVBn7hEDikonnoK+MCTL3/wYATxwv\nvbdLgIEU6GX2X4HOzs7WqIndQpcd9wxOc1iNu5QAwM4gkHjB4DQH3DkLYCdaW1vffvvtJ598cuDA\ngYmJibdu3eK6RtxAlx1f0GkOOkNKq3LLd52rwpASCEuAt5sQO5zHDvB6+E7WMWvWrKKioo0bN7q7\nu6elpcXGxpaUlLi52d3/egQSj2BICcQhMVJmDx/XK1csMye2vLz8u+++++GHH2JiYgghcrk8NDT0\nzJkzUVFRFnl/AUGXHe9gSAlAoNra2p544onCwsLnn3/e29v78ccf/+KLLx561PXr1yUSyYQJE+iP\ncrnc2dm5srLSypXlIwQSH2FICUCINBrNhQsX5syZM2XKlMOHD0+aNCk1NbW0tNT0UVFRUTU1NY6O\nf3wbnzx5UqVShYeHW7++vINA4iksxgogUImJifPmzRs2bNiGDRvc3d1LSkpM7+/i4uLt7U1f//zz\nzwkJCbNnz7bP+20xhsRfJoaUMM0BwGY++OADtVqt/9qY//qv/6Iv3NzcevXq1dzc3Jn36ejo+OST\nTz744IP58+dv2rTJkg0QDgQS361UBAV4u7GXFKKZlDhcRuMKAKzKxcXF4OvO7N/J97l79+7kyZMb\nGxtzc3PHjBnT1ZoKHrrsBCAxUoYhJQCx0mq1U6ZMkclkJSUl9pxGBFdIQkGHlFbnlrNXZMFdSgAi\nkJOT8/PPPxcVFdXU1DCFffv2dXd357BWnMAVkmAEerutjDH8fD9cKgEI18mTJx88ePDEE08EsRw+\nfJjrenHAQcTPgQ8JCbl27RrXtbC8/LL7UZ/pzttZFROEISWwMbH+F+O5bv6z8/m3hisk4dG/S4lg\nRjgACB8CSZAM3qVEu++EuOw/AABBIAkXvUtJP5OwyBAACBQCSdhWKoIMzgiP+uw8h7UCAOgCBJLg\nGey+yy+7hyElABAWBJIYGOu+w4xwABAQBJJ4GOu+QyYBgCBgpQZRod13UdtK2J11WNABrCckJITr\nKoB44MZYEVLWt+k8Cp0QEujthvVYAYDPX4zoshMhY0NK6L4DAD5DIImW/pASwYIOAMBjCCQxM7Gg\nAy6VAIBvEEgih+47ABAKBJJdMNF9l192n6taAQCwIZDshbHuu6Ssy7hUAgA+QCDZERPdd5jpAACc\nQyDZHdp9Ny7Yi12ImQ4AwDkEkj0K9HbbGR+KmQ4AwCsIJDtFu+8MznRIyrrCVa0AwJ4hkOyawZkO\nu4qrgtae5qpKAGC3EEj2jl4q7YwPZRcq69swzQEAbAyBBIQQkhgpK08byS6h0xyQSQBgMwgk+EOg\nt5v+45SitpXgzlkAsA0EEvwHHVLSySTcOQsAtoFAgj+hmcS+SwnTwQHANhBIoMvgXUrIJACwNgQS\nGBDo7TYnUqafSZgODgDWg0ACw4wtfIfp4ABgJQgkMMXgLUqYegcA1oBAgoegtyhh6h0AWBsCCR7O\n4HRwTHMAAMtCIEGnGFz1Dg9SAgALQiBBZxmceocVhgDAUhBIYAZjU+/wcD8A6D4EEphN/0FKGFIC\ngO5DIEFX6E9zIHi4HwB0DwIJusjEw/0wpAQAXeDEdQU6q6OjY+HChRqNhilxdHT84osvOKwS0CEl\nQsiq3P901tEhJf3rJwAA0wRzhVRWVpafn+/v7z+AhetKASH/HlJil9BM2lVcxVWVAECIBHOFdPXq\nVU9Pz/fff9/BwYHruoAu+nA/9vxvZX3b6tzym/Vt9BIKAOChBHOFdPXq1dDQUKQRb+FBSgDQTUIK\nJC8vryVLlowdO3bKlCmffvppa2sr15WCPzH2ICVMvQOAzhBSIOXn5z/66KOLFi0aNmzYl19+mZSU\nxJ7jAHxg8M5ZTL0DgM7g3RhSU1NTdnY286OPj8+ECRPUanVKSkpERMSQIUMIIdOnTw8PD3/33Xdz\nc3NjY2O5qywYZmzqXeJwGYaUAMAYB61Wa433bWhoiIuLmzt3bkJCgv7WvLy87Ozs0tJSqVQql8uT\nk5MDAgLopps3b/7lL39h9oyIiNixY4f+O6jV6qeeemr27NnLly83VoeQkJBr1651uynQRfll96M+\nK9EpHBfstTM+FDPCAbjC5y9Ga3XZpaenV1ZWNjc3629au3btwoULjx49KpVKGxsbd+/ePXXq1IKC\nAro1ICDgPAtNo7q6uqKiIrVazbyJRCJxcnIy+P7AE+OCe+usMEQIyS+7hxnhAGCQJQNJrVYrlcrj\nx48vXbp0+/btBvcpLCzMzMzs37//oUOHsrKyDh48uHnzZpVKlZaW1tZmdIxBqVTOnj37zJkzTMkv\nv/zS3NwcGhpq7BDgA4OrOSjr25KyrmD2HQDosGQg3bhxQ6FQzJ8//+DBg8b22bt3LyFk2bJlfn5+\ntEShUMTGxlZXV586dcrYUUOHDpXL5WlpaYcPH66qqjpx4sSyZcsCAgLi4uIsWH+wBjrNQf9SaVVu\nedRn57mqFQDwkCUnNchksi1bttDXFy5cSE9P19+nsLBQIpGMHTuWXRgdHZ2Tk3P27Nno6GiD7+zi\n4vLFF1+sWbPm9ddf12q1zs7OY8aMWb16tYuLi+kqhYSEMK95221qD+ilUkZxFXumQ37ZvaC1p7HI\nEIC1sb8J+cySgSSVShUKxR/v62TgnVUqVW1trb+/v5vbn76AgoODCSEVFRUm3tzX1/fzzz9vbW2t\nra318/Mz+P76EEL8gYXvALjC/ibkczjZ9D6kxsZGjUbTq1cvnXJa0tDQ8NB3cHd3DwgI6GQaAQ8Z\nfJZS1LaS/LL7HNYKAPjApoGkUqmIoYsnWkK3gujpP0tJWd+WlHUZ0xwA7JxNA0kikRBDwUNL6Faw\nBwYzade5KmQSgD2zaSC5u7sTQvTXoKMldCvYCf0Z4ViMFcDO2TSQpFKpp6fn7du32be4EkKUSiUh\nRCaT2bIywLlAb7c5kTIsxgoAlK0XVw0NDW1vb7906RK7sKSkhBASFhZm48oA57AYKwAwbB1IMTEx\nhJANGzYwJVVVVV9//bVEIhk/fryNKwM8oZ9JdOodMgnArtg6kOLj4wcMGFBUVDRt2rQdO3Z89NFH\nM2bMaGlpmTt3rq+vr40rA/xh7DnoGFICsB+2vqHH2dk5IyNjzZo1eXl5tOPOw8PjzTffTE5OtnFN\ngG8MPged3kWLh1YA2ANrPX7ioVQqlVKp9PDwkMlkVnowOZ9XWQdjlPVtq3PLdZYDXxUThEwCsAg+\nfzFy9sRYZ2fngQMH+vr6WimNQKACvd1WxugOKa3KLcc0BwDRE8wjzMF+GJx6h2kOAKKHQAKeWqkI\nOpYawS7BNAcAcUMgAX/pP3MWKwwBiBgCCXjN4Kp3WGEIQJQQSMB3Bp+DjmkOAOKDQAIBwDQHAHuA\nQALBWKkI2hkfyi7BNAcAMUEggZAkRsr0pzlgSAlAHBBIIDB0SCkx8k8PK8GQEoAIIJBAeAyu5oDu\nOwChQyCBINFpDsdSMSMcQDwQSCBg44J769ylRNB9ByBYCCQQNoN3KaH7DkCIEEggeMbuUkL3HYCw\nIJBAJOgzZ9F9ByBcCCQQDxPdd/ll97mqFQB0EgIJRMVY911S1mV03wHwHAIJREi/+44OKaH7DoDP\nEEggTgYXdMDsOwA+QyCBaBlb0AGz7wD4CYEEYkaHlDD7DkAQEEggfrh5FkAQEEhgF3DzLAD/IZDA\njpjovsONSgCcQyCBfTHWfYcblQA4h0ACu0O77/Sfho6ZDgDcQiCBnaJPQ8eNSgD8gUAC+xXo7bYz\nPlS/+27XuSpkEoDtIZDA3hlbZwiZBGBjCCQAwzMdVuWWR312nqsqAdghBBIAIUZuVMovu4dpDgA2\ng0AC+A+Ds++itpUgkwBsAIEE8Cd09p3OkBKm3gHYAAIJQBcdUtLJJEy9A7A2BBKAAfqPU8LUOwBr\nQyABGGbwcUqrcsuTsq5wVSUAcUMgARhlcOrdruKqoLWnuaoSgIghkAAewuBzKzAdHMDiEEgAD7dS\nEXQsNYJdgungABaHQALolHHBvTEdHMCqEEgAnWVwOjim3gFYCgIJwAz608HJv585i+47gG5CIAGY\nx+B0cNp9t6u4iqtaAYgAAgnAbAangyvr21aj+w6gGxBIAF1k7EFKeGgFQNcgkAC6zuCDlOhDK/LL\n7nNVKwCBQiABdIux7rukrMvovgMwCwIJwALonbMGu+8w+w6gkxBIAJYxLri3/ozw/LJ7mH0H0EkI\nJACLMTYjHLPvADoDgQRgSXRICbPvALoAgQRgeSZm32FICcAYBBKAVRibfYf1WAGMQSABWJGx7jtk\nEoA+BBKAdemvEU4IwZASgD4EEoDVYUgJoDMQSAC2gCElgIdCIAHYDoaUAExAIAHYFO2+GxfsxS7E\nI/4ACAIJwPYCvd12xoei+w5ABwIJgAN0SOlYagS7EN13YOcQSACcGRfcW2dIiaD7DuwYAgmASwZn\nhNPuO2QS2BsEEgDHTMwIx2Nnwa4gkAB4weCMcDx2FuwKAgmAL/RnhGOaA9gVnSwTRwAAGcVJREFU\nBBIAjxicEY5MAjuBQALgl0BvtzmRMv1MwmKsIHoIJADeMTjNAYuxgughkAB4aqUiaGd8KLsEqzmA\nuCGQAPgrMVKGxVjBfiCQAHgNi7GC/UAgAfCdicVYdxVXcVUrAItDIAEIgLHVHFb+gO47EA8EEoBg\n6K/mcOs+hpRAPBBIAEJicDFWDCmBOCCQAATGxGKsGFICQeNpIF24cGHEiBEtLS3swuLi4gULFkyY\nMOGFF17YsmXLgwcPuKoeAOcMPt9vNbrvQMj4GEg1NTUffvhhfX29VqtlCk+ePDl79uyWlpbExMTI\nyMj09PTU1FT2DgD2Rv/5frhLCQTNiesK/Mm3336bmZlZWlqq0Wh0Nq1fvz4sLCwjI8PR0ZEQMmjQ\noGXLlp0+fXrUqFFc1BSAF+iQUkZx1arc/4TQqtzyXeeqjqVE6DyLFoDn+HWFFB4ePn/+/PXr18+Y\nMYNdXldXd/369SlTptA0IoTExsa6uLgUFBRwUU0AHsGQEogGv66QBg0aNGjQIEJIW1vbnj17mPLb\nt28TQgIDA5kSZ2dnPz+/yspKm9cRgI9WKoICvN2Ssq4wJXRI6WZ920pFkIkDAfiDX1dIxtDZDR4e\nHuzCHj16NDc3c1QjAN7BwncgdNxcITU1NWVnZzM/+vj4TJgwwcT+dPKCzhQGrVaLSQ0AbBhSAkHr\nViA1NDTExcXNnTs3ISFBf2teXl52dnZpaalUKpXL5cnJyQEBAXRTfX39+vXrmT0jIiJMB5K7uzsh\npLW1lV3Y2toqk8m6U38A8aHP9yOEsDOJDintjA8bF9ybu6oBPES3Aik9Pb2ystJgv9natWszMzNd\nXV3DwsIaGxt37969f//+rVu3jh49mhASEBBw/rwZj7+kwaNUKpkStVpdUVGBKXYA+ug0B/0hpaSs\ny4nDZRhSAt4yewxJrVYrlcrjx48vXbp0+/btBvcpLCzMzMzs37//oUOHsrKyDh48uHnzZpVKlZaW\n1tbWldVNfH19ZTLZ4cOHmZL8/HyVSjV8+PAuvBuAPcCQEgiO2YF048YNhUIxf/78gwcPGttn7969\nhJBly5b5+fnREoVCERsbW11dferUqa5V9OWXXz537tzatWuvX7+enZ394Ycf9uvXLyoqqmvvBmAP\nsPAdCIvZXXYymWzLli309YULF9LT0/X3KSwslEgkY8eOZRdGR0fn5OScPXs2Ojq6CxWdM2dOQ0PD\nzp07MzMzCSHh4eHr1q1zc3vIIG1ISAjz+tq1a104L4CgmRhSQved/WB/E/KZ2YEklUoVCsUfBzsZ\nOFylUtXW1vr7++ukRXBwMCGkoqKiM2eZOXPmzJkz2SUSiWTp0qWLFi2qrKz08vLy9PTszPsghADo\nkNLYAV5Rn5UwhbT7jhCCTLIH7G9CPoeT5e9Damxs1Gg0vXr10imnJQ0NDd15c4lE0r9//06mEQAw\n9Be+I4Ssyi2P+syM6UUAVmX5QFKpVMTQxRMtoVsBwPYMDinll93DkBLwhOUDSSKREEPBQ0voVgDg\nhImF75BJwDnLB5LBm1iZEroVADik/yh0mkmYEQ7csnwgSaVST0/P27dvq9Vqdjm9rRVrKwDwAe2+\nw11KwCtWWVw1NDS0vb390qVL7MKSkhJCSFhYmDXOCADmopmUGPmnvxGRScAhqwRSTEwMIWTDhg1M\nSVVV1ddffy2RSMaPH2+NMwJAFwR6u62M0R1SonfOclUlsGdWWe07Pj4+KyurqKho2rRpkyZNqqmp\nycnJaWlpeeWVV3x9fa1xRgDoGjrNgejdORu09jQWCAcbs0ogOTs7Z2RkrFmzJi8vj3bceXh4vPnm\nm8nJydY4HQB0k8HFWLGaA9iYg1UfKaRSqZRKpYeHh0wmc3BwsN6JDAoJCcFKDQCdZ3D+96qYIGSS\nmPD5i9G6T4x1dnYeOHCgr6+v7dMIAMylP/WOYDFWsCFhPMIcAGzD4GoO9Mopv+w+V7UCO4FAAoA/\nMbaaQ1LWZcwIB6tCIAGAAQZXc6CLsaL7DqzEupMauMXnsTsAoVDWt2UUV7EnhRNCAr3dMAFPoPj8\nxYgrJAAwhfbgGbxawmQHsCwEEgA8nMF1huhkh13FVVzVCkQGgQQAnRLo7bYzPlR/ssNqLH8HFoJA\nAgAzGOu+QyZB9yGQAMA8Bu9VwpASdB8CCQDMhifPgjUgkACgi1YqgnbGh7JL8ORZ6A4EEgB0XWKk\nDENKYCkIJADoFjqkNC7Yi12INR2gCxBIANBdBmeE55fdw5ASmAWBBAAWYGKaA7rvoJMQSABgMbhL\nCboDgQQAlmRwkSHcpQSdgUACAAsL9HZbGYPuOzAbAgkALI8OKR1LjWAXovsOTEMgAYC1jAvurTOk\nRNB9B8YhkADAigwufIdFhsAgBBIAWBdmhEMnIZAAwBYwIxweCoEEADaCGeFgGgIJAGwHM8LBBAQS\nANiUsSEldN8BAgkAOKA/pETQfWf3EEgAwA0TM8JxqWSfEEgAwBlj3Xe7zlUhk+wQAgkAOGZsRnjU\nZ+c5rBXYHgIJALhnsPsuv+wehpTsCgIJAHgBCzoAAgkAeAQLOtgzBBIA8AvtvtOfEY4hJdFDIAEA\n72BIyT4hkACAj0wMKe0qruKqVmBVCCQA4C+DQ0qrMaQkUggkAOA1/SElTHMQKwQSAPCdwSElTHMQ\nHwQSAAiAwSElTHMQGQQSAAiGsWkOyCRxQCABgJAYnOaA1RzEAYEEAAKDaQ5ihUACAOExNs0BQ0qC\nhkACAEEK9HabEynDkJKYIJAAQKiwQLjIIJAAQNiwQLhoIJAAQPDokFJipIxdiDtnBQeBBABiEOjt\ntjIGd84KGwIJAEQCC4QLHQIJAEQFC4QLFwIJAMQGd84KFAIJAETIxJ2zXFUJHgqBBADiZGxICdMc\neAuBBABihjtnBQSBBAAihztnhQKBBADipz/NgeDOWf5BIAGAXTA4zQF3zvIKAgkA7AUWY+U5BBIA\n2BcMKfEWAgkA7I6xISV033ELgQQA9sjgkBK677iFQAIAO2VsSAndd1xBIAGAXdMfUiLovuMIAgkA\n7B2673gCgQQAYKr7Luqz87hUsg0EEgDAHwx23+WX3cMj/mwDgQQA8B/Guu/wiD8bQCABAPwJ7b47\nlmrgEX9Y+86qEEgAAAaMC+6tf6mk05sHloVAAgAwTH+mQ37ZPQ7rI3pOXFcAAIDXViqC5kTKVueW\nLxjhG9ZXynV1xAyBBADwEIHebjvjQ7muhfihyw4AAHgBgQQAALyAQAIAAF5AIAEAAC/wdFLDhQsX\nFixYcOTIkR49etCSjo6OhQsXajQaZh9HR8cvvviCowoCAICF8TGQampqPvzww/r6eq1WyxSWlZXl\n5+f/9a9/dXP748Y0R0dc3gEAiAe/Aunbb7/NzMwsLS1lXwlRV69e9fT0fP/99x0cHDipGwAAWBW/\nAik8PHz+/PmEkLNnz+7Zs4e96erVq6GhoUgjAACx4lcgDRo0aNCgQYSQtrY2/UDy8vJasmTJ+fPn\ne/bsOXbs2IULF7q7u3NUUwAAsDDBDMNcvXo1Pz//0UcfXbRo0bBhw7788sukpCT9nj0AABAobq6Q\nmpqasrOzmR99fHwmTJhgYn+1Wp2SkhIRETFkyBBCyPTp08PDw999993c3NzY2FirVxcAAKyvW4HU\n0NAQFxc3d+7chIQE/a15eXnZ2dmlpaVSqVQulycnJwcEBNBN9fX169evZ/aMiIgwHUgSieTll19m\nl/zlL39ZvXr1xYsXEUgAAOLQrUBKT0+vrKxsbm7W37R27drMzExXV9ewsLDGxsbdu3fv379/69at\no0ePJoQEBAScP2/Gc67q6upu3LgxbNgwiURCSyQSiZOTk8FTAwCAEJk9hqRWq5VK5fHjx5cuXbp9\n+3aD+xQWFmZmZvbv3//QoUNZWVkHDx7cvHmzSqVKS0tra2vrQi2VSuXs2bPPnDnDlPzyyy/Nzc2h\noVh/FwBAJMwOpBs3bigUivnz5x88eNDYPnv37iWELFu2zM/Pj5YoFIrY2Njq6upTp051oZZDhw6V\ny+VpaWmHDx+uqqo6ceLEsmXLAgIC4uLiuvBuAADAQ2Z32clksi1bttDXFy5cSE9P19+nsLBQIpGM\nHTuWXRgdHZ2Tk3P27Nno6GhzT+ri4vLFF1+sWbPm9ddf12q1zs7OY8aMWb16tYuLi+kDQ0JCmNfX\nrl0z97wAACLA/ibkM7MDSSqVKhSKPw52MnC4SqWqra319/dn1vihgoODCSEVFRWdOcvMmTNnzpzJ\nLvH19f38889bW1tra2v9/PwMnlofQggAgP1NyOdwsvy078bGRo1G06tXL51yWtLQ0NCdN3d3d2em\n6gEAgJhY/sZYlUpFDF080RK6FQAAQIflA4nOzNYPHlrCzNsGAABgs3wg0fXlWltbdcppCVafAwAA\ngywfSFKp1NPT8/bt22q1ml2uVCoJITKZzOJnBAAAEbDK4qqhoaHt7e2XLl1iF5aUlBBCwsLCrHFG\nAAAQOqsEUkxMDCFkw4YNTElVVdXXX38tkUjGjx9vjTMCAIDQWWW17/j4+KysrKKiomnTpk2aNKmm\npiYnJ6elpeWVV17x9fW1xhkBAEDorBJIzs7OGRkZa9asycvLox13Hh4eb775ZnJysjVOBwAAIuCg\n1Wqt9+4qlUqpVHp4eMhkMts/fTwkJAQrNQAAsPH5i9G6D+hzdnYeOHCgVU8BAADiIJhHmAMAgLgh\nkAAAgBcQSAAAwAsIJEHi8wLy3SHWdhHxNk2s7SKibhpvIZAAAIAXEEgAAMALCCQAAOAF694Yyy10\nAQMA6OPtjbFiDiQAABAQdNkBAAAvIJAAAIAXEEgAAMALCCQAAOAFBBIAAPACAgkAAHgBgQQAALyA\nQAIAAF5AIAEAAC8gkAAAgBcQSAAAwAtOXFfA8vLy8rKzs0tLS6VSqVwuT05ODggI4LpSXdfQ0BAX\nFzd37tyEhAT9rUJs7HfffXf8+PHr16/37NlzyJAh8+bN69u3r84+gmuXWq3+6quvTp8+XV5e/sgj\njzzxxBNJSUk+Pj46uwmuXTpaW1tnz57t7e29fft2nU2Ca9rBgwePHTumU9inT5933nmHXSK4dlF3\n7tzJy8s7e/ZseXl5v379ZsyYERMTo7MPD5smtsVV165dm5mZ6erqGhYW1tjYeOPGDVdX161bt44e\nPZrrqnXRxo0bt2/f/sYbbyxYsEBnk+Aaq1arX3/99SNHjkil0sGDB1dXV9+8edPNzW3Hjh3Dhw9n\ndhNcuzo6OhITE4uLi/v06RMcHHzx4sXW1lYPD4+9e/cGBQUxuwmuXfree++93bt3Dx48eN++fexy\nITZt4cKFR48e9fDwYBf6+voeOHCA+VGI7SKE3Lp166WXXqqrq+vbt2+fPn0uX75MCElNTV28eDGz\nD0+bphWRs2fPyuXyCRMmVFRU0JIffvghNDR0zJgxra2t3NbNLB0dHeXl5fn5+UuWLJHL5XK5/PPP\nP9fZR4iN/eqrr+gfYkwNd+/eLZfLn3nmmfb2dloixHbt2rVLLpe//fbbKpVKq9W2tbWtW7dOLpe/\n/PLLzD5CbJeOI0eOhIWFDR06NC4ujl0u0KZFR0f/9a9/NbGDQNv1+++/R0VFDRkypKCgQKPRaLXa\nGzduDB8+fNCgQUxDeNs0UY0h7d27lxCybNkyPz8/WqJQKGJjY6urq0+dOsVp1cxz48YNhUIxf/78\ngwcPGttHiI3dtWuXk5PTRx995ObmRktmzpz5zDPP1NTUXL9+nZYIsV379+/v0aPHihUrnJycCCGu\nrq7z5893cXH55ZdfNBoN3UeI7WKrq6tLS0tLSUnx9vbW2STEprW0tFRUVISGhprYR4jtIoTk5uZW\nVlYuX7581KhRDg4OhJCgoKDk5GQvL6+LFy/SfXjbNFEFUmFhoUQiGTt2LLswOjqaEHL27FmOKtUV\nMplsy7/NmzfP4D6Ca6xWq71z505ISMijjz7KLqedWhUVFfRHwbVLo9H89ttvTz31lFQqZQq9vb17\n9uzp4uKi/XeXuODapWPFihUymSwlJUV/kxCbdv36da1WO2jQIBP7CLFdhJDs7GyJRDJt2jR24auv\nvnr69OnY2Fj6I2+bJp5JDSqVqra21t/fn/nrmwoODias7ztBkEqlCoWCvqZ/dOsQYmPVavXq1av1\n5y+UlpYSQgIDA4kw2+Xo6Hj8+HGdwu+//76urm769OkSiYQIs11sWVlZp06d2rdvH20Om0CbRh+Z\n+thjj23YsOHXX391d3cPCQmZPXv2I488QncQaLsIIUVFRUOHDpVKpf/617/Onz9fWVk5YMCAESNG\nCKJp4gmkxsZGjUbTq1cvnXJa0tDQwEWlrEWIjXVycpo+fbpO4enTp8+cOTNgwIABAwYQYbaL7ddf\nfz1y5MilS5dOnDgRExOzYsUKWi7odimVyk8++WTx4sUDBw7U3yrQptFAWrJkyYMHDwIDA2/fvn3s\n2LGvv/568+bNI0eOJIJtV1NTU3t7u0wmS09P37BhA1Pes2fPdevWTZgwgfC7aeLpslOpVMTQ9QQt\noVtFQxyNPXz4cGpqqqur6wcffMBcSRAht+uXX37ZsWPH8ePHtVqtl5eXWq2m5cJtl1qtXr58+aBB\ng+bOnWtwB4E27dq1aw4ODvPmzTt37tz333//008/LVmypLGx8Z133vn999+JYNv122+/EULOnj37\nt7/97b333isoKCgoKEhLS2tra3vjjTcqKysJv5smnkBif6Ox0RL9rgZBE3pjGxsb/9//+3+vvfaa\nh4fHjh07IiIiaLnQ2zVr1qwLFy4UFxe/9dZbe/bsmThxYn19PRFyu7Zu3VpaWvrxxx87Ohr+rhBo\n0/73f//36NGjr776qqurKyFEIpGkpKRMnTq1pqbmxx9/JIJtF61efX398uXLExISfHx8fHx8Xn75\n5cWLFz948CA9PZ3wu2niCSR3d3dCSGtrq045LaFbRUPQjT169Ohzzz23b9++SZMmZWdns+9AEnS7\nGD179kxOTk5KSrp79+7JkyeJYNt18eLFzz///NVXX/X29m76N41Go9FompqaWlpaiGCb5uPj4+vr\nq1NIbx3917/+RQTbrt69e9MXcXFx7PLJkycTQq5cuUL43TTxBJJUKvX09Lx9+zbTT0IplUpCiEwm\n46Za1iHcxu7cuTMlJcXFxSUjI2PTpk1eXl7srUJs1+XLl99///28vDydcjqF6fTp00SY7SKE/Prr\nr2q1etOmTcNZqqurr1y5Mnz48JdeeokItmn0Bh2dQjpPkvZ6CbRdvXv3dnZ27tGjh84Nv/Q/Gr1e\n53PTxBNIhJDQ0ND29vZLly6xC0tKSgghYWFhHFXKWoTY2KNHj3788cfDhw8/cODA008/bXAfwbVL\nq9V+8803W7du1Smn//n79OlDfxRcuwghgwcPTtUjlUp9fHxSU1NfeOEFupvgmqZUKkNDQ9944w2d\ncjrTgc6vIQJsFyHE2dn56aefbmlpocNFDDqXlVkZiL9N4/CmXIvLzMyUy+UJCQlMyZ07d5588snQ\n0NDKykoOK9YdeXl5BldqEFxjOzo6nn322WHDhjU1NZnYTXDtam1tjYyMlMvlly5dYgrb2tqmTZsm\nl8vz8vJoieDaZcy4ceN0VmoQXNM0Gs3TTz89dOjQK1euMIX37t0bPXp0WFhYWVkZLRFcu6hvvvlG\nLpcvX76cXgVqtVq1Wp2SkiKXy7/77jtawtumiWfaNyEkPj4+KyurqKho2rRpkyZNqqmpycnJaWlp\neeWVV/T7i4VOcI0tLS29efNm//79/+d//kd/a1JSkr+/PxFgu9zc3FasWLF8+fLZs2cnJCQ8/vjj\n1dXVu3fvrqiomDhxIr3ZkAiwXZ0nuKY5ODisWrVq0aJF8fHxL730UkhISFVV1VdffXX37t3Fixc/\n/vjjdDfBtYuaPn16fn7+/v37q6ur6ajYDz/8UFxcPGzYsClTptB9eNs0sS2uWldXt2bNmry8PNo9\n6uHh8eqrryYnJ/N2VsxDHTlyJDU11eDiqsJq7Lfffsvcl6MvKyvrqaeeoq+F1S7q6NGjn3zyCe2F\nJ4R4eXnNmzdv9uzZdBIXJcR26YuKivLy8tJZXFWITTt69Oinn35aVlZGCHFwcOjfv//y5cufffZZ\n9j5CbBch5MGDB+vWrTtx4gTtuPP29p48efJbb73l7OzM7MPPpoktkCiVSqVUKj08PGQyGV3NScTE\n2lghtquxsbGiooLOtTW2jxDb1UlCbNr9+/erqqr69++vMwuATYjtou7cuePg4GBingLfmibOQAIA\nAMER1Sw7AAAQLgQSAADwAgIJAAB4AYEEAAC8gEACAABeQCABAAAvIJAAAIAXEEgAAMALCCQAAOAF\nBBIAAPACAgkAAHgBgQQAALyAQAIAAF5AIAEAAC8gkAAAgBcQSAC20N7eXl1djcePAZiAQAKwiu3b\nt48ZM6agoKCkpGTWrFnDhw8fO3ZsRETEihUrWlpauK4dAB8hkACsoqmpqaamprCwcN68eXV1dQqF\nYvjw4a2trd9+++27777Lde0A+MiJ6woAiNn27dtTUlIWLVrk6OhICPn5559ffPHFH3/8sba29tFH\nH+W6dgD8giskACsaNGjQ4sWLaRoRQp588smwsDCNRnPr1i1uKwbAQwgkACsaPXq0g4MDu2TAgAGE\nkPv373NUIwD+QiABWFHfvn25rgKAYCCQAKxI5/IIAExAIAEAAC8gkAAAgBcQSAAAwAsIJAAA4AUE\nEgAA8IIDVnsEAAA+wBUSAADwAgIJAAB4AYEEAAC8gEACAABe+P9rXjanqMT7ZwAAAABJRU5ErkJg\ngg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "m = 500;\n", "kappa = 10;\n", "lam = linspace(1,kappa,m)';\n", "A = sparse(diag(lam));\n", "b = randn(m,1);\n", "b = b/norm(b);\n", "[xCG,~,~,~,resnorm] = pcg(A,b,1e-14,100);\n", "semilogy(resnorm,'.-')\n", "hold on\n", "ylim([1e-16 1])\n", "xlabel('n')\n", "title('Convergence of CG')\n", "legend('||r_n||_2')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this simple problem, we can calculate the true solution and therefore the error as well." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAABcSAAAXEgFnn9JSAAAA\nB3RJTUUH4AsHFRA5P71DSAAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJ\nbmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAwNy1Ob3YtMjAxNiAxNjoxNjo1N7HH9vsAACAA\nSURBVHic7d15XNTV4j/+A8OwDqgY6LCIJoLgVqKWSwpiorldXIg+4u6lxNRM63ajMi1aNC29P78q\nD03gVpdrammIqVgouICGpSnqRRgVAlFRQBhkmJnfH8fevZ2NAWZ4L/N6/vWeM+/3vM+JaV7OmbPY\nabVaAgAAwDV7risAAABACAIJAAB4AoEEAAC8gEACAABeQCABAAAvIJAAAIAXEEgAAMALDlxXAISn\nvLw8Nze3oKCgoKBAq9UGBwcPHDhw5syZjo6OXFcNrOLs2bNXrlypq6vr06fP8OHDTZ/cireHziW9\ne/ceOnRoTEyMRCKxdFOA37QALbFjxw43Nzf9N1JAQMB//vMfrmsHFnb79u3w8HDmr/ziiy+aPr8V\nbw9jlwQHBx84cMAKbQL+QiCBuR4+fPjiiy+a+MeNnZ0dMklk5syZw/4TmwikVrw9mr1EKpVmZWVZ\nv5XAF3ZaLB0E5lm/fv3KlSvpcZcuXWbOnDls2DB7e/vMzMydO3eq1WpCiLOz85kzZ/r27ctpTcFi\nOnfuXFVVRQhJSEh4//33HR0dO3ToYPDMVrw9DF5y796977///vDhwyqVihDSoUOHwsJCuVxu7ZYC\nL3CdiCAMt27d8vDwoO+ZJ5988tq1a+xn9+/fz3T3v/nmmwZfQaFQlJSUaDQaM29XVFSkVqvbWO2m\npqYrV67U1NSYc3JpaemVK1f0C8vKyix+LzMbWFxcXFRUpFKpLHVffSb+LnV1dcwHxbFjx0y8SCve\nHqYv+fjjj5lbb968uXVNA8FBIIFZmH/JEkKOHDmif0J0dDR9Njg4mF1++fLl8ePHM/+sdnd3f/75\n5y9cuMA+Z9q0ab6+vr6+vt9///2hQ4d69epFT3ZxcZk3b96DBw/oabGxsfS0hQsXsi9PSkqi5fPm\nzWMKT58+/dxzz7m6uhJC7O3t+/Xr9+GHH7I/dtk3LSgoGDx4MCFELpfTZzUazQcffODn50dr4ufn\nt2vXrh07dtBL4uPj2RVo0b1MNJBRWFgYFRXF/EeTSCRTpkw5c+aMzmnN3teEZv8uvr6+Pj4+zB/d\ny8vL19f3X//6l8FXa8Xbw/QlDx8+DAkJcXBwcHBweOGFF8xpEYgAAgnMEhERYTBvGLW1tX/8iSlM\nTU11dnbW/14ulUq3bNnCnPbcc8/R8piYGP2xWFOnTqWnbdmyhfn0fPjwIXN5nz59aHlKSgot+eST\nTwwO0Jo+fXp9fb3OTdesWePu7k6PmUCKiYnRudbBwWH06NH0eNq0aczdW3Qv0w2kvv76a5ox+hU4\nePBgi+5rjDl/F/1nCSEfffSRwRdsxduj2UvABiGQwCxdunShHx8zZ84085Lr168zo6dGjx6dkpKS\nlpYWFRVFS5ycnC5fvkzPZD6vCSHu7u4JCQlLly6VyWS0xN7e/saNG1qttrKy0sHh0USFH3/8kV5b\nVFRESxwdHe/fv6/Vak+cOGFnZ0dvkZiY+N133y1fvtzLy4uetmbNGp2bMpUMCAgYPHiwVqv94Ycf\nmPrMmDFj9+7dW7ZsCQkJYQqZQGrpvUw3UKvV/vHHH0z5qFGj0tLSNm3aFBYWRktkMll1dbX5923L\n3yUjI2P37t1MtT/77LOMjIyioiJLvT1acQmIHgIJmnfv3j3mg+mtt94y86q4uDh6ybBhwxobG2lh\nU1PT888/T8snT55MC5nP686dO9NQ0Wq1P//8M3NT5pvB2LFjacmiRYtoyfr162nJpEmTaMnTTz9N\nS77++mumMvv27aOFXl5e9AsEOyQmTZp0+/Zt5mTmKaaGWq32zp07HTt21Amklt6r2QYuXLiQljz9\n9NMNDQ20sLy83MXFhZZ/99135t+3jX+XBw8eMDU8deqUsRdsxdvDxCXR0dEBekpLS815WRA6rNQA\nzXN2dmZ6hxoaGsy86tSpU/QgISFBKpXSY4lE8uqrr+qcwHjhhReYXzUGDx7M3PTWrVv0IDY2lh7s\n379fq9USQphPYTqAuLq6+ty5c4QQe3t7f3//03/q3Llzp06dCCG3b98+evQo+6ZPPPHEf//73yee\neIIp+fXXX+nByy+/zBR27tx52rRp7Atbca9mG3js2DF6sHjxYicnJ3rctWvXn3/+OSMjIyMjo0+f\nPq24L1sr/i6mteLtYeKSioqK63qamppaVCUQKKzUAM1zdnbu0aMH7RxTKBQGzykpKSkoKKDH48eP\npyX0IbuzixASGhpKD27fvn379m2mo4kQwh7d6+bm5ufnd/36dUJIY2MjLYyOjn7llVcaGxvLysp+\n+eWXHj16nDhxgtZw8uTJhJDff/+dnqnRaEaOHGmwquXl5eyH48ePZ75/0FrV1tbS44CAAPaZPXv2\nZD9sxb1MN7CxsbG4uJg++9RTT7EvfOaZZ5hj2uQW3ZdRX1/fir+Laa14e7i6uhq7pFevXvX19YSQ\nhoaGK1eumFkHEAcEEpilT58+9OPj6NGjNTU1zIBdxsqVK/fu3UsI8fDwqKqqqqur02g09Cmd397Z\nD+lcEwbznYCyt9f9Bt+xY8eoqCj6G8++ffsCAwPpBJfx48fTgQkPHz5krg0MDDTYFqZilK+vL/uh\nh4eHnd2j+XnV1dXsp9jDoFt3L9MNrKuro80hhDDfXfS14r6MpqamVvxdmtXSt4eJS1JTU+nB/v37\np0yZ0qJqgNAhkMAsEydOpJ1jtbW1mzZteuedd9jP3rhxg+k6GzlypEQi8fDw8Pf3v3nzJiGkqKho\nwIABzMnMMISOHTuyBxabKTY2lh1ItJCZ8M+MuNNqtb/88gszQMAEnY9+Jyenrl270m8YWVlZw4YN\nY546cuQI+8xW3Mu0Tp06yeVyeuurV6/279+feSo1NfX8+fOEkDFjxgwcOLDV97XS36Wlb49mLyGE\npKent6gOIAL4DQnMMn/+fOZz8L333lu7di3TrX///v0lS5Yw/7Rnfmhhxobt2LGD/VLbt2/XOaFF\nJk+eTHvYLly4cODAAUKIq6vrxIkT6bNdunShk4e0Wu3+/fuZqzQazbJly6ZOnTp16tSKiopmb0EP\n1q9fn5ubSwhRq9WrVq3Kz89nn2aRe+kYNGgQPfjyyy+Zwnv37i1evHjDhg0bNmzQarVtvK81/i6t\neHuYuKSpqWnt2rX/+c9/WloNEDzOhlO0XH5+fnx8fGRk5IwZMzZu3MiMQYL2kZubS4caU66uruHh\n4aNHj2bGnhFCJkyYwJz/+++/M18+4uLijh49+vPPPy9YsICW2Nvb5+fn0zOZQWjvvvsu+449evSg\n5cnJyezy6dOns9/DM2bMYD/71Vdf0XIPD4+XX345Kytr9+7dTMZERESYvqlWq7127RozvpwQ4ufn\nR8dJM1OImFF2bbyXfgN//fVX5tZz5849fvx4SkrK0KFDaUnXrl2VSqX59zXI/L+LmaPsqJa+PYxd\nMnbsWP3vZwqFwvTdQRwEE0jHjx8PDg6Oi4v797//vXbt2n79+s2fP9/MSelgKUePHtX5xYWtX79+\nlZWV7PPXrVun/zsQtWrVKua0lgYSe34MIeTbb7/Vqaf+tFbK29v70qVLpm/KtLRz587sa8eMGZOU\nlESP2RNj23Ivgw385JNP2B/TDKlUmpub26L7GmPm36VFgaRt+dvD9CVBQUHz5s2jxwgkGyGYQJo4\ncWJ0dDSz9tf+/fuDgoLY/39C+7h79+7cuXN1/g3r6em5YcMGg0uu5ebmPvXUU8wnrJ2dXd++fX/6\n6Sf2OS0NpPr6euaHE5lMZnDOzc6dO9mfdPb29tHR0cxUXBM3ZZSXl3/11Vevvfbam2++uXfv3ocP\nH65Zs4ZeEhcXZ5F7GWvgsWPH+vfvz46lyMjIs2fPtrSNJpjzd2lpIGlb/vYweIm9vf3kyZPv379/\n6NAhBJJNEcZq33fu3Bk+fPhbb73F/ItJpVINHDgwLi7uH//4B7d1s1mVlZXnzp2zs7MLDQ1l1nwz\npq6u7tKlSxqNpk+fPm3/8d98d+7cuXjxopub25NPPunp6WnmVSdPniwrKyOEdO/ena5xR40cOTIn\nJ4cQ8vbbbzPfltp4LxNqa2svXLgglUp79erF7viy4H2t93dp0duDfYmTk1NYWBizmBPYFGEE0rlz\n52JjY7du3cqsf0UIGTduXFBQ0KZNmzisGIjSp59++tZbbxFCZDJZQUEBXQv1zJkzw4YNa2pqsrOz\n++233/r168d1NQHERhij7OhEOZ1tJV1dXXXmhQBYxMyZM+kM1gcPHoSEhISFhQ0ePHj48OF0GNjC\nhQuRRgDWIIxAol/jdL7M0T5HjmoEYubn53f8+PERI0YQQtRqdUFBwdmzZ1Uqlaur65YtW5KTk7mu\nIIA4tevE2Orq6ujo6Pnz5zPLO7JlZWXR5YRlMllQUNCCBQuYhVvovBOlUsk+X6lUYh9JsJLAwMCc\nnJxLly4VFBSUl5fLZLKAgIDw8HCDG0MAgEW0ayBt3769rKzMYD9bUlJSWlqak5NTaGhoTU3Nrl27\n9u3bt3nzZvqvVBo87DWv1Gp1aWnp8OHD26vuYItCQ0OZFd4AwNqs3mWnVqsVCsWxY8eWL19urK8j\nLy8vLS2tW7duBw8eTE9Pz8zM3Lhxo0qlSkxMpCsB+/j4yOVy9sIt2dnZKpWKmdYOAABCZ/VAKi4u\njoqKio+Pz8zMNHbOnj17CCErV65kJlVERUWNGzeuoqKCWdh49uzZZ8+eTUpKunr1akZGxocffujv\n788edGfM97/ftlvx09BNv/RIOmmJBgEAgFVYvctOLpczI7PPnz/PrJfFlpeXJ5FIRo0axS6MjIw8\ncODA6dOnIyMjCSFz5syprq7euXNnWloaIaRfv34ff/yxwW2YGcHBwYSQu73Gk+DJp69XE0KemLSi\n89UMC7UMAECQeLuvh9UDSSaTMbsjs9cHY6hUqsrKSj8/P510oXvPlJaW0ocSiWT58uVLly4tKyvr\n1KmTmfPmrly5EvH/zmVfe7Q9pfuz0w6lJnX3NBVjghAcHMzbt1RbiLVdRLxNE2u7iHibRv+lzk/c\nD/uuqanRaDTMNpoMWqKzIY1EIunWrVuLZnGviurBHCuqGualF7ahsgAAYC3cBxLdCkz/yxMtaelG\nYfrCe3Z8f+xfmZR97V7KGcObaQIAAIe436CP7talHzy0RGdTy1Z4qLj40tfTvLRBtRLXtzrHE0JW\nHy4J79lJBB13AABiwv03JIOTXpkS+mxb1F88SeprxijPRj84vuT+XkKIoqohYktBG1+WW6Ls2ibi\nbRcRb9PE2i4i6qbxFveBJJPJ3N3db968yewpSdFpsG1fi6Eu/xBzHF2XM6ShkODHJAAA/uE+kAgh\nISEhjY2NFy9eZBcWFBQQQto+T75zzArm2Lfp9id3k32b7hBCUs6Urz5U0sYXBwAAS+FFII0dO5YQ\nsm7dOqakvLz8m2++kUgko0ePbuOLu/QZqpNJ/771aCeblLPl2dfut/H1AQDAIngRSLGxsYGBgfn5\n+VOnTt2xY8dHH300ffr0+vr6+fPn62w92TqdY1a49hnGPPRtus38mDQv/VLbXx8AANqOF4EklUpT\nU1OjoqIuX768du3a1NRUpVK5YsWK5cuXW+oWXRZ/LvXyZx7ixyQAAL7h146xKpVKoVC4ubnJ5XI7\nO7s2vprORGtV5c2ShGeYh2UOXqN9P6fHO2ND5g7GThYAIH58XoGCF9+QGFKptFevXj4+Pm1PIwMv\n7u2vP8CBHq8+XKKoarD4HQEAwHz8CiRr0/kxiT0zCR13AADcsq1AIoZ+TKKjwLOv3cMocAAADtlc\nIBnquNtGj1POlqPjDgCAKzYXSIQQj4gYdiYNaShExx0AAOdsMZAIIR7hMQZHgaPjDgCAKzYaSFJv\n/66vfsE8ZI+4Q8cdAAAnbDSQiKElhdBxBwDAIdsNJEKIR3jMY6PAWSPusIkfAEA749dKDZZlzoRk\nneUbvpONpJv4dfd0LkkcZvw6ABsVHBzMdRWgTXs18XmlBu53jOWW1NvfIzymJnsXfRj94Phet+fy\nnUMUVQ2rD5Wsiuph+nIAG8TbjzMbIeJ/E9h0lx3F/iWJEILRDQAAnEAgEam3f9fFj424Y0Y3rD6M\nIeAAAO0EgUQIIS59hhoc3ZByBjv4AQC0EwQSIYRIvf27LP6ceejbdHtJ9V56jB38AADaBwLpETq6\ngXkY/eA4s4MfhoADALQDBNJfOsesYK8nhN2SAADaEwLpL/oLgWN0A4Btev/99ydOnKh/zNXr2AgE\n0mM8Igyv3YDRDQA2RalU1tbW6h9z9To2AoGkC6MbAAA4gUDSZWJ0AxZdBQCwHgSSAcZGN2Rfu4eO\nOwBgiHgtUE4gkAwwtM15MqGjG7B9H4ANmzJlyvLly+vr6xctWhQcHLx3716uayQqCCTDdEY3DGko\nZLaUxbQkAJt19erV69evz507NzU1tUOHDj4+PlzXSFQQSEbpjG5gT0viqEYAwL2DBw9ev379ypUr\n+fn5Q4cO5bo6ooJAMsrEoqsY3QBgsxoaGpKTk/39/Zs/FVoIgWSKsUVXMboBwGbJ5fIBAwZwXQtx\nsvUN+kyji64yW8rSaUlvdY6noxvCE57mtnoA/JF97X6qdX5ebVBpmjRamZPEGi9OCBnVs+PcwXLz\nz/f19bVSTQCB1AxjW8rS0Q0teh8DiNixIqGO9+neyblF50ulUivVBNBl1zwsugoA0A4QSM3DoqsA\nYHFKpfIf//jHU0891atXr7lz5964cYPrGnEPXXZm8YiIqcneVX/xJH0YXZeT5xyS7xyScqZ8zmB5\neM+O3FYPgHMBns4C7cEeFdiJk/vOnDkzPz9//fr1Li4uiYmJ48aNKygocHZuWf+hyCCQzKU/umGW\ncyIhZF76pZLEYSYvBRC/uYPlAg2kFikstMyUj5KSku+///7HH38cO3YsISQoKCgkJOTUqVMREREW\neX2BQpeduXQ67oY0FEY/yCFYTwjAtjU0NAwYMCAvL2/KlCmenp5PPvnktm3bmr3q6tWrEolkzJgx\n9GFQUJBUKi0rK7NyZfkOgdQCHuEx7NENzM4UKWexWxKAjdJoNOfPn58zZ86kSZOOHDkyYcKEhISE\noqIi01dFRETcunXL3v7RJ3BOTo5KperXr5/168trCKQWkHr7d331sbUbsOgqABBC5s6du3DhwrCw\nsHXr1rm4uBQUFJg+39HR0dPTkx7/+uuvcXFxs2bNwnxb/IbUMnTthr9GN2BaEoAYffDBB2q1Wv/Y\nmCFDhtADZ2fnDh061NXVmfM6TU1Nn3766QcffBAfH79hwwZLNkCY8A2pxdiLrpLHpyVxUR0AsDxH\nR0cXFxf9YxPnt/R1bt++PXz48K+++urw4cObNm1ycMDXAwRSy2HRVQBoI61WO2nSJLlcXlBQMHLk\nSK6rwxfI5NbQ7bjDtCQAaIkDBw78+uuv+fn5t27dYgq7dOnS7FcxccM3pNagi64yD+m0JHo8L/0S\nR5UCAMHIycl5+PDhgAEDerAcOXKE63pxDIHUSiamJaHjDsB2uLq6arXaESNGMCVlZWXz5s0zfdWn\nn36q1TN58mQrV5bvEEitpz8tie6WlHIG05IAAFoMgdR6+tOS0HEHANBqCKQ20d1S9sFxdNwBALQO\nAqmtdKYlMR132OYcAKBFEEhtZWBaUjUzLQkddwAA5kIgWYBHRIyxjjuscQcAYCYEkmXod9zRg5Sz\n5djmHADAHAgky9DvuGMWAsfoBgDBef/99ydOnKh/zIdXEzEEksXod9wNaSgkhGRfu4eOOwBhUSqV\ntbW1+sd8eDURQyBZkrGFwNFxBwDQLASSJaHjDgCg1RBIFqYzVXZIQyHTcZdyppy7egEA8B0CycL0\nFwJn7+CHjjsAAGMQSJansxA4ewe/iC0F3NULACxJo9GcP38+JyeH64qIBwLJKjzCHx9xV5dDO+4w\nVRZABGpqal5//fXnn38+Kyvr6tWrNTU1XNdIJBBIVmGi4y7lLDanABCwe/fu0d2Pjh49+vrrry9Y\nsMDDw4PrSokEAslaTHTcYY07AOFatWrV3bt3P/vsM64rIkIOXFdAzDzCY5QXT9VfPEkfRtfl5DmH\n5DuH0FHgO2NDuK0egAUpL56q/vm/XNeiNVz7DPOIiDHzZI1Gs3Xr1u7du8+cOZOW/POf/+zfv7/V\namdbEEhWRDvuShKeoQ9px91o388JISlnyucMlof37MhpBQEspv7iyZrsXVzXojWk3v7Nn/Sn+vp6\nlUr17rvvzpo1y3pVslnosrMuYx13hBCMbgAQHGdnZ4lE4ubmxnVFxAmBZHX6I+6wgx+AQDk4OMyc\nOTM7O9uyL3vjxg2JRNK7d2/LvqzgoMvO6gx13G2b1SWREDIv/VJJ4jCTVwMIg9TL3yPc3F9ieIX9\n70VzbN68efbs2UlJSQMGDDh16lRAQEB8fHwb65CamhoaGlpYWHjq1KmhQ4e28dWEC4HUHmjH3d1d\n6+nDIQ2F0Q9yvpM9p6hqSDlTPnewnNvqAbSdR0SM+UMDBE0mk+3du/fu3buXL19evXq1g0NbP0W1\nWm1KSsqKFSu+++67L7/80pYDCV127cQjPEbq9ddvp8wOfqsP45ckAOHp3Lnz8OHDmTRqaGgYMGBA\nXl7elClTPD09n3zyyW3btpn5UsePH79+/fqMGTNefPHF//73v/X19VarNd8hkNqJ1Nu/66tYCBxA\nnOgyQnPmzJk0adKRI0cmTJiQkJBQVFRkzrU7d+6MjIz08vKaOnVqQ0PD7t27rV1b3kIgtR/9hcAx\nugFATObOnbtw4cKwsLB169a5uLgUFDS/duWDBw9279790ksvEUI8PT3HjBmzc+dO69eUp/AbUrvS\nGd2wpHrvW53j6QJ34QlPc1s3AGB88MEHarVa/9i0IUOG0ANnZ+cOHTrU1dU1+2q7du2qq6srKiqi\nSz9IJJJjx44VFxc/+eSTlmqLgOAbUruSej82Eom9zTm+JAHwh6Ojo4uLi/5xs1e19NV27tzZrVu3\nnJycH3744YcffqipqXFyckpJSWlT7QULgdTe2PNkCWubcyxwB2Br/ve//+Xm5m7ZsuUYy4wZM1JT\nUzUaDde14wACqb3pb3Me/SCHYGcKANuTkpLStWvXqKgodmFcXNyNGzd++uknrmrFIQQSB3RGNzBD\nwFPOlmNLWQAbodFo0tLSZs6cKZFI2OWRkZFdu3b98ssvuaoYhzCogQN0nmz9qkergNMh4HR0w7z0\nwp8xugFAaFxdXbVaLbukrKzM9CX29vY3b97UL5dIJOXl5ZasnHDgGxI39IeAY3QDANg4BBJndLaU\nZTruMLoBAGwTAokzOjtTMF+SsHYDANgmBBKXdBa4Y4aAo+MOAGwQAolL+tv3MQvcYQg4ANgawYyy\na2pqWrx4MXuymL29vfnr6fKWR0RMTfau+ouPRtzRjrt855Dsa/ewMwUA2BTBfEO6du1adna2n59f\nIAvXlbIMY6MbsDMFANgUwXxDunz5sru7+3vvvWdnZ8d1XSzMxPZ9qw+VrIrqwW31AHQEBwdzXQUQ\nJyEFUkhIiPjSiPIIj6n5eZfq9qNZckuq934ne44QknK2fM5geXdPZ05rB/CXK1eucF0FEC3BdNld\nvny5U6dOr7322qhRoyZNmvTZZ58plUquK2UxJrbvQ8cdANgIIQVSdna2t7f30qVLw8LCvvzyy3nz\n5olpQVyXPkMN7kyRcqYcQ8ABwBbY6ay/xLna2tqMjAzmoZeX15gxY9Rq9ddffz1w4MC+ffvS8j17\n9rz99tsbN24cN26csZcKDg4WVveCqvIms30fIaTMwWu07+eEkO6eziWJw4xfBwBgLj5/MForkKqr\nq6Ojo+fPnx8XF6f/bFZWVkZGRlFRkUwmCwoKWrBgQUBAAH3q+vXrf/vb35gzBw4cuGPHDv1XUKvV\nTz/99KxZs9544w1jdeDzf3djan7eVbH5Nebh/9dh6r86TiWEzB0s3xkbwl29AEAk+PzBaK1BDdu3\nby8rK2N28GVLSkpKS0tzcnIKDQ2tqanZtWvXvn37Nm/ePGLECEJIQEDAuXPndC65c+dOcXFxWFgY\ns067RCJxcHAw+PqCRhddZaYlRdfl7JWNLHN4IuVM+aieHTEtCQBEzJK/IanVaoVCcezYseXLlycn\nJxs8Jy8vLy0trVu3bgcPHkxPT8/MzNy4caNKpUpMTGxoMLoVkEKhmDVr1qlTp5iS3377ra6uLiRE\nbF8apN7+JqYlYbckABAxSwZScXFxVFRUfHx8ZmamsXP27NlDCFm5cqWvry8tiYqKGjduXEVFxYkT\nJ4xd1b9//6CgoMTExCNHjpSXlx8/fnzlypUBAQHR0dEWrD9PSL39DY5uwKKrACBuluyyk8vlmzZt\nosfnz5/fvn27/jl5eXkSiWTUqFHswsjIyAMHDpw+fToyMtLgKzs6Om7btm3NmjVLlizRarVSqXTk\nyJGrV692dHQ0XSX2DD7edpvq6xyzQnnxFDMt6ZO7yXR0A110NbxnR05rBwACI5S5zJYMJJlMxmwO\n7+Bg4JVVKlVlZaWfn5+z82MzPXv27EkIKS0tNfHiPj4+W7duVSqVlZWVvr6+Bl9fn4BCiI2u3cCM\nbmC2lCWEzEu/hBF3ANAi7E9CPodTu85Dqqmp0Wg0HTp00CmnJdXV1c2+gouLS0BAgJlpJGgeETEG\nt5RVVDWknLHR7Y0BQNzaNZBUKhUx9OWJltBngaEzuoHZLQlrNwCAKLVrINFB2/rBQ0uYId1ASb39\nuy5+bD2hJff3EoxuAACRatdAcnFxIYTor0FHS+izwEanJTEPo+tyfJvuEEKyr93DEHAAEJl2DSSZ\nTObu7n7z5k21Ws0uVygUhBC5HLM+dRmbloRFVwFAfNp7cdWQkJDGxsaLFy+yCwsKCgghoaGh7VwZ\nQTA2LQmLrgKAyLR3II0dO5YQsm7dOqakvLz8m2++kUgko0ePbufKCEXnmBXsTytFtQAAHsFJREFU\nh8zohnnpl7ioDgCAVbR3IMXGxgYGBubn50+dOnXHjh0fffTR9OnT6+vr58+f7+Pj086VEQr90Q3R\nD3II7bg7hI47ABCJ9g4kqVSampoaFRV1+fLltWvXpqamKpXKFStWLF++vJ1rIiw6oxuYBe5SzpZj\ndAMAiANn+yGpVCqFQuHm5iaXy620MTmfV1lvBeXFUzdXTWMeficbSdduwM4UAGA+Pn8wcrZjrFQq\n7dWrl4+Pj5XSSHx0viQxazdgdAMAiINgtjAHord2A9Nxh9ENACACCCQhoYuuMg+xwB0AiAkCSWA8\nwmOkXv7MQ/YCdxjdAACChkASGJ0vSewF7rB2AwAIGgJJeHR2pmAWuMPoBgAQNASSIOntTLGNHmN0\nAwAIFwJJkPRHN2DtBgAQOgSSUOmMbmCv3YCOOwAQIgSSUOmPbqAj7hRVDei4AwAhQiAJmO7ohj93\npkDHHQAIEQJJ2NijGwhrWhIWXQUAwUEgCZv+zhSsjrtC7uoFANBiCCTBM7boava1e1hPCAAEBIEk\neFJvf71pSVhPCACEB4EkBvodd1hPCAAEB4EkEjodd9PqjjG7JaHjDgAEAYEkEjodd/Kmu+i4AwBh\nQSCJBxYCBwBBQyCJike47kLg6LgDAKFAIIkKRtwBgHAhkMTGRMcdpsoCAJ8hkETIWMdd9rV7WOMO\nAHgLgSRCJjrusMYdAPAWAkmcTHTcRWwp4K5eAABGIZBEy1jHHTanAAB+QiCJlumOO+wqCwB8g0AS\nM5Mj7rCrLADwCwJJ5Ex03GEUOADwCgJJ5Ex03GVfu4eOOwDgDwSS+Jla4w6jGwCANxBINsHEVFl8\nSQIAnkAg2QQTHXcY3QAAPIFAshX6HXfRD3IIOu4AgDcQSDZEp+NuSfVeeoD1hACADxBINkT/SxLt\nuMMOfgDABwgk2+LSZyj7S9KQhkLfpjsEoxsAgAcQSDZHZ3QD7bjDL0kAwDkEks2Revt7hMcwD6Mf\nHMcQcADgAwSSLWL/kkQIwRBwAOADBJItknr7d138BfMQQ8ABgA8QSDZKZ3TDkuq9dHQDhoADAFcQ\nSDbKwAJ31czOFFgFHAA4gECyXS59hhob3ZByppy7egGAjUIg2TRjoxtWHy5Bxx0AtDMEkk3TH93w\n184UWLsBANoXAsnW6YxuiK7LeTS64Uw5Ou4AoD0hkGydoQXuttFjdNwBQHtCIAFx6TOUnUlDGgr/\nmpaEjjsAaC8IJCCEEI/wGKmXP/Pwr2lJ6LgDgPaCQAJC6OiGVx8b3YCOOwBoZwgkeERnWhK74w5T\nZQGgHSCQ4C+dY1YY7LjDQuAA0A4QSPAXEx13WAgcAKwNgQSPQccdAHAFgQS60HEHAJxAIIEuYx13\niqoGdNwBgPUgkMAA/Y47uhA4dvADAOtBIIFhOh13zELgKWfL0XEHANaAQALDDK1xl0zQcQcAVoNA\nAqM8ImLYC4Gj4w4ArAqBBKZ0Wfw5c8x8SSKEpJwtx3pCAGBZCCQwxcQOfpiWBACWhUCCZhjbwS/7\n2j0sBA4AFoRAgmZIvf31Ou7+Wgico0oBgAghkKB5OiPuHtvBD6MbAMBCEEhgFqM7+GF0AwBYCAIJ\nzKK/ntCS6kejG9BxBwAWgUACc+msJxT94DidlpRyBms3AIAFIJCgBdi/JBHWekJYuwEA2g6BBC2g\nPy0JoxsAwFIQSNAyOtOS6C9JBKMbAKDNEEjQMiYXXcXaDQDQeggkaDFjuyVh7QYAaAsEErSGwS9J\nhJDVh0vQcQcArYNAgtYwsegqpiUBQOsgkKCVjC26imlJANA6CCRoJUOjGx4tuoppSQDQCjwNpPPn\nzw8dOrS+vp5deObMmZdffnnMmDExMTGbNm16+PAhV9UDytjoBkxLAoBW4GMg3bp168MPP6yqqtJq\ntUxhTk7OrFmz6uvr586dO3jw4O3btyckJLBPAE50jlnBXnSVvaUsOu4AoEX4FUjffvvtpEmTwsPD\nf/vtN52n1q5dGxoampqaGhcX98YbbyQlJeXm5p48eZKTegLDxLQkfEkCgBbhVyD169cvPj5+7dq1\n06dPZ5ffuXPn6tWrkyZNsrd/VOFx48Y5Ojrm5uZyUU14jEdEDHt0A6YlAUDr8CuQevfuPWnSpEmT\nJj311FPs8ps3bxJCunfvzpRIpVJfX9+ysrJ2riEYpLel7F/TkjiqEQAID78CyRg6usHNzY1d6Orq\nWldXx1GN4DEmpiVhPSEAMJMDJ3etra3NyMhgHnp5eY0ZM8bE+XTwgs4QBq1Wi0EN/EGnJdVffPSr\nXnRdzl7ZyDKHJ1LOlM8ZLA/v2ZHb6gEA/7UpkKqrq6Ojo+fPnx8XF6f/bFZWVkZGRlFRkUwmCwoK\nWrBgQUBAAH2qqqpq7dq1zJkDBw40HUguLi6EEKVSyS5UKpVyubwt9QcLoqMb6lc9CiS6pexbneMJ\nIfPSL5UkDjN5NQBA27rstm/fXlZWZrDfLCkpafHixT/99JNMJqupqdm1a9fkyZOZMQgBAQHnWHbs\n2GH6RjR4FAoFU6JWq0tLSxFIvGJsS1mMuAMAc7Q4kNRqtUKhOHbs2PLly5OTkw2ek5eXl5aW1q1b\nt4MHD6anp2dmZm7cuFGlUiUmJjY0tGblTR8fH7lcfuTIEaYkOztbpVINGjSoFa8G1oNpSQDQai0O\npOLi4qioqPj4+MzMTGPn7NmzhxCycuVKX19fWhIVFTVu3LiKiooTJ060rqKzZ88+e/ZsUlLS1atX\nMzIyPvzwQ39//4iIiNa9GlgJpiUBQKu1+DckuVy+adMmenz+/Pnt27frn5OXlyeRSEaNGsUujIyM\nPHDgwOnTpyMjI1tR0Tlz5lRXV+/cuTMtLY0Q0q9fv48//tjZ2dn0VcHBwczxlStXWnFfaCmPiJia\n7F3M6AY6LSnfOYROS5o7GL2sAO2N/UnIZy0OJJlMFhUV9ehiBwOXq1SqyspKPz8/nbTo2bMnIaS0\ntNScu8yYMWPGjBnsEolEsnz58qVLl5aVlXXq1Mnd3d2c10EIcaLL4s9LEp6hx/RL0mjfzwkhqw+X\nhPfs1N2zmX9GAIBlsT8J+RxOlp+HVFNTo9FoOnTooFNOS6qrq9vy4hKJpFu3bmamEXAFuyUBQCtY\nPpBUKhUx9OWJltBnQfSwWxIAtJTlA0kikRBDwUNL6LMgelJvf531hJZU76XH2C0JAAyyfCAZnMTK\nlNBnwRZIvf0xLQkAzGf5QJLJZO7u7jdv3lSr1exyOq0VU1ltiolpSYqq1sxIAwARs8riqiEhIY2N\njRcvXmQXFhQUEEJCQ0OtcUfgJ/1pSRjdAADGWCWQxo4dSwhZt24dU1JeXv7NN99IJJLRo0db447A\nWzq7JUXX5dCOO4xuAAAdVgmk2NjYwMDA/Pz8qVOn7tix46OPPpo+fXp9ff38+fN9fHyscUfgM4xu\nAABzWCWQpFJpampqVFTU5cuX165dm5qaqlQqV6xYsXz5cmvcDnhOp+NuSENh9IMcgtENAPA4O6tu\nKaRSqRQKhZubm1wut7Ozs96NDAoODsZKDTyhqrxZumq66vZN+rDMwYuu3dDd0/nnRQOxdgNAu+Hz\nB6N1d4yVSqW9evXy8fFp/zQCXpF6+3d99bG1G5hFV7GlLABQwtjCHERAd+2GP6cl0UVXuasXAPAF\nAgnaD3t0A2FNS1p9uATTkgAAgQTtB4uuAoAJCCRoVx4RMT3+Xx7z8NXqvf++lUQISTlT3iPpJL4n\nAdgyBBK0N1OjwPE9CcCGIZCAAx7hMew17pZU72U2p8AABwCbhUACDhgaBb6NHuNLEoDNQiABN1z6\nDGVvTsHuuMPMJADbhEACzuhsTsHuuMO6qwA2CIEEnDGwOQXWXQWwYQgk4JJHRIzOrrJYdxXAZiGQ\ngGP6HXf0ALvKAtgaBBJwTL/jjll3FSPuAGwKAgm4p7Or7JCGQmZXWUxLArAdCCTgBZ1dZdnrrnJU\nIwBobwgk4AX9dVcxLQnA1iCQgC90NkzCtCQAW4NAAr7AtCQAG4dAAh7RWU+IPS0JHXcAoodAAn4x\ntp5Q9rV76LgDEDcEEvCLsY47RVUDOu4AxA2BBLyD9YQAbBMCCfjI6ELgWE8IQLwQSMBHJjvuMLoB\nQJwQSMBT+h13dD2h7Gv3sJ4QgCghkIC/2F+SCCHs9YTQcQcgPggk4C/99YSW3EfHHYBoIZCA13QW\nAo+uy0HHHYBYIZCA73QWAmfWE0LHHYDIIJCA73RG3A1pKPxrWhI2pwAQEQQSCIBHeIzhaUlYCBxA\nRBBIIABSb/+urz4+ugELgQOIDgIJhEF/IXA6ugEj7gBEA4EEgqGznhAzLQkLgQOIAwIJBEN/PSGa\nSVgIHEAcEEggJDrTkoY0FDIdd5iWBCB0CCQQGJ1pSez1hDiqEQBYBgIJBEZ/PSFscw4gDggkEB6X\nPkPZHXfY5hxAHBBIIDwmdkvClrIAwoVAAkEyNi0Ji64CCBcCCYTKxG5JXFQHANoKgQRChd2SAEQG\ngQQCpjO6IbouB4uuAggXAgkETOrtb2y3JKzdACA4CCQQNqm3v7FFVzHiDkBYEEggeMZGN6ScLceW\nsgACgkACwcPoBgBxQCCBGOgsusqMbsDaDQACgkACkdBbdHUbPcboBgChQCCBSOisJ8TemQKjGwAE\nAYEE4uERHmNwS1mMbgAQBAQSiIfJLWUxugGA7xBIICrGtpTF6AYA/kMggdgY21IWoxsAeA6BBGKD\naUkAAoVAAhEytugqOu4A+AyBBCJkYktZdNwB8BYCCcTJ2JaymJYEwFsIJBAtLLoKICwIJBAt/dEN\nmJYEwGcIJBAz3UVX/+y4y752Dx13AHyDQAKRY09LIui4A+AxBBKIHDruAIQCgQTih447AEFAIIFN\nMNFxh6myADyBQAKbYLLjDlNlAXgBgQS2wljHHabKAvAEAglsCDruAPgMgQQ2xFjHXWWtqnsnZ+7q\nBQCEIJDA1hjsuHN0sOvuiUAC4BgCCWyOTsddmirlwspnuKoMADAQSGBzdDru7KrKHL76B4f1AQDK\ngesKAHDAIyKmJntX/cWT9KHy4qnSVdO5rRLwVpfFn0u9/bmuhU1AIIGN6rL485KERz11qts3Vbdv\nclsf4K1b2970e/c/XNfCJqDLDmyUzq6yAMbU/3ZMVaHguhY2AYEEtssjPAaZBM3ymv2etGt3rmth\nE3jaZXf+/PmXX3756NGjrq6utKSpqWnx4sUajYY5x97eftu2bRxVEMSAfknyCI9RXjzFdV2Avzwi\nYriugq3gYyDdunXrww8/rKqq0mq1TOG1a9eys7P/7//+z9n50XwRe3t8vQMLkHr74ydrAD7gVyB9\n++23aWlpRUVF7G9C1OXLl93d3d977z07OztO6gYAAFbFr0Dq169ffHw8IeT06dO7d+9mP3X58uWQ\nkBCkEQCAWPErkHr37t27d29CSENDg34gderU6bXXXjt37pyHh8eoUaMWL17s4uLCUU0BAMDCBPMz\nzOXLl7Ozs729vZcuXRoWFvbll1/OmzdPv2cPAAAEiptvSLW1tRkZGcxDLy+vMWPGmDhfrVYvWrRo\n4MCBffv2JYRMmzatX79+b7/99uHDh8eNG2f16gIAgPW1KZCqq6ujo6Pnz58fFxen/2xWVlZGRkZR\nUZFMJgsKClqwYEFAQAB9qqqqau3atcyZAwcONB1IEolk9uzZ7JK//e1vq1evvnDhAgIJAEAc2hRI\n27dvLysrq6ur038qKSkpLS3NyckpNDS0pqZm165d+/bt27x584gRIwghAQEB586dM/9Gd+7cKS4u\nDgsLk0gktEQikTg4OBi8NQAACFGLf0NSq9UKheLYsWPLly9PTk42eE5eXl5aWlq3bt0OHjyYnp6e\nmZm5ceNGlUqVmJjY0NDQiloqFIpZs2adOvXX7MXffvutrq4uJCSkFa8GAAA81OJAKi4ujoqKio+P\nz8zMNHbOnj17CCErV6709fWlJVFRUePGjauoqDhx4kQratm/f/+goKDExMQjR46Ul5cfP3585cqV\nAQEB0dHRrXg1AADgoRZ32cnl8k2bNtHj8+fPb9++Xf+cvLw8iUQyatQodmFkZOSBAwdOnz4dGRnZ\n0ps6Ojpu27ZtzZo1S5Ys0Wq1Uql05MiRq1evdnR0NH1hcHAwc3zlypWW3hcAQATYn4R81uJAkslk\nUVFRjy52MHC5SqWqrKz08/Nj1vihevbsSQgpLS015y4zZsyYMWMGu8THx2fr1q1KpbKystLX19fg\nrfUhhAAA2J+EfA4nyw/7rqmp0Wg0HTp00CmnJdXV1W15cRcXF2aoHgAAiInlJ8aqVCpi6MsTLaHP\nAgAA6LB8INGR2frBQ0uYcdsAAABslg8kur6cUqnUKaclWH0OAAAMsnwgyWQyd3f3mzdvqtVqdrlC\noSCEyOVyi98RAABEwCqLq4aEhDQ2Nl68eJFdWFBQQAgJDQ21xh0BAEDorBJIY8eOJYSsW7eOKSkv\nL//mm28kEsno0aOtcUcAABA6q6z2HRsbm56enp+fP3Xq1AkTJty6devAgQP19fV///vffXx8rHFH\nAAAQOqsEklQqTU1NXbNmTVZWFu24c3NzW7FixYIFC6xxOwAAEAE7rVZrvVdXqVQKhcLNzU0ul7f/\n7uPBwcFYqQEAgI3PH4zW3aBPKpX26tXLqrcAAABxEMwW5gAAIG4IJAAA4AUEEgAA8AICSZD4vIB8\nW4i1XUS8TRNru4iom8ZbCCQAAOAFBBIAAPACAgkAAHjBuhNjuYUuYAAAfbydGCvmQAIAAAFBlx0A\nAPACAgkAAHgBgQQAALyAQAIAAF5AIAEAAC8gkAAAgBcQSAAAwAsIJAAA4AUEEgAA8AICCQAAeAGB\nBAAAvODAdQUsLysrKyMjo6ioSCaTBQUFLViwICAggOtKtcb3339/7Nixq1evenh49O3bd+HChV26\ndNE5R9CNVSqVs2bN8vT0TE5O1nlKiO36448/srKyTp8+XVJS4u/vP3369LFjx+qcI8R2abXaw4cP\nf/fdd9evX+/YsWNgYODcuXN79uypc5pQmlZdXR0dHT1//vy4uDj9Z5ttBZ+babppzX6e8KFpYltc\nNSkpKS0tzcnJKTQ0tKampri42MnJafPmzSNGjOC6ai2gVquXLFly9OhRmUzWp0+fioqK69evOzs7\n79ixY9CgQcxpQm/su+++u2vXrj59+uzdu5ddLsR23bhx46WXXrpz506XLl06d+586dIlQkhCQsKy\nZcuYc4TYLkLIm2++uW/fPjc3t759+16/fr2iosLBweGLL754/vnnmXME1LT169cnJye//vrrL7/8\nss5TzbaC58001jRzPk/40jStiJw+fTooKGjMmDGlpaW05McffwwJCRk5cqRSqeS2bi3y9ddf03+h\nMNXetWtXUFDQc88919jYSEuE3tijR4+Ghob2798/OjqaXS7Edj148CAiIqJv3765ubkajUar1RYX\nFw8aNKh3795MK4TYLq1We/To0aCgoClTpty7d0+r1Wo0mszMzKCgoKFDhzY1NdFz+N+0pqamkpKS\n7Ozs1157LSgoKCgoaOvWrTrnNNsKfjbTnKY1+3nCn6aJ6jekPXv2EEJWrlzp6+tLS6KiosaNG1dR\nUXHixAlOq9YyKSkpDg4OH330kbOzMy2ZMWPGc889d+vWratXr9ISQTf2zp07iYmJixYt8vT01HlK\niO06fPhwWVnZG2+8MXz4cDs7O0JIjx49FixY0KlTpwsXLtBzhNguQsgvv/xCCFmwYEHHjh0JIXZ2\nduPHjw8MDLx7925xcTE9h/9NKy4ujoqKio+Pz8zMNHZOs63gZzPNaVqznyf8aZqoAikvL08ikYwa\nNYpdGBkZSQg5ffo0R5VqMa1W+8cffwQHB3t7e7PLe/ToQQgpLS2lDwXd2HfeeUculy9atEj/KSG2\nKyMjQyKRTJ06lV34yiuvnDx5cty4cfShENtFCNH/FwMhRK1W29vbP/HEE/Qh/5sml8s3/WnhwoUG\nz2m2FfxsZrNNM+fzhD9NE8+gBpVKVVlZ6efnx/wrgKK/vjKf4/ynVqtXr16tP36hqKiIENK9e3ci\n8Mamp6efOHFi7969EolE5ymBtis/P79///4ymex///vfuXPnysrKAgMDhw4dynxkC7RdhJAxY8as\nX78+JSVl9OjRbm5uhJDDhw+XlJQ8++yznTp1IgJpmkwmi4qKoscODgY+9JptBW+b2WzTmv084VXT\nxBNINTU1Go2mQ4cOOuW0pLq6motKtYaDg8O0adN0Ck+ePHnq1KnAwMDAwEAi5MYqFIpPP/102bJl\nvXr10n9WiO2qra1tbGyUy+Xbt29ft24dU+7h4fHxxx+PGTOGCLNdVEBAQEpKytKlS0ePHh0WFlZc\nXFxSUjJ8+PAvvviCniDcprE12wrhNrPZz5P79+/zp2ni6bJTqVTE0L8RaAl9VqCOHDmSkJDg5OT0\nwQcf0G8VAm2sWq1+4403evfuPX/+fIMnCLFdd+/eJYScPn36X//617vvvpubm5ubm5uYmNjQ0PD6\n66+XlZURYbaLUVtbq9Fo7t+/n52dXVJSQgipr6+/d+8efVbQTWM02wpxNJPS+TzhVdPEE0jsT2o2\nWqLfOyQINTU1//znP1999VU3N7cdO3YMHDiQlgu0sZs3by4qKvrkk0/s7Q2/8YTYLlq3qqqqN954\nIy4uzsvLy8vLa/bs2cuWLXv48OH27duJMNtFffXVVwkJCV5eXmlpaefPnz9z5kxiYuKlS5emTJly\n5coVIuSmsTXbCnE00+DnCa+aJp5AcnFxIYQolUqdclpCnxWWn3766YUXXti7d++ECRMyMjLYM5CE\n2NgLFy5s3br1lVde8fT0rP2TRqPRaDS1tbX19fVEmO2iw88IIdHR0ezyiRMnEkIKCwuJMNtF7dmz\nRyqVJicnP/PMMw4ODh4eHrNnz3799deVSuWRI0eIkJvG1mwrRNBMY58nvGqaeAJJJpO5u7vfvHlT\nrVazyxUKBSFELpdzU63W2rlz56JFixwdHVNTUzds2EB/QGYIsbG///67Wq3esGHDIJaKiorCwsJB\ngwa99NJLRJjt6tixo1QqdXV1pb/5M+ifrKqqigizXYSQe/fuFRYWBgUFMaOBKboCBR1/JdCm6Wi2\nFUJvponPE141TTyBRAgJCQlpbGy8ePEiu7CgoIAQEhoaylGlWuOnn3765JNPBg0atH///meffdbg\nOYJrbJ8+fRL0yGQyLy+vhISEmJgYeprg2iWVSp999tn6+nr6cxGDjmJiFl8RXLsIIW5ubhKJpLa2\nVqf8/v375M/EJcJsmr5mWyHcZjb7ecKfpokqkOg/3NgjncrLy7/55huJRDJ69Gju6tUyarX6k08+\ncXd337Ztm0wmM3aa4Brbv3//ZXo8PDy8vb2XLVs2c+ZMeprg2kX+rPPGjRu1fy7EpdFoNm/eTP7s\nuCPCbJejo2NoaOiNGzcOHTrEFGq12q1btxJCwsLCaIkQm6av2VYItJnmfJ7wp2niGfZNCImNjU1P\nT8/Pz586deqECRNu3bp14MCB+vr6v//97z4+PlzXzlxFRUXXr1/v1q3b559/rv/svHnz/Pz8iFga\nq0+I7Zo2bVp2dva+ffsqKiro/9s//vjjmTNnwsLCJk2aRM8RYrsIIe+//35sbOzSpUv79u07fvx4\nR0fHzMzMc+fO9evXj/k3hECbpqPZVgi0meZ8nvCnaWJbXPXOnTtr1qzJysqi/aFubm6vvPLKggUL\nhDIMhhDy7bffvvPOO8aeTU9Pf/rpp+mxCBobERHRqVMnncVVhdiuhw8ffvzxx8ePH6cdd56enhMn\nTnzzzTelUilzjhDbRQj57bff1q9fn5eXRx86OjrGxsa++uqr7JkrAmra0aNHExISDC6u2mwreN5M\ng00z8/OEJ00TWyBRKpVKoVC4ubnJ5XK6tpiIibWxAm3XH3/8YWdnZ+KnYIG2S6lUlpaWOjk5+fr6\nGvuQEmjTdDTbCnE00yDOmybOQAIAAMER1aAGAAAQLgQSAADwAgIJAAB4AYEEAAC8gEACAABeQCAB\nAAAvIJAAAIAXEEgAAMALCCQAAOAFBBIAAPACAgkAAHgBgQQAALyAQAIAAF5AIAEAAC8gkAAAgBcQ\nSADtobGxsaKiAtuPAZiAQAKwiuTk5JEjR+bm5hYUFMycOXPQoEGjRo0aOHDgO++8U19fz3XtAPgI\ngQRgFbW1tbdu3crLy1u4cOGdO3eioqIGDRqkVCq//fbbt99+m+vaAfCRA9cVABCz5OTkRYsWLV26\n1N7enhDy66+/vvjii4cOHaqsrPT29ua6dgD8gm9IAFbUu3fvZcuW0TQihDz11FOhoaEajebGjRvc\nVgyAhxBIAFY0YsQIOzs7dklgYCAh5P79+xzVCIC/EEgAVtSlSxeuqwAgGAgkACvS+XoEACYgkAAA\ngBcQSAAAwAsIJAAA4AUEEgAA8AICCQAAeMEOqz0CAAAf4BsSAADwAgIJAAB4AYEEAAC8gEACAABe\n+P8BP8AqSzJYSCkAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "x = b./diag(A);\n", "Anorm = @(x) sqrt(x'*A*x);\n", "Aerr = Anorm(x);\n", "for n = 1:100\n", " [xCG,~] = pcg(A,b,1e-14,n);\n", " Aerr(n+1) = Anorm(x-xCG);\n", "end\n", "semilogy(resnorm,'.-')\n", "hold on\n", "semilogy(Aerr,'.-')\n", "ylim([1e-16 1])\n", "xlabel('n')\n", "title('Convergence of CG')\n", "legend('||r_n||_2','||\\epsilon_n||_A')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Effect of condition number\n", "\n", "We already know that a large $\\kappa(A)$ creates error in the solution of $Ax=b$. But it carries another penalty in Krylov methods: slower convergence. We have essentially the same result for CG as for MINRES, but using the $A$-norm of the errors:\n", "\n", "$$\\frac{\\|\\epsilon_n\\|}{\\|\\epsilon_0\\|} \\le 2 \\left( \\frac{\\sqrt{\\kappa}-1}{\\sqrt{\\kappa}+1} \\right)^n,$$\n", "\n", "where $\\kappa=\\kappa_2(A)$ equals the ratio of max to min eigenvalues." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAABcSAAAXEgFnn9JSAAAA\nB3RJTUUH4AsHFRcsHSExZAAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJ\nbmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAwNy1Ob3YtMjAxNiAxNjoyMzo0NEGrHOIAACAA\nSURBVHic7d15XFTl/gfwB4ZhHUbFMIcdURQQTZHKJRc00et2cUG8uYK5kLmkZKX91JL06rVbtmgm\nF6RbImouF0lNBb2KAoZXCzVlDRBFRZYAZWD4/fHYcZqNMzDLmTOf9x++Zp7zzJznPMj58pzvc55j\n0dLSQgAAAIzN0tgNAAAAIAQBCQAAOAIBCQAAOAEBCQAAOAEBCQAAOAEBCQAAOAEBCQAAOMHK2A0w\nml9++SU9PV1znQ4dOsyaNYt5e/ny5V9//bWuri4gIGDw4MGtlpsEU2z8tWvXzp07RwhxcnL629/+\nRgg5fvx4Xl4eIaR3797Dhw/X8Fn2NdvMALsA4KcWc/Xll1+22jne3t608v379+XPLNOnT9dcrls7\nd+7cunXr1q1bf/75Zx1+rWEarw+ffvopbbO/vz8tmTZtGi1ZsGABU01lv6msqVsG2AUAL5nvCEkr\nq1atUjmcUleuW7GxsSUlJYSQzp079+7dW1dfa5jGG5Ge+g0A9AQBiRBCJk2a1LlzZ+VyZ2dn+uI/\n//kPfREdHb1+/Xpra2vN5SbBpBuv4P3331+wYAEhxNXVVVc1DdAYAJCHgEQIIf/3f//Xv39/dVvr\n6+srKyvp6+nTpzNRSl25vObm5vz8fIlE4ujoqLkNhYWFMpnM09PTykpnP5Ti4uKWlhZPT08LCwuF\nTWwa3552ati1soqKitraWm9vb0tLtbNsysvLW1paXFxcVG4NDAwMDAxsdUdsara/5Wwaw+aQy8rK\nLCws1B1yq3T+I1CJ/f9bw7RHHvvfPkJIWVlZXV2dr6+vQiHLH4FW+yLsDpBl32q7a04z5vVCo5LP\nIf3000/qqrm6usr/d3R2dnZ1df3ss8/UlTMfvHTp0iuvvGJvb08IsbS0DAwM3Lhxo0wmU/j+Gzdu\nhIaGdujQgX6PQCCYNGlSdnY23Tp48GBXV1eBQEC3durUydXVddKkSZoP7ebNm2PHjmW+09HR8dVX\nX5XPo7TaeGWa28l+1y0tLVOmTHF1dXV1dT18+PCJEyd69OhBK9vZ2c2bN+/333+XryyTydatW+fp\n6UnreHh47N27VzmHtGjRIvqdq1ev1txvCjX10XLlXWh7yB9++KGbmxut4+bmlpycHBcXR7+h1byU\nzn8EKrH8/6DD9kRERNA68+fPl/9sbGwsLZ83bx5T2Opvn/wec3JygoODCSESiaQNPwI2v+ladTj7\nvmV5kjEhCEitBCSiykcffaSunH5q8+bNzNlQ3tSpU+vr65kv//bbb+l/JgVWVlY//PBDS0uLj4+P\n8taBAwdqOK49e/bY2toqf0ooFO7YsUPzQan7zlbbyX7XLS0tr7zyCi0PDw9Xvk44efJkpmZzc/Nf\n//pXhQoWFhbMcFbdpAYN/aZyxoFuW668C/aH3NLSEh4erlDBysoqJCSEvp4yZUo7f/ratqfN/x90\n254dO3bQt46Ojk+ePGE+GxAQQMsTEhJoCZvfPmaPH3zwATOwYAIS+x8By9909h3Ovm9Z7tq0ICAR\nQoirq2s3VfLz81NSUg4cOMDU/Mc//pGSkpKXl6euvKWl5cKFC/SihI2NzZo1aw4dOrRixQrmmtgH\nH3xAG3Dnzh2RSEQLhw0blpiYuH379qCgIFoiEomqq6vT0tJSUlKee+45Wrhs2bKUlJQLFy6oO6ji\n4mIHBwdaOSQkJCEhITExMTQ0lJbY2NjcvHmzpaVFQ+OVsWkn+123yP1yEkIcHR2jo6OXLl3K7MLS\n0vK3336jNf/973/L/5olJSV99dVXL7zwAlOoLiBp6DflaKHzlmsISK0eMpPYI4RMmzbtwIEDO3bs\n8PPzYwo1BCR9/Aja/P9B5+2pqKhgrlwdP36cfpBOryeEWFtbV1VVtbD+7WP2yLTQ09MzODhYqx8B\ny32x73D2fct+16YFAUkT+gvz+++/MyUXL15kvkFdeb9+/Wjht99+yxQeOXKEFjo7O9O/X+bPn09L\n+vXr9/jxY1qtvLzczs6Olh86dIgWuru705J//etfmg9q5syZtOagQYMaGxtpYVNT06uvvkrLJ06c\nqLnxyli2k/2umV/Ozp0705NIS0tLWloa0x7mL0HmLCA/Jb2qqorpEM3TvlX2m3JNnbdcQ0Bq9ZCZ\nmsxOW1paHjx40LFjR1quISDp40egjP3/W523Z/To0fTt4sWLaZ1t27bRkgkTJtASlr998hFiwoQJ\n9+/fZyqz/xGw3Bf7A2Tft+x3bVqwUgMhhIhEog6qtCGhWl1dfeXKFUKIpaWlu7v7pT907ty5U6dO\nhJD79++fPn2aEHL27Fn6kTfeeMPGxoa+7tq1K/3rPiUlhbkWwd7Fixfpi+joaKFQSF8LBIIlS5Yo\nVGCPZTvbsOu//OUvzIXy4OBg5vrDvXv3CCE1NTU3btygJcuWLWM+pXC3cvvpvOUatPrB//3vf/TF\nwoULmU917tx5ypQpHDkQ9v9vdd6eiIgI+vbo0aMtLS2EEOYUPH36dKLNbx/jueee27dvHzOYJqx/\nBG3YV6sHyLJv27Zrk4BZdoQQcvbsWQ2z7LTyyy+/0BcymWzo0KEq65SXlzc2NhYUFNC38tegCCEv\nvfRS23ZdX19fWFhIX8tfYSCE+Pv70xf379+/f/8++zl1LNvZtl1LJBLmtYODg5ubW3FxMd0pIYS5\nGiP/JVTPnj1Ztr9V+mi5Bpo/eP/+/draWrqVmcdBqcyKGf5A2P+/1Ud7wsLCFi1a1NjYWFZW9tNP\nP3l7e1+4cIEQYmtrO3HiRML6t0/+7dixY5nxB9HmR9CGfWk+QPZ927ZdmwQEJB178uQJfWFpadm9\ne3eVdWQyWV1dXXNzM33L/PHYTk1NTTKZjL5WyHbKv5VKpey/k2U727Zr5m9ASmE8Kv9WYaIwm3nD\nLOmj5Rpo/qBYLLawsKB/+1dXV8tvqqur0/zNhjkQ9v9v9dGejh07hoaG0hzPkSNHunfvThszduxY\nOjGB5W+f/FuFe8XY/wjasC/NB8i+b9u2a5OAgKRjzPWKlpaWn376iUlRKpNIJPSvmFu3bvXp04cp\n37Nnz7Vr1wgho0aNGjt2LPtdi8Vid3d3ujZBXl5e3759mU3MaKNjx45a3dfSqVMnlu3U+a7lf9Ny\nc3MHDhzIvM3Pz2f/PZrpo9PazMbGpmvXrrS3T506NWjQIGbTjz/+qPmzhjkQ9v8f9NSeiIgI+YBE\nC+n1OqLNbx9D4dTP/kfQhn1pxr5vdb5r7kAOiRBCHj9+3KCGtl/1/PPP09sXWlpajh49ypTLZLJl\ny5ZNnjx58uTJd+/eJYQMGDCAbvrXv/7FVHv06NEbb7zx8ccff/zxxy1K87MrKio0752ZkBMXFydf\nvnv3boUK7LFsp853LRKJmGs7n332GVP+5MmT7777Tquv0txv+ui0NqOXnggh27ZtO3/+PCGkubl5\n3bp1WVlZrX7WMAfC/v+tPtozceJEeoXt559/PnbsGCHE3t5+/PjxdCv73z7Nu6AvNP8IdLIvBSz7\nVh+75gojTKRoq6ysrAULFowcOXLatGmffvopMwulbdjMsiOE1NTUaDvLjpmsLBaLFy5ceOrUqQMH\nDjD/y0eMGEGr/e9//2Omsc6dO/fcuXMJCQnMOKBr164NDQ20JnNe9vT0XLt2LXO/hbJffvmF+Ytv\n5syZp0+fTktLi4qKoiWWlpZZWVmaG6+MZTvZ75qZcfT+++/L78jb25uW79q1i5bs2bOHaWRUVNS5\nc+dSU1OZOVqktVl2KvtNuabOW65hll2rh5yfny9/W76bmxudl8zcvKJhlp0+fgTK2P+/1VN7pk6d\nSuRMmzZN/iMsf/vU7VGrHwHLfbE/QPZ9y37XpsVkAtK5c+d69uw5c+bMb775ZsuWLYGBgZGRke25\nJ1l/AalF1Y11VJcuXa5fv85U27x5s8p0iFAoPH/+PFMtMjJSfqvmG2O3bt2qLhOwbt06No1XxrKd\nLHfN/uwjlUpHjRql/G3MJQvNAUllv6msqduWtycgtbS0nD59WmFlxVGjRsXGxtLXmm+M1fmPQCWW\n/x/01B75W+gIIfv371doHpvfPg0BqUWbHwHL33T2B8i+b1nu2rSYzCW7LVu2+Pv779mzZ+bMmTEx\nMbGxsefPn8/IyDB2u1Tbt29ffHy8fL7U0tIyLCzs3Llz8jOOVq9enZ6e3qdPH/n/giNHjrx48aL8\no4m2bdsWFhbGcvHTVatWnTt37oUXXmC+08LConfv3mfOnFm/fn3bDodlO3W+a3p3+tKlS8ViMS0R\nCoXLli1j7j7RjH2/6aPT2iwkJOSXX37597//vXz58rfffvv7778/duwYk+6WnxKmzDAHwvL/g57a\n85e//IXJmohEonHjxilUYPnbpwH7H0H796WAfd/qfNdc8HQ+Ccc9ePBg8ODB77zzzrx582iJVCrt\n37//zJkzV69ebdy2afbgwYPc3FwHB4du3bo5OTmpq1ZbW/vzzz8LhcIePXow998paGpqqqysbG5u\ndnJyUpiro1JdXd3169dlMllAQICu0p5s2qmPXctksuvXr9fX1/ft25fNscvTqt/00WlaycjIKCsr\nI4R4eXnRBdaooUOH/ve//yWEvPfee8yf6hoY5kBY/n8wWHsUsPztU9C2H0Hb9qUB+77V+a6NydhD\nNFZycnJ8fX3PnDkjXxgaGvrmm28aq0kA+rB582b6iykSiW7dukULs7KyaGrBwsLi2rVrxm0h7+FH\nYESmccmuvr6eyK06Rdnb27d6cwaAaXnttdfovZO///67n59fUFBQcHDw4MGDm5qaCCHz589n+ZQN\naDP8CIzINAJSS0sL8698YYspXG8EYM/Nze3cuXNDhgwhhDQ3N+fk5Fy+fFkqldrb2+/YsWPXrl3G\nbiD/4UdgRAa9Mba6ujosLCwyMpJZdVHeqVOn6JrTIpHI19c3KiqKWbqDZhEV7gpqaGiQX4cDgB+6\nd+/+3//+9/r16zk5OeXl5SKRyNPTc/jw4SqfSgD6gB+BsRg0IO3evZs+llF5U2xsbGJioo2Njb+/\nf01NTXJy8pEjR7744gv6dwoNPEVFRUz95ubm0tJShWknulV1p6aji1h/3w+ggb+/v8IKfmBg+BEY\nnt4DUnNzc0lJSXFx8eHDh1NTU1XWyczMTExM9PDwSEhIoLMYT5w4sWLFijVr1pw4ccLW1tbFxUUi\nkfz4449z586lH0lPT5dKpcyNzRpc+z6xy50c+4BB4hHhTQ8u1WTMsLR3E9i5Wdq7EUKsnnuJvqZv\n5f3zL89uL5ePTMxr5kUnlw6t1kFsAwDQTO8BqaCggFnYQ52DBw8SQlatWsXMqQ8NDR0zZsyxY8cu\nXLgwcuRIQsjs2bP//ve/x8bGTps27datW9u2bXN3dx8xYoSGry2qfDwv6UZUxp4XH9+oSU9+mLxN\nNNBT0JHI6ktl9aXkISGEPCl5dpMdE6gs7V0V4lPVnRqVr9tAc6xqNbYhsAEAX+k9IEkkku3bt9PX\n165dY5axkpeZmSkQCIYNGyZfOHLkyGPHjl26dIkGpDlz5lRXV8fHxycmJhJCAgMDN23apPLpyFRR\n5WPfNf+R2nf+6vHTZ+pI75c05N4RDVZ7yPKBihBCyFItjpM1Jp61M7DVtdQqv1bxQqZYwryWLwEA\n8/Hrr78auwmq6T0giUQi5qHF8itEMaRSaUVFhZubm0J0oU8fKS0tpW8FAsGKFSuWLl1aVlbWqVMn\nutq8BvOSbkjtO29++KcpMcKubB9bUFPF9YGIg4WjytdtoDwg02rQRjg/buvZsydnfwMNBp1A0AmE\nEJ0+TkznjP/4iZqaGplMxjxFkUFLFB5JIhAIPDw82HxtfISfd2zGZx0mZ9r4Ta479+LjG4SQ+qvN\n9Vebhc9bkj+Ck6WDhbCrirnv4o41895KeNrCR0/PtjVVz877tVVi8ue4xbyueeSoUMJxuCAJAFxg\n/IBEn9ClPHiiJVo9TU6el5Nt16uJZX1nHxK9ckj0imvTg8m/n1tS/T0hpPn3ZkLIY7lH6ghEFoQQ\nu0CPzuErmx5kyupLCSFicoluFXds1zmaCUuWdm61VY50AkXNI7GlvZvA3o38EQMe3XkaepXDQzuD\nhCHp6oIkm1kkKjcRxDYAk2X8gEQfH6kceGiJwrMmtSIuubgoZt36k4WEkDKr5z7rOLmy78QvA6tq\n0pPrc/+0Kmvz7y2EkN8vFj/J22oXMLDDiKV2AQPJH1ml5oZSQggNVM0NpTRcadGMZ/Hsuuh5Qkgm\nIcT5OUIIoVMnBIHPZlII7Nws7QOUp/xRKocyyi/kY1tWVtaLL76oqzhhMLodtE2wfY3OmeTlBUkA\n3jB+QFJ50ytTonlt41bNCZak51el5z+ib78tturRM3jdhnBpRcnD5G016ckK9aX3S6TpJTXpyUJn\n987hK8Ujwi3t3Wgf2bg/ewoLE6hk9aWy+rK2BSr6PU//ffin8qeBys7N0t7N6rmX6Gur515Wed7U\nbB6ZprJc8zhMYdBGNMY/jqu6U+Ng4aiTsaZWFyS5NmhD7oSgEzjP+AFJJBI5OjqWlJQ0NzfLj4fo\nbbDtXIvBy8mWJpOYkoTL5cO6dxru4951ySedw1fWpCc35F5UGDARQqT3S+5+sfxh8ja7gIEdRkyn\nAyYGveym0HdMdGluKG16kElfSx9eakOz5QOVwsR0Qoiw88vkjzuoCCFWz73chl20IbappHnQRnBB\nUpU2X5DEoA34zfgBiRDi5+eXlZWVm5sr/xj5nJwcQkj775T2crJNi+4/4ssc+rao8vG8pOuFawYR\nQoRd3DuHrySENORerE7bx2bApGFHNFrQQKUwnCKE0Mikk0D1pP4AUXMHFdF4q68+6CqwEXYXJHkz\naFP5ug14P0MSzAonAtLo0aOzsrK2bt36zTff0JLy8vLvvvtOIBCEhIS0//uH+3RcP9qbJpPIHzfM\nxkc8e4aVXcBAu4CBdMBUk5YsvV+i8A3MgEk8Ilw8PFzYxZ393mlgsLGfSjQGKvm32mJ5q68hA1Ub\n6CO2sRm0tVqfy4w1Q5JrFySBHzgRkCIiIpKSkrKysiZPnjxu3Lh79+4dO3asvr7+9ddfd3Fx0cku\nFJJJCdnlXp1s14V6y9ehAybx8HB6EU/lgOlh8raHydvsAwaJh4drHjC1SkOgav9MCvmvYpOg4nKg\nagNd/fnfhlkkCvVNJbARXJAEHTlw4EB0dHROTo6bm9anFIM+Mfb06dPR0dFvvfXWwoULFTY9ePDg\ngw8+OHXqFH1OsIODw6JFi6Kiotozy07hJriiyscjduQUVT6mb72cbOMj/If7qH0Uo7SiRN2AiRI6\nu7dhwNQ28gmqds6kUEfdTApdfb+ZM58LkjqEC5L6oNe7g0tLS/v06fPo0aOioiLmcQ3scesR5lKp\ntKioyMHBQSKRyD9Svm2U+z09v4pJJhFCvJxsaTJJU5MqShpyLyrPFGcInd1VTnwwDB3OpFBHhzMp\nQCfaeUHSDAMbwQVJOfoLSDKZbOTIkU5OTt9//z0fApJuqez3DScKmWQSIWRusEQ+maQBHTA9TN6m\nrgIdMNFZEkanwwSVOkacSQE6YVYXJHWFBxck9ReQNm/efPLkyfXr1w8bNgwBSZHKfqczGphkEiFk\n/WhvhWSSBhwfMLVKhwkqdUxrJgXoBC5ItoGxLkjqKSBdvnx57NixP/3002+//fbKK68gIClS1+/a\nJpNUMq0BU6t0dauvOkhQAUu4IKktzUGLKMW26TFhOg9IdXV1/fv3f++99+bMmXP+/HkEJBU0/CHQ\nhmSSSqY+YNJMZYIKgQpMAi5IqpPUsFPnAWn+/PmVlZXff/89IQQBSTXNI9M2J5NUUrcWEcO0Bkya\nYSYFmBWerbOlISBVVlYGBQWtXLlyyZIlCpsOHz68d+/e69evOzo6BgYGrlq1qkePHnTT0aNHZ8+e\nfe7cOWdnZ0JIZmZmWFhYVlZWt27dOnfurFXbzDcgtTOZpBK9jqdyLSLKpAdMrcJMCgDNuHBBUkNA\nevfddzdv3vzRRx+9++678uXLli3bvn27ra1t//79Hz16dPPmTVtb20OHDtFn3b333nubNm1S/rah\nQ4eePXtWq7aZb0AiOkomqcRmwNTqWkS8gZkUALrVnguSm/+7Vv7E2NTUVFBQkJeXl5iYuG/fPkKI\nQkBKS0sLCQnx8fE5deqUl5cXIeTgwYPTp0+XSCS3bt2ys7O7d+9eRUUFUz8nJ2fu3LnHjx/v0aNH\nt27dtDousw5IRFUyKW1xfy8ntU9G1wrLW2vtAwbxcsDUKsykADA8hRPjL7/8EhgYKF9BISDNnj37\nm2++2b9//9SpzxaUmTFjRlJS0uHDhydNmqTw/e3JIXFi6SAjUrnMXVp0P518OV2LqHP4ypq0ZE1r\nERHdrEVkctgvmt7+p3soL5qOQAVACPHw8Ni/fz99nZ2dvWXLFoUKaWlpVlZW48aNky+cNGlSUlJS\nWlqackBqD3MPSERpmbv0/EcbThS2M5mkQDwinM5oUDdgqs/NqM/NaNvirTyjedF0U3m6B4C+FRcX\nE0LaMApRIBaLmaGPUChU2NrY2Hjnzh1vb2+FR9P5+fkRQgoLC4mSIUOGtPnCGwLS02cmySeT/nhm\nkg6SSfLkF29VOVOcDphq0pLtAgbaBwwytwGTZiyf7tHORdMJx57uAaBOSEjI0KFD4+PjCSH37t1b\nu3ZtfHz8jBkzmGcm3L59W+UHGxsbb9++3aFDhy5durS6l6qqKplM5uTkpFBOSyorK9t1DEoQkAj5\nYzqDwjOTdJhMkifs4i7s4i4eEa7u1lrmIUwYMLFh4EXTTffpHsAnJSUlBQUFGzduJIRcuHAhOzt7\n7dq1Z86ckX9eT69evWQymcqP+/r6Llq0aMeOHa3uqLGxkRBiZaUYKWgJ3apDCEhP6TWZpJJWAya+\nzhTXH/kElT4eP0/+HKjk90uQoAI9O3/+PCFkyJAhcXFxYrF4+fLlhJD8/Hz5OqmpqSo/GxUVFRcX\n5+7O6s9c+rwFqVSqUE5L2vM0BpUQkJ4xQDJJGfsBE59urTUijjx+niBQQTukpqba2NjMnDlzwoQJ\n48ePV1mH3iSkzMHBQd0mlZUJIXV1dQrltIRu1SEEpGcMlkxSCQMmIzLYTAqVCSqCmRSgpYyMjJCQ\nkOjo6NLS0rFjx9rZ2SUlJbEc9GhFLBZ36NChoKCgqalJ/sLdrVu3CCEeHh663R0C0p8YMpmkEgZM\nnIKZFMBB5eXlBQUFH330ER0bLViwwN/f/5133vn2228JIcXFxXTqnYuLi8rZbg8ePJBIJLNnz/77\n3//OZnf9+vVLT0/Pycl58cUXmcKMjAy6SSdHxEBAUmT4ZJJKGDBxGWZSgBGdO3eOEDJkyBD61tLS\n0tLSkr6+ePHib7/9RgNSZGSkyoC0c+fOyMhI+eii2eTJk9PT01evXp2WlkZLSkpKPv/8c4FAMHHi\nxHYeiwIEJBXWhXobPpmkkvyASeVaRPIDJvNZi4izMJMCDOD48eNeXl6urq5MSV1dXXBwMCEkJSWF\nWWSBzsFTduDAgdjYWPa7W7hw4c6dO9PT04OCgmbMmFFWVrZ37976+vrVq1fjkp2BGDGZpJKwi3vX\nJZ9ouLVWer/k7hfL6Uxxs12LiLMwkwJ06MqVK2PHjpUviYmJSU1NtbW17datm0gk0u3urK2tT58+\n/eabbx46dCgnJ4cQ4ujouGnTppiYGN3uiGAtOw2KKh97xz67UKbbZe7aqSH3YnXaPizeykt4ugdo\nlpeX5+HhYW1tLV9YXFwslUq7d+/e6sfbfGKkd9SKRCIPDw8LC4s2fEOrEJA0UXhm0nCfToZPJmnA\ncvFW3FrLD3i6B+iEnh5hrhMISK0Y8eUV3T4zSedafWotIUQ8PBxrEfEVnu4BWkFAMg6d9Lv+npmk\ncxgwgTw83QNUQkAyDl31O5eTScpaHTDx+6m1oBkSVICAZBw67HeOJ5NUUndrLQO31gJDIUFF9BOo\nkKDiAgQk49Btv3M/maQSBkzQHoacSYEElWEgIBmHbvvdhJJJKmHABDqk75kU8gkqBCrdQkAyDp33\nu2klk1TCgAn0CjMpuA8ByTj00e8J2eXzkm4wb00imaSSurWIGBgwga5gJgWnICAZh576fV7SjYTs\ncuatqSSTVKLX8RpyL2LABIaHW32NAgHJOPTU7wrJJEJIWnR/E0omqYS1iIA7MJNCrxCQjEN//c6D\nZJJKuLUWuMyQMymYBBX/AhUCknHotd8VkkleTraFawbpaV8GRic+1OdmaBgw2QcMEg8Px4AJuAAz\nKbSCgGQc+u53PiWTVMKACUwXZlKog4BkHPrud14mk5RhpjjwCWZSICAZhwH6na/JJJVway3wmPks\nmo6AZByG6XceJ5NUwoAJzIpRElR6DVQISMZhsH7nfTJJJQyYwGzJJ6hMbiYFApJxGKzfzSSZpBIG\nTAAMk5hJoacTY1JS0u7du/Py8iQSyfTp06OjoxUess4GApJumFUySSWsRQSgDqdmUujjxPjJJ5+8\n9dZbkZGR/fr1O3/+/P79++fNm/f1119r+z0ISDpjbskklbAWEQB7RplJETB4jm5PjE1NTZ06dVq8\nePGWLVtoSXh4+KFDh2pqauzs7LT6KisdNsvMzQ2WnM2vYpJJRZWP5yXdiI/wM26rDEzYxZ2OgdQN\nmKT3S6TpJTXpyViLCIAOYuhZ2MZ9KlOuw5kU9KvIw2clkk7CdjZbQUlJiYeHR1hYGFPy4osvHjp0\n6MmTJ9oGJIyQdEk5mRQf4Tc3WGLINnAKy1tr7QMGYcAEoJnKBFUbAtXAGGu9nhhv3LgxduzYgQMH\n7t27V9vPIiDpGJJJKrW6eCvWIgJoG21nUugvIJ05c+a11167e/duUFDQxYsXhUKth2IISLqHZJI6\nWIsIwGBUzqQghAyIzNHTifHu3btZWVm3b9/etm1bjx49Dh065OTkpNU3ICDpnpMmDgAAIABJREFU\nhcKdSXODJeaWTNKA5Uxx+4BBGDAB6JyGE2NlZWVQUNDKlSuXLFmisOnw4cN79+69fv26o6NjYGDg\nqlWrevTooW4XV69efeGFF+Lj4+fOnatV2yy1qg0srRvtLX+ZLiG7XD4+mTlhF3fxiHC3DQe8v8xU\nOQtcer+kJj357hfLCxe/9DB5m7RC9XAKAHRr69atRUVFtbW1CuXLli0LCws7evSoWCyuqqr6+uuv\n+/bte+LECbr18OHDY8aMqaqqYur37dvX2dn5woUL2jYAAUkvaOpIvmTDyUL5yQ5A/piS5/1lZtc3\nPrEPUHFVU3q/5GHyttJ1U+9+vrwh96LhWwjAe01NTbdu3UpNTY2IiNi8ebNyhbS0tO3bt/v4+Ny4\ncePChQvXr1/fv39/Y2Pj/PnzGxoaCCFCofDEiRMHDx5kPlJaWvrgwYPu3btr2xhM+9YXLyfb+Ag/\nJplEJ+AhmaRM2MWdjpnUrUUkP1Mct9YC6NbNmzcDAwM1VIiPjyeEbN682cvLi5ZMmTJl2rRpSUlJ\nJ0+enDRp0qhRo/z9/d9///379+9PnTq1sLBw1apVTk5Of/vb37RtDEZIejQ3WLJ+9LNF7eidSUZs\nD8exHDAVLn4JAyYAXfHw8Nj/h7ffflu5QlpampWV1bhx4+QLJ02aRDcRQmxsbP7zn//07dv3vffe\n69Gjx+jRozt06HDixAl3d62nJmGEpF9zgiXp+VXp+Y/o24Ts8mE+Hc35zqRWYcAEwEZxcTEhxNPT\ns53fIxaLp059ek+u8kTtxsbGO3fueHt7K9zi6ufnRwgpLCykb7t16/bDDz80NDSUlJS4uLiIRKK2\nNQYjJP2iF+7kS5BMYgkDJgANQkJC1q9fT1/fu3fv9ddft7KymjVrFlPhthqNjY23b9+uqKhgs5eq\nqiqZTKY8e5uWVFZWyhfa2dn5+vq2ORoRjJAMAMmk9pAfMGEtIgCqpKSkoKBg48aNhJALFy5kZ2ev\nXbv2zJkzISEhTJ1evXrJZDKVH/f19V20aNGOHTta3VFjYyMhxMpKMVLQErpVhxCQDGFusKS48vH6\nk0+Ht+a5zF07Cbu4d13ySefwlepurZXeL7n7xfKHyduwFhHw3vnz5wkhQ4YMiYuLE4vFy5cvJ4Tk\n5+fL10lNTVX52aioqLi4OJYJHoFAQAiRSqUK5bSEbtUhBCQDQTJJJ+h1vM7hK9WtRUSv4z0k2zBg\nAh5LTU21sbGZOXPmhAkTxo8fr7JOaGioynIHBwd1m1RWJoTU1dUplNMSulWHkEMyECSTdMsuYGDX\nJZ/QW2uFzir+1qMDJtxaC7yUkZEREhISExMjEonGjh07ePDgkhK9/CcXi8UdOnQoKChoamqSL791\n6xYhxMPDQ7e7Q0AyHIWYRJNJRmwPD9ABk9uGA61MfIh+qXTd1Jo0tUu7ApiQ8vLygoKCOXPmjB8/\nftGiRZcvX3748OE777xDt9LZd4QQFxcXiSoFBQUSiWT16tUsd9evX78nT57k5PzpZJWRkUE36e6w\nCMElOwNDMkkfFGaKq8ww1edm1Odm0AwTFm8Fk3bu3DlCyJAhQ+hbS0tLS8unQ4uLFy/+9ttvdC54\nZGSkyqVKd+7cGRkZ+eKLL7Lc3eTJk9PT01evXk3vOiKElJSUfP755wKBYOLEie08FgUISIaGZJL+\n0AGTeHi4usVb6YCpJi0ZT60F03X8+HEvLy9XV1empK6uLjg4mBCSkpLy7rvv0kI6B0/ZgQMHYmNj\n2e9u4cKFO3fuTE9PDwoKmjFjRllZ2d69e+vr61evXq3zS3YISIZGL9zJPzNpw8nC4T6d8MwkXcGt\ntcBvV65cGTt2rHxJTExMamqqra1tt27d2nMbkErW1tanT59+8803Dx06RC/cOTo6btq0KSYmRrc7\nInj8hLHgmUkGw/JpFxgwganIy8vz8PCwtraWLywuLpZKpWzWM23ziZHeUSsSiTw8PCwsLNrwDa1C\nQDKaDScKmWQSwTOT9E/dgImBAROYAy6fGBGQjIbOaGCSSYSQ9aO914V6a/gItB8GTGDmuHxiREAy\npqLKx/LJJC8n2/gI/+E+HY3YJPOhbi0iBgZMwEtcPjEiIBkZkknGRa/jNeRexIAJzASXT4wISMaH\nZBIXsBkwYS0i4AEunxgRkIwPySTu0HBrLUWv4+HWWjBdXD4xIiBxApJJnEInPtTnZmgYMNkHDBIP\nD8eACUwOl0+MCEhcgWQSB2HABPzD5RMjAhKHIJnETZgpDnzC5RMjAhKHIJnEcSxvrcWACbiMyydG\nBCRuQTKJ+zBgApPG5RMjAhLnIJlkKrAWEZgiLp8YEZC4CMkkE4IBE5gWLp8YEZC4CMkkU4S1iMAk\ncPnEiIDEUUgmmSisRQQcx+UTIwISdyGZZNKwFhFwE5dPjAhInIZkkqljeWutfcAgDJjAMLh8YkRA\n4jQkk3ijIfdiddo+DJjA6Lh8YjSZgNTU1PTGG2/IZDKmxNLS8quvvtLwES73O3tIJvEJ1iICo+Py\nidFkAtKvv/46ceLEv/3tb7a2trTE0tIyJiZGw0e43O9aSc+vGvFlDvMWySRTx3KmuH3AIAyYQOe4\nfGI0mYB05MiRDz/8MDs728LCguVHuNzv2kIyiZewFhEYHpdPjJbGbgBbN2/e9PPzYx+NeGZOsGS4\nTyfmbUJ2+YYThRrqg0kQdnHvHL7S+8vMrm98Yh+gYtQrvV/yMHlb6bqpdz9f3pB70fAtBDAkkxkh\nzZs3z9HR0dLS8sqVK2KxeNiwYW+88YadnZ2Gj3D5D4E2KKp8PGJHTlHlY/oWyST+wVpEYABcPjGa\nTEAaOHBgXV1dREREz549f/755+Tk5D59+nz33XeWlmoHeVzu97ZBMskcYC0i0Csunxg5F5Bqa2tT\nUlKYt87OzqNGjWpubv7222/79+/fu3dvWn7w4MH33nvv008/HTNmjLqv4nK/txmSSeYDAybQBy6f\nGPUVkKqrq8PCwiIjI2fOnKm89dSpUykpKXl5eSKRyNfXNyoqytPTk24qLi7+61//ytTs379/XFyc\n8jc0Nzf369dv1qxZGibacbnf2wx3JpkbDJhAt7h8YrTS0/fu3r27rKysrq5OeVNsbGxiYqKNjY2/\nv39NTU1ycvKRI0e++OKLIUOGEEI8PT2vXLmi8JEHDx4UFBQEBQUJBAJaIhAIrKysVH4/v3k52cZH\n+MknkxIulw/r3gnJJL4SdnEXdnEXjwhXtxaR9H6JNL2kJj0Zt9aCqdPlLLvm5uaioqKzZ8+uWLFi\n165dKutkZmYmJiZ6eHj88MMPSUlJqampn376qVQqXbNmzePHj9V9c1FR0axZsy5efDbL6OrVq3V1\ndX5+5ni1ik5nYN4WVT6el3TdiO0BwxB2ce+65BPvLzM7h68UOquYBS69X3L3i+WFi196mLwNU/LA\nFOkyIBUUFISGhi5YsCA1NVVdnYMHDxJCVq1a5erqSktCQ0PHjBlz9+7dCxcuqPtUnz59fH1916xZ\n8+OPP5aXl587d27VqlWenp5hYWE6bL8JGe7Tcf3oZ5fp6HU8I7YHDObpTPEdme4bDoqHqxgM0Zni\nJeumFC5+qSZN7UpFABykyxzS77//zgSVa9eu7d69+6233lq4cKF8nWHDht2/fz8nJ4dZcIEQcuzY\nsbfeemv27Nlr1qxR9+V37tz54IMP0tPTW1pahELh0KFDN2zY4OzsrKE9PXv2lH/L2cumbYNkEhCs\nRQTsmMrJUJc5JJFIFBoa+vR7rVR8s1QqraiocHNzk49GhBAfHx9CSGlpqYYvd3Fx2blzZ0NDQ0VF\nhaurq8rvV8bZfm8/JJOA/DFgEg8PVzfxgQ6YHiZvsw8YJB4ejgyTeZI/EyoEJ04x6EoNNTU1Mpms\nQ4cOCuW0pLq6utVvsLOz8/T0ZBmNeE9lMomJT2A+6KwHtw0HNGSY6nMzmAyTtEL1cArAuAwakKRS\nKVE1eKIldCtoBckkkEcHTG4bDmAtIjBFBg1IdNK2cuChJcyUbtCKwjJ36fmPsMydmVMYMClXkN4v\nqUlPphMfNNx4C2BgBg1IdOm5hoYGhXJaonlhOlCHJpO8nJ6l5RIul6fnVxmxScARLBdvLVz8EgZM\nwAUGDUgikcjR0bGkpKS5uVm+vKioiBAikUgM2Rg+QTIJNMCACUyFoR8/4efn19jYmJubK1+Yk5ND\nCPH391fzIWgdkknQKgyYgOMMHZBGjx5NCNm6dStTUl5e/t133wkEgpCQEAM3hmeQTAI25AdM6m6t\nxYAJjMLQASkiIqJ79+5ZWVmTJ0+Oi4v76KOPpk6dWl9fHxkZ6eLiYuDG8AySSaAV+bWIMGACLjD0\nDT1CoXDPnj0ffPDBqVOn6IU7BweHlStXRkVFGbglvESTScwzk2gyKW1xf/koBSCPXscjhGDxVjA6\noz0PSSqVFhUVOTg4SCQSPT2YnMurrOuVwjOThvt0SovuZ8T2gAnBWkS8x+UTI+ce0KdDXO53fRvx\n5RUscwdtRh/CVJ+boTxgYmAtIhPF5RMjAhI/FVU+ll/mjl7KwzJ3oC0MmPiHyydGBCTeKqp87B37\nbKlNLydbJJOgbfDUWj7h8okRAYnPkEwC3aIDJg1zwTFg4j4unxgRkHgOySTQOQyYTBqXT4wISDyH\nZBLoD8sBk8r1isBYuHxiREDiPySTQK8wYDItXD4xIiCZBSSTwADU3VrLwICJC7h8YkRAMhdIJoFh\nYMDEcVw+MSIgmQskk8DA2AyYsBaR4XH5xIiAZEaQTALDY3lrrX3AIAyYDIPLJ0YEJPOSkF0u/5wk\nLyfbwjUqlnkG0LmG3IvVafswYDI6Lp8YEZDMzrykGwnZ5cxbJJPAkLAWkdFx+cSIgGR2FJJJhJC0\n6P5IJoEhsZz4YB8wCAMmnePyiREByRwhmQQcgbWIDI/LJ0YEJDOFZBJwB2aKGxKXT4wISOYLySTg\nGqxFZABcPjEiIJkvJJOAmzBg0isunxgRkMwakknAZRgw6QOXT4wISOYOySTgOAyYdIvLJ0YEJEAy\nCUwD1iLSCS6fGBGQQEUyKT7Cb26wxIhNAlAHaxG1E5dPjAhIQAiSSWCCsBZR23D5xIiABE8hmQSm\nCGsRaYvLJ0YEJHhGIZk0N1gSH+FnxPYAsNTqxAdCiH3AIPHwcAyYuHxiRECCZ5BMAlOHAVOruHxi\nRECCP0EyCXgAM8U14PKJEQEJFCGZBLyBW2uVcfnEiIAEKiCZBHyCAZM8Lp8YEZBABSSTgJcwYCLc\nPjEiIIFqSCYBX5n5gInLJ0YEJFALySTgNzZrEfFvwMTlEyMCEmiy4UTh+pOFzFskk4B/6HW8htyL\nZjJg4vKJEQEJNCmqfDwv6UZ6/iOmBMkk4CszWbyVyydGBCRoBZJJYFZ4f2stl0+MCEjQOiSTwNzQ\niQ/1uRkaBkwmuhYRl0+MCEjACpJJYJ74N2Di8okRAQlYQTIJzBnLmeL2AYO4P2DS04kxNTV1x44d\n169f79Kly+TJk5cuXWpjY6PtlyAgAVtIJgGwvLWWywMmfZwYd+zYER0dPX369OHDh1+/fn3Xrl2v\nvvrq0aNHLSwstPoeBCTQApJJAMTEb63V+Ymxrq7OxcVl6tSpcXFxtITGp+zs7AEDBmj1VZY6bBbw\n3txgyfrR3sxbeh3PiO0BMAphF3fxiHC3DQe8v8xUedus9H5JTXpyybophYtf0jCc4ofCwsKampoZ\nM2YwJYMHDyaE3L59W9uvQkAC7cwJlgz36cS8Tcgul1+GFcCsCLu4dw5f6f1lZtc3PrEPUHG1QHq/\n5GHytsLFL939fHlD7kXDt9AAvLy8srKyBg16dvjZ2dmEEF9fX22/CpfsQGtIJgGoZBJrEen7xPjL\nL7+MHDkyICDgzJkz2n4WIyTQmpeTrfycb7o0uBHbA8ARwi7uXZd8YrYDppaWlq+//nrQoEGenp7J\nyWqjsgYISNAWSCYBqCOfYRIPVzELXD7DVJPWlhO3nlRWVnp7e3/++efKmw4fPjx9+vTAwMBBgwYt\nXLhQOT9UXFz86quvLlmyZNmyZefPn3/uuefa0AAEJGgj5WTShhOFGuoDmBtmwNQ5fKXQWcUscOn9\nkrtfLKcTH7gwYNq6dWtRUVFtba1C+bJly8LCwo4ePSoWi6uqqr7++uu+ffueOHGCqZCbm/vSSy/V\n1tZevXr1ww8/tLa2blsDkEOCtlNOJsVH+A/36WjEJgFwVkPuxeq0fUZfvFXhxNjU1FRQUJCXl5eY\nmLhv3z5CyEcfffTuu+8yFdLS0kJCQnx8fE6dOuXl5UUIOXjw4PTp0yUSya1bt+zs7GQyWf/+/Tt3\n7pyamtqGm2HlWbXnw2DmaDKJuVhXVPl4XtJ13JkEoJJdwEC7gIGdw1eqW4uIDpgeJm8z5K21N2/e\nDAwM1FAhPj6eELJ582YajQghU6ZMmTZtWlJS0smTJydNmnT69OmrV69u3br17Nmz8h8MCgrq3Lmz\nVo1BQIJ2mRssKa58zCxzR5NJWOYOQB06U1w8PFzdrbV04kNNWrJh1iLy8PDYv38/fZ2dnb1lyxaF\nCmlpaVZWVuPGjZMvnDRpUlJSUlpa2qRJkzIzMwkhMTExCh88fvx4aGioVo1BQIL2mhMsSc+vYpa5\nS8gu9+pkuy7UW/OnAMyZsIs7nfugbi0i6f0SafrTTSoHTMXFxYQQT0/PdrZELBZPnTr1aauEQoWt\njY2Nd+7c8fb2trOzky/38/MjhBQWFhJC1q5du3bt2nY2g8KkBmgvhVnghJCEy+Xp+VXGag+ACWF5\na23puqkKM8VDQkLWr19PX9+7d+/111+3srKaNWsWU+G2Go2Njbdv366oqGDTvKqqKplM5uTkpFBO\nSyorK7U83FZghAQ6gGQSQHuwHzDRW2vrB4YXFBRs3LiREHLhwoXs7Oy1a9eeOXMmJCSE+UivXr1k\nMpnK3fn6+i5atGjHjh2tNqyxsZEQYmWlGCloCd2qQwhIoBtIJgG0H8sM0w9f7yCEBHWxi4uLE4vF\ny5cvJ4Tk5+fL10xNTVW5i6ioqLi4OHd3VjMmBAIBIUQqlSo2QypltuoQAhLoDJJJADrR6oDpzK07\n1pYWMyOmj+rRdfa8SJVfom5CgYODA/u5Bg4ODoSQuro6hXJaQrfqEHJIoDNIJgHolroM05WHjS93\nsY7yFdk+rp7ydmyQRJyzMUoft9aKxeIOHToUFBQ0NTXJl9+6dYsQ4uHhodvdISCBLikvczcv6boR\n2wPAAwprEd1/LCupaw7ztBshsYnoZn8w5LlHdQ3rd+2laxHlJn1BP+Xi4iJRpaCgQCKRrF69muXe\n+/Xr9+TJk5ycP61XmZGRQTfp9khxyQ50DMkkAD2haxGdsfElx94YPHAgKfkfIcTS4tnAIvtm/p0r\na+3PJtoFDJw1YbTVc67KX7Jz587IyMgXX3yR5U4nT56cnp6+evXqtLQ0WlJSUvL5558LBIKJEye2\n/6DkISCB7iGZBKA/P2Zke3l5Bf8zlVmLqKG5JdBJSAhJu/tkYU8RnZIXRYjQUsVaRAcOHIiNjWW/\nu4ULF+7cuTM9PT0oKGjGjBllZWV79+6tr69fvXo1LtmBCUAyCUB/rly5MnbsWEKIXcBAunjr0ukT\nLjyySCqod3cQ2FtZMDXlF2+VViiuVMSStbX16dOnp06devXq1ZiYmE8++aS+vn7Tpk1aRTWWsLgq\n6EtCdrn8Mym8nGxxZxJA++Xl5Xl4eCisqJ2fc6n2RlbXRwUaFm+1DxgkHh4evOj9tp0Y6R21IpHI\nw8PDwsKi9Q9oDwEJ9GjDiUImmUQImRssQTIJQK/oTHGVi7dSpx8JF58uNnCrWMIlO9AjPDMJwMDo\nTHG3DQfUrUV0+pHignXcgYAEeoRkEoBRyM8U7xy+8lm5s/vPdTpeXkGHEJBAv3BnEoARKdxaKx+c\nOAgBCfRubrBk/ehnc77pnUlGbA+AuWEGTPp+ulI7ISCBISCZBACtQkACQ6AX7rycbJkSJJMAQAEC\nEhiIl5NtfIQ/8xbJJABQgIAEhjPcpyOSSQCgDgISGBSSSQCgDgISGBSSSQCgDgISGBqSSQCgEgIS\nGAGSSQCgDAEJjAPJJABQgIAExoFkEgAoQEACo1GZTCqqfGzEJgGAESEggTEhmQQADAQkMDKFZFJ6\n/iMkkwDME0cD0rVr1wYOHFhfXy9fmJ2dvXDhwlGjRoWHh2/fvv3JkyfGah7oEJJJAEBxMSDdu3dv\n48aNlZWV8o9X/+9//ztr1qz6+vq5c+cGBwfv3r07Ojqax89fNytIJgEAIcTK2A34k/379ycmJubl\n5clkMoVNW7Zs8ff337Nnj6WlJSGkV69eq1atysjIGDx4sDFaCjpGk0nrTz69WEeTSWnR/YzbKgAw\nJG6NkAIDAxcsWLBly5apU6fKlz948ODWrVsTJkyg0YgQMmbMGGtr6/PnzxujmaAXSCYBmDluBaRe\nvXpNmDBhwoQJL7zwgnx5SUkJIcTLy4spEQqFrq6uZWVlBm4h6A+SSQBmjlsBSR06u8HBwUG+0N7e\nvq6uzkgtAr1AMgnAnBknh1RbW5uSksK8dXZ2HjVqlIb6dPKCwhSGlpYWTGrgHySTAMxWuwJSdXV1\nWFhYZGTkzJkzlbeeOnUqJSUlLy9PJBL5+vpGRUV5enrSTZWVlVu2bGFq9u/fX3NAsrOzI4Q0NDTI\nFzY0NEgkkva0H7hpXah3en5Vev4j+pYmk9aFemv+FACYunYFpN27d5eVlam8bhYbG5uYmGhjY+Pv\n719TU5OcnHzkyJEvvvhiyJAhhBBPT88rV66w3xENPEVFRUxJc3NzaWkpptjxVXyE34gdOczFuoTL\n5cO6dxru09G4rQIAvdI6h9Tc3FxUVHT27NkVK1bs2rVLZZ3MzMzExEQPD48ffvghKSkpNTX1008/\nlUqla9asefy4LfkAFxcXiUTy448/MiXp6elSqXTAgAFt+DbgPi8n27TF/Zm3SCYBmAOtA1JBQUFo\naOiCBQtSU1PV1Tl48CAhZNWqVa6urrQkNDR0zJgxd+/evXDhQtsaOnv27MuXL8fGxt66dSslJWXj\nxo3u7u4jRoxo27cB93k52WKZOwCzovUlO4lEsn37dvr62rVru3fvVq6TmZkpEAiGDRsmXzhy5Mhj\nx45dunRp5MiRbWjonDlzqqur4+PjExMTCSGBgYGbNm2ytbXV/KmePXsyr3/99dc27BeMCMkkAJ2Q\nPxNymdYBSSQShYaGPv2wlYqPS6XSiooKNzc3hWjh4+NDCCktLWWzl2nTpk2bNk2+RCAQrFixYunS\npWVlZZ06dXJ0dGTzPQhCpg7JJID2kz8Tcjk46f4+pJqaGplM1qFDB4VyWlJdXd2eLxcIBB4eHiyj\nEfAAkkkA5kP3AUkqlRJVgydaQrcCsIdkEoCZ0H1AEggERFXgoSV0K4BW1oV6Y5k7AN7TfUBSeRMr\nU0K3AmgLy9wB8J7uA5JIJHJ0dCwpKWlubpYvp7e1Ym0FaBskkwB4Ty+Lq/r5+TU2Nubm5soX5uTk\nEEL8/f3VfAigFXQ5cOZtUeXjETtyjNgeANAtvQSk0aNHE0K2bt3KlJSXl3/33XcCgSAkJEQfewQz\nMTdYMjf42SC7qPIxkkkAvKGXgBQREdG9e/esrKzJkyfHxcV99NFHU6dOra+vj4yMdHFx0ccewXys\nG+0tn0xaf7IQySQAftBLQBIKhXv27AkNDb158+aWLVv27NnT0NCwcuXKFStW6GN3YFYUkkmEECST\nAPjBQq+PFJJKpUVFRQ4ODhKJxMLCQn87Uqlnz55YqYGvErLL5e9G8nKyLVwzyIjtATAVXD4x6veJ\nsUKhsEePHi4uLoaPRsBvSCYB8I9pPMIcQBmSSQA8g4AEpgrJJACeQUACE4Y7kwD4BAEJTNvcYEnh\nmkHMtbuiyscWK88kZJcbt1UA0AYISGDylK/dbThZiGt3ACYHAQn4ANfuAHgAAQl4QnkiOB6bBGBa\nEJCAPxQmgidklyOZBGBCEJCAP5BMAjBpCEjAK0gmAZguBCTgGySTAEwUAhLwEJJJAKYIAQl4CMkk\nAFOEgAT8hGQSgMlBQALemhssWT/am3mLZBIAxyEgAZ/NCZYM9+nEvEUyCYDLEJCAzxQu3BEkkwA4\nDAEJeA7JJABTgYAE/IdkEoBJQEACs4BkEgD3ISCBWUAyCYD7EJDAXCCZBMBxCEhgRpBMAuAyBCQw\nL0gmAXAWAhKYF5XJpPT8KmO1BwAYCEhgdpSTSfOSrhuxPQBAISCBOUIyCYCDEJDATCknkzacKDRi\newAAAQnMlHIyKeFyOZJJAEaEgATmC8kkAE5BQAKzhmQSAHcgIIG5QzIJgCMQkMDcIZkEwBEISABI\nJgFwAgISACFIJgFwAAISwFNIJgEYFwISwFNIJgEYFwISwDNIJgEYEQISwJ8gmQRgLAhIAIqQTAIw\nCgQkAEVIJgEYBQISgApIJgEYHgISgGpIJgEYGAISgFpIJgEYEgISgFpIJgEYEgISgCZIJgEYDAIS\nQCuQTAIwDAQkgNYhmQRgAAhIAK2jF+68nGyZEiSTAHQOAQmAFS8n2/gIf+YtkkkAOoeABMDWcJ+O\nSCYB6A8CEoAWkEwC0B8EJAAtIJkEoD8ISADaQTIJQE8QkAC0hmQSgD4gIAG0BZJJADqHgATQFkgm\nAegcAhJAG6lMJhVVPjZikwBMGgISQNshmQSgQwhIAO2ikExKz3+EZBJA2yAgAbQLkkkAuoKABNBe\nSCYB6AQCEoAOIJkE0H4ISAC6gWQSQDshIAHoBpJJAO2EgASgM0gmAbQHAhKALiGZBNBmCEgAOrYu\n1BvJJIA2QEAC0D0kkwDaAAEJQPe8nGzTFvdn3iKZBMAGAhKAXng52SKZBKAVBCQAfUEyCUArCEgA\neoRkEgB7CEgAeoRkEgB7CEgA+oVkEgBLCEgAeodkEgAbCEgAhoBkEkB7IEo8AAATjUlEQVSrrIzd\nANWuXbu2cOHC06dP29vb05KmpqY33nhDJpMxdSwtLb/66isjNRBAOzSZ5B2bQd/SZFLa4v7yUQrA\nzHFxhHTv3r2NGzdWVla2tLQwhfn5+enp6W5ubt3lGLGRANqiy4Ezb4sqH4/YkWPE9gBwDbdGSPv3\n709MTMzLy5MfCVE3b950dHT8v//7PwsLC6O0DaD95gZLzuZXJWSX07dFlY83nChcF+qt+VMAZoJb\nI6TAwMAFCxZs2bJl6tSpCptu3rzp5+eHaASmbt1ob/nLdOtPFiKZBEBxKyD16tVrwoQJEyZMeOGF\nFxQ23bx5s1OnTsuXLx82bNiECRP+8Y9/NDQ0GKWRAO2hcGcSIYQHdyZVVlZ6e3t//vnnypsOHz48\nffr0wMDAQYMGLVy48Pbt24ZvHpgKbgUkDW7evJment6lS5elS5cGBQX961//mjdvnvKVPQDu418y\naevWrUVFRbW1tQrly5YtCwsLO3r0qFgsrqqq+vrrr/v27XvixAmjNBK4zzg5pNra2pSUFOats7Pz\nqFGjNNRvbm5evHhx//79e/fuTQiZMmVKYGDge++9d/LkyTFjxui9uQC6xoNkUlNTU0FBQV5eXmJi\n4r59+5QrpKWlbd++3cfH59SpU15eXoSQgwcPTp8+ff78+bdu3bKzszN0i4Hz2hWQqqurw8LCIiMj\nZ86cqbz11KlTKSkpeXl5IpHI19c3KirK09OTbqqsrNyyZQtTs3///poDkkAgmD17tnzJX//61w0b\nNvz8888ISGCi1o32Ts9/xFysW3+ycFj3TsN9Ohq3VezdvHkzMDBQQ4X4+HhCyObNm2k0IoRMmTJl\n2rRpSUlJJ0+enDRpkgEaCaalXZfsdu/eXVZWVldXp7wpNjb2jTfeOHPmjEgkqqmpSU5Onjhx4vnz\n5+lWT0/PK3Li4uI07+jBgwdZWVnNzc1MiUAgsLKyUrlrAJNg6skkDw+P/X94++23lSukpaVZWVmN\nGzdOvpDGobS0NAO1EkyK1iOk5ubmkpKS4uLiw4cPp6amqqyTmZmZmJjo4eGRkJDg6upKCDlx4sSK\nFSvWrFlz4sQJW1ut7wQsKiqaNWtWXFzckCFDaMnVq1fr6ur8/Pw0fxCAy2gyiVnajiaTCtcM0utO\ni4uLCSHM5Yo2E4vFzGxYoVCosLWxsfHOnTve3t4Kl+bo72xhIVZOAhW0HiEVFBSEhoYuWLBAXTQi\nhBw8eJAQsmrVKhqNCCGhoaFjxoy5e/fuhQsX2tDKPn36+Pr6rlmz5scffywvLz937tyqVas8PT3D\nwsLa8G0A3DE3WDI3WMK8pckkve4xJCRk/fr19PW9e/def/11KyurWbNmMRVua1RRUcFmL1VVVTKZ\nzMnJSaGcllRWVurmYIBftB4hSSSS7du309fXrl3bvXu3cp3MzEyBQDBs2DD5wpEjRx47duzSpUsj\nR47UdqfW1tZfffXVBx988Oabb7a0tAiFwqFDh27YsMHa2lrzB3v27Mm8/vXXX7XdL4ABKCeTPJ1s\n5aOUDpWUlBQUFGzcuJEQcuHChezs7LVr1545cyYkJISp06tXLw3zVxctWrRjx45Wd9TY2EgIsbJS\nPMPQEroVDEb+TMhlWgckkUgUGhr69MNK/9sIIVKptKKiws3NTeHSnI+PDyGktLSUzV6mTZs2bdo0\n+RIXF5edO3c2NDRUVFS4urqq3LUyBCHgPoVl7gghG04WDvfppI9l7mged8iQIXFxcWKxePny5YSQ\n/Px8+ToaLn4QQtzd3dnsSCAQEEKkUqlCOS2hW8Fg5M+EXA5Oup/2XVNTI5PJOnTooFBOS6qrq9vz\n5XZ2du2/9g3ANQZLJqWmptrY2MycOXPChAnjx49XWYf5i7M9HBwcCCHK045oCd0KoED3N8bSv4DU\nDdWV/2ICAKIqmaSP5/hlZGSEhITExMSIRKKxY8cOHjy4pKRE53shhIjF4g4dOhQUFDQ1NcmX37p1\nixDi4eGhj52CqdN9QMJQHaBtFJa5S8guZ+6c1Yny8vKCgoI5c+aMHz9+0aJFly9ffvjw4TvvvEO3\n0tl3hBAXFxeJeqtXr2a5u379+j158iQn50+LUGRkZNBNujss4A/dX7KjszyVF5qjJbg9G0AdfSeT\nzp07Rwhh7p2wtLS0tHz6J+nFixd/++03ej08MjJS/skvCl588UWWu5s8eXJ6evrq1auZu45KSko+\n//xzgUAwceLENh8F8JjuA5JIJHJ0dCwpKWlubpYfDxUVFRFCJBK9zB0C4Ae9JpOOHz/u5eXF3IxB\nCKmrqwsODiaEpKSkvPvuu7SQzsFrv4ULF+7cuTM9PT0oKGjGjBllZWV79+6tr69fvXo1LtmBSnpZ\ny87Pzy8rKys3N7dPnz5MIR25+/v762OPALyhvMzdvKQb8ouxttmVK1fGjh0rXxITE5Oammpra9ut\nWzeRSNT+XciztrY+ffr0m2++eejQIfrr7+jouGnTppiYGN3uCHhDLwFp9OjRWVlZW7du/eabb2hJ\neXn5d999JxAI5G93AACVFO5MSsguH+bTsf13Jh04cEBhaLJkyZIJEyZIpdJ2Pn950qRJKq/yde3a\ndf/+/Y2Njbdv3xaJRB4eHnikGWigl4AUERGRlJSUlZU1efLkcePG3bt379ixY/X19a+//rqLi4s+\n9gjAJ3pKJqmMOga4j8La2jogIEDfewEe0MvzkIRC4Z49e0JDQ2/evLlly5Y9e/Y0NDSsXLlyxYoV\n+tgdAP/w75lJAK2y0DCdpv2kUmlRUZGDg4NEIjH8UL1nz55YqQFM2oYThetPPlvabm6wRCfJJDBn\nXD4x6veJsUKhsEePHi4uLrhwDNAGc4Ilw306MW91fmcSAKeYzCPMAcyQwoU7QsiGk4Um9MwkAK0g\nIAFwGpJJYD4QkAC4bm6wZP1ob+atnpa5AzA6BCSe4/JS8wbDg05ofzKJB53QfugEjkNAAjABSCaB\nOUBAAjANSCYB7yEgAZgM+WSSl5NtfARWhgRe0e+NscaF68XAP1L7zvf6ziaEuF38p7HbAqaKszfG\n8jkgAfBSUeVjXT0hCYBTEJAAAIATkEMCAABOQEACAABOQEACAABOQEACAABOQEACAABOQEACAABO\nQEACAABOQEACAABOQEACAABOQEACAABOsDJ2A3Tv1KlTKSkpeXl5IpHI19c3KirK09PT2I3Sr8OH\nD589e/bWrVtisbh3797z589//vnnFeqYT7c0NDTMmjXLyclp165dCpt43wl37tw5derUpUuXCgsL\n3d3dp06dOnr0aIU6/O6ElpaWkydPHjp0qLi4uGPHjt27d587d66Pj49CNV52QnV1dVhYWGRk5MyZ\nM5W3tnrIXOgTvq1lFxsbm5iYaGNj4+/vX1NTU1BQYGNj88UXXwwZMsTYTdOL5ubmN9988/Tp0yKR\nKCAg4O7du8XFxba2tnFxcQMGDGCqmVW3vP/++8nJyQEBAd9//718Oe874bfffpsxY8aDBw+ef/75\nzp07X79+nRASHR29bNkypg7vO+Htt98+cuSIg4ND7969i4uL7969a2Vl9cknn7z66qtMHb52wrZt\n23bt2vXWW28tXLhQYVOrh8yVPmnhkUuXLvn6+o4aNaq0tJSWHD9+3M/Pb+jQoQ0NDcZtm558++23\n9G8Z5gCTk5N9fX1feeWVxsZGWmJW3XL69Gl/f/8+ffqEhYXJl/O+E37//fcRI0b07t37/PnzMpms\npaWloKBgwIABvXr1Yg6Z951w+vRpX1/fSZMmPXr0qKWlRSaTpaam+vr6Dhw4sKmpidbhWSc0NTUV\nFhamp6cvX77c19fX19d3586dCnVaPWTu9AmvckgHDx4khKxatcrV1ZWWhIaGjhkz5u7duxcuXDBq\n0/QlISHBysrqo48+srV9+jyCadOmvfLKK/fu3bt16xYtMZ9uefDgwZo1axYvXuzk5KSwifedcPLk\nybKyspiYmMGDB1tYWBBCvL29o6KiOnXq9PPPP9M6vO+En376iRASFRXVsWNHQoiFhcXYsWO7d+/+\n8OHDgoICWodnnVBQUBAaGrpgwYLU1FR1dVo9ZO70Ca8CUmZmpkAgGDZsmHzhyJEjCSGXLl0yUqP0\nqKWl5c6dOz179uzSpYt8ube3NyGktLSUvjWfblm7dq1EIlm8eLHyJt53QkpKikAgmDx5snzhokWL\nMjIyxowZQ9/yvhOU/xAhhDQ3N1taWj733HP0Lc86QSKRbP/D/PnzVdZp9ZC50yf8mdQglUorKirc\n3NyYsQJF85nM2ZlPmpubN2zYoDx/IS8vjxDi5eVFzKlbkpKSLly48P333wsEAoVN5tAJWVlZffr0\nEYlEt2/fvnLlSllZWffu3QcOHMiciM2hE0aNGrVt27aEhISQkBAHBwdCyMmTJwsLC19++eVOnToR\nPnaCSCQKDQ2lr62sVJzPWz1kTvUJfwJSTU2NTCbr0KGDQjktqa6uNkaj9MvKymrKlCkKhRkZGRcv\nXuzevXv37t2J2XRLUVHR3//+92XLlvXo0UN5K+87oba2trGxUSKR7N69e+vWrUy5WCzetGnTqFGj\niBl0AiHE09MzISFh6dKlISEhQUFBBQUFhYWFgwcP/uSTT2gFc+gEBa0eMqf6hD+X7KRSKVH1NwIt\noVt578cff4yOjraxsfnwww/pQMEcuqW5uTkmJqZXr16RkZEqK/C+Ex4+fEgIuXTp0mefffb++++f\nP3/+/Pnza9asefz48VtvvVVWVkbMoBOo2tpamUxWVVWVnp5eWFhICKmvr3/06BHdaiadIK/VQ+ZU\nn/AnIMmff+XREuXLODxTU1Pz7rvvLlmyxMHBIS4urn///rTcHLrliy++yMvL27x5s6Wl6v/PvO8E\neiCVlZUxMTEzZ850dnZ2dnaePXv2smXLnjx5snv3bmIGnUAI+fe//x0dHe3s7JyYmHjt2rXs7Ow1\na9Zcv3590qRJv/76KzGPTlDQ6iFzqk/4E5Ds7OwIIQ0NDQrltIRu5aszZ8785S9/+f7778eNG5eS\nkiJ/BxLvu+Xnn3/euXPnokWLnJycav8gk8lkMlltbW19fT0xg06gk8oIIWFhYfLl48ePJ4TcuHGD\nmEEnEEIOHjwoFAp37dr10ksvWVlZicXi2bNnv/XWWw0NDT/++CMxj05Q0Oohc6pP+BOQRCKRo6Nj\nSUlJc3OzfHlRUREhRCKRGKdZ+hcfH7948WJra+s9e/Z8/PHHNHnL4H23/PLLL83NzR9//PEAOXfv\n3r1x48aAAQNmzJhBzKATOnbsKBQK7e3taSafQf8zVFZWEjPohEePHt24ccPX15eZu0zRtSrobDHe\nd4KyVg+ZU33Cn4BECPHz82tsbMzNzZUvzMnJIYT4+/sbqVH6debMmc2bNw8YMODo0aMvv/yyyjr8\n7paAgIBoJSKRyNnZOTo6Ojw8nFbjdycIhcKXX365vr6eposYdL4lswAMvzvBwcFBIBDU1tYqlFdV\nVZE/YjPheyeo1Oohc6hPDHkXrr4lJib6+vrOnDmTKblz584LL7zg5+dXVlZmxIbpSVNT06uvvhoU\nFFRbW6uhmrl1S0tLy/DhwxVWauB9J+zbt8/X1zcmJoYu09DS0tLc3Lx48WJfX9/Dhw/TEt53wtSp\nU319fY8fP86UyGSyN99809fXNz4+npbwuBNOnTqlcqWGVg+ZO33Cn2nfhJCIiIikpKSsrKzJkyeP\nGzfu3r17x44dq6+vf/31111cXIzdOt3Ly8srLi728PD45z//qbx13rx5bm5uxPy6RSXed8KUKVPS\n09OPHDly9+5depHq+PHj2dnZQUFBEyZMoHV43wnr16+PiIhYunRp7969x44da21tnZqaeuXKlcDA\nwNdee43W4X0nKGv1kLnTJ3xbXPXBgwcffPDBqVOn6PVQBweHRYsWRUVF8XL+zP79+9euXatua1JS\nUr9+/ehrs+oWQsiIESM6deqksLgq7zvhyZMnmzZtOnfuHL1w5+TkNH78+LffflsoFDJ1eN8JV69e\n3bZtW2ZmJn1rbW0dERGxZMkS+fts+NoJp0+fjo6OVrm4aquHzJE+4VtAoqRSaVFRkYODg0Qioet6\nAUG3EELMoxPu3LljYWGhIR3N+05oaGgoLS21sbFxdXVVd0rlfScoa/WQjd4n/AxIAABgcng1yw4A\nAEwXAhIAAHACAhIAAHACAhIAAHACAhIAAHACAhIAAHACAhIAAHACAhIAAHACAhIAAHACAhIAAHAC\nAhIAAHACAhIAAHACAhIAAHACAhIAAHACAhIAAHACAhKAITQ2Nt69exePHwPQAAEJQC927do1dOjQ\n8+fP5+TkvPbaawMGDBg2bFj//v3Xrl1bX19v7NYBcBECEoBe1NbW3rt3LzMzc/78+Q8ePAgNDR0w\nYEBDQ8P+/fvfe+89Y7cOgIusjN0AAD7btWvX4sWLly5damlpSQj53//+N3369BMnTlRUVHTp0sXY\nrQPgFoyQAPSoV69ey5Yto9GIEPLCCy/4+/vLZLLffvvNuA0D4CAEJAA9GjJkiIWFhXxJ9+7dCSFV\nVVVGahEAdyEgAejR888/b+wmAJgMBCQAPVIYHgGABghIAADACQhIAADACQhIAADACQhIAADACQhI\nAADACRZY7REAALgAIyQAAOAEBCQAAOAEBCQAAOAEBCQAAOCE/weWH+dtSFQ4LgAAAABJRU5ErkJg\ngg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "for kappa = [10 100 1000 1e4]\n", " bound = 2*( (sqrt(kappa)-1)/(sqrt(kappa)+1) ).^(0:100);\n", " semilogy(0:100,bound*Aerr(1))\n", " hold on\n", "end\n", "ylim([1e-16 1])\n", "text(60,1e-15,'\\kappa=10')\n", "text(95,1e-8,'\\kappa=10^2')\n", "text(95,1e-3,'\\kappa=10^3')\n", "text(95,1e-1,'\\kappa=10^4')\n", "title('Effect of conditioning on convergence')\n", "xlabel('n')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The effect is the same for MINRES and CG, but measured in different terms. In particular, CG does *not* guarantee a nonincreasing residual norm." ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAABcSAAAXEgFnn9JSAAAA\nB3RJTUUH4AsHFR0WIcIAXAAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJ\nbmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAwNy1Ob3YtMjAxNiAxNjoyOToyMpEv/jUAACAA\nSURBVHic7d15XFTV/z/wwzI6LCKgoCMgIAgq2idRMjGTxcTS8oca0scNRS2xXEL9+sn6oBZWlqaW\nZaYJuES4lH1A01yg3FDDstBUllFEEBQEFNAR5vfH0ev1zsIAs9x75/V89MfMmTsz54DdN+ec9znH\nQqlUEgAAAFOzNHUFAAAACEFAAgAAnkBAAgAAXkBAAgAAXkBAAgAAXkBAAgAAXkBAAgAxaGhomDRp\n0uDBgz/99FNT1wVayNrUFQAA0IOEhIQtW7YQQnr16mXqukALoYcEAIL3888/L1++3NS1gNZCQAIA\nYSsqKpowYUKnTp0GDBhg6rpAqyAgAYCAKRSKcePGVVZWbt++3dXV1dTVgVZBQAIAAVu0aNGJEyf+\n+9//hoaGmrou0FpIajAvJSUlR48ezcnJycnJUSqV/v7+gYGB48ePb9OmjamrBs1w5syZixcv3r17\nNyAgYNCgQaauTqtUVVXV1tYSQmQymdoLbt68+eeff167dq1Hjx69e/e2s7NjXvrxxx9XrVoVFhb2\n3nvvGam6YFBKMBubNm1i/8/M8PT0/O6770xdO9BJeXl5SEgI87sbN26cqWvUWn379qVtUSgUnJfu\n3bs3Z84c9r9VW1vbb775hrnghRdeIIT07dt36NChQ4cO7dixIyHEzc1t6NChBQUFxm0H6AF6SGbh\n/v37kyZN+v7779W+euXKlX//+9+EkOjoaOPWC5pt/vz5mZmZpq6F3uzdu/fs2bNqX1IqlWFhYceO\nHWMX1tbWTp8+/Z9//mEvNuJ8QnFxcXFxcU1NjSEqDAZlocR5SGZg5cqV8+fPp487deo0fvz44OBg\nS0vLvXv3bt68uaGhgRAilUpPnz7du3dvk9YUmtChQ4eKigpCSFxc3JIlS9q0adO+fXtTV6olrl+/\nnpKSsmTJknv37tEShUJhbf34T+SdO3e++uqrhBArK6sFCxY8/fTTu3fvTktLI4RYW1tfunTJ29v7\n3Llz9KdBLV68+Pjx4yNHjoyPjw8KClI7HgC8ZuouGhjcjRs3HBwc6K+7W7du+fn57Fd/+uknKysr\n+urChQvVfoJcLi8sLGxsbNTx6/Ly8hoaGlpZ7QcPHly8eLG6ulqXi69du3bx4kXVwuLiYr1/l44N\nLCgoyMvLUx2GavH3KpXKu3fvMv/nZmVlqb2mWb8stT83Qztw4ICNjY3qvYj9s2psbOzTpw8tX7Zs\nGVP+4osv0sKpU6eqfvLLL79MCJkxY4YxmgEGgIAkfkzfiBDyyy+/qF4QGRlJX/X392eX//PPPy++\n+CLzB3i7du1eeOGFv/76i33NmDFj3Nzc3Nzcfvzxx/3793fv3p1ebGNjM2XKlDt37tDLoqOj6WXT\npk1jvz0xMZGWT5kyhSk8efLk4MGDbW1tCSGWlpZ9+vT54IMP2HdY9pfm5OQEBQURQmQyGX21sbHx\n/fffd3d3pzVxd3dPS0vbtGkTfQvnbtWs79LSQMaFCxciIiKYH5qVldWoUaNOnz7NuazJ71Xl5ubW\npUsX5lfp4uLi5ub2+eeft+yXpfpzM5r//e9/av84Zgekq1evMuVlZWVM+Z49e2ih2mojIAkdApL4\nMemwnHjDqKmpuf4IU5icnCyVSlXvGhKJ5KuvvmIuGzx4MC2PiopSTdUbPXo0veyrr75ibpT37t1j\n3h4QEEDLk5KSaMlHH33E9NjYxo4dW1tby/nSZcuWtWvXjnOHioqK4rzX2to6LCyMPh4zZgzz7c36\nLu0NpLZt20ZjjGoF9u3b16zvVaV6PSFk+fLlLfhlqf25aVFXV9f6Li8jNzd33iNjx45lassOSEeP\nHqWFXbt2Zb+3uLiYlltYWLD/IYE4ICCJX6dOnej/w+PHj9fxLVeuXGHG38PCwpKSklJSUiIiImhJ\n27Zt//nnH3olc48jhLRr1y4uLm727Nn29va0xNLS8urVq0qlsqysjJke+Pnnn+l78/LyaEmbNm1u\n376tVCqPHTtmYWFBv2Lx4sU//PDDvHnzXFxcmNso50uZSnp6egYFBSmf/Ov71Vdf3blz51dffdWz\nZ0+mkAlIzf0u7Q1UKpXXr19nyocMGZKSkrJ27dp+/frREnt7+6qqKt2/V1V6evrOnTuZynz66afp\n6el5eXkt+GWp/ty0OH78OH0Lu2TEiBGEEC8vL6Zw5cqVS7TaunWr6ocfOXJEbUDatm0bLXzmmWfY\n1zc0NFhaPlw9mZeXp73mIDgISCJXWVnJ/A+/aNEiHd81YcIE+pbg4OD79+/TwgcPHtAsW0LIK6+8\nQguZe1yHDh1oUFE+eZdhegbDhg2jJTNnzqQlK1eupCUvv/wyLWEygLdt28ZUhhmlcXFxoR0IdpB4\n+eWXy8vLmYuZl5gaKpXKmzdvOjo6cgJSc7+ryQZOmzaNlvTt27e+vp4WlpSUMPMlP/zwg+7fq9ad\nO3eY7z1x4kSLf1mqPzctFixYQAiJjo6mTwkhX3/9NW0sO0qxhxPVioiIUP1wTQHpww8/pIXPP/88\n5y1MX/Dw4cO61B8EBDs1iJxUKmVGh+rr63V814kTJ+iDuLg4iURCH1tZWb355pucCxgvvfQSM4ER\nFBTEfOmNGzfoAyan/KefflIqlYQQ5i48btw4QkhVVRXN37W0tPTw8Dj5SIcOHZycnAgh5eXlhw4d\nYn9px44dv//+e7r6hPrjjz/og9dff50p7NChw5gxY9hvbMF3NdnArKws+mDWrFlt27aljzt37nzk\nyJH09PT09PSAgIAWfG+TWvDLUv25aXH48GFCSGJiYnFxsYWFxddffz1jxox169adOXOGHSBtmsL8\nTHTBXExTQNmYEqzmFh+sQxI5qVTq7e1NB8fkcrnaawoLC3NycuhjmsVUWFhIn7IHuwhrY//y8vLy\n8nJmoIk8uczezs7O3d39ypUrhJD79+/TwsjIyDfeeOP+/fvFxcW///67t7c3XWIilUpfeeUVQsjf\nf/9Nr2xsbHz++efVVrWkpIT99MUXX2Tna5WXlzOrTzw9PdlX+vj4sJ+24Lu0N/D+/fsFBQX01aef\nfpr9RvaOn8yqGt2/V7va2toW/LI4Pzctqqurf//9d0JIenr6nDlzSktL6QhwmzZtmNFIihmA1QsP\nDw/6gG7iwFAoFAqFgj7u2rWrHr8R+AABSfwCAgLozeLQoUPV1dVMCjhj/vz5u3fvJoQ4ODhUVFTc\nvXu3sbGRvsSZe2c/Ze4LFOfvX2agn+Ho6BgREUHnePbs2ePr60v/1H3xxRfpBDuzHsXS0tLX11dt\nW5iKUW5ubuynDg4OFhYPl9ZVVVWxX2InTLfsu7Q38O7du8xf7kw3RVULvle7Bw8etOCXxfm5abF/\n/3764MCBA6NGjVq4cGFYWNjEiRNVf7/6xQQbJsxTTPS1srJqcpAQBAcBSfxGjhxJB8dqamrWrl37\n7rvvsl+9evUqM3T2/PPPW1lZOTg4eHh4FBUVEULy8vL+9a9/MRczfwU7Ojq24HYQHR3NDki0kI7X\nEUKYjDulUvn7778zCQJacG79bdu27dy5M+1hHDx4MDg4mHnpl19+YV/Zgu/SzsnJSSaT0a++dOnS\nU089xbyUnJx87tw5QsjQoUMDAwP1+70t+2VpCZkcP//8MyFk+fLl//nPfxobG3/88ccxY8bExMTc\nu3ePM2L27rvvcv4I4OjZs2dcXJyO38t0cKuqqi5fvsxk2586dYo+cHNzU5upCMJm0hksMIaGhgbm\nPmhhYfHxxx8zs8eVlZV0uIzavHkzLf9//+//0ZIXX3yR/VFMRnV4eDgtYebJ33vvPfaV3t7etHzD\nhg1MYU1NDTNSRG9ntra27KU8zOIh9oR/Q0PD7NmzIyMjIyMjS0pKtHypUqlkpo4cHBx+++03pVL5\n4MGD//73v0wbmaSGVn6XagPpIhjOD62iooJJacvIyND9e9VSm9TQ+l+WFt26dSOE/Pnnn0zJO++8\nQx5lKNy8eZO5h+g3qYFd29jYWFrS2NjI7CQyf/58HZsAAiKkgHTq1KkZM2aEh4e/+uqra9asYRKZ\noElHjx6lqcaUra1tSEhIWFgYk3tGCBkxYgRz/d9//838ET1hwoRDhw4dOXIkNjaWllhaWp46dYpe\n2ayApFQq2etOCCGvvvoq+9WtW7fScgcHh9dff/3gwYM7d+5kQmZoaKj2L1Uqlfn5+eztZ9zd3Wk8\nYP6cZwJSK79LtYF//PEH89UxMTG//vprUlLSwIEDaUnnzp3r6up0/1611Aak1v+yNMnPz6fXswsn\nTZpEHq1/6ty5M5MFFxQU5KHVxIkTVb9CS0A6ePAg81JYWNgHH3zAzIrZ2treuHFDlyaAsAgmIP36\n66/+/v4TJkzYsmXLihUr+vTpM3XqVB33RwGlUnno0CEtMwd9+vRhr4dXKpWffPKJpnmChIQE5rLm\nBiT2ShpCyI4dOzj1VF3WSrm6up4/f177lzIt7dChA/u9Q4cOTUxMpI/ZC2Nb811qG/jRRx+xAz9D\nIpEcPXq0Wd+rltqA1PpfliYbN24khPTt25dd+NZbbxFC/ve//ymVSkJIamqqLh+liZaApFQq6V52\nHBYWFitXrmzNlwJvCSYgjRw5MjIyklku/tNPP/n5+bH/J4cm3bp1KyYmhjO04uzsvGrVKrVbrh09\nevTpp59m7rAWFha9e/fmLP5obkCqra1lJk7s7e3VrrnZvHkzO3ZaWlpGRkYyqzu1fCmjpKRk69at\nc+fOXbhw4e7du+/du7ds2TL6lgkTJujluzQ1MCsr66mnnmKHpfDw8DNnzjS3jWppCkjK1v2yNKEr\nnL799lt24eXLlwkhCxYskEgkH3zwgS6fo4X2gKRUKj/77DNmUwlCSNeuXffv39/KLwXeEsZu3zdv\n3hw0aNCiRYumTJlCSxQKRWBg4IQJE/7v//7PtHUTorKysrNnz1pYWPTq1YuZ0tDk7t2758+fb2xs\nDAgIaP0kvO5u3ryZm5trZ2fXrVs3Z2dnHd91/PhxuruMl5cX3auNev7553/77TdCyDvvvMP0llr5\nXVrU1NT89ddfEomke/fu7HFRQ3+vfn9Z5eXlVVVVqtmADx48OHjwYL9+/dip5IajVCoLCgqKiop6\n9OjRuXNnI3wjmIowAtLZs2ejo6PXr1/PPqV4+PDhfn5+a9euNWHFgG8+/vjjRYsWEULs7e1zcnJo\ndtbp06eDg4MfPHhgYWHx559/MttIAwCvCGOnBro4jnO6ia2tLWdxCcD48ePpCtY7d+707NmzX79+\nQUFBgwYNevDgASFk2rRpiEYAvCWMgES7cZzOHB1zNFGNgKfc3d1//fXX5557jhDS0NCQk5Nz5swZ\nhUJha2v71VdfbdiwwdQVBACNjLowtqqqKjIycurUqcx2kGwHDx6kuxfb29v7+fnFxsYyi+Po4pW6\nujr29XV1dezdXAAoX1/f33777fz58zk5OSUlJfb29p6eniEhIWoPhgAA/jBqQNq4cWNxcbHacbbE\nxMSUlJS2bdv26tWruro6LS1tz54969ato3/q0sDD3oqtoaHh2rVrgwYNMlbdQWB69erFLFsBAEEw\n+JBdQ0ODXC7PysqaN2+epgGT7OzslJSUrl277tu3LzU1de/evWvWrFEoFIsXL6YbVHfp0kUmk7F3\nf8nMzFQoFP379zd0/QEAwDgMHpAKCgoiIiJmzJixd+9eTdfs2rWLEDJ//nxmZUZERMTw4cNLS0uZ\n3ZEnTZp05syZxMTES5cupaenf/DBBx4eHuykO00O7Uy7NLbLT4umXfjtsD4aBAAABmHwITuZTMZk\nZp87d46u/ebIzs62srIaMmQIuzA8PDwjI+PkyZPh4eGEkMmTJ1dVVW3evDklJYUQ0qdPnw8//FDt\nsc0Mf3//PnYNC70UxJL0yNt79rMz/29anZbrAQDMwcWLF01dBfUMHpDs7e2Z05TZm4wxFApFWVmZ\nu7s7J7rQA2yuXbtGn1pZWc2bN2/27NnFxcVOTk7sxdtabNmypSjh4clsz3Rqe3H3Hy1uCK/4+/vz\n9p9Ua4i1XUS8TRNru4h4m+bv72/qKmhk+rTv6urqxsZG5ixOBi3hbGhvZWXVtWtXHaMRIcTa5fE2\nBIryossXL7eusgAAYCimD0j06DDVzhMt4Rws1lwSVw/2U/cHN1vzaQAAYDimD0j0lC3VwENLWnkG\nl7+//193H38CekgAYLb8/f35PF5H+BCQ1C56ZUqY89xa5uLFiwOeGcA8FU1AEuXQNhFvu4h4mybW\ndhExNu3ixYs8b5TpA5K9vX27du2KiooaGhrY5XQZbOv3YrAJGMg87nQ9p5WfBgAABmL6gEQI6dmz\n5/3793Nzc9mFOTk5hJDWL7aXuDyeRpJX1rfy0wAAwEB4EZCGDRtGCPnkk0+YkpKSku3bt1tZWYWF\nhbXmk/39/afGv8M8dXtwMzP/dms+EABAoPg/h2TUvew0iY6OTk1NPXXq1OjRo0eMGHHjxo2MjIza\n2trp06dzjjdtrosXLyrKigrjHk4juT0o//Xi5RCfIO3vAgAQHzqBxOeYxIuAJJFIkpOTly1bdvDg\nQTpwZ2dnFx8fHxsbq4cPfzLz2+1Bees/EwAA9M6oASk8PFxTjkfHjh3Xrl2rUCjkcrmdnZ1MJrOw\nsNDX99oGBNfmHqePL1+8/KK+PhcAAPSHFz0khkQioWdO6wvtnB6KfnxIqKKsSI+fDwAgFHwerKN4\nkdRgODTv/onM75Kz8grk2gGA2cE6JF5gZ367PbiJ5G8AAB4yj4D0ZF5DVl6lqWoCAACamEVAYu/5\n7fagXDQbCAEAiIlZBCROD0lRjrwGAADe4VeWnd7RrJKLFy+yM7/dcAgFAJgf/mfZiTwgqU0peVBW\nlJl/O8TH0fj1AQAwFf7v1GAWQ3bkyT2/n7l3QV7BPe0CAABMy1wCEifzO/l0qQkrAwAAqswmID2Z\n15CZX4ltvwEAeMVcAhIn89vtwc3k0yUmrA8AAHCIPKmBybJT3fM7Mx/LYwHAjPA5nYESeQ+JvXeT\nbUAwUz6g/oK8oh6jdgBgPrCXHY9wEu0IIUv3F5quOgAA8AQzCkjsHhJdGyuvRPI3AABfmFFA4uQ1\nPFN/QV5Rn4TUBgAAfjCjgCRx9eBMIxFCsCAJAIAnzCggcdBpJCxIAgDgCfMKSB2i4pnHzBarWJAE\nAMAH5rIOiT5VnUY6Je2JBUkAYA6wDsnEOHn3aqeRkNoAAOaA/+uQRN5D0o5OIxFClh4ojAmSmbYy\nYFZqampycnLq6+t9fX29vLysrKy0XFxbW5uTk1NVVeXl5dWtWzcbGxuj1RPAmETeQ1LlEBLFPGam\nkbBrAxjNxo0be/To0b59+5CQkOHDh/v6+nbr1m39+vUPHjxQvfi7777r27dv+/btBw8ePHLkyN69\ne8tkskWLFt28iUMmQYTMrofE3q+BmUYihCSfLsGRfWBQtbW1b7zxxpYtWzp27Dhu3Lhnn332zp07\nGRkZ2dnZM2fOPHbsWEpKioWFBb24oaFh4cKFq1atcnBwGDFixHPPPefk5HT58uUtW7Z8/PHH6enp\nv/76q7Ozs2lbBKBfZheQJK4eEhcPRXkRfcp0kpDaAIa2aNGiLVu2BAYG/u9//+vSpQstXLx4cVFR\nUXh4+NatWwMDA+fNm0fLlyxZsmrVqm7duu3fv9/X15f5kGXLlk2cODEtLW3ixIkZGRkmaAaAwZjd\nkB158myk0Xd/pQ/kFfXY2g4Mp6CgYP369R07dszMzGSiEeXh4bFlyxZCyKZNm2hJeXn56tWr7ezs\njh49yo5GhJA2bdps3bq1c+fOP//8c3FxsdHqD2AE5hiQ1E4jEUKSzpRgJgkMZMWKFQqFYvbs2e3a\ntVN9dcCAAYmJiaNGjbp79y4hZOXKlXfu3Jk6dapMpibXRiKRfPDBBxMnTiwoKDB4vQGMyEKpVJq6\nDobi7++vNsdRUVZUGDeAeTqx02I6jUQI8XKWFi4OVn0LQCs988wzp0+fvnTpUvfu3Zu8ePDgwUeP\nHs3Jyenbt68R6gZmRdONkQ/MsYfEWY0UZft4pE5eUT8l9YIpKgUid+nSJUtLS09PT10uvnz5MiFE\nl9AFICYiD0j+/v5qFyezc+0Ci34J8XFiniadxsAd6FlVVVVVVVX79u3btGnT5MV37ty5ceOGo6Oj\nvb09u7xXr15uT5o9e7bBqgwipOl+yB8iD0iaViY/eTZS+Zd9nkixm5J63uA1A3Pi4OBga2tbWVlZ\nX1/f5MW2trZSqbS6uvrevXvscjc3N49HOnbseP369cpKpIZCM/B/pwaRByRN2JvaEUJcS85uju7J\nPMXAHeiXhYWFj48PIUQul2u6JjEx8d///vfFixctLS179uzZ2NiYn5/PvuCXX345+ciaNWsMXWcA\n4zPTgMSZRqrLPRETJGMP3OFYCtCvgIAAQkhqaqqmCz777LMdO3a4ubkRQvr160cI2b59u6aLL1zA\nH0wgQmYakMiTR1EoyooIIZxOEpYlgR7Fx8cTQtatW3fnzh3VV3fv3n3r1q2goCA6b7Ro0aK2bduu\nWbPmxo0bqhfX19d/+umnhq4wgPGZb0Bij9opyovqck94OUuXDPNmCtFJAj3q379/ZGTkzZs3hwwZ\ncv36dfZL586de/PNNy0sLBITE2mJj4/P22+/fefOnSFDhly6dIl98e3bt7ECCcTK7LYOYtBRu9rc\n4/TprbSV7kt3Tg6SJZ0pkVc8nHmeknoey5JAX9avX3/79u0jR47861//euGFF5599llra+vff/99\n69at9+/f/89//hMaGspcvGTJkjt37nz++edBQUEREREDBw60sbE5d+7cnj17bt269emnny5atMiE\nbQEwBHNcGMuoPpJWum4ufSxx8fD+KpsQknS6hJ3RsGSYd0KEt/r3AzRTQ0PDsmXLkpOTr1y5whT6\n+fmtWbNm+PDhqtf/8MMPy5Yt+/vvv+le4BYWFr179/7222/79+8/ePBgLy8vuucQgO74vDDWrAMS\nZ8sGj6W7bAIG0hQ79l6riEkAIBp8DkjmO4dEVHLtbqWtJIR4OUs54WfJgULvxOPMOB4AABiCWQck\n8uRGq4qyIppuF+LjyM5uIITIK+pDv8pBTAIAMBxzD0jsPYQU5UXVmWn0cUKENzsLnDyKSUatHACA\nOTH3gCRx9WB3kqqPpDGPY4JkhYuDvZylTIm8oj70y7NGrR8AgNkw94BEOCtky4vqck8wT72cpUdm\nBnJ2cMCuQgAAhiDyLDv6oMmUkmsJY5kFSbYBwe5Ld7JfVZ1AQt4dqCWvqJdX1mflVWbm32YnatJ+\ndsIw75ggNQfuARiH7rdEUxF5QNLx585ekEQe5X+zL5BX1HsnHmeeejlLN0f3CvFx1FdVQQQy82+H\nftnELKOXszSmvwx/zYAJIe2b7zjhp+rI95wLvJylKtuBn0fSHTB0iUaEEHlFPVYRAGiCgEQIIRJX\nD/ZMEnsaiRETJGPnguOICmDoGI0YWEUAoJb57mXH4RASRRfGEpr/fSTNITSKc83kIBl7biAzv3Lp\n/sLJQTLyaJIATEheUZ98ukRe+fAu7+Uk9XSWejnbGGJklfNdSadL2K/Scbkhvk5eTg//VWTmVyaf\nLmXPKtGYdGRmIP7lADAwh/QYO7WB2dqOQ9Pftl7O0hAfp4Rh3ri/GN/S/YXsLXE5WpZQQNMTCCFM\nUKEfJa+oX3qgkBOB2GKCZJwVbIzM/NuckV5MRoLx8XkOCQHpsbrcE0UJY5inqqkNFCfBgQ1hyTiS\nTpckny4lhLD7HNrpkk1A+z2cBLlm0RKNmK/g/EHj5SzFjvJgTHwOSJhDeswmYCB7a7vSL+aqvYxz\nbBKbvKI+6XRJ6Fc5ONzPcKakXqC73zYrbNBsAi2/l8z8296Jx5ccKDRcNCLqVrbhKEgABuaQntAh\nKr424WHvhy6SVdtJSojw9nSWcmYFGPTeRy8zaG3Faun+QnllvbyifnJQZ844W+iXZ7UEDNoNYqZ2\n5BX1nIs1/V6am5XA+S4vJ+kQXycdR95oxiZ7R/mkMyWTg2ToVQNgyO4JirKiG+vmaVkkqxYd6qE3\nOzasn20u1bM/6CyLl5NUNS+ALcTHSTV6EQ2/GtpNYQJAc6ORXtYScQZ+Q3ycjsT1bc0HAuiIz0N2\nCEhcTS6S1UTtvQ8xSXfNTYZeMsx7iK8TIaTJrsnS/YWqMSmm/8PoxXkpJkg2OUim+pm0Ynrsx3Bq\ndSQuENkNYAR8DkgYsuOyCRgocfFQlBfRp1VHvtcxIDEHKbHvMhi70xHnoF7tmpucpvp7YYZVObTM\nA+l9SG1ykIydHDgl9TyyG8DMIamBi7NItjozTe06WU0SIrw5KQ9JZ0oy82/rrX5CJq+o5/yXmX97\n6f5Ci/jDnGjk5SzVFABoWlpzOxMJEd6cvdtV6ZKVoEdeztKEJ5daYy95MHPCG7I7d+7c66+/fujQ\nIVtbW+1XtrhnyjnaXNOaJC04ozFmu9yEDmMSQpqVS00nVDhDoF7OUi8nmxAfx9Z0NzXN9hHTDa5y\n0jSMHBTBDPF5yE5gAenGjRtvvfXWn3/+mZOTY2dnp/3i1vzcOTNJnWetVt24QTvVmMRMWlA0QWuI\nj6OBdhMwLS23fu1Up/flFfX6HS6jy5gy8yvpujEvJ6kJk9ywlzwYGQKSHuzYsSMlJSUvL6+xsZEQ\nYuiAxEm3a0EniaibS9dEZCtqdW84Gx3CMsMDGlSXWiMmgeHwOSAJZg6pT58+M2bMWLFixdixY43w\ndZyZJEV5kaZ1slpMDpKxl0BqwayoFfqGm3QipFnRiC40LlwcXLg42AyjEVHZS55g3hHMlWCy7Hr0\n6NGjRw9CSH19/c6dTa8Naj26cQPTSar7+5imdbKaqC6B1I6uwhHuehRNC3pigmR06ShnXzgjVo3v\naCRmMjvo+SbYehXMjWACkkl0mvUZk92guFlc+sXc5g7c0ZiUmV+ZpfIHr5eTUksOlwAAIABJREFU\nlO5HwA5XdAdxHYdrOO8N8XEy4f1L7TBdiI/T5uieuKvqIiZIdoWVjI7twMEM8S4g1dTUpKenM09d\nXFyGDh1qqsrQgTv2sRS30layh/J04eUsjXGWaRmM4kxrLzlQ6OksvVJRTwhhxm3klXVeTjZeztIh\nPo5XKurllfVqN5zmSaoYMePEwtbgrJei/zCwOAnMh6ECUlVVVWRk5NSpUydMmKD66sGDB9PT0/Py\n8uzt7f38/GJjYz09PelLFRUVK1asYK4MDAw0YUAihDiERNXlnmAG7qqPpDmERElcPfT4FXQnG/a0\nttolovKKepLPPXqHY8mBwqQzJcb8s1p1sx+CjXBagXPmFv3xIhEczIShkho2btxYXFx89+5d1ZcS\nExNnzZp1+PBhe3v76urqtLS0V1555ejRo/RVT0/PsyybNm0yUA11pJrdcGPdPL1/i+q0dovRP6uN\ns4G02mi0ZJg3olGL0X8J7L8nkk6XaP8rBEA09BmQGhoa5HJ5VlbWvHnzNmzYoPaa7OzslJSUrl27\n7tu3LzU1de/evWvWrFEoFIsXL66v52mCmU3AQIeQx4uQanOPVx9J0/u3xARpG9bThO5owOkP0X1x\nQr88a9CcPRr5OBuhHokLRL5yK9EeM7tk6YFCoadfAuhCn0N2BQUFI0eO1H7Nrl27CCHz5893c3Oj\nJREREcOHD8/IyDh27Fh4eLge66NHHaLi63JPMBvc3UpbaRMwUL8Dd4SQhGHeXk5Sur8Zs2aTEOL5\nKN5k5d+WV9TLK+voS+xbv+r6ysz8ytCvclq/L7VaqktnMGmkR0x+Jn2KySQwE/oMSDKZbO3atfTx\nuXPnNm7cqHpNdna2lZXVkCFD2IXh4eEZGRknT57Ue0Dy9/dnHrdmLRgduGP2bqDZDZ3fXN3a+j2J\nbs+qZdcALV0o+mc1Z3ME2lWSV9brdxJCNb0bx57qXUyQLCv/NjNYR8/xQ+8TWoZ9J+QzfQYke3v7\niIiIh59rreaTFQpFWVmZu7u7VPrEDdfHx4cQcu3aNV2+5dVXX3311Vd1rJIeFyQ7hEbV5h6vznw4\nWFedmWYbENzc/YR00eJ8BBrPhvg6TUk9z+4qJZ0uycyvbE3AoFugalpNhWhkIAnDvDPzK9nplwQ7\nx0OLsO+EfA5ORt2pobq6urGxsX379pxyWlJVVWXMyrRAh6h4icvjYbpbaSsVZUUmrI9aIT6OR2YG\ncnYcpyNs7ChF/+Jeur9wSuqFpNMa9wWgOy+Efpmj5WQ8RCMDUZ1M0n4KO4DQGXUdkkKhIOo6T7SE\nvspnxhm4az2mq8QeWKPzEIQQLycbQgg7wNBxIdozo6udyKPz6LTvMYH0bkPjTCYR9JNA1IwakKys\nrIi6wENL6Kv6RTunBh244+SF80eIj2Ph4mB2pgN9oClf6/Gr+Tp9Pg5KMA7ODg7k0Woz89yIFlqD\nz4N1lFEDko2NDSGkrq6OU05L6Kv6ZYhNbTkZd9VH0mwDgpu1x53R0DEfvezZSo/PmBwko6dmIJvO\nmBIivD2dpex+El0Blny6FDszge7o/ZDPYcmoc0j29vbt2rUrKipqaGhgl8vlckKITCaMP/fUbgTO\nw8kkisYktZuO08zyGB2OAqK7cSdEeHs5S0N8HBGNjC8mSHYkLpBTSDP7sS84iIax97Lr2bPnqVOn\ncnNzn3rqKaYwJyeHENKrVy8jV6bFHEKj6AQSfaooL7qWMNZ96U69r0zSCy9nKT2AVV5ZL6+oI4Rc\nqagf4uvEjivsV7PybzPrn0R5eKBA0TFYztYY2BccxMTYAWnYsGGnTp365JNPtmzZQktKSkq2b99u\nZWUVFham96/T+xwSg+7dwIlJLTjEz2gebuigIbqwX8XMBG/Rvy2STpew927AslnQEZ8H6yhjH9AX\nHR3t6+t76tSp0aNHb9q0afny5WPHjq2trZ06dWqXLl30/nUXL1400NmIElcPh5Ao24DHdwFFeVHh\nzAGG+C4AtpggGWcYlk4pmbBKIAiGux/qi7EDkkQiSU5OjoiI+Oeff1asWJGcnFxXVxcfHz9vnv53\nLDU0iatHp1mfsVcmISaBcajdgxVLlEDoLJRKpUm+WKFQyOVyOzs7mUxmYWFhiK8wztHxirKiawlj\nmaQ7QojExYPPY3cgGpwdBbGdIOjCODfGljF2D4khkUi6d+/epUsXA0Ujyt/f39DDphJXD/elO9FP\nAuPjnFoir6hPxkEVoJkR7oetZLKAZBzGGTNVG5NKv5hr6O8F4Jxawt77DoADc0jmgsYkdkl1Zhpi\nEhhBwrAnTiFZegAzSSBUCEh6I3H18P7yiakjxCQwAi9nKTpJIA4ISPqEmAQmgU4SiAMCkp4hJoHx\nqXaSTFgZgBYTeUAySVaJ2pjE7OkAYAicThLWyYIqZNmZmKmySlRj0q20lYhJYDiYSYImIcvOfElc\nPTrPeuLsvltpKzF2B4aDmSQQOgQkA3IIjeLEpOrMtMKZA3h7VgUIGqeTpOVkegB+QkAyLNWYRPcF\nr8s9YaoqgYixO0mEEOxuB8KCgGRwDqFR3l9mq+7jgCkl0DvOZkKZ+ZVJ2EwIhEPkAYknWSVq9xa6\nlbYSw3egdyE+TuyTKTCTBAye3A+1EHlA4k9WCY1J9Fg/Bh2+Q1cJ9MjLWZoQ8WR2AwbugBDCp/uh\nJiIPSLwicfXoEBXfISqeXYiuEuhdiI8ju5OUdKYEKeAgCAhIRkVjEmdKiaCrBPrGOZYCA3cgCAhI\nJkCH79R2lRCTQC+QAg5ChIBkGrSr5LF0F6erhOE70BdOCviU1POmqgmAjhCQTMkmYKDarhKG76D1\nVM+TRQo48BwCkokxXSV2IZPpgPWz0BpIAQdhEXlA4n/ePWUTMND7y2zbgGB2oaK8qChhTOkXczGC\nBy2jmgKOXcDNGf/vhxZKpdLUdTAUf39/nifdcyjKitQeVCFx8XAIjXIIiZK4eqh9I4AWoV+eZU5I\n8nKWbo7uFeLjaNoqgQnx+cYo8h6SsDBJ4apdpVtpKzGxBC3DTQHHOlngKwQk3lGbFE5YE0sIS9As\nXs7SJayMO2xwB7yFgMRTtKukJSyVfjEXKQ+go8lBMi9nKfMU2Q3AT5hD4jtNE0uUbUCwQ0iUQ2iU\n2lcBGJn5t0O/zGGexgTJ2EN5YD74fGNEQBIG7WEJWQ+gC3Z2AyHkSFwgshvMEJ9vjAhIQqI9LBFC\nHEKi2oeOswkYaMxagVDIK+q9E48zT72cpUdmBrKH8sAc8PnGiDkkIWHS8FTnlqjqzLSihDFIfAC1\nVPduwLIk4BX0kIRKUVZUl3uiOjOtNve42gswjgdqTUm9wM6yWzLMm714FkSPzzdGkQck+oC3P329\naHJ6ySZgYIeoeIQloOQV9aFf5TAnJGGprPng/y1R5AGJtz93vdPeYaJhCdNLQHEy7rycpYWLg7Vc\nD2LC5xsjApLYKMqKbqWtrM5MU/uqbUBwh6h4hCVYur9wCWs1ErLAzQefb4xIahAbiatH5zdXa0p8\nqM09jqwHIIRMDpI9ccz56RJsKQQmhx6SmOkyvYRxPLOlOpmUMMybfc4siBKfb4wISOJHw1L1kTRF\nufpjLJCPZ7ZUJ5Ni+suQdCdufL4xIiCZiybTxAk2IjJLnMkkQkiIj9ORuL6mqg8YGp9vjAhIZqf6\nSFpt7nFNWQ/k0VCebUAwIpOZUI1J2MRBxPh8Y0RAMlN0HK8u94SWDhMTmSSuHphnEjF5RX3y6RLV\nmIThO1Hi840RAcncNbk/HoVuk+hxchwo7OMgPny+MSIgASGPZpi0D+Ux0HMSK7UxCUuURIbPN0YE\nJHiCLrkPbBIXD4mrh7WLO41P1i7uSNUTNLXDd5ujeyIdXDT4fGNEQAL1mtVn4mCiFCGEBipCCPpS\nAsJJc8DeQmLC5xsjAhI0rbndJk0kLh6EEHasYp6iX8U3OGFWrPh8YxR5QKIPePvTFxxFWdGD8mu1\nuce1p+e1DA1XhBWxJK4eTAwjhCBuGRn7hFlsCi4C/L8lijwg8fbnLg40PinKimpzj9NAZYQvZcct\nQggTugirB0YvQABrJc4Js1gwKw58vjEiIIE+KcqKCCF1uScIITQ+0Yiladci49Cl70VYkQwYnMmk\nI3GB6CQJHZ9vjAhIYCRMd4qwYhXzmD+a7IGpfUmsVDdgRXaD0PH5xoiABKZHoxSNT/SxoryIKTR5\nB0tHTLgiTw4bcgoFF9WSTpdMSb3APEUnSej4fGNEQALBUBu3OOX0qSACGAefUxA52Q3oJAkan2+M\n1qauAICuHvYwdF7VRAMV0dz3InyKXg+D66PKsJd/0VhFm2ySDTImB3VmApK8oj4z/zY6SWAI6CEB\nEKKu+0U098AIK3KYBBOijBOf5BX1U1IvMDEJ6XaCxucbIwISQMtxOmHsEsIKWpou01dUo/HJITTK\nNiDYQMEJM0miwecbIwISgImx0zfIo0HFFqcgGuj8X066HTpJwsXnGyMCEgCvcVYf6zjpRXdkbx86\nTo8dJnSSxIHPN0YEJACBadYGGRIXj85vrtZLWOJ0krC7nUDx+caIgAQgbMwGg9VH0jR1nmwDgjtE\nxbc+LLE7Scj/Fig+3xiFFJCysrK+++67vLy8Dh06vPDCC5MmTWrTpo2W6/n8cwcwBC3n/9K5pQ5R\n8a35fM7udjgnSYj4fGO0NHUFdPXdd9/NmDHDxsYmNja2T58+a9asmT17toCiKYARSFw9OkTFe3+Z\n3XnWavbOEYQQRXnRrbSVhTMH0J0GW8bLWRri48Q8TT5d2vK6AqgQxsLYurq6Tz/9dMyYMcuXL6cl\nPj4+S5Ysyc3N7d27t2nrBsA3ElcPujipLvfErbSV7HE8RXlR6RdzW9NVSojwzvzy4YKkzPxKLJIF\nPRJGD+natWt37twZOXIkUxIYGEgIkcvlJqsTAL9JXD0cQqPcl+7kxB6mq8ReMqU7LydOJ6mktRUF\neEQYAcnNzW3nzp19+z5e9/DXX38RQry9vU1XKQABYAbx6OZ4DEV50bWEsS2ISV7OUnaXiNm+AaD1\nhBGQbG1t+/TpY2NjQ59evnx51apVAwYMCAgIMG3FAARB4urRadZnql2lawlj1WZAaDeZlchAt7bT\nQxUBeJhlV1NTk56ezjx1cXEZOnQo81SpVO7YseOjjz7q1q3bN9984+TkpO4zHuJzMgmASSjKiq4l\njOVkh3eIim/ulBJ7/2/s2iAsfL4xGiqpoaqqKjIycurUqRMmTFB99eDBg+np6Xl5efb29n5+frGx\nsZ6envSlioqKFStWMFcGBgYyAen69evvvPPOmTNnpk+fHhcXJ5FIDFR5ALGSuHq4L915K20lezfx\nW2krFWVFnd9crfvnILUBDMFQAWnjxo3FxcV3795VfSkxMTElJaVt27a9evWqrq5OS0vbs2fPunXr\nnnvuOUKIp6fn2bNnVd91+fLlmJiYLl26/PTTT926dTNQtQFEj84qSVw92IN11Zlpdbkn3Jfu1HEH\nPC8nKftpVl4lAhK0nj7nkBoaGuRyeVZW1rx58zZs2KD2muzs7JSUlK5du+7bty81NXXv3r1r1qxR\nKBSLFy+ur6/X9MmNjY3z58/39fXdtm0bohFAK9GYpHZKScc0By9nKXtJbNIZ5NqBHugzIBUUFERE\nRMyYMWPv3r2artm1axchZP78+W5ubrQkIiJi+PDhpaWlx44d0/SuEydO/PPPP88///ypU6eOsty+\njdlUgBai2XfsEhqTqo+kaXoLW8KwxzmuSG0AvdDnkJ1MJlu7di19fO7cuY0bN6pek52dbWVlNWTI\nEHZheHh4RkbGyZMnw8PD1X7yn3/+SQhhzy1RmzZtogN9mvj7+zOPeTuPB2AqElcP7y+z2WkOdJWS\noryoyTQHumsDk9qQfLoEo3a8xb4T8pk+A5K9vX1ERMTDz7VW88kKhaKsrMzd3V0qfWIA2sfHhxBy\n7do1TZ8cFxcXFxfXgiohCAFoR9McVGMSIaTJmMQ+2hwLkviMfSfkc3Ay6jqk6urqxsbG9u3bc8pp\nSVVVlTErA2AEd+/ezcnJ+fHHH7OysrT8yWVaNCZxws+ttJXXEsZqfyN7ywZ5RX0Sdm2A1jFqQFIo\nFERd54mW0FfBzP3www+WlpaWlpYzZ87UdM1PP/1Er5k+fTotOXDggKWl5UsvvcRcs3nzZktLS0dH\nx6IiNbP0u3fvtrS0nDt3Ln26detWSxXW1tZdu3Z96aWXfvvtN+aNaq9kv4W5sqqqKj4+XiaT9evX\nLzIyMiQkxMPDIzw8/NSpU638ERmC2jSH2tzj2ncYwl6roF9G3VzVysqKqAs8tIS+ql+0c4qBO2Gh\ni7VTU1M/++wzzugulZSURK9hlnUrH2F/iFKprKqqmjlzJnupNedb2E89PDz69evHlNTX1//666/7\n9u3bt29fSkrKxIkTNV3JsLCwoA/u3bs3atSorKysfv36jRs3rlOnTidPnvz5558PHz4cHh5+8uRJ\nfm4yQgMSOx2cpjloSQfHgiQB4fNgHWXUgET3/qmrq+OU0xJmZyA9QigSqPbt29++ffvHH3+Mjo7m\nvFReXp6enu7o6KhjjmVGRsa2bdvGjx/f5JUhISEpKSnskrt3765cuTIhIWHWrFmvvfYa0wFSvZIj\nJSUlKytrxIgRe/bsoX9pTZo0iRDywQcfvPfeezNnzvz11191qbzxdYiKdwiJ4kwpaYlJnAVJSG3g\nM3o/5HNYMuqQnb29fbt27YqKihoaGtjldNNumQwnfcFDo0aNsrGxSUpKUn1p+/btCoVCNVCp9dpr\nr9nY2MyZM6e8vLwF1bCzs5s/f76dnV1NTU2z/riho3wxMTGcfv8777zToUOH7Ozse/futaA+xkGn\nlNjHKWlZouTlLF3Cyv9GagO0hrE3V+3Zs+f9+/dzc3PZhTk5OYSQXr16GbkywFuOjo6RkZG//PJL\ncXEx56WkpCQXF5cRI0bo8jndu3d///33b9269dZbb7WsJra2tr6+voQQ7Rsncjx48ICoOx7F0tIy\nOTn5iy++0LIMnA9oTGJvEK4lJnH2WkVqA7SYsQPSsGHDCCGffPIJU1JSUrJ9+3YrK6uwsDC9f52/\nvz+f+6egRUxMTGNj45YtW9iFf/zxxx9//DF+/HjddzKcO3duUFDQ999//9NPP7WgGqWlpefPn/f3\n9+/SpYvu76Ir6pYsWbJq1aqKigr2SyNGjJg+fbpqrinf0A3CdYlJnNSGpQcKjVRFaCb+3w+NfWJs\ndHR0amrqqVOnRo8ePWLEiBs3bmRkZNTW1k6fPr1Z/8PrCHNIwhUeHu7h4ZGUlLRo0SKmkA7ixcTE\nlJbqmtBlZWX17bffBgYGxsXFDRkyREskqK+vZ0b2lEpldXX12bNnly5dam1t/e6777Kv/P333+fP\nn6/6CZGRkYMGDSKETJ069dSpUxs2bIiPj1+4cGH//v0HDhwYEhISFhbWrl07HWtucjQm3Vg3rzb3\nOC1RlBfdWDfPfelOzpXs1Aa6awNmkniI/3NIxg5IEokkOTl52bJlBw8epAN3dnZ28fHxsbGxRq4J\n8JylpeXEiROXL19+8uTJZ599lhCiUCi2bdv29NNP/+tf/9I9IBFCevfu/c477yxdunTBggWadlkk\nhOzYsWPHjh2cQhcXl1OnTvXu3ZtdeP78+fPnz6t+gpeXFw1IFhYWX3/99fTp0zds2HDgwIHs7Ozs\n7OzVq1e3adNmypQp//3vfw3x55chqMak2tzjpV/M5WwNHuLjyN61Yen+whAcSAHNZ6iAFB4erql3\n0rFjx7Vr1yoUCrlcbmdnJ5PJmGRZALbJkycvX748KSmJBqT09PSbN29yOis6euedd3bt2vXNN9+8\n9tproaGhaq/x8/Njv3T//v2zZ8/+8ccfY8aMOXz4MLP7IiEkMjLy008/Vf2Ejh07sp/279+/f//+\nhJCrV68eP3587969e/bs+frrrw8dOnT8+HEXF5cWNMT4aExi591VZ6bRdUvsyzi7NqCTBC1g7B4S\nQyKRdO/e3dDfgnVIgubn5xccHPz999+vXr1aKpUmJSVJJJJ///vfLfioNm3abNq0KTg4ePr06efO\nnVN7zYABA9avX88uUSqVsbGxmzdvTkxM/PLLL5lye3v7Zu0637Vr165du0ZHR5eUlERGRmZnZ3/2\n2WfLly9vQUNMguY4FMYNYEqqD6faBgTbBAxkStjTSAT537zE58E6ShhHmLfYxYsXEY0ELSYmhi5I\nunHjxt69e0eMGNHijsUzzzwzd+7c/Pz89957T8e3WFhYxMfHE0KatWzo6tWrTk5Ow4cPV31JJpP9\n5z//IYScPHlS9w/kA4mrR+dZj4fpFDeLS7+Yy75ANf9bXsHrTEIzxP/7ocgDEgjduHHj6IKkbdu2\nPXjwICYmpjWf9v777/v6+q5ZsyY7O7vpqwkhhPTs2bNt27bXr1/X/Vu6du1qbW2dlZVVU1Oj+qql\npSUhpEOHDrp/IE84hEaxh+mYPVgZnPxvpNtBcyEgAa85ODiMHj36l19++fzzz11cXNi71bWAjY3N\nN99809jYuGrVKh3fYmlpKZPJqqqqOKu5tXvllVfq6+vHjRvH2Zfk3r17X3zxBXmUFy44DiFR7ETw\n6iNpdbknmKfcU/tOl2BNEjSLyeaQAHQUExOzbds2uVw+d+5c3ZcfaRISEjJjxoyvv/5a97c4ODg0\nNjZevny5R48etCQrK2v06NFqL962bZuNjc3nn39+9uzZffv2eXp6Tp482d/fv23btnl5eampqZcu\nXRoxYsTUqVNb2RCToAkOzGSSoryo9Iu53l897m4mDPNmD9YtPVAY4uPk5axmQ0IAVQhIwHdhYWEe\nHh5FRUWtHK9jrFixIiMjQ/fDIAICAs6dO/fdd98tXbqUlly9evXq1atqL6Z7NNja2h4+fHjt2rVr\n1qxh5+O5u7svW7Zs4cKFbdq0aV0jTIbm1zGDdXTgjhnK83KWJgzznpJ6gT6VV9RPSb1wBCngoBsL\nzp7HYsKklPB8Hg9AWBRlReyVSRIXD3YniRAyJfUCe7Buc3RP9lAemAr/b4kin0Pif1YJgOBwFiHR\ngTv2BQnDvNnDdMhu4An+3w9FHpAAwBBsAgY6hEQxT+tyT7D3uPNylm6OfrxXMnZcBR0hIAFAS2hP\nAaebCTFP0UkCXSAgAUBLSFw9OJ0kdgo4ISQh4vE6WXSSQBcISADQQugkgX6JPCDx//wPAOHi7CdU\nm3scnSQ+4//9UOQBif9ZJQCCZhMwkH3YOTpJfMb/+6HIAxIAGBQnBbzJTlJm/m3jVQ6EBgEJAFql\nWZ2kKalqDjYEoBCQAKBVmttJwkwSaIKABACthZkk0AsEJABoLXSSQC8QkABAD9BJgtYTeUDif949\ngDg0t5O0dD9ikrHx/34o8oDE/7x7ANFoVicp6QxG7YyN//dDkQckADCaJjtJm6N7Mo/p2X3GqxwI\nAQISAOiNTcBA24Bg5innnCQvZym7k5SZX2m8moEQICABgN5wtgBXlBdVH0ljX4BOEmiBgAQA+uQQ\nGsXuJHFmkrycpezjzNFJAjYEJADQsyYPOGceo5MEbAhIAKBnnJkkTmoDOkmgCQISAOhfp1mfMY/R\nSQIdISABgP6pHnCuKCtinqKTBGqJPCDxf2UygFhpP+Cc00nC7nZGwP/7ocgDEv9XJgOIlWoniT2Z\nxFmThN3tjID/90ORByQAMCFOJ6nqyPfsV3GYLHAgIAGAoXA2E+J0krhbgGO7VbOHgAQABsTZuIEz\nkzQ5qDPzODO/Ep0kM4eABAAGJHH16DxrNfOUs+NqiI+Tl7OUeZqM1AbzhoAEAIal5VgKL2dpTH/k\nf8NDCEgAYFjaj6WYzFqQhIP7zBwCEgAYHKeTxN64gbNIFgf3mTMEJAAwOE4nSVFepKWThNQGs4WA\nBADGoOXsPk7+N1IbzBYCEgAYg/az+0J8HJnHSG0wWwhIAGAkDqFRHkt3MU9L181lYlJChDfTSZJX\n1Id+edYE9QNTE3lA4v9mggBmhTNwx04BZy+SlVfWySvqjVozM8D/+6HIAxL/NxMEMDec7Aamk8Se\nRpJX1GMmSe/4fz8UeUACAL7R1ElC/jcgIAGAsWnqJCH/28whIAGAsWnqJGH/bzOHgAQAJtBp1mfM\n4yc7Sdj/23whIAGACUhcPTR0kpzYlyG1wawgIAGAaaidSfJyli4Z9vgkWSySNSsISABgGppmkrD/\nt9lCQAIAk9HUSUL+t3lCQAIAk9HUSUpgjdqhk2Q+EJAAwJTUHkvh5SxlZzcsOVCImGQOEJAAwJQ0\nHUuREOHNvgwxyRwIKSBlZGTExMSEhoaOGzcuOTlZoVCYukYAoAdqO0khPo7sdDtCyJIDhVNSLxi7\ncmBEgglIycnJ8fHxbm5u06ZNc3Nz+/jjj5cuXWrqSgGAHli7uKufSYrw5sSkpNMl3onHjVo5MCJh\nBKSGhobVq1fHxsYmJiaOHz9+1apVL7zwwu7du+vrsUE9gOBxz+4re3zAeUKE9+bonuyL5RX13onH\ncTiFKAkjIJWUlHTp0mXo0KFMyVNPPUUIuX//vukqBQB6YxMwkHmsKC+qOvI98zQmSFa4OJh9sbyi\nPvSrHMQk8RFGQHJ3d8/IyOjbty99mp+fv3Xr1uHDhzs4OJi2YgCgF5xOEtNDorycpYWLg72cpUyJ\nvKIe80niI4yAxDh58uSgQYNeeuklJyenFStWmLo6AKA3ms6koLycpUdmBrJjErZeFR8LpVJp6jo8\noaamJj09nXnq4uLCHqm7efPmuXPn5HL5t99+6+XltW7duvbt22v6KH9/f54fjwgAbNcSxtbmPsxZ\nkLh4eH+VzbmAM1gXEyTjzDBBk/h8Y7Q20OdWVVVFRkZOnTp1woQJqq8ePHgwPT09Ly/P3t7ez88v\nNjbW09OTvlRRUcHu+gQGBrIDUseOHcPCwgghwcHBo0aNOnTo0OjRow0/ogT4AAAWE0lEQVTUBAAw\nsg5R8bUJDwMSzf9mzy0RuqtQf9mSAw8XJGHrVZEx1JDdxo0bi4uL7969q/pSYmLirFmzDh8+bG9v\nX11dnZaW9sorrxw9epS+6unpeZZl06ZNhJCDBw/GxsZWV1czH9KjRw9nZ+ecnBwD1R8AjE/TTkJs\nOFVWxPQZkBoaGuRyeVZW1rx58zZs2KD2muzs7JSUlK5du+7bty81NXXv3r1r1qxRKBSLFy/WksNt\nbW199OjRAwcOMCWlpaWVlZVMvwoAxIGd2lCbe5yT3UBUdhXC9g1ios8hu4KCgpEjR2q/ZteuXYSQ\n+fPnu7m50ZKIiIjhw4dnZGQcO3YsPDxc7buCg4N9fX1Xr15dUVERERFx7dq1jz/+uH379k1+HQAI\ni03AQImLh6K8iD6tOvI9Z9SOEDI5qDMzWCevrDNq/cCQ9NlDkslkax+ZNm2a2muys7OtrKyGDBnC\nLqRx6OTJk5o+uU2bNuvXr+/Ro8eqVauGDRs2derUdu3affvttzKZTNNbKH+W5jcIAIxN4urBTrer\nyz2hKCviXMPuIckr6pNwqmxThHIn1GcPyd7ePiIi4uHnWqv5ZIVCUVZW5u7uLpVK2eU+Pj6EkGvX\nrmn5cA8Pj40bN9bX15eWlrq6utra2upSJd4mkwCAJpxFsrfSVnZ+czX7Ajpqx3SSlh4oZJ+fBKrY\nd0I+xySjrkOqrq5ubGxUTdSmJVVVVU1+glQq9fLy0jEaAYAQaV8kS7H3Akdqg2gYNSDR/blVO0+0\nBLt3AwClfZEsIcTLCakNImTUgGRlZUXUBR5aQl/VL/6PmQKAKomrh/b8by9naYiPI/M0M78SMalJ\n/L8fGjUg2djYEELq6rhZMbSEvqpfFy9exDQSgBBxOkmqMWlykIy9k1DSmRJst6od/++HRg1I9vb2\n7dq1KyoqamhoYJfL5XJCSJMpcwBgPjiLZKuPpHHS7bycpZujezFP6a5CxqsfGICxN1ft2bPn/fv3\nc3Nz2YV0w4VevXppeBMAmKNOsz5jHqvtJHFOlZVX1GPgTtCMHZCGDRtGCPnkk0+YkpKSku3bt1tZ\nWdFN6vSL/2OmAKCJarqdasbd5CAZO7sh6UwJMu404f/90NgBKTo62tfX99SpU6NHj960adPy5cvH\njh1bW1s7derULl266P3r+D9mCgBaNDmT5OUsZW/4La+oT8Y6WQ34fz80dkCSSCTJyckRERH//PPP\nihUrkpOT6+rq4uPj582bZ+SaAAD/SVw9Os96vCpW0+527IWx2AJcuEx2HpJCoZDL5XZ2djKZzMLC\nwhBfwedjPwBAR4qyomsJY5nd7TSdk+SdeJx5ujm6J/Zu0ITPN0aTnRgrkUi6d+/epUsXA0Ujiv9j\npgCgHWd3O0V5UekXcznXcLcAP4DUBjX4fz8U2BHmzcX/MVMAaJJDaNQTKeCZaaoDd9hMqEn8vx+K\nPCABgDiwU8AJIaqdpBAfR2wmJHQISAAgAJzsBrUDd5ODOjOP5ZV16CQJDgISAAgDZ+BOdVkS55wk\n5H8LDgISAAgGZ++GqiPfs19F/rfQiTwg8T+rBAB01+TeDQnDkNqgEf/vhyIPSPzPKgGAZtG+dwM3\n/xupDSz8vx+KPCABgMg0uXcDJ7XBeDWDVkNAAgCBsQkYKHHxYJ5yOkmc1IYkpDYIBwISAAgMZ+8G\nTieJM2qXfLrUqJWDVkBAAgDh4Rzfx1mTxN61ITO/EqkNQiHygMT/rBIAaAFOuh0nu8HLScq+OCsP\n+d+ECOF+KPKAxP+sEgBoGe4Gd6wzzjkLkpLOYBqJECHcD0UekABAxLSccT6ZFZCwIEkoEJAAQKi0\nrJPFXqtChIAEAAKmZZ0s9loVHAQkABAwLetksdeq4CAgAYCwcdbJMingXs7SJcOeyP+WV9Qbu3LQ\nHAhIACBsqmecM50kTmoDjjbnOZEHJP7n3QNA62laJ4sDKdj4fz8UeUDif949ALSeaiep+kgafcw5\nkMKc0+34fz8UeUACADPB6SQx6Xacre2wSJbPEJAAQCQ0dZI2R/dkyuUV9VNSLxi7ZqAbBCQAEAld\nO0mnS7AmiZ8QkABAPDSl27E7SYSQKannjVot0A0CEgCIh5Z0Owzc8R8CEgCIiqZOUkyQjD1wh3OS\neAgBCQBExdrFXe1MElHJbjDnFHB+QkACAFHhbAHO3t2OM3CHThLfiDwg8X9lMgDonU3AQPbTqiPf\nM49DfJy8nB+fJ2tWnST+3w9FHpD4vzIZAPSOs3FDXe4J9mGyCU/uuGo+nST+3w9FHpAAwDyxR+04\n5ySZcyeJ5xCQAECEVA+TZR6bcyeJ5xCQAECcNO0kRJ48u4+gk8QbCEgAIE4SVw9N+d/cdbI44Jwf\nEJAAQLQ0LZIlOOCclxCQAEC0NO0kRLAmiZcQkABAzNBJEhAEJAAQM01nUhB1B5zLK+qNWjl4EgIS\nAIicpp2EiOoB5weQbmdKCEgAIHI2AQMlLh7MU3SSeAsBCQBEjrOTEDpJvIWABADih06SIIg8IPF/\nd1sAMAJ0kogQ7ociD0j8390WAIxDeydpCSsmJZ0uEeWaJP7fD0UekAAAKO2dpMlBMvYW4FNSzxu1\nckAIQUACAPOhfU0SZ+AuCetkjQ4BCQDMhZbTzQkhIT5O7L0bxDqTxGcISABgRproJEU80UmaknrB\nqJUzewhIAGBGmuokObI7Sdhx1cgQkADAvGjpJBFCnjgnqaIeZ/cZEwISAJgX7Z0kTgp4Zn4lshuM\nBgEJAMyO9k4SJwUc2Q1Gg4AEAGanyU4SJwUc2Q3GgYAEAOZIeycpJkjGZDdYW1oM8XE0auXMFQIS\nAJgj7Z0kwspueNCofM4bAckYBBmQ9u/fP3DgwNLSUlNXBAAETHsnyctZWrg4uOi9QbkLB/h2tDF6\n7cyR8AJSaWnpu+++W1FR0djYaOq6AICANdlJ8nKWuju27dXJzuhVM1MCC0iNjY0LFix49tlnTV0R\nABAD7Z0kMDKBBaSNGzdaWFhMmjTJ1BUBADFospMExiSkgPT3339/++23H330kaWlkKoNAHyGThJ/\nCObOXldXFx8f/3//939dunQxdV0AQDzQSeIPa1NXgKumpiY9PZ156uLiMnToUELIBx980L1798jI\nSNNVDQDEiXaSanOP06elX8z1/irbtFUyT4YKSFVVVZGRkVOnTp0wYYLqqwcPHkxPT8/Ly7O3t/fz\n84uNjfX09KQvVVRUrFixgrkyMDBw6NChhw8f3r9//7Zt28rLywkhlZWVhJBbt27Z2to6Oprj+gB/\nf3+eH0XcMmJtFxFv08TRLtpJYgKSoryoLvfE06NjRNA0YTFUQNq4cWNxcfHdu3dVX0pMTExJSWnb\ntm2vXr2qq6vT0tL27Nmzbt265557jhDi6el59uxZzlv++OOPmpqaV155hV04duzYoKCgrVu3GqgJ\nAGA+VDtJpq2PedJnQGpoaCgqKrpy5cqPP/64d+9etddkZ2enpKR07do1KSnJzc2NELJ///558+Yt\nXrx4//79UqlU7bsmTZo0cuRI5mlubu6iRYs2bdrE9KsAAFpD4urRISq+NuFxJyncSf3tCAxHnwGp\noKCAHTbU2rVrFyFk/vz5NBoRQiIiIoYPH56RkXHs2LHw8HC17+rYsWPHjh2ZpzU1NYSQbt26IcEB\nAPSF3UmyDQj+6+Kfpq6R2dFnQJLJZGvXrqWPz507t3HjRtVrsrOzrayshgwZwi4MDw/PyMg4efKk\npoDUYv7+/vr9QP4Qa9PE2i4i3qaJqV197Bpe62S1+pq07K+/CLEUU9MEQZ8Byd7ePiIi4uHnWqv5\nZIVCUVZW5u7uzhma8/HxIYRcu3ZNxy/q16+fLpONmJAEgBYYY+oKmC2jrkOqrq5ubGxs3749p5yW\nVFVVGbMyAADAK0YNSAqFgqjrPNES+ioAAJgnowYkKysroi7w0BL6KgAAmCejBiQbGxtCSF1dHaec\nltBXAQDAPBk1INnb27dr166oqKihoYFdLpfLCSEymcyYlQEAAF4x9uaqPXv2vH//fm5uLrswJyeH\nENKrVy8jVwYAAPjD2AFp2LBhhJBPPvmEKSkpKdm+fbuVlVVYWJiRKwMAAPxh7N2+o6OjU1NTT506\nNXr06BEjRty4cSMjI6O2tnb69OnYdgEAwJwZOyBJJJLk5ORly5YdPHiQDtzZ2dnFx8fHxsYauSYA\nAMArFkql0iRfrFAo5HK5nZ2dTCazsLAwSR0AAIA/TBaQAAAA2Hh3YmzraTn9T1h+/PHHrKysS5cu\nOTg49O7de9q0aZ06deJcI+jG1tXVTZw40dnZecOGDZyXhNiu69evHzx48OTJk4WFhR4eHmPHjqUp\nPGxCbJdSqTxw4MAPP/xw5coVR0dHX1/fmJgYuv8km1Ca1uKzQ3W8wIS0N63J+wkfmia2HhLn9L+C\ngoK2bdsyp/8JRUNDw1tvvXXo0CF7e/uAgIDS0tIrV65IpdJNmzb179+fuUzojX3vvffS0tICAgJ2\n797NLhdiu65evfraa6/dvHmzU6dOHTp0OH/+PCEkLi5uzpw5zDVCbBchZOHChXv27LGzs+vdu/eV\nK1dKS0utra1Xr179wgsvMNcIqGkrV67csGHD22+//frrr3NearIVPG+mpqbpcj/hS9OUInLy5Ek/\nP7+hQ4deu3aNlvz88889e/Z8/vnn6+rqTFu3Ztm2bRv9C4Wpdlpamp+f3+DBg+/fv09LhN7YQ4cO\n9erV66mnnoqMjGSXC7Fdd+7cCQ0N7d2799GjRxsbG5VKZUFBQf/+/Xv06MG0QojtUiqVhw4d8vPz\nGzVqVGVlpVKpbGxs3Lt3r5+f38CBAx88eECv4X/THjx4UFhYmJmZOXfuXD8/Pz8/v/Xr13OuabIV\n/GymLk1r8n7Cn6YZex2SQWk6/a+0tPTYsWMmrVrzJCUlWVtbL1++nDmn49VXXx08ePCNGzcuXbpE\nSwTd2Js3by5evHjmzJnOzs6cl4TYrgMHDhQXFy9YsGDQoEE0Q8fb2zs2NtbJyemvv/6i1wixXYSQ\n33//nRASGxvr6OhICLGwsHjxxRd9fX1v3bpVUFBAr+F/0woKCiIiImbMmKHpJGuiQyv42Uxdmtbk\n/YQ/TRNVQNJ0+h8h5OTJkyaqVLMplcrr16/7+/u7urqyy729vQnr1ChBN/bdd9+VyWQzZ85UfUmI\n7UpPT7eysho9ejS78I033jh+/Pjw4cPpUyG2ixCi+hcDIaShocHS0pI5xJn/TaNnh1LTpk1Te02T\nreBnM5tsmi73E/40TTxJDfo6/c/kGhoali5dqpq/kJeXRwjx8vIiAm9samrqsWPHdu/erbq/u0Db\nderUqaeeesre3v7y5ctnz54tLi729fUdOHAgc8sWaLsIIUOHDl25cmVSUlJYWJidnR0h5MCBA4WF\nhc8++6yTkxMRSNNaf3Yob5vZZNOavJ/wqmniCUiiOf3P2tp6zBjukZXHjx8/ceKEr6+vr68vEXJj\n5XL5xx9/PGfOnO7du6u+KsR21dTU3L9/XyaTbdy4kb0nloODw4cffjh06FAizHZRnp6eSUlJs2fP\nDgsL69evX0FBQWFh4aBBg1avXk0vEG7T2JpshXCb2eT95Pbt2/xpmniG7ER8+t8vv/wSFxfXtm3b\n999/n32mlOAa29DQsGDBgh49ekydOlXtBUJs161btwghJ0+e/Pzzz997772jR48ePXp08eLF9fX1\nb7/9dnFxMRFmuxg1NTWNjY23b9/OzMwsLCwkhNTW1lZWVtJXBd00RpOtEEczKc79hFdNE09AEuXp\nf9XV1f/5z3/efPNNOzu7TZs2BQYG0nKBNnbdunV5eXkfffSRpaX6f3hCbBetW0VFxYIFCyZMmODi\n4uLi4jJp0qQ5c+bcu3dv48aNRJjtorZu3RoXF+fi4pKSknLu3LnTp08vXrz4/Pnzo0aNunjxIhFy\n09iabIU4mqn2fsKrpoknIInv9L/Dhw+/9NJLu3fvHjFiRHp6OnsFkhAb+9dff61fv/6NN95wdnau\neaSxsbGxsbGmpqa2tpYIs100/YwQEhkZyS4fOXIkIeTChQtEmO2idu3aJZFINmzYMGDAAGtrawcH\nh0mTJr399tt1dXW//PILEXLT2JpshQiaqel+wqumiScgiez0v82bN8+cObNNmzbJycmrVq2iE8gM\nITb277//bmhoWLVqVX+W0tLSCxcu9O/f/7XXXiPCbJejo6NEIrG1taVz/gz6K6uoqCDCbBchpLKy\n8sKFC35+fkw2MEV3oKD5VwJtGkeTrRB6M7XcT3jVNPEEJCKi0/8OHz780Ucf9e/f/6effnr22WfV\nXiO4xgYEBMSpsLe3d3FxiYuLi4qKopcJrl0SieTZZ5+tra2l00UMmsXEbL4iuHYRQuzs7KysrGpq\najjlt2/fJo8iLhFm01Q12QrhNrPJ+wl/miaqgCSO0/8aGho++uijdu3aff311/b29pouE1xjn3rq\nqTkqHBwcXF1d58yZM378eHqZ4NpFHtV5zZo1ykcbcTU2Nq5bt448GrgjwmxXmzZtevXqdfXq1f37\n9zOFSqVy/fr1hJB+/frREiE2TVWTrRBoM3W5n/CnaeJJ+yZiOf0vLy/vypUrXbt2/eyzz1RfnTJl\niru7OxFLY1UJsV1jxozJzMzcs2dPaWkp/X/7559/Pn36dL9+/V5++WV6jRDbRQhZsmRJdHT07Nmz\ne/fu/eKLL7Zp02bv3r1nz57t06cP8zeEQJvG0WQrBNpMXe4n/Gma2DZXvXnzJj39j46H2tnZvfHG\nG7GxsUJJgyGE7Nix491339X0ampqat++feljETQ2NDTUycmJs7mqENt17969Dz/88Ndff6UDd87O\nziNHjly4cKFEImGuEWK7CCF//vnnypUrs7Oz6dM2bdpER0e/+eab7JUrAmraoUOH4uLi1G6u2mQr\neN5MtU3T8X7Ck6aJLSBRZnX6n1gbK9B2Xb9+3cLCQstUsEDbVVdXd+3atbZt27q5uWm6SQm0aRxN\ntkIczVTL5E0TZ0ACAADBEVVSAwAACBcCEgAA8AICEgAA8AICEgAA8AICEgAA8AICEgAA8AICEgAA\n8AICEgAA8AICEgAA8AICEgAA8AICEgAA8AICEgAA8AICEgAA8AICEgAA8AICEgAA8AICEoAx3L9/\nv7S0FMePAWiBgARgEBs2bHj++eePHj2ak5Mzfvz4/v37DxkyJDAw8N13362trTV17QD4CAEJwCBq\nampu3LiRnZ09bdq0mzdvRkRE9O/fv66ubseOHe+8846pawfAR9amrgCAmG3YsGHmzJmzZ8+2tLQk\nhPzxxx/jxo3bv39/WVmZq6urqWsHwC/oIQEYUI8ePebMmUOjESHk6aef7tWrV2Nj49WrV01bMQAe\nQkACMKDnnnvOwsKCXeLr60sIuX37tolqBMBfCEgABtSpUydTVwFAMBCQAAyI0z0CAC0QkAAAgBcQ\nkAAAgBcQkAAAgBcQkAAAgBcQkAAAgBcssNsjAADwAXpIAADACwhIAADACwhIAADACwhIAADAC/8f\nHKtogcetVgoAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "kappa = 1e4;\n", "lam = linspace(1,kappa,m)';\n", "A = sparse(diag(lam));\n", "b = randn(m,1);\n", "b = b/norm(b);\n", "[xCG,~,~,~,resnorm] = pcg(A,b,1e-14,100);\n", "semilogy(resnorm,'.-')\n", "hold on\n", "text(60,0.3,'CG')\n", "[xMR,~,~,~,resnorm] = gmres(A,b,[],1e-14,100);\n", "semilogy(resnorm,'.-')\n", "text(40,.0075,'MINRES')\n", "ylim([1e-4 1])\n", "xlabel('n')\n", "title('Convergence for \\kappa=10^4 ')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "CG and MINRES are often not much different, practically speaking. MINRES also works for indefinite problems, while CG has considerably more (or at least better-known) theory behind it." ] } ], "metadata": { "kernelspec": { "display_name": "Matlab", "language": "matlab", "name": "matlab" }, "language_info": { "codemirror_mode": "octave", "file_extension": ".m", "help_links": [ { "text": "MetaKernel Magics", "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" } ], "mimetype": "text/x-octave", "name": "matlab", "version": "0.11.0" } }, "nbformat": 4, "nbformat_minor": 1 }