M348-Solving Linear Systems by Fixed Point algorithms for linear systems.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fixed point iterations\n",
    "\n",
    "by [Richard Tsai](https://www.researchgate.net/profile/Richard_Tsai3)\n",
    "\n",
    "$f:\\mathbb{R}^d\\mapsto\\mathbb{R}^d$\n",
    "\n",
    "* Need to come up with $g(x)$ such that \n",
    "$$ f(x^*)=0 \\iff x^* = g(x^*) $$\n",
    "\n",
    "* Introduce the fixed point iteration for $g$: $$ x_{k+1} := g(x_k).$$\n",
    "\n",
    "**What should $g(x)$ be?** \n",
    "So that $x_k\\longrightarrow x^*$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear systems of equations\n",
    "Let $$f(x)=Ax-b,$$\n",
    "where $A$ is an $d\\times d$ invertible matrix, and $x,b$ are two vectors in $\\mathbb{R}^d$.\n",
    "We want to find the solution of $f(x)=0$.\n",
    "\n",
    "Typically, we write the above linear system of equations as \n",
    "$$ Ax =b $$\n",
    "\n",
    "We will consider fixed point iterations with $g(x)=Gx+c$.\n",
    "\n",
    "### General strategy for creating fixed-point iterations for systems by splitting the matrix\n",
    "\n",
    "A general strategy in coming up with fixed point iterations is to *split* $A$ into the sum of two matrices, separate the equations accordingly by the equal sign\n",
    "$$\n",
    "A=M-N \\implies Ax=b \\iff Mx=Nx+b\n",
    "$$\n",
    "and define the iterations\n",
    "$$\n",
    "Mx^{(k+1)} = Nx^{(k)} +b, \n",
    "$$\n",
    "\n",
    "leading to the general form\n",
    "\n",
    "$$x^{(k+1)}=M^{-1}Nx^{(k)}+M^{-1}b \\iff x^{(k+1)}=G x^{(k)}+ c $$\n",
    "\n",
    "**Of course, $M$ has to be very simple and easy to invert.**\n",
    "\n",
    "**Both iterative schemes can be used to solve *nonlinear* equations**\n",
    "\n",
    "** Learniing objectives **:\n",
    "* Know how convergence of a fixed point iterations is related to the eigenvalues of the matrices involved\n",
    "* Know in which situtations does one expect that iterative methods be better than direct methods\n",
    "\n",
    "\n",
    "## Three fixed point iterations to solve linear systems of equations\n",
    "\n",
    "We first discuss three particular fixed-point iterations (two different $g$'s):\n",
    "\n",
    "* Jacobi iterations (J)\n",
    "* Gauss-Seidel iterations (G-S)\n",
    "* Successive Over Relaxations (SOR): this method builds on top of (G-S)\n",
    "\n",
    "**Learning objectives**:\n",
    "\n",
    "* Know how to write Julia code for Jacobi and for Gauss-Seidel methods\n",
    "* Know how convergence of the iterative methods is related to the spectral radius of $G$\n",
    "* Know how to use the Jacobi and Gauss-Seidel **ideas** to decouple nonlinear systems of equations into solving smaller ones at a time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Jacobi iterations\n",
    "\n",
    "In the Jacobi iterations, $M=D$ is the diagonal of $A$ (and $N$ is $A$ with its diagonal culled out): $$G=G_J:=D^{-1}(L+U)$$ \n",
    "\n",
    "$$x^{(k+1)}=G_{J} x^{(k)}+ c_{J}$$\n",
    "\n",
    "* **Very easy parallelization**\n",
    "\n",
    "\n",
    "$$ \n",
    "a_{11} \\mathbf{x_1^{(k+1)}} + a_{12} x_2^{(k)}+a_{13} x_3^{(k)}+...+a_{1n} x_n^{(k)}=b_1 \n",
    "$$\n",
    "\n",
    "$$ \n",
    "a_{21} x_1^{(k)} + a_{22} \\mathbf{x_2^{(k+1)}}+a_{13} x_3^{(k)}+...+a_{1n} x_n^{(k)}=b_2 \n",
    "$$\n",
    "\n",
    "$$\\vdots$$\n",
    "\n",
    "**Jacobi method converges for diagonally dominant matrices.**\n",
    "\n",
    "\n",
    "**Jacobi iterations for solving nonlinear equations:** $ f_1(x) = 0,f_2(x)=0,...,f_n(x)=0, x=(x_1,x_2,...,x_n) $\n",
    "\n",
    "$$ f_1(\\mathbf{x_1^{(k+1)}}, x_2^{(k)}, ..., x_n^{(k)})=0 $$\n",
    "$$ f_2({x_1^{(k)}}, \\mathbf{x_2^{(k+1)}}, ..., x_n^{(k)})=0 $$\n",
    "$$\\vdots$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.05053332648312581, -0.17956771997209509, -0.2529689571345711]\n",
      "[0.21180423113479036, 0.15437869579166155, 0.10962149865236209]\n",
      "[0.03623781124595906, -0.03366306109688892, -0.09040311496193816]\n",
      "[0.13385354670795924, 0.09316449631918618, 0.05345472684491567]\n",
      "[0.06533068268268807, 0.009590496806577939, -0.03591004197416863]\n",
      "[0.10885491323093499, 0.06360916807428396, 0.024524276145713682]\n",
      "[0.07992088354142911, 0.027828979346798524, -0.01425567037951667]\n",
      "[0.0987109052444953, 0.05105669920742557, 0.011469577917663582]\n",
      "[0.08634778678321582, 0.035755927434402486, -0.005233154230045877]\n",
      "[0.09441876078213328, 0.0457338518831608, 0.0057641091074878605]\n",
      "[0.08912399689112149, 0.039182686815102646, -0.0014106389758010238]\n",
      "[0.09258665432971978, 0.04346452308602135, 0.0032984488021751777]\n",
      "[0.09031756074214319, 0.04065751265787132, 0.0002198544738837235]\n",
      "[0.09180254112626063, 0.042494042012541665, 0.0022377354495803736]\n",
      "[0.09082987096261756, 0.04129089485316502, 0.0009173618955706435]\n",
      "[0.09146661246069582, 0.04207842298507372, 0.0017823023387556246]\n",
      "[0.09104962470135858, 0.04156264907497887, 0.0012161225733968131]\n",
      "[0.09132263341298519, 0.041900317716945655, 0.0015869176047715286]\n",
      "[0.09114386117517942, 0.04167919738129534, 0.0013441591775546578]\n",
      "[0.09126091277047454, 0.04182397334893029, 0.0015031259328095459]\n",
      "[0.09118426755037108, 0.04172917288808297, 0.0013990433586966253]\n",
      "[0.09123445241477442, 0.04179124464308906, 0.0014671970548809875]\n",
      "[0.0912015919549179, 0.04175060053410776, 0.0014225724185192973]\n",
      "[0.09122310816762266, 0.04177721317794751, 0.0014517922127743493]\n",
      "[0.09120901970057821, 0.041759787603524216, 0.0014326598633190154]\n"
     ]
    },
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: x not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: x not defined",
      "",
      "Stacktrace:",
      " [1] top-level scope at In[2]:29"
     ]
    }
   ],
   "source": [
    "#Example: Jacobi iterations\n",
    "\n",
    "A=[10. 2. 3.; 4. 15. 6.; 7. 8. 19.]\n",
    "\n",
    "#What happens if the diagonal elements are smaller, but A is still invertible?\n",
    "#A=[1. 2. 3.; 4. 5. 6.; 7. 8. 19.]\n",
    "\n",
    "\n",
    "M=zeros(3,3)\n",
    "\n",
    "M=diagm( 0=> diag(A))\n",
    "\n",
    "N=M-A\n",
    "\n",
    "Minv = inv(M)\n",
    "\n",
    "b=ones(3)\n",
    "bb=ones(3)\n",
    "xx=rand(3)\n",
    "\n",
    "#the following for-loop is only for illustrating teh concept of Jacobi iteration\n",
    "#in practice, one should never allocate memory for M or Minv\n",
    "for k=1:25\n",
    "   \n",
    "   global xx,bb\n",
    "   xx = Minv*(N*xx+bb) \n",
    "   println(xx)\n",
    "end\n",
    "\n",
    "println(\"\\nChecking the residual: \", A*x-b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The eigenvalues of Minv*N is important. \n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       " -0.654780355409762  \n",
       "  0.22381548783575786\n",
       "  0.4309648675740044 "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "println(\"The eigenvalues of Minv*N is important. \")\n",
    "eigvals(Minv*N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "JacobiIterations! (generic function with 1 method)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#a more appropirate way to implementing the JAcobi iterations\n",
    "function JacobiIterations!(x, A, b, numitrs)\n",
    "   \n",
    "   xnext=copy(x)\n",
    "   sum=0.0\n",
    "    \n",
    "   for itr=1:numitrs \n",
    "   for i=1:size(A,1)\n",
    "        sum=0.0\n",
    "        for j in [1:i-1; i+1:size(A,2)]\n",
    "            sum+=A[i,j]*x[j]\n",
    "        end\n",
    "        xnext[i] = (b[i]-sum)/A[i,i]\n",
    "   end\n",
    "    \n",
    "   x.=xnext\n",
    "   #println(itr,\": \", x)\n",
    "   end \n",
    "end\n",
    "\n",
    "\n",
    "function JacobiIterations!(x, A, b, tol; max_itrs=100000)\n",
    "   \n",
    "   xnext=copy(x)\n",
    "   res=1.0\n",
    "\n",
    "   itr=0\n",
    "    \n",
    "   while res>tol && itr<max_itrs\n",
    "        x.=xnext        \n",
    "    for i=1:size(A,1)\n",
    "        sum=0.0\n",
    "        for j in [1:i-1; i+1:size(A,2)]\n",
    "            sum+=A[i,j]*x[j]\n",
    "        end\n",
    "        xnext[i] = (b[i]-sum)/A[i,i]\n",
    "    end\n",
    "    \n",
    "    res=norm(xnext-x)\n",
    "    itr+=1\n",
    "   end \n",
    "    \n",
    "   println(\"Jacobi finishes in \",itr,\" iterations.\")\n",
    "\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Jacobi finishes in 0 iterations.\n",
      "[0.323553813262617, 0.289911152918378, 0.5783193054253175]\n",
      "[4.550318354738879, 8.112798379378042, 14.572232719266376]\n"
     ]
    }
   ],
   "source": [
    "x=rand(3)\n",
    "\n",
    "JacobiIterations!(x, A, b, 25)\n",
    "\n",
    "println(x)\n",
    "println(A*x-b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gauss-Seidel iterations\n",
    "\n",
    "In the Gauss-Seidel iterations, $M$ is the lower triangular part of $A$ (including the diagonal)\n",
    "\n",
    "$$ A= D-L-U, M=D-L, N=U \\implies G_{GS}:=(D-L)^{-1}U$$\n",
    "\n",
    "$$x^{(k+1)}=G_{GS} x^{(k)}+ c_{GS}$$\n",
    "\n",
    "* This algorithm is of **sequential nature**:\n",
    "    * Solve for first variable, $x_1$ in the first equation, using the previous guess for all other variables    \n",
    "$$ \n",
    "a_{11} \\mathbf{x_1} + a_{12} x_2^{(k)}+a_{13} x_3^{(k)}+...+a_{1n} x_n^{(k)}=b_1 \n",
    "$$\n",
    "    Obtain the value of $\\mathbf{x_1^{(k+1)}}$.\n",
    "    \n",
    "    * Solve for the second variable, $x_2$, in the second equation, using the value $\\mathbf{x_1^{(k+1)}}$, and the previous guess for all other variables\n",
    "$$ \n",
    "a_{21} \\mathbf{x_1^{(k+1)}} + a_{22} \\mathbf{x_2}+a_{13} x_3^{(k)}+...+a_{1n} x_n^{(k)}=b_2 \n",
    "$$\n",
    "    Obtain the value of $\\mathbf{x_2^{(k+1)}}$.\n",
    "$$ \\vdots$$\n",
    "\n",
    "Solving nonlinear equations: $ f_1(x) = 0,f_2(x)=0,...,f_n(x)=0, x=(x_1,x_2,...,x_n) $\n",
    "\n",
    "$$ f_1(\\mathbf{x_1^{(k+1)}}, x_2^{(k)}, ..., x_n^{(k)})=0 $$\n",
    "$$ f_2(\\mathbf{x_1^{(k+1)}}, \\mathbf{x_2^{(k+1)}}, ..., x_n^{(k)})=0 $$\n",
    "$$\\vdots$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GaussSeidelIterations! (generic function with 1 method)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function GaussSeidelIterations!(x, A, b, tol; max_itrs=100000)\n",
    "   \n",
    "   itr=0\n",
    "   res=1.0\n",
    "    \n",
    "   while res>tol && itr<max_itrs #for itr=1:numitrs \n",
    "    xold=copy(x)\n",
    "        \n",
    "    for i=1:size(A,1)\n",
    "        sum=0.0\n",
    "            for j in [1:i-1; i+1:size(A,2)]\n",
    "            sum+=A[i,j]*x[j]\n",
    "        end\n",
    "        x[i] = (b[i]-sum)/A[i,i]\n",
    "    end\n",
    "        \n",
    "    res=norm(xold-x)\n",
    "    itr+=1\n",
    "   end\n",
    "   \n",
    "   println(\"Gauss-Seidel finishes in \",itr,\" iterations.\")\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gauss-Seidel is faster than Jacobi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "100×100 Array{Float64,2}:\n",
       " 1.00157      0.000250114  0.00149694   …  0.000152362  0.00096891 \n",
       " 0.00112331   1.00066      0.0012599       0.00196979   0.00120999 \n",
       " 0.000949708  0.00144849   1.00193         5.50479e-5   0.00126919 \n",
       " 0.00174919   0.000610772  0.00188507      0.000116819  0.00174736 \n",
       " 0.00152282   0.000825299  0.00126561      2.8276e-5    0.000647735\n",
       " 0.00160587   0.000913379  0.00177287   …  0.000129018  0.00138372 \n",
       " 0.00163859   0.000439454  0.00134131      0.00194249   0.0016576  \n",
       " 0.00161961   0.000643386  0.00181091      0.000847802  0.000494094\n",
       " 3.86519e-6   0.00194103   0.00183352      0.00089188   0.00116436 \n",
       " 0.000924208  0.0012276    0.000769738     0.000449646  0.00095254 \n",
       " 0.000112491  0.00119056   9.38858e-7   …  0.00105343   0.000418245\n",
       " 0.00049524   0.00181685   0.001806        0.00169527   0.0014174  \n",
       " 0.00116172   4.9957e-5    0.000502303     0.000196907  0.00138762 \n",
       " ⋮                                      ⋱                          \n",
       " 0.00168877   0.000459332  0.00146636      0.000277095  0.000335432\n",
       " 0.000315832  0.00156697   0.000515194     0.00126211   0.000382564\n",
       " 0.00136167   0.00122778   0.000679116  …  0.000682485  0.000185331\n",
       " 0.000144466  0.000493564  0.000782751     0.000706767  0.000434623\n",
       " 0.00139753   0.00173488   0.0019791       0.000472515  0.000488367\n",
       " 0.0013592    0.000438725  0.00112949      0.000358258  0.00176541 \n",
       " 0.000157781  0.00196137   0.000891789     0.00169076   1.57269e-5 \n",
       " 0.000859979  0.000128232  0.000673996  …  0.000540568  0.00158413 \n",
       " 0.00056923   0.00103615   0.0016076       0.00187581   0.00181114 \n",
       " 0.00118505   0.00138502   0.00120809      0.000487718  0.00179536 \n",
       " 0.000675486  0.000603159  0.0010186       1.00156      0.0013642  \n",
       " 0.00181384   0.000540957  0.00193187      0.00154747   1.00118    "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gauss-Seidel finishes in 6 iterations.\n",
      "[0.9120309529511862, 0.9134136435738935, 0.912948635440938, 0.9060039177705486, 0.9026038099255504, 0.9110451184619583, 0.9074284029174863, 0.918765349590702, 0.9105146950900884, 0.9065923163744541, 0.913001234240074, 0.9047249312221112, 0.9147624665697199, 0.9111162661508538, 0.9176317230479875, 0.903364453387345, 0.9002610627036045, 0.9030185409403939, 0.9136841881246258, 0.9137164495444415, 0.9098211951329184, 0.9162138112671334, 0.9149218177176396, 0.9093738389861514, 0.9123114313041721, 0.909444685494364, 0.9145052067293258, 0.9133833387846739, 0.9030115836091868, 0.9140474115365999, 0.9086369072804279, 0.8978479852950875, 0.8977784629087492, 0.9096124025316028, 0.9164716800126204, 0.9027854538927532, 0.9080350853705987, 0.9083876019313846, 0.9097797497958655, 0.9174994252940214, 0.9206257643015576, 0.9035102397103149, 0.9145669444492912, 0.9022490434999164, 0.9095006574870261, 0.9006906638241964, 0.9081087668519177, 0.903050671083964, 0.9198334017704316, 0.9153917072439197, 0.9094357648688994, 0.9156196128094234, 0.9096584626308989, 0.9141302224787026, 0.9075916196045632, 0.9082501899132176, 0.9032961063941768, 0.9098181806351797, 0.9116582609214445, 0.9144526118072659, 0.9099005773937218, 0.9063616621657441, 0.9245321030034884, 0.9106920351492253, 0.9148254727469501, 0.9180704207349027, 0.9166669324716543, 0.9179061116898678, 0.9082264324903524, 0.9113155527228088, 0.910714489302068, 0.9029004602077825, 0.9032761565843909, 0.9120088857363203, 0.9075061941190623, 0.9059699914081191, 0.9055204680942488, 0.9185040506509351, 0.9083719307717008, 0.9151056912311001, 0.9065919141631896, 0.9118301906085289, 0.9209549021072783, 0.910641668695975, 0.9074346689665992, 0.910168255329857, 0.9092514853132451, 0.9082381830287243, 0.9165391860655773, 0.9097258563123216, 0.9074361402164047, 0.9172849742813564, 0.9075502638295138, 0.9147864455787139, 0.9222939411852139, 0.9121159357841137, 0.9012397390579431, 0.9082822429142173, 0.9115058013784992, 0.9102344741191543]\n",
      "[7.075806607304003e-11, 6.931255569497807e-11, 7.492539921827301e-11, 9.288947389052282e-11, 9.706457859692819e-11, 1.1197465177303911e-10, 1.0333467415080122e-10, 9.092149255707227e-11, 8.647038640674509e-11, 1.1713252590084267e-10, 1.0310929887680231e-10, 1.2066014853928664e-10, 1.1037459834994934e-10, 1.1665246546499475e-10, 1.0854450671615723e-10, 1.3989009950421405e-10, 1.2679635119638988e-10, 1.1928014131967757e-10, 1.1424572399221233e-10, 1.2190204401463234e-10, 1.037248065216545e-10, 1.0986700438309072e-10, 1.1624967655166074e-10, 1.1073719718979191e-10, 1.062130383644444e-10, 1.0421397078630434e-10, 8.633516124234575e-11, 9.733214234586285e-11, 1.012125938615327e-10, 9.14459619139052e-11, 9.544542933781486e-11, 1.0213829781946515e-10, 9.033640502309481e-11, 8.655542949043138e-11, 6.706790678379093e-11, 8.285638841698528e-11, 7.689204828409402e-11, 7.966716175644706e-11, 6.140465913517801e-11, 6.615152869926533e-11, 5.531619606813365e-11, 6.307909750091767e-11, 5.6036730811115376e-11, 6.536060581652237e-11, 4.361044858569585e-11, 5.543099312887989e-11, 4.537104025814642e-11, 4.287348254194967e-11, 3.3726355042063005e-11, 3.20958815080985e-11, 2.794697806507429e-11, 3.154898564616815e-11, 3.4796610037801656e-11, 3.090017131057721e-11, 2.645927921207658e-11, 2.4158008926633556e-11, 2.8917090943991752e-11, 1.689093309664713e-11, 2.0632162645028984e-11, 1.766231605415669e-11, 1.5758283566924547e-11, 1.6981971384666394e-11, 1.1391998455678731e-11, 1.1153744594594173e-11, 1.191757803553628e-11, 7.588818462522795e-12, 6.7912342416320826e-12, 5.375255796025158e-12, 6.280975739514361e-12, 5.999645225074346e-12, 6.902478588699523e-12, 6.5571992280410996e-12, 4.690470234436361e-12, 3.610223231476084e-12, 4.247491247610924e-12, 2.830846668189224e-12, 1.7152945730458669e-12, 1.7625900738948985e-12, 2.3450130726132556e-12, 1.4013235016818726e-12, 1.2867484855405564e-12, 7.345235530920036e-13, 7.72937269744034e-13, 4.973799150320701e-13, 1.9251267247000214e-13, 4.221067939624845e-13, 1.7497114868092467e-13, 7.416289804496046e-14, 1.2323475573339238e-13, -9.026113190202523e-14, -6.361577931102147e-14, -1.0635936575909e-13, -1.0591527654923993e-13, -2.162714451969805e-13, -1.2367884494324244e-13, -1.9828583219805296e-13, -3.0164759579065503e-13, -1.5054624213917123e-13, -9.037215420448774e-14, 0.0]\n",
      "\n",
      "Jacobi finishes in 9 iterations.\n",
      "[0.9120309495445641, 0.9134136402153573, 0.9129486320869242, 0.906003914093999, 0.9026038061390056, 0.9110451149403699, 0.9074283992503518, 0.9187653463798391, 0.9105146916309084, 0.9065923126675279, 0.9130012308282196, 0.9047249274307608, 0.914762463251356, 0.9111162626472583, 0.9176317197777655, 0.9033644495907563, 0.900261058788583, 0.9030185371547628, 0.9136841847253557, 0.9137164461207431, 0.9098211916210694, 0.9162138079993898, 0.9149218143404056, 0.909373835450122, 0.9123114278792728, 0.9094446819481674, 0.9145052034153334, 0.9133833353747078, 0.9030115798363003, 0.9140474081952786, 0.9086369036854717, 0.897847981285272, 0.8977784588982276, 0.9096123990271113, 0.916471676808686, 0.902785450114791, 0.9080350818222046, 0.9083875983386684, 0.9097797463137001, 0.9174994220841801, 0.9206257612301828, 0.9035102359291952, 0.9145669411058042, 0.9022490397597931, 0.9095006539598878, 0.9006906599590283, 0.9081087633453573, 0.9030506673085731, 0.9198333986391798, 0.9153917039888196, 0.9094357614014991, 0.91561960956299, 0.909658459138321, 0.9141302191367606, 0.9075916160765096, 0.9082501864004598, 0.9032961027116578, 0.9098181771619462, 0.9116582575528255, 0.9144526085230724, 0.9099005739592368, 0.9063616585195027, 0.9245321001293133, 0.9106920317653293, 0.9148254695315124, 0.9180704175866601, 0.9166669292516898, 0.9179061085414724, 0.9082264289860028, 0.9113155493274695, 0.9107144858782037, 0.9029004565229446, 0.9032761528509005, 0.9120088823941436, 0.9075061906267361, 0.905969987850436, 0.9055204644871563, 0.9185040475687953, 0.9083719272721061, 0.915105688007237, 0.906591910573825, 0.9118301872843476, 0.9209548990790987, 0.9106416653106479, 0.9074346654734091, 0.9101682519263338, 0.909251481838302, 0.9082381795709833, 0.9165391828513386, 0.9097258529024068, 0.907436136709689, 0.9172849711330436, 0.9075502603081099, 0.9147864423071675, 0.9222939382551639, 0.9121159323971616, 0.9012397352883125, 0.9082822394257438, 0.9115057980291166, 0.910234470720035]\n",
      "[-3.670757475759956e-9, -3.6154575999702843e-9, -3.609662679870951e-9, -3.940836212912302e-9, -4.058058333811232e-9, -3.749733523683574e-9, -3.915619273264781e-9, -3.428328065879782e-9, -3.712229190711014e-9, -3.944096715891021e-9, -3.6391273328106877e-9, -4.033442246864638e-9, -3.531618220087296e-9, -3.7247508410942487e-9, -3.474564080896414e-9, -4.024646504952045e-9, -4.168467016008037e-9, -4.03422228956174e-9, -3.614801014073521e-9, -3.628847666803381e-9, -3.750675769964573e-9, -3.472608200993932e-9, -3.5833822575881413e-9, -3.7696673560461136e-9, -3.6517642243438786e-9, -3.787348767936294e-9, -3.5523444186225106e-9, -3.642447676810434e-9, -4.039491963148123e-9, -3.575722939963555e-9, -3.847961171921099e-9, -4.294005262472922e-9, -4.308401635455539e-9, -3.761647882072339e-9, -3.4519302971602883e-9, -4.062707503749152e-9, -3.820245786378962e-9, -3.861320818643321e-9, -3.7641032513136e-9, -3.457979458332261e-9, -3.3178365610453397e-9, -4.0853750382652265e-9, -3.612433463473508e-9, -4.044430346183958e-9, -3.827826833280312e-9, -4.189719682301529e-9, -3.8079772668453415e-9, -4.102124528948536e-9, -3.402572223976108e-9, -3.545323368214781e-9, -3.782654744988179e-9, -3.5341443105352255e-9, -3.800630254957582e-9, -3.6376787138081568e-9, -3.849544238931912e-9, -3.8346769093422495e-9, -4.021177058000092e-9, -3.799842107632401e-9, -3.6810450243507375e-9, -3.589760488864613e-9, -3.761001288182797e-9, -3.988106178631767e-9, -3.147854865659383e-9, -3.7110361450487517e-9, -3.5242937457269363e-9, -3.4533323978180874e-9, -3.531226311359603e-9, -3.4548287564106772e-9, -3.847162255432579e-9, -3.72558406347423e-9, -3.756728816917132e-9, -4.047894575087696e-9, -4.096865180436282e-9, -3.672907089580235e-9, -3.838160678171221e-9, -3.91034249425104e-9, -3.968554596056606e-9, -3.3876513816366582e-9, -3.846544527341678e-9, -3.546183902081168e-9, -3.942255633049285e-9, -3.658617742097192e-9, -3.328717634865086e-9, -3.7242650074986727e-9, -3.844866425239957e-9, -3.74590114482487e-9, -3.817943294848192e-9, -3.8053310502661475e-9, -3.5312993640346235e-9, -3.75355069248684e-9, -3.858145802837498e-9, -3.4632134937595538e-9, -3.871414633316306e-9, -3.5955702859524763e-9, -3.2250915271703207e-9, -3.721958741209619e-9, -4.147586274427795e-9, -3.8363411336561626e-9, -3.6857147334146134e-9, -3.738665821373388e-9]\n"
     ]
    }
   ],
   "source": [
    "A=[10. 2. 3.; -4. 5. 1.3; 2. 8. 19.]\n",
    "\n",
    "N=100\n",
    "\n",
    "A = I(N)+ 0.002*rand(N,N)\n",
    "\n",
    "display(A)\n",
    "\n",
    "b=ones(N)\n",
    "\n",
    "x0=rand(N)\n",
    "\n",
    "x=copy(x0)\n",
    "GaussSeidelIterations!(x, A, b, 1e-7)\n",
    "\n",
    "println(x)\n",
    "println(A*x-b,\"\\n\")\n",
    "\n",
    "x=copy(x0)\n",
    "JacobiIterations!(x, A, b, 1e-7)\n",
    "\n",
    "println(x)\n",
    "println(A*x-b)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Successive over relaxation (SOR)\n",
    "\n",
    "SOR is obtained by *averaging* the Gauss-Seidel update, $x^{(k+1)}_{GS}$, with the approximation computed in the previous iteration, $x^{(k)}$:\n",
    "\n",
    "$$ x^{(k+1)}_{SOR}:= \\omega x^{(k+1)}_{GS}+(1-\\omega) x^{(k)}. $$\n",
    "\n",
    "\n",
    "Three important theorems about SOR:\n",
    "\n",
    "Theorem (Kahan): Let $a_{ii}\\neq 0, i=1,2,\\cdots,n$. **If SOR converges, then $0<\\omega<2$.**\n",
    "\n",
    "Theorem (Ostrowski-Reich): **SOR converges for all $0<\\omega<2$ iff $A$ is invertible, symmetric, positive definite.**\n",
    "\n",
    "**Choose $\\omega$ to make $\\rho(G_{SOR})$ as small as possible.**\n",
    "\n",
    "* This was pioneered by David Young, here at UT Ausitn in 1950'\n",
    "\n",
    "Theorem (Young): $\\omega_{opt}$ as a function of $\\rho(G_J)$.\n",
    "\n",
    "For more in-depth exposition about the choice of $\\omega$, check out GGK, Sec 11.3.4.\n",
    "\n",
    "## Symmetric Successive over relaxation (SSOR)\n",
    "\n",
    "Involves two SOR updates:\n",
    "\n",
    "* One SOR step in which the G-S update involves solving each row of the linear system in **increasing order** (row 1, than row 2,..., row n). \n",
    "    * In other words, the G-S iteration is defined by $A=D-L-U=M-N$, with $M=D-L$.\n",
    "\n",
    "* One SOR step in which the G-S update involves solving each row of the linear system in **decreasing order** (row n, than row (n-1),..., row 1).\n",
    "    * In other words, the G-S iteration is defined by $A=D-L-U=M-N$, with $M=D-U$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SOR! (generic function with 1 method)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function SOR!(x, A, b, ω, tol; max_itrs=100000)\n",
    "   \n",
    "   sum=0.0\n",
    "   itr=0\n",
    "   xold=copy(x)\n",
    "    \n",
    "   res=1.0 \n",
    "    \n",
    "   while res>tol && itr<max_itrs \n",
    "    xold.=x\n",
    "    for i=1:size(A,1)\n",
    "        sum=0.0\n",
    "        for j=1:i-1 \n",
    "            sum+=A[i,j]*x[j]\n",
    "        end\n",
    "        for j=i+1:size(A,2) \n",
    "            sum+=A[i,j]*x[j]\n",
    "        end\n",
    "\n",
    "        x[i] = ω*(b[i]-sum)/A[i,i]+(1.0-ω)*xold[i]\n",
    "    end\n",
    "    res=norm(xold-x)\n",
    "    itr+=1\n",
    "   end\n",
    "   \n",
    "   println(\"SOR finishes in \",itr,\" iterations.\")\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gauss-Seidel finishes in 25 iterations.\n",
      "[0.17124735729379462, 0.052854122621400344, -0.04439746300198322]\n",
      "[-3.8125058665627876e-13, -1.3171685964152857e-12, 0.0]\n",
      "\n",
      "Jacobi finishes in 481 iterations.\n",
      "[0.17124735729413046, 0.052854122621843184, -0.04439746300180666]\n",
      "[3.04911651483053e-12, 5.956790616323815e-12, 8.365530490550555e-12]\n",
      "\n",
      "SOR finishes in 59 iterations.\n",
      "  0.000050 seconds (137 allocations: 8.188 KiB)\n",
      "[0.17124735729340707, 0.0528541226203849, -0.04439746300088907]\n",
      "[-1.4548362514688051e-12, -7.472356067239616e-12, 4.481748305806832e-12]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#Example: Jacobi iterations\n",
    "\n",
    "A=[6. 2. 3.; 4. 11. 6.; 7. 8. 14.]\n",
    "\n",
    "\n",
    "M=zeros(3,3)\n",
    "\n",
    "M=diagm( 0=> diag(A,0))\n",
    "\n",
    "N=M-A\n",
    "\n",
    "Minv = inv(M)\n",
    "\n",
    "tol=1e-12\n",
    "\n",
    "b=ones(3)\n",
    "x=rand(3)\n",
    "\n",
    "x0=rand(3)\n",
    "\n",
    "x=copy(x0)\n",
    "GaussSeidelIterations!(x, A, b, tol)\n",
    "\n",
    "println(x)\n",
    "println(A*x-b,\"\\n\")\n",
    "\n",
    "x=copy(x0)\n",
    "JacobiIterations!(x, A, b, tol)\n",
    "\n",
    "println(x)\n",
    "println(A*x-b,\"\\n\")\n",
    "\n",
    "#SOR with non-optimal weight\n",
    "x=copy(x0)\n",
    "@time SOR!(x, A, b, 0.68, tol)\n",
    "\n",
    "println(x)\n",
    "println(A*x-b,\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: solving a two-point boundary value problem for Poisson's equation\n",
    "\n",
    "$$ -u_{xx} =f(x),~~~0<x<1,~~~u(0)=u(1)=0 $$\n",
    "\n",
    "* Discretize the unit interval: $x_j = j\\Delta x, j=0,1,2,\\cdots,N,\\Delta x =1/N$\n",
    "\n",
    "* Approximate the differential equation on the grid nodes:\n",
    "    $$ -u_{xx} (x_j)=f(x_j) \\approx -(u_{j+1}-2u_j+u_{j-1})/\\Delta x^2 = f(x_j),~~~j=1,2,\\cdots,N-1 $$\n",
    "\n",
    "* Set the boundary conditions: $u_0=0=u_N$.\n",
    "\n",
    "This leads to a linear system to solve: $AU=F$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.000204 seconds (11 allocations: 38.938 KiB)\n",
      "\n",
      "\n",
      "Gauss-Seidel finishes in 4016 iterations.\n",
      "  0.033931 seconds (204.90 k allocations: 96.885 MiB, 18.40% gc time)\n",
      "\n",
      "\n",
      "Jacobi finishes in 8021 iterations.\n",
      "  0.050704 seconds (401.13 k allocations: 189.708 MiB, 9.77% gc time)\n"
     ]
    }
   ],
   "source": [
    "using LinearAlgebra\n",
    "\n",
    "N=50\n",
    "xgrid=range(0.0, stop=1.0, length=N+1) #linspace(0,1,N+1)\n",
    "Δx=1/N\n",
    "\n",
    "f=x-> x.^2+0.0*sin(20*x)\n",
    "\n",
    "#The vector U stores (u1,u2,u3,cdots, u_(N-1))\n",
    "U=zeros(N-1)\n",
    "\n",
    "#N-2 because we exclude x=0 and X=1.\n",
    "#A=spdiagm( (-ones(N-2), 2.0*ones(N-1),-ones(N-2)), [-1,0,1])\n",
    "A=-2.0 .*diagm(0=>ones(N-1)) .+ diagm(1=>ones(N-2)) .+ diagm(-1=>ones(N-2))\n",
    "F=collect(f.(xgrid[2:end-1]))\n",
    "\n",
    "b=F.*Δx^2\n",
    "\n",
    "@time U=-A\\b\n",
    "\n",
    "Ugs=copy(b)\n",
    "Uj=copy(U)\n",
    "Usor=copy(U)\n",
    "\n",
    "println(\"\\n\")\n",
    "@time GaussSeidelIterations!(Ugs, A, b, 1e-10)\n",
    "\n",
    "println(\"\\n\")\n",
    "@time JacobiIterations!(Uj, A, b, 1e-10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The timings above for Gauss-Seidel and Jacobi iterations are not fair, because in real implementation, one would skip a lof of multiplications involving the zero entries in the matrix $A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAGdCAYAAADqsoKGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzde1iUZeI+8HtmgEE5DCcBQRTEEwcRQUVwPVWLWKZWCmU77rarm7WVZCdNLfVXkdWW25qWm1ttmbKlppknaK1QAU+ARzwLKGfEGQ4yM8y8vz/Q6UugMgi8c7g/1zXXtQ3PzNy0Brfv8z7PIxEEQQARERGRhZOKHYCIiIioI7DUEBERkVVgqSEiIiKrwFJDREREVoGlhoiIiKwCSw0RERFZBZYaIiIisgosNURERGQV7MQO0FUMBgOKi4vh4uICiUQidhwiIiJqA0EQUFNTAz8/P0ilt78WYzOlpri4GAEBAWLHICIionYoKipCr169bjvGZkqNi4sLgKZ/Ka6uriKnISIiorZQq9UICAgw/h6/HZspNTennFxdXVlqiIiILExbbh3hjcJERERkFVhqiIiIyCqw1BAREZFVaFepWbVqFYKCguDo6Ijo6GhkZGTcdvzGjRsRGhoKuVyO0NBQbN68+ZZjn3zySUgkEqxYsaLZ89XV1VAqlVAoFFAoFFAqlbh27Vp74hMREZEVMrnUpKamIjk5GQsXLkROTg5Gjx6NiRMnorCwsNXxmZmZSEpKglKpRF5eHpRKJRITE5Gdnd1i7HfffYfs7Gz4+fm1+NqMGTOQm5uLnTt3YufOncjNzYVSqTQ1PhEREVkpiSAIgikviImJQVRUFFavXm18LiQkBFOnTkVKSkqL8UlJSVCr1dixY4fxuYSEBLi7u2P9+vXG565cuYKYmBjs2rULDzzwAJKTk5GcnAwAOHXqFEJDQ5GVlYWYmBgAQFZWFmJjY5Gfn4+BAwfeMbdarYZCoYBKpeLqJyIiIgthyu9vk67UaLVaHD58GPHx8c2ej4+Px/79+1t9TWZmZovxEyZMaDbeYDBAqVTipZdeQlhYWKvvoVAojIUGAEaOHAmFQnHLzyUiIiLbYtI+NZWVldDr9fDx8Wn2vI+PD0pLS1t9TWlp6R3HL1++HHZ2dnjuuedu+R7e3t4tnvf29r7l52o0Gmg0GuM/q9Xq1r8pIiIisgrtulH4txvgCIJw201xbjf+8OHD+Mc//oHPP//cpPe40+empKQYbypWKBQ8IoGIiMjKmVRqvLy8IJPJWlwdKS8vb3E15iZfX9/bjs/IyEB5eTl69+4NOzs72NnZoaCgAC+88AICAwON71FWVtbivSsqKm75uQsWLIBKpTI+ioqKTPlWiYiIyMKYVGocHBwQHR2NtLS0Zs+npaUhLi6u1dfExsa2GL97927jeKVSiaNHjyI3N9f48PPzw0svvYRdu3YZ30OlUuHAgQPG98jOzoZKpbrl58rlcuORCDwagYiIyPqZfPbTvHnzoFQqMWzYMMTGxmLNmjUoLCzEnDlzAAAzZ86Ev7+/cSXU3LlzMWbMGCxfvhxTpkzBli1bkJ6ejr179wIAPD094enp2ewz7O3t4evra1zVFBISgoSEBMyePRuffPIJAOCvf/0rJk2a1KaVT0RERGT9TC41SUlJqKqqwrJly1BSUoLw8HBs374dffr0AQAUFhZCKv31AlBcXBw2bNiARYsWYfHixQgODkZqamqzlUxtsW7dOjz33HPGlVSTJ0/GypUrTY1PRGQRDAYBpeoGXKqsw8WqOhRW1UOrN9zxdS5yO/Tt4YzgHs7o28MJTnKbObeYyPR9aiwV96khInNUr23E8StqXKysxcXKelyqrMOlqqZHg+7OJeZOeioc0beHE4JvFJ3gHs7o5+0MX4VjB6Qn6nym/P5mhSci6kKCIOBceS1+Ol2Bn89U4MDFq7e8AmMnlSDAozsCPbujj6cTnOSyO77/1TotzlfU4UJFLSprtShRNaBE1YB956qajRvo44IJ4b5ICPNFSE+X264+JbIUvFJDRNTJahp02HeuCj+fKcfPpytQrGpo9nVfV0f093FGkJcTAj2dENTDCUGeTvB37wZ7WfvPHVbV63C+shbny2txvqIO5ytqcb6iFgVV9dAbfv3R38ezOxLCfDEh3BeRvdwglbLgkPkw5fc3Sw0RUSeoqtVg05ErSDtVhiMF1Wj8PyXCwU6KkX09MW5AD4wd2AN9vZy69EqJql6HH/PLsON4KX45UwFN469XinxdHTEhzAcTwn0xItADdndRqog6AktNK1hqiKizCYKAQwXVWJdVgO3HSptNK/X1csLYgT0wdkAPxAR5opvDnaeSukKdphE/n6nAjuOl+N+pMtRp9cav+bt1w9/G98O06F5wsGO5IXGw1LSCpYaIOktNgw7f5VzBV1mFOF1WY3w+opcC06J7YdwAb/T27C5iwrZp0Omx/3wldh4vxe6TZbhWrwPQVG6eGheM6cN6QW5nHmWMbAdLTStYaoioo50oVmFddiG+y7mC+htXOBztpZg8xA9/GNkHEb3cRE7Yfg06Pb7OLsTHP59HeU3TOXp+Ckc8Nb4fElluqAux1LSCpYaIOoIgCNhzuhwr/3cORwqvGZ8P7uGEx2P64JGoXlB0txcxYcdq0Omx/kAhVv/0a7npqXDEU+OCkTgsAI72LDfUuVhqWsFSQ0R368DFq3hnZz4OFVQDaFpyPSHcF3+I6YORfT2sell0g06P1INFWP3TeZSqm1Zv+bo2lZvHY3rzhmLqNCw1rWCpIaL2OlGswnu7TmPP6QoAgNxOij/FBeIvo4Pg7WJbm9g16PT476EirNrza7kZEuCG9xOHILiHs8jpyBqx1LSCpYaITHWpsg7vp53B1rxiAIBMKkHS8AA8d09/m9+RV9PYdOXm3V2nUdPQCLmdFPMnDsIfYwO5zw11KJaaVrDUEFFblakb8OGPZ5F6sMi4v8yDQ/ww7/cDEOTlJHI681J87Tpe2XgUGWcrAQCxfT3x7vQI9HI3/9VeZBlYalrBUkNEd9Kg0+OjPefwr4wLxnOXxg3sgRfjByLcXyFyOvMlCAK+yirAW9vzcV2nh7PcDq9NCsX0Yb2s+j4j6hosNa1gqSGi28m6UIVXNx3Dhco6AEB0H3e8PGEgYvp6ipzMclyqrMML3+Th8I0bqe8d5I2URwbb3H1H1LFYalrBUkNErVFd1+HtHaew/kARAMDbRY6lk8OQEO7LqwztoDcI+FfGBby/+wy0egPcu9vjjamD8UBET7GjkYViqWkFSw0R/dbO4yVYvOUEKm7svzIjpjdeSRgERTfr2WdGLPmlasxLzcPJEjUA4LERAVg2JfyuDugk28RS0wqWGiK6qVTVgNe2HMfuk2UAgL49nJDy0GBONXUwbaMB//zfWXy05xwMAjB2QA989HgUnOV2YkcjC8JS0wqWGiIyGASsP1iIt7fno0bTCDupBE+NC8bfxvfjzridKP1kGZ5dn4PrOj3C/Fzx2Z+Gw9uV99lQ27DUtIKlhsi2FV2txwv/zcOBS1cBAJEBbnj7kcEY5MufB10hr+ga/vz5QVTVaeHv1g1f/Hk4+nm7iB2LLIApv785uUlEVm/n8RI88GEGDly6iu4OMrz+YCg2PhXHQtOFhgS4YdPTcQjycsKVa9fx8Kr9yL5QJXYssjIsNURktRp0ery+5TjmfHUE6oZGRAa4YVfyGDwxKggy7nrb5fp4OmHjU3GI6u0GdUMjlGsP4PsbuzUTdQSWGiKyShcr6/DI6v34IrMAAPDkmL74Zk4sAjy4062YPJwc8PXskZgQ5gOt3oBn1+dgzS/nYSN3QlAnY6khIquzJfcKJn2YgRPFang4OeCzJ4Zjwf0hXE5sJhztZVj1eDT+FBcIAHhrez6WbD0BvYHFhu4O19URkdW4rtVj2bYTxo30RgR54MNHh9r84ZPmSCaV4PUHQ9HLvRve+OEUvsgsQKm6Af98LAoOdiyf1D78k0NEVuFsWQ2mfrQP6w8UQSIBnru3P76eFcNCY8YkEglmje6LlTOGwkEmxa4TZXjxmzwYeMWG2olXaojI4n2XcwULNh3DdZ0ePVzkWJEUiVH9vMSORW00KcIPznI7zPriELbmFcPDyQGvPxjKYyrIZLxSQ0QWy2AQ8M7OfCSn5uK6To/f9fPC9udGs9BYoHEDvfH3xCEAgM/3X8I//3dO5ERkiXilhogsUr22EfNS87DzRCkA4OlxwXgxfiCkXKptsaZE+qO6Tosl35/E+2ln4O7kAOXIPmLHIgvCUkNEFqdEdR2zvjiEE8VqOMikePuRwXg4qpfYsagD/GlUEK7W6/Dhj2fx2pbjcO9uj0kRfmLHIgvBUkNEFiWv6Bpm/ecQKmo08HRywCfKaAwL9BA7FnWg5+/rj+o6Lb7MKsDzqblwdbTHmAE9xI5FFoD31BCRxdh2tBiJn2SiokaDgT4u+O5vo1horJBEIsGSyWGYFNETOr2AJ788jJzCarFjkQVgqSEisycIAlakn8EzX+dA02jAvYO8sfHpOO4ObMVkUgneT4zE6P5euK7T44nPD+JsWY3YscjMsdQQkVlr0Onx7PocrEg/CwD465i+WDNzGJzlnD23dg52Unz8h2hEBrjhWr0OyrUHcLm6XuxYZMZYaojIbF2t0yJpTRa2HS2BvUyCdx6JwKv3h/AwShviJLfDZ38ajn7ezihVN2Dm2gOoqtWIHYvMFEsNEZml4mvXMf3j/cgruga37vb48i8xSBweIHYsEoG7kwO+/MsI+Lt1w4XKOjzzdQ7PiaJWtavUrFq1CkFBQXB0dER0dDQyMjJuO37jxo0IDQ2FXC5HaGgoNm/e3OzrS5YswaBBg+Dk5AR3d3fcd999yM7ObjYmMDAQEomk2WP+/PntiU9EZu58RS2mrd6P8xV18FM4YuNTcRjZ11PsWCSinopu+PyJ4ejuIEPmhSp8kHZG7EhkhkwuNampqUhOTsbChQuRk5OD0aNHY+LEiSgsLGx1fGZmJpKSkqBUKpGXlwelUonExMRmpWXAgAFYuXIljh07hr179yIwMBDx8fGoqKho9l7Lli1DSUmJ8bFo0SJT4xORmTt2WYXpH2eiWNWA4B5O+PapOAT3cBY7FpmB/j4uSHl4MABg5Z5z2HO6XOREZG4kgiCYdA0vJiYGUVFRWL16tfG5kJAQTJ06FSkpKS3GJyUlQa1WY8eOHcbnEhIS4O7ujvXr17f6GWq1GgqFAunp6bj33nsBNF2pSU5ORnJysilxW7ynSqWCq6tru96DiDpX5vkqzP7PIdRqGhHRS4HP/jQcns5ysWORmVn83XF8mVUAt+722Pbs79DLnavgrJkpv79NulKj1Wpx+PBhxMfHN3s+Pj4e+/fvb/U1mZmZLcZPmDDhluO1Wi3WrFkDhUKBIUOGNPva8uXL4enpicjISLz55pvQarW3zKrRaKBWq5s9iMh87TpRij9+dgC1mkbE9vXE17NHstBQqxZNCkFELwWu1evwt69zoGnUix2JzIRJpaayshJ6vR4+Pj7Nnvfx8UFpaWmrryktLW3T+G3btsHZ2RmOjo744IMPkJaWBi+vXw+lmzt3LjZs2IA9e/bgmWeewYoVK/D000/fMmtKSgoUCoXxERDAGwyJzNU3h4rw1FeHoW00ID7UB589MZxLtumW5HYyfDQjCopu9sgruoa3fjgldiQyE+26Ufi3x8ELgnDbI+LbMn78+PHIzc3F/v37kZCQgMTERJSX/zpf+vzzz2Ps2LGIiIjArFmz8PHHH2Pt2rWoqqpq9TMXLFgAlUplfBQVFZn6bRJRF/g04wJe+vYoDAIwPboXVj0eBUd7mdixyMwFeHTHB0lNV/O/yCzA1rxikROROTCp1Hh5eUEmk7W4ylJeXt7iasxNvr6+bRrv5OSEfv36YeTIkVi7di3s7Oywdu3aW2YZOXIkAODcudaPp5fL5XB1dW32ICLzIQgC3t2Vjzdu/C179uggvDMtAnYy7jRBbXPPIB88PS4YADB/41GcK68VORGJzaSfHg4ODoiOjkZaWlqz59PS0hAXF9fqa2JjY1uM37179y3H3yQIAjSaW2+wlJOTAwDo2bNnW6ITkRkRBAFLvz+Jj/acBwC8nDAQr94fctsrvkStmff7ARjZ1wP1Wj2eXncY9dpGsSORiEyetJ43bx6USiWGDRuG2NhYrFmzBoWFhZgzZw4AYObMmfD39zeuhJo7dy7GjBmD5cuXY8qUKdiyZQvS09Oxd+9eAEBdXR3efPNNTJ48GT179kRVVRVWrVqFy5cvY/r06QCabjbOysrC+PHjoVAocPDgQTz//POYPHkyevfu3VH/LoioCwiCgGXbTuLz/ZcAAG8+FI7HY/qIG4oslp1Mig8fG4oHPtyLM2W1WLj5ON5PHMKCbKNMLjVJSUmoqqoy7hkTHh6O7du3o0+fph9KhYWFkEp/vQAUFxeHDRs2YNGiRVi8eDGCg4ORmpqKmJgYAIBMJkN+fj6++OILVFZWwtPTE8OHD0dGRgbCwsIANE0lpaamYunSpdBoNOjTpw9mz56Nl19+uSP+HRBRFxEEAf9v2yl8tu8SAGD5I4ORNJx/MaG74+3iiJWPDcWMT7OxOecKhgd6YEYM/1zZIpP3qbFU3KeGSFyCIOCNH05h7d6LAICUhwfjsRH8xUMd5+Ofz+PtHflwkEmx8ak4DO6lEDsSdYBO26eGiKg9BEHAm/+n0Lz1EAsNdby/ju6L+0J8oNUb8Mz6I7iu5f41toalhog6lSAISNmRj09vFJo3Hwrn1AB1CqlUgr9PH4KeCkcUVNXj/bTTYkeiLsZSQ0SdRhAEvL0jH2t+uQAA+H9TeVMwdS5Fd3u8+VA4AGDt3ovILbomciLqSiw1RNQpBEHA2zvz8cnNQjMlDMqRLDTU+e4Z5IOpkX4wCMAr3x6FttEgdiTqIiw1RNThBEHAO7tO45OfmwrNsilhUMYGihuKbMprD4bB08kBp8tqsOqn1jdpJevDUkNEHe7vu89g9U9NG+stnRyGmSw01MU8nBywZHLTtiAf7TmH06U1IieirsBSQ0Qd6pOfz2Plnqa/GS95MBR/jAsUNxDZrEkRPXFfiA90egEvf5sHvcEmdjCxaSw1RNRhUg8WImVHPgDglYRB+NOoIJETkS2TSCR486FwuDjaIe+yCp/tuyh2JOpkLDVE1CF2Hi/Bgk3HAABPju2Lp24cNEgkJh9XRyy8PwQA8N7u07hUWSdyIupMLDVEdNf2nq3Ec+tzYRCAR4cHYH7CILEjERklDQ9AXLAnGnQGzN90FDaykb5NYqkhoruSU1iNv355CFq9ARPDffHmQ4N5mCCZFYlEgrcfjoCjvRRZF65i/YEisSNRJ2GpIaJ2O1NWgyc+P4h6rR6/6+eFFY9GQiZloSHz09uzO16MHwgASNl+CiWq6yInos7AUkNE7VJ0tR7Ktdm4Vq9DZIAbPlFGQ24nEzsW0S09MSoIkQFuqNE0YtHm45yGskIsNURksooaDZRrs1Gm1qC/tzM++9NwOMntxI5FdFsyqQTvTIuAvUyCH/PLsTWvWOxI1MFYaojIJKrrOsz89wFcqqpHL/du+PIvMXB3chA7FlGbDPBxwbP39AcALP3+JK7WaUVORB2JpYaI2qxBp8esLw7iVIkaXs5yfPWXGPgqHMWORWSSOWODMcjXBVfrtFiRfkbsONSBWGqIqE30BgHJG3Jx8FI1XBzt8J8/j0Cgl5PYsYhM5mAnxWsPhgIA1mUX4lw5j1CwFiw1RNQmb20/hZ0nSuEgk+LTmcMQ6ucqdiSidosL9sJ9IT7QGwS8tT1f7DjUQVhqiOiOPtt3EWv3Nm0x/17iEMT09RQ5EdHdW3D/INhJJfhffjn2nq0UOw51AJYaIrqtXSdKsWzbSQBN5zlNHuInciKijhHcwxl/GNkHAPDGDyd54KUVYKkholvKKazG3A05EARgRkxvzBnbV+xIRB1q7r394epoh/zSGnx7mDsNWzqWGiJqVUFVHWZ9cQgNOgPGD+yBZZPDePwBWR13Jwc8d2/TEu/3dp9BraZR5ER0N1hqiKiF6jotnvjsIKrqtAj3d8XKGVGwk/HHBVknZWwf9PHsjooaDT75+bzYcegu8KcUETXToNNj9n8O4UJlHfzduuHff+RuwWTd5HYyLJjYdLL8vzIuoPgaz4WyVCw1RGRkMAh44Zs8HCpo2ovmsyeGw9uVm+uR9ZsQ5osRgR5o0Bnw3q7TYsehdmKpISKj5Tvz8cPREtjLJPhEGY0BPi5iRyLqEhKJBAsfCAEAbMq5gqOXr4mciNqDpYaIAABfZRXgk18uAACWPxKBuGAvkRMRda0hAW54aKg/AOCNbad4ircFYqkhIuw/V4nXt54AAMz7/QA8HNVL5ERE4nhpwkDI7aQ4cOkqdp0oEzsOmYilhsjGXaqsw1PrjkBvEPDQUH88e08/sSMRicbPrRtmj27ajyllxyloGw0iJyJTsNQQ2TDVdR3+8sVBqK7rEBnghpSHB3MvGrJ5c8YFw8tZjoKqevwn85LYccgELDVENqpRb8Cz63NwvqIOPRWOWDMzGo72MrFjEYnOWW6HF+MHAAA+/PEsquu0IieitmKpIbJRKTvy8cuZCjjaS/GvmcPg7cKl20Q3TR8WgEG+LlA3NOKf/zsndhxqI5YaIhuUerDQeOr2+4mRCPdXiJyIyLzIpBIsuL9pife67AKU1zSInIjagqWGyMYcuHgVi747DgB4/r4BuH9wT5ETEZmnMf29MLS3GzSNBvzrxnYHZN7aVWpWrVqFoKAgODo6Ijo6GhkZGbcdv3HjRoSGhkIulyM0NBSbN29u9vUlS5Zg0KBBcHJygru7O+677z5kZ2c3G1NdXQ2lUgmFQgGFQgGlUolr17g5EpEpiq7WY85Xh6HTC3ggoieeu5crnYhuRSKR4Ll7mg67/CqrEFW1GpET0Z2YXGpSU1ORnJyMhQsXIicnB6NHj8bEiRNRWFjY6vjMzEwkJSVBqVQiLy8PSqUSiYmJzUrLgAEDsHLlShw7dgx79+5FYGAg4uPjUVFRYRwzY8YM5ObmYufOndi5cydyc3OhVCrb8S0T2aZaTSNmfXEIV28cUvnetCFc6UR0B+MG9sBgfwWu6/T49MaULZkviWDilokxMTGIiorC6tWrjc+FhIRg6tSpSElJaTE+KSkJarUaO3bsMD6XkJAAd3d3rF+/vtXPUKvVUCgUSE9Px7333otTp04hNDQUWVlZiImJAQBkZWUhNjYW+fn5GDhw4B1z33xPlUoFV1dXU75lIounNwh48stDSD9VDm8XObY8Mwo9Fd3EjkVkEdJOlmH2fw7ByUGGva/cA3cnB7Ej2RRTfn+bdKVGq9Xi8OHDiI+Pb/Z8fHw89u/f3+prMjMzW4yfMGHCLcdrtVqsWbMGCoUCQ4YMMb6HQqEwFhoAGDlyJBQKxS3fR6PRQK1WN3sQ2ar3dp9G+qlyONhJsWbmMBYaIhPcF+KNkJ6uqNPq8e99vFpjzkwqNZWVldDr9fDx8Wn2vI+PD0pLS1t9TWlpaZvGb9u2Dc7OznB0dMQHH3yAtLQ0eHl5Gd/D29u7xXt7e3vf8nNTUlKM998oFAoEBAS0+fsksibf5xVj9U/nAQDvTotAZICbyImILEvTvTVN9599vu8SVNd1IieiW2nXjcK/nYcXBOG2c/NtGT9+/Hjk5uZi//79SEhIQGJiIsrLy2/5Hnf63AULFkClUhkfRUVFd/y+iKxNfqkaL397FADw5Ni+mBLpL3IiIss0IcwXA3ycUaNpxOf7Lokdh27BpFLj5eUFmUzW4upIeXl5i6sxN/n6+rZpvJOTE/r164eRI0di7dq1sLOzw9q1a43vUVbW8mCxioqKW36uXC6Hq6trsweRLVHV6/Dkl4dxXafH6P5eeHnCILEjEVksqVSCZ2+shFq79wJqGni1xhyZVGocHBwQHR2NtLS0Zs+npaUhLi6u1dfExsa2GL979+5bjr9JEARoNBrje6hUKhw4cMD49ezsbKhUqju+D5EtMhgEJKfmoKCqHr3cu+HDR4dCJuVKJ6K7cf/gnujbwwnqhkb8J7NA7DjUCpOnn+bNm4dPP/0U//73v3Hq1Ck8//zzKCwsxJw5cwAAM2fOxIIFC4zj586di927d2P58uXIz8/H8uXLkZ6ejuTkZABAXV0dXn31VWRlZaGgoABHjhzBrFmzcPnyZUyfPh1A0+qqhIQEzJ49G1lZWcjKysLs2bMxadKkNq18IrI1K9LPYM/pCsjtpPj4D9FcrUHUAWRSifEU+08zLqBO0yhyIvotk0tNUlISVqxYgWXLliEyMhK//PILtm/fjj59+gAACgsLUVJSYhwfFxeHDRs24LPPPkNERAQ+//xzpKamGlcyyWQy5Ofn45FHHsGAAQMwadIkVFRUICMjA2FhYcb3WbduHQYPHoz4+HjEx8cjIiICX3755d1+/0RWZ/eJUnx446yatx8ZzCMQiDrQgxF+CPTsjup6Hb7K4tUac2PyPjWWivvUkC04V16LqR/tQ62mEX+KC8SSyWF3fhERmeS/h4rw8rdH4eXsgIyX70E3B55u35k6bZ8aIjJfNQ06PPnlIdRqGjEiyAMLHwgROxKRVXpoqD96uXdDZa0WXx9ofTd9EgdLDZEVMBgEvPhNHs5X1MHX1REfzYiCvYz/eRN1BnuZFH8b33RvzSc/n0eDTi9yIrqJP/WIrMDqn89j14kyOMikWP2HKPRwkYsdiciqPRLVC34KR5TXaPDfQ9wHzVyw1BBZuJ9Ol+O93acBAEunhGFob3eRExFZPwc7KZ4aFwwAWP3TeWgaebXGHLDUEFmwwqp6zN2QC0EAHhsRgMdG9BY7EpHNmD4sAD6ucpSoGrDx8BWx4xBYaogsVoNOj6fWHYbqug6RAW5c6UTUxRztZZgztulqzUd7zkGnN4iciFhqiCzUsm0ncaJYDffu9lj1eBTkdlxWStTVHhvRG17ODrhy7Tp2nWj9gGXqOiw1RBZoc85lfJ1dCIkEWPHoUPi5dRM7EpFNcrSXYUZM0+azn/GgS9Gx1BBZmDNlNXh10+i0BMcAACAASURBVHEAwLP39MfYAT1ETkRk2/4wsjfsZRIcLqhGXtE1sePYNJYaIgtSp2nE0+uO4LpOj9/188Lce/uLHYnI5nm7OGJShB8A4LN9F0VOY9tYaogshCAIWLDpGM6V18LHVY4Vj0by5G0iM/HnUUEAgB+OlaBc3SByGtvFUkNkIb7KLsTWvGLIpBKsnBEFL2dusEdkLgb3UmBYH3fo9AIPuhQRSw2RBTh6+Rr+3/cnAQDzEwZheKCHyImI6LeeuHG1Zl12IY9OEAlLDZGZU9Xr8PS6I9DqDYgP9cGs0UFiRyKiVkwI84GfwhFVdVp8n1csdhybxFJDZMYMBgEvfJOLy9XX0dujO96dPgQSCe+jITJHdjIplLGBAJqWdwuCIG4gG8RSQ2TG1mRcQPqpcjjYSbHq8SgoutmLHYmIbuOxEQFwtJfiZIkaBy5eFTuOzWGpITJT2Req8O6upoMqlzwYhnB/hciJiOhO3Lo74KGhvQBwMz4xsNQQmaGqWg2e25ADvUHAQ0P98diIALEjEVEbPTEqEACw+2Qpiq7WixvGxrDUEJmZpvto8lCm1qCftzPefCic99EQWZABPi74XT8vGATgP5mXxI5jU1hqiMzMp3sv4KfTFZDbSbFyxlB0d7ATOxIRmejm1ZoNB4tQp2kUN4wNYakhMiM5hdV4Z2fTfTSvPRiKQb6uIiciovYYP9AbgZ7dUdPQiE1HLosdx2aw1BCZCdV1HZ5dn4NGg4AHBvfEjBG9xY5ERO0klUrwx7hAAMBn+y/BYODy7q7AUkNkBprOdTqKy9XXEeDRDSmPDOZ9NEQWbvqwALjI7XChog6/nK0QO45NYKkhMgPrsgux/Vgp7KQS/POxKLg6cj8aIkvnLLfD9GFNKxe5vLtrsNQQiexUiRrLtjWd6/RKwiBEBriJnIiIOsqf4gIhkQA/n6nAufJaseNYPZYaIhHVaxvxzNdHoG00YPzAHvjL73iuE5E16e3ZHfcO8gEAfLH/krhhbABLDZGIXttyAucr6uDjKsffEyMhlfI+GiJr8+cby7u/PXwZqnqduGGsHEsNkUg251zGt4cvQyoB/vHoUHg4OYgdiYg6QWywJwb6uOC6To9vDheJHceqsdQQieBCRS0Wbj4OAHju3v4Y2ddT5ERE1FkkEgmUsX0AAKkHi3h6dydiqSHqYppGPZ75Ogf1Wj1G9vXAs/f0FzsSEXWyyZF+cLSX4mx5LXKKrokdx2qx1BB1seU7TuNkiRoeTg74x6NDIeN9NERWz9XRHveH9wQA/Pcgp6A6C0sNURf66XQ5/r3vIgDg3WkR8HF1FDkREXWVxOFNe9Z8n1fM86A6CUsNURepqNHgxW/yADTtXXFviI/IiYioK8UEeSDQszvqtHr8cKxE7DhWqV2lZtWqVQgKCoKjoyOio6ORkZFx2/EbN25EaGgo5HI5QkNDsXnzZuPXdDodXnnlFQwePBhOTk7w8/PDzJkzUVxc3Ow9AgMDIZFImj3mz5/fnvhEXU4QBLz0bR4qa7UY6OOC+RMHiR2JiLqYRCIx7jDMKajOYXKpSU1NRXJyMhYuXIicnByMHj0aEydORGFhYavjMzMzkZSUBKVSiby8PCiVSiQmJiI7OxsAUF9fjyNHjmDx4sU4cuQINm3ahDNnzmDy5Mkt3mvZsmUoKSkxPhYtWmRqfCJRfL7/En46XQEHOyk+fGwoHO1lYkciIhFMi+4FqQQ4VFDNHYY7gUQwcW1ZTEwMoqKisHr1auNzISEhmDp1KlJSUlqMT0pKglqtxo4dO4zPJSQkwN3dHevXr2/1Mw4ePIgRI0agoKAAvXs3nVQcGBiI5ORkJCcnmxLXSK1WQ6FQQKVSwdXVtV3vQdQep0rUmLJyH7R6A5ZNCcPM2ECxIxGRiGZ9cRDpp8rx5Ji+WHB/iNhxzJ4pv79NulKj1Wpx+PBhxMfHN3s+Pj4e+/fvb/U1mZmZLcZPmDDhluMBQKVSQSKRwM2t+Rk4y5cvh6enJyIjI/Hmm29Cq9Xe8j00Gg3UanWzB1FXa9Dp8dz6HGj1Btw7yBvKkX3EjkREIku8MQW18chl6PQGkdNYFztTBldWVkKv18PHp/kNjj4+PigtLW31NaWlpSaNb2howPz58zFjxoxmjWzu3LmIioqCu7s7Dhw4gAULFuDixYv49NNPW32flJQULF261JRvj6jDvfnDKZwtr0UPFznemRYBiYTLt4ls3fhB3vBylqOyVoP/5ZdjQpiv2JGsRrtuFP7tD2ZBEG77w7qt43U6HR599FEYDAasWrWq2deef/55jB07FhEREZg1axY+/vhjrF27FlVVVa1+5oIFC6BSqYyPoiLelEVdK+1kGb7MKgAAvJ84BJ7OcpETEZE5sJdJ8Ui0PwDeMNzRTCo1Xl5ekMlkLa6ylJeXt7gac5Ovr2+bxut0OiQmJuLixYtIS0u747zZyJEjAQDnzp1r9etyuRyurq7NHkRdpUzdgJe/bVq+PXt0EEb37yFyIiIyJzenoPacLkeZukHkNNbDpFLj4OCA6OhopKWlNXs+LS0NcXFxrb4mNja2xfjdu3c3G3+z0Jw9exbp6enw9LzzOTg5OTkAgJ49e5ryLRB1OoNBwAv/zUN1vQ5hfq54ccJAsSMRkZkJ7uGM4YHuMAhNp3dTxzDpnhoAmDdvHpRKJYYNG4bY2FisWbMGhYWFmDNnDgBg5syZ8Pf3N66Emjt3LsaMGYPly5djypQp2LJlC9LT07F3714AQGNjI6ZNm4YjR45g27Zt0Ov1xis7Hh4ecHBwQGZmJrKysjB+/HgoFAocPHgQzz//PCZPnmxcHUVkLj7dewF7z1Wim70MHz42FHI7Lt8mopYShwXg4KVqfHOoCE+PC+Y9dx3A5FKTlJSEqqoq454x4eHh2L59O/r0aVrVUVhYCKn01wtAcXFx2LBhAxYtWoTFixcjODgYqampiImJAQBcvnwZW7duBQBERkY2+6w9e/Zg3LhxkMvlSE1NxdKlS6HRaNCnTx/Mnj0bL7/8cru/caLOcPyKCu/uOg0AeP3BUAT3cBY5ERGZqwciemLp9ydxqaoe2RevYmTfO89S0O2ZvE+NpeI+NdTZrmv1eOCfGbhQUYeEMF+s/kMU/+ZFRLe1YNNRrD9QhIeH+uP9pMg7v8AGddo+NUR0a29tP4ULFXXwcZUj5eHBLDREdEc3bxjefrwE6gadyGksH0sNUQfYc7rcuHz7velD4O7kIHIiIrIEkQFuGODjjAadAVtzi+/8Arotlhqiu1RVq8HL3x4FADwxKpDLt4mozSQSifFqzX8Pcc+au8VSQ3QXBEHAgk3HUFGjQX9vZ7ySwNO3icg0D0f1gr1MgqOXVThVwiN97gZLDdFd+ObQZew+WQZ7mQQrHo3k6dtEZDIPJwf8PrRpQ9pU7jB8V1hqiNqpoKoOS74/AQB4IX4gwvwUIiciIkt1cwrqu9wr0DTqRU5juVhqiNqhUW/A86m5qNfqMSLIA7NH9xU7EhFZsNH9e8BP4Yhr9TrsPlEmdhyLxVJD1A6rfzqPI4XX4CK3w/uJQyCTcvk2EbWfTCrBtOheAHjD8N1gqSEyUV7RNaz48SwAYOmUMPRy7y5yIiKyBtNvTEHtPVeJch5y2S4sNUQmqNc24vnUXOgNAh6I6ImHhvqLHYmIrESAR3dE9XaDIAA/HCsRO45FYqkhMsFb20/hQmUdfF0d8ebUcO4aTEQd6sEhfgCArXnciK89WGqI2uh/+WX4KqsQQNOuwW7duWswEXWsByJ6QioBcgqvoehqvdhxLA5LDVEbNO0afAwA8OdRQfhdfy+RExGRNfJ2cURscNNp3bxaYzqWGqI7EAQBCzcfR2Vt067BLycMFDsSEVmxyTemoL5nqTEZSw3RHXyXewU7T5TCTirBB0ncNZiIOldCWE/YyyTIL63BmbIaseNYFJYaotsoUV3Ha1uadg2ee29/hPtz12Ai6lyK7vYYO8AbAHhyt4lYaohuwWAQ8NI3R1HT0IghAW54alyw2JGIyEZMjvx1FZQgCCKnsRwsNUS38FV2Afaeq4TcTor3E4fATsb/XIioa9wX4o1u9jIUXq1H3mWV2HEsBn9KE7XiYmUd3tp+CgCwYOIgBPdwFjkREdmS7g52xpO7OQXVdiw1RL/RqDdg3n9z0aAzYFQ/T8yMDRQ7EhHZoJuroLYdLYbewCmotmCpIfqNT365gJwbh1W+O20IpDyskohEMGZAD7g62qG8RoMDF6+KHccisNQQ/R8nilVYkX4GALBkchj83LqJnIiIbJWDnRQTw3sC4EZ8bcVSQ3SDplGPeal50OkFTAjzwcNRPKySiMR1cxXUjuMl0DYaRE5j/lhqiG54P+0MTpfVwNPJAW89NJiHVRKR6Eb29UQPFzmu1euw91yF2HHMHksNEYCDl65izS8XAAApDw+Gp7Nc5ERERIBMKsEDg29MQXEV1B2x1JDNq9M04oX/5kEQgGnRvRAf5it2JCIio5tTULtPluG6Vi9yGvPGUkM2763tp1B4tR7+bt3w2oOhYschImpmaIAbAjy6oV6rx4/5ZWLHMWssNWTTMs5WYF12IQDg3WkRcHW0FzkREVFzEokED0bcODaBU1C3xVJDNkvdoMMr3x4FAMyM7YO4fl4iJyIiat3NKaifTldAdV0nchrzxVJDNuuNbSdRrGpAb4/umD9xkNhxiIhuaZCvKwb4OEOrN2DXiVKx45gtlhqySXvyy/HfQ5chkQDvTR+C7g52YkciIrqtm8cmfM+N+G6JpYZsjqpeh/mbmqad/jwqCCOCPERORER0Zw/eKDX7zlWiokYjchrzxFJDNmfJ9ydQptagr5cTXpowUOw4RERt0sfTCUMC3GAQmnYYppZYasim7DpRis05VyCVAO8lDoGjvUzsSEREbfZgBDfiu512lZpVq1YhKCgIjo6OiI6ORkZGxm3Hb9y4EaGhoZDL5QgNDcXmzZuNX9PpdHjllVcwePBgODk5wc/PDzNnzkRxcfP/w6qrq6FUKqFQKKBQKKBUKnHt2rX2xCcbdbVOi4WbjwEA/jomGFG93UVORERkmgeH+EEiAQ4VVONydb3YccyOyaUmNTUVycnJWLhwIXJycjB69GhMnDgRhYWFrY7PzMxEUlISlEol8vLyoFQqkZiYiOzsbABAfX09jhw5gsWLF+PIkSPYtGkTzpw5g8mTJzd7nxkzZiA3Nxc7d+7Ezp07kZubC6VS2Y5vmWzV4i3HUVmrxQAfZzz/+/5ixyEiMpmPqyNGBDbdB7jrBDfi+y2JIAiCKS+IiYlBVFQUVq9ebXwuJCQEU6dORUpKSovxSUlJUKvV2LFjh/G5hIQEuLu7Y/369a1+xsGDBzFixAgUFBSgd+/eOHXqFEJDQ5GVlYWYmBgAQFZWFmJjY5Gfn4+BA+98X4RarYZCoYBKpYKrq6sp3zJZgW1Hi/HM1zmQSSX47ulRGNxLIXYkIqJ2+ffei1i27SRigjyQ+mSs2HE6nSm/v026UqPVanH48GHEx8c3ez4+Ph779+9v9TWZmZktxk+YMOGW4wFApVJBIpHAzc3N+B4KhcJYaABg5MiRUCgUt3wfjUYDtVrd7EG2qaJGg8XfHQcA/G1cMAsNEVm0+DAfAE0H8V6t04qcxryYVGoqKyuh1+vh4+PT7HkfHx+Ulra+GVBpaalJ4xsaGjB//nzMmDHD2MhKS0vh7e3dYqy3t/ct3yclJcV4/41CoUBAQMAdvz+yPoIgYOHmY6iu1yGkpyueuYfTTkRk2Xq5d0eYnysMApB+ilNQ/1e7bhSWSCTN/lkQhBbPtWe8TqfDo48+CoPBgFWrVt32Pe70uQsWLIBKpTI+ioqKbpmPrNeW3GLsPlkGe5kEf58+BA52XPBHRJYvPtQXALCb99U0Y9JPeC8vL8hkshZXR8rLy1tcjbnJ19e3TeN1Oh0SExNx8eJFpKWlNZs38/X1RVlZy//jKioqbvm5crkcrq6uzR5kW8rVDXh96wkAwHP39EeoH/8MEJF1mBDe9Lsv42wF6rWNIqcxHyaVGgcHB0RHRyMtLa3Z82lpaYiLi2v1NbGxsS3G7969u9n4m4Xm7NmzSE9Ph6enZ4v3UKlUOHDggPG57OxsqFSqW34u2TZBEPDq5uNQXdch3N8Vc8YFix2JiKjDDPRxQW+P7tA0GvDLmQqx45gNkw+8mTdvHpRKJYYNG4bY2FisWbMGhYWFmDNnDgBg5syZ8Pf3N66Emjt3LsaMGYPly5djypQp2LJlC9LT07F3714AQGNjI6ZNm4YjR45g27Zt0Ov1xis7Hh4ecHBwQEhICBISEjB79mx88sknAIC//vWvmDRpUptWPpHt2ZpXjPRTTdNO700fAnsZp52IyHpIJBLEh/rg070XsftEGRLCe4odySyYXGqSkpJQVVWFZcuWoaSkBOHh4di+fTv69OkDACgsLIRU+usvkLi4OGzYsAGLFi3C4sWLERwcjNTUVONKpsuXL2Pr1q0AgMjIyGaftWfPHowbNw4AsG7dOjz33HPGlVSTJ0/GypUrTf+OyeqV1zSfdhrky2knIrI+E8J98enei/gxvxw6vYF/eUM79qmxVNynxjYIgoAnvzyM3SfLEObniu/+Nor/oRORVdIbBIx4Mx1VdVqsmxWDUf28xI7UKTptnxoic7c1r2m1k52U005EZN1kUgnuC2m6YXj3ida3N7E1/IlPVqOiRmOcdnr2nv4I6ckrckRk3W5uxLf7ZBlsZOLltlhqyCoIgoBF3x3DtXodQnu64unxXO1ERNZvVD8vdHeQoUTVgGNXVGLHER1LDVmFbUdLsOtE07TTu9MjOO1ERDbB0V6GcQN7AOBGfABLDVmByloNXtty42yn8f0Q5seznYjIdhh3Fz7J+2pYasjivbbluPFsp7+N7yd2HCKiLjV+kDfspBKcKavFxco6seOIiqWGLNoPR0uw/VjpjdVOETzbiYhsjqKbPWKDm3bit/VVUPwNQBarqlaDxTemnZ7mtBMR2bD40F9XQdkylhqyWK9tPYGrdVoM8nXBM5x2IiIb9vsb99UcKaxGeU2DyGnEw1JDFmnn8RL8cLQEshub7HHaiYhsma/CEUMC3CAIQPrJcrHjiIa/CcjiVNdpsei7pmmnp8YGI9yf005ERL9OQdnufTUsNWRxlm07icpaLfp7O+PZezntREQEABPCmqag9p+rQk2DTuQ04mCpIYvy46kybM65AqkEeGdaBOR2MrEjERGZhX7ezujbwwlavQE/na4QO44oWGrIYqiu6/Dq5mMAgFmj+2Job3eRExERmZdfN+KzzVVQLDVkMVK2n0KZWoMgLyfM+/0AseMQEZmdmwdc7skvh6ZRL3KarsdSQxYh42wFNhwsAgAsfyQCjvacdiIi+q3IXm7wdpGjVtOIzPNVYsfpciw1ZPbqNI2Yv7Fp2umPsX0wIshD5EREROZJKpXg9za8ER9LDZm9d3bm48q16+jl3g0vJwwSOw4RkVmLv7EKKu1kGQwGQeQ0XYulhszagYtX8UVmAYCmaScnuZ3IiYiIzFtsX0+4yO1QUaNBTtE1seN0KZYaMlvXtXq8/G0eAOCxEQEY1c9L5ERERObPwU6K8YO8AdjeAZcsNWS23k87jUtV9fB1dcSC+0PEjkNEZDFu3lez57RtHZnAUkNm6UhhNdbuvQgAeOvhcLg62ouciIjIcozu7wWpBDhTVovia9fFjtNlWGrI7Gga9Xj526MwCMDDQ/1xzyAfsSMREVkUt+4Oxg1KbWl3YZYaMjv//PEczpXXwstZjtceDBU7DhGRRRo3oAcA4CcbmoJiqSGzcqJYhY9/Pg8AeGNqGNy6O4iciIjIMo0b2HSz8L5zldA2GkRO0zVYashsNOoNeGXjUTQaBNw/2BcJ4T3FjkREZLHC/Fzh5eyAOq0ehwquih2nS7DUkNn4V8ZFHL+ihqKbPZZMDhM7DhGRRZNKJRhzYwrqZxu5r4alhszC+YpafJB+BgDw2qRQeLs4ipyIiMjy3ZyCspWbhVlqSHQGg4D5G49C22jA2AE98HCUv9iRiIiswpgbS7tPl9XYxNJulhoS3VfZBTh4qRpODjK8+VA4JBKJ2JGIiKyCW3cHRAa4AQB+PmP9V2tYakhUl6vrsXxHPgDglYmD0Mu9u8iJiIisy69TUNa/tJulhkQjCAJe3XwcdVo9hge64w8xfcSORERkdcYNbLpZeN+5Kqtf2s1SQ6LZdOQKfjlTAQc7Kd5+JAJSKaediIg6WrifAl7ODqjVNOJwQbXYcTpVu0rNqlWrEBQUBEdHR0RHRyMjI+O24zdu3IjQ0FDI5XKEhoZi8+bNzb6+adMmTJgwAV5eXpBIJMjNzW3xHuPGjYNEImn2ePTRR9sTn8xAeU0Dlm07CQBIvq8/gns4i5yIiMg6SaUSjOl/Y3fhM9Y9BWVyqUlNTUVycjIWLlyInJwcjB49GhMnTkRhYWGr4zMzM5GUlASlUom8vDwolUokJiYiOzvbOKaurg6jRo3C22+/fdvPnj17NkpKSoyPTz75xNT4ZCaWbD0B1XUdwvxcMXt0X7HjEBFZtbEDbWO/GokgCIIpL4iJiUFUVBRWr15tfC4kJARTp05FSkpKi/FJSUlQq9XYsWOH8bmEhAS4u7tj/fr1zcZeunQJQUFByMnJQWRkZLOvjRs3DpGRkVixYoUpcY3UajUUCgVUKhVcXV3b9R7UMXYeL8Gcr45AJpVg6zOjEOanEDsSEZFVu1qnRfQbaRAEIHPBPeip6CZ2pDYz5fe3SVdqtFotDh8+jPj4+GbPx8fHY//+/a2+JjMzs8X4CRMm3HL87axbtw5eXl4ICwvDiy++iJqaGpPfg8Slqtdh0XcnAABzxvZloSEi6gIeTg4Y0uvG0m4rvlpjZ8rgyspK6PV6+Pj4NHvex8cHpaWlrb6mtLTUpPG38vjjjyMoKAi+vr44fvw4FixYgLy8PKSlpbU6XqPRQKPRGP9ZrVab9HnUOd744SQqazUI7uGEZ+/pL3YcIiKbMW5gD+QWXcNPpyvw6IjeYsfpFCaVmpt+uzmaIAi33TDN1PGtmT17tvF/h4eHo3///hg2bBiOHDmCqKioFuNTUlKwdOlSkz6DOlfG2Qp8c/gyJBJg+SMRcLSXiR2JiMhmjBvojRXpZ7HvXCV0egPsZda3ANqk78jLywsymazFVZby8vIWV2Nu8vX1NWl8W0VFRcHe3h5nz55t9esLFiyASqUyPoqKiu7q8+ju1GsbsWDTMQDAzJF9MCzQQ+RERES2JcJfAQ8nB9RY8dJuk0qNg4MDoqOjW0z5pKWlIS4urtXXxMbGthi/e/fuW45vqxMnTkCn06Fnz56tfl0ul8PV1bXZg8Tz991ncLn6OvwUjngpYZDYcYiIbE7T0m4vANZ7wKXJ00/z5s2DUqnEsGHDEBsbizVr1qCwsBBz5swBAMycORP+/v7GlVBz587FmDFjsHz5ckyZMgVbtmxBeno69u7da3zPq1evorCwEMXFxQCA06dPA2i6yuPr64vz589j3bp1uP/+++Hl5YWTJ0/ihRdewNChQzFq1Ki7/pdAnSunsBr/3ncRAPDmw4PhLG/XrCcREd2lcQO98V1uMX46XY75E63vL5gm/3ZJSkpCVVUVli1bhpKSEoSHh2P79u3o06dpi/vCwkJIpb9eAIqLi8OGDRuwaNEiLF68GMHBwUhNTUVMTIxxzNatW/HEE08Y//nmpnqvv/46lixZAgcHB/z444/4xz/+gdraWgQEBOCBBx7A66+/DpmM92WYM22jAfM3HoMgAA8N9cf4G2eQEBFR1xszoAckEiC/tAalqgb4KhzFjtShTN6nxlJxnxpx/CP9LD5IPwMPJwekzxsLDycHsSMREdm0KR/tQ17RNSx/ZDCShpv/KqhO26eGyBRny2qwck/TjdyvPxjKQkNEZAbGDbhxZIIV3lfDUkOdQm8Q8MrGo9DpBdwzyBuTh/iJHYmIiPDrqd17zzYt7bYmLDXUKb7MvIQjhdfgLLfDG1PDTd6XiIiIOkdELze4d7dHjaYRR6xsaTdLDXW4y9X1eGdX0wq2VxIGws/Ncs4YISKydjKpBGNuTkGdsa4pKJYa6lCCIGDh5uOo1+oxPNAdj8f0ETsSERH9xs0pKGu7r4alhjrUd7lX8POZCjjIpEh5OAJSKaediIjMzZj+TUu7T5WoUaZuEDtOh2GpoQ5TWavB0u9PAgCeu7cf+nk7i5yIiIha4+ksR4S/AoB1ndrNUkMdZun3J3GtXodBvi54cmyw2HGIiOg2xt7YDPWnM+UiJ+k4LDXUIX48VYbv84ohlQDvTIuwytNfiYisyc37ajLOVqLRSpZ28zcP3bVaTSMWfXccAPCX3wUhopebyImIiOhOhvRyg6KbPWoaGnH0ikrsOB2CpYbu2rs781GiakBvj+6Y9/uBYschIqI2kEklGNnXAwCQeb5K5DQdg6WG7srhgqv4T1YBAOCthwajmwMPGCUishSj+nkBAPadqxQ5ScdgqaF20zTq8cqNE7inRffC7/p7iR2JiIhMEBfc9HP7UEE1GnR6kdPcPZYaarfVP53HufJaeDk7YNEDIWLHISIiEwX3cIK3ixzaRoNVHJnAUkPtcrasBh/tOQcAWDI5DG7deQI3EZGlkUgkv05Bnbf8KSiWGjLZ/z2B+74QbzwwuKfYkYiIqJ1igz0BAPut4GZhlhoy2VdZBcYTuP8fT+AmIrJoN6/UHL2sQk2DTuQ0d4elhkxy5dp1vLMzH0DTCdw9FTyBm4jIkvm7dUOgZ3foDQIOXLwqdpy7wlJDbSYIAhZtPoY6rR7RfXgCNxGRtYgNvrm027KnoFhqqM2+P1qCPaebTuB+++HBPIGbiMhKjOp3874ay75ZmKWG2qS6ToulW08AAP42vh/6+7iInIiIiDpKbN+mPahTqgAAHp9JREFUUpNfWoPKWo3IadqPpYba5I0fTqGqTosBPs54ahxP4CYisiaeznIM8m36y6olH5nAUkN3lHG2AhuPXIZEArz9SAQc7PjHhojI2txcBWXJS7v524luq17biFc3HwMA/DE2EFG93UVOREREnSEu2PLvq2GpodtakX4WRVevw0/hiBcn8ARuIiJrNSLIAzKpBAVV9bhcXS92nHZhqaFbOnZZhU8zLgAA3ngoHM5yO5ETERFRZ3FxtMeQXgoAljsFxVJDrWrUGzB/01EYBODBIX64Z5CP2JGIiKiT3Ty1e/85y5yCYqmhVq3dexEnitVQdLPHa5NCxY5DRERdIK7fr+dACYIgchrTsdRQCwVVdfgg/QwAYOEDIejhIhc5ERERdYWo3u6Q20lRXqPB+YpaseOYjKWGmhEEAa9uPoYGnQFxwZ6YHt1L7EhERNRFHO1lGBbYtMrVEo9MYKmhZjYeuYJ956ogt5PirYcG8wRuIiIbY7yvxgKXdrPUkFFlrQZv/HASAJB83wAEejmJnIiIiLrazf1qMs9XQW+wrPtqWGrIaNn3J3GtXoeQnq6YNTpI7DhERCSCwf4KuMjtoG5oxMlitdhxTMJSQwCAPafLsTWvGFIJsPyRwbCX8Y8GEZEtspNJEdPXAwCwz8KmoNr1m2vVqlUICgqCo6MjoqOjkZGRcdvxGzduRGhoKORyOUJDQ7F58+ZmX9+0aRMmTJgALy8vSCQS5ObmtngPjUaDZ599Fl5eXnBycsLkyZNx+fLl9sSn36jTNOL/t3f/UVGWef/A38PADKIwIggDhIhUCoHfBBNBTStDLX9gu19wbSe/PbWrqz3raPsYBqZZ8sN23Z7nmO7RfE6d3Y7wLUPbFl3wKU0FzVgwDNIMEU0mRXEGfzHMcD1/sIw7AcqMwj1z836dM8e4/dz3fOYjp/mc677u68oqPA4A+LcJkRh932CJMyIiIindmlfjXpOFHW5qCgoKoNfrkZmZiYqKCkyaNAkzZsxAfX19l/FlZWVIT0+HTqfDsWPHoNPpkJaWhiNHjthirl27hgkTJiA3N7fb99Xr9SgsLER+fj4OHjyIq1evYubMmbBarY5+BPqJ3xefwA9XbuA+/wFYnvKg1OkQEZHEOja3PHr6MsyWNomz6TmFcHB1ncTERMTHx2Pz5s22Y9HR0UhNTUVOTk6n+PT0dJhMJuzevdt2bPr06fD398f27dvtYuvq6hAZGYmKigo8/PDDtuNGoxFDhw7Fn//8Z6SnpwMAzp8/j/DwcBQVFWHatGl3zNtkMkGj0cBoNMLPz8+RjyxrlWevYO6mQxACeP/fxmHyg0OlTomIiCQmhMAj6/ai8aoZBb8ej8QRAZLl4sj3t0MjNWazGeXl5UhJSbE7npKSgtLS0i7PKSsr6xQ/bdq0buO7Ul5ejtbWVrvrhIaGIjY2ttvrtLS0wGQy2b3IXqu1DRk7voYQwNwxYWxoiIgIAKBQKJDkhregHGpqGhsbYbVaERxsvw9QcHAwDAZDl+cYDAaH4ru7hkqlgr+/f4+vk5OTA41GY3uFh4f3+P36iy1f1OJbQzP8fbyQ9XS01OkQEZELmRDVsWWC+0wWdmqi8E8XZBNC3HaRNkfje+p211m5ciWMRqPtdfbs2bt+Pzk53XgN//k/3wEAVs2MQcAgboVARES3dEwWrqi/gutmi8TZ9IxDTU1gYCCUSmWn0ZELFy50Go3poNVqHYrv7hpmsxlNTU09vo5arYafn5/di9oJIfDqx1UwW9ow6YFAzB0TJnVKRETkYoYF+OA+/wGwtAl8efqy1On0iENNjUqlQkJCAkpKSuyOl5SUIDk5uctzkpKSOsUXFxd3G9+VhIQEeHl52V2noaEBx48fd+g61O7D8nMoq70Eby8PrEvlVghERNS15Khbu3a7A09HT1i+fDl0Oh3Gjh2LpKQkbNmyBfX19Vi0aBEA4LnnnkNYWJjtSailS5fi0UcfRV5eHubMmYNdu3Zh7969OHjwoO2aly9fRn19Pc6fPw8AOHHiBID2ERqtVguNRoMXXngBL7/8MgICAjBkyBD87ne/Q1xcHKZOnXrXRehPLja3YN3fagAAy6Y+iGEBPhJnRERErmrC/YH4/1+dc5t5NQ43Nenp6bh06RLWrl2LhoYGxMbGoqioCBEREQCA+vp6eHjcGgBKTk5Gfn4+srKysGrVKkRFRaGgoACJiYm2mE8++QTPP/+87ed58+YBAFavXo01a9YAAP74xz/C09MTaWlpuHHjBp544gm89957UCqVTn3w/mrtp9Uw3mjFQ6F+eGEit0IgIqLuJf1zpOab8yZcuW7GYB+VxBndnsPr1LgrrlMDfP7tBTz/3lF4KIBdSyYi7j6N1CkREZGLe3LDfnx34So2PxuPGXEhff7+vbZODbmvay0WZO28tRUCGxoiIuqJjnk1R9xgsjCbmn7iD8UnuRUCERE57JHI9s0t3eEJKDY1/cCxs1fwXulpAMCbqbHwUTk8lYqIiPqpccPbm5oagwmmm60SZ3N7bGpkrtXahoyPq9AmgDkPh2LKyCCpUyIiIjcS5OeNiAAfCAGUn2m68wkSYlMjc+8eOI2aBhMG+3hh1cwYqdMhIiI39Mg/R2uOuvgtKDY1MlbXeA1v7z0JAMh6OgaB3AqBiIic0HEL6mgdmxqSgBACmTur0GJpw4T7A/CzeG6FQEREzumYLHzsrBE3W60SZ9M9NjUyteMfP+DQqUtQe3IrBCIiujvDA3wQOEgNs7UNX58zSp1Ot9jUyNClqy1482/VAAD91AcxPHCgxBkREZE7UygUGBfpD8C1b0GxqZGhNz6txpXrrYgO8cOLk7gVAhER3b2OycKuvF4NmxqZ2X/yInZWnoeHAsh9Jg5eSv4TExHR3etoav5xpgnWNtfcYYnfeDJy3WxBZmEVAOD/JUfi/4QPljgjIiKSi+gQP/iqPdHcYkFNg0nqdLrEpkZG3t77Hc413UDY4AF4mVshEBHRPaT0UCA+wrXn1bCpkYnjPxjx7oFaAMAbqQ9hoJpbIRAR0b01LtK116thUyMDFmsbVv5zK4SZo0Pw+KhgqVMiIiIZujVZuAlCuN68GjY1MvBeaR2qfjDCz9sTr83iVghERNQ7Rt+ngcrTA41XW1B36brU6XTCpsbNnb18HX8obt8K4dWnohHk6y1xRkREJFfeXko8fF/7QyiuuA8Umxo3JoRA1s7juNFqRWLkEKQ/Ei51SkREJHOP/HMRvi9dcF4Nmxo39smx89h/8iJUnh7IfoZbIRARUe97xIU3t2RT46auXDdj7V/bt0L498fuR9TQQRJnRERE/UFChD88FMCZS9dxwXRT6nTssKlxU+v+VoNL18x4IGgQFk6OkjodIiLqJ3y9vRAd4gfA9W5BsalxQ6WnGvFh+TkoFEDuz+Kg8uQ/IxER9R3bLSgXmyzMb0M3c7PVisydxwEAv0yMQELEEIkzIiKi/qZjEb4v65okzsQemxo3s/GzUzjdeA3Bfmr8x/SRUqdDRET9UMdIzbcGE4w3WiXO5hY2NW7khKEZf9r/PQDg9dmx8PP2kjgjIiLqj4b6qhEZOBBCtO/a7SrY1LgJa5tAxsdfw9ImkBITjOmxWqlTIiKifuyR4e3r1RxxoXk1bGrcxAdHzqCi/goGqT3x+pyHpE6HiIj6OVdcr4ZNjRtoMN7A+j0nAAArpo9EiGaAxBkREVF/1zFZ+OtzV3Cz1SpxNu3Y1LiB1bu+wdUWC8YMG4xnEyOkToeIiAjDhvggyFeNVqtA5dkrUqcDgE2Ny9tz3IDi6h/h6aFAzjNxUHpwKwQiIpKeQqHAI5GutV4NmxoXZrrZitWftK9Js3DyCIzS+kmcERER0S3jhnesV8Omhu7grT0n8KOpBcMDfPDvjz8gdTpERER2OiYL/+NMEyzWNomzYVPjssrPXMZfjpwBAGTPjYO3l1LijIiIiOyN1PrC19sT18xW1DQ0S50OmxpXZLa0YeXHVRAC+L8J9yH5/kCpUyIiIupE6aHA2Ij29Wpc4RaUU03Npk2bEBkZCW9vbyQkJODAgQO3jd+xYwdiYmKgVqsRExODwsJCu78XQmDNmjUIDQ3FgAEDMGXKFHzzzTd2McOHD4dCobB7ZWRkOJO+y9vyxfc4+eNVBAxU4dWnoqVOh4iIqFuuNFnY4aamoKAAer0emZmZqKiowKRJkzBjxgzU19d3GV9WVob09HTodDocO3YMOp0OaWlpOHLkiC1m/fr12LBhAzZu3IijR49Cq9XiySefRHOz/VDW2rVr0dDQYHtlZWU5mr7Lq714Ff/12SkAwGuzYuA/UCVxRkRERN0b9y+L8AkhJM1FIRzMIDExEfHx8di8ebPtWHR0NFJTU5GTk9MpPj09HSaTCbt377Ydmz59Ovz9/bF9+3YIIRAaGgq9Xo9XXnkFANDS0oLg4GDk5eVh4cKFANpHavR6PfR6vVMf1GQyQaPRwGg0ws/PNZ8iEkJg3pbDOHL6Mh59cCjef/4RKBR8hJuIiFxXi8WKuDXFMFva8D8vT0bU0EH39PqOfH87NFJjNptRXl6OlJQUu+MpKSkoLS3t8pyysrJO8dOmTbPFnz59GgaDwS5GrVZj8uTJna6Zl5eHgIAAPPzww1i3bh3MZnO3uba0tMBkMtm9XN2HX53DkdOX4e3lgXWpsWxoiIjI5ak9lXg4fDAA6W9BOdTUNDY2wmq1Ijg42O54cHAwDAZDl+cYDIbbxnf8eadrLl26FPn5+fj888/x0ksv4e2338bixYu7zTUnJwcajcb2Cg8P7/kHlcDF5hasK6oBACyb+iDCh/hInBEREVHPuMp6NZ7OnPTTEQQhxG1HFXoSf6eYZcuW2f579OjR8Pf3x89//nPb6M1PrVy5EsuXL7f9bDKZXLqxeePTahhvtCImxA8vTIyUOh0iIqIeGxc5BPhc+s0tHWpqAgMDoVQqO43KXLhwodNISwetVnvbeK1WC6B9xCYkJKRH1wSA8ePHAwBOnTrVZVOjVquhVqt78Kmkt+/EBXxy7Dw8FEDuz+LgqeST9kRE5D7iI/wxc3QIxkUOQVubgIdEW/o49O2pUqmQkJCAkpISu+MlJSVITk7u8pykpKRO8cXFxbb4yMhIaLVauxiz2Yz9+/d3e00AqKioAAC7RsgdXTdbkLWzfSuE5ydEYvR9gyXOiIiIyDGD1J7YOD8ezyUNl6yhAZy4/bR8+XLodDqMHTsWSUlJ2LJlC+rr67Fo0SIAwHPPPYewsDDbk1BLly7Fo48+iry8PMyZMwe7du3C3r17cfDgQQDtt530ej2ys7PxwAMP4IEHHkB2djZ8fHwwf/58AO2TjQ8fPozHHnsMGo0GR48exbJlyzB79mwMGzbsXtVCEm/v/Q7nmm4gbPAALH/yQanTISIiclsONzXp6em4dOmSbc2Y2NhYFBUVISIiAgBQX18PD49bA0DJycnIz89HVlYWVq1ahaioKBQUFCAxMdEWs2LFCty4cQOLFy9GU1MTEhMTUVxcDF9fXwDtt5IKCgrw+uuvo6WlBREREfjVr36FFStW3O3nl9TxH4x490AtAODN1FgMVDs1xYmIiIjgxDo17srV1qmxWNswd1Mpqn4wYuboEGycHy91SkRERC6n19apoXvnvdI6VP1ghJ+3J16bFSN1OkRERG6PTY0Ezl6+jj8UnwQAvPpUNIJ8vSXOiIiIyP2xqeljQgis2nUcN1qtGBc5BGljXXftHCIiInfCpqaPffp1A/aduAiV0gPZc+MkffSNiIhITtjU9CHj9Va8/tdvAACLH4vC/UH3dtMvIiKi/oxNTR/K2V2DxqtmRA0diN9MiZI6HSIiIllhU9NHDtdeQv7RswCAnGdGQ+2plDgjIiIieWFT0wdutlrxamEVAOAX44a1b/xFRERE9xSbmj6wad/3qL14DUN91ciYMUrqdIiIiGSJTU0v++7HZmzedwoAsGbWQ9AM8JI4IyIiInliU9OL2toEXi2sQqtV4IlRQXgqTit1SkRERLLFpqYX5R89i6N1TfBRKbE2NRYKBdekISIi6i1sanrJBdNN5OyuAQD8LmUkwgYPkDgjIiIieWNT00te/2s1mm9aMPo+DRYkD5c6HSIiItljU9ML9lb/iL9VNUDpoUDOM3FQcisEIiKiXsem5h672mLBa7uOAwBenBSJh0I1EmdERETUP7Cpucf+UHwC5403ET5kAPRPPCh1OkRERP0Gm5p76NjZK3i/tA4AsC41DgNU3AqBiIior7CpuUdarW3I+LgKbQKYOyYMjz44VOqUiIiI+hU2NffIfx88jZoGEwb7eCHr6Wip0yEiIup32NTcA/WXruOPe08CALKejkHAILXEGREREfU/bGrukhACmTurcLO1DclRAfhZfJjUKREREfVLbGruUnH1jzjwXSNUnh5YNzeOWyEQERFJxFPqBNzd46OCkDFjFJQKBSIDB0qdDhERUb/FpuYueSk9sGhylNRpEBER9Xu8/URERESywKaGiIiIZIFNDREREckCmxoiIiKSBTY1REREJAtsaoiIiEgW2NQQERGRLLCpISIiIllwqqnZtGkTIiMj4e3tjYSEBBw4cOC28Tt27EBMTAzUajViYmJQWFho9/dCCKxZswahoaEYMGAApkyZgm+++cYupqmpCTqdDhqNBhqNBjqdDleuXHEmfSIiIpIhh5uagoIC6PV6ZGZmoqKiApMmTcKMGTNQX1/fZXxZWRnS09Oh0+lw7Ngx6HQ6pKWl4ciRI7aY9evXY8OGDdi4cSOOHj0KrVaLJ598Es3NzbaY+fPno7KyEnv27MGePXtQWVkJnU7nxEcmIiIiOVIIIYQjJyQmJiI+Ph6bN2+2HYuOjkZqaipycnI6xaenp8NkMmH37t22Y9OnT4e/vz+2b98OIQRCQ0Oh1+vxyiuvAABaWloQHByMvLw8LFy4EDU1NYiJicHhw4eRmJgIADh8+DCSkpLw7bffYuTIkXfM22QyQaPRwGg0ws/Pz5GPTERERBJx5PvboZEas9mM8vJypKSk2B1PSUlBaWlpl+eUlZV1ip82bZot/vTp0zAYDHYxarUakydPtsWUlZVBo9HYGhoAGD9+PDQaTbfv29LSApPJZPciIiIi+XKoqWlsbITVakVwcLDd8eDgYBgMhi7PMRgMt43v+PNOMUFBQZ2uHRQU1O375uTk2ObfaDQahIeH9+ATEhERkbtyapduhUJh97MQotMxR+PvFNPV9W/3vitXrsTy5cttPxuNRgwbNowjNkRERG6k43u7J7NlHGpqAgMDoVQqO42OXLhwodNISwetVnvbeK1WC6B9NCYkJKTbmB9//LHTtS9evNjt+6rVaqjVatvPHUXhiA0REZH7aW5uhkajuW2MQ02NSqVCQkICSkpKMHfuXNvxkpISzJkzp8tzkpKSUFJSgmXLltmOFRcXIzk5GQAQGRkJrVaLkpISjBkzBkD73J39+/cjLy/Pdg2j0Ygvv/wS48aNAwAcOXIERqPRdp07CQ0NxdmzZ+Hr63vbUaV/ZTKZEB4ejrNnz3JycR9i3aXBukuDdZcG6y4NZ+ouhEBzczNCQ0N7FOyQ/Px84eXlJbZt2yaqq6uFXq8XAwcOFHV1dUIIIXQ6ncjIyLDFHzp0SCiVSpGbmytqampEbm6u8PT0FIcPH7bF5ObmCo1GIz7++GNRVVUlfvGLX4iQkBBhMplsMdOnTxejR48WZWVloqysTMTFxYmZM2c6mr5DjEajACCMRmOvvg/ZY92lwbpLg3WXBusujd6uu8NzatLT03Hp0iWsXbsWDQ0NiI2NRVFRESIiIgAA9fX18PC4Nf84OTkZ+fn5yMrKwqpVqxAVFYWCggK7J5lWrFiBGzduYPHixWhqakJiYiKKi4vh6+tri/nggw/w29/+1vaU1OzZs7Fx40ZH0yciIiKZcnidmv6Ea9tIg3WXBusuDdZdGqy7NHq77so1a9asuedXlRGlUokpU6bA09OpB8XISay7NFh3abDu0mDdpdGbdedIDREREckCd+kmIiIiWWBTQ0RERLLApoaIiIhkgU0NERERyUK/b2o2bdqEyMhIeHt7IyEhAQcOHLht/I4dOxATEwO1Wo2YmBgUFhb2Uaby4kjdt27dikmTJsHf3x/+/v6YOnUqvvzyyz7MVj4c/X3vkJ+fD4VCgdTU1F7OUJ4crfuVK1ewZMkShISEwNvbG9HR0SgqKuqjbOXD0bq//fbbGDlyJAYMGIDw8HAsW7YMN2/e7KNs5eGLL77ArFmzEBoaCoVCgZ07d97xnP379yMhIQHe3t4YMWIE/vSnPzmfQK8s6ecmOlZH3rp1q6iurhZLly4VAwcOFGfOnOkyvrS0VCiVSpGdnS1qampEdnZ2p9WR6c4crfv8+fPFO++8IyoqKkRNTY14/vnnhUajEefOnevjzN2bo3XvUFdXJ8LCwsSkSZPEnDlz+ihb+XC07i0tLWLs2LHiqaeeEgcPHhR1dXXiwIEDorKyso8zd2+O1v0vf/mLUKvV4oMPPhCnT58Wf//730VISIjQ6/V9nLl7KyoqEpmZmWLHjh0CgCgsLLxtfG1trfDx8RFLly4V1dXVYuvWrcLLy0t89NFHTr1/v25qxo0bJxYtWmR3bNSoUXbbPPyrtLQ0MX36dLtj06ZNE/Pmzeu1HOXI0br/lMViEb6+vuL999/vjfRky5m6WywWMWHCBPHuu++KBQsWsKlxgqN137x5sxgxYoQwm819kZ5sOVr3JUuWiMcff9zu2PLly8XEiRN7LUe560lTs2LFCjFq1Ci7YwsXLhTjx4936j377e0ns9mM8vJy27YLHVJSUlBaWtrlOWVlZZ3ip02b1m08deZM3X/q+vXraG1txZAhQ3ojRVlytu5r167F0KFD8cILL/R2irLkTN0/+eQTJCUlYcmSJQgODkZsbCyys7NhtVr7ImVZcKbuEydORHl5ue3Wdm1tLYqKivD000/3er79WXffq1999RVaW1sdvl6/XUaxsbERVqsVwcHBdseDg4NhMBi6PMdgMDgUT505U/efysjIQFhYGKZOndobKcqSM3U/dOgQtm3bhsrKyr5IUZacqXttbS0+++wzPPvssygqKsJ3332HJUuWwGKx4LXXXuuLtN2eM3WfN28eLl68iIkTJ0IIAYvFgt/85jfIyMjoi5T7re6+Vy0WCxobGxESEuLQ9fptU9NBoVDY/SyE6HTsbuKpa87Wcf369di+fTv27dsHb2/v3kpPtnpa9+bmZvzyl7/E1q1bERgY2FfpyZYjv+9tbW0ICgrCli1boFQqkZCQgPPnz+Ott95iU+MgR+q+b98+rFu3Dps2bUJiYiJOnTqFpUuXIiQkBKtWreqLdPutrv6dujreE/22qQkMDIRSqezUtV+4cKFT19hBq9U6FE+dOVP3Dr///e+RnZ2NvXv3YvTo0b2Zpuw4Wvfvv/8edXV1mDVrlu1YW1sbAMDT0xMnTpxAVFRU7yYtA878voeEhMDLywtKpdJ2LDo6GgaDAWazGSqVqldzlgNn6r5q1SrodDq8+OKLAIC4uDhcu3YNv/71r5GZmQkPj347W6NXdfe96unpiYCAAIev12//lVQqFRISElBSUmJ3vKSkBMnJyV2ek5SU1Cm+uLi423jqzJm6A8Bbb72FN954A3v27MHYsWN7O03ZcbTuo0aNQlVVFSorK22v2bNn47HHHkNlZSXCw8P7KnW35szv+4QJE3Dq1ClbEwkAJ0+eREhICBuaHnKm7tevX+/UuCiVSoj2B2p6Ldf+rrvv1bFjx8LLy8vxCzo1vVgmOh7527Ztm6iurhZ6vV4MHDhQ1NXVCSGE0Ol0djPlDx06JJRKpcjNzRU1NTUiNzeXj3Q7wdG65+XlCZVKJT766CPR0NBgezU3N0v1EdySo3X/KT795BxH615fXy8GDRokXnrpJXHixAnx6aefiqCgIPHmm29K9RHckqN1X716tfD19RXbt28XtbW1ori4WERFRYm0tDSpPoJbam5uFhUVFaKiokIAEBs2bBAVFRW2R+kzMjKETqezxXc80r1s2TJRXV0ttm3bxke678Y777wjIiIihEqlEvHx8WL//v22v5s8ebJYsGCBXfyHH34oRo4cKby8vMSoUaPEjh07+jhjeXCk7hEREQJAp9fq1av7PnE35+jv+79iU+M8R+teWloqEhMThVqtFiNGjBDr1q0TFoulj7N2f47UvbW1VaxZs0ZERUUJb29vER4eLhYvXiyampokyNx9ff75513+/7qj1gsWLBCTJ0+2O2ffvn1izJgxQqVSieHDh4vNmzc7/f4KITiuRkRERO6v386pISIiInlhU0NERESywKaGiIiIZIFNDREREckCmxoiIiKSBTY1REREJAtsaoiIiEgW2NQQERGRLLCpISIiIllgU0NERESywKaGiIiIZIFNDREREcnC/wJnKmecV0MsggAAAABJRU5ErkJggg==",
      "text/plain": [
       "Figure(PyObject <Figure size 640x480 with 1 Axes>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "1-element Array{PyCall.PyObject,1}:\n",
       " PyObject <matplotlib.lines.Line2D object at 0x13ad57b00>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using PyPlot\n",
    "\n",
    "plot(xgrid[2:end-1], U)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Poisson equation in three dimensions\n",
    "\n",
    "$$-(u_{xx}+u_{yy}+u_{zz})=f(x,y,z),~~~0<x,y,z<1$$\n",
    "\n",
    "with zero Dirichlet boundary conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAGiCAYAAACh/hJSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAcoUlEQVR4nO3dX3BU9f3/8dchMcufSdYmDolbEwgzjCggKlFHQIVR00EGSx21iEaqNzKNQoxDIUXbaIdssS1fbFNw4oXQYVBuJFJbR1PFAIN/QkLU2o5IzUCUZlI7zm6AukByfhcd9teYkD9wzr73hOdjZsfZs4dzXnM2uy8/ez571nFd1xUAACk2yjoAAODCRAEBAExQQAAAExQQAMAEBQQAMEEBAQBMUEAAABMUEADABAUEADBBAQEATAS6gDZu3Kji4mKNHj1aM2fO1J49e6wjDSoajeq6665Tdna2xo8fr0WLFunTTz+1jjVs0WhUjuOooqLCOsqgvvzySz3wwAPKy8vT2LFjdfXVV6u5udk61oBOnz6tJ598UsXFxRozZowmTZqkZ555Rj09PdbR+ti9e7cWLlyoSCQix3FUX1/f63HXdVVdXa1IJKIxY8Zo7ty5+uSTT4zS/n8D5T516pRWrVql6dOna9y4cYpEInrwwQd19OhRw8T/Ndjx/l+PPPKIHMfRhg0bUphw6AJbQNu3b1dFRYXWrFmjAwcO6KabbtL8+fN15MgR62gDamxsVHl5ud577z01NDTo9OnTKi0t1fHjx62jDVlTU5Pq6up01VVXWUcZ1Ndff63Zs2froosu0uuvv66//e1v+s1vfqOLL77YOtqA1q1bp+eff161tbX6+9//rmeffVa/+tWv9Lvf/c46Wh/Hjx/XjBkzVFtb2+/jzz77rNavX6/a2lo1NTWpoKBAt99+u7q6ulKctLeBcp84cUItLS166qmn1NLSoldeeUUHDx7UnXfeaZC0t8GO9xn19fV6//33FYlEUpTsHLgBdf3117vLli3rtWzKlCnu6tWrjRKdm87OTleS29jYaB1lSLq6utzJkye7DQ0N7i233OKuWLHCOtKAVq1a5c6ZM8c6xrAtWLDAffjhh3stu+uuu9wHHnjAKNHQSHJ37NiRvN/T0+MWFBS4v/zlL5PLvvnmGzccDrvPP/+8RcR+fTt3fz744ANXknv48OEUpRrc2XJ/8cUX7ne/+133r3/9qzthwgT3//7v/wzSDS6QI6CTJ0+qublZpaWlvZaXlpZq3759RqnOTSwWkyTl5uYaJxma8vJyLViwQLfddpt1lCHZuXOnSkpKdM8992j8+PG65ppr9MILL1jHGtScOXP01ltv6eDBg5KkDz/8UHv37tUdd9xhnGx42tra1NHR0eu1GgqFdMsttwTyteo4TtqPnnt6elRWVqaVK1dq6tSp1nEGlGkd4Fx89dVX6u7uVn5+fq/l+fn56ujoMEo1fK7rqrKyUnPmzNG0adOs4wzq5ZdfVktLi5qamqyjDNnnn3+uTZs2qbKyUj/96U/1wQcfaPny5QqFQnrwwQet453VqlWrFIvFNGXKFGVkZKi7u1tr167VfffdZx1tWM68Hvt7rR4+fNgi0jn55ptvtHr1ai1ZskQ5OTnWcQa0bt06ZWZmavny5dZRBhXIAjrDcZxe913X7bMsnT366KP66KOPtHfvXusog2pvb9eKFSv05ptvavTo0dZxhqynp0clJSWqqamRJF1zzTX65JNPtGnTprQuoO3bt2vr1q3atm2bpk6dqtbWVlVUVCgSiWjp0qXW8YYtyK/VU6dOafHixerp6dHGjRut4wyoublZzz33nFpaWgJxfAP5Edwll1yijIyMPqOdzs7OPv+nla4ee+wx7dy5U7t27dJll11mHWdQzc3N6uzs1MyZM5WZmanMzEw1Njbqt7/9rTIzM9Xd3W0dsV+XXnqprrzyyl7LrrjiirSfrLJy5UqtXr1aixcv1vTp01VWVqbHH39c0WjUOtqwFBQUSFJgX6unTp3Svffeq7a2NjU0NKT96GfPnj3q7OxUUVFR8nV6+PBhPfHEE5o4caJ1vD4CWUBZWVmaOXOmGhoaei1vaGjQrFmzjFINjeu6evTRR/XKK6/o7bffVnFxsXWkIbn11lv18ccfq7W1NXkrKSnR/fffr9bWVmVkZFhH7Nfs2bP7THM/ePCgJkyYYJRoaE6cOKFRo3q/PDMyMtJyGvZAiouLVVBQ0Ou1evLkSTU2Nqb9a/VM+Xz22Wf6y1/+ory8POtIgyorK9NHH33U63UaiUS0cuVKvfHGG9bx+gjsR3CVlZUqKytTSUmJbrzxRtXV1enIkSNatmyZdbQBlZeXa9u2bXr11VeVnZ2d/D/DcDisMWPGGKc7u+zs7D7nqcaNG6e8vLy0Pn/1+OOPa9asWaqpqdG9996rDz74QHV1daqrq7OONqCFCxdq7dq1Kioq0tSpU3XgwAGtX79eDz/8sHW0Po4dO6ZDhw4l77e1tam1tVW5ubkqKipSRUWFampqNHnyZE2ePFk1NTUaO3aslixZYph64NyRSER33323Wlpa9Nprr6m7uzv5Ws3NzVVWVpZV7EGP97eL8qKLLlJBQYEuv/zyVEcdnO0kvPPz+9//3p0wYYKblZXlXnvttYGYyiyp39uLL75oHW3YgjAN23Vd949//KM7bdo0NxQKuVOmTHHr6uqsIw0qHo+7K1ascIuKitzRo0e7kyZNctesWeMmEgnraH3s2rWr37/ppUuXuq7736nYP//5z92CggI3FAq5N998s/vxxx/bhnYHzt3W1nbW1+quXbvSNnd/0nkatuO6rpuirgMAICmQ54AAAMFHAQEATFBAAAATFBAAwAQFBAAwQQEBAExQQAAAE4EvoEQioerqaiUSCesow0Lu1Apqbim42cmdWkHMHfgvosbjcYXDYcVisbS/UOD/IndqBTW3FNzs5E6tIOYO/AgIABBMFBAAwETaXQ27p6dHR48eVXZ29pB+UCkej/f6b1CQO7WCmlsKbnZyp1Y65XZdV11dXYpEIn1+VuR/pd05oC+++EKFhYXWMQAA56m9vX3AH9xMuxFQdna2dYRzFovFrCOck3A47Ov2g3pcAJybeDyuwsLCQd/P066AgvA75mcTlJknqcZxAS5Mg72fMwkBAGCCAgIAmKCAAAAmKCAAgAnfCmjjxo0qLi7W6NGjNXPmTO3Zs8evXQEAAsiXAtq+fbsqKiq0Zs0aHThwQDfddJPmz5+vI0eO+LE7AEAA+fJF1BtuuEHXXnutNm3alFx2xRVXaNGiRYpGowP+2zMX1AuiNPtO75D5PfU9qMcFwLkZ6oVRPR8BnTx5Us3NzSotLe21vLS0VPv27euzfiKRUDwe73UDAIx8nhfQV199pe7ubuXn5/danp+fr46Ojj7rR6NRhcPh5I3L8ADAhcG3SQjf/ljHdd1+P+qpqqpSLBZL3trb2/2KBABII55fiueSSy5RRkZGn9FOZ2dnn1GRJIVCIYVCIa9jAADSnOcjoKysLM2cOVMNDQ29ljc0NGjWrFle7w4AEFC+XIy0srJSZWVlKikp0Y033qi6ujodOXJEy5Yt82N3AIAA8qWAfvjDH+rf//63nnnmGf3zn//UtGnT9Oc//1kTJkzwY3cAgABKux+k43tAqcf3gAB4yex7QAAADAUFBAAwQQEBAExQQAAAE77MgrtQ+Xky388T+X5PEgjqcQHgL0ZAAAATFBAAwAQFBAAwQQEBAExQQAAAExQQAMAEBQQAMEEBAQBMUEAAABMUEADABAUEADBBAQEATFBAAAATFBAAwAQFBAAwQQEBAExQQAAAExQQAMAEBQQAMEEBAQBMUEAAABMUEADARKZ1gLOJxWLKycnxfLuO43i+zVTwM7frur5t2+/tB/m4ABc6RkAAABMUEADABAUEADBBAQEATFBAAAATFBAAwAQFBAAw4XkBRaNRXXfddcrOztb48eO1aNEiffrpp17vBgAQcJ4XUGNjo8rLy/Xee++poaFBp0+fVmlpqY4fP+71rgAAAea4Pn/d+1//+pfGjx+vxsZG3XzzzYOuH4/HFQ6HuRJCCgX5G/9cCQFIP0N9H/f9UjyxWEySlJub2+/jiURCiUQieT8ej/sdCQCQBnydhOC6riorKzVnzhxNmzat33Wi0ajC4XDyVlhY6GckAECa8PUjuPLycv3pT3/S3r17ddlll/W7Tn8joMLCQj6CS6Egf9TER3BA+jH/CO6xxx7Tzp07tXv37rOWjySFQiGFQiG/YgAA0pTnBeS6rh577DHt2LFD77zzjoqLi73eBQBgBPC8gMrLy7Vt2za9+uqrys7OVkdHhyQpHA5rzJgxXu8OABBQnp8DOttn8i+++KJ+9KMfDfrvmYadekE+18E5ICD9mJ0D4kULABgKrgUHADBBAQEATFBAAAATFBAAwITv14JLN35OkmCGXeoF9flksg7ACAgAYIQCAgCYoIAAACYoIACACQoIAGCCAgIAmKCAAAAmKCAAgAkKCABgggICAJiggAAAJiggAIAJCggAYIICAgCYoIAAACYoIACACQoIAGCCAgIAmKCAAAAmKCAAgAkKCABgggICAJjItA5wNuFw2Jftuq7ry3b93rafHMexjnDOgvp8+nnMg/p3iAsPIyAAgAkKCABgggICAJiggAAAJiggAIAJCggAYIICAgCY8L2AotGoHMdRRUWF37sCAASIrwXU1NSkuro6XXXVVX7uBgAQQL4V0LFjx3T//ffrhRde0He+8x2/dgMACCjfCqi8vFwLFizQbbfdNuB6iURC8Xi81w0AMPL5ci24l19+WS0tLWpqahp03Wg0qqefftqPGACANOb5CKi9vV0rVqzQ1q1bNXr06EHXr6qqUiwWS97a29u9jgQASEOO6/Glc+vr6/WDH/xAGRkZyWXd3d1yHEejRo1SIpHo9di3xeNx366ELXGl4P5wNezU42rYGMnOvI/HYjHl5OScdT3PP4K79dZb9fHHH/da9tBDD2nKlClatWrVgOUDALhweF5A2dnZmjZtWq9l48aNU15eXp/lAIALF1dCAACYSMkvor7zzjup2A0AIEAYAQEATFBAAAATFBAAwAQFBAAwQQEBAEykZBbcuRjsG7Tnim+gjyxBfT793HZQjwkuPIyAAAAmKCAAgAkKCABgggICAJiggAAAJiggAIAJCggAYIICAgCYoIAAACYoIACACQoIAGCCAgIAmKCAAAAmKCAAgAkKCABgggICAJiggAAAJiggAIAJCggAYIICAgCYoIAAACYoIACAiUzrAKnmuq5v23Ycx7dt+5nbz21L/h4XP/F89hXUY4L0xAgIAGCCAgIAmKCAAAAmKCAAgAkKCABgggICAJjwpYC+/PJLPfDAA8rLy9PYsWN19dVXq7m52Y9dAQACyvPvAX399deaPXu25s2bp9dff13jx4/XP/7xD1188cVe7woAEGCeF9C6detUWFioF198Mbls4sSJXu8GABBwnn8Et3PnTpWUlOiee+7R+PHjdc011+iFF17wejcAgIDzvIA+//xzbdq0SZMnT9Ybb7yhZcuWafny5frDH/7Q7/qJRELxeLzXDQAw8jmuxxdgysrKUklJifbt25dctnz5cjU1Nendd9/ts351dbWefvrpPstjsZhycnK8jOY7rpPVv6BeC85PQX0++RvHUMTjcYXD4UHfxz0fAV166aW68sorey274oordOTIkX7Xr6qqUiwWS97a29u9jgQASEOeT0KYPXu2Pv30017LDh48qAkTJvS7figUUigU8joGACDNeT4Cevzxx/Xee++ppqZGhw4d0rZt21RXV6fy8nKvdwUACDDPC+i6667Tjh079NJLL2natGn6xS9+oQ0bNuj+++/3elcAgADzfBLC+Rrqyat0xAna/jEJoa+gPp/8jWMozCYhAAAwFBQQAMAEBQQAMEEBAQBMUEAAABOefxH1QubnLJ4gzz4K6nFBX0F9Lplhl54YAQEATFBAAAATFBAAwAQFBAAwQQEBAExQQAAAExQQAMAEBQQAMEEBAQBMUEAAABMUEADABAUEADBBAQEATFBAAAATFBAAwAQFBAAwQQEBAExQQAAAExQQAMAEBQQAMEEBAQBMUEAAABOZ1gEwNK7r+rZtx3F827bkb3Y/t+0nv4+5X4L6XAb5b3wkYwQEADBBAQEATFBAAAATFBAAwAQFBAAwQQEBAExQQAAAE54X0OnTp/Xkk0+quLhYY8aM0aRJk/TMM8+op6fH610BAALM8y+irlu3Ts8//7y2bNmiqVOnav/+/XrooYcUDoe1YsUKr3cHAAgozwvo3Xff1fe//30tWLBAkjRx4kS99NJL2r9/v9e7AgAEmOcfwc2ZM0dvvfWWDh48KEn68MMPtXfvXt1xxx39rp9IJBSPx3vdAAAjn+cjoFWrVikWi2nKlCnKyMhQd3e31q5dq/vuu6/f9aPRqJ5++mmvYwAA0pznI6Dt27dr69at2rZtm1paWrRlyxb9+te/1pYtW/pdv6qqSrFYLHlrb2/3OhIAIA15PgJauXKlVq9ercWLF0uSpk+frsOHDysajWrp0qV91g+FQgqFQl7HAACkOc9HQCdOnNCoUb03m5GRwTRsAEAvno+AFi5cqLVr16qoqEhTp07VgQMHtH79ej388MNe7woAEGCO6/EvKXV1dempp57Sjh071NnZqUgkovvuu08/+9nPlJWVNei/j8fjCofDisViysnJ8TIazoIf60o9fpAutfgbT62hvo97XkDniwJKPV6cqUcBpRZ/46k11PdxrgUHADBBAQEATFBAAAATFBAAwITn07ARPH6fQPXzBDAnf1MrqM8lf+PpiREQAMAEBQQAMEEBAQBMUEAAABMUEADABAUEADBBAQEATFBAAAATFBAAwAQFBAAwQQEBAExQQAAAExQQAMAEBQQAMEEBAQBMUEAAABMUEADABAUEADBBAQEATFBAAAATFBAAwAQFBAAwkWkdACOf67q+bdtxHN+27WfuoB4TPwX1ufR7+0E+LoNhBAQAMEEBAQBMUEAAABMUEADABAUEADBBAQEATFBAAAATwy6g3bt3a+HChYpEInIcR/X19b0ed11X1dXVikQiGjNmjObOnatPPvnEs8AAgJFh2AV0/PhxzZgxQ7W1tf0+/uyzz2r9+vWqra1VU1OTCgoKdPvtt6urq+u8wwIARo5hXwlh/vz5mj9/fr+Pua6rDRs2aM2aNbrrrrskSVu2bFF+fr62bdumRx555PzSAgBGDE/PAbW1tamjo0OlpaXJZaFQSLfccov27dvX779JJBKKx+O9bgCAkc/TAuro6JAk5efn91qen5+ffOzbotGowuFw8lZYWOhlJABAmvJlFty3L57nuu5ZL6hXVVWlWCyWvLW3t/sRCQCQZjy9GnZBQYGk/46ELr300uTyzs7OPqOiM0KhkEKhkJcxAAAB4OkIqLi4WAUFBWpoaEguO3nypBobGzVr1iwvdwUACLhhj4COHTumQ4cOJe+3tbWptbVVubm5KioqUkVFhWpqajR58mRNnjxZNTU1Gjt2rJYsWeJpcABAsA27gPbv36958+Yl71dWVkqSli5dqs2bN+snP/mJ/vOf/+jHP/6xvv76a91www168803lZ2d7V1qAEDgOa71T+J9SzweVzgcViwWU05OjnUcpLmR/GuR5yqov4jqp6A+l1Iw/8aH+j7OteAAACYoIACACQoIAGCCAgIAmPD0i6hAqvl5cjmIJ3/93jYTHFJvJD+fjIAAACYoIACACQoIAGCCAgIAmKCAAAAmKCAAgAkKCABgggICAJiggAAAJiggAIAJCggAYIICAgCYoIAAACYoIACACQoIAGCCAgIAmKCAAAAmKCAAgAkKCABgggICAJiggAAAJiggAICJTOsAQLpyXde3bTuO49u2/czt57b95Ofx9lsQn894PK5wODzoeoyAAAAmKCAAgAkKCABgggICAJiggAAAJiggAIAJCggAYGLYBbR7924tXLhQkUhEjuOovr4++dipU6e0atUqTZ8+XePGjVMkEtGDDz6oo0ePehoaABB8wy6g48ePa8aMGaqtre3z2IkTJ9TS0qKnnnpKLS0teuWVV3Tw4EHdeeednoQFAIwcjnseX4V1HEc7duzQokWLzrpOU1OTrr/+eh0+fFhFRUWDbvPMN2hjsZhycnLONRqQ1oJ6JYSg4koIqTXU93HfL8UTi8XkOI4uvvjifh9PJBJKJBLJ+/F43O9IAIA04OskhG+++UarV6/WkiVLztqC0WhU4XA4eSssLPQzEgAgTfhWQKdOndLixYvV09OjjRs3nnW9qqoqxWKx5K29vd2vSACANOLLR3CnTp3Svffeq7a2Nr399tsDfgYYCoUUCoX8iAEASGOeF9CZ8vnss8+0a9cu5eXleb0LAMAIMOwCOnbsmA4dOpS839bWptbWVuXm5ioSiejuu+9WS0uLXnvtNXV3d6ujo0OSlJubq6ysLO+SAwACbdjTsN955x3Nmzevz/KlS5equrpaxcXF/f67Xbt2ae7cuYNun2nYuBAwDTu1mIadWr5Nw547d+6ABySIBwsAkHpcCw4AYIICAgCYoIAAACYoIACACQoIAGDC94uRAujLz9miTPEeWUby88kICABgggICAJiggAAAJiggAIAJCggAYIICAgCYoIAAACYoIACACQoIAGCCAgIAmKCAAAAmKCAAgAkKCABgggICAJiggAAAJiggAIAJCggAYIICAgCYoIAAACYoIACACQoIAGCCAgIAmMi0DgDAW67r+rZtx3F827afuf3ctuTvcfGTdW5GQAAAExQQAMAEBQQAMEEBAQBMUEAAABMUEADAxLALaPfu3Vq4cKEikYgcx1F9ff1Z133kkUfkOI42bNhwXiEBACPPsAvo+PHjmjFjhmprawdcr76+Xu+//74ikcg5hwMAjFzD/iLq/PnzNX/+/AHX+fLLL/Xoo4/qjTfe0IIFC845HABg5PL8Sgg9PT0qKyvTypUrNXXq1EHXTyQSSiQSyfvxeNzrSACANOT5JIR169YpMzNTy5cvH9L60WhU4XA4eSssLPQ6EgAgDXlaQM3NzXruuee0efPmIV9jqKqqSrFYLHlrb2/3MhIAIE15WkB79uxRZ2enioqKlJmZqczMTB0+fFhPPPGEJk6c2O+/CYVCysnJ6XUDAIx8np4DKisr02233dZr2fe+9z2VlZXpoYce8nJXAICAG3YBHTt2TIcOHUreb2trU2trq3Jzc1VUVKS8vLxe61900UUqKCjQ5Zdffv5pAQAjxrALaP/+/Zo3b17yfmVlpSRp6dKl2rx5s2fBAAAjm+P6/UtNwxSPxxUOhxWLxTgfBKSZoP4gnd+sf9gtXQ32Ps614AAAJiggAIAJCggAYIICAgCYoIAAACY8vxgpgJHLz5lqQZ5hF9TjYo0READABAUEADBBAQEATFBAAAATFBAAwAQFBAAwQQEBAExQQAAAExQQAMAEBQQAMEEBAQBMUEAAABMUEADABAUEADBBAQEATFBAAAATFBAAwAQFBAAwQQEBAExQQAAAExQQAMBEpnWAb3NdV5IUj8eNkwAYKXg/sXHm/fxs0q6Aurq6JEmFhYXGSQCMFOFw2DrCBamrq2vAY++4g1VUivX09Ojo0aPKzs6W4ziDrh+Px1VYWKj29nbl5OSkIKE3yJ1aQc0tBTc7uVMrnXK7rquuri5FIhGNGnX2Mz1pNwIaNWqULrvssmH/u5ycHPODfi7InVpBzS0FNzu5Uytdcg9l1MkkBACACQoIAGAio7q6uto6xPnKyMjQ3LlzlZmZdp8oDojcqRXU3FJws5M7tYKWO+0mIQAALgx8BAcAMEEBAQBMUEAAABMUEADABAUEADBBAQEATFBAAAATFBAAwMT/A+FF/l6Ww3JBAAAAAElFTkSuQmCC",
      "text/plain": [
       "Figure(PyObject <Figure size 640x480 with 1 Axes>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "PyObject <matplotlib.image.AxesImage object at 0x13b070e10>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using LinearAlgebra, SparseArrays\n",
    "\n",
    "N=5\n",
    "xgrid=range(0,stop=1,length=N+1)\n",
    "Δx=1/N\n",
    "\n",
    "#assemling the finite different operator for the left hand side of the PDE\n",
    "T=spdiagm(0=>-2.0 .*ones(N-1), 1=>ones(N-2), -1=>ones(N-2)) \n",
    "Ii=spdiagm( 0=> ones(N-1))\n",
    "\n",
    "A=kron(T,Ii)+kron(Ii,T)\n",
    "\n",
    "using PyPlot\n",
    "spy(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAGiCAYAAACh/hJSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dfXBU1f3H8c9CkuXBZCU87LJDZDI1o2IAMVEkWKECcRgRHVtBRUunTkcU0BQYJfqH0WkTakd8GCot1BHQ2vhHjaVTEeKoQZphhAhjAEfpQEuQpKlO3A2IGwzn94fD/rpkg9zkLmcf3q+ZOyP33tx7zr2bfDy73z3XY4wxAgDgAhtguwEAgMxEAAEArCCAAABWEEAAACsIIACAFQQQAMAKAggAYAUBBACwggACAFhBAAEArEjqAHrxxRdVWFioQYMGqaSkRB988IHtJvXL9u3bdcsttygYDMrj8ejNN9+M2W6MUVVVlYLBoAYPHqzp06dr//79llrbNzU1NbrmmmuUm5urUaNG6bbbbtOnn34as08kEtHSpUs1YsQIDR06VHPnztXRo0cttbhv1q5dqwkTJigvL095eXmaMmWKtmzZEt2eDn08W01NjTwejyoqKqLr0qGfVVVV8ng8MUsgEIhuT4ffyzM+//xz3XPPPRo+fLiGDBmiq666Sk1NTdHtF7qvSRtAr7/+uioqKvT4449rz549+uEPf6jZs2fryJEjtpvWZydOnNDEiRO1Zs2auNuffvpprV69WmvWrNGuXbsUCAQ0a9YsdXZ2XuCW9l1DQ4MWL16snTt3qr6+Xt9++63Ky8t14sSJ6D4VFRWqq6tTbW2tduzYoePHj2vOnDnq7u622HJnxowZo1WrVmn37t3avXu3brzxRt16663RX9Z06OP/2rVrl9atW6cJEybErE+Xfl555ZVqbW2NLs3NzdFt6fB7KUkdHR2aOnWqsrOztWXLFh04cEDPPPOMLr744ug+F7yvJklde+21ZtGiRTHrLr/8crNy5UpLLXKXJFNXVxf99+nTp00gEDCrVq2Krvvmm2+Mz+czv//972000RXt7e1GkmloaDDGGPPVV1+Z7OxsU1tbG93n888/NwMGDDBvv/22rWa6YtiwYeaPf/xj2vWxs7PTFBUVmfr6ejNt2jTz8MMPG2PS514+8cQTZuLEiXG3pdPv5aOPPmquv/76Xrfb6GtSjoC6urrU1NSk8vLymPXl5eVqbGy01KrEOnz4sNra2mL67PV6NW3atJTucygUkiTl5+dLkpqamnTq1KmYfgaDQRUXF6dsP7u7u1VbW6sTJ05oypQpadfHxYsX6+abb9bMmTNj1qdTPw8ePKhgMKjCwkLdeeedOnTokKT0+r3cvHmzSktLdccdd2jUqFGaNGmS1q9fH91uo69JGUBffPGFuru75ff7Y9b7/X61tbVZalVinelXOvXZGKNly5bp+uuvV3FxsaTv+pmTk6Nhw4bF7JuK/WxubtZFF10kr9erRYsWqa6uTuPGjUurPtbW1uqjjz5STU1Nj23p0s/Jkydr06ZN2rp1q9avX6+2tjaVlZXpyy+/TKvfy0OHDmnt2rUqKirS1q1btWjRIj300EPatGmTJDt/g7ISclSXeDyemH8bY3qsSzfp1OclS5bo448/1o4dO75331Ts52WXXaa9e/fqq6++0l/+8hctXLhQDQ0Nve6fan1saWnRww8/rG3btmnQoEHn/XOp1s/Zs2dH/3v8+PGaMmWKfvCDH2jjxo267rrrJKXH7+Xp06dVWlqq6upqSdKkSZO0f/9+rV27Vj/96U+j+13IviblCGjEiBEaOHBgj9Rtb2/vkc7p4kzVTbr0eenSpdq8ebPee+89jRkzJro+EAioq6tLHR0dMfunYj9zcnJ06aWXqrS0VDU1NZo4caKef/75tOljU1OT2tvbVVJSoqysLGVlZamhoUEvvPCCsrKy5Pf706KfZxs6dKjGjx+vgwcPptXv5ejRozVu3LiYdVdccUW0sMtGX5MygHJyclRSUqL6+vqY9fX19SorK7PUqsQqLCxUIBCI6XNXV5caGhpSqs/GGC1ZskRvvPGG3n33XRUWFsZsLykpUXZ2dkw/W1tbtW/fvpTqZzzGGEUikbTp44wZM9Tc3Ky9e/dGl9LSUi1YsCD63+nQz7NFIhF98sknGj16dNr8XkrS1KlTe3wl4rPPPtPYsWMlWfoblJDSBhfU1taa7Oxs89JLL5kDBw6YiooKM3ToUPOvf/3LdtP6rLOz0+zZs8fs2bPHSDKrV682e/bsMf/+97+NMcasWrXK+Hw+88Ybb5jm5mZz1113mdGjR5twOGy55efvgQceMD6fz7z//vumtbU1unz99dfRfRYtWmTGjBlj3nnnHfPRRx+ZG2+80UycONF8++23FlvuTGVlpdm+fbs5fPiw+fjjj81jjz1mBgwYYLZt22aMSY8+xvO/VXDGpEc/ly9fbt5//31z6NAhs3PnTjNnzhyTm5sb/VuTDr+Xxhjz4YcfmqysLPPrX//aHDx40PzpT38yQ4YMMa+++mp0nwvd16QNIGOM+d3vfmfGjh1rcnJyzNVXXx0t5U1V7733npHUY1m4cKEx5rsyyCeeeMIEAgHj9XrNDTfcYJqbm+022qF4/ZNkXn755eg+J0+eNEuWLDH5+flm8ODBZs6cOebIkSP2Gt0HP//5z6OvzZEjR5oZM2ZEw8eY9OhjPGcHUDr0c/78+Wb06NEmOzvbBINBc/vtt5v9+/dHt6fD7+UZf/vb30xxcbHxer3m8ssvN+vWrYvZfqH76jHGmMSMrQAA6F1SfgYEAEh/BBAAwAoCCABgBQEEALCCAAIAWEEAAQCsIIAAAFYkdQBFIhFVVVUpEonYbkpC0c/0kQl9lOhnurHVz6T+Imo4HJbP51MoFFJeXp7t5iQM/UwfmdBHiX6mG1v9TOoREAAgfSUsgF588UUVFhZq0KBBKikp0QcffJCoUwEAUlBCHkj3+uuvq6KiQi+++KKmTp2qP/zhD5o9e7YOHDigSy655Jw/e/r0aR07dky5ubnq7OyU9N3wMJ2d6R/9TH2Z0EeJfqYbt/tpjFFnZ6eCwaAGDDjHOCcRM5xee+21ZtGiRTHrLr/8crNy5crv/dmWlpZeZ1RmYWFhYUmdpaWl5Zx/710fAXV1dampqUkrV66MWV9eXq7GxsYe+0cikZjKC9OHmohQKHTe+/p8vn4fI91wTQC4KRwOq6CgQLm5uefcz/UA+uKLL9Td3d3jEa5+v7/Ho14lqaamRk8++WS/zulG1UY6V7j0FdcEQH94PJ5zbk9YEcLZJzbGxG1MZWWlQqFQdGlpaUlUkwAAScT1EdCIESM0cODAHqOd9vb2HqMiSfJ6vfJ6vW43AwCQ5FwfAeXk5KikpET19fUx6+vr61VWVub26SR9N9o6e+mN+e4x5D2WeMf4vuFjuuCaALAhIWXYy5Yt07333qvS0lJNmTJF69at05EjR7Ro0aJEnA4AkIISEkDz58/Xl19+qaeeekqtra0qLi7WW2+9pbFjxybidACAFJR0c8GdmZOov5x2q7e3lpLs8lxQXBMAfXG+c8sxFxwAwAoCCABgRUI+A3JDvKGbkwosp28f9bY+3nEy5S0oJ9fkXPsDQDyMgAAAVhBAAAArCCAAgBUEEADAiqQtQojH6YfiiTpnpn8IT3ECADcwAgIAWEEAAQCsIIAAAFYQQAAAKwggAIAVSVsFF29GbKfT6MTjtGIu3rGpAouP6YwAOMEICABgBQEEALCCAAIAWEEAAQCsIIAAAFYkbRWckwfSJbKiykkFF9Vx8TGfHoB4GAEBAKwggAAAVhBAAAArCCAAgBUEEADAiqStgovHjbnG3HiqqtMKLqrjeuKaAGAEBACwggACAFhBAAEArCCAAABWpFQRQm/cmOrFjeIEp3iAW08UJwCZgxEQAMAKAggAYAUBBACwggACAFhBAAEArEiLKrh43KqmclJl5bRizsk5M70KjOo4IP0wAgIAWEEAAQCsIIAAAFYQQAAAKwggAIAVaVsF15tkqqZycs5kancyYT49IHUxAgIAWEEAAQCsIIAAAFYQQAAAKxwH0Pbt23XLLbcoGAzK4/HozTffjNlujFFVVZWCwaAGDx6s6dOna//+/a41GACQHhwH0IkTJzRx4kStWbMm7vann35aq1ev1po1a7Rr1y4FAgHNmjVLnZ2d/W5sIhlj4i4ej6fH4vQYTrlxznjHSOTTXZMN1wRIfo7LsGfPnq3Zs2fH3WaM0XPPPafHH39ct99+uyRp48aN8vv9eu2113T//ff3r7UAgLTh6mdAhw8fVltbm8rLy6PrvF6vpk2bpsbGxrg/E4lEFA6HYxYAQPpzNYDa2tokSX6/P2a93++PbjtbTU2NfD5fdCkoKHCzSQCAJJWQKriz31c/8/57PJWVlQqFQtGlpaUlEU0CACQZV6fiCQQCkr4bCY0ePTq6vr29vceo6Ayv1yuv1+tmMwAAKcDVEVBhYaECgYDq6+uj67q6utTQ0KCysjI3T3XBuFFN5VZ1XH/b7bSqL91QMQgkF8cjoOPHj+uf//xn9N+HDx/W3r17lZ+fr0suuUQVFRWqrq5WUVGRioqKVF1drSFDhujuu+92teEAgNTmOIB2796tH/3oR9F/L1u2TJK0cOFCbdiwQY888ohOnjypBx98UB0dHZo8ebK2bdum3Nxc91oNAEh5HpNkc9SHw2H5fD6FQiHl5eXZbk5cbj0CwcnbPG7dJh5T0FOmP9ICcNv5/h1nLjgAgBUZ90A6N7j1cDgn/4ft9ENxN86ZKXjYH2AHIyAAgBUEEADACgIIAGAFAQQAsIIAAgBYQRWci5KpmooKrv5zcj+5roBzjIAAAFYQQAAAKwggAIAVBBAAwAoCCABgBVVwF4Ab1XFOj9EbKrj6L971ouoQcI4READACgIIAGAFAQQAsIIAAgBYQQABAKygCs4iN+Yac6s6Dv2TTPMAAqmCERAAwAoCCABgBQEEALCCAAIAWJG0RQg+n6/Hukz54NaNfjo5htOCBafty+TpfyhOAHrHCAgAYAUBBACwggACAFhBAAEArCCAAABWJG0VXCgUUl5eXsw6KoeSg9P7wAPceqI6DmAEBACwhAACAFhBAAEArCCAAABWEEAAACuStgouHjce4Iae3HqonZMKLqrA4uM1jkzCCAgAYAUBBACwggACAFhBAAEArCCAAABWpFQVXG+Yaywx3KqOS+Q5M+V+8hpHOmIEBACwggACAFhBAAEArCCAAABWOAqgmpoaXXPNNcrNzdWoUaN022236dNPP43ZJxKJaOnSpRoxYoSGDh2quXPn6ujRo642GgCQ+hwFUENDgxYvXqydO3eqvr5e3377rcrLy3XixInoPhUVFaqrq1Ntba127Nih48ePa86cOeru7na98edijIm7eDyeuEsmc3pNeru2Tq+5G+fM5PvJNUGq85h+1Gz+97//1ahRo9TQ0KAbbrhBoVBII0eO1CuvvKL58+dLko4dO6aCggK99dZbuummm773mOFwWD6fL+4jud1A6WpPib4mTv74OT0n97MnrglsO9+/4/36DCgUCkmS8vPzJUlNTU06deqUysvLo/sEg0EVFxersbEx7jEikYjC4XDMAgBIf30OIGOMli1bpuuvv17FxcWSpLa2NuXk5GjYsGEx+/r9frW1tcU9Tk1NjXw+X3QpKCjoa5MAACmkzwG0ZMkSffzxx/rzn//8vfueeV86nsrKSoVCoejS0tLS1yYBAFJIn6biWbp0qTZv3qzt27drzJgx0fWBQEBdXV3q6OiIGQW1t7errKws7rG8Xq+8Xm9fmtEnTPXSUzJdE6fn5AFuPSXT/QTOxdEIyBijJUuW6I033tC7776rwsLCmO0lJSXKzs5WfX19dF1ra6v27dvXawABADKToxHQ4sWL9dprr+mvf/2rcnNzo5/r+Hw+DR48WD6fT/fdd5+WL1+u4cOHKz8/XytWrND48eM1c+bMhHQAAJCaHAXQ2rVrJUnTp0+PWf/yyy/rZz/7mSTp2WefVVZWlubNm6eTJ09qxowZ2rBhgwYOHOhKgwEA6aFf3wNKhER/D6g3vD/ek1vXxI0vQbpxzky+lxKvcVw4F+R7QAAA9FVaPJDODVRT9eRWNZWTh6n1JpHn5H5m7mscdjECAgBYQQABAKwggAAAVhBAAAArCCAAgBVUwX0Pqql6cqM6zukx3MAcafHxGoctjIAAAFYQQAAAKwggAIAVBBAAwAoCCABgBVVwfUA1VXxuzDWWyJm23Xiq6rn2TydcE1wIjIAAAFYQQAAAKwggAIAVBBAAwAoCCABgBVVwLqJyKL5kmWvMjaeq9nacTL6XEq9x9A0jIACAFQQQAMAKAggAYAUBBACwgiKEC4APbnty65q48WA7N86ZyfdS4jWOvmEEBACwggACAFhBAAEArCCAAABWEEAAACuogrOIqV56SmR1nJPKOKeoAouP1zjOhREQAMAKAggAYAUBBACwggACAFhBAAEArKAKLgkx11hPblSZOb1WTqrmnM5Vl+l4jUNiBAQAsIQAAgBYQQABAKwggAAAVhBAAAArqIJLEcw1Fl+yzDWW6ffBDbzGMw8jIACAFQQQAMAKAggAYIWjAFq7dq0mTJigvLw85eXlacqUKdqyZUt0eyQS0dKlSzVixAgNHTpUc+fO1dGjR11vNAAg9TkKoDFjxmjVqlXavXu3du/erRtvvFG33nqr9u/fL0mqqKhQXV2damtrtWPHDh0/flxz5sxRd3d3QhqP7z6Ijbd4PJ64S6Zw45r0dm2dyPT74AZe4+nLY/pZSpKfn6/f/va3+slPfqKRI0fqlVde0fz58yVJx44dU0FBgd566y3ddNNN53W8cDgsn8+nUCikvLy8/jQto1E51JNb18SNP3KZfB/cwms8eZ3v3/E+fwbU3d2t2tpanThxQlOmTFFTU5NOnTql8vLy6D7BYFDFxcVqbGzs62kAAGnK8feAmpubNWXKFH3zzTe66KKLVFdXp3Hjxmnv3r3KycnRsGHDYvb3+/1qa2vr9XiRSESRSCT673A47LRJAIAU5HgEdNlll2nv3r3auXOnHnjgAS1cuFAHDhzodf8z79X2pqamRj6fL7oUFBQ4bRIAIAU5DqCcnBxdeumlKi0tVU1NjSZOnKjnn39egUBAXV1d6ujoiNm/vb1dfr+/1+NVVlYqFApFl5aWFue9AACknH5/D8gYo0gkopKSEmVnZ6u+vj66rbW1Vfv27VNZWVmvP+/1eqNl3WcW9J+TyqFM4VY1VX8r4+AOXuOpz9FnQI899phmz56tgoICdXZ2qra2Vu+//77efvtt+Xw+3XfffVq+fLmGDx+u/Px8rVixQuPHj9fMmTMT1X4AQIpyFED/+c9/dO+996q1tVU+n08TJkzQ22+/rVmzZkmSnn32WWVlZWnevHk6efKkZsyYoQ0bNmjgwIEJaTwAIHX1+3tAbuN7QIl1oWeJTgVufJ/E6ds8mX7NE4nXuH0J/x4QAAD9QQABAKzggXQZJt5bEZk+pYkb/Uz0dD5uvB2Yyfcz069JsmIEBACwggACAFhBAAEArCCAAABWEEAAACuogkOvlUBUDiUPJ1+u5H72xDVJToyAAABWEEAAACsIIACAFQQQAMAKAggAYAVVcOgVlUOJ4/TaOtnXjeq4TLmXvMbtYgQEALCCAAIAWEEAAQCsIIAAAFZQhADH+OA2cdwoTnDjnJl+L3mNXxiMgAAAVhBAAAArCCAAgBUEEADACgIIAGAFVXBwDVO9nD+3ptFxcuzexDs2VWDx8Rp3FyMgAIAVBBAAwAoCCABgBQEEALCCAAIAWEEVHBKOucZ6SqYqMycVXMnU7mTCa7xvGAEBAKwggAAAVhBAAAArCCAAgBUEEADACqrgYAXVVPG5MdeYG09VdWuuuky+n1yT78cICABgBQEEALCCAAIAWEEAAQCsoAgBSYUPbuNzY6oXN4oTnOIBbj3xGv9/jIAAAFYQQAAAKwggAIAVBBAAwIp+BVBNTY08Ho8qKiqi6yKRiJYuXaoRI0Zo6NChmjt3ro4ePdrvhgIA0kufA2jXrl1at26dJkyYELO+oqJCdXV1qq2t1Y4dO3T8+HHNmTNH3d3d/W4sMpcxJu7i8XjiLpnArWvS23GcHNuNc2byvZQy8zXepwA6fvy4FixYoPXr12vYsGHR9aFQSC+99JKeeeYZzZw5U5MmTdKrr76q5uZmvfPOO641GgCQ+voUQIsXL9bNN9+smTNnxqxvamrSqVOnVF5eHl0XDAZVXFysxsbGuMeKRCIKh8MxCwAg/Tn+Imptba0++ugj7dq1q8e2trY25eTkxIyKJMnv96utrS3u8WpqavTkk086bQYAIMU5GgG1tLTo4Ycf1quvvqpBgwad98+deR8znsrKSoVCoejS0tLipEkAgBTlKICamprU3t6ukpISZWVlKSsrSw0NDXrhhReUlZUlv9+vrq4udXR0xPxce3u7/H5/3GN6vV7l5eXFLACA9OcogGbMmKHm5mbt3bs3upSWlmrBggXR/87OzlZ9fX30Z1pbW7Vv3z6VlZW53njASeVQpkimaion50ymdieTdH6NO/oMKDc3V8XFxTHrhg4dquHDh0fX33fffVq+fLmGDx+u/Px8rVixQuPHj+9RsAAAyGyuz4b97LPPKisrS/PmzdPJkyc1Y8YMbdiwQQMHDnT7VACAFOYxSTbXdzgcls/nUygU4vMg9FkmT/ffGzem+3frbR43zsn9TN7X+Pn+HWcuOACAFQQQAMAKnoiKtOTGE0TTjRtPJ3XrqaqJPGcm389UuyaMgAAAVhBAAAArCCAAgBUEEADACgIIAGAFVXDIGFRTxedGNZVb1XFOuFHVl25S7TXOCAgAYAUBBACwggACAFhBAAEArKAIARkv1T64vRDcuiaJnGnbjXNmimR9jTMCAgBYQQABAKwggAAAVhBAAAArCCAAgBVUwQG9YKqXnpKpmsp2BVc6sP0aZwQEALCCAAIAWEEAAQCsIIAAAFYQQAAAK6iCAxxy4wFu6caN6ji3HmqXyVWKbrlQr3FGQAAAKwggAIAVBBAAwAoCCABgBQEEALCCKjjABck0R1oycWOuMbeq49A/ibgPjIAAAFYQQAAAKwggAIAVBBAAwAqKEIAEojghPjf66eQYTj8od9q+TJ7+J14/w+GwfD7f9/4sIyAAgBUEEADACgIIAGAFAQQAsIIAAgBYQRUcYAHVccnN6X3gIYV9wwgIAGAFAQQAsIIAAgBYQQABAKxwFEBVVVXyeDwxSyAQiG43xqiqqkrBYFCDBw/W9OnTtX//ftcbDQBIfY5HQFdeeaVaW1ujS3Nzc3Tb008/rdWrV2vNmjXatWuXAoGAZs2apc7OTlcbDaQrY0zc5ez/8eNhbM70dl2dincfzlXtdr73MlPvp+MAysrKUiAQiC4jR46U9N3Ffu655/T444/r9ttvV3FxsTZu3Kivv/5ar732musNBwCkNscBdPDgQQWDQRUWFurOO+/UoUOHJEmHDx9WW1ubysvLo/t6vV5NmzZNjY2NvR4vEokoHA7HLACA9OcogCZPnqxNmzZp69atWr9+vdra2lRWVqYvv/xSbW1tkiS/3x/zM36/P7otnpqaGvl8vuhSUFDQh24AAFKNowCaPXu2fvzjH2v8+PGaOXOm/v73v0uSNm7cGN3n7Pcyz7zn2ZvKykqFQqHo0tLS4qRJAIAU1a8y7KFDh2r8+PE6ePBgtBru7NFOe3t7j1HR//J6vcrLy4tZAADpr18BFIlE9Mknn2j06NEqLCxUIBBQfX19dHtXV5caGhpUVlbW74YCmYxqqsRwqzrOjXNm4v10NBnpihUrdMstt+iSSy5Re3u7fvWrXykcDmvhwoXyeDyqqKhQdXW1ioqKVFRUpOrqag0ZMkR33313otoPAEhRjgLo6NGjuuuuu/TFF19o5MiRuu6667Rz506NHTtWkvTII4/o5MmTevDBB9XR0aHJkydr27Ztys3NTUjjAQCpy2OSbG7wcDgsn8+nUCjE50HAOTDdf+I4fevLjWueTvfzfP+OMxccAMAKAggAYAVPRAVSFE9VPX9uPOG0L8d345zpfD8ZAQEArCCAAABWEEAAACsIIACAFRQhAGkmEz/M/j7JdE3cKoiId5xUu5eMgAAAVhBAAAArCCAAgBUEEADACgIIAGAFVXBAhkjnaqq+cqs6Lt56pzNqJ/KcyXo/GQEBAKwggAAAVhBAAAArCCAAgBUEEADACqrggAyXDtVUbnOjOs7pMdyQTHPenQ9GQAAAKwggAIAVBBAAwAoCCABgBQEEALCCKjgAPaRaNdWF4sZ8ek6vlZOqOTeeqnqu/d3GCAgAYAUBBACwggACAFhBAAEArCCAAABWUAUH4LwlazWVbckyn54bT1Xt7TiJaDcjIACAFQQQAMAKAggAYAUBBACwgiIEAP1GcUJPbl0TNx5s58Y5E3EvGQEBAKwggAAAVhBAAAArCCAAgBUEEADACqrgACSM7aleklEiq+OcVMY55UY13tkYAQEArCCAAABWEEAAACsIIACAFY4D6PPPP9c999yj4cOHa8iQIbrqqqvU1NQU3W6MUVVVlYLBoAYPHqzp06dr//79rjYaAJD6HAVQR0eHpk6dquzsbG3ZskUHDhzQM888o4svvji6z9NPP63Vq1drzZo12rVrlwKBgGbNmqXOzk7XGw8gNRljeiwejyfukiniXROn16W3Yzg9tpP7EO+4oVDovPrsMQ5qH1euXKl//OMf+uCDD3ptSDAYVEVFhR599FFJUiQSkd/v129+8xvdf//933uOcDgsn8+nUCikvLy8820agBSXyROXnksir4uTgHdyvvP9O+5oBLR582aVlpbqjjvu0KhRozRp0iStX78+uv3w4cNqa2tTeXl5dJ3X69W0adPU2NgY95iRSEThcDhmAQCkP0cBdOjQIa1du1ZFRUXaunWrFi1apIceekibNm2SJLW1tUmS/H5/zM/5/f7otrPV1NTI5/NFl4KCgr70AwCQYhwF0OnTp3X11QLuSIoAAAKMSURBVFerurpakyZN0v33369f/OIXWrt2bcx+Zw/rzrzXGE9lZaVCoVB0aWlpcdgFAEAqchRAo0eP1rhx42LWXXHFFTpy5IgkKRAISFKP0U57e3uPUdEZXq9XeXl5MQsAIP05CqCpU6fq008/jVn32WefaezYsZKkwsJCBQIB1dfXR7d3dXWpoaFBZWVlLjQXQLpyowosHTm5LomUiPvgaDLSX/7ylyorK1N1dbXmzZunDz/8UOvWrdO6deuiDayoqFB1dbWKiopUVFSk6upqDRkyRHfffXe/GgoASC+OAuiaa65RXV2dKisr9dRTT6mwsFDPPfecFixYEN3nkUce0cmTJ/Xggw+qo6NDkydP1rZt25Sbm+t64wEAqcvR94AuBL4HBOB/8f2g+Nx4pIUbb9vFO2dCvgcEAIBbeCAdgKTm1gPc0o2TB9I5eajduY7T333PxggIAGAFAQQAsIIAAgBYQQABAKwggAAAVlAFByAlOangyuTKOMmd6rhETPXDCAgAYAUBBACwggACAFhBAAEArEi6IoQzH36Fw2HLLQGQLvh7El+ir8v3FX8kXQB1dnZKkgoKCiy3BEC68Pl8tpuQlBJ9XTo7O895jqR7HMPp06d17Ngx5ebmqrOzUwUFBWppaUnrRzOEw2H6mSYyoY8S/Uw3bvfTGKPOzk4Fg0ENGND7Jz1JNwIaMGCAxowZI+n/687z8vLS+uafQT/TRyb0UaKf6cbNfp7P6IoiBACAFQQQAMCKgVVVVVW2G3EuAwcO1PTp05WVlXTvFrqKfqaPTOijRD/TjY1+Jl0RAgAgM/AWHADACgIIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAArCCAAgBX/B9dfAMowZNxOAAAAAElFTkSuQmCC",
      "text/plain": [
       "Figure(PyObject <Figure size 640x480 with 1 Axes>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "PyObject <matplotlib.image.AxesImage object at 0x13b899dd8>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A=kron(kron(T,Ii),Ii)+ kron(Ii,kron(Ii,T))+ kron(kron(Ii,T),Ii)\n",
    "\n",
    "size(A)\n",
    "\n",
    "using PyPlot\n",
    "spy(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "size of the matrix:(13824, 13824)\n",
      "Optimal direct solver: \n",
      "  0.548099 seconds (74 allocations: 55.256 MiB, 0.67% gc time)\n",
      "LU:\n",
      "  0.276399 seconds (70 allocations: 88.185 MiB, 31.96% gc time)\n",
      "  0.364837 seconds (82 allocations: 173.523 MiB, 31.28% gc time)\n",
      "  0.142359 seconds (42 allocations: 86.761 MiB, 20.14% gc time)\n",
      "\n",
      "\n",
      "Gauss-Seidel finishes in 172 iterations.\n",
      "302.226181 seconds (4.76 M allocations: 245.111 GiB, 11.84% gc time)\n",
      "\n",
      "\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "N=25\n",
    "xgrid=range(0,stop=1,length=N+1)\n",
    "Δx=1/N\n",
    "\n",
    "f=(x,y)-> x.^2+y+z.^2\n",
    "\n",
    "#assemling the finite different operator for the left hand side of the PDE\n",
    "T=spdiagm(0=>-2.0 .*ones(N-1), 1=>ones(N-2), -1=>ones(N-2)) \n",
    "Ii=spdiagm( 0=> ones(N-1))\n",
    "\n",
    "#discretize the Laplace operator in 2D\n",
    "#A=kron(T,I)+kron(I,T)\n",
    "\n",
    "#discretize the Laplace operator in 3D\n",
    "A=kron(kron(T,Ii),Ii)+ kron(Ii,kron(Ii,T))+ kron(kron(Ii,T),Ii)\n",
    "\n",
    "F=ones(size(A,1))\n",
    "\n",
    "b=F.*Δx^2\n",
    "\n",
    "U=copy(F)\n",
    "Uj=copy(F)\n",
    "Ugs=copy(F)\n",
    "Usor=copy(F)\n",
    "\n",
    "println(\"size of the matrix:\", size(A))\n",
    "\n",
    "println(\"Optimal direct solver: \")\n",
    "@time Ud=A\\b\n",
    "\n",
    "println(\"LU:\")\n",
    "@time fac=lu(A)\n",
    "#Step 1:\n",
    "@time y=fac.L\\b[fac.p]\n",
    "#Step 2:\n",
    "@time x=fac.U\\y\n",
    "\n",
    "#@time U=inv(fac.U)*inv(fac.L)*b[fac.p]\n",
    "\n",
    "println(\"\\n\")\n",
    "@time GaussSeidelIterations!(Ugs, A, b, 1e-1)\n",
    "\n",
    "println(\"\\n\")\n",
    "#@time JacobiIterations!(Uj, A, b, 1e-12)\n",
    "\n",
    "println(\"\\n\")\n",
    "#@time SOR!(Usor, A, b, 1.5, 1e-12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Convergence of iterative algorithms\n",
    "\n",
    "Here we discuss the convergence issue for linear systems of equations and *linear* methods:\n",
    "\n",
    "Solve $Ax=b$ by the fixed point iterations defined below:\n",
    "\n",
    "Fixed Point Method: $$ x^{(k+1)} = Gx^{(k)}+c $$ where $G$ is a matrix and $x^{(k)}$ and $c$ are vectors.\n",
    "\n",
    "Solution vector $x^*$ is a fixed point for $Gx+c$: $$x^* = Gx^*+c$$\n",
    "\n",
    "Derive an equation for the errors:\n",
    "\n",
    "$$x^{(k+1)}-x^* =G(x^{(k)}-x^*)$$ \n",
    "\n",
    "* linearly convergent algorithms: (**the error in the next step is a linear function of the current error**)\n",
    "$$ e_{n+1} = a e_n $$\n",
    "\n",
    "### What you need to know\n",
    "The number $a$ above is related to the Eigenvalues of $G$.\n",
    "\n",
    "Consider the following simplied situation, where $x^{(0)}-x^*$ happens to be an Eigenvector of $G$, with the associated Eigenvalue denoted by $\\lambda$.\n",
    "\n",
    "Then, we notice that\n",
    "* $e^{(1)} =\\lambda e^{(0)}$ ... \n",
    "* $e^{(k)} =\\lambda^k e^{(0)}$\n",
    "\n",
    "So if $|\\lambda|<1$, $e^{(k)} \\rightarrow 0$ as $k\\rightarrow \\infty$!\n",
    "\n",
    "The next step is to assume that $G$ has a complete set of Eigenvalues so we can use the Eigenvectors to form a basis.\n",
    "$$Gv_j=\\lambda_j v_j$$\n",
    "\n",
    "Write the initial error vector using the new basis (as a linear combination of $v_j$'s):\n",
    "$$ x^{(0)}-x^* = c_1 v_1+c_2v_2+\\cdots+c_n v_n$$\n",
    "\n",
    "Then we have..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The conjugate gradient (CG) method of Hestenes and Stiefel (1952)\n",
    "\n",
    "$$Ax=b$$\n",
    "* It is optimal for symmetric positive definite systems: $A^T =A,~~~x^TAx>0,~\\forall x\\neq 0$\n",
    "* It is a type of **Krylov methods**\n",
    "* In each iteration, say iteration $k$, the algorithm chooses the next guess, $x^{(k+1)}$, by \n",
    "  * minimizes the residual, $r^{(k+1)}:=b-Ax^{(k+1)}$, \n",
    "  * in a special norm built from $A$ and $x^{(j)},j=1,2,\\cdots,k$.\n",
    "* Very simple, elegant and fast algorithm.\n",
    "\n",
    "* We have seen that solving teh model Poisson's problem in 1 and 2 dimenions, the backlash operator is faster than CG. But the situation changes when we solve the Poisson's equation in 3D, for which the matrix has more strides away from teh diagonal.\n",
    "\n",
    "* See Algorithm 11.15 in GGK or this workbook (instead of BF) for a more comprehensive view of teh CG algorithm.\n",
    "\n",
    "**Learning objectives**:\n",
    "\n",
    "* Know when CG is applicable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "CG (generic function with 2 methods)"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CG(A,b;kwargs...) = CG(x -> A*x,b;kwargs...)\n",
    "    \n",
    "function CG(A::Function, b; tol=1e-2, maxItr=1000)\n",
    "\n",
    "    n = length(b)\n",
    "\n",
    "    #initial guess\n",
    "    x=zeros(n)\n",
    "    \n",
    "    #initial residual\n",
    "    r=b-A(x)\n",
    "\n",
    "    p = copy(r)\n",
    "    \n",
    "    res0  = norm(r)\n",
    "\n",
    "    itr   = 1 \n",
    "    \n",
    "    for itr=1:maxItr\n",
    "        Ap = A(p)\n",
    "        pAp= dot(p,Ap)\n",
    "        α  = dot(r,p)/pAp\n",
    "\n",
    "        #update the approximate solution\n",
    "        x.+=α.*p\n",
    "        \n",
    "        #update the residual vector\n",
    "        r.-=α.*Ap\n",
    "            \n",
    "        #break the loop when tolerance is met\n",
    "        if norm(r)/res0 <= tol\n",
    "            break\n",
    "        end\n",
    "\n",
    "        #update search direction\n",
    "        β = -dot(Ap,r)/pAp           \n",
    "        p.=r.+β.*p\n",
    "    end\n",
    "\n",
    "    println(\"cg terminates after \", itr, \" steps, with residual=\", norm(r))\n",
    "    \n",
    "    x\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cg terminates after 1 steps, with residual=1.7008927367038972e-15\n",
      "  0.020896 seconds (353 allocations: 9.081 MiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "3.609376587170471e-14"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@time uu=CG(A,b, tol=1e-14)\n",
    "\n",
    "norm(Ud-uu)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "using KrylovMethods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cg iterated maxIter (=1500) times but reached only residual norm 1.05e+00 instead of tol=1.00e-10.\n",
      "  0.392077 seconds (877.73 k allocations: 44.803 MiB, 5.25% gc time)\n",
      "cg iterated maxIter (=1500) times but reached only residual norm 1.05e+00 instead of tol=1.00e-10.\n",
      "  0.000845 seconds (14.55 k allocations: 393.406 KiB)\n",
      "cg iterated maxIter (=1500) times but reached only residual norm 1.05e+00 instead of tol=1.00e-14.\n",
      "  0.000858 seconds (14.55 k allocations: 394.094 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(3, 3)"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@time uu=KrylovMethods.cg(A, b, tol=1e-10, maxIter=1500, out=1);\n",
    "\n",
    "@time uu=KrylovMethods.cg(A, b, tol=1e-10, maxIter=1500, out=1);\n",
    "\n",
    "@time uu=KrylovMethods.cg(A, b, tol=1e-14, maxIter=1500, out=1);\n",
    "\n",
    "size(A)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Julia 1.2.0",
   "language": "julia",
   "name": "julia-1.2"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}