lectures/lecture9/lecture-9.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 9: From dense to sparse linear algebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Recap of the previous lectures\n",
    "- Algorithms for the symmetric eigenvalue (SEV) problem\n",
    "    - QR algorithm with shifts\n",
    "    - Divide-and-Conquer\n",
    "    - Bisection\n",
    "    - Jacobi method\n",
    "- Reduction problem of computing SVD to SEV problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Large scale dense matrices\n",
    "\n",
    "- If the size of the dense matrix is huge, then it can not be stored in memory\n",
    "- Possible options\n",
    "    - This matrix is **structured**, e.g. block Toeplitz with Toeptitz blocks (next lectures). Then the compressed storage is possible\n",
    "    - For unstructured dense matrices **distributed memory** helps\n",
    "    - MPI for processing distributed storing matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Distributed memory and MPI\n",
    "\n",
    "- Split matrix into blocks and store them on different machines\n",
    "- Every machine has its own address space and can not damage data on other machines\n",
    "- In this case machines communicate with each other to aggregate results of computations\n",
    "- [MPI (Message Passing Interface)](https://en.wikipedia.org/wiki/Message_Passing_Interface) is a standard for parallel computing in distributed memory\n",
    "\n",
    "- Skoltech owns supercomputer, which is on the 7th place in the list of TOP50 supercomputers from Commonwealth of Independent States (see the rating [here](http://top50.supercomputers.ru/?page=rating), sorry for only in Russian)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Example: matrix-by-vector product\n",
    "\n",
    "- Assume you want to compute $Ax$ and matrix $A$ can not be stored in available memory\n",
    "- Then you can split it on blocks and distribute blocks on separate machines\n",
    "- Possible strategies\n",
    "    - 1D blocking splits only rows on blocks\n",
    "    - 2D blocking splits both rows and columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### 1D blocking scheme\n",
    "<img src=\"1d_block.jpg\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Total time of computing matvec with 1D blocking\n",
    "\n",
    "- Each machine has $n / p $ complete rows and $n / p$ elements of vector\n",
    "- Total operations are $n^2 / p$\n",
    "- Total time for sending and writing data are $t_s \\log p + t_w n$, where $t_s$ time unit for sending and $t_w$ time unit for writing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### 2D blocking scheme\n",
    "\n",
    "<img src=\"2d_block.png\" width=400>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Total time of computing matvec with 2D blocking\n",
    "\n",
    "- Each machine has $n / \\sqrt{p} $ size block and $n / \\sqrt{p}$ elements of vector\n",
    "- Total operations are $n^2 / p$\n",
    "- Total time for sending and writing data are approximately $t_s \\log p + t_w (n/\\sqrt{p}) \\log p$, where $t_s$ time unit for sending and $t_w$ time unit for writing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Packages supported distributed storage\n",
    "\n",
    "- [ScaLAPACK](http://www.netlib.org/scalapack/)\n",
    "- [Trilinos](https://trilinos.org/)\n",
    "\n",
    "In Python you can use [mpi4py](https://mpi4py.readthedocs.io/en/stable/) for parallel programming of your algorithm.\n",
    "\n",
    "- PyTorch supports distributed training and data storage, see details [here](https://pytorch.org/tutorials/intermediate/dist_tuto.html) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Summary on large unstructered matrix processing\n",
    "\n",
    "- Distributed manner of storage\n",
    "- MPI\n",
    "- Packages that use parallel computrations\n",
    "- Different blocking strategies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Sparse matrices intro\n",
    "\n",
    "- For dense linear algebra problems, we are limited by the memory to store the full matrix, it is $N^2$ parameters.\n",
    "\n",
    "- The class of **sparse** matrices where most of the elements are zero, allows us **at least** to store such matrices.\n",
    "\n",
    "The question if we can:\n",
    "\n",
    "- solve linear systems\n",
    "- solve eigenvalue problems\n",
    "\n",
    "with sparse matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Plan for the next part of the lecture\n",
    "\n",
    "Now we will talk about **sparse matrices**, where they arise, how we store them, how we operate with them.\n",
    "\n",
    "\n",
    "- Formats: list of lists and compressed sparse row format, relation to graphs\n",
    "- Matrix-by-vector product\n",
    "- Parallell processing of sparse matrices \n",
    "- Fast direct solvers for Gaussian elimination (start)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Applications of sparse matrices\n",
    "\n",
    "Sparse matrices arise in:\n",
    "\n",
    "- partial differential equations (PDE), mathematical modelling\n",
    "- graphs mining, e.g. social networks analysis\n",
    "- recommender systems\n",
    "- wherever relations between objects are \"sparse\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Sparse matrices are ubiquitous in PDE\n",
    "\n",
    "The simplest partial differential equation (PDE), called \n",
    "\n",
    "**Laplace equation**:  \n",
    "\n",
    "$$\n",
    "   \\Delta T = \\frac{\\partial^2 T}{\\partial x^2} + \\frac{\\partial^2 T}{\\partial y^2} = f(x,y), \\quad x,y\\in \\Omega\\equiv[0,1]^2,\n",
    "$$\n",
    "\n",
    "$$\n",
    "    T_{\\partial\\Omega} = 0.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Discretization\n",
    "\n",
    "$$\\frac{\\partial^2 T}{\\partial x^2} \\approx \\frac{T(x+h) + T(x-h) - 2T(x)}{h^2} + \\mathcal{O}(h^2),$$\n",
    "\n",
    "same for $\\frac{\\partial^2 T}{\\partial y^2},$\n",
    "and we get a linear system.  \n",
    "    First, let us consider **one-dimensional case**:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "After the discretization of the one-dimensional Laplace equation with Dirichlet boundary conditions we have\n",
    "\n",
    "$$\\frac{u_{i+1} + u_{i-1} - 2u_i}{h^2} = f_i,\\quad i=1,\\dots,N-1$$\n",
    "\n",
    "$$ u_{0} = u_N = 0$$\n",
    "or in the matrix form\n",
    "\n",
    "$$ A u = f,$$\n",
    "and (for $n = 5$)\n",
    "$$A=-\\frac{1}{h^2}\\begin{bmatrix} 2 & -1 & 0 & 0 & 0\\\\ -1 & 2 & -1 & 0 &0 \\\\ 0 & -1 & 2& -1 & 0 \\\\ 0 & 0 & -1 & 2  &-1\\\\ 0 & 0 & 0 & -1 & 2 \\end{bmatrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The matrix is **triadiagonal** and **sparse**  \n",
    "(and also **Toeplitz**: all elements on the diagonal are the same)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Block structure in 2D\n",
    "In two dimensions, we get equation of the form\n",
    "\n",
    "$$-\\frac{4u_{ij} -u_{(i-1)j} - u_{(i+1)j} - u_{i(j-1)}-u_{i(j+1)}}{h^2} = f_{ij},$$\n",
    "\n",
    "or in the **Kronecker product form**  \n",
    "\n",
    "$$\\Delta_{2D} = \\Delta_{1D} \\otimes I + I \\otimes \\Delta_{1D},$$\n",
    "\n",
    "where $\\Delta_{1D}$ is a 1D Laplacian, and $\\otimes$ is a **Kronecker product** of matrices. \n",
    "\n",
    "For matrices $A\\in\\mathbb{R}^{n\\times m}$ and $B\\in\\mathbb{R}^{l\\times k}$ its Kronecker product is defined as a block matrix of the form \n",
    "\n",
    "$$\n",
    "   A\\otimes B = \\begin{bmatrix}a_{11}B & \\dots & a_{1m}B \\\\ \\vdots & \\ddots & \\vdots \\\\ a_{n1}B & \\dots & a_{nm}B\\end{bmatrix}\\in\\mathbb{R}^{nl\\times mk}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "In the block matrix form the 2D-Laplace matrix can be written in the following form:  \n",
    "\n",
    "$$A = -\\frac{1}{h^2}\\begin{bmatrix} \\Delta_1 + 2I & -I & 0 & 0 & 0\\\\ -I & \\Delta_1 + 2I  & -I & 0 &0 \\\\ 0 & -I & \\Delta_1 + 2I & -I & 0 \\\\ 0 & 0 & -I & \\Delta_1 + 2I   &-I\\\\ 0 & 0 & 0 & -I & \\Delta_1 + 2I \\end{bmatrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "##### Short list of Kronecker product properties\n",
    "\n",
    "- It is bilinear\n",
    "- $(A\\otimes B) (C\\otimes D) = AC \\otimes BD$\n",
    "- Let $\\mathrm{vec}(X)$ be vectorization of matrix $X$ columnwise. Then \n",
    "$\\mathrm{vec}(AXB) = (B^T \\otimes A) \\mathrm{vec}(X).$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Sparse matrices help in computational graph theory\n",
    "\n",
    "- Graphs are represented with adjacency matrix, which is usually sparse\n",
    "- Numerical solution of graph theory problems are based on processing of this sparse matrix\n",
    "    - Community detection and graph clustering\n",
    "    - Learning to rank\n",
    "    - Random walks\n",
    "    - Others\n",
    "- Example: probably the largest publicly available hyperlink graph consists of 3.5 billion web pages and 128 billion hyperlinks, more details see [here](http://webdatacommons.org/hyperlinkgraph/) \n",
    "- More medium scale graphs to test your algorithms are available in [Stanford Large Network Dataset Collection](https://snap.stanford.edu/data/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Florida sparse matrix collection\n",
    "\n",
    "More sparse matrices you can find in  [Florida sparse matrix collection](http://www.cise.ufl.edu/research/sparse/matrices/) which contains all sorts of matrices for different applications."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"700\"\n",
       "            height=\"450\"\n",
       "            src=\"http://yifanhu.net/GALLERY/GRAPHS/search.html\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.IFrame at 0x10d501f70>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import IFrame\n",
    "IFrame(\"http://yifanhu.net/GALLERY/GRAPHS/search.html\", width=700, height=450)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Sparse matrices and deep learning\n",
    "\n",
    "- DNN has a lot of parameters\n",
    "- Some of them may be redundant\n",
    "- How to prune the parameters without significantly accuracy reduction?\n",
    "- [Sparse variational dropout method](https://github.com/ars-ashuha/variational-dropout-sparsifies-dnn) leads to significantly sparse filters in DNN almost without accuracy decreasing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Sparse matrix: construction\n",
    "\n",
    "- We can create sparse matrix using **scipy.sparse** package (actually this is **not the best** sparse matrix package)\n",
    "\n",
    "- We can go to really large sizes (at least, to store this matrix in the memory)\n",
    "\n",
    "Please note the following functions\n",
    "- Create sparse matrices with given diagonals ```spdiags```\n",
    "- Kronecker product of sparse matrices ```kron```\n",
    "- There is also overloaded arithmectics for sparse matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x129c668b0>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD8CAYAAACxd9IeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAaDklEQVR4nO3dbWxc93Xn8e+ZB3JGDyPRjWTZsaQkrWMlMhArYrdKt0CN2rvIZmUrL5o6RdOoD6nedFG3tVHYu3CAwin8ZtfrAC4KeJMWSVOkENxgZSiLtrbSBIWBGpUgpY0tOW6b2vTKspSGpCRySA5nTl/Mg0YShxzOfZh75/4+ADGcIfmfv2ke/e/933PONXdHREZfbtgTEJF4KNhFMkLBLpIRCnaRjFCwi2SEgl0kI2ILdjP7uJm9bmb/ZGaPxfW+/TKznWb2N2Z21sxeNbOHW6/fYmYvmtkbrceJYc+1zczyZnbazI63nid5rlvN7HkzO9f6HX8sqfM1s99p/Q18z8y+bmalpM51PWIJdjPLA38I/Bfgw8AvmtmH43jvdVgGHnH3DwEHgN9szfEx4IS73wmcaD1PioeBs13PkzzXLwJ/6e57gI/QnHfi5mtm7wV+C5h097uBPPBpEjjXdXP3yD+AjwF/1fX8ceDxON47wJyPAf8JeB24rfXabcDrw55bay530Pyj+zngeOu1pM61AvwAsBteT9x8gfcCU8AtQAE4DvznJM51vR9xHca3f4Ftb7deSyQzex+wD3gFuNXd3wFoPW4f4tS6PQP8HtDoei2pc/0AcAn4k9Zpx5fMbCMJnK+7/3/gfwJvAe8As+7+1yRwrusVV7DbCq8lMk/XzDYBfwH8trtfHvZ8VmJmB4GL7n5q2HPpUwH4KPBH7r4PmCOhh8Gtc/FDwPuB24GNZvaZ4c4qHHEF+9vAzq7ndwDnY3rvvplZkWag/5m7f6P18rtmdlvr67cBF4c1vy7/EXjQzP4V+HPg58zsayRzrtD8//+2u7/Sev48zeBP4nzvB37g7pfcvQZ8A/hpkjnXdYkr2P8euNPM3m9mYzQ3PF6I6b37YmYGfBk46+5Pd33pBeBw6/PDNM/lh8rdH3f3O9z9fTR/l99y98+QwLkCuPsFYMrM7mq9dB/wGsmc71vAATPb0PqbuI/mZmIS57o+MW58fAL4PvDPwP8Y9mbFCvP7GZqnFv8AnGl9fAL4MZobYW+0Hm8Z9lxvmPe9XNugS+xcgXuAk63f7/8FJpI6X+D3gXPA94A/BcaTOtf1fFjrP05ERpwy6EQyQsEukhEKdpGMULCLZISCXSQjAgX7IJVsZnYkyHvGLU3zTdNcIV3zTdNcexk42ANUsqXtl5am+aZprpCu+aZprisKsrL/B+Cf3P1f3H2JZtrmoXCmJSJhGzipxsx+Hvi4u3+u9fyXgZ9y9//W62eKG7d4vrIdB3Jm5HPGB2/dRM5WqpPpT8Od7797lXrDQxnvRpcuXWLbtm2hjRelNM0V0jXftMz11KlTP3T3FSdaCDBuX5VsrXOdIwDjO36CHYef6XytVDAefeBuDt1zOxvHB5vKqTen+eyXX2FuqU65mOPRg3sDjSeSZmb2Zq+vBTmM76uSzd2fc/dJd5+88V+ChWXnyeOvcf/T32FucXmgSezZsZlKuUi5mGNxucEXvhlsPJFRFSTY113JttLhdbVWZ2Z+idNTM5x6c3rdQbpxvMBLv/uzPHFwL6Vinvml5njHzpxXwIt0CVQIY2afoNkxJQ/8sbv/wWrfv+mOu/yOX/nfLCxf/54GbK+Mc3VhmUq5yEu/+7PrPgyfW1zm/qe/w8z8EovLDUrFPFsGHEskrczslLtPrvS1QNfZ3f3/ufsH3f3H1wp0gA/euonPP3A35WL+uteLeWNmfom5pTrTc4sDrcpa4UVWF2sGXc6MQ/fcztYNRUqFa4f0S3WnVm+u9kHO4zeOFzh0z+1s0Tm8yE1iT5dtr8DdK/x4IcdY/tpUqrU6s9Ua5y5cGXj87hV+0LFERslQcuPbK/DWDUU2jOXZuqFIpVxkLG8YzeDfOJ5noVYfaEXuXuE3jOXZXCoMPJbIqIi1U83k5KSfPHmy83xucZlzF66wc6LMg8++zPT8IrW6U8wZtYYH3mSbW1zm9NQMjx79LlcWagNv/omkRWQbdEFtHC+wf/cEU9NVrizUWFx2Gg6L9ebjfIANu/b45WKeKws15lqH84Ne4hNJu0SUuLYTY7o37dqCJt60x24fzj969Lt89suvaNNOMicRwb7Spl23ai34Jbk//fWf4ulfuKezyuuynGRNIoIdrt+0C3uFb58u3LNzq1JrJbMSE+wQ7QrfPb4SbySLEhXsEO0K3z2+Em8kaxIX7NDfCq/iGZH1Gep19rW0i1um5xZVPCPSh8ReZ1/Laiv8WCHHbPXa9XOl1oqsLtHBDtefw5eLuU46baVcoFIqMl7IsblUYM+OzYHGV2qtjLrEBztcW4H/z+GfZHtlHDMwrPnY+jyM8Z/77CSGceSrJ7VhJyMn1mBvuA+cqtpOfb26sMxCrcFstcZstdb6vLm5dvHyQuDxlVoroyrWDbpNd9zlu3/ti4E31GarNTaXChjGbLW5uTZeyAXeZFtpfBXQSJqstkEX619vveHMLdVpuHPszPl1d4FtH26fu3Clc45+7Mx5njz+KtVaA2gWzzQajcDjL9TqHPnqyUDzFUmS2Ff2nb/6TKiXuS5eXuDAUydo3PCfUS426+SDlMfqspykTWIuvX3w1k2hJ7JMTVcpKbVWZE1D6UEXZqrqnh2b2RJReaxSa2WUDK0H3Y2JLEFTX1U8I7K6oaXLhr3zvVpqLegcXrIhMefs3cJuKqHiGZHVJaIQJsxVU8UzkmWJXNm7hVmQouIZkZUlItgh3IKUfopn1JdesiYRh/Hdwuz13h7rkaNnmK3WqJSKmMHM/BK1uvrSy+hJ/GF8txsLUoJsgq1UPNM811ZfesmexAU7XOv1HkYiS3ff+Eq50LmBZJv60ktWJDLYw7zM1X2J75mH9oWeWqu+9JIWiQx2CDdVtbtvfBSptepLL2mQ2GCH8BNZlForWZboYIfwi1HUl16yKvHBDsMpnpmZXxo4UUYrvCTRmsFuZjvN7G/M7KyZvWpmD7dev8XMXjSzN1qPE1FOdKVEliA732ut8IvLDXZOlEOZr1Z4SYJ+VvZl4BF3/xBwAPhNM/sw8Bhwwt3vBE60nkcqzuKZUjHP9y9e5dSb0wM3stQKL0my7gw6MzsGPNv6uNfd3zGz24Bvu/tdq/3seu8Is5ooime6x2qX3V6uLrEQwfgqnpEohNZw0szeB+wDXgFudfd3AFoBv73HzxwBjgDs2rVrPW+3qvaqeezMeb7wzdeYX6oDcO7CFfbvXt8ZRXejyZ0TZaamq52Gk/NdjSyDjh/GXEUG1fcGnZltAv4C+G13v9zvz7n7c+4+6e6T27ZtG2SOPYVdPLN/9wTbK6XONfNNpQJjeesU0qh4RtKsr8N4MysCx4G/cvenW6+9zhAP47tFUZAyt7jMff/rO0zPL1KrO8WcUWuoeEaSLVAhjJkZ8GXgbDvQW14ADrc+PwwcCzrRQYVZPNN27sIVrizUOkUzi3UVz0i6rbmym9nPAH8L/CPQaL3832metx8FdgFvAZ9y9x+tNlZUKzuEvwkWZU873XlGorLayp64evYg5haXr9sEKxdzPHFw78B3cmmP9+Tx16jW6jd9vVQwPv/A3QONP7e4fNOdZ4LOVyRV9exBpCm1VsUzEreRCnZQ8YxILyMX7JCuFT6K+YqsZCSDHYZXPKO+9JJUI7VBt5I47zyjvvQybJnZoFtJnMUzxbwxM7/EXIBr8VrhJSojH+wQ/s53r3P4pbp3GloG3anXObyELRPB3hb1nWfGCznG8td+pdWa7jwjyZGpYIfo7jyzYayZUVcpF1U8I4k08ht0vYR955l2eeyDz76s4hkZmkxv0PUS9p1n9u+eYGq6quIZSazMBjuEe+eZ7vGiSLzRnWckqEwH+zBSawftWqs7z0hQmQ52iD+1NkjXWhXPSBCZD3aIN7V2vJDjpbMXB+5Y22u+WuFlLZndjV9JVKm17dTX8ULuusewGmwotVbaQusuO+q6u8x2N5VouHPszPl1N5W4sWvtS2cv8uTxV6l2dax1d05PzVAu5tmzY/NA43c37PAB5yqjTyt7D1GsmhcvL3DgqRM0un7lKp6RMOk6+wCiSFedmq7edH/4sUKO2eq16+dKrZWoKNhXEXa66p4dmzu7/u102kq5QKVUVGqtRE6H8X0IO7X29NQMjxw9w2y1RqVUxAxm5peo1ZVaK8HoMD6gsFNry8U8VxeWWag1mK3WWufaSq2VaCnY+xRmam136mulXOjUwLcptVaioGDvU5iJLN2pr888tO+mTTsI1rVWqbWyEgX7OoSZWtud+rolguIZpdbKjRTs66S+9JJWCvYBqC+9pJGCfUDD6EuvFV6C0HX2gOLsSw/B7h7bPb5Sa0eTrrNHKM6+9KA7z8jgtLKHKMxVU3eekUFoZY9J1H3p21Q8I4NQsIcsqr70Kp6RoHQYHxEVz8gwhHIYb2Z5MzttZsdbz28xsxfN7I3W40RYEx4FcRfPDNq1dqW5qnhmNK3nMP5h4GzX88eAE+5+J3Ci9Vy6xFk8E6Rr7Y3jq3hmNPUV7GZ2B/BfgS91vXwI+Err868Anwx1ZiMgzuKZoF1rVTwz+vo6Zzez54GngM3Ao+5+0Mxm3H1r1/dMu/uqh/JZOmfvFvZlLnWtlV4CnbOb2UHgorufGvDNj5jZSTM7eenSpUGGSL2oime+9rkD/N3j9/HEwb2MF3JUa43QGmAo8Wb0rLmym9lTwC8Dy0AJqADfAH4SuNfd3zGz24Bvu/tdq42V1ZW9LaoVc6WutaDU2iwKtLK7++Pufoe7vw/4NPAtd/8M8AJwuPVth4FjIc13ZIVdPNO2UtdaUGqtXG9d19nN7F6unbP/GHAU2AW8BXzK3X+02s9nfWVvi7N4Rqm12RJauqy7f9vdD7Y+/zd3v8/d72w9rhrock2cxTPFvDEzv8RcgPN4rfCjQemyQxJ226heDTCW6t65Jh+0zZUaYKSbgn3Ioi6eGS/kGMtf+99cral4JqsU7AkQVfHMhrHmbnylXGQsb51CGhXPZJMKYRIk7OKZ9t1jH3z2ZabnF6nVnWLOqDVUPDOqVM+eEmEXz+zfPcHUdJUrC7VO0cxiXXeeySoFe8KEWTzTPV4UXWtVPJMuCvaESVNfehXPpIuCPYHS1Jded55JDwV7QqVphY9ivhI+BXuCpWmFj2K+Ei4Fe8IN484zKp4ZTbrOnhIqnpF+6Dr7CBhW8cygjSy1wiePgj1FhlE8E6SRpc7hk0XBnkJxFs+UinmmpquJmKsEo2BPqaiLZ7ZuGLtu3EG71q401y3lIjsnykqtjZk26FIuiuKZPTs2A3TGvVxdYiGkrrXt4pxDf/gyl6sqoAmbNuhGWBTFMxvHC51xL1eXmA+xa227OOdyVam1cVOwj4Cwi2fadk6UWVhuXPda0MSbKOcrq1Owj4CoLnOt1rVWqbXpo2AfEVFc5tqzYzNbIiqP1WW5+CnYR8gwUmu1wqeHduNHUJyptaA7zySJduMzJs7UWlDxTFpoZR9xYa6aKp5JPq3sGRZ1am3bWCHHbPVa80n1pU8eBXsGRJVaWy7mOr3oK+UClVJRfekTTIfxGRJ2au3pqRkeOXqG2WqNSqmIGczML1Grqy/9sOgwXoDwU2vLxTxXF5ZZqDWYrdZa59rqS59UCvaMCTNVtbtvfKVc6NTAt6kvfbIo2DMmzMtc3Zf4nnloX+iptepLHy4FewaFmara3T0nitRa9aUPj4I9o9SXPnsU7BmmvvTZomDPuGH1pR80UUYr/OD6CnYz22pmz5vZOTM7a2YfM7NbzOxFM3uj9TgR9WQlGislsgTZ+V5rhQ/SsfbG+WqF71+/K/sXgb909z3AR4CzwGPACXe/EzjRei4pFWfxTKmY5/sXr3LqzemBG1lqhV+/NTPozKwCfBf4gHd9s5m9Dtzr7u+Y2W3At939rtXGUgZdOkRRPNM9VrvsNoxGliqeud5qGXT9/DY+AFwC/sTMPgKcAh4GbnX3dwBaAb+9x5sfAY4A7Nq1a4DpS9zaq+axM+f5wjdfY36pDsC5C1fYv3t9Z2vtsdpdZaemqyzU6hz56knma83+dmGMH8ZcR10/h/EF4KPAH7n7PmCOdRyyu/tz7j7p7pPbtm0bcJoSt7CLZ/bvnmB7pdS5Zr6pVGAsb51CGhXPRK+fYH8beNvdX2k9f55m8L/bOnyn9XgxminKsLRXzec+O4lhHPnqydA2wQzDDMwAd/7t6hK/EWD8KOc6KtYMdne/AEyZWft8/D7gNeAF4HDrtcPAsUhmKEMVZvFM27kLV7iyUOsUzSzWVTwTh75KXM3sHuBLwBjwL8Cv0vyH4iiwC3gL+JS7/2i1cbRBl05hb4JF2dMu7P57abPaBp3q2aUvc4vL122ClYs5nji4l0P33D5wwB87c54nj79GtVa/6eulgvH5B+4eaPz2babaG4FzIcw3LVTPLoGlKbVWxTMrU7BL31Q8k24KdlmXNK3wUcw3zRTssm7DKp5RX/pgtEEnA4vzzjPqS98fbdBJJOIsninmjZn5JeYCXIvP+gqvYJdAwt757nUOv1T3TkPLoDv1WT2HV7BLKKK+88x4IcdY/tqfa7WmO8+sl4JdQhPVnWc2jDUz6irloopnAtAGnYQu7DvPtMtjH3z2ZabnF6nVnWLOqDV055kbaYNOYhX2nWf2755garqq4pmAFOwSiTDvPNM9XhSJN1m584yCXSIxjNTaQbvWZuXOMwp2iUzcqbVButZmoXhGwS6RijO1dryQ46WzFwfuWNtrvqOywms3XmIRVWptO/V1vJC77jGsBhtpS60N2l1WJLDuLrPdTSUa7hw7c37dTSVu7Fr70tmLPHn8VapdHWvdndNTM5SLefbs2DzQ+N0NO3zAuSaFVnaJXRSr5sXLCxx46gSNrj/nLBbP6Dq7JEoU6apT09Wb7g8/VsgxW712/TzrqbUKdhmKsNNV9+zY3Nn1b6fTVsoFKqWiUmtbdBgvQxV2au3pqRkeOXqG2WqNSqmIGczML1GrZyO1Vofxklhhp9aWi3muLiyzUGswW621zrWVWgsKdkmAMFNru1NfK+VCpwa+LcuptQp2GbowE1m6U1+feWjfTZt2EKxrbZpTaxXskghhptZ2p75uiaB4Jq2ptQp2SQz1pY+Wgl0SRX3po6Ngl8QZRl/6LKzwus4uiRVnX3oIdvfY7vGHmVqr6+ySSnH2pYfRv/OMVnZJhTBXzVG+84xWdkm9qPvSt41y8YyCXVIjqr70WSme0WG8pI6KZ3oLfBhvZr9jZq+a2ffM7OtmVjKzW8zsRTN7o/U4Ee60RVYWd/HMoF1rV5rrMItn1gx2M3sv8FvApLvfDeSBTwOPASfc/U7gROu5SCziLJ4J0rX2xvGHWTzT7zl7ASibWQHYAJwHDgFfaX39K8AnQ5+dSA9xFs8E7VqblOKZvs7Zzexh4A+AKvDX7v5LZjbj7lu7vmfa3Vc9lNc5u4Qt7Mtcae9aG+icvXUufgh4P3A7sNHMPrOONz9iZifN7OSlS5f6/TGRvkRVPPO1zx3g7x6/jycO7mW8kKNaa4TWAGNYiTf9HMbfD/zA3S+5ew34BvDTwLtmdhtA6/HiSj/s7s+5+6S7T27bti2seYt0RFE8s3/3BNsrJe7/0HYWlxvXfT2txTP9BPtbwAEz22BmBtwHnAVeAA63vucwcCyaKYqsLezimbaVutZCOlNr+z1n/33gIWAZOA18DtgEHAV20fwH4VPu/qPVxtE5u0QtzuKZJKbWrnbOrqQaGTlzi8s33XmmXMzxxMG9A93NZW5xmWNnzvPk8deo1uqd18fyhhksLjulgvH5B+4ONH77zjNB5qrceMmUsNtG9WqAsVT3zjX5oG2u4jiHV7DLyIq6eGa8kGMsfy2EqrVkF88o2GWkRVU8s2Gs2eiiUi42D+ch8cUzOmeXTAi7eKZ999gHn32Z6flFanWnmDNqjeEWz+icXTIv7OKZ/bsnmJqucmWh1imaWawn+84zCnbJjDCLZ7rHi6JrbRTFMwp2yYw09aWPonhGwS6Zkqa+9GFfQlSwS+akaYUPc74KdsmkNK3wYc1XwS6ZNYw7zwyzeEbX2SXzRql4ZlOpqOvsIr3EeeeZYt6YmV/qjB9Wam17rliuZ0wr2EUYTvFMkEaWvc7hi+/ZvbfXzyjYRbrEWTxTKuaZmq6GOlfL5XqeEyjYRW4QdfHM1g1j1407aNfalebqjUbPQbRBJ9JDFMUze3ZsBuiMe7m6xEJIXWvPXbjC5Pvfc9ob9Y+u9D1a2UV6iKJ4ZuN4oTPu5eoS8yF2rd2/ewK80ej1PQp2kVWEXTzTtnOizELIXWvXomAXWUVUXWBX61obZIVfjYJdZA1R9Ijbs2MzWyIqj+1FwS7Sh2Gk1oa9wms3XmQd4kytBSgXm5fr1JZKJGZxptZC8OKZblrZRQYU5t1cwiqe0couEoGoU2vbxgo5ZqvXmk8O2ktewS4SQFSpteVirtOLvlIuUCkVA/el12G8SAjCTq09PTXDI0fPMFutUSkVMYOZ+SVq9dX70uswXiRiYafWlot5ri4ss1BrMFuttfYFgvWlV7CLhCTM1NruvvGVcqFTA982SOKNgl0kJGGm1nZf4nvmoX2hpNYq2EVCFGZqbXf3nDBSaxXsIiEbZl969aATidmw+tKrB53IEAyjeEY96ESGZKWkmyB3ZF1rhQfTYbzIsPQqnokitRa8Z1uqWDPozOwSMAf8MLY3De49pGe+aZorpGu+4czVcrnie3bvtVyu4I3Gcu2Hb766Wt+4fsdqrujeqM2+S6N65eZjfGIOdgAzO9krnS+J0jTfNM0V0jXfNM21Fx3Gi2SEgl0kI4YR7M8N4T2DSNN80zRXSNd80zTXFcV+zi4iw6HDeJGMULCLZISCXSQjFOwiGaFgF8mIfwdjkjO8nqqopwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.sparse\n",
    "from scipy.sparse import csc_matrix, csr_matrix\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.linalg\n",
    "import scipy.sparse.linalg\n",
    "%matplotlib inline\n",
    "n = 10\n",
    "ex = np.ones(n);\n",
    "lp1 = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr'); \n",
    "e = sp.sparse.eye(n)\n",
    "A = sp.sparse.kron(lp1, e) + sp.sparse.kron(e, lp1)\n",
    "A = csc_matrix(A)\n",
    "plt.spy(A, aspect='equal', marker='.', markersize=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Sparsity pattern \n",
    "\n",
    "- The ```spy``` command plots the sparsity pattern of the matrix: the $(i, j)$ pixel is drawn, if the corresponding matrix element is non-zero.\n",
    "\n",
    "- Sparsity pattern is really important for the understanding the complexity of the sparse linear algebra algorithms.  \n",
    "\n",
    "- Often, only the sparsity pattern is needed for the analysis of \"how complex\" the matrix is."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Sparse matrix: definition\n",
    "\n",
    "- The definition of \"sparse matrix\" is that the number of **non-zero elements** is much less than the total number of elements. \n",
    "\n",
    "- You can do the basic linear algebra operations (like solving linear systems at the first place) faster, than if working for with the full matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## What we need to find out to see how it actually works\n",
    "\n",
    "- **Question 1:** How to store the sparse matrix in memory?\n",
    "\n",
    "- **Question 2:** How to multiply sparse matrix by vector fast?\n",
    "\n",
    "- **Question 3:** How to solve linear systems with sparse matrices fast?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Sparse matrix storage\n",
    "\n",
    "There are many storage formats, important ones:\n",
    "\n",
    "- COO (Coordinate format)\n",
    "- LIL (Lists of lists)\n",
    "- CSR (compressed sparse row)\n",
    "- CSC (compressed sparse column)\n",
    "- Block variants\n",
    "\n",
    "In ```scipy``` there are constructors for each of these formats, e.g. \n",
    "\n",
    "```scipy.sparse.lil_matrix(A)```."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Coordinate format (COO)\n",
    "\n",
    "The simplest format is to use **coordinate format** to represent the sparse matrix as positions and values of non-zero elements.\n",
    "\n",
    "```python\n",
    "i, j, val\n",
    "```\n",
    "\n",
    "where ```i, j``` are integer array of indices, ```val``` is the real array of matrix elements. <br>\n",
    "So we need to store $3\\cdot$ **nnz** elements, where **nnz** denotes number of nonzeroes in the matrix.\n",
    "\n",
    "**Q:** What is good and what is bad in this format?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Main disadvantages\n",
    "\n",
    "- It is not optimal in storage (why?)\n",
    "- It is not optimal for matrix-by-vector product (why?)\n",
    "- It is not optimal for removing elements as you must make **nnz** operations to find one element (this is good in LIL format)\n",
    "\n",
    "First two disadvantages are solved by **compressed sparse row** (CSR) format."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Compressed sparse row (CSR)\n",
    "\n",
    "In the CSR format a matrix is stored as 3 different arrays:\n",
    "\n",
    "```python\n",
    "ia, ja, sa\n",
    "```\n",
    "\n",
    "where:\n",
    "\n",
    "- **ia** (row start) is an integer array of length $n+1$ \n",
    "- **ja** (column indices) is an integer array of length **nnz** \n",
    "-  **sa** (values) is an real-value array of length **nnz**\n",
    "\n",
    "<img src=\"https://www.karlrupp.net/wp-content/uploads/2016/02/csr_storage_sparse_marix.png\" width=60% />\n",
    "\n",
    "So, we got $2\\cdot{\\bf nnz} + n+1$ elements."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Sparse matrices in PyTorch and Tensorflow\n",
    "\n",
    "- PyTorch supports sparse matrices stored in COO format\n",
    "- Incompletre backward operation for these matrices, see summary [here](https://github.com/pytorch/pytorch/issues/9674)\n",
    "- Tensorflow also supports sparse matrices in COO format\n",
    "- The list of supported operations is [here](https://www.tensorflow.org/api_docs/python/tf/sparse) and gradient support is also limited"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### CSR helps in sparse matrix by vector product (SpMV)\n",
    "\n",
    "```python\n",
    "\n",
    "   for i in range(n):\n",
    "        \n",
    "        for k in range(ia[i]:ia[i+1]):\n",
    "            \n",
    "            y[i] += sa[k] * x[ja[k]]\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Let us do a short timing test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7.91 ms ± 453 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "10.8 ms ± 536 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.sparse\n",
    "import scipy.sparse.linalg\n",
    "from scipy.sparse import csc_matrix, csr_matrix, coo_matrix\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "n = 1000\n",
    "ex = np.ones(n);\n",
    "lp1 = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr'); \n",
    "e = sp.sparse.eye(n)\n",
    "A = sp.sparse.kron(lp1, e) + sp.sparse.kron(e, lp1)\n",
    "A = csr_matrix(A)\n",
    "rhs = np.ones(n * n)\n",
    "B = coo_matrix(A)\n",
    "%timeit A.dot(rhs)\n",
    "%timeit B.dot(rhs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "As you see, **CSR** is faster, and for more **unstructured patterns** the gain will be larger."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Sparse matrices and efficiency\n",
    "\n",
    "- Sparse matrices give complexity reduction.  \n",
    "- But they are **not very good** for parallel/GPU implementation.  \n",
    "- They do not give maximal efficiency due to **random data access**.  \n",
    "- Typically, peak efficiency of $10\\%-15\\%$ is considered good.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Recall how we measure efficiency of linear algebra operations\n",
    "\n",
    "The standard way to measure the efficiency of linear algebra operations on a particular computing architecture is to use **flops** (number of floating point operations per second)\n",
    "\n",
    "We can measure peak efficiency of an ordinary matrix-by-vector product."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time:  1.2e-02, Efficiency:  2.6e+00 Gflops\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import time\n",
    "n = 4000\n",
    "a = np.random.randn(n, n)\n",
    "v = np.random.randn(n)\n",
    "t = time.time()\n",
    "np.dot(a, v)\n",
    "t = time.time() - t\n",
    "print('Time: {0: 3.1e}, Efficiency: {1: 3.1e} Gflops'.\\\n",
    "      format(t,  ((2 * n ** 2)/t) / 10 ** 9))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time:  5.8e-04, Efficiency:  2.1e-02 Gflops\n"
     ]
    }
   ],
   "source": [
    "n = 4000\n",
    "ex = np.ones(n);\n",
    "a = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr'); \n",
    "rhs = np.random.randn(n)\n",
    "t = time.time()\n",
    "a.dot(rhs)\n",
    "t = time.time() - t\n",
    "print('Time: {0: 3.1e}, Efficiency: {1: 3.1e} Gflops'.\\\n",
    "      format(t,  (3 * n) / t / 10 ** 9))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Random data access and cache misses\n",
    "\n",
    "- Initially all matrix and vector entries are stored in RAM (Random Access Memory)\n",
    "- If you want to compute matvec, the part of matrix and vector elements are moved to **cache** (fast and small capacity memory, see [lecture about Strassen algorithm](https://nbviewer.jupyter.org/github/oseledets/nla2018/blob/master/lectures/lecture-3.ipynb))\n",
    "- After that, CPU takes data from cache to proccess it and then returns result in cache, too\n",
    "- If CPU needs some data that is not already in cache, this situation is called **cache miss**\n",
    "- If cache miss happens, the required data is moved from RAM to cache\n",
    "\n",
    "**Q**: what if cache does not have free space?\n",
    "\n",
    "\n",
    "- The larger number of cache misses, the slower computations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Cache scheme and LRU (Least recently used)\n",
    "\n",
    "<img src=\"cache_scheme.png\" width=\"500\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### CSR sparse matrix by vector product\n",
    "\n",
    "```python\n",
    "\n",
    "   for i in range(n):\n",
    "        \n",
    "        for k in range(ia[i]:ia[i+1]):\n",
    "            \n",
    "            y[i] += sa[k] * x[ja[k]]\n",
    "            \n",
    "```\n",
    "\n",
    "- What part of operands is strongly led cache misses?  \n",
    "- How this issue can be solved?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Reordering reduces cache misses\n",
    "\n",
    "- If ```ja``` stores sequential elements, then they will be moved to cache altogether and number of cache misses decreases\n",
    "- This happens when sparse matrix is banded or at least block diagonal\n",
    "- We can convert given sparse matrix to banded or block diagonal with *permutations* \n",
    "\n",
    "- Let $P$ be row permutation matrix and $Q$ be column permutation matrix\n",
    "- $A_1 = PAQ$ is a matrix, which has less bandwith than $A$\n",
    "- $y = Ax \\to \\tilde{y} = A_1 \\tilde{x}$, where $\\tilde{x} = Q^{\\top}x$ and $\\tilde{y} = Py$\n",
    "- [Separated block diagonal form](http://albert-jan.yzelman.net/PDFs/yzelman09-rev.pdf) is a cache-oblivious format for sparse matrix by vector product\n",
    "- It can be extended for 2D, where separated not only rows, but also columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Example\n",
    "\n",
    "- SBD in 1D \n",
    "<img src=\"sbd.png\" width=\"400\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Sparse *transpose* matrix-by-vector product\n",
    "\n",
    "- In some cases it is important to compute not only $Ax$ for sparse $A$, but also $A^{\\top}x$\n",
    "- Mort details will be discussed in the lecture about Krylov methods for non-symmetric linear systems\n",
    "- Transposing is computationally expensive\n",
    "- [Here](https://people.eecs.berkeley.edu/~aydin/csb2009.pdf) is proposed **compressed sparse block** format of storage appropriate for this case"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Compressed sparse block (CSB)\n",
    "\n",
    "- Split matrix in blocks\n",
    "- Store block indices and indices of data **inside** each block\n",
    "- Thus, feasible number of bits to store indices\n",
    "- Ordering of the blocks and inside blocks is impoprtant for parallel implementation\n",
    "- Switching between blockrow to blockcolumn makes this format appropriate to transpose matrix by vector product"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Solve linear systems with sparse matrices\n",
    "\n",
    "- Direct methods (start today and continue in the next lecture)\n",
    "    - LU decomposition\n",
    "    - Number of reordering techniques to minimize fill-in\n",
    "- Krylov methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let us start with small demo of solving sparse linear system..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Reshaped solution')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEICAYAAABLdt/UAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJztnXm8JGdZ779P72efJbNmMpkskAARAoxhUZRLAgLGoF5UEBH0cnNdUPCqKBdF3L3IVfDjwo2IXmURZZWoQCIGXNiSEEJghpBAkpnMdmY9a+/P/aOquut019pd3V3d5/1+Pudzzumu6nq7uvr3PvV7n/d5RVUxGAwGw+SQGXUDDAaDwZAsRtgNBoNhwjDCbjAYDBOGEXaDwWCYMIywGwwGw4RhhN1gMBgmDCPsCSIirxSRfx91OwBEREXkygRf74D9mrke93+ZiHwiqfYY+kdEHhKRG1LQjjeJyLsSfs2/EpHf6mP/r4jIsxNs0lBJvbDbF99JEZlxPfYqEbkj4v53iMirBtZAQxdenYCqvltVnzfKdo0r9ndgXURWROSELVqzo27XpODVCajqE1T1jhE1qW9SL+w2OeA1wzhQrxGpwTBgvkdVZ4FrgScDrx9xewwpZlyE/feBXxCRLV5PisgzReQLInLB/v1M+/HfBp4F/LEd7fyxx75OdPnfROQR4JP2408Xkf8UkfMi8iX3bZltuXxDRJZF5Jsi8rKO13yLiJyzn3uB6/EfE5FD9n7fEJH/4Xru2SJyVET+l4ictqO0l7meL9qv+4h9B/N2EZlyPf+LInJcRI6JyI8HnUy/9otIRkR+RUQeFpFTIvLXIrLg8xobbuM7bqc/bf8+b5/3Z3TaVH6fmf3cHSLymyLyH3YbPyEiFwW9p82Cqp4APo4l8EDwtSEiF4nIrfZ1fFZE/k1E3N/7a0XkXvtzeJ+IlOz9ttr7LdrX8q0iss91zDtE5HdF5PP2vh8RkW2u54O+P5eJyKfsz/Y2wPezDWq/iDzObsd5sayTm3xeo8sitb/zV4rIzcDLgNfZ1+pH7edb17d9ft9qf7eO2X8X7eec7+3P29+Z4yLyY4Ef4jBQ1VT/AA8BNwAfBH7LfuxVwB3239uAc8DLsSL7l9r/b7efvwN4VcDrHwAU+GtgBpgCLgbOAC/E6vyea/+/w95mCbjK3n8P8AT771cCNeC/A1ngJ4FjgNjPfzdwBSDAdwJrwFPs554N1IE/AIr286uu47wV+Af7/c4BHwV+137u+cBJ4Bq7fe+x39OVHu83qP0/DjwAXA7M2uf8bzrOU879ubhe903Au7y2dZ2bf4/xmT0IPNb+PO4Afm/U1+KovwP23/uALwNvcz0fdG38LvB2IG//PMt1PT4EfB7Ya+97CPgJ+7ntwH8Fpu3X/Hvgw65j3gE86rrmPuD6/H2/P/bzn3Fd598BLDv7erx3z/bbfz8A/C+gADzHfh3nuv4r2nrRuvZcr9v6fri39TnnvwF8FtiJpQH/Cfxmx/f2N+w2vRDre711pNfMqC/aqBe1fQFdsE+sW9hfDny+Y5/PAK90XYBRhP1y12O/hC1orsc+DrzCvojP2xf9VMc2rwQecP0/bb/2bp9jfxh4TccFMuN6/u+AX7Uv5FXgCtdzzwC+af/9TlzChyWIQcLu1/5/AX7K9f9VWB1VjmSFPcpn9iuu534K+Nior8URfwdWsIRL7c9pi/1c2LXxG8BHfK6Fh4Afcf3/ZuDtPm24Fjjn+v+Ojmvu8UAVK6AJ+v7s97jO34O/sHu2H0vgTwAZ12PvBd5k//1XJCfsDwIvdD33XcBD9t/PBtY7rvVTwNNHec2MixWDqt4H3Ar8csdTe4GHOx57GCtqiMMR19+XAj9g3+KdF5HzwLcDe1R1Ffgh4CeA4yLyjyJytWvfE642r9l/zgKIyAtE5LP2LeV5rN7dfRt6zn599/vYi9WZTQN3udrzMftx5xwc6djPk5D2d57Lh7FEfZff6/VIlM/shOvvNexzuIn5XlWdwxKSq2lfN2HXxu9jRbafsO23zu+P53kWkWkR+b+2LbeEZa9tEZGsa/vOay5vt8v3+4P12Xtd5374tX8vcERVmx2vE/d7HwWv78Ve1/9nVLXu+n/k1+vYCLvNr2HZHO4P7xjWheRmP9ZtIlg9cxTc2x3Biji2uH5mVPX3AFT146r6XKwL9TDw52EvbntyHwDeAuxS1S3AP2FFXA5bxZX9Y7+PY8BprKjgCa72LKg1mAZwHLikYz//N+rf/s5z6URXJz1eZhVLUBx2uw8RdHyP4zjHetRjW4MLVf0UVoT5FvuhwGtDVZdV9edV9XLge4D/KSLXRzjUz2PdsT1NVeexLBPYeL12XnM1uz1B35/jeF/nfu/Xr/3HgEs6xgv8rqEN16qI7O54Pu716nwvU8tYCbuqPgC8D/hZ18P/BDxWRH5YRHIi8kNYt4W32s+fxPKM4/Au4HtE5LtEJCsiJXuQZJ+I7BKRm+wLs4J1i9yI8JoFLE9xEaiLNajqlf736yJSEJFnATcCf29HJX8O/KGI7AQQkYtF5Lvsff4OeKWIPF5EprE6QE9C2v9e4Ofswa1Z4HeA93VEIw73AC8RkbyIHARe7HpuEWjif97DPjNDMG8Fnisi14ZdGyJyoz1IKFhjKw2iXa9zWB3GeXtQ1Oua+hHXNfcbwPtVtUHA90dVHwbupH2dfzuWYHsS0P7PYQn26+xr8Nn26/ytx8t8CXiCiFwr1uDwmzqeD9OI9wK/IiI7xBrEf6P9HlPLWAm7zW9g+cQAqOoZLAH8eawBmtcBN6rqaXuTtwEvFmtk/4+iHEBVjwAvwhqYWcSKQH4R63xl7GMdA85iDXL+VITXXMbqkP4Oa6Dwh7EGvNycsJ87BrwbayDrsP3cL2Hdkn7WvjW+HSuiQlX/GevL/kl7m08GNCWo/e8E/gbrtvubQBn4GZ/X+VWsgeBzwK9j+aTOe10Dfhv4D/tW/Okd5yLsMzMEoKqLWIP9v2o/5HttAI+x/1/BGsf4U42Wn/1WrIHr01gDhx/z2OZvsO4eTgAl7IAr5PsD1rX/NKzr79fs9+KHZ/tVtQrcBLzAbuOfAj/q+r60UNX7sXTjduDrQOckwr8AHm9fqx/2aMNvYXVG92INXN9tP5ZanNFxw4ixI453qeq+sG0NhlEj1gTBd6nqO0bdFkM34xixGwwGgyEAI+wGQwAi8hoRuU+sCTCvHXV7DIYo9GXFiMgPYA1EPA64TlXvTKhdBsPIEZFrsAbjrsPK0f4Y8JOq+vWRNsxgCKHfiP0+4PtpTyE3GCaJxwGfVdU1OzPoU8D3jbhNBkMofRW8UtVDAFYmUnQuuugiPXDgQD+HNhh8ueuuu06r6o7wLUO5D/htEdmOlfr3QqzsiA2IVW/kZoDpaXnqZVfE+1qda0yHbzQkluqlUTdhYMznyqNuwga2ZtfCN+rgK1+uRbq2h1bJ0H3x79+/nzvvNK6NYTCISNBMxsio6iER+d/AbVjpdl/CmrDVud0twC0A1zyxoB/4x3j1yt6/9JT+G9sHt524OnyjCeO5u7uyIofOi+fvjr3P1fuPR7q2Q4VdRG5n46xChzeo6keiNsh98R88eNDkWBrGAlX9C6w8Z0Tkd4CjSb7+KER9Mwp5J53nYBRC//6lp/Qk7lEIFXZVHfkKKwbDqBCRnap6SkT2Y40nPSOp1x6WqBshD2dUQj8ocTeLShgMwXzA9thrwE+r6rlRNygMI+T9k4aIvh/6yooRke8TkaNYUcw/isjHk2mWwZAOVPVZqvp4VX2Sqv5LUq87qGjdiPpguO3E1QM7t4O4FvoSdlX9kKruU9Wiqu5S1e8K38tg2NwM4os8SOExtBkXcTdWjMEwxhgxHz7OOU+zPWNKChgMQyTJyMyI+mhJ+i4pyWsjVRH7O/7tG+zbOs3zr/HKruymUm/wkXuO8QNP3Rd7kpQhPXz4i49y/eN2MlfKR9r+0/cvcs+R8/zMc64cq889qS9umgX9oaNJzAvz5sC+xYG9dj8kGcEnlSWTKmH/f595iG+9dFtkYf/Xw4u87v338rjd83zLvoXBNs4wEI6eW+O177uH3/3+b+Gl1wUu/NTiU/cv8r4vHOFnr3/MgFuXLtIg6IMU7n6PPWrhv+3E1amxZ1Il7HPFPEvlWuTtL6xXAThv/zaMH+fXaht+R2FpvcZcKVWXbij9ROujEPRRCnivBLV5WKKfRPSeRNSeqm/H/FSOpXWvVdi8cbaNs48hXTgdeZwOfalcYz6ibZMGxkXUx1HMo+J+b8MQ+X4Fvl9xT5ewl/I8cjZ6YZxeRMGQLtqdc5yIvc78VKou3cQZlqBPspj7MUyRH5U9k6pvx/xUPuYXvLbht2H8aHfOMe7UyjV2z49HFcJeovVBi/pmFHM/hiHyvUbv/UTt6RL2Uj7mF9yO9kzEPrb00jkvlWs8dtfcoJqUGGkSdSPm4QzbrolCr+Keqjz2+akcK5U6jWa04o9tUTAe+7jSS+e8XK4zn/LB0zSI+kNHd7R+DPEYxLnr9fPt5VpKlbA7ecwrEaP2ZROxjz1O57wc8TNXVTsrZnwGT4eNEfNkSfJ8DmvsJFXC7kRhUYXa2S6qKBjSR6tzjmjFrFYbNJVUD56OKlo3gj5YRnlu415Tqfp2zE9ZUdiF9RqXRNjeDJ6OP3Ezm5zPOq3pjr0sc5eUqKeJ4iOFxF6rsj8981QeOrqjb/+910wZS9z/MdK26RJ2+8saPWI3Vsy44wh1udakUm9QzGWDt7c/aycIGHfGXdSTFPC4xxiV4Dvnux+BH3QaZLqE3b69jjIYWm80WamYCUrjjjsLarlcpzgbIuz2Zz2siF1Efg54FaDAl4EfU9VUrIo8TEEfhoDHZdSC32/0PkhxT5nHHj1id0Q9nxUTsY8xS+s18llp/R1lexiOxy4iFwM/CxxU1WuALPCSpF6/n2h90KJefKSw4Wec6Gz7INufNgvMIV3Cbt9eRxkMdba5eMsUa9UGtUZzoG0zDIalco2Lt0wBET/3ytA99hwwJSI5YBo4lsSLplHUx1XIozDI99XPoPWgsmRSJeyzRceKCY/cLtjbXLw1uigY0kWzqaxU6q3PMMqdl2PFDKMImKo+CrwFeAQ4DlxQ1U8M/MA+DCrrZVLF3ItBdl5pEvd+1zz9fRE5LCL3isiHRGRLP6+XzQhzxVy0L7i9zb4tVhaCyYwZP5YrdVTdn2F45+x8zsPIYxeRrcCLgMuAvcCMiPyIx3Y3i8idInLn6rlwf7eXL3LSgj7J0XlUBvH+02LN9Bux3wZco6pPBO4HXt9vg6x6MVG+4LYVYyL2sWW5vPGuK2qHPpXPUsgN5WbzBuCbqrqoqjXgg8AzOzdS1VtU9aCqHpzZmp5I0I/NLOZeTKK497uY9SdU1VHUzwL7+m3QXClmxB5DFAzpwumcW59hpMHToVZ2fAR4uohMi7VU0/XAoWEdPGk2e4QeRNLnJq64J23HJBn2/Djwz35Pum9XFxf9U4SiVnh0ttm31Vgx44rTGe+eL5HNRMtuGmYtdlX9HPB+4G6sVMcMcEs/rxn3C5xU9GcEPRqTcp5CQx8RuR3wWqvuDar6EXubNwB14N1+r6Oqt2B/KQ4ePOhb5Wu+lOfR8+thzWKpXEcE9m4p2f8bYR832qmLeeZL0RZZWSrXhjo5SVV/Dfi1oR3QRRKiPilCNUycc9ZvPnwSs1R7JVTYVfWGoOdF5BXAjcD1qhqtLGMA86UchyOI9HK5xmwxx5Zp60Mwk5TGD2dy0sJUnvmpfMtzD2K5XGfbjBGrKBhR74/iI4WhinuSE5b6zYp5PvBLwE2qGn3powCiWzF15kt5ZgpZMmIi9nHEXfclai3+pfXxWhbPTRwbpt9o3Yh6MozruES/o1B/DBSB26yxJT6rqj/RzwvOl3Is2zXZsxnx3e6CvaCxiDBXirfy0qHjS7z6PXfz/p94JltN9JcI5VqD7/vT/+SNNz6eZ1yxPdI+Tmc8W8oxP5VrzU0IwvncJ5lxFvW5h/u+afdl+VJ/PUgzo7Bk+s2KuVJVL1HVa+2fvkQdYPtsEVU4txZ8C3RmtcKOuSIA04Usa9VG5GMcOr7Eg4ur3H9yua+2GtocObvGoeNLfOXYhcj7rFcbFHMZshlh+0yRMyuVwO1rjSbn12tcNFvst7kTyyhEfe5hbf1MwnG86Pe8DjsFMlUzT4GWWJ9aCv6Sn1qqsMP+gpfyWcr16CUFyjVr21PLwccwRMc5l5VYn0ODUt4q+rVjrhj6eZxZqaLavkYmkTTkQEdhlCLbefxhtWEYnWZSaY+pu6fd6Qj7cpnHM++5jaqyuFxhx7y1bTGXoVyLHrE72y4aYU8M51zG+xyalPJWbLFzrshatcFKpd4qLdHJqeVya9txYxgr5wxSeEYl4FFxt2+Qlk0SA6rDIIXCbqUvBkVvF9ZrVBvN1ralfDaeoNQboccwxMMR3bifgxOx75x37tTKzO6Y9T6GfRe3c77UT1MnkkGIetrF3I/OdqfFmx+m1546K8b5ggdF044gO5FbKZ+hUuvFiklFWe2JwBHdcqzPoUHJXlgjSofe+blPGr3aMEmL+igtlkGQ9HsZtCWTxN1d6oS9lM8yV8pxaslfdB0R2THn9tijR4oVY8Ukzqk+rZjW2EqgsFvXxLgNng5rAeMkmCRBdzNu4t4vqRN2sCKyxYAMicWVjV5rKRfTirG3DRugNUSnZcXEHDwtOlbMXPid2uJyhW0zhWEVABsLkhSYSRV1hzTciQxrcDyV35Cdc6VA0e30Wkv5TCwLwMncMFZMcrSyYmJ57M2Wx74wlaeQywR+JqeWK8aGcWFEvTeSeq9pjtpTKexhqW+nlitMF7Kt7InYg6f2tufWalRjRJgGfxYdjz3G+azUGpTs6FtE2DFbbL2OF6eWKxOd6jgK0hDFjoK0i3u/9l0qhX3nXJFTy2X8Ss90Rm7xhb0tPkGWjyEa69UGy/YatHE7WCdiB2vgPKhDX1wqj52wD8pfT0JQNqOguxnV+x+GHZNOYZ8vUq41W2LRyamOL3gxn4nn7boGWs0Aav+4z2EsK8Y1eAqwY7boa8WoKosrlVb2zCQR94tuRD05krhjSaMlk05hd1LffG7LF5c3fsFLuSzVepNmM9oHVK41mLIjxaDsG0M0HDGeymfjpTvWo0fs59Zq1Bo6sR77sNis1ksYaRT3fu72UirswRkSix1eqyMOUaezl2tN9m+zFugwk5T6xzmH+7dNx0o77bJi5kqcX6tR8XgN51pw5jmMA1G+mMOM1o2gBzNMcR+0HZO6mafgmoXocVvu+LnuL7hzO1+uNZgqZLv26aRca3Dx1inuP7UcS9gvrNf40N1H+dFnHCATUHlynFFV/vYLR7j+6p2RZ3g6dz2XbJvmy4+ej3yccq3ZGjyFdod+eqXKxVumNh6jVU5geFaMiFwFvM/10OXAG1X1rUNrhItxE/WFB/sPmi5cMdyO3DlPaZmt2iupjNh3L1hf6qPnuldSOnrOKvu+Z8FlxdhRX9RosVJvMlPMMV/KcyGkiqSbTx4+yZs++lXui1HBcNx49Pw6r//gl/novccj73PeLre7e6EY2Ypx7q6Kroh9jy3mR892l/Z3rgX35z5oVPVrTuVS4KnAGvChoTXARa+iPmzrZeHBSusnja8XlV7PWVqi9lQK+2wxx+75Eg+eWul67gH7sSt3zLUea0fsUa0YK82ukMtQbcQvRfCAR7smha/b7y1Odku13iSXEWYKucj7OSUg3FbMFTtmAHhg0ftzL+UzXZH8ELkeeFBVH07qBQd9Oz4sQR+W+A5b5Ich7mH06rOn0ooBeMyuWd8vOMAVO2dajzn1RqKKiuPtFrKZWGVmnZz3SRZ2pzONk99frTcp5DIU81kq9Saqir3wii/O3ZU7K2bvwhTThazn+X3g1ApX7JgdpQX2EuC9Xk+IyM3AzQALe6yOJ8k0x16EYhiiPuwo2u/Yg7Rr5h7WsbRlUhmxA1yxY5YHTq105bI/sLjCxVummC60+6SWFRNZ2K00u2IuE1vAoB3VTiJfP2kLe4w7mWrDEnZHpKN0ls5n5XTKAJmMtD73Th44tcKVO72rPg4aESkANwF/7/W8qt6iqgdV9eDM1vSlviXNKKyRIAbdlnEcdO53zdPfFJF7ReQeEfmEiOxNqmFX7pxlrdrg2IWNA6hfP9n9BS/GsGJUtZVmV4gr7LbYeVlEk4JzlxQ7Ys9mYt05lT2sGLA+905hX63UefT8Olf6lPMdAi8A7lbVk6NqQBwGJURpE3Q3aWtb1LusQdlx/Ubsv6+qT7QHl24F3phAmwB4jC3e7i95s6l84/RK6zmHOIOn1UYTVWuffDaex+5Eog+fXZvIUgSq2jrfcYU9n8247pxiROz5jZfglTtnOX6hzIprcto3FlcBy54bES/Fx4YZNGmY/JI20QxiUG0dt6i93zVPl1z/zgCJvfsrPYT90fPrlGvNrojdiRSjzHp0RKfoDJ72YMU0mspDZ1Yj7zcunF6pthaUjnNeKo0mRZcVEy1id4S9O2KHjXdFDywub3humIjINPBc4INJvu6gIrWkBWhcBL2TceqMBkHfHruI/LaIHAFeRoIR+/bZIlun8xuEvZUR0xWxR7diKi5BKWR7E3Zoe9GTxNdPtRf3juWx1x2PPfqdk1MCwitiB7o+91xGuHT7DMNGVddUdbuqTm6OqweTIoyjfA+jvNsKFXYRuV1E7vP4eRGAqr5BVS8B3g28OuB1bhaRO0XkzsXFaMtDXblzlvtPtsXmaye9I7dijMFTt7cbN92x2mgwXcgiAg96ZOz4cd+jF3jJLZ9hqVyLvE+/VOoNfvSdn+c/HzgdeZ8HbctjtpjrKSsmTgfrfFbF3MaI/dJt0+SzsuFzv//kCpdunyafTe1Y/0BIgw0zCSQl7qOyY3rJsAr9pqjqDap6jcfPRzo2fQ/wXwNep5U5sGNHtNvQp122nXuOnOe0XYHxXw6d5Kpdc2yZ3njBO7MXIwm7K80urhVTqyvzpTwzhRxL69FF+t6jF/jsN87ysftORN6nX/7t/tN8+v5Fvngk2kxQoPWedi+U4kfssQdPvT32XDbDUy/dyu2HTqKqrFcb/McDp3na5dsjt2ezkpTwTEKkvtnpNyvmMa5/bwIO99ecjdz4pD00mso/33eC4xfW+cJD57jxiXu6tmtbAPHS7OJH7FZkms9KTOGzjnlrjNmc/XLrvcfsY8e3mmYK2XgdXqOdxw7RhL3SGuvoLgFx4xP38uDiKodPLPPJw6dYqzY8P/e0Mk5L4XUyqaI+qvc1qruufu9tf8+2Ze4Fnge8JoE2tbhq1xxX7pzl1i8d4x9tUbzxSd0ZlXHy2N1WTLEHj90S9gy1GMJea1iR1H88cJozQ6j/Xq41uO2rJ+1jxxPofFYo5uIJe2ceeyQrpu49eArwgmt2k80IH/3SMW699xg75oo87bLNFbEbGyZ5khD3QdgxgxhI72vmqar6Wi9JICLc+MQ9vO1fvs7hE8tcc/E8l13UPYCWzQj5rMROs4trxVRsyyGfzVCtR/+Anei+0VQ+/pWT/PDT9kfa7yvHLvDGj3yFRlN5/QuujmxH3PG1RVar1vuML+zWeVmPWVKg4Ep39KrO2ImfFQPWwPkzr9jOX3/mYSr1Bj983X6yE1p0LSmSEJxJjdY3I6ktKeDw0uv2c/j4sjUY+IwDvttFXdDanWbXixWTt9Mk4womWB3QkXPdBa78uPOhc9z18DkA/vPBM5GF/YhdRCuXkZiWipLPWlbThfU+smL6mKDk8DPPeQz/91MPkskIr3jmgchtGQeGtaBxHDaLqC88WOm7BEHcMgPFRwpU9kcvNpgEqRf2XfMl3v7yp4ZuZ9UpiZdmFz/dsUHRFr64wp7NCMVchlpM79rr79B22ttOFbJUG/HuLJyIPfadTC7jGsSOdudk3Wl5u4HXXbaN6y7bFrkNk0RcG2bcJs+MmiTEPe1MTP5YKZ+JnWbXywSlXj32fNYSsXrEVZ6c/cCK9OPsV7f3my5k47Wz3qSQFQq5bOw7mWIPEbu7FrthdGyWaH0zMTHfrKgLWrsnKDlCG3VJvXZWTCZeJGxPuY9bwsAR5el8/CyVbEYo5WMKu2M19TCoHLukQMeyeIbRsFlFvd/3Hfcuqd/B8LiZVhMk7JmY3q5lOUD0WZbOIGEhG99SKTgWTsz9MmIVOYsr0Dnb5oh/Z2Gdl7jljAvZTHsQO+LgqRH2/unHhtmsou6Qpvef9LjL5Ah7LtpCyu7B02IueplZcFkxufgeuxOxx7ViWvvFuEOoNbTn7J181ilnHCMrxr6TgeiD2JVas1WV09DGpDkOlzSJe5JMzDerlM9GrFHSHrRrRewRhb3W0N499pzEntjkdAi5HgZrLUsl/n6Wx55p+fthNJpKo6mtc1nMR+9gSx6Tk8adpXrw0n1pyYiZVEEbJmketJ4gYY86eNoetCvYGRlRxa9S79Fjd0Xsca0YZ9A1bofQmxVjZ8XEOJ7z+q2IPZ+JVmWz3vDMYTdEp1dhMaK+kUk8H6lPd4xKMZ+NWLa37e3Gjdir9UbbY4+dbdK7FVPowYppdSSx2tn22J1IPGxikGNjOZ1k5DunWtN47BNO4fDRRF6nevW+RF4niF5TINO6dN7EhEzRJyg1u4U96uCpq1ZMbx57b/v1YsVYYwE93FnYdyQQrcNztim6IvbIVowR9p5Jc7ReOHw0MVF3Xm8SGOb4yeQIez4TrQhYvdEatCvEEDBnu0JPloqVx56LmUbYrxVT6CELx/HYIaKwd1oxMWYAj4MVIyJbROT9InJYRA6JyDMGdaxBf/GHJeqDet1BC/yoLZmw8Zc4KY/p/2ZFJE4euzNoV4iRFVNvNGkqvUfCjqUSc6JRr1kxPVkxrpmnAJVGhKUG650ee4w7p/EYPH0b8DFVvRp4EnBoxO1JJcMQXuc4g6QXcU/jIOrEeOxTtqCoKiL+npdlxdgRe4+RaWyPvdFktpgjlxGWyvFmdOZsC6cSwd5wHy/fa/ZONkOxByumkLVEupTPcnY1vC5Gpd5olflNKyIyD3wH8EoAVa0CPRf9SCojphchGVRhkPWzAAAgAElEQVQ0OgqbxDnmMLz3cWViIvbpYpamhs96XKnUmSla/VkxhsdedQ0S9u6xx7diCtneslva+8WfIRurw+uI2GeKWVar9aBdAOtzmC2mW9iBy4FF4C9F5Isi8g4R6Sov6l4drHY+epG3YTFJoj6M44/akkmCiRH2OVus3avbe7FaqTNXsrbtZZCwpzz2em+1YuobLJV4++UyGQq5+HnzhZzEGlSu2nZNPmvdJc0Wc6yGfAb1RpNyrclsMR+5bSMiBzwF+DNVfTKwCvxy50bu1cHyW6aH3caRMGpRd0hLO9LGxAj7TERhX6nUmSlY28aJTCtdwq6oRhNbSzCzsSN9y4qRnvbL92gZOWMBEP+8gCXsy+WwztXqDGbSH7EfBY6q6ufs/9+PJfSGFDEIcY8btUe1x4aVGTMxwj5rC3tYtLhSqTNrR+wtAYsySNhop/U5IhY1iq66slviT1DqbRC00NPxtGcrpugS9kq9GdjeFduqce6c0oqqngCOiMhV9kPXA18dYZNiMwhbwUTJ6WfihD0oWlRVViv11rZxBKw1w9L22N2PRdm3kM2Qy2aoxbRiCn1YMfkYpQFgYz125//QfToGT2cidLAr9mfkbJtyfgZ4t73847XA74yyMaPOwEirqKe1XaMiEWEXkV8QERWRi5J4vV5wovAgQVmvNWgq3cIeQfwcAXMiaIgj7M4M0viDrr1YMe6smGqjGckyUtUNkT5E7fCs125ZMaVwS8x5bnYMhF1V77H98yeq6veq6rlBHMcU/+qfpMV9nAdR+xZ2EbkEeC7wSP/N6Z0oHntnpFi0o8xeBk8h+ozVWr23WjHVHq0Yx/op2HcWUQZsG01FlVZ1R4iaBmrZWG6PHSZH2A1tTFQ8PiQRsf8h8DpgpPeIUbJinOccb7fXtL528bAYHnvOmnkaxxpp13GPv18+E+/Ownn9fK5HK6ZD2IPunJznZlPusY87SUac4yLq49LOQdOXsIvITcCjqvqlCNu2cn0XFxf7OawnkSJ2+7mesmJcE5TyOdtjj1zu11mgQ6g1o1kjYHvlWYldK6beKhNsC3uEmuyOiMfNinHn90P7cwga62jdORWMsI8D4yaW49ZeN0lNYgsVdhG5XUTu8/h5EfAG4I1RDuTO9d2xI/ma1NOFLCIhg3YdkWI2I2QzEi0rZsMEpeiRcKOpNLXtzataj0XBnX4Y34qxBk+d/6McC4hdK6Yz3XGuNdbhf04775wMhqRJStzj3PWMemDbTeg3S1Vv8HpcRL4FuAz4kj2Ffx9wt4hcZ6eJDRURYbYQnEPtRIpubzfq+p7utL44HnvNFQm3H1OilElxr5XatDuEsDK6zjGdwVp3G6K2M94EpY3pju07p5rvPq07J+OxD4ykbJhxjn43Mz1bMar6ZVXdqaoHVPUA1mSOp4xC1B1mS7nYg3aFXDxhj+uxty0OaadJNqNF3/WmUxUyXnplvdGe6Rp1P8eu6TWP3TlWlLTTlUp9QwdpiMawI8JxF/Vxb38/TNQ3K2w6u9egXdSSuO4iYPEE09UhOBObYnjzbs87yvGaTbU7hHjtbHVAro4rStXLar3ZsrTAPXgabMUYG8YwDDaruCcm7Hbkfjqp1+uFmWJwxL7sEbEXc5nIAgYdE5Ri5HnnsxlymeiRvpVX3inQ4fs5dwPu/aIsaL3BY485eFpwRd7ZjDCVzwZbMeX6prNh0rLWaVQ2qyB6MY757BMVsc+FWDGrlTq5jLT8YIhvxcRNB3R713FmrDq5524rph5lv0Z7v0IuvsdeyGXIZIRcJloBMWdVKTeWJeYfsbtn/xqSZxyFaDMxjMloEyXsM4Vca4DUCydSdNdrj5px0rJiYkbQbo+9XWMmbofQa0fS2/HAEvhodyQewh7hzmmzReyG0bEZ7z4mSthnS8Ee+0ql0RUpxh487bVeeUwrZsNgZg8dSS5mh1Cpdwt71P0K2W5hDxvrmDPCnlo2oxAmRVpSHifq2zVbzLV8dC9WKjVvYY8YseezQiYTM9vEFQk37YlJkfZrtiP9XqyYwoasmCgee7sjsfaP3uEVOyL2mWI2+M7JROwGw0CZrIjdjhT9ZnauVhpd09jjCJgTmTpeeSxrJNeeMDQsK6YV6UexVDpmkMa5k+m2YvKhYx2mnECbKJ5rWiJBw3gwUcI+U8zRVKuKoxde3m4vAhYn/bBabw9m5odgxdTcVkwPg6fOPoVcplVGIQjPwdNiNthjL5vBU8Nw6ddeGrcB6YkS9lbJWB8bYKVc6/J2CzHSHR0By8eJhD3quEexVBwrJhfTiql5WDFx8vR7sWK6PPaA7KRao0ml3jTCPiD6FSDjr08GE/XtchZHXqnU2enx/Gql0bUcWxyPvSXsMVZQclsjmYy2XivqfoUkrJgYHruzT7EPKyZoPsHqmJUTEJGHgGWgAdRV9eBoW2QwhDNZEbu9OLKfqKxU6l0LKBeH5bFnMz1bMfEGQbuzYnpNd4xWj71b2OeKOar1puf+TqmBMcuK+S+qeq0R9fFmWHcjaRgPmShhn3FF7J00m8pqtd6K6h3yEYW9YhfkAlwCHUX4nBWGJJbn3U5bdNWYidQBuTz9noqAif074p2M67w4BC2Pt1odr4h9M2FsmMlhooR9zonYPTz2tVoD1e7FHeJYMU5anzMzM06tmLgThuoeVkw9QvEwZ5uCu2xvzJm10G9WjH9t/FaFzfHJilHgEyJyl4jc7LWBe62B2vm1ITfPMGkkUX5ibL5dUQiK2P283agCVusQsKirGrktDsHJY4+RV76hTHA8K6Yfjz3y4GmjSdFjghJ4V3hsV9iMULc4HXybqh4TkZ3AbSJyWFU/7d5AVW8BbgGYu2r36O/DDZueiYrYt0xb+cAX1rsLUJ1fsx5bmNrosZfyVlZM2KpG5XqDoquIej4rERd7dkXsPaQf5jL9WDG9e+ylfJZKPXwBknKtSTG/8TJamLbOsdfn4DzW+TmkFVU9Zv8+BXwIuG60LTIYwpksYZ/KkxE4s1Lteu7MipUGdtFsccPjU/ksjaaGRrXr1QZThbawF3JRa8y0I+GWpdLjBKW4Vkw2I2QkurBnhFb53al81nc+gJtyrcFUfuNdkHOOz6x2p96dtj+bzs8hjYjIjIjMOX8DzwPuG22rDIZwJsqKyWSEbTNFb0FZdQRl4yy/Ut4S6/Vao8srdmMJmDtijybs7qqJ2DoZzVJxBl17s2Li1nxxltNzmCpkWa8GC7uqsl5rMFXYeN62z1jn2K+DzWWE+dJYROy7gA/ZReNywHtU9WOjbZLBEM5ECTtYwn3aQ1BOL1tiv32mI2K3o/ByrRFoD6x7CnuUtMV2tokyeCvGSZPMubJboixmXavrholGpXyWci34eLWG0mjqhvMCliVm3Tl1d7BnVqpsmymQibDE36hR1W8ATxp1OwyGuEyUFQOwfbbgLSirFbIZ6RJvR5TCotNOKyafjVavvNZoIrbFMQwrpuayYpzfUTuSvOuOZSqfpdpoBrbVsWpKHcKezQjbZgqtuyQ3Z1YrbB8DG8ZgGGcmT9hnipzxEhSfSHG60LZigijXmh3CHq1eedVeBUnESpF0HgvDy4qJd4fQLn8QWdiz7XPjnJdywHss2+dsutB947d9pujZwZ5eqXbZYQaDIVn6EnYReZOIPCoi99g/L0yqYb1iReweVsxKteX9unF77H7UG02qjeYGyyHq4Gmt0Z6xKmLZKr1aMdGycDqsmFz0lZDcHnupEH4n4zzX6bGD/+dwZrUyFgOnBsM4k0TE/of2dOtrVfWfEni9vrhotshKpd6KJh38BMUR63KQgNmv1ZPH3hEJ57OZeFZMrh3tx7FiNkbs0e4Q3B57FIvK67w4bJ/1v3Py6mANBkNyTKAVY4nG6Q4b4MxKle0eFsBUBCum5SX36LG7I+E4Qgu46tNEHaztXjAjaj32vJewRzkvXsI+U2gNWDusVeusVRvGYzcYBkwSwv5qEblXRN4pIlv9NnJPu15cXEzgsN60cqg7bIAzK8ER+1pAZFquNjdsC9G962pdO4Q9eocAtHz5XIwJUe589Fgee659Z+HYK4HCXvWP2HfMFVnuuHNyPhOvDtZgMCRHqLCLyO0icp/Hz4uAPwOuAK4FjgP/x+91VPUWVT2oqgd37Oi/FoIfjmi4c9nXqw1Wqw1PQYnisXtZDnGyTTpLEURb4q6dTeMcL6oV00tH0uWxR7FiWh67d8QOcNZlx5xuTRIzwm4wDJLQPHZVvSHKC4nInwO39t2iPnGicncuuyPyF814ROyuPHY/WsJe6LBUouSHe3js0RalbmfTxDpeRz56vKyYbism0nnx8djBitL3bplq/Q3dcwkMBkOy9JsVs8f17/eRgunWrYh9xR0p+lsA03GyP1xT5/MxsmLcgpmL4827UjNzMbJpcq6OxMreiT946qQwRrqT8YrY7XN92nXn1Opg54ywGwyDpN+Zp28WkWuxSps+BPyPvlvUJ9OFHFP57IYcaudvr0G7Ui7ciil7CFh0i2Ojxx55MLNjwlAhYn1078Ha7mJcXvvNu0rpRsmKKQdE7M7dkWcHa7JiDIaB0pewq+rLk2pIkmyfLWxItTsTICiZjFDMZQKFfc1jkDCyx96xJmg+m6HejBZBdwp0PWLk3eWxR6yrvjGP3fp7LcrgaUDEvrGDrTJbzHlm0RgMhuSYuHRHsCJzd7qjYwf4ZWNMFbKDzWPP9Wap9GrFdHv6vZUUgGj5/aVct1BPF7KU8pmNHexqxWTEGAxDYCKFfc98iUfPrbf+f/TcOvOlnOfUdwgvUdvOY3dZI7kMlQhlbTuzTaIuxdcptFGXquu0Ygo5q958lHZ2FgGDcI+9mMt4FvQSEXbPlzh6rr2i0KPn1tk1Xwpti8Fg6I+JFPards/x0JnVllVw6PgSV++e993eEvaAmigeVkwxomCWa40NEW0hohVT9/Dme7FiirlsxHY2KeU3diT5rASPPXQURuvkqt1zHD6+DFhrzh4+sczVu+dC22IwGPpjIoX9cXvmaSrcf3KZZlP52ollHrfHX1BK+Szr1e5l3By8ZliW8lnqTQ3NSe8UzKiWSrXRbE1OirNfpxVTymcCUxbb7dy4QpS1b3BN9rVqw3Pg1OFxe+b5pt3BPnp+nZVKncft8e9gDQZDMkxcPXagJeKHji+xdbrAarXB1QGCMlUIt2LcS80BLbEu15vMZv37x3KtsaFDiGPFdE5sWvVYy9Vrv66IPaSuOkCl1uwa1JzKZ0Pz2IOE/erd86jC104uc3KpbD82fhG7iGSBO4FHVfXGUbfHYAhjIoX9kq3TzBSyHD6x3FoHNShSnC5kA0Vz3SMydaLbSq3RWrzZi0q9SbEjbbEXKybqYG3nfqW85c03m+q7uEWzqdai1B0rSE2HdHjlWrAV83j7nB8+vsTJpQoilj0zhrwGOASY2w3DWDCRVkwmI1y1e46vHl/i8IklROCxu2Z9ty+FeeweAuaO2IPojNh7tWKilvutdkxQco4d5LM7z3VG7GFWTFjEvm/rFDOFLIeOL3Ho+BIHts/4DmCnFRHZB3w38I5Rt8VgiMpECjtYEfqh40t89dgSl4UISpjl4OUllyJMuVdVK2LvTJPs0YqJW/8daEXhQe0st8YQNl4OoRZVyOBpJiNcvWeeQ8eXOXRiKXCcI8W8FXgd4Hvy3QXuaufX/DYzGCJxYF//RRInVtiv3jPPcrnOHV9b5OoQQZmKEJl2RrOOYAb5104kXOxKW4xmqWyM2KNbMb1G7J2Dp+HnpduX7+Tq3XN88cg5Hj6zFpiZlEZE5EbglKreFbSdu8Bdfsv0kFqXPNWr9426CYaEGK/74hh81+N3cddDZ6k2mvzoMw4EbhsWmXpZMU4UXq777+eIfqkjTTJaud/mBqEt5KLlsVteuTt7p4+IPZ/lwrp/OYJyiBUD8NLr9nN+vUZGhO950t7Q9qeMbwNuslcGKwHzIvIuVf2REbfL0APD6riWLx39Qu0TK+w750u89SVPjrRtKWyCkpcVkwu3YhzRdwtmMRcx/bDeoNjLfvakoa52BnRA7XZ2vMcoVkyIsF9z8QJ/8sNPCW13GlHV1wOvBxCRZwO/YETdMA5MrBUTh6l8lmq9ScMnW8VrkNAR3UArptZtcTgTm1SDbZVKbWOWSjEfbUKU5elv3A+sfHo/yrVuywjssYewwdMAj90wfC5cYSpnGoywA+06634R8XqtsWFZPGhHwpVIkbBbaK39wmyVSr3R0SFYnU94h7Bxv5IrLTNoH6udHh57SB77Zinopap3bIYcduOzTwZG2IEpO2PGb3m8crXBdFdWTJRIuLtIVmvQNST6rtQ3zlh1/o67X3ssIKCdrXTH6HnsjaZSrTdb9ewNvVPZ373odydp8G3HmX47rHG7EzLCTvhqQWueeezhHrtXfrgjtGGzQSsdg6etCVEBAl1vNKk3tcv6CWun85xXSYFyzZrc1EnQ6kkGg2G0GGHHtaiEnxXjOfM0PIJuCWa+O688yMKpNyy/f4PHHmE/x97ZMHgaa4JSdx47eA+8OmmQnRaVYfwxdsz4Y4SdtsfulbPdbKptb8SP2B2bJq4V41gjnh1CpMFaj5o2PUTsQasoBa2eZDCkiWF2VGmxzIywE1x73IlWu/LYc+Eee8Uz3THciql4CG3LwgnKm291CB4WTo+Dp+B9XowVk17GzQ/ebEQZU+mXvoVdRH5GRL4mIl8RkTcn0ahhEyhgHrXYAXLZDLmMBOeHe6U7tmrMRBBoDysmymCtd8QeId2xw4pxbBavaL+9LJ6JDaKQxDTxYWLsmPGmrwlKIvJfgBcBT1TViojsTKZZw6XlJXtYDkGRaSkfXBLXa0ZnrFIEnt58vNIApUiRvvcSd20rpvuYXjXqDYa0kUQHNY53QP2GWz8J/J6qVgBU9VT/TRo+03mrf1vxKN27FrBgcymfCYnYncHT/oU21n7uhT0yQkbCI/aMsGGBDqCVyuh9Xur2NhM7eTl1DNu/Hfeofdzb3w/9CvtjgWeJyOdE5FMi8q1+G7or4C0upuu2dH7KEqflcreALdm1Uuan8l3PFXPBVSG9sk0iRd4e1kivEbuI2GmLwR1QKZ9FZKNwzJes97xc7q4Xs7Ret7cxwp5GkooyN7M4jjOh30oRuR3Y7fHUG+z9twJPB74V+DsRuVw9pkeq6i3ALQAHDx4ML1M4ROZsAVvyEjD7sQUvYQ+Z5l+pNRCho4xub9UWow26+njl+eB1TzsXA3Fw3vOSV4cXcF4MhlGTVIc0jjYMRBB2Vb3B7zkR+Ungg7aQf15EmsBFQLpC8hCyGWGumGtFoW6CItNSLhuYbVK2BdMdCbc99ijph901X6JYMZ0iHVZArHMxEAfnTmbJo8Kj85jTKRoml+rV+ygcPjrqZkRmVHcZSVhlSQ2y92vFfBh4DoCIPBYoAKf7bdQomJ/Ke5aovRBgxVgLRQdH7F113COUBvCO2OPUf/eYQRqSN+8l7M6Sf37nZSqf3bAYiGFyMZbMeNHvt/KdwOUich/wt8ArvGyYcWCulPO2Yhxh94hMwzz2cq3b4ohmxXjNWI0+eNpLxO5lxeSyGWaLfuel3oroDelkXG2EfjAdkEVf30xVrQITUZ96YSrvbTmU/SPTUj7jmTHiUK77r7wUOOjqMYM0SqTvl49ejOKx+6QtWufF22M3/npyVPZXKT5SGHUzAkm7JZO0qI9zx2juo23mp/Leg4QBkWlYtkml1uzKDY+T3eLuFEpRIn2PapLW/xE8dh9LxfdOplzzvIsxDJZRT1lPa0Sc1naNCiPsNvMl74j9wrq/gBVzIXnsHasggZV+aC22Ec9SyWcFkZDSAB4Tm8CZSBVcUsAvYg8ae/Aad5gkRKQkIp8XkS/ZM6t/fdRtissgok4jounHCLvN/FTO14rxE7AoM087o2ewV1GKOQja7hDC93OnV1rtDN/PL2L36/CW1uubIYe9AjxHVZ8EXAs8X0SePuI2GToYREcTt0Mc9Z1UJ0bYbRam8ixX6l3L4wV5yaFWTMcydQ5hnrfzmp2+fviEqAa5jJDLxtvPL90R7PPik8c+6R67WqzY/+btn7FMDkiatETtaWlHVIZRAAyMsLfwm2VpWTHekWkxnwlOI6x5pxGGWzFN8lkhm9kYBYRG7B5ZOBCelmm10ydin8p1WTHNprK0CawYABHJisg9wCngNlX9nMc2rVnVtfNrw29kCIMaBBy1qI76+GnGCLuNI1KdGSDW4Kmfx26tQ+q1whA46496ROyRBNqjQ4hgqXh55VYee8hYgMfxwOrwVip16q41WlerdZrqnQI6aahqQ1WvBfYB14nINR7b3KKqB1X1YH7LtO9rJTX5JE23/dWr941EYAd5zHHOhnEwwm7Tnj7fjk6bTWU50IqxTp/fwtReC3RY+4UMZtYbnhF0KZcNXjzbJ7sl1NMPiNid9+5O63SyhybdinGjqueBO4Dnj7gpPTFosRqWuA+6I+nlPKWpo3Uwwm7j2C3ugcKwyNQZGPXzry3vuoeIvR4QsYcMugZF7F5zx1TVM9/ewetOpl0YbbIHT0Vkh4hssf+eAm4ADg/qeIP2X4ch7oMS3WHcGYw6Uk+yZr8RdhtHwNx+8oUQASuGLGJhzej08tiDs2n8inIVc2ETjbytn1I+i6r3nUW10US1e7aqg9PheZ6Xybdi9gD/KiL3Al/A8thvHXGbgN6jxGGIV9ICbLx0i+fujh5TTHbIFYN5DyumXQAsOGL3s0esGixeWTEZzq36R2flWsNzpmtYaQD/DiHjej7btQ/4L5jhfV786+dMEqp6L/DkUbdjHBk3Me61w4vTwQ4rIwZMxN5iwctyCClN217QujsSrjeaNJrqn8fewyBo74Ou/paR12IgbtrnxSXsm9BjnyRGbTmkjUk8H0bYbWYKWTLiZ8UED556CmZAJBxqqfhm0wQPnlY8ZroCrQFVL/vHecx3glKQRTX5Vkxs5nPlwOeT9FH7GbSbRDEbNmkcNHUwwm4jIna9GA/LwbekQJRI2C9LpQdLJUq6o+d+vUfsrUFlj/MyO/kzT4fOMG/XDZPbwRlhd9E5fd6xHPyLgPkX9Gp51x7WSJQVjTzTJKMMunruF6Wd3pfCTCFHRrotqrlirmsClWH4mKi9d9L0/pO8kwMj7BtY6KjwGLZKUClSJNzLIKiPFZMPnrHqV1c9Sjv9Bk8zGa87Gf9JW4bxIk3iNk7E7VD7vROLkxEDRtg3sGU6z5mVSuv/M6sVFqbyvpGpE7GvewjmetVfMMPXSvUZBO0x/91pg2c7Q4QdYOt0gTMr7QvzzGqFrTNG2A3jy6R3aEbYXVyybZoj59Zb/x85u84l26Z8t58tWuK2WukWzFV7puZcsdvGKeay1Ju6YZq+G79BUKeYl98iVX6Drs4Sd6sei4I4j816tNNh39Ypjpxr10A5cnaNS7b6T5039Efc6K7fQbxJF7lONsP77UvYReR9InKP/fOQXSxpbNm/bZqzq9VWIbAjZ9fYv81fwJzBw5VKd1lbZwq+1wCjI76+pQh8inkVcxmaCnW/2jQ+1STn7DZ4VWl0HpsLGAjdv22aR85awt5sKkfOrQeeF8P4sRnEDpJ5n8O2YXqhL2FX1R9S1WvtIkkfAD6YTLNGgyNWR86u02gqR8+tc0mAgE3ns4jAiodgrgREwmELUweVFHCe70RVffdz2uC1jF9QOx32b5vm/FqNC+s1Ti1XqNabgefFEEzSA2VJsVnEfTOQiBUjIgL8IPDeJF5vVDjC/sjZNU4ulak2moGRaSYjzBZyLHsIphMJe0bsef9l7ppNpdrwLykA3qsoOdG/134zjrB7dUD2YzMhwg7WHYwTuZuIfbAM247ZDIwiWh8VSXnszwJOqurXE3q9kXBJDwI2W8oFRuxzxe5BxnaapL9Ae1eFjJC26LFfIZehmPNeeHulUqeYy3iWMHDo5bxsZuJmMKSJSY3aR/W+onTQg7iDC51hIiK3A7s9nnqDqn7E/vulhETrInIzcDPA/v37YzZzOCxM5VmYyvPI2TUWpi1BDhX2Ys5bMMt1shnxqe7oH7E79kxgxB5zP7A8dM87i0o90F8H2L+9fSezWm2QEdi7xX9Q2TAali8V5h7uf4GnC1cUWXiwEr6hIbWECruq3hD0vIjkgO8HnhryOrcAtwAcPHgwtcuLOQOFW6bzkQRsxk/YK3VmClksl2ojjvh6lyIIzn/33c/Jm/cR9tlizjcrJsiGAWvi1pZpq8NbqzbYszAVGOEbDGkiqWh9XGwYSMaKuQE4rKpHE3itkbN/2zRHzq5x5Owae7dMkc8Gn6K5krewL5frvhObokXe8QZPWwtg++Sj+1pG5XrgwKnDfjsV9JGQTCFDcoyyvMCkWDKjfB9JfX69WHtJCPtLGPNBUzeXbJvm6Ll1vnkmmoDNFv089pqvYLYE2qsoV90/8g4aPA3az2mnnxUTRdgvsTs8I+zJMKjMmCSjygtXFMdW4Me57UnQt7Cr6itV9e1JNCYN7N82TbXR5N6j56MLu6fF0fAtkhVkqbQi74C66l4LaAft12qnT8Qe5rGDdV6+eXqVxeVKy3OfdETkEhH5VxE5JCJfEZHXjLpNo2CcRHJQbR0nGwbMQhtdPO8Ju7j/5DK1RpOXPe3S0O39LI7lSp0tPvVUpgtW5L3mIeyODz5d6P5opuz91qvdx1uzZ7967QcBg7wRI/YfPHgJK+U6GYEXXbs3dPsJoQ78vKreLSJzwF0icpuqfnVYDajsr1J8pDCswwXiCGZaB1bHpfNxM6g7NyPsHVw0W+RNNz0h8vZzxRwr1TrNppJx1ZRZKdfYt9V74HUmaIp/1ckrD5po1N0hOKLttR/YHZCfsEeI2C+7aIbf/N5rQrebJFT1OHDc/ntZRA4BFwNDE/ZeSCo7xo+0Zc0MWtB7idZHXX7ZpDb0yQ3HgQgAABFtSURBVGwph2p39L1SqXvWiYFgYXdE2yuKDuwQQmaQzhbzAYOnpqBXGCJyAGuZvM95PHeziNwpIneunrO+0EnmsvciEoO2DhzLY5RR8jCOP24WjIMR9j7xE9uVsn8a4UzBf4r/WsV/JmirmJeHFdOO9L2POVfKUW00N0yKqtQbVBtNZn2ifIOFiMxilcx4raoudT6vqreo6kFVPTizNbptMujSAsMSpWGK/DCP1ev5SzJa7zVAMMLeJ47YugtsNZrKarXhGz1nM8JUPusTsfsLdDGXIZuRwP38I3anA2oL+2rA3YHBQkTyWKL+blUdWS2kXsVi2BHnIER3FHcHwxD1QXbs5hvdJ3Ol7ujbiZ6Dsk2siU1eg6fWYzOF7ihaRJgpZH3LBGczEpgVA9adxLaZQutvgFmzdqkndg2kvwAOqeofjLo9vQ6kOiI1SN+9k3EcyHQYV/vFjYnY+8Txp93+dUswAyLh2WKWNR9LpZTPkPOZGBWUXuk30xXaxciWXSWGnb9NxO7LtwEvB57jKk/9wqg7R7mNHmalx0kQrEGyfKn0dY5GPWDqxnyj+6SdqdIWzKBa7A7TBe8p/lYpgoD9eiwNMOeK2FvHilCLfTOjqv8OpEoN+01/HHTGzLjSb6cXV9SjdOj9DMCbiL1PvBaxiFLj3D/yDhZov9o0q9Xg/WY9LKMo7TSkj34jw34j00ljEs+FEfY+8VrEIkokPFP098oDBbroN+jaCO0QutoZ4c7C0B9ps2PcTKKgxSGpDm4Q0Xq/GGHvE690x3Yk7D8oOeNjqVgzQf3TD2cKOd8OIWi/uQBh98u3N6SXpPzczSruSb3vQfnq/c6DMMLeJwV7kYplj4jdbxYohAyC9mrhBHjzLSvGw2MPK9trGDy9RHFJivtmEvi0i3oSGGFPgLmOAlvLrUg4fsQexWP3mqAUVvNlKp8lI90Ru0i7do1h/EhSXCZd3NPQgQ3LdjPCngCddViiROyWQDdoNjdmKKxU6swGRN69dggiYpXudXdAdi12vxRJQzKM01J5oxa+QZH0+0pztA5G2BOhsyTuSqXGVD7rm4tu7eNd4XGtGmbFZKk1lGpH6d4wCwdgrpTvitiNv54eeo3mkhaZNES2STJuop5EIGC+1QmwdbrAmdX2h31mtcqW6eDZnO5BV8dCUVVWqyGDp679Cjkrn7lab0aq+bJlOs9ZVzvPrlZZmE5HSVhDfwyivK9bEMct9z2NHdMws59MxJ4AuxdKnLhQbv1/4kKZPQulwH280iTXqg1UgwczvdIWHWsmLGLfs1DiuKudxy+U2RvSTkMyjJMd44UTxadRMGFj+wbZxrRbMA5G2BNg70KJU8tl6g3LHjl+ocyesEWwC91pklEE2qvCY1DhMDd7FqY4fmG99f/xC+vs2WKEPU30E9UNS3TSIvLDbscwzm9SAUBfVoyIXAu8HShhrTbzU6r6+SQaNk7s2TJFU+HkcoW9CyWOnV/nhsftDNwncMJQhIh9Q4dQjTaDdM+WEufXaqxXLV///FqNPQvBHZBhvBj2ikudojpIy2bUHUk/DHsSWr8e+5uBX1fVf7aLI70ZeHbfrRozHNvl+Pl1pvNZKvUmu0ME08mY8SqjG5R+6FR9XNmwn7OcXrDH7rTzmCtq32si9tRxYN8iDx3dMepm9ETSvnxaxLyfaH0UM4v7FXYF5u2/F4Bjfb7eWOJEvcculFvrkoZ518EzVuNF7EGrLnm18/j5ts++e95E7MPiubsPc9uJqyNt24+4OyI06rVS0yLK/TIuvrqbfoX9tcDHReQtWH79M/02FJGbgZsB9u/f3+dh04XjU5+4YEXs1mPBguk1eBrHY+9l8HSvI+wmYt8UpGkh7HEkCUGPE60nOcAeKuwicjuw2+OpNwDXAz+nqh8QkR/EWpTgBq/XUdVbgFsADh48OF65UyHMl/LMFnMcO19mKh8vYnfXZA9b3s79nNega1jEvmvBWvzAnRmz22TFpJYkLJm0RO/jxrBFPWlChV1VPYUaQET+GniN/e/fA+9IqF1jh5VKuM50IUsuI2yfDV5Bxons3V55NCvG8ebjR+zFXJaLZgutiP2i2QLFnCkn4IeIvBO4ETilqtck8Zpx7BhIzm830Xs0xtF28aLfdMdjwHfafz8H+Hqfrze27NkyxfELZY5fKLNrvkQ2E+wvZjLOMneuPHZnWbyAiUbFXJZ8VlitugZPq+H7tdq5MMWx82WOnS+bjJhw/gp4/qgbkRSV/dWJEa5BkOS5iRutJz3PoV+P/b8DbxORHFDG9tA3I3sXSnz12BJT+Wxk37qz7ksrHz2gVozffvmsRIq+9yyUeOjMKgAHts9EaudmRVU/LSIHRt2OpLNkTPS+kaQ7u1FaMA59Reyq+u+q+lRVfZKqPk1V70qqYePGnoUpTq9UePjMWuRIuLME72qlznQhSyYk2p8pdO8XtfTu3i1WxH78fJm9IQO8hsHQS3SWtFiY6H0w5yANog6mVkxiHLhoGoATS2UObJ+OtM9sybvaYhhzPe4HcOn26VansH9btHYagnFnfC3sGVxnOYj8drewbZYoflAdWq+iPohyE0bYE+KF37KH7TNFas0m1x3YFmmfHbNFjrkyVE4tl9kxFzzoCrBjrsip5Urs/QBeet1+Lt8xC8DTLovWTkMw7oyvi5+wJVLGV9xBVAdHPAYxgWmSM2gGfXeSJlEHI+yJkc9m+PbHXBRrn10LJb545Hzr/xNLlUhFuXbNl/j6ydPt/S6UeeyuuUjHLOWzfOdjx3NW46TRq7jDYGendorgOAr9sGymtFgvnRhhHyG750ucXa1SqTco5rKcWirz5P1bIu23uFKh0VSyGeHUUoVnPcaIddKIyHuxSmRcJCJHgV9T1b8YbavaDKv0gJ9IpkHwRzlO0K+oD7LipxH2EbJr3rJPTi1V2Dlf5MxqlV1zUSL2Io2mcmalwkwxx3Klzq55M9EoaVT1pYM+Rj9ROwzWmgljmIKftoHeNIs6GGEfKY4Yn1xyzwQN98qd/U4slVvZMFH2Mwyerdm12Pv0K+6QrsJhaRPhJBml9fLi+bv5lYjbmnrsI8SZzn9iqdwS9yiRd2u/C2VOXoi+n2E4vHj+7pEc98C+xdR6vpNAUue2l2g97jVlhH2E7G5F7BVOLllZLlFqt7T2W65wcrm84THDeJLkrbkR92RJssMchqiDsWJGysJUnkIus8GKieKxb58tks0IJy+UWbGtGBOxp4sXz9/N+5eeEmufJCwZB7cQpcWiGSfS0jn2evdnhH2EiAi759vrpRZymdBFsAGyGWHHbJETS2VmiznmirnIM08Nw2PU4u5gRD4agxbzYa57a9RgxOyeL3FiqW2niERbnGDXQomTS2VWCjl2mdK7qSUt4u5gRH4jw4rMh2XBOBhhHzE754vc9+gFhHb6YxR2zRV56MwqM8VcrP0M48Egxd1hs4r8MG2WXqP0fgfgjbCPmN3zJW4/dBKAay5eiL7fQonPfuMMs8UcT79i+6CaZ0iAXqJ2aIvCoAUeJl/kh+2ZD9N28cII+4jZvVCiXGvy8Nk1bnjcrsj77ZovsVSus1ypm4yYMaBXcYfhCjz4i+A4CH4aBj37FfUk0mWNsI+Y5z1+N19+9AKNpvK9T7448n4vuGY3XzuxTFOVm67dO8AWGtLCsAW+kzDRHJbwp0G8vUgiSk9qDoQR9hGzf/s0b3vJk2Pvd/mOWf7opfH3M4yOfqJ2N8Pw33shrYI7aEZtu3hhJigZDEMkqYjsubsPp1JQNhNJfwZJzlg2wm4wjDFG4EdD0uc86TIUfQm7iDxJRD4jIl8WkY+KyHxSDTMYJpVB1JIx4j4cBtGRDuJ66Ddifwfwy6r6LcCHgF/sv0kGw+QzKHE3Aj8Yxu3c9jt4ehXwafvv24CPA7/a52saDIY+cAtQGgdZx4VhCPmgKoH2K+z3ATcBHwF+ALjEb0P3gr/79+/v87AGw3AQkecDbwOywDtU9feSeu2ksmSC8BInI/beDDsiH2R551BhF5Hbgd0eT70B+HHgj0TkjcA/AL4V9t0L/h48eDDSgr8GwygRkSzwJ8BzgaPAF0TkH1T1q0kdYxji3kmngG1WoR8nayUuocKuqjeEbPI8ABF5LPDdSTTKYEgJ1wEPqOo3AETkb4EXAYkJexrYLNZNmoR80IuxiGrvwbOI7FTVUyKSAf4KuENV3xlhv0XgYZ+nLwJO99yoZDFt8SbtbblUVfueBikiLwaer6qvsv9/OfA0VX11x3YtmxG4BsuiHCVp+XxMOzaSRDsiXdv9euwvFZGftv/+IPCXUXYKapiI3KmqB/tsVyKYtnizidriVUO5KxJy24xpODdpaINpx2jb0Zewq+rbsAaWDIZJ5CgbEwL2AcdG1BaDITJm5qnB4M8XgMeIyGUiUgBegpUkYDCkmjQWAbtl1A1wYdrizaZoi6rWReTVWPMzssA7VfUro2pPDNLQBjDt6GRo7ehr8NRgMBgM6cNYMQaDwTBhGGE3GAyGCSM1wi4izxeRr4nIAyLyy0M+9iUi8q8ickhEviIir7Eff5OIPCoi99g/LxxSex6yK2beIyJ32o9tE5HbROTr9u+tQ2jHVa73fo+ILInIa4d1XkTknSJySkTucz3meR7E4o/s6+deERnqdM5RXr+uNnhex6NCRLIi8kURuXWEbdgiIu8XkcP2eXnGCNrwc/bncZ+IvFdEBr+WpaqO/AdrYOpB4HKgAHwJePwQj78HeIr99xxwP/B44E3AL4zgfDwEXNTx2JuxKmkC/DLwv0fwGZ0ALh3WeQG+A3gKcF/YeQBeCPwzVu7504HPDfncjOz6dbXD8zoedjtc7fmfwHuAW0fYhv8HvMr+uwBsGfLxLwa+CUzZ//8d8MpBHzctEXtr6raqVgFn6vZQUNXjqnq3/fcycAjrA0kTL8K6SLF/f++Qj3898KCq+s0YThxV/TRwtuNhv/PwIuCv1eKzwBYR2TOclo72+nVI03UsIvuwSoy8YxTHt9swjxUc/AWAqlZV9fwImpIDpkQkB0wzhLkQaRH2i4Ejrv+PMroL8gDwZOBz9kOvtm/t3zkM+8NGgU+IyF32dHWAXap6HKwvMLBzSG1xeAnwXtf/ozgv4H8eRnkNpeb6dfC4jofNW4HXAc0RHR+sO6hF4C9tS+gdIjIzzAao6qPAW4BHgOPABVX9xKCPmxZhjzR1e+CNEJkFPgC8VlWXgD8DrgCuxfpQ/s+QmvJtqvoU4AXAT4vIdwzpuJ7Yk3NuAv7efmhU5yWIUV5Dqbh+HTyu42Ef/0bglKreNexjd5DDsvL+TFWfDKxi2XdDww56XgRcBuwFZkTkRwZ93LQI+8inbotIHuvL8G5V/SCAqp5U1YaqNoE/x7rlHjiqesz+fQprZarrgJOOtWD/PjWMtti8ALhbVU/a7RrJebHxOw+jvIZGfv06eF3HI+DbgJtE5CEsW+o5IvKuEbTjKHBUVZ27lvdjCf0wuQH4pqouqmoNq6bWMwd90LQI+0inbouIYPlwh1T1D1yPuz3a72MIVftEZEZE5py/scoi34d1Pl5hb/YKrMVNhsVLcdkwozgvLvzOwz8AP2pnxzwd65b3+JDalIrSA37X8bBR1der6j5VPYB1Lj6pqgOPUj3acQI4IiJX2Q9dz/BLLj8CPF1Epu3P53qssY/BMswR4pDR4xdijeI/CLxhyMf+dqxb53uBe+yfFwJ/A3zZfvwfgD1DaMvlWFkVXwK+4pwLYDvwL8DX7d/bhnRupoEzwILrsaGcF6zO5DhQw4q+/pvfecCyQ/7Evn6+DBzcLNevqw2e1/Eo2uJq07MZbVbMtcCd9jn5MLB1BG34deAwVgD0N0Bx0Mc0JQUMBoNhwkiLFWMwGAyGhDDCbjAYDBOGEXaDwWCYMIywGwwGw4RhhN1gMBgmDCPsBoPBMGEYYTcYDIYJ4/8DaqynBLky6L4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n = 10\n",
    "ex = np.ones(n);\n",
    "lp1 = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr'); \n",
    "e = sp.sparse.eye(n)\n",
    "A = sp.sparse.kron(lp1, e) + sp.sparse.kron(e, lp1)\n",
    "A = csr_matrix(A)\n",
    "rhs = np.ones(n * n)\n",
    "sol = sp.sparse.linalg.spsolve(A, rhs)\n",
    "_, (ax1, ax2) = plt.subplots(1, 2)\n",
    "ax1.plot(sol)\n",
    "ax1.set_title('Not reshaped solution')\n",
    "ax2.contourf(sol.reshape((n, n), order='f'))\n",
    "ax2.set_title('Reshaped solution')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Take home message\n",
    "\n",
    "- About parallel matrix-by-vector product and different blocking.\n",
    "- CSR format for storage\n",
    "- Cache and parallel issues in sparse matrix processing\n",
    "- Reordering and blocking as a way to solve these issues"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Questions?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"./styles/custom.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  },
  "nav_menu": {},
  "toc": {
   "navigate_menu": true,
   "number_sections": false,
   "sideBar": true,
   "threshold": 6,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}