lectures/lecture3/lecture-3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "code_folding": [
     1
    ],
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 3: Matvecs and matmuls, memory hierarchy, Strassen algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Recap of the previous lectures\n",
    "\n",
    "- Floating point arithmetics and related issues\n",
    "- Stable algorithms: backward and forward stability\n",
    "- Most important matrix norms: spectral and Frobenius\n",
    "- Unitary matrices preserve these norms\n",
    "- There are two \"basic\" classes of unitary matrices: Householder and Givens matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Examples of peak performance\n",
    "\n",
    "**Flops** –– floating point operations per second.\n",
    "\n",
    "Giga = $2^{30} \\approx 10^9$,  \n",
    "Tera = $2^{40} \\approx 10^{12}$,      \n",
    "Peta = $2^{50} \\approx 10^{15}$,       \n",
    "Exa = $2^{60} \\approx 10^{18}$ \n",
    "\n",
    "What is the **peak perfomance** of:\n",
    "\n",
    "1. Modern CPU\n",
    "2. Modern GPU\n",
    "3. Largest supercomputer of the world? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Clock frequency of CPU vs. performance in flops\n",
    "\n",
    "FLOPS = sockets * (cores per socket) * (number of clock cycles per second) * (number of floating point operations per cycle).\n",
    "\n",
    "- Typically sockets = 1\n",
    "- Number of cores is typically 2 or 4\n",
    "- Number of ticks per second is familiar clock frequency\n",
    "- Number of floating point operations per tick depends on the particular CPU"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "\n",
    "1. Modern CPU (Intel Core i7)   –– 400 Gflops\n",
    "2. Modern GPU (Nvidia Quadro RTX 8000) –– 16.3 Tflops single precision\n",
    "3. [Largest supercomputer in the world](https://www.top500.org/lists/2019/06/) –– 513.85 Pflops –– peak performanse \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Matrix-by-vector multiplication (matvec)\n",
    "\n",
    "Multiplication of an $n\\times n$ matrix $A$ by a vector $x$ of size $n\\times 1$ ($y=Ax$):\n",
    "\n",
    "$$\n",
    "y_{i} = \\sum_{i=1}^n a_{ij} x_j\n",
    "$$\n",
    "\n",
    "requires $n^2$ mutliplications and $n(n-1)$ additions. Thus, the overall complexity is $2n^2 - n =$ <font color='red'> $\\mathcal{O}(n^2)$ </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## How bad is $\\mathcal{O}(n^2)$?\n",
    "\n",
    "- Let $A$ be the matrix of pairwise gravitational interaction between planets in a galaxy.\n",
    "\n",
    "- The number of planets in an average galaxy is $10^{11}$, so the size of this matrix is $10^{11} \\times 10^{11}$.\n",
    "\n",
    "- To model evolution in time we have to multiply this matrix by vector at each time step.\n",
    "\n",
    "- Top supercomputers do around $10^{16}$ floating point operations per second (flops), so the time required to multiply the matrix $A$ by a vector is approximately\n",
    "\n",
    "\\begin{align*}\n",
    "\\frac{(10^{11})^2 \\text{ operations}}{10^{16} \\text{ flops}} = 10^6 \\text{ sec} \\approx 11.5 \\text{ days} \n",
    "\\end{align*}\n",
    "\n",
    "for one time step. If we could multiply it with $\\mathcal{O}(n)$ complexity, we would get\n",
    "\n",
    "\\begin{align*}\n",
    "\\frac{10^{11} \\text{ operations}}{10^{16} \\text{ flops}} = 10^{-5} \\text{ sec}.\n",
    "\\end{align*}\n",
    "\n",
    "Here is the YouTube video that illustrates collision of two galaxisies which was modelled by $\\mathcal{O}(n \\log n)$ algorithm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/jpeg": "/9j/4AAQSkZJRgABAQAAAQABAAD/2wCEABALDA4MChAODQ4SERATGCgaGBYWGDEjJR0oOjM9PDkzODdASFxOQERXRTc4UG1RV19iZ2hnPk1xeXBkeFxlZ2MBERISGBUYLxoaL2NCOEJjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY//AABEIAWgB4AMBIgACEQEDEQH/xAAbAAADAQEBAQEAAAAAAAAAAAAAAgMEAQUGB//EADoQAAICAQMCBAQFBAIBBAMBAQECAxEABBIhMUETIlFhMnGBkQUUI6HBQrHR8OHxUgYkM2JDcoKSFf/EABgBAQEBAQEAAAAAAAAAAAAAAAABAgME/8QAHxEBAQEAAwADAQEBAAAAAAAAAAERAiExAxJBUXET/9oADAMBAAIRAxEAPwD8/wAMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDGeMpVlTYvhgc74T+F4teS9t++AmGUeF0NOKNXzibTV9hgcwxvDfrtNXV5qj/DZZViaN4iJBfxVt5IAN/I9MDHhlW07oDvpWFeU3ZvOTwS6eVopkKSL1U9RlwTwxxGSm/ctXVXz9s4qFmoEdCbJxgXDHSGSRXZFLBBbV2HrilSK5Bv0OMHMMfwzRsgMDW3vll0Urabx7ULZCgnlq618rH3xgzYZdtJMqsWSivLL3HTqO3UYkMEk2/wAMAhFLtZqhjBPDOgXlPAI8Tc6KU7X8XyrrjBLDG2HdtPB732xvCGwt4iWBe3m+vTp9cgnhjKhbp3NDjrlPysplESqTIb8tUQfTnvgRwyskDxcOVVv/ABvn/ecRo2UgEdRY5wFwypgYLuFMtgFh0BI6YPp5EWNmAAkFrZ6j1/30wJYZRIJHfYilm54HOBhNMVIdVq2HTn54E8MpFC8pbZXlUsbNcDJ0cAwy6aSR4vEFUTSDm3NjgfcZIRki+3rgLhlIIvGcr4iJxduaBzr6d028obJFK4PT5f6cCWGPDGZpVjUqCxoFjQ++VGjl8JJTtWN7piw5rr9fbAz4YwWx8QHzyjaaRRdqRuK2G4se/fr2wI4ZUwENtLKCQCLsdcHgeOXw2rdxQ9bwJYZol0ckLyI5TdFW8XRH+c7JoZo3o7a8vmB48wsH16e2BmwxihDFTVg11y8OiaWZY2lij3VTMTR+wOMGbDKvp5ELgj4DRN/P/GDwMkKyEimPTLglhjmMiIPYosRXfj/vFC2fQep7ZMHMM6RRIu/cYAEmhyfTA5hjhD4mxyE5olr8vzxavA5hnSKzmAYYYYBhhhgGOuzxQTZS/kaxMMCnlYO+4A3wpHUc/wC/XO7JIwhkV1RxuW+AwurH2OJ5dgonfZsVxXb+c67BttbuBzZwHdjJM7IrEOxoMdx+/c43gtJGnhoN22+DZbk9v96ZKN5EYeG7Kea2mu3P7ZXSahoGKl5UjJtvDNNYBA5+pwL6VvALIys+9WG3bZHHNWDXzHpnYH2ahI4pJYy52Wz7doPS/TqcjrNQZZSUjWJeAFQUBQ/nrnA4aMvsLSDlmZuD9Mv+CvibSFAE4iJYq4sN69O31wiDTxmKVQHch/FdSzbQCOPb/HtmZY91+bou7gXjRvJv3h2QDylhdAG+OPXnKOboyTfQLQpaN4hQbNwcXdbeb+fplxHp2iZvEYPRbbssdeBf84sCx0JHkUOHFIykgj1J9Mg7F4UUDGS2aQEDa1FPmO4OcCeFEkrIwZjce5bVhyD+9YgagSWHJoqOtZoWBnkuFWnjjptotqB5rj638jgcKJBtLSMJ9u+weFPBUggmzXyrKwNCFSSSRwS1Oi8llsEn2+3bJnTKBppJnCJMeWvdXPJI+ubHhk1P4uy6ZBLM6u9FjRXbxXfp059MoYkPUcCzMCoLK+7c5XmrHQbSKHPQe2ZFgkhK/lpmLTcCKKywN0A3TteNpVkCuYZ4ImSMkjqWvnj37ds7L+IMhbwgiyWGEsYKFWu+OeR3+vbjAy7zpZGSOXcNysGQkUR0Pz5OFSNE0bybE5kUPxuPTj34/bPT/JrrIo4NMIJXADTOOZA180b5HPy49cxnSjTzb1kJVW8j/DuoWSPb045wMewgbldTZrg89MtrJBqtU8m0RFhupu/Hy74k6op2xkOo53gHpiqyM36pfpVjntxkDRS7GQ7Sdt2N3U52Z5G2AOGCg7dgPlFk/wDOPCyRSlCGfcpDbD04PSjz2xYJWjVkoANe4H+v2/br7nGBFkXwdjhiRZHI6mvbOy6lptnim9oC+5A6ZUeDKfFUOZN1uJGDbhXXt75NWgLRo6uFVvMy9WW/T1wEO3xSIixW/KTwfa8pp2ndjHE5LSARbfUE9Pvkze5vCsqtkECjXqcVt6WOV3AWOljrgViZ1n3LTODdDoft2xtKNMb8cupU2Nq3u9uoxE2AyMkhWk8oYXZPBH2J59sRStci2B6VwRgNLtDnbtAbkBOg9uc0ePqIkgYufDiBWMgDy3dj2PXrzmVYnetqk810zZDu0iuZlVlLbexKnuQD3wDTSbJG1ErI5QqQXfcwII6C+friPqpFeQyR/HH4Y3D4QCPuaFfXMrAEtsHlvi+uCi1ohiT8PPGBojk4ZdPGfhAYkAk89enA56fLLeCFbwtR5Y3+Cbb5Qfb255yKKs0jyLUZQbtpPArtyecQpIpUagSCM8/f0xBREQSspZxp7HiBRTDn0J9ffIyoEkCK4Y0LNigfmDWaDFphESDJv27q6bea5vr26euPOy6jVuNM4ZS1iVowpPua6ZbBEb0MpLK0pPIBDAjrdjj/ALxRztqJVC2CTZDGv75pZpH1SwinKsGD+Epdxx26N68nnJeJHCrpHI62PiF04sg8fI/tjAk88pBDjgn/AMeQQK61nYkSScP4YeMbSyIxv0rnv8s5vYSKot3DbdrAMK4qh68ZRWkTVyK0Txvu3LCOAG7cH2OA8moJbxJZFDORtKKCwXzA2etj0PXLSq1xB0/TcWSzBS/a6B5o39PbMcEJ1DMo+MjeKWxXez2+2aHmfWyxiWQSmQWwQW116n0r1rEHJTG8Eokapl2LHHz5QAbst0+V9c0IEXRmGV1QWGahtYk+g2jtXeucjp01OqcMmn8XcLO6+R3PFenzymo1SJIsbSu9gpKSeavoOOOAOt5rP0ZtQjQzpBKkcKWpNeYDiiffvx63kp1iSRwkwkEZpeDTetemUVotRKmmRUhRjXiuT6mj7dh9Mz6qIwybHP6g+Na+E3++ZtHVAi1AMUgZVNq5WgSOemLM6u1qu0EcgevrksdkaNtrqVPWjkAUKtXBPsbzqkR71ZLfoDdbTeKSLsCvbHHhmJiWKyA8Cvi/xgI4YOd97r5vreL0PGM5JYljZPJN3ecIr3wLNvkhWSQ2qeRel+v8nIYw2lT13Xi4BhhhkBhhhgGUVAjoZlbw2F8cEjpY+oyeaHhhWFH8fczrYUKbU+9/XpeAp0xKb1eMryRbgGgfT+MVWRNrL5yQQwZeB1HGdYMYk3CqHl8tWL/fm8doJZHd1hpQwB2DyqT0GBZIoCBOzEqSpbZQ2c+Ybbsjkc2BjamQTytHHBFHyTuJqx26nj784aiRpg0TBYyjEWTW42OD2HS/pmdovD3iUN1Kq68gkdR6dxmvBpklR9IEvxPLSswoqbs9/c/PDT6dpnWBWYPKKViaulPlv1uuMhIFQBVkYwkml3Cwfl9vnWVZdTq9NvZ3kMPI3Pe1as/x3yhBLpxLKQgVGBAVWPHHr/nOJOdN4kcbB4yVJBWw9G6Pt1yYTeqGNAKUlrbrXN/8ZZ4vDMUgA8VqPhslryODd83kGeSdpXLHau7ghFCg/QZZhpnCIjmNUss7ry3IHQX88bTLDLPJ+ZtCb8qhVF8nv05A++VliiOuEc0pi0rLuBHJ20Stj1Pv0vJghE0KwyGaK3NeG23jjt1H85aARcflJ3jk/qDHYWHl4DDjrf2x9BqWDeHxuYhWRlCxlAK5Pb7ZKQGKF4ZEeMhyUXubrqe48vbKF0//ALhh+YaRlBCgluAfTvXA9D0xkeXTwbHhXa+4KSpVuR/5dxz0vHmGyRjGUSN6HitHtN96q6969cg+oduNiUFCBjZ6EeYX34wOyRxxtGEceFLTbpALXqD0JPW/8Z1plTSJCQspDFgeeO1UePfIm4H5KuzAggqCB8v847TSSLKQ8ahqYr688Affpga9VqoilQ+FHMPIfCUkEbaJ3E/tVDtlHj000h8ZpAkZEYLSXSAEkggc9qv1zFpncbnj5kfyGMIfOD1HGb9CGbReHLFLW/bvACKhv+pqvj60PTAzTxvFAu9I301gBlPLdb2k88kHEOlhjmRZN6xszKZG+H6Fbvm+mX1MUss0unSZ3mknNxKTsbn1ocjk/UZhm08kS3JSEGhGzeYfTtgTjleIkoxBIo/LHM7yqEfzcALfVRd0PTrkcomzqy2evxUMgvIi/l1ji3Fm87UwIoAmuO/XMpI2qAtEdT65og1E0MRCbKUN1PTcNpr6ZGKRlDopAEg2tY7WD/AwH0/5glk0+8kqQwTuve/bAMojVXo2QTQFgC+M6kQZSFYbxZIPSgL4/fF2eGisw+LkHrgHh7FHiI9sLHauvtjvCYwBuTYxPnu7oA/z++NNN4jKyFmCjaAy9h07n5/XIoQfKVLE9K63lFlSNIy+2ZiCvmApa7g/XpnJNjr41lCTyu39x2yqJsZY01KsjC2XaSC3oR3+eKUSLUESx7lQ7GVWNfO/3wJRNEJCWU7OhuievbO+Nbqnww2DsYkKaHWr6n+c4qNewDhyKBarrOmEq7oFZnTcChvgAdbzIVjIHEvl3PZ7V3HTtjorSbuDJL0IIsAV6+uJRjiSRXUndyhHSuhIPXvnZW/SjVEQE2Symyb7H06dPn64EmD3b3bcgnvloJVSN1ZmjY0AyDnvf05zjFfFUs5r/wAgLFD0GVTdFNDqUSSm6H1aua++WB9MwYbYWlEhBXZGCWquT8uvHHXFMTQ6gb3KvuIV0NcjvzVc16ZCFDvU+U9zfO0XXOaRIRKDCrbyfJAoLKt9+evy980M6RuZjsewT8SWf+cpHLu1RmCGZY23VILsdBft07+mPpYG1MjpGfAoAbAT5mAPr7j9+Bllih3NFvfTjnqhJcigAe3W8gxRwbxakcck7v4/jNXgs8bgS+WLaQIlHJb1IyxVIEiTw/CkiqRiwPnB47j/AI653Th59O8MSQjewC21MDxRrv8AT1y4Fg3pLLHp5Zk0yIzXJXWiPkOb++S/FGK6lURo3VACrCjfHzPpm2Wb8rE8eoSGYtzbSbnVrsqfXkc55xU02qeWMP8AEu0rybA+Ht69MtQu5VjkeNYzuQAooYhPfnv/AJzHlXkeR2a+ZOCq8fSsDFwQgdnF7l29AO/98yqkgjZkUuijbdovTjv9ciKBDHkXyL5ziMFYEqGA7Hoc4Mg5nbJFdscqmxmDEEMAqkdRz3+33zkSF32ir6+YgD98Bo1j+KQMVrkKQDfNYoKkneT7HC2Wmsi8AR6D64HFNG6B9ji4wLRuCCVZTwRwQcXAMMMMgMMMMAyhG47VJI7Xxxk8vsMjEqRZF0BQ+XGBrR3nB008ccUUS25VQG8oPf1Nj9vTKQPLBovEiURl5+NU0p5//nvXW674ks8+mijikkZo2YTOu4gncP24v53nIZkkjXSAfpl2KGb4QD/USOQRXuMKyUSxBVWZjYfn/ec1aieTUTLFLEI5RJbMOt9O59K+2B0rpEmyVYy7BCrBhusfFZ4rmv4ycmhfThDO6rvXcAOaF1z9j0yooYteXEDQmNmYku3lB3UOSeK4zCQVHPQ++VkMgPgCS1/+reU/7x9s66eVEOxLtixB4voP2/fAbVERiOIReGyLyeLa+eo69T9MSQxNRRJEXn4m3X6dh2rJk7nBboTz6ZSkkcJGQq/+Tn25wNEqPGzaadkgF7hYDV7Eiz9PXOalBCYWOm2oY7Tn47sg36i/2wmK6XUkLJHIyDzOtukxu+b+nbtiwyGOdGhDr4YoOTtonofbKE02okgnd0epCCAx9T3wjLz6hfEd93JBjUFr6j/e2Ui088oUxbAjtvW+Bu7Lfr7Zb8SaZ9d+iOfD8vh3u2HkA8Dsa4FVgZIJpEnjcyBirAjfZAN9cNRE8MjQvQUMSrFaDdeR7Zf8jqZ5lkY0k5J8UjjuTddOmWP4e67IpxtduI9xsN9ug5/fvgYET8xIx3Rx9+eB6cYyeCCVYUwVuWs2e3pWepqPwPVwxxhNOjLNTKQ67rF8Lz3+5rMsuhl8VZNQm1WQHmxQqh15NWOPbBiEPioUBkqNG3EiyFvi+Oe3bH0pRPOmo2zpbKuwurexFe3PUYphQpN4ayPz5GQ2oUWTfrxjfllRmaN1kjQByLAZh3+XfA9DwdRAT48kchZSsqlyAqcAdB0pR+2eXqaWUhpnkKAKlgEe/c8enrj6aWQagSyTuEcbWkveQPcX7dDhrNOsGpG6QvE6b45CPjHIHHbkVgTSSKGd2hAKFKCy+bkij0rpZ+2cnYPqGIREVgDSAGgef9+2CzRoIyrSblUnigQ3aj6dM7HAA/6kjQ0t/DzdEgAXz06++BKQqzhVJ2qKB9ffKQaVpmpeb8qnoC3pZyIjIm8NiEN0S3QY4m8mxrI3E9e5yAbeZC7MFffVDykH5dsYlZA43KFHIJHP/eUSNXDsrqyP1Z7tOepruf5wLaeSR2fxQp5HIY3X8nKJnTtHzKrBD0av99RxlmKPOh2QqqggKo5oeoJ5P1xJ9Q8rGV5y8tDzEf7yMcPK0JCadH5L+IASwoc/Ic/t7YHJTNHN4dyEA7qJ6j1vvYzjgSRyshKgNaoVsgX3OTd2j4VgDtF1zff9rrOnasiF0LXy4JoHk9CO1VkFkgeTTyTSh2RR5XAuq7fLpmVt3Elnnocs0wk04QiJFj/8RTPf/XtkldFYuvDBgVUqCv1v6cc4HZVNhpAxa/O274r5w8QRqqrRptxFD6f3ONO8jyHxW3bRx/TQ9AMaQxtqFewdw84ck81ye2BJHTvGBxQI/vjvxEH8VzLv8o2kcet5xQikkeZd3AfgGu5rCSeQxCEylohTKgJ2qfkfmcBtM5hckiMFLNOL3HpR++PpJp1pdNGfE3AbkB3MTfF/x3rMrtubdQ568ZoR9qhFDRiQg+U2QOn+co2NMx0Y5RZt53sliVgT0Pbr/HGdYLHCAd/iqz7lccqLo2fofXr65mWYLEheIWG3bnJJbnnnt09O5zVoTrA0gVWYyyRl7NndyfUEG7/fKA7VmeXSJLEpVuLU7+RXHoTXHOU0X4fqG01ssyyNIJFAj6V8+nQ+3GIhLMjod+0l2KEkBjZoiq7/ALdRlIpZRpXOplkJVCIjVlfN/V8rzSMhZyTFEVBQ3tI3eIS4quvtmWRDLqXaSRAS3mZel3/pzdqDJHqIlg16yPGtmz5VJ5IF8H/OeWCRweRzQPTM1XVR55wkSFmc0qqOpxkk2srhQSOt8g/TCKEOpJlC11BBxzGFFbgOBwRyflmQkojCokYBYXucEnd6cVxWIFI4PHzy0hkikYFDFIOCDYII69e+dZ1kjrbfmsyN8RPfn0wIutIpLWxPT0HGEkfhybGK2O4Nj9sqwPhoyRbUawGI+Kvc/wAYkiMI45CqhWsCj6e3bAuNG0kqoZFIK2vPUE0P75GJljMobqUKi1B5+vT5jOr4sjvIHVWLUeQvW/p2yQUs1AEk9hlFDFtlQSsKkAO4EEgHv16+xyGUKHjaGvvY9OuTyAwwwyAwwwwDNMDLBIxfxQydDG22vrWZs0rOsTo6+c15g3r88C0Krq0jhYbK3HxAhZmoXX7YkUpMqwugalMf6hvaD6X0oknFMqagSSS2JWYU26h3skUScPEnsssa3MaBWOrI/wDHj+2FNPN4esDwuzSCwxoAXZFADtVZteb89qG0+rXwtSzKgd1rYb83AHrz07nnPHvLQLvkCtvuSgpC3Zv98qPQ17UkKeK7MqK1qKAU7rodu3Hv63mOaTxoo13W1szFjz/x9MvpdUNPqG8QJMNQAGNE0CeeOOc16j8nqHcJKigR0I2JG07qFk8Ei/tlGNSfy2kVZBtLlWLk7VN8jp0qiavIEsUaCOJGcOWMkYJJHSvllm0Ekf4eZXkB3P5UUbjQBs30rkdD/bBJmZJWG8xouxHdQaFihfyB6YHZBLPpdNHuZmcFhuvmrAA+3b1wd4Wh20hcCtoJXqOO3bnqe+aoiuoAWGUf/HYXhSpF8ihxyAfcZCVAkhXqyRuzsihRdVwe68DnvZwH0f4dP5mePxIlNGND8fUWD6e4zuniTTTf+5qCt22Q0/mHavTt8/lm3TOdHGiDzxlDIGFqJDxdE10F9PliQaabUQsIwjsZ9zGlJdqsKCeg/nLiCbVRJD4bacRyBQkcqgg7wR5gKHp9PrWSYLK0us1okkjWTw6qtzcny0aHRQR75WT8O1H4hrIzLpniidypES/ATbcA1xyLP/GUf8N1EJg02ybzceIVujt6GjY46WLrDUINa5laLTxxx6OySJAVr18tsbF0K55+1YQur0zazwakQMy2eHa+vFURwOncYkcUmhcieAtzYaQUfKRW4dgeB9cJ4nkng2y0ULONvNkE810PSv7jAy6oSpp6mjV96lo5RYQA1bHi75A7AHtishi0zHYQUAYusm1gTRB47cj7Z7aJFq/CEe2GNhsDzD4bq/KbBBNX6Z5cpYMrFk8Eg7mhWmU9Txx8698FjG0XhS/rOs3WjYIY9PW/qe4x4odRLpkVtOsioS6ySNVL0rk1Vj/bzTEYZkhjlltoUcjfZ6Ub8t8cEcHteTmZ5LgSRyjvSutKpJFkEDtzxz64RjfTI2qJEkARiwNtW2vbt7f8ZlO+IBlIAfkGwWFH7jNbGKGRohvdHQbwP/K+3qP9rJtqU/LukTSqegFitp7evr88gikZljkfcirEovoCea49TnNiBmtqUA1wLPpxiott8QX39Mo0O2GM7lZpOnPA5qjkAkZcFUYUTZFcgDv8ssgZI42B3qSQUU1fuPpiMrmNZQ5VGBSyaHHb19Pvh+XfwTIDdKCPqR/zlCLMNojfd4QBO0f+VGjnfGbdGA7R7EK3uJ9ftd4OIzZtQbACAHyg+/t7+udk1Mm+y7bxVMDR45B479ecCe8hHj3KwIBurN+l+2dZBe7baVfkuu1iz88UFQFDg0Oa29fr6ZRvN5oNyeIxXYpNVxx/vpkCru8U+ArG/wCjqau6Pr0yrxqum3vsU7j5G4a/b7Y4LLp13bxIo/TfftAs306+v3GC6dtRE7MzO8QogUdo6X6kfTKJRyLtdQEsqFXco9et4kcMrh2RNwQWx7AZZoPC0mn1NAhy60O5Ff5HpktztHNK3iWxCkg8G+aP2wGmZZQzq9VRCMD1oAhe1fPsBkTIpdm8JQpuls0vy57e+dkjaOhuBB5FG+ozu0UeQzMOAF6f85A40+9f06JC7m56Dt9cpAZogkyIyO4bbKf6q60cy7W6gGrq82BZEcToyyqLCtxdKKuuw+fplg6sm2LYVZwOpFFR/vm79/bK6WaGZG/REb0sYEdixZJJPyrreRlENsI2UK4Xyg7juyw1XhCN4iC8S8goANzdQBx/ozWjQEim1C6dFCbq5BtOa6n147dcprkij0b+Lp2SQgqhV7pbFX6cjuOec7A408ZklWORJglMy0EYgqODzxXX07ZOaWVXkmldGjW/CJjoMyqBe3/ec1+JrFJI3gxmZTuC+GgZAPJXW+//ADkJ9mzw4fPGnmZ66np6AgdOMtrZZt407eKFUKUVuTfXqOvU0ffM43ooJVjG4v8A/YX/AJzFUzmMQmFZCTv3WD5Kr0q7yVM/CW9Cya6DL+GdwjZ4kJXcAQKHHc/L55w6UyuPyh3gi9u8FuBZNffIEZZTCpdT4dmjXc++PTwo5EZ2P8N01c3V11+2SDSTsiB9xHILUPnz9MsFY7mIkscuV7Hn/nECT6mSaNVY0iE7F7KCb4GIYJBFvKHbV7h0rNEMcE6CPwpVkCMxZTu3EcjjsM5pJlVWhlVpIm/pD7a733547jAzQlRIvifBfmF1ebGoufDRNNtoGyQ/e+f7/wBshpQBqlpo0DWtyDcosd8keCACbHpgXklLu+9gdv8AST8XarA5P+MyZom1UkuommFRma9yx+UUe1emZ8lBhhhkBhhhgGa9Hql05dim48FTS8EfMHjrxmTNSQSTxyvCRsDC4lPPJ4od+TXrgRldHkZkj2Aknbd0Owyyu7xNvehxtLD07A9u2dVIX0oBZElD10O5gf2oZ1A0U8+n8YralP0wG30enB6cdRgZOhysDuhYooJrqR8PPX2xVjZydq2QLOajOjIv6MKAvbAWLHHHy4PT1wJSwNtsmRptxDDbY477r5xW8WIxS7NhABU7eo9ctGniOZGJjTzDgkdugJxdTzHD4krmXbRVx8AHAr6ZRr0WufQwyoyyRyOGAkKhh2NUR+/uMvLF4iJp42X4RJYjIVSFAtue4I7elVZzzJQJYkMe47LQkgC+TX+2ctC0834iChXTOSWDXtCD2+2RWmDSNNqLiIcrIrBlj2gK3NkDt2rtnrxfg0CyNqtZP5ev6pFOBzyL57fXNEUk2k00rI8QUVJHNLQFEFuBzzfPHHHbNGn1AOjiIjmYuxklIJ3NY4IHACgiufTvedOOJTw/h+ngiiEUkcYc+VWG4Cxe4A8g9Bd40X4dqZGkhDtp4UY1EqdV3Xe7dwTfX2zVotJEdNDLqiFLhVVI34kAJIJHTkD9s0aX8VkkfaNOPCPwMh6DoL9TwenGZ5WTx148f6XTaBoCP0fDRyA90zsK4B+RODfh8kU51CMxdGIVQoBVWIujQ9Ce/U9ck/4vJDImp1DOsArcUClTd13voQPnmjSayfVKZOQP/MqyLt4N2a9/rk1chPyz6mOVGhUdg8h5J5Bvnn0zy5Pw0zbdVG0YV/jKOBtSgauuaIPIz1pPxIRMQpUxrVgclifTjnkjpfbLaZg7/oQJCdzC6A83uOvocdwyV85Pp59N+ISSoiSxyOqh432tfNqPNz0vk+mebrYoJdLGs0bvOWbcUbdz3aulXXTpVXn1s34cZtSr3B4Km1Vhuo1W726A5i1H4bHpPwsQVHDuNKjD4jfci/Qf8Y1PrXxczRQTzKpUS+KVDR8LtI6/zWaCp1bAGNH1Pxu6tuD2OLvvyf2zT+J/h/5PS6eVyyzTRVIykkWbIFV5en2v55lSAw+DLqZYkFKtIQG2+p79OcsYxklJih8GZls7SpChuPZuoqhxz36ZC2Yu/CNt5ApbHTgfxno6145iW00ZIdjtRF4G0AfDXuTfucwrHNqfOIwfDTdyDyo4xYy5FpJDGZFljVlYAAvR+YvtffKu/jgRzEssTMXZWvcT6HoLrJuwdy08nnJ3AR1QujwBx68cY04gcctJHIqfAVtbvoDdgVkEJFKKo2lb5o9a6jGnl6KFUAAdjY+//WCsvl8RWZh8RujXSucJJX1Uu6VhYAs9OB7YHRK8rEWoscqOFavbJuaNo5N3ZPHH/WciDF/0zTdBXU3gN0S2ARZq8CsjnVy7kiSMgBVjQHnoOPfvnq/h8Jm/BZJFWptMS6td9xxX1P2HvnlBVMYSK3ZqO3mweb4+3OfU/wDpeQz/AIVrIVkAahQK8dOlevHpz75rjOx8wHSaV/EUuGoLtG2j7AcfTLQOI9cwjj8QSGkscNZ45NViyK0SiJgRt3eGU+K778df97ZSFtRKrRxJumdgRGqilPABI7da7dcwORh0TcTIunD+baA6pu68evA+2RZoiW2mLoWFJ170b6V0/wA5s0kUc8EnjFvERQ6BCFJN/v19emJHEqxKsJcySqUAKKoJ4/qPXr25rjCsocRTEFfGtfNTdD6g/TIOjBUpywbkKLsf7/GejP8AhEuxpkR2jX4mUhj1oUOD9cxySTvpo4pfEEasdgri+Lr36YHEbcT4fkVvi5NAd/U1mnTNs3SOVYNHRet4jHQWvr0r+clpHAmRNqbGbzjmyp/p/wCuc6UUIy7SiMaV3sD6/vhHY22yRzUu4g7l581dz7981QoAISVicOCzCSNlc/F6WSO9j0zu5vNDGjOH53igpUCuhAs/z6YRzHwhEqiOVXYOGPnYHjab79u/0zcHZZAySkyxh/MBG7cV0FfQng9MlBqQyFAxBRSsayHdS0LPNDtwPXOL4jaje9GNQEZd5AC/U3V+/f3xFceJud5A0l2IyAdh9fQ/xl0TfTxqrpIZBMhavLxQ9fe/plNhpAEVPFamdgNnPT/9e947RtPpKSZTMoAKbgpcduL5P0J5+WZzBIm2B94UPtbal0e47WR6ZBOF3MxNK+4U28XY+fbp1GMzoI2aOEqC1B+ePl7+vzyKKxRyoPl5J9B0yg1JMJh2qIz1Fd/X5/LJopHE0OsieGNpAGBWxwxFEjj3yy6khJ7ESeYWteYqeqg5jWaTxC0YALDaQFscisafxN8aygoCq9SSCK4P2xotLO6G0l86t/8AIhN8dOcyrs2neWDDpXfGQu42ILq2r5Dr+2dmVXLSRbqJ5FfD6WRxgTlZWclFKrxQJvOiF2haUfCvX2xZE2uwBsA9arOxhnbaGA4/qNZBy+DYxcckeGAB5r5Pau385PAMMMMAwwwyAysdB/0g5c1sI4Ib6ZLPV08Y2b0/MTvDbKw4VYx39Ryfp9cDH4O1t17lC7hfF81VX64ixGU/prbk8IoJJ+WPEpOojWI73JDFgpNd+nevl65rd1hXekpapW/W8MbyNvHB6A36++FQg8bTs9SiNtvKlviFX/vzyMYSaUiaTwxXBokD2oZ19LONQsLJUzkUnF88j++dKSmTwZi36ak1128XhFJH8SIyEtIwO12Y8dKWvoMoUSdkOoZI2QKGU8eTiq9T1+mSicww0GjJJJA2gm/n1ykezxCZpRFIqkhtpPP/AI8ZrBnkZwAWQrExJUHoe1j1z6L8F0GlMMP5uN3lsScHdSDkVV9eBRzxPC/NSrDG5LKwjUV5T15FfIet59J+MSo0LwwKUXag3c7lF/1ccEGqA9cT+iDFpyulm05ePTgGt3mItrbkcXQ4PY++ex+HypqERZY1jCRhXRqvxRVD1vv0zwPxF5oE0rPM/wCYC02rViwI5Fbh1r/enPqyK/8A/wAwCJBDtlDlkUbyQKJ282au6Pbms1hr3YphIrlJNscMphXa9iyQCCPa+PlmGO9S50ohl/SrexO4csQR63ye4+2X0EIOmERhuNT4oA8rgnnv8N897++erBo1hLSIAm424vr7/PnONuO88R/JOFaoT5QqBVYgV7V6c9rzPrQ7BBAioieUln2gG66d/YHjPfj+AUK9sWWBZEKm+epBxOeJ+vjZt+n8LfpVqmQuFoK3Q9DQBPHPtm5pSqIYz4Pi8EuxsdSOKsdx17982/iOmJ07gqu5j6bgvboOenbnPE1jsuomU6jYdSKSQEqNwAFc/XvnWTYz49SVmBlraw8y3Jwwvmr9Onzypji1CxMrGJ08ii9qyiuhrqOe2ZoHGqgTxCGKkBXHqeN3Xkc/LplhqolnWAali6hRtG0URXSvY3075Lx/jcqGrp3KyxOokjKUW3KQCavy8nn3zwPxX8F1GmR/KrLflkAXavJrgAewJOfSawJD4YLHxdw5C0FFiz6dzz74uqZ9XAFRSjbbQHy0eldKvjpnScXPlXyiaKX8lMu3czU6stqRRPt8jR9M8nSQTM7RHxIgwCsoLA11/joc+o/MadNO8TI8LBuUPm8oPJJ5u7PX16d8hqJ3gKiPTsA4vcpILcHt0zV4SuP2fLsEeVEkgKP/AFbTV8cUD0zszhNMIwzP1FMopOh49+P3z0S0Y8TzIVV2k2E1ts9r56D/AHpndV4O4mRI2D8yKPMwA73ZI7c8XnP6Na8yfezLslDNKPMqmhXBAP8AOLMngNsdU3cfCdwoi+vT/vNr6eJtpWQvaH4FJrg116D/ABg2hRY0iaQt5jtCkckgcXkvFdebSjsenByh+BGK3Tc10Htl5NBPuRRtk4sIJLIBxUhRiBIWjX128HtdZMCRsIZFnCsPNa9Og9D69M93/wBLT+Jq9Sku8K8JPiDkqe9ffMEGnCKxlTeWSk4sMfQGjz/jPU/CQYvxAaVt0niqNwVLDD5npwc1xnYx/jOnkh/EEMbl9xBMidWIHWgfQ365mjjmOrkaCEWJCeAQRZHv7jPf/EPw+SbVyu6GPTNRJJ2npwOL45IOJqQkU0cihZ3UBiGHNqCOvU9fQ/zlvEePp/0oijiUBx4nhyXTpQPYi7KgYaaSZ1DwwJvARV2rZTktam/mMacziUyxuYlJJWXb5CLPPHK8Xx3yDySnYyTp8fLht1sAOT7fTOdiqq0R2SQ8gsSym9zNz79OlV1x5JnnPgaqTx1hO6IONrtzRAPXkf2GSeSZvCWExkQm1k8PaFF2S3tdeo9Mnpg8sibRPK5FcAUeLIJo32rJhrTqtHoFjLaJZvGb/wDFIAQv/wDV9qxF0+oGn8XbYfh0a2LUOSe9Ak9Mizy+B+oDAyMTtZSCwJ2kAdOKN5s08XlbUalSFVzaGCzJZ6g2KvkcZ0zUZJ31EEHhuY4lB3COzbbjzXoOBxhq3g1GqlNBEUBdxO4IK4Arr3vrm9vwzXUhWAzHaW2SkArxzfoBXw+30zNJpJTNBckP6hBfZqEO49eewHzxgzLqPBkmiEREospIoJaq9z/45jEdMqUXfqwU8VV/5z2J9JqdUfFm06NMpEZWJ1sjnrR6npfaumZzpptEhk8B5QSVk3jkEUasdunQ3konLt1EccaCGN+rHheK4A7nv72cWXwBqjpxvVQSD4zcI3Ivgc9uwyTzn8uVdVpzuIA5Boc2bPPH74RQHVufy8TbyQAi2bFcm/8AeuT0RMjLK58puwQB5T9MrCsCjc0TzxhQZCBt2H5+l/6MadYGUxxhkeJqrbdjuSf+MzMkhk897n81sav3s4FFd14WRohXAU97sX/zj75tPJHIJVd9tjnds5Io334/fJqoDKqSqzk8ODQHX1AxfJ4NkjeSeOcCsMy0FnJ2RqdgWgST71/fGVNOIAyzPIzmjGBtKm+LPIPGZvLtHW8CQtgWR7jAs3OoB1RkCn4iBbfvkytxb7HHFVlQZ5oJJnDSKGALsLo8984oMa1MjMrpuSnodav36HjAmAAy0Ax4PrkcqQpjJWwR1565LIDDDDAMMMMgMvBFPIrPCCwTqAeefQd+mQzbpNixyu0u1thChAdxP8cA84EmjaWVkjgZZbYtGoPlA7V14o9cGn2uu1eFYGmN388aTUM0RDSu8sjbpGY3YoVz1vk39MgEsMdyiu3r8sDTDERPMijxGVCbjTdVc32r54LEJ5JJYk8NA3wA7io/nK6LZNLNIyogAB2rakDvt7du/rlXnEUoXRSELGdwZiPMTzf9rBJ6ZqQTaV9OImRRHPx5gCCw4IvmqP8AGWeaZ/w5JBMfFiBB8PqoJIO4gcWOOp4AxE8J4leZXkQ8AqwuPk9ugBvoe4yr+bTv+XeR0UnaSvIHFhu1c/2zQ3fgqTF4J9G67mdElCEIDyeDZPtzWNOsXi6gQ6eWZ5rL0SNo3Xu6H5cZi/DF1P4fqIvF0rV46FieQoDA8Zu/E4tms1CqZXcksvmAZCRZFD26/XLIjLp2m2xRapGlhBPh+EvKA+oA6H/GevFsSMakbo+TGrE0gXtxfuOvPf0Oedpp59ryEBJUhK7QaACjgiunU+3XPR06vqNCmlJjMaEEMVrxAKLUT1HFf5yyD2fwomCk1KxQzSWNgtg2zgke3tnuRmlUEjaejXd54P4dOx1JfVRbWN+dR5aoWPbt69DnrwONxSM2B8Vt0Iq/75i/Hrvx5THoKeOM6zccZGN6Sna6HJ9cWWUCHcCOelnMf8u02esP4oGaI+FTFlKbTwL9b9hefN6rTz+BqY4ndfCO5NvLEURt5v0u/fPc10krJS8kHoF3FeD+3zzJp3/MsoZtzbAUkQkFbteneiDnp+uTGN1z8H0pOjWNxGRIgZiqV057VfXsf+apG6PHIi+Pyd0pbaL56dT6C77e2DLBFCI0ZNNZZ1VIjvBPBoH5kZCLUiPUlEjDb1BTkBQpPJvjsRx2vMT+t+NkymORRI5O+Pb4fUcd9xHoenTMv5nfrGjkULuU7Sp3buP4J6/5Gc/ENRHGkcW1JiVFXRZz2r1sfLvmXRzw6hDJDt0wvabBtfSwev8Ax7HNudJr4pEmO1NjSJu8zA9bI5Pwn4ug988XURNJqdsaOwZAPDVzIavn+BQ7Z7n/AKgjin0UazLHHJuIMldaK9L5PF8+xzwNdHENYo3MdwUhbO74eVvjJfGbEGVneaOMoYzudVAAKHsbr+2R1O3YsqkgOKsAjbzYbrz39uuWnkGnmAjdNzoBtU7ioA6buP8Ae+E08qRRiSAU4UbCvlHcAG/Qg/O8wy5ppFFbSZFRSqxrx17+vc9sTTjwpGkhlKyJHXPxDmjQ/wA0cWNZ4tYhjXw35k3duBZAP0+mSZZp9RbKp+I2OQ3+nIsa96CUaiJSgDKdhfk8ghe3v7cZoTUh5wrNHPpz5V3xn4rJHP8Aozzo3J0pi2m+SjBQPW7Pftj6M3KsMdM1dTxz9ehvjt0GXVezIIizRR6WxssbgAAOx49fbtWex+GbxJ4oCjaGDAgdvSufTr7+ufPRTqizBNyxzCgzOKBHUX365634ZNFqI/GjICxKUVKKbuB378Anqc6S9pfFvxiXWPooodK20sakkskKKHQ9/S88gaNINIsupjLPv2xKW27uhFEe1/8Aee/48ckZ8MsincTYALDkXfPoOflmWWZ5Pw9wkUMkxcMqqQ4APG7pRy8okrxoHjGkfTGijShLPXaehF17dfv2y0P4POJN+kXxPN+pHIoAFHgHnjkHpzz88+iGm0irHJKFldfNbEHoQLHOCAbApUxwyXYjajfNnj19uePnnO8WtfOS/gcksf5piARW5YVtWJ993PPtloo4oA26EvEo3MSpXcCeKUEWRyP+89EyRGAyzKVRADuav0rrkd+4P985q3jLsY2XVRydmYMTyBYAHTr64xXnP4Ess0v6TR7WCN8Hyvg324HzxPzc2p0c0OpayCrbmPxC/iHPFDNmp/CmWWSXSQKsJQNsV7UkEkWDVj+1cZhhgmhLIYCjEMofZYF2aANAHng+/wAhjtCaR/zE8kE0T+G54KsGNE2bJPPDD3+WDaQyawNDDG+1C6IAobcLAFjr0uu9Zh1UE8LkCRkdBt2M92P7Ch73jQ+JG+wKyM6lVJCsbHIHT1BGTRmkkdIy20buQK8gG7vQo+o9MaLVt4cMQkcGJWPlbYa6gX359svqovzu508SKYhQ0RAWxRJ4HXkA3mKJpvGV4wWnvaK6oRVEfvmFaP8A28aSSvpEbyCrYslnryOh9MjprcXE5jJsWFJLX2Ha+wrNuqdJoWVNyqkm5VdiFY0SwrpfHb781mVBA0iymNwVYSSCLoq8cDueo78fvlGaNDIZXLEMDzuBJb55TUJ4hMsTSPEvALAnYOwPb7ZLaJvEnkbaCxJA9z6fX9jnIg8gVFiaTqxUDr9vTIOuTJtVQPJZrdwfkMWQAHlUUubqzac9M458m0oqk9T0P+8/tnY0V5IwbUn4i3A9iMDdPoFOnRoNNqN6IviMzCrIuwKvpZ+QvM8ywyRboUoonn22VHNCyfp7Zwy7FKJqN/hyEoaPm9+fX09hk2dhCxXlCeSVHB9sASKXykozF72i+px98k8YuRjHHXkZ6AF9B9z0yHisTv3Hfd33vGEbh1jWzI54Uc3dVgcdlcny0CxN9aGRzQhtkWP9Nx1bdVn+Mz5AYYYYBhhhkBmiMLJHtS1lUEmzw/8Ajj75ny/hKunEhkG9jQQCzXr9x0wCFkUOSqk9genXGkIaNY4md23EVXUdv7njOuZEjWIJsRgG5q29wf4zkkbQmSCWJhqA9EEcisCulnkj0+pj3kCVALs9u3HtmqVmLlKd1kgUM6Rgc10+XT3zPC8w/D44496CSVlJV+G4HBHXv+5xonhjmV1nkqNgWAJUuK81Ht6e95uBoZv0FEb7WsDygAivTvfvmmASzReEPKzH+rm75PAHTpmRJ0iklmSJBEwqMSJdkeld/X54vj+ItM6Ag2QN3nr19B16euXRul1iDUaWOVjPHAPLb0D35Pp/t57mp16alY9SINrEMztGykqAK3A3dWfY8Z8rponKFqkUqp2uANv/ANrvr5b4z1vw14p9E+lQEBCSkiLRcWCwa7HQX6Zr7VFtedPLN+WLy0w3I8a7y3/kPfkX9MaEaLRTJJIZWeDagBPxXxXPoK6emRSAzqNVGxh8IMeWO0gDp1rk7eh75eBDGFEnhbiLRmC8sQD5jVEX/wA85Z2Pb08rLL+nujElBhsBuxfz49e+erG4HmdthPNet9PrnjxF45bjSSRz8TGrFX09vTNXiKkqq0UjODxwLYjofsevtnaRNv49UvsiPIJB59sRHXa0Zuuasf7eedHqlhCOVZexcgV09fnxjCbfuR1Aojbt8pPN3wenH7HL9VnJo1E0JKK8iKS2277n+3/WIRp4wiLII2iF2LCivX/ByGqDXETW0MWclQBQ5Hp3A5vPNZ2VyWjUhBu3P5mUCwG56cgf85jl0s5SvW1G+TSIqzygr5VkWlLHiqB4o+uefCkCadQXVSo+JuQ4uxz6WT86+WeKPxPUq6TPuKciPdQtOQRY6/8AGWP4tAiv4WmJjDVuFivMDSj/AHrfGc9i60KA2qTUzFfCIott5Zj0AF+2SgmMOt8KPwzG1Ei2IjVeeRxY5/fPPPieHEJXqRvOiizZu6Bvr1/cY0Mcg1MUjqjeMKl2yVwB/V6cf72xpa9n8U1Mg0MckTEqhYUq3QraR7f854wRpdLM8khsgOEeS+a5HW+/U5f8TEx0egicyQxmV90pG2jSnjn5855SkSa0yltzsdyptskewHHbJeXbHp/BE8kzL4cPlBZnO4qeKpunP84r6WcRaeV0EgY+RF4F+hPr19+MaVnkjkMEzxrOiPIXYhSbpuvXnvmWKeUhlM2yYWUugtA3wexJ/wB5zNq4pPGYNSoXayox3y/0Nd8Djjjj55mdJAWJ8nPKHg3xx++PB+JSwO8jIkpkq9w4NV/HH1y35toWikYK5FsoL7gqn+nvk6CpCiswkis+J5Q1rQHY9heV0ChdQojcM+34evfpx7+uVj1KHTSzQFoj25bqSL555HJ+WX0xEerM6FRDKG8MsdhFc3f+fbNSKjTwrtSCNyNw8SjywN8Xweg9eD756UMkZ0aSJp1Qgg/Jh1FduPvmIqdRuEQjaSQsFZF4NGyB7kfxmrR6XUawTL+ZHiKu7byvN9f2GWI9TSsjzyxOikhwWUkbSCOneuvXjtiSxrpnaeFEWa+CbAZaFr7CyDmaUCHWBvEK+KhJDRkmzwT7dP7Z6ukg0f4hoC0m/wAYGyLor7DNpZ3rD+ZTSQQny7pJGUmQEqbHarNdOPvnofpPon04oIi3+mQps83x0+hzI/4fO0Saf/5YUfzh6Q1W0dPW8fU+FNpvEmXxAtxEOaJoUSAPXrkgyavQ/noSVlBkXyLuY+WuwANX7+/fISt4LbJWKaeJdt0L3c88nqOOL4s568UKNIzKNoSMbfDski+KPFjriapI54yk+1UbzSWt7TXIBPSx3xYry9FqZdFDSSsDKQok3KVXnr7A+g9fbPUn08OvcSXtlji2lo2NtdcV8v755y6GODVhInkkjJPDJSkfPrYNVnIvJI0QMqNM6GTctBbry8jg9wfTIrL+IaI74VjjdmFh4RR2k/az05+/TMG6dCzineM+Gp4B+Ej7f73z359PDrEkj1Lzrscokx/pY80D9vrnm6uBItHGJVk3pFQEo8oUnlh73dA316ZmxK85dTqopfGedzuHEhtt98lb6+gOV1khSNjFPteRVeVNtMrLwDY9Sfn/ADwwxrJtjk8OOIEusimwelel1XTF/EJtRqZGVIilnbsPNjij8+9+2ZGfxnjeIO7eJE7bUcfD14Pe7Az0tx8E6WcrGYrKpGtsoJuz148w47jMmrSLxmE+xAoBCEEO3oL7cfbG0viyzBA4TUIuxA5sEAeUc8c2efTJ+qERZYzpPDEkgBDPW3a10BQ68/71w1BKaeGXTQbQU2A0rFTxYI59b7HnOaATxzy/lpOHADSDoGPFXfQE5ldpzp1iVVAjqyhPVvX/AHt7Y0H4iiw63a6kxEbgopeSOoHYfxk32wxPEhV95DB65oA1/fHc3JH4kAARgfDK9VPN2OSMys5YeayRQBJ6D0yAJ2AI0QV1Jsm7Pz/3vjox86+U7xR/vnN9LbQrTrQJvj3H2xJCQ3l+HoOKJHvgHp7Y0y3tlLJ+pZ2r/Tz+2OsyjSsnhJuYim5sVf8Ae/2yBBB+WA8dNIqFhGGNFiLoeuSzvfOYBhhhgGGGGQGa9IjgTshIKxG2UngdCDQ79PTnMmaEvb4IVw5J3AtQPHH84G2MR6vYdVOwMShP1GsHnovIqh75njjb84NOHSQNQ3DaffqeBzk5NS4gTTGNEEZs+Xknnk386+gzjVG6OACerAfCDfTCtCuDHI/9Sf0IoVVoijd8nk4k06OZVjNoVAUyqNwquldOmMZQhLSMZGlkbxFjeuOD6Vyfn0yyLBMsjSCViOVQABmJFk7ugA29DffNIyRzK8KxyGth8tL8zz68n98eeNEhap4yd3KbKIINf84RLJpozLJCqrInk3oDu9xf9xj2NJFLBJFIk7qvxgUCGv05FV9cB4opGjRp0YRAD+rasgXmr9aPbLoJo2fSRtKjq+0wOwCkc8X3PPXMsMYj/WmgfUwgeerVVa+lj6ffH1GoIli/Lv4ipGCygEi6O6769T98D6LSGMxImgkLssW8x2TT+l9jwB19PbIwlVCaeeAE6VrZZP8AyLAdr/cZi02peSGRhDL+Xmlbc6GizdRYB4A49qzWICqyFmVldFKSWLCkkDi+fQkH3zrKPQ0TRSyCRpHdd1DzHlhQ6AUO+bYWLU8yCOTaV7AsAevftWeBIx0lrHqYXJAaNtwuwfMe/wDfPd0MpngR99llsqB71dE9fXO/x3XPl0aMrLpIJSJJCrDk9R/9un+98sur07TSRO/htGLIkWuL7eozNp5GGpaOSAuGUFmLCgQ3f6V0/fCXXaZ5pPEgZOkRciwxFGvXoc1bnrO6bXxOYpCCJFIa1dqA+R7emeR+JaswaSKIeHMZWILKDvXpzXyGerqPDYsIlWeJVtlDbi9+nPtnganSszGYwzkJ0QqKA5sdD065jnZfF40hiEUC6iURywqKUlvibrdDnpQ9sgGmh0yTpKUUSFlIILKDY+fbqPXNE8SeFKsglWcgWCGJNdaFcf7zkIYvCVULBmEg3DgN7UT296zz1vTxGCYJGJHU14jOwJugOCL9T6dvfORiJ9ZteJnDUDGj0rH7+gIvF1TSPGocKFZghYEbeeSD9v8ARnpaAt48DLIj6eypmjU2RzQBPTmx8u+GlPxqSBdBpdRHJIix8GOuCfS+5FHkH0z5uSY6hgynwk3bmofCfb6AfbPsP/VUcdIqyBBItSSbd4Udel8fMDPklfTHTWHUMQUAN7h05PYjrjkji6lpo5Vld/DYrStZAPaz6dcyCHeoPiqWY0Bff39M16eaOLTqgqRpCLDKDtAPb9/t75DUGnVdl/8AgQKJFnr785mq4AqKhCglRuIkawfpiEqhNBW6EH55ynQNbCmv616YyxrMAf8A46IBcglRxxfXrWQXikcQSLVAgLVjqBd1jfhzyxzMrSPEQCFtbF+lEYgZSZipJsWOL3c9/T98XSSyRjcpe0awFsV9fplGmKR9PPsYp4iv0onnjk9jnqfh2smjYl3XcTtIsktzd1fHUDn2GecryCSdYUbeq8ECiPcn68fTvWaURk0xlkEhFqUVeAeB9+n7dso3KYrSSCetQVUryCQovrxyenHas9bRapfAvUsOeGcnhgL7+vI+h9s8XTOXR0lgIkCkCReRt9Rz637i++bhLHDHBAV2LaFw3Ia+SfQ80ePTNyo9xpJ9jsr+U/ASbDf9V9sNMUllCO8ZnXzDZXC/vwbr65KNwY40AikcLe6NqXpVjr8voPXPPDzGdRKJkAmRLcbty7j7dOma1i9PTlikinMqBnMo8wPAQe/3Ocb8PWOVpAX2k06lqHt8h04yyTLRjVaEfk8hsivf1zdFAwkdxIxVj0NYXjteNqtMUBYBYyu0u4AN0enua/j6TZDGzyNtjXgRhqp+LBo87u2eydKBOAF4o2a5I4sf2yOqjTTL5TKS1oABZBPNk5nlyyuk42zXiEzLqKaLwzKygFgCBddAOn3PzzHNBD+ITsdZK9qxJXcVZQL3cel8D5Z7bQokTSM0iyraoLHl28DjoDX85i88Wn2rA8KAnxGDhwRzZvnoeDx2yM3Y8gqYWkj1WkikRZFpiwDrxVm+1V1rn55gaGB9ezImyM92JLEDgjgntfOerrUgljMu9JClkkxBWF0CPkL4I9s83VxwzToulWVHKl9483VaKgd+OfqcyG/EJIGm0j6V2liQBGVmpgLFm+18/wCnMZiaWFtQJEYpGVcsLHI6dOvP7XnW1XiaNoTp02R0isG2WaNE+poH55nD1qd7FURk6JR4rob6nA2apB+HxiNYIf8A4QJjZc7jfY8c+3TMxSORAEglTULGtJV7+nm+xv7ZfSpIkMmo8SWIAK4VUrxkDdeBQrjH14/L0fAvVMCd4boOD8PrVj0rpeYqxkmGqhnXUTeXdW7oaHoR2+uZpmillLICu48igoGCanaNjjfHRtR5dx7EkdaJ74eLu0/hISoBBIr4j2v1qzgIyOu0OrKCLXd6eozskTtO0aqCVJ4jO4cdaIuxjaeMSODJIig9AWon/GTEm1Su0Wf6+bGAlcE2OM6srjbR+E2PY51kAiD+IpJNbObHv6YSRvCy7vKxAYUex5GQNIdz+IwHPXbxeRzQChUySnxWKnjcQVPYnjnM+AYYYYBhhhgGVmCiQpFIZI1PlJFftkst4QtU3oXY11oL8z0wFl8zlvLZ5IUUBlpoHg0sLMbWfzL5uBRIPGV2ghZpWKSBiWZwGVz14FfL75xYDqJ+PCjURCRtjBhQHPfr3r58YU04WZdPFphcYQtTMLBs3dd+PbgDjOzyL+a8jCceCELBdoHFA/avTnI6chJdu5l3jkx+Yjof9GLIhj0yHwyhbjeGsP3r+2VGjUbHRYVk4VmIX4j7KCBz8/nwO/Vmi0s0Rb9ZzFTeJZEZNiqI5rg4sZ1L6eGDTg+Vi11XmPa7/wDrxkAQdx2SO+0l269T16ehy6KNJIIVeRl8Jm/+FfLuHc8cVxk2aOewEERBJFWb54BJPYYeJ+iF+E9jsFEc9/r/ALWE0PgKm10cSJv3e3pz3sHA1aMy6fW7dLNKpZDbKKYr14v1AGW0WoWPf4xRFMo2shINdwD12/8AGYYzPPPEWllBUKquLJC9P2zW8cRjKNKkkSUARSsq9SRz16ijiDdAiPqEbVONtWjk81yvryPr24vPU/CZIodTLp5JBulIZT8K7qog3/P88/Lx6pEKmZG1Sx0oDsQoHPHyzfPu8LdHAJG3BtiMxjjNkdB14UG7++b48s7Svq9WTEqssTbhY44A47enT9s8rURqzLEupKhuhBO4mum455+o10vhRwagyF1JYvACV5rqOOf7dOt5iOrXUqyHesKm/M+5i37Z15fNKx9HsSzEM/h8yFwFAbiKu/oel1eZ9Rr5UeUptDBt/iuoBYgg0P3+eZS8UEcjQylo1YGJbJXr1PT/AG8eOZbJkBKsxqWrBJ6jkkXVZm8tWRZ5ZFMKx6YyzSHcRs7f/U9up7d+cZSHijmZEgVJCW3UaPNDzAX/ABmJpZNvgRqRFIQGDrdNfX+Pvl01A0enBRgPIyMCNx/wAenGYadkn36oxKpiknk3Dbt4BNEWO1d6656X4THHu06Rrtmj51O9qDCyqqPWunYcZ4SakjURpO8kblirBhQUfL52KPGW034iI9XHMuxiSqglaAI5J6dOvHXm8zasbf8A1KXk/FzuZhH4Q2CPkE0LHXpz1zxixdE0rgMovbKBdLx09R3+ucmSXU6u1IlklJcr8PqehzskJfT6UMKDK4UBu+7/AH0x6jsB06ahSdzKp2hl47jn++JryiSH9MDcpBAriuBXtx/fEUPGxkWPyIKa/MvpfpycvroxFodOS5aRhta+QB2ojg8fPKrLEhfTOUR2ZSb43ACuvt/OSRjtKrXPXjKStGUuMtwa5Fbh68fTOQb28SOMhd6ksD6KL6/TMjunikKlkCtuIVUBtmN9h9Mos8pR1AcMtAtXTtROT0zMNm2rV7FKCxsfv0yi7pXeNkfhtxQDriD1Ifyw1IHlcKQBZulNnqOvOUYs0nihgkanxAjMeF6gD6t2zP8AhTbBKZFHg1XAHFnk38v75tljj0jJ+X0sszBNzK/JU87SPpmxp/DkMcciB2JeXzkuL29aB6Hv++NrCNS3jzy7GR967WN7rO1aPSrB+QOJ+FPGlaeQsX3BrZh1Is3fXt9s48q/mUg3eApPAU2D3JofU3xlQ7TGOCPUQV4aKi/pnkckE9Olj/nPR0mnl1EialifF3Esd5FgdK9b/tnnKgCB36qCrqooMBybHHYf7Weh+HTfmIS6ORpnQhU2Enqwu/QdPpmoWNmmnn1LssyOkZIA6EC6JFDuKI79c9+NfIBd++fNwxypOrRm0D7xW0mRR7315P8At59DoRSeag3UqDe2+2c/k6jp8eRUjiubr1xZVO3coBdRwCavLbQOc7tNHpnG8t7dNeN+Ib44/D/MGIu1blWz3oV2+d581rA1w6f8zGDDHy3Rzz056txf1759jrI1LAoo8Q92Xiuh5r0J/fPnvxGDTtOv6TAox+In4a5rrzyevNXxm+N67ZrzTE+s08kjFTIFFb0IIWlBF+nXjMEpjUDwXeYsCHJHKMOFI9L/AJz2n18c0QhQbowdrSr5SR6UT25vt7dM8qGaPTO8e9NkyhibK04HxC+eoGdHJ47aWWXVmNNmxWahI4Fd6JHfFeQCaaPU1axgDbXUdOR7YJHJp5Cs3ma7A3D4q4J//wBWD0zRKIxq0uFonUKzsfNXYsaFDt2zKJXDOB5pFWgaHlCkmia53dAeKujmbUySzvKz6kSAUwo1fyH+1mj8U1Jn1WxI1jjitQoIYX8xx1sj55knSWGNY5Nu09K7G+frkqobG27tpocE1nCrIRYIJFjL7/HcoXKK5Bosdu7pZ/fEZzt2uA/IAJNkAdh7ZBLqcZ+XJ/tlNOB4qgI0jkEBVW7J6fPJkHdXQ5Ahx2rizZ/tnK+w9MrqI5FnKyEFqskG+2BNzuNmrPoKxMqkZlkEaFbPdjQ++SwDDDDAMMMMAy8wdpXadgH4JodePbjIZWGHxuPFjQ7gvnNde/yGBSeXxSGkleV1ARSRxtA4/wCsXUKilCsiNuXkKCNvz9/lnESN0YFmEgraALDG/Xtxnd0ckLAqiyDbRur7dP8APpgNqHj8VX0zMvlA4G3naAe/reNvjMKq6HhSE8lWTfN3zzWZxwptR16+mMsrRyiSFmjI6FTz0yjTDCJopVLSRhFssxteOKPp198hMAzssCkopPN2SPnQ/sMGncsCeV379h+G+/GWZZlhfcPDjJ8yhqBarHH1yhtQ35fdptQgdkTaqhiBGxHJr19frmeTxHUTGzVAkrwD2H2GPHCzRKTFal+ZFG5gOO1+4+uEYgZmd1cRqK8vUk3R9q/jA6sktrKCrVaKFIsE30X69hlog2nHhTRyhAhZ1awDzxYHaxXbIxpYVIWeQsbIAPUDgD1PJzkCyzPI6R7yqkuD0rpf3/fA5BJIinwZGViVIoc7geOe2a9Fq5Y921V1JG1vDlPAC83X7DnpfHpEtFHJuvYGIJiVTwOa6n5Hr3yM/iAK0jAhxahWHHboPlgal1Eza2PxJ9ylt7IrMAD6fP3za1Qx8JLEJrlXawFUrA3QJ68fLr1zyohOknA8JuCWZTwprk+3P74NM7jdtCjkblFde198B4VImMIgLahmG0E976VmiLUgukLKvnPn822z73xixwxrAqaeWWWeYEFY+h56VVn1xGlDRXAqKAPMNgFfWzfXrlReOZPGZ5FbeAUpK63d+/8A1nZJDro2d9O8kl7V2UoFDjt/GZ5H8LaJBE4LE8dbIo21WayTMFkdRGlNQFc/Uel1ltqqakmYopCjwR5m2gbiT6j+99sTzeEzIpVAKJPrd8YMjyRK6p4cYG1ueLA6+v8A3mnSyRpC+4Ekx8b2Pr/T/vHOZkEI5o4lkVUsmwHvzDtmjRASzgElYerR7uW6ccZncxC45SdysbKCx975oDjNMepaJ45Eb9JTxHyOOa5++WDNI77PD3KRHwP/ALC+OM3s51X4YQ/iNsr3JoG+ew75PWXLqJYowwZmBG5to9xyPUf95rSTWRfh0kxjAUIFIlHXaa5HfrWaiPMKl9kCCUop3eHKwXb1uh8szujxyq7KxFAruHBFf2ykjeLCOF4O5iByBz09Bzmn8z0u2oD9SxZPQHkemZVhgcxTKyimu0IF856IgVJ5XbmkJYKp62OF6V173kt2mdDOgZFvzLXAayQB1r/jPS08MWp1CPvXxXQNyprkc2Tfz+mJEQ08WnGgldlnPnCiuBI3NX9M26mYlo4pzIwZQhWiLIqtoHue/XnPSTQlPw7cskMaEmiOhINbhRzCN0s7BnLuAWUA/a/SyAc3hHnDTw6WdxJqE8RuVUi7PHc/Xqec9ZJ5wkqxRN+pTeWQOSD6D6j7+2eUWUaFv0pPDcgABd6hxwDY/wBN+2N+YehIqNJIvDKqFQhvjnv+37ZNV7EB8XULcjhHG/cF+FunB+19/TKQhtHJ4G5QCLRWFBQx5IIPBPPvnjaf8WOnEg8W3fkALYuz1PUHn1Obl1bflYWUoJXHhtFFyavqL6/v19Muj29PNCGMIart1I9Ls8k9Pl889XQysVCsrKUobS18XV33B/jPkY5DPGjws++PyBVSgw9/Uc3djqc9qLVrpPIZX3Ci5eUiwB1vn36emW9xrj6+oU3jjkdc82HUoCB4yy7gGsHgDt9DWaE1aFtrMqm6osOc81mOlhtUFEbbq6el585+IO+o0p02tDJvKeJtBIoNtvkULOe1LMks/gllB4K2bsfL04zzNW3jaloiXUuhI3mkKjqKHPUg+v73vjC+PGkp9DMqMjFN+4EADaCTZruDQF9eczaVIpJFVNM8Sk7wEJ6bq635h8XIzuu/EI9XppEWRYSa8LYpUWCPLYNHk+ne8i8bliNRNLDqFa6d/IAtWWBPHPHA7dOc2514iaavxR4pZFhUmml6hen+/XPQ1eqhOg1GkkkiYbd0cwX4mWhQ69gOeOtdsYmtaDBPPGHC75CfMzbqJr5gcc8Wc8idN35dIiC6gh9qEVzdn719Ml6ZA1Mm2QK6IGPiEbaAYdAB2+mUk0niMiwMH3BA6ILZWrnih0/nrmOYh9QdrXuPVjlY0k00tGTwZOzbj0o+nUHIH1LaaHVEQKzRqNh7EmqJ/wB/bFnlhMR2BmeRr8zXtUcAfPIyuHJO5mNk7m/q+mWm08SabxI3MgHBcIQu49ufYXkGW2Vr5U9R2zVYi1MqagK5ZipaWwQT/Ua598iYGSISSeWwGRWHxjnkfbGjijJV5Zht3AEAEk319Onz74E3K/0WBQsE3zXOIceZBHIyWDtJFggg+/GcjdkDVXmFGwDkCqwAYbQb7nti44VqYgWALJrpiYBhhhgGGGGAZ3OY8kbRMVcU3cdxgaJNQv5eJYQYnA/UoHzHsbJ/wMgqIQtvRb1HC8984pZWVx26cXlZZjJKJixeXdduoogVXGAPpHRYTuT9VC3xDiif8YCWtimR/JW1x1TuQBfqciu2m3Eg15aF2c6iPIwRFLMeAALJywP47hpCXYmQUxJ+Lm+friom5HbzWoB6WKuuTljHGXdSvhBjaFmJ2j06c5BFJKg7qJrjk/bKNELPPNaSGOUgsTdbiOQBQ9h9clKDzRKg8hGNmu398RVY7iisdoskDoPXNKxqsUkgJZ2XyKgsD/yv5D++ApnKgxacsUDWhI5A5/v/ABlxJCdUXpZfEULsIPUjk0Pe+B6jM6SEeCYxuKXann3PHpjwTiNvEO5KfeoVAVDfX2vA7N+X/NhpY2iWgWjRfhPoLPpz9cVpn1bwIxRPCSgSCRQs89f8ZWbVCbTzpHGxDOHZmonjgdvcfXK/lm0usjjedfGYUfICoUdOpHWh6cHAwEyAOwYj+lqNXfavplImESBnSN0exV88f2yYZmmQqb21yR0r1zXqFnGojk1IMfimmkKbewvj2vAmpXTTWszVXleM888cj5XldNPpQ2x4PKSxtnIoVx9RkNMVeSOJzEELHluAL7k+gxIi7hlQXY83pV/9YFtTuYs8yinvw/D+EkULxVDRjxD4ZFgstjd17XnIZ3RgrMFUjm+eLuhhC5KB5EaVIugJ8q2caKM6W8YRhG3lBJBoj3AGK5j8JgsjSKeIxQFVV2PlnHXcu6NnIQAMJFqieOOeeMAgKeCsjsASUCrwSfn04GNoCwAjJO8MbY3f06dcsShkLINun3EBTyQDkV0+2IOZRZHKp5iOtf29cAIlEXiMwAXeSg5u+nJ9BlGldNLqdG0m79OIcsSbF+vtd5qd3OhljUAmGMrKR57J9/SyOnphp1eBSZEmEddJKAAPqPkc06zTDSab8tpFWWTVIPKALUUSb9z269M3nWo+cDnYy92788Z0ySgD9Rum0ebt6fLOskaCNt4cm9yixX1/3pnNxkQKFWkFXQBPP7nnOSxP759J/wCn08Z4ZZFtFUoWZwVsHiwenUZ4IdbLKvlTlQaPfv65r0Lpp5QskcrLIgPw117jrx75qXsx7/45PNqdsaSVpwRsiC7lYVzdd7v9s8QlzI80O2QFiaLA8etX7jjNckLync0gZSjEjq1cksePTvmV9XHw6fEzMNyqAzcjknsavjpmqmYpL4jwF5CsMMnPhNfUdPcXhpppImFTReSjvcE7feu/2vOzSDUQxAGNwW3bWoMvzN/L7YK/h60EwxF99gqNwNjpR5PT19ciieTUCGYRlzGqgqzdSpNE8cHmu4Pzx9DMkWpI1MkgWNbRym9dt2eOP574sy6kar/3cYPi010OD2A9+P8AOJHqIgsfjBiwNlA61Jz37jm+PYYHsaTfr9IsaSqI5Qsclqw8TaTZBvuPvWb9JqNJp4nVIJEWVvDjZqKsOhoeg559T758/HNPApeysQpwViJAs8A9OOua21U5uZaVY1KhVUtXrag+Ucf95WpX0GnI/SijnkkRI1QRgDdVdGHXkf2OUj1CnVSDeUnLFSiIWBNWDZNdL9M8zSSfm1i3uYmAsbB5nIN+3P1J5ys+sjGnkAvU6lXvwy4Yo3W+b4qq/wA4sjX2WMimKMxzeEFYqfEBB2DgC/TqeO/yzbrYdN4cjEFyPjqhfINnjn7Z5KzvNrpIyhUBE3oXsg1xsH0JOehI8GnK6nUSgsf0149O5PJ9epOa4eMcvXipNDDrVEmmeNbQRlQTuFEWRx7j5nMOrGo/EJY5/wAwUMKjyyWtAkjqfcf7zmz8YWI/iC6rV6pJE6Roj/B7HvmP8Q/GPFhlWFXXTsSDICCC+08fL/vJRhikfVQahVE0jRstSoeEUHgnizXPU98o0CuzIZdqkUWfZQBtL9uQPX1yWl1ZaPWLDplWPwa8u0lVsWeevF++ShnXWTr+b1ciqjUHoDYvQEV/YZhGIKunndZVLUCB259aOVMMk6tLDGSqWSelDntZ4AGFxQ6WmhR5GkBtmPCjtXXm+uN4xPgK8pkicANGrUQAx8t/v9cg7KnjeGU0+x0dYnVmPJoAWSfZvlktRKfgWwnYdgPQeovG1crLrWeMspV9ysQA19b4xZbmdFkuOVVO9pWPPUj5ccfbARtpZWa2TiyOPoMRFs/EF7i8InCPbIHH/iemK1A0GseuQVDRmGTdGzOSNr7qr1sd8S1KqqrTC7N9cA5WNk4pq7C+PftnYYzIWANMFJHNdBeAodgCATTcEeuJlfKIyGUhwRX/ADksAwwwyAwwwwDLTKhSN0MYteUUkkEcc369ePXI5YKsbPslDFSQDtNMOecAJmRPLI5jW1BUnaLHI+oziyAwsjmq5WlBN+l9QOudj8Roiga0vcV3fvWd1BRmZ4VEcb8iOy23npZwEmgkg2eIteIodebsHviqSp3K1EdCDzltWY/0fDP/AOMWAxO02eORx+/zyJUhQbHPa+RgdeR5GLO7Ox6ljZxbINg1lbRS+11NoAPJdni+vTvzkiKA5vLos5WVGkJVXFDaq0CK6/2++TjveNo3H09cqVeQKsk4MgIjVHbovzPAAyBBViO49DjQ8iGJqPDd1ogqfQ5SHVGNk3rvjW9yXW4EV9675F23sWIFk2aFD7YOjRuyOKZTRHocaNEcsmm1gbb4JB+Eg0t+2dmaJWG1jJw3nXykm+Cbvt298ykkgWbrKwal4LG1XQkEo4sEgED+5xop4X5fUPBNSsCVfoa+WAlEuoLahvF30N7k8VX8CszZUOhkZVUIjgLbeYr05/b98aKIW8+mRw8YJcEcBqB559ryUbMtlOCePasbTqp8QsAUUWSd3r7evTn1ykkXhsdN+oHJDMGHC0CexN9euNCLsAVnolT8BB5y0yyeHF+m9ISEBFACzf7+/rmd4kDRqJQbALHst4/iyBGkQBRe21JHUen0ONHI1JkC2iFxW5qI/wCOmWj07ooJkTzgNsvlh6+nXMiyMDe4g+uP4p2MnfcDus2MQazHAJUBnoILIoWPUc8E9MPzMYVlEbyREgW3HA5IoH3GYGYt1/tjwK0kuyOMyFrCqOvTLo1mSW6hm2oR8ZsXXv1zkglleXVKGpK3sW3cni7PPNHMod6rxG54q84zkoF/pHPTJoC21rVjjTy+K7OQAzEk7RQ59sVirKtUpA5685ytyk7uR2PU40VijjVkOpMiobsIPN7HnjrmibbHqI4mMDqovdCevUiy3z/bIIoEbPMCwdDtNkc38uf+cmeUThuLA3NxXtjVehB+LahG/SNScES15gAOnBoj2Od1Oo02o3SeAF83KjihtAvpwSc854zGzKxG4VyrAivpjTshcKhJAAtq+I9zl+w0PAPDLfAt0j9QSOovv1Ga42eQ/qmgg2hmJDVdWOe/es8uIhDvam2kEIwsN88uoYNSszxhv6eA1d/Uf85NRs1Ok8DUE7nMNXu+KvQ/fi8jLHpliYvIFk3k+CBYA549etd+mNDrkRpLRacXtfzi93c3fAJzTp4/w16dG3yk8RP5enYnpz0+ua9EWLPGLNhgKKndtAXdXr39+ntlfzLR6YiAnTFRv4UXusVZ69vf/FD+HW2998I3hXijAscUa557+mcj0jDUIAhMEtgDoxs1yR0INdRl7GqL80syajyCZ3Dbyw5DdRXQ97P3rNn5mVtRBpoHhjYncQi7Yz0o2Lsf4zDHFJo0EGoYJFJfhuAKDevNAcCiaOanbSxaOOeRYzuBZHXzL5RXmUd778f51DWn8PQ6cRNGr6rVMG2liBVX3bpzfHoL71k9TqdRHIjqIAefOXDKhsgBiQKF3869s8l9dql1ceqlUPIwtWLVvs1R56CzXGSQAhY/hIPnsi1c0RXfj39cfb8Ro1suqOhVG08QSvj2X1A9PUjj3rMc0ksSxanTxpBp3PiLEpuq8p56n/n557H4g+nXQkCPUsefDaHgL0on081/PtWfPa0tA6QkblRdtsDRN2a9uR98nLpTT/l5ZNM8TuXY1IrNZFHjkgDp/bF0QlLxeHCjOGJjJ4JIo/X5e+UjhnMLybY1CDcB5Q4A54vn5+2ZoUtkmZg/nt16tQqzR69f2OYCy7HnYRkupPlJG379f752TcWjWWMRoPQbb9efXGlRPBR96iRhe1RxXT78Y+7cEm3ionFqW3AAnsCbPN39MDOI2kidwTsjoUbPX6V2OHiHZtNG+rd/llZQZl8SNBFGeBzwSAN3J59DXvmwxyLL4gjJUwndEpAKJdc2D2PX3wPOJMsaoEjXw1PmHBbnv6400BgeLeyuzqHKDsDyB9sdy2jmRVZHKBXBoMLIB+uRAd5LsWebJoZAzCN/DUBYqFMxJNnEaaRpRKXbxBRDXzx0x9SzySGaUqXlO/ykevt0zkiqbaEOYwBZYdDXP75R2aVp5HeR7YnvkMoY3EQkKkIxIDEcEjr/AHGTyAwwwwDDDDIDHTd5iGogetYmPWzr8X3GAQ7fGQPQWxd3VfTGSpZVWaXYpIBYiwo+WdWOH8o7tJ+ruAVBfTuTxX74jNv8zEluAOOwFf4wFU1fAN8fLKKVbchYKhJYMVs8A0L98qJn0UpEDsrq3JNeldOfU5mqut+2A0cbyvsjUs3YAc4KSoBIBUnkHvg0boqMykBxuX3F1/GJgWmczSSzPQeRt1c9znXQNI4aaO04BANNXpxgil9K5GxFj6tzbk9B+xPbv7Yumdo9QjKiObra/wAJ9jgcj/TdJGVWW72nm69R6ZXVyRz+HIgVDsAZBfUcemQdClX35GW1YkURIZPEiAPhNVWu48/e8BYlSUBTSbQSWvk/c4m5mVVL+VboE9MTDAZRv2oiMZCa45v0AGdUtDKCV8ynlWH7EZUrHI6LplkEhIUKTZY/P+M74v5mdQYI97LsFeQFj/Uffn5YEEZlPkYqTxYNcd8ssQdCyte1S0nlPk5r9+PvkAWQmiQenBxkKBXDKWYilINUb6+/F4HY3VRRjDHcDzfQXx9f4w8RgAjFjHu3FL4vDaIw4csr/wBI29f8YNciNIWFggEE8m/+v7YE8vP4eyNUUiQWJDdgm+K+mKQil0QeIAb30egyYBJoYFfDvc5JEQsB64J5r71jxIk8e0RhTEpZ3B5I+RPuOmRW9hABIHJ9v9vLRalI9xMKljVUSBQ6ggdQcBBIYvCZUCup3but+nHTE6KfN19MLXZRB3diDwBnHUpQN0RY98Bt0SvGyqzAUXVj1PzHbHiCySxotqW8rG+tn6VxizpsZUBjalHKG7vnn35r6Yu1mDNx5eTzgOzsm5CpUjykdCOeh/59MmGJ43UL79MGVlrcCNwsWOox9qElU5BNKWNEfPAU3Xb6ZwN09RlzvjTnUdBt2qSaFmx6e/1yDqoNo1rZAvg17jAd3LuXcA2PllRE5ZAJoyCu4ktQX2578ZFGUK9miRwKsZxQzqas1zgViJEn6aF37cXf0zm8GioClR09TlyqbFnknBkZT5GUhruhRHB9e3Q/XPIrIEZl2hlta7j1wLJqpUiCwOyBfM1MRuNjNf5zWtCqjWSEPYssVUccizwf+cyTRwmT9DcR6Hj5XjSrt3BElOmI8hdaF9L/AL981o5I0sjok2qaQkbgS9hSRfJPT3y0eu5v8tEjqQ0ZArYfXm7yUUMTAb5NjMNyq3w187+2dLePO0cEa+GAwQMOAOT3Jo5No9KPWzJDJp28NQoAVXUDevHl6Gj9uvbGXWLrYWSoY5nbxEKtySAFAa/v3+VZl0UCT6WdGXd4Z3K6Dl+Rxz0uutcffKRaVIgsO8KGLFJZaAtTRFc30HH+c3NFvxZpIBJppDLOsC8c0FJ/qI9BYrjvnkvpzEQruu+QblkDeUrRsdOt8Zb8WlEskbeYEIAysK83PP7ZiLO8QHibgv8ASew/0nM29i0ULyzKkgbcQCoui46UvHJ9MvDAkBMzPNELI2AVIEPG4dL9PociIzpUXxV2yHzxypJyp69j/wBZyKaaHURmUtKAgATxDRU87eD055GQciaOSVQ7hR4gJYrZruSOnHpj7Smmm8Lew3+dl/8AjYDp9bxYolmLCVyhCnYp4B56fc44Qadzp5HjY7wGQsdo9Tx9vXIGrSPoI4ogx1LOSXKnrXwCut+uQkIinjBjMeygwJ3c9+P4ysEkulZFk3qI5N6FQKviyD8qrKSaeOV5JHiaFRHaIHUFvTr14HNd/nmhk1BHiyFwSW8yEUAAfYX27dsSYp4txs7WASWFG65/fFjCmRRISFvk1dDNGoTbudi0iuxMcjIR4gsi7OQRYwgGjIzFRyQAA3f1sfbGHitp2p6j3glN3U0ea+/PvkmC0NpJ45sY8bBY3smyKABwEZicXNOoiIGnVRG29LHh2SbJ4Pv2zNgGGGGQGGGGAYYZWNPFXZGjvMW4Ci+KwOBA6rXlrgsx4JzkiqpGx94oEmq5rkfTEs1WXZQLVJhsUFhu4s+leuBHKRhp3Cu7EIhqz0ABNCzlTDJJvaV0UmMzL5lAbn59evHX2yemVGnTxZRGl8sVLV9O+Am07Qyup4si6I7f7WCpaksQK7HqfljzLEFDRylyzNalNpA7H6/xiyqU2r4qyDaD5STV9vngTyskEkUaPItLJe031rEjYK6syBwDZU9D7cYWXYmRiSe55wAuz7VdyQo2rZ4UXf25OLnWILEqNoJ4HpgaoVd98DqI0jhEUsx6AYuUkMe8+EGC3xuPNX/jKJHPGsTGAlJWBS0veR2B+vT5YEmkYvu463QUAfbplNUYjOWg4QgGh2NC/wB7yTcufLRvp6ZyhtJ3C76d8Dp8wLEqCK49cZI5AysFIsFlJ6ED5/I480UcSCpN7sFYVRABHIPPBvJhwQdws0ACO2BUeLrtQxdpJJWtia3E9/8ATkxE3h+IVOzkX7/6RjQJ4rrGjhHaxudgq1Xr2xNn6Re+jURR/wCsADFVIG4MeDR7emNCJHJijDtv6ov9VYsrBmsEkUByK7ZxXZVIWhfeufvgOoEbNvNOnQEWCQen985ErlmMabtoJPlsAeudkjXc5ifei15iNpP0xWk3Rom1RtvkDk364HKG3rz2GMhaVgrbnCqaG7oACf8AnFV9o4VboiyL6/POoAysK8/UEsAAADf8f6cDQGEQ2FJEkBKyyfEQDYIA6dMgFRoyxkp76EdR8/XG087xLIgZtki0wVqvJHljQodhgNF4e9fGLbO+3r++DRugXchXcNwJFWPXF68AY8yyRt4UyFHTgqy0R35wGWWRohAhamPKr/V8/XJBiEK0KJBuufv9cGYs1mr9hWcwKPEUhSQhqckA7eOK79+uPAxI8NVjDE3vf09OeKyNGrrj1xtrKocqQjWAa61gcU0bIDAdjjmWzaoq8k8X37ZxNqsC6lh3F1nDtoUDffnLAzyb9pIsjg+4za0olbThdQLHTcTSff16/U5grdQF3lVldI/CLXGr7qHr064CE7y1njrfrmjZHGrMJP0nW08wLWDXIH89szBwSocUgPO0C6zVt0xiYK7BfiG4eY9enp1F/L6YAZoIUOxFmVrrxFplIPoDwK+ebNVOkn4bpoxpxBA7My+eyTY79u/7Z5oRTCpZJNxet9+WgOR069MdRGyvHsBZB5WQE7vnzxx7dhllwWOmiOrjEsyFXW1G/cCboAm+L9T65J9NsjU+HOrKSspZfLuvp9scmCVli8ExNvCli9heTZ9L6fbNKyvO0umj3eE1KIx5nY+q1V/49clgwRKaZRsLobAILFvYdvfnFV9kJoAsTVljYHyy0kJeSMaMSPuJ2nbRPT0xIg8LtMxjYowDRyG93079P7ZAzssyhVeMPs83lAU0BVf/AG98lI+1TG8ZWRSBfSgLsEet/wBs6oMRKoyEt/WDVCuev+8ZWGBJ5ubIkBCEvR3V8if85Qb9oEUEiyCRQZNwAArmrP8ABxdRqmlYMyncBSncTto9F9vbIlHD+CVYODtC979McoYy+neNN4ei/PFcH2r6ZRo1G/8AKwnxkkEqmkoWnm/ayL4zNNIxZCoYIgG1XO6v2rOwQrLqAjzpGn/m3QZ16iRtjoyyLwDyVF/seP3wIyXus1zz5emdRf03by8V1POEa7iU2qSehJqv4xUUuwVbJJoDIOFiSDZ4zmUQuAYgQA5F/TJ4BhhhkBhhhgGadHrDpBOFRW8aMxk82AfTM2dFA8gHAs2nfxNsZSU7PEPhm6FWfqB1yTVwACOObOaIjIiI4MSqrFbYA9RySOpHvWZmFMRYNdx3wHQJ4bs5sigqg0b9enTj98pBqHVBDEqKXtWY15gSOCT24zsIV0LFF/THQA294TrEXEqBBF5QURuenPXn69MCU8kkszPM5eQnlibvEFX5r+mcykD+HIslXtIPKhh9jxgEieBM6ExybbFqbU+4OTs5aSWK/wBOLgpR3m/N3IqvtkxR4rzE8G6GB2NA6t5qYCwPXFvygbRYPXuc26ceHq4zNANiAK+1N930NXV85le5WeWlFtZCgDr6D0+WBPLKrppzKI2reAsovykdvS/8Y2kjilaQTOEAQlSWobu3Y39MTTx+NKELBV5JYgkKK68YCxsok3Sqzjno1G/nittvykke4rGDfolNi3e7d3+WJgGMiqZFDttUkWauhi5dRGqRyEq+1huj6WOvXA5qBGHpJN5HFhNqkDgEfP3AyRFAc9eox5gN1rto80pJ232xCxNWSa6X2wGWPdIEZ1S/6ieB9sCihN3iKTuraL6euLWa9E0AjlM6k1GQu1dxUn+o9vb6jAzAll5fhRQBP9s6BDZ3NIRx0Uc+vfEJvtlZWUlbVg4HmYtdn1+3bAWSP9UhFYAnyqw8xB6ftk8dUdgSAaUcnsMWuLwKI/iTqZi7gkA0fMR6A4skbxtTqQfcZ2FyrBTRUsCQ3Q1hJIJDdENZ73gcQFjsBA3EcnHiRZGLTOdikbyOWr2B64kT+HKjlVcKQdrdD7HFJ4A7YDu5ZR5UBJJsCj2/bBmG0qadiQ2/m/cf76Z12jdAQCrAAVd7vU+3bjFijMsgRWUE93IUfc4HfFIg8IKACbJBPm9L7cfzisWoJvLKpNc8YuOxLqoCjyLRIHa+p+/9sB9SGErEyeJ0827d26XioSqE+UhrHPXEVtvv88qXjAULFXw7ra7rr8rywIgYhiAxCjkjt8/rWcN0LFWOPfKCWIJIPB5atp3fDz++KzBkVQqjb1I74CkjaFoWO4zVqWqdj4vjOJT/APYHn173zmUDi7GcUkMCDVc3gaZige4nPUmuwziEGRY9rqhoPtNlvU5YafUS6aGZmX8vuMYJahfWj6HrkUDSyny2DZO3kihyf5yiumiV9gEnl6uP/FrNcVzwB98cKEbnUbmXciEmigXuVI7i6Hri6TV+Dqy6RrKRu8MfDXvlGgkeCXWMQAygksu0kmun3yQYSx542g9h0zTp9V4cis5YQnd5IiAflz2+eR8zKsJYlg3A3WvP+jLwwgyeDHJG7vtKtRvvYHbv39MCPjEyvMLL7rBaiB+2EzvLEhMiGmICKKN8Wfr/ABlZymmBSCZikq+ZWWuL4/teZ453h3GJihdSjV3B7YFgNOyRiJpI5wSWlkal9qAF5xyJGVId/Ml7B8PQdLPXrk5ZTKVdiKUBQqiqA9v3+ZOLC6rKpddy3ZXAZkV0aVSABQIZhuJN8genH9soHdol07EtEhLeRRYJ6m++SknaV3JpVdtxVRxfP+cdS8aIzGRQ3AAFblv1+eUKG26lXI3hWFK46gdARk3bc5agtm6HQZd5ZArsiKsctp8IugQevW+nOQClug565BXwkXSrLuYOxIAI4Ye3yzPlztaNFQDcBbMchgGGGGQGGGGAZ1iWYkjkn0rOYYDEs1szE+5xcbxH8Pw97bL3bb4v1rOE2b4+mBWHUSw+WKVor+JkJBIIrt7X98mil3Cgi2NCzQ++LhgU3SQPIiuVJtH2NwR6cdRk8M6TZJoD2GBzLHTuFra/i0G2Bf6SLu8jjRyyRNujdkNVamuMCk+nm0kirKuxyoYCweD06YPK07hiP1DZZrNueSSb74k0niyb9irYApb9OvPr1zkcjwyLJE7I6mwymiPrgXWn0ZBQM4+FrNqOpFdK6/fIWNgAFNZs31w3naFHA/vi4GnTzQrS6lGKoCU8MKG3GupINjjpmcg9xV84/i/omMxxnpT15hV/57+gxXkeTbvdm2jaLN0PQYC5aVIxpoWWt53bvPfy4rjELqd1RoLqqvj5c4pNn/jA6pAVgRyRxg7Fm3EAE+goZy+ntgDRsgH2OAzBQiMrglhyvNrjSxGCRo2NtQ6Aijwa5H0yWdDEG+93gVeNoYl8SMqZKZCbFrz+x/jJNV8En5jBmLNZ65zAqrbIHX9QFyBYNAgckEd+2OglmAgVyVALBS3HAJ/zmcm+uGAyxs6uyiwgtvYXX9yMcp4k4GnjcqzUiHzMfQcDn7Yiu6BgrFQwpgD1HWj9hgZHZVVnYqgpQTwO/GBztRPHpgASCQCQOT7ZzOhiAQCQGFGu+BzNPhoukaSr3y7UJ60Bz39xmbDAYnBWZQwViAwogHqM4TZvAnpwBWAFiVC3wCT/AL9soZHXTrFutCd9V0PTJZ2zVYDFlZFGxVK3bC7b548DkK8RZVWTqWHAr+e31xPEbw/DvyXde+LeUGaJYQHG+oQUVlB53Cu1DvmbC8g2wyMsjjTvUCliBLXoa49ayIEkiyTsC4UjcTZ5PqfvkM6HYAgE0eoy6NSMYoiyrCUloFGpjQ/cftm7TukuqWFFRJZBsVo7cAEVt5PX3GePfGPNL4shfYiX2QUMaNMenYzlGje/DLMBxsFXZ9RXOPPpfBgRrtlTz1RAJ6f3zPBrZYBtUIy2SQyg3Yo89cl4r3YNdOmXYLTQmCRUbYxZQ13wLGL4axS7XKupsBlPHpf/AHnZdbNOqLM3iKl0D7/LOHVSb0dfIycAqSDVVXXIORJI8hgh8xkYKBXxc8Y0kQgYo6sWRyrlW4Pt/fEk1EkhTcQNgoUoGcEzCNk4pjZ45v55BQ7E3q0JtgChLWVHX6msCrPp0IYsENEc+W+nt65C8eWbxCCI0T/9b5/fLovNFUzmGJti8H+oAgc8jj1+mG9HQEhnlMhd02gAj2I5Hf8AbIeM3hCMUBd2ByfrnWndjuNb7J3AUecaHWVRHMFYpuqlC3f17Znxg4CkbVJPQm+MXAMMMMgMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDAMMMMAwwwwDDDDA//9k=\n",
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"400\"\n",
       "            height=\"300\"\n",
       "            src=\"https://www.youtube.com/embed/7HF5Oy8IMoM\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.YouTubeVideo at 0x1094f5940>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import YouTubeVideo\n",
    "YouTubeVideo(\"7HF5Oy8IMoM\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Can we beat $\\mathcal{O}(n^2)$?\n",
    "\n",
    "- Generally speaking **NO**. \n",
    "- The point is that we have $\\mathcal{O}(n^2)$ input data, so there is no way to be faster for a general matrix.\n",
    "- Fortunately, we can be faster <font color='red'>for certain types of matrices</font>.\n",
    "Here are some examples:\n",
    "\n",
    "    * The simplest example may be a matrix of all ones, which can be easily multiplied with only $n-1$ additions. This matrix is of rank one. More generally we can multiply fast by <font color='red'>low-rank </font> matrices (or by matrices that have low-rank blocks)\n",
    "\n",
    "    * <font color='red'>Sparse</font> matrices (contain $\\mathcal{O}(n)$ nonzero elements)\n",
    "\n",
    "    * <font color='red'>Structured</font> matrices:\n",
    "      * Fourier\n",
    "      * Circulant\n",
    "      * Toeplitz \n",
    "      * Hankel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Matrix-by-matrix product\n",
    "\n",
    "Consider composition of two linear operators:\n",
    "\n",
    "1. $y = Bx$\n",
    "2. $z = Ay$\n",
    "\n",
    "Then, $z = Ay =  A B x = C x$, where $C$ is the **matrix-by-matrix product**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Matrix-by-matrix product (MM): classics\n",
    "\n",
    "**Definition**. A product of an $n \\times k$ matrix $A$ and a $k \\times m$ matrix $B$ is a $n \\times m$ matrix $C$ with the elements  \n",
    "$$\n",
    "   c_{ij} = \\sum_{s=1}^k a_{is} b_{sj}, \\quad i = 1, \\ldots, n, \\quad j = 1, \\ldots, m \n",
    "$$\n",
    "\n",
    "For $m=k=n$ complexity of a na&iuml;ve algorithm is $2n^3 - n^2 =$ <font color='red'>$\\mathcal{O}(n^3)$</font>.   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Discussion of MM\n",
    "\n",
    "- Matrix-by-matrix product is the **core** for almost all efficient algorithms in numerical linear algebra.  \n",
    "\n",
    "- Basically, all the dense NLA algorithms are reduced to a sequence of matrix-by-matrix products.\n",
    "\n",
    "- Efficient implementation of MM reduces the complexity of numerical algorithms by the same factor.  \n",
    "\n",
    "- However, implementing MM is not easy at all!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Efficient implementation for MM\n",
    "**Q1**: Is it easy to multiply a matrix by a matrix in the most efficient way?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Answer: no, it is not easy\n",
    "\n",
    "If you want it as fast as possible, using the computers that are at hand."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Demo\n",
    "Let us do a short demo and compare a `np.dot()` procedure which in my case uses MKL with a hand-written matrix-by-matrix routine in Python and also its numba version."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "code_folding": [],
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "def matmul(a, b):\n",
    "    n = a.shape[0]\n",
    "    k = a.shape[1]\n",
    "    m = b.shape[1]  \n",
    "    c = np.zeros((n, m))\n",
    "    for i in range(n):\n",
    "        for j in range(m):\n",
    "            for s in range(k):\n",
    "                c[i, j] += a[i, s] * b[s, j]\n",
    "                \n",
    "    return c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numba import jit # Just-in-time compiler for Python, see http://numba.pydata.org \n",
    "\n",
    "@jit(nopython=True)\n",
    "def numba_matmul(a, b):\n",
    "    n = a.shape[0]\n",
    "    k = a.shape[1]\n",
    "    m = b.shape[1]\n",
    "    c = np.zeros((n, m))\n",
    "    for i in range(n):\n",
    "        for j in range(m):\n",
    "            for s in range(k):\n",
    "                c[i, j] += a[i, s] * b[s, j]\n",
    "    return c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Then we just compare computational times.\n",
    "\n",
    "Guess the answer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/anaconda3/envs/pytorch/lib/python3.6/site-packages/jax/lib/xla_bridge.py:122: UserWarning: No GPU/TPU found, falling back to CPU.\n",
      "  warnings.warn('No GPU/TPU found, falling back to CPU.')\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "755 ms ± 12.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "1.21 ms ± 65.9 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "21.4 µs ± 499 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n",
      "204 µs ± 6.93 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
     ]
    }
   ],
   "source": [
    "import jax.numpy as jnp\n",
    "from jax.config import config\n",
    "config.update(\"jax_enable_x64\", True)\n",
    "\n",
    "n = 100\n",
    "a = np.random.randn(n, n)\n",
    "b = np.random.randn(n, n)\n",
    "\n",
    "a_jax = jnp.array(a)\n",
    "b_jax = jnp.array(b)\n",
    "\n",
    "%timeit matmul(a, b)\n",
    "%timeit numba_matmul(a, b)\n",
    "%timeit a @ b\n",
    "%timeit (a_jax @ b_jax).block_until_ready()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Is this answer correct for any dimensions of matrices?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dimension = 10\n",
      "Dimension = 20\n",
      "Dimension = 30\n",
      "Dimension = 40\n",
      "Dimension = 50\n",
      "Dimension = 60\n",
      "Dimension = 70\n",
      "Dimension = 80\n",
      "Dimension = 90\n",
      "Dimension = 100\n"
     ]
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "dim_range = [10*i for i in range(1, 11)]\n",
    "time_range_matmul = []\n",
    "time_range_numba_matmul = []\n",
    "time_range_np = []\n",
    "for n in dim_range:\n",
    "    print(\"Dimension = {}\".format(n))\n",
    "    a = np.random.randn(n, n)\n",
    "    b = np.random.randn(n, n)\n",
    "\n",
    "    t = %timeit -o -q matmul(a, b)\n",
    "    time_range_matmul.append(t.best)\n",
    "    t = %timeit -o -q numba_matmul(a, b)\n",
    "    time_range_numba_matmul.append(t.best)\n",
    "    t = %timeit -o -q np.dot(a, b)\n",
    "    time_range_np.append(t.best)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEPCAYAAABoekJnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzde3yT5f3/8dfVNGnTND0XWnqgLSc5yknACcpAxXkA9KsTJ1PwwJiH6XRjDnXoNkU2N3FT5xSwE0/j52FT5zwCOjdFUJByEDm1UI490XOb0/X7405Dj0Bp2qTN5/l49EFz586dqyHNu9dZaa0RQgghTldYoAsghBCie5MgEUII0SESJEIIITpEgkQIIUSHSJAIIYToEAkSIYQQHRIe6AIEQlJSks7Kygp0MYQQolv58ssvi7XWyc2Ph2SQZGVlsWHDhkAXQwghuhWlVEFrx7t9kCilbMBTgANYq7V+McBFEkKIkBKUfSRKqRVKqaNKqS3Njl+klNqhlNqllLrHe/gK4FWt9c3A9C4vrBBChLigDBIgF7io8QGllAl4EvgeMAS4Rik1BEgH9ntPc3dhGYUQQhCkQaK1/gQobXZ4HLBLa71Ha+0AXgFmAIUYYQJB+vMIIURP1p0+eNM4XvMAI0DSgNeB/1NK/QV4q60HK6XmKaU2KKU2FBUVdW5JhRAihHSnznbVyjGtta4G5p7swVrrZ5RSh4DLLBbLGL+XTgghQlR3qpEUAhmNbqcDB9tzAa31W1rrebGxsX4tmBBChLLuVCNZDwxQSmUDB4BZwA/acwGl1GXAZf379z+l8+vr6yktLaWyshK3W/rxhWhgMpmw2+0kJCQQERER6OKINjjdHg4dq2N/WQ37S2vYX1bDDydkkRIb6dfnCcogUUq9DEwGkpRShcAirfVypdRtwHuACVihtd7anutqrd8C3ho7duzNJzu3vr6effv2ER8fT1ZWFmazGaVaa10TIrRorXE6nVRUVLBv3z4yMzMlTALE49Ecraw/HhSltb7vC8tqOVRei6fR3oWmMMU5/ZP8HiQqlHZIbFQjuXnnzp0nPPfQoUOYzWaSkpK6pnBCdEPFxcU4nU5SU1MDXZQeSWtNWY3TV5toHhQHympxuD1NHtM7JoL0+Cgy4q1kJESRER9FeoKVjPgoUmMjCTedfo+GUupLrfXY5seDskbSWdpTI6msrETW4xLixGJiYsjPz5cg6YCqepe3NlHD/rJab0gYQbG/tIZqR9Nm9fgoM+nxUQxOtXPhkN6kJxwPjbQ4K5FmU5f/DCEVJO3pI3G73ZjN5s4vlBDdmNlslv7Dk3B7NPtKaygoqWZ/WS2F3tpEQ82irMbZ5Pwoi4mM+CgyEqxMyEn01iqMoEiPt2KPDL7PpZAKkvbUSADpExHiJOR3pKmiynp2HK7km8MVfHO4kh2HK9l5tJI65/HmJ4spjLR4K+nxVoYNT/WFhvFvFPFR3a8/NqSCRAgh/KHW4ebbI0ZQbD9cwQ5vaJRUO3znJEVHcEaKnWvH92VQip3sJBsZ8VH0skcQFta9guJkQipI2jv8VwgR2tweTUFJtTcwKtnhDY2C0hoaxilZzSYG9o5m6uBeDEqJ4YwUO4NS7CRFh85ItpAKkvY2bYmebfLkyeTn55Ofnx/oooggcLJmKaUgO9HG4NQYZo5K8wZGDJkJUZh6WA2jvUIqSETr1q5dy3e/+10Abr31Vp544okW5xw9epT09HScTifnnXcea9eubffz5ObmcuzYMe68886OFlmI09bQLNU4MFo2S1k4IyXG1yx1RoqdAb3sWC1dPyKqOwipIJGmrROLjIzkpZde4g9/+EOLCWYrV65Ea014+Om/ZXJzc8nPz5cgEV3C4fKwr7TGGxqtN0tFmsMY2Nse0s1S/hBSQSJNWyd2+eWX8/LLL/PPf/6T73//+03ue+6557j44ov56KOPAlQ6IVpyezQHj9Wyp7ia/OJq9jb6Kiyr8c3qVgqyEm2ckSLNUp0hpIJEnNjo0aPZtm0bzz33XJMg+eKLL9i6dSu//e1vWwTJ+++/z/Lly1m/fj2HDh0iIiKCcePGce+993Leeef5zsvKyqKgwNjuufHQxjVr1jB58mRff8XatWv56U9/ypo1a1BKMWPGDJ544gmioqJ45JFHePbZZzl06BBDhgzhz3/+M+ecc47vWrm5ucydO9d3zcakP6T70tpYBmRPUTX5JU3DYl9JTZOZ3VEWE9lJNoanxzJjZB+yEm0M6B0tzVKdTIJENDF37lzuuusuCgsLSU839gtbsWIFvXr14tJLL21xfm5uLqWlpVx33XWkp6dz4MABli1bxtSpU1mzZg2TJk0CYOnSpfzyl7+kuLiYxx57zPf4wYMH+76vrq5mypQpnHvuuTzyyCOsX7+eFStWUFdXR2JiIuvWreP222/H6XTy6KOPctlll1FQUIDdbu/kV0V0hbJqR9OaRUk1e73hUdNodrfFFEbfxChykmxMPaMX2Uk2spJs5CTZSLZHdLs5GD2BBEk7PfjWVrYdrAh0MZoY0ieGRZcN9cu1Zs+ezYIFC3j++edZuHAhtbW1vPLKK9x0002t9o88++yz2Gy2Jsfmz5/P0KFDWbx4sS9IZs6cydKlS6mtrWX27NmtPndxcTELFizg5z//ue86ZWVlrFq1itGjR/PZZ5/5VhsYPHgwM2bM4KWXXuJHP/qRX3520fmq6l1NmqDyi6uN8Cip5lijGd6mMEV6vJXsJBvjshPISbaRlWgjO8lGnzirNEcFmZAKEulsP7nExESmT59Obm4uCxcu5PXXX6e8vJwbbrih1fMbh0hVVRX19fWYTCbGjx/P559/3q7nNplM3H777U2OTZo0iTfeeIP58+c3WbKmIaBOtvim6Hpaa/aX1rL9cEXTsCiu5mhlfZNzU2MjyU6yccnwVLKTbL7aRUZ8FJbw7rRdUmgLqSDxR2e7v/7yD2Zz587lkksu4dNPP2XFihWMGzeOIUOGtHru7t27uffee3nvvfc4duxYk/va28SQmppKZGTT5a3j4+MByM7ObvV4SUlJu55D+J/Wmg35pXxZUMaXBWV8ta+M4qqmQ2mzEm2cNzDZ1wSVlWTUMKTfomcIqSARp2batGmkpaXx4IMPsmbNGv7yl7+0el5VVRXnnnsu1dXV3HnnnQwfPhy73U5YWBiLFy9m9erV7Xpek6ntD5W27mu8DcKJgsvlcrWrLKJtTpeHaoeLGoebGoebg+V13PT8ZwBkJUZx7sBkRmfGMywtluwkG7HW4FtkUPiXBIlowWQycd1117F48WKsViuzZs1q9byPPvqIgwcPsmLFCubOndvkvvvuu6/F+Z3dCZqQkABAaWlpi/v27t0rqzmfBo/W1DmNwKipN8KjYZRUmFJYzSaiI8J59rqxjMqMk/kXIUqCRLRq/vz5WCwWcnJyaGuP+4ZaQvPN0d5//33WrVvX4vzo6GjKysrQWndKqAwcOBCADz/8kCuuuMJ3/OWXX+bgwYP07dvX78/Z07jcHmocbl+No9bhxuP9/zWbwoiymEiyRBAVYSLSbCJMKRwlZiYM7h3gkotAkiARrcrMzOSBBx444TkTJ04kJSWFu+++m/z8fNLT09m0aRMrV65k+PDh5OXlNTl/woQJvP3229x222185zvfwWQyMWXKFHr16uWXMg8aNIjzzz+fv/71r2itGTlyJJs2beKNN96gf//+OJ3Ok18khGitqXd5qK4/3kxV7zKG2SoUVouJBJuFKIuJKEu4dH6LNoVUkMioLf+Ki4vjvffeY8GCBfz5z3/G5XIxZswY3nnnHZYvX94iSO6880727NnDq6++ytNPP43H42HNmjV+CxIwlnK5/fbbefHFF1m5ciWTJk1izZo1/PjHPw75yYhuj8cXGNX1Lmodbtze2kZ4mFHbiLeZsVnCsZpNPW6pc9F5QmrP9gZjx47VGzZsOOE527dvbzJZTojuRGuNw+Wh2uGmxttMVedsqG1AhNmEzVvTiIowYTGFnXZzo/yuhA7Zs12IHszj0dQ4vaFRb9Q6XB6jU9wUpoiyhBNrNXubqUyYwqSZSviPBIkQ3YzWGqevU9wYTVXn9KAxWhciwk3YI8OxRRg1jojw069tCHEqJEiECHIej6bWebyJqtrhxtVoCG6UxUSyPcJX2wg3SW1DdC0JEiGCjMPl8YVGjcNNrdPtG2JtCQ/DHhHuC41Is0lqGyLgJEiECCCP1tR6A6MhPJyNJ/xZTCRFW4yRVBYTZqltiCDU7YNEKZUD3AvEaq2vDHR5hDgRp9vjmyFe3by2YQrD5h1F1VDbCJPahugGAhokSqkVwKXAUa31sEbHLwIeB0zAMq31I21dQ2u9B7hRKfVqZ5dXiPZouryIUeNoWF5EeZcXSYo+PuFPahuiuwp0jSQXeAJ4vuGAUsoEPAlcABQC65VSb2KEyuJmj79Ba320a4oqxIk1jKRqGIJb62y5vEiiJQKbxUSkRWoboucIaJBorT9RSmU1OzwO2OWtaaCUegWYobVejFF7OS1KqXnAPDCW/xDCH+qcbipqnZTXOqltmPDnrW3I8iIiVAS6RtKaNGB/o9uFwPi2TlZKJQIPAaOUUr/0Bk4LWutngGfAmNnuv+KKUKK1MRS3vNZJRa3LtzZVlCWclJhIbBGyvIgIPcEYJK39Brb5wa+1LgHmn9KFZa0tcRq01lTXuyivc1FR68Tp9qBQ2CJMJEVbiYk0Y5YahwhhwRgkhUBGo9vpwEF/XNgfOySK0ODxaKrqXZTXOqmsc+LyaMKUIjrCqHnYI8Nl4p8QXsH4m7AeGKCUylZKWYBZwJv+uLBS6jKl1DPl5eX+uJzoASZPnkxWVhZgrI57rMZBQUk12w5VkF9STUWdE3ukmb6JUQxOjSEryUa8zRLUITJnzhyZpCi6VEB/G5RSLwOfAYOUUoVKqRu11i7gNuA9YDuwSmu91R/Pp7V+S2s9r62NmkLV2rVrUUqhlOK2225r9ZyjR49isVhQSjF58uTTfq7c3FyWLl162o/3N62NYbp7i6vZdqiSfaU1VDvcxEeZyU6yMTg1hoyEKGKtFvrlZKOUYuLEia1eq+EDvLi4uIt/CiECK6BBorW+RmudqrU2a63TtdbLvcff0VoP1Fr301o/5K/nkxrJiUVGRvLSSy9RX1/f4r6VK1eitSY8vGOtocEQJA6Xm6LKenYfraLG4cLt0dS73CRFW+iXHM3gFDtp8VHYI82tDtH973//yz//+c8AlFyI4BS89fNOIDWSE7v88sspKytr9UPyueee4+KLLyYiovvtyd0w0upIRR3fHqnkm8OVHCqvxa015vAwwk1hDOptJzXWii0i/ITNQn379iU5OZlf/vKXuN3uLvwphAheIRUkUiM5sdGjR3PmmWfy3HPPNTn+xRdfsHXrVubOndvq495//32uvvpqcnJysFqtxMXFceGFF/Lxxx83OS8rK4uPP/6YgoICX1OaUoq1a9cCx/sr8vPzufzyy4mLiyM+Pp45c+ZQVVWFx+Ph4YcfJjs7m8jISEaPHs1///vfJs+Rm5uLUoo1a9ZQXe/iUHktO45UsvNIJZdOO58pZw0jNdbKoBQ7A3vbjQ2d4JT7FKKjo7nvvvvYvn07ubm5Jz3/RP0VSinmzJnju52fn49SigceeIBVq1YxcuRIrFYr/fv39/2f7Nu3jyuvvJKEhATsdjuzZ8+msrKy1esXFRVx3XXXkZiYiM1mY+rUqWzcuLHFeU899RQXXnghaWlpWCwWUlNTmT17dsjvKClOXTCO2uo0Mmrr5ObOnctdd91FYWEh6enpAKxYsYJevXpx6aWtzwfNzc2ltLSU6667jvT0dA4cOMCyZcuYOnUqa9asYdKkSQAsXbqUX/7ylxQXF/PYY4/5Ht94d73q6mqmTJnCueeeyyOPPML69etZsWIFdXV1JCYmsm7dOm6//XacTiePPvool112GQUFBdjtdt+SJAD7SmtILKpCeUdaJUdHEGUJJzxMkWzvWK1q/vz5PP744yxatIgf/OAHWK3WDl2vubfffpunn36aW265hYSEBJYvX84NN9yAxWJh4cKFTJkyhYcfftj32kRGRrJs2bIW17noootISEjggQce4PDhwzzxxBOce+65fPbZZwwb5luRiEcffZQJEybwk5/8hISEBLZs2cKyZctYvXo1eXl5JCYm+vXnEz1PSAWJX/z7Hjicd/LzulLKcPhem8uRtcvs2bNZsGABzz//PAsXLqS2tpZXXnmFm266qc3+kWeffRabzdbk2Pz58xk6dCiLFy/2BcnMmTNZunQptbW1zJ49u9VrFRcXs2DBAn7+85/7rlNWVsaqVasYPXo0n332GWazGTACaMaMGSzPXcn/XTuHijonRZVG/05EeBiZCVFER4YT7t0N0F8DmSwWC7/5zW+49tprefzxx7nnnnv8c2Gv7du3s23bNvr27QvA1VdfTUZGBj/84Q959NFHueuuu4Djr83zzz/P0qVLiY6ObnKdvn378tprr/lqRFdccQVnnXUWP/vZz3j33Xd95+Xl5bX4/5s+fTrnn38+y5cvZ8GCBX79+UTPI01boonExESmT5/ua7Z5/fXXKS8v54YbbmjzMY0/hKqqqigpKcFkMjF+/HjWrVvXruc3mUzcfvvtTY5NmjQJrTXz58/HbDajtTHHI2f4GAA2btlORZ2TmEgzSdFGbSMl1kpclMUXIv52zTXXMHr0aJYsWUJpaalfrz1z5kxfiAAkJyczaNAgwsLCuPXWW5ucO2nSJJxOZ6vNUAsWLGjSrDZmzBguuOACPvzwQ6qqqnzHG/7/PB4P5eXlFBcXc+aZZxIbG9vu/z8RmkKqRuKXpi0//eUfzObOncsll1zCp59+yooVKxg3bhxDhgxp8/zdu3dz77338t5773Hs2LEm97V3PkNqaiqRkZFNjsXHxwOQlpHJ4fJajtU4cbg9hFmMv8DdtRUMSY0x1riymNr1fKdLKcUjjzzChRdeyEMPPcQf/vAHv107JyenxbH4+HhSU1NbDHZoeG1KSkpaPKZxk2GDIUOG8P7771NQUMDQoUMBWL16Nb/+9a9Zt24ddXV1Tc4vKys77Z9DhI6QChJxaqZNm0ZaWhoPPvgga9as4S9/+Uub51ZVVXHuuedSXV3NnXfeyfDhw7Hb7YSFhbF48WJWr17druc2mZoGgdPtobLOCcDBCgfplQ6iI8NJiY3EHmk0cYWHKV9gnSi4XC5Xu8pyMhdccAHnn38+Tz75JHfccUer57RVnhOVpflrcLLjgG9Pk5Npft769eu58MIL6d+/P4888gjZ2dlYrVaUUsyaNQuPx3NK1xWhLaSCRNbaOjUmk4nrrruOxYsXY7VamTVrVpvnfvTRRxw8eJAVK1a0GNV13333tTj/VGoobo+mos5JWbWD6noXx2qMIEmyWTgj1X7CfTsSEhIAWm1u2rt3r69/xV+WLFnC2LFjuf/++1v92RqXp+F7gD179vi1HK3Zvn07EyZMaHHMZDL5ms5eeukl3G43//73v8nOzvadV11dLbURccpCqo9E5pGcuvnz57No0SKefvppTvR6NfyV3Pwv3ffff7/V9vXo6GjKyspanK+1xu0xvrYfqmB/aQ0Ol4dkeyQpsUZTV2yU5aSbPw0cOBCADz/8sMnxl19+mYMH/bJkWxOjR49m1qxZvPDCC+TltRyE0VZ5/NkU1pbf/e53TV7nr776ig8//JCpU6f6Oubb+v97+OGHpTYiTllI1UjEqcvMzOSBBx446XkTJ04kJSWFu+++m/z8fNLT09m0aRMrV65k+PDhLT5cJ0yYwNtvv81tt93G2WefjQfFyPETCbfFU+fdCCouykx8lLGXh1KqXTsHDho0iPPPP5+//vWvaK0ZOXIkmzZt4o033qB///44nc72vhQn9dBDD/Haa6/x1VdftbjvmmuuYeHChcybN49vvvmGxMRE/v3vf3fJMioFBQVMmzaN6dOnc+jQIZ544gmsViu///3vfedcfvnlPPbYY1x88cXMmzcPi8XCBx98wObNm0lKSur0MoqeIaRqJML/4uLieO+99xg/fjx//vOfufvuu9m2bRvvvPMOo0ePbnH+nXfeyZw5c1n1/17l+uuv5/ofzmbDpjxsESYizCbMpjDS46NOOsP8RFauXMkVV1zBiy++6Au4NWvWkJaW1tEft1XZ2dnMn9/6TgYxMTG88847DB06lIcffphFixbRp0+fJsNvO8u7775L7969WbRoEY899hhjxozh448/ZsSIEb5zzjnnHF577TVsNhv3338/DzzwAFarlY8//rjFkGAh2qJOtZOuJ2jUR3Lzzp07T3ju9u3bWx31Ik6Py+2hvM7JsRon1fVGR7PNEk5clJlYqzmoV9MVJya/K6FDKfWl1nps8+Mh1bQlM9u7lkdrKutcHKtxUFHnQmtNRLiJlJhI4qLMWMK7ZqiuEKJzhVSQiM6ntabG4aasxkF5rRO3RxMeFkaizUJclBmr2SR7ZQjRw0iQCL+oc7o5VuPkWK0Dh8tDmFLEWM3ERZmxd6C/QwgR/CRIxGlzuj2U1zo5VuOgxuFGAbaIcHrbI4mxmjGFSXgIEQokSES7udweiqrqKa5yoLXGajaRGmslLsrcrqG6QoieIaSCRGa2d4zHoymprudoZT1ujyY+ykKyPYJIs3SaCxHKQurPR5nZfnq01pRVO/j2SCWHyuuIsoQzoJedjIQoCREhRGjVSET7aK2prHdxuLyOOqcbq9lETryN6Ej/rlclhOjeJEhEq2ocRoBU1buweDeJirWaZfSVEKIFCRLRRL3TzeGKOsprnYSHhdEnzkqCzUKYBIgQog0SJAIwhvIerayntMqBUtDLHkmy3YKpk3YYFEL0HBIkIc7t0RRX1VNUWY/WEG8z0zsmUobxCiFOWY8IEqXUTOASoBfwpNb6/QAXKeh5vCOxjlTU4/J4iLUaASKjsIQQ7RXwPzuVUiuUUkeVUluaHb9IKbVDKbVLKXXPia6htf6H1vpmYA5wdScWt9vTWlNe42DnkSoOHKslIjyMA9u+JCspGqslnGXLlrX6OKUUl156aReXVgjRHQQ8SIBc4KLGB5RSJuBJ4HvAEOAapdQQpdRwpdTbzb56NXrofd7HiVZU1bvYXVRNQWkNCshKtJGTbMNqOV4LWbRoEbW1tYErpBCi2wl4kGitPwGab7A9Dtiltd6jtXYArwAztNZ5WutLm30dVYYlwL+11i23qQtxdU43+cXV7Cmqwun2kB5vZUDvaGKaDecdO3YsBw8eZOnSpQEsrRCiuwl4kLQhDdjf6Hah91hbbgfOB65USrW6VZ1Sap5SaoNSakNRUZH/ShrEHC4P+0tr2Hmkkup6FymxkQzqbSfBFtHqfJDvf//7jBkzhiVLllBSUnLCayulmDNnTovjubm5KKVYu3at79gDDzyAUopt27Zx5513kpqais1mY+rUqezYsQOA119/ndGjR2O1WsnKyuKZZ55p8zk//PBDJkyYQFRUFCkpKdxxxx1UV1f7zvvjH/+IUqrFPukA9fX1JCQkMHXq1BP+fEKIUxesQdLapIU2t3LUWv9Jaz1Gaz1fa/10G+c8AzwIfGWxWPxUzODk8ng4VF7Lt0cqOVbrJDE6gkEpdnrZIwk7wYq8SimWLFlCeXk5Dz30kN/Ldf311/P111+zcOFC7r77bj7//HOmTZvGypUrufXWW5k5cya///3viY+P50c/+hGffvppi2t89dVXzJw5k7PPPptHH32USZMm8ac//Ynp06fj8Xh8zxMREcHy5ctbPP6NN96grKyMG2+80e8/nxChKlhHbRUCGY1upwMHO3pRf+yQuOSLJXxT+k1Hi+JXZyScwS/G/cK7qKKDo5V1vkUVe8dEtGsnwqlTp3LBBRfw1FNPcccdd9C3b1+/lTMlJYU333zTVxtKSkrijjvu4JZbbmHr1q1kZmYCcPXVV5ORkcGTTz7JxIkTm1wjLy+PN954g5kzZwJwyy23cMcdd/CnP/2JVatWMWvWLBITE7niiit4/fXXKS0tJSEhwff45cuXEx8fzxVXXOG3n0uIUBesNZL1wAClVLZSygLMAt7s6EWVUpcppZ4pLy/vcAGDiabxooq13kUVo8lIiDqt7WyXLFmCw+Hg/vvv92s5f/KTnzRpUps0aRIAM2bM8IUIQHJyMoMGDWLnzp0trjFo0CBfiDS45x5jUN8bb7zhOzZv3jzq6+t58cUXfcfy8/P56KOPuPbaa4mMjPTPDyWECHyNRCn1MjAZSFJKFQKLtNbLlVK3Ae8BJmCF1nprR5/LHzWSX4z7RUeL4Tdaa6rqXRwqr2N/WY3fFlUcNWoU11xzDS+++CI/+9nPGDFihF/Km5OT0+R2fHw8ANnZ2S3OjY+Pp6CgoMXxwYMHtziWmppKXFwce/bs8R2bPHkyAwcOZPny5dx+++0APPfcc2ituemmmzr0cwghmgp4jURrfY3WOlVrbdZap2utl3uPv6O1Hqi17qe19kuDfU+qkdQ4XOwtrmZvcTUerclMiKJ/r2i/rcz729/+lvDwcH7xi/YFp8vlavM+k6n12lFbx7Vu2S3W1qKRrZ1788038/XXX/Pll1/i8XjIzc1l7NixnHnmmW2WUQjRfgEPkq7UE/Yj8WhNYVkNu45WUef00CfOysDeduKiLH5dmTc7O5sf//jHvPvuu6xZs6bF/QkJCZSWNh+1TZNaQWfYtm1bi2OHDh2ivLy8RY1nzpw5WCwWli9fzgcffMC+ffukk12IThBSQdLdeTyafSU1lFY7SI6OYFBKNEnREZ22Mu99991HTExMq7WSgQMH8tlnn1FTU+M7VlZWxnPPPdcpZWmwY8cO/vGPfzQ5tmTJEoAWfSdJSUnMnDmTl156iSeeeIKoqCh+8IMfdGr5hAhFIRUk3blpy+3R5JdUU1HnpE+cldQ4a6evzJuUlMTPf/5z1q9f3+K+2267jeLiYqZMmcKTTz7JQw89xKhRo/w6yqs1w4cPZ/bs2fz0pz/lqaee4qqrruLxxx/nvPPO4+qrW66OM2/ePMrLy3n77be56qqriImJ6dTyCRGKQipIumvTltvjIZ36oq0AACAASURBVL+4mup6F+nxUSRFR3TZc991112kpqa2OH7ttdfyu9/9jkOHDnHXXXfxwgsv8Ktf/Yr581udD+o3o0eP5h//+Af/+9//uPvuu/nkk0+47bbbeOuttwhrJVinTJlC//79AaRZS4hOolrrpOyplFKXAZf179//5taGlja2ffv2VkcIdTWX28Pe4mrqXB4y463ERvXsyZQnopTi+uuvJzc3t12PGzp0KG63m2++Ca75Pz1FsPyuiM6nlPpSaz22+XGpkQQxp9vDnqJq6l0e+iZEhXSInK7Vq1ezbds25s2bF+iiCNFjBXweiWidw+VmT3E1LrcmK9FGdKT8V7XH6tWr2b17N4sXLyY5OZmbbz7tqUNCiJOQT6cgVOd0++aH5CTZiIqQ/6b2+vWvf82nn37KkCFD+Nvf/obdbg90kYTosdr1CaWUsgM/BS4EegPXaa0/U0olAbcAq7TWQdsQ3aiPJNBFaVOtwwgRgJyk6CZ7hYS69vTnNV59WAjRuU65j0QplQxsAO4HEoEcwAqgtS4GrgeCuiE62PtIqutd7CmuQilabDglhBDBqj01kt8CKcB4YB9wtNn9/wRkk4fTVFXnJL+khnCTIifJdlqLLQohRCC0Z9TWpcBT3h0IW2tj2EPTpd+7va4aGl1Ra4SIxRRGv+RoCRHRbYTS9AHRtvYESRKw6wT3e4CgXpu7PTPbTSYTTqez08t0rMZBQWkNEeFh5CTbMJtCakS26OacTmebi26K0NGeT63DQL8T3D8Ko8kraLWnj8Rut1NRUdGp5SmtdrC/tIYos4mcZBvhEiKim6moqJARcaJdQfIOcKNSqsV6GUqp8cB1GP0kPUJCQgJlZWUUFxfjcDj8XoUvrqqnsKwGW0Q4WUm2Tl83Swh/0VrjcDgoLi6mrKysyQ6UIjSd8hIpSqkU4EuMjabeBG4EXgAswBUYW+GO0Vq3XFs8yIwdO1Zv2LDhpOfV19dTWlpKZWUlbrfbb89fWeekvNaF1RxGgs2/y78L0RVMJhN2u52EhAQiIrpu7TcRWG0tkXLKo7a01oeVUhOAJ4AbAAX8EKPj/R3gx90hRNojIiKC1NTUVhctPB1aa37/3g6eWrufGSP78OhVI6RPRAjR7bVrQqLWej8wQykVAwzCCJNdPS1AOoPHo/n129vI/V8+14zL5Lczh2EKk5qIEKL7O621N7TWFUDLTSqCXKBmtrs9ml+8tplXvyzk5knZLLx4sDRnCSF6jHYHiVIqCsjCmN3e4tNQa/1Jx4vVObTWbwFvjR07tstW8HO4PPz075v4V94h7jx/AHdMHSAhIoToUU45SJRSNuAxjNFZ5tZOwegvkUHlXnVONz9+4UvW7CjivksGc9OknJM/SAghupn21EieBq4F3gD+A5R1Sol6iKp6Fzf9bT3r9pby8OXD+cH4zEAXSQghOkV7gmQGsFxrLRs7nER5jZPrn/uCvAPlPPb9kcwclRboIgkhRKdpT5A46YYd7F2tqLKeHy5fx56iap66djTThqYEukhCCNGp2jOJYTXGyr+iDQeP1XL1Xz+joKSG5XPGSogIIUJCe4LkbmCqUuoOpVRrne0BoZQarJR6Win1qlLqx4EqR0FJNVc9/RlFlfU8f+M4Jg1IDlRRhBCiS51ykGit9wELgT8A1UqpAqXUnmZfu9vz5EqpFUqpo0qpLc2OX6SU2qGU2qWUuuck5dqutZ4PfB9oMXW/K+w8UslVT39GjcPFSzdP4KwsWXtICBE62jP8dw6wHHAAO/DPqK1cjCVXnm/0PCbgSeACoBBYr5R6E2NY8eJmj79Ba31UKTUduMd7rS615UA5P1y+jnBTGH//0dkM7C0roQohQkt7OtvvBTYB07xb63aY1voTpVRWs8PjMJZd2QOglHoFmKG1XoyxuVZr13kTeFMp9S/gJX+U7VRsyC9l7nPribGaefGm8WQl2brqqYUQ4sTcTigrgNLdULILSnYb31/8KCQN8OtTtSdI0oDH/BUiJ3me/Y1uF3KCTn6l1GSM1YcjMBaPbOu8eXj3lM/M7Picjv/sLGLe81+SGhvJCzeNp0+ctcPXFEKIdvG4obzQGxa7j4dFyS4jRHSjVcsjYyGxP9RX+r0Y7QmSHUBXNP63tn5Im2vda63XAmtPdlGt9TNKqUPAZRaLZcxplw54f+thbntpIznJNlbeOJ5kuyyjLYToJFpD5eFmNYs9xvele8Fdf/xccxQk9oOUETD0CuP7hH5GgEQlQCctz9SeIHkYeFwplau1LuyU0hgKabr3ezrGXicd5o+1tv656QB3rfqaYWmx/G3uWcRFWfxRNCFEKNMaakqbhkXJLu/tPeCsPn6uyQIJOUZADLiwaVjYUzotLE6kPUEyGDgAbFdKvQHsBZrv9qS11r/pYJnWAwOUUtne55sF/KCD1wQ6vvrvy1/sY+EbeYzPTmDZ9WcRHXFaiycLIUJVXUXTZihfWOyGumPHz1MmiO9rBETfid6wyDHCIjYdwoJrScP27JDoOYXTtNb6lH9CpdTLwGQgCTgCLNJaL1dKXQwsxRiptUJr/dCpXvNUnOoOic09sXonGwrKeHr2GCLNwfUfKYQIIs46KNoOh/Pg8BY4sgWKd0L10abnxWYcD4jGNYv4vmAKmul6Pm3tkNieIOl7KudprQvaWbYu06hGcvPOnTvb/XitNW6PJlx2NRRCNKgugcObjdA4ssX4t2jH8Y5usw16D4XkQU3DIiEbzN1rkE6Hg6QnOd0aiRAihHk8ULbXW8to9FXZqAvX3gdShjf9is+GsJ7xx2eH92zvCQK1Q6IQoptx1MDR7XCkUWAc2QqOKuN+ZTJqGNmTjgdG7+FgSwxsuQOkzRqJUupXGMNuH9Jae7y3T8Yfne2dTmokQgifqqMtaxklO0F7u4Ut9ma1jGGQPBjMkYEtdwC0u2nL27muAavW2tEZne2BIkEiRAjyuI35Fw39GQ1fVUeOnxOb0aiGMcz4N65vj2ma6qjTadrKBtBaOxrf7s6kaUuIEOGoMZqiGofG0W3grDHuDws3ahX9pho1jIbgiJIFV0/HCTvblVKZQJHWurbritT5pEYiRA9SXQyHvvYGhjc4SnYdb5qKjDVmejfUMFKGG/0b4bIiRXudbmf7XuCHdOFCiEII0aoWo6a8oVF56Pg5DU1TQ684HhpxmQGZ7R1KThYkPerVl6YtIboJV70xaqpxYBzeAg7vgoO+UVPnGrWNhtCQpqmACKnhv/5Ya0sI4We1x1rWMoq+AY/LuN9sM/oxzpx1PDB6DQnJUVPBKqSCRAgRQFobS543CY3NcGzf8XOiextBMeBC49/UM3vUhL6e6lSCZJJS6pQDR2v9/MnPEkL0aG4XFH/bNDAO50Ftw8aqylguJG0sjJl7vHnK3jugxRan51QCwrch1EkojHknEiRChJq6cihcD/vWwf7PofDL40ufmyKg9xAYPP14LaPXEIiIDmyZhd+cSpA8A3ze2QXpCtLZLoQfaA1l+bB/nfG1b50xRwMNKswYZjvqWqO2kToCEgeASVrRe7JT+d/9j9a6Rwz/lc52IU6Dy2E0Te373Kht7P/i+Gxwix0yzoIhMyBzPKSNgQh7YMsrupz8mSCEaKqmtGlt4+BX4Koz7ovrC9nnGaGRMQF6DQ66TZZE15MgESKUaW1suLR/3fHaRvG3xn1h4UZ/xtgbIWMcZIyHmNTAllcEJQkSIUKJsxYObvQ2U31hBEhtqXFfZJwRFmfOMmobfUaBJSqw5RXdwgmDRGstg7eF6M4qjzRqpvrcWJPK4zTuS+wPg75nhEfmBKNTXOZriNMQUjUSGbUleryS3bD34+PDcMvyjeOmCKOGcfYtRnBkjAdbUkCLKnoO2WpXiO7MUQP5n8KuD2DnB8aihgBRSUYto6G2kXqmrHYrOky22hWipyjZbYTGrg+MEHHVQbjV2PZ1wi3Qb4oxa1xWvBVdRIJEiGDXVq0joR+MmQMDLoC+54DZGtBiitAlQSJEMDpZrWPA+ZCQE+hSCgFIkAgRHKTWIboxCRIhAkVqHaKH6BFBopSyAZ8Ai7TWbwe6PEK0qq1aR2J/Yyn1AedLrUN0SwENEqXUCuBS4KjWelij4xcBjwMmYJnW+pGTXOoXwKpOK6gQp0tqHSIEBLpGkgs8QaM9TJRSJuBJ4AKgEFivlHoTI1QWN3v8DcAIYBsg+26KwJNahwhBAQ0SrfUnSqmsZofHAbu01nsAlFKvADO01osxai9NKKW+C9iAIUCtUuodrbWnlfN8G3RlZmb688cQoa6sAL59D3a+17LWcfat0P98SMgOdCmF6DSBrpG0Jg3Y3+h2ITC+rZO11vcCKKXmAMWthYj3vGcwNuli7NixoTedX/iPx23sBvjtu0aAHN1mHE/o16jWMRHMUkkWoSEYg6S16bgn/eDXWuee9MKy1pY4XXXlsOsjb83jfWPF3LBwyDwbLnwIBl4ESfK+EqEpGIOkEMhodDsdOOiPC8sOiaJdSnZ7ax3vQsH/wOMCazwMuBAGToN+U8EaF+hSChFwwRgk64EBSqls4AAwC/iBPy4sNRJxQm6nsdR6Q5NVyU7jePJgOPs2Y8n19LNkR0Ahmgn08N+XgclAklKqEGMeyHKl1G3AexgjtVZorbf64/mkRiJaqCmFXR8a4bHrQ6MJy2SBrIkwbh4MvBDiswJdSiGCWqBHbV3TxvF3gHf8/XxSIxFoDUU7jtc69n8O2gO2ZDjjMhh0EeRMhgh7oEsqRLch+5GIns9VDwX/NYLj23ePb/aUMhwGfs/oKO8zSnYHFOIkZD8SpEYSUqqKjNFV3/4bdq8BRxWER0L2eXDOHTBgGsSmBbqUQvQIUiMRPYPWcGSLUePY8S4c+BLQYO9jjLAaeBFknwuWqECXVIhuS2okoufxeIw+ji2vGeFRUWgcTxsD311oBEjKCNkpUIhOFlJBIk1bPUTRDti8CvJWwbF9YI4ytpedfI8xx8PeO9AlFCKkSNOW6B4qDxs1j82r4NAmUGGQ810YcTWccQlERAe6hEL0eNK0Jbqf+ir45m3Y/HfYs9YYpps6EqYthmH/JzUPIYJESAWJNG11A26XERqb/26EiLMGYjNh4l0w4vuQPCjQJRRCNBNSQSIz24OU1nBwo9FsteVVqC6CyDij2WrE1ZAxXuZ4CBHEQipIRJApy4fN/8+ofZTsNJYmGXiRER4DLoDwiECXUAhxCiRIRNeqKYWtbxi1j/2fG8f6ToTv3A5Dphur6wohupWQChLpIwkQZ50xUXDzKmO2uccJyWfA1EUw/EqIkx0rhejOZPiv6Bwej7G+1ea/w7Y3ob4colOM4BhxtbHOlUwUFKJbkeG/omsc2WaER96rxkxzSzQMvswYcZV9nuzlIUQPJEEiOq7ioBEcm1fBkTxQJug/FS540NgMymILdAmFCBke7aHKWUWlo9L3VeGo8H0/LWsavaJ6+fU5JUjE6fG4YcvrsOkF2PMxoI01rr73Oxh6BUQnB7qEQnRLWmtqXDVUOiopry8/HgjOpqFQUV/R5HjDfVWOKjRtd1kMih8kQSICTGvY9RF88Cs4utXYPfC8BTD8+5AkgxiEaMyjPRyuPsze8r2U1pU2qSX4gqG+4niNwXvMoz0nvK7NbMNusRtfZjspUSkMjB/Y5JjdYifGEuM7FhNhfB9t9v9yQiEVJDJqq4MOfQ3v3w97PzYC5MrnYOjl0mkuQp7b4+ZA1QF2H9vN7vLd7Dm2h93lu9lbvpdaV22L863h1iYf9MlRyeTE5fhuNw4AXxBYYoixxGAz2wgPC66Pbhm1JU7u2H5Y/VujE90aB+f9AsbeCOGWQJdMiC7l9DjZX7Gf3eW72X3seGDkl+fj8Dh85/WO6k2/uH7kxObQL64f2bHZJFmTiLHEEG2JxhxmDuBPcfpk1JZov9pj8Okf4fOnjdvn3AETf2qEiRA9WL27nvzyfPaU7zECw/vvvop9uLTLd15adBr94vpxTp9zyInLISfW+Iq2hNZq1BIkoiWXAzYsh49/B7VlxryPKfdBXEagSyaEX9U4a9hbsdeoWXibpfaW72V/5X5fP0WYCiPTnkl2bDZTMqf4ahlZMVlEmWXHTZAgEY1pbSxf8tGDxjpY2efBhb+B1DMDXTIhOqTSUcme8j0tAuNA1QHfOeEqnL4xfRkYP5DvZX+PfrH9yInLoW9MXyJMsu7biUiQCEPBZ/D+fXBgA/QaCte+ZswFkY500Q6Vjkq2FG9hc9FmimqL0Fqj0b6/7j3a0/ZtDR48vsdorX33N75O49taazw0fWzja3swRk0drTnqK6MlzEJ2bDYjkkdwef/LfX0ZGTEZ3bbvItAkSEJd8U74YBHs+BfYU2HGk3DmNTIDXZyUy+Ni97HdbC7ezOaizeQV5bGnfI9vDkNcRBxhKgyFQilFGGGg8B0LU8bWAG3dVsr75b3Pd6yN2w2PNYWZfM+lUIxPGU9OXA79YvvRL64fadFpmOT97VfdPkiUUpOB3wBbgVe01msDWqDuouoorH0Evsw19jyfch9MuBUs0uYrWldUU8Tmos2+4NhastU3tDUuIo7hScO5KPsiRiSNYFjyMGIsMQEusegqAQ0SpdQK4FLgqNZ6WKPjFwGPAyZgmdb6kRNcRgNVQCRQ2InF7Rkc1fDZk/Dfx8FVB2NvMIbzykx00Uidq47tpduN4PCGx+HqwwCEh4VzRvwZzOw/kxHJIxiRNIIMewZKmkFDVqBrJLnAE8DzDQeUUibgSeACjGBYr5R6EyNUFjd7/A3Af7TWHyulegN/BK7tgnJ3Px43bHwB1jwMVYfhjEvh/AdlNrpAa01BRQF5xXl8XfQ1ecV5fFv6rW+Yax9bH0Ymj2T44OGMSB7B4MTB0vksmghokGitP1FKZTU7PA7YpbXeA6CUegWYobVejFF7aUsZIO/u5rSGnR8YS5oUbYf0s+CqXOh7dqBLJgKkvL6cvOI88ory+Lr4a7YUb6G8vhyAqPAohiUNY86wOQxPMoIjyZoU4BKLYBfoGklr0oD9jW4XAuPbOlkpdQUwDYjDqN20dd48YB5AZmaIbKR0cBN8cD/s/QTis+Gqv8GQGTISK4Q4PU52lu00OsOL89hctJn8inzA6IjuF9ePqZlTGZE0guHJw+kX2086okW7BWOQtPYp1+Y6Llrr14HXT3ZRrfUzSqlDwGUWi2VMB8oX/I7tg49+A3mrwJpgrMg7Zq4sadIDaa2pd9f7lg2vclRxsPogeUV55BXnsa1kG3XuOgASIhMYkTyC6f2mMyJ5BEMTh4bcDGzROYIxSAqBxlOo04GD/riw1vot4K2xY8fe7I/rBZ3aMvjPH2HdX41ax8SfGl+RsYEumWiF1ppaV60RAM4q48tRRaXTCIQm3zcEhfecxue6PK4W1zaHmRmcOJgrB15pdIgnj6CPrY90iItOEYxBsh4YoJTKBg4As4Af+OPCPXb1X1c9rF8Gn/zeWB/rzGtgyr0Qmx7okoUEp9vJ/qr9HKw6SEV9RZMP/UpHJdXO6lZDodpZjVu7T3hthSLaHE20xfiym42VYrPMWdjNduOYd2nwhvuTrEkMiB+AxSQ1UNE1Aj3892VgMpCklCoEFmmtlyulbgPewxiptUJrvdUfz9fjaiRaw9bX4cMH4VgB5HwXLvg1pI4IdMl6HI/2cKT6CPkV+RRUFFBQUeD7/kDVgVb3jzApkxEA5uMf9qm21KbHmt1vt9h9e01Em6OJMkf5JtsJEaxCahn5RjWSm3fu3Bno4nRM/n+NJU0OfgW9hxkB0n9qoEvVrWmtKasvM0KivGlg7K/cT7273neuNdxKVkwWfWP6+r7S7enEWmJ94WANt0pTkuhR2lpGPqSCpEG33o+krADeWwjfvA32PsaM9DNnyZIm7VDtrG5RqygoN25XOit954WHhZNhz6BvTN8WoZFsTZaQECFH9iPp7lz18L8/wSd/MDrSZUmTE3K4HRRWFrZoitpXsY+i2iLfeQpFqi2VvjF9uSTnErJis8i0Z5IVk0VqdGrQ7UQnRDAKqd+SbtvZvns1vPNzKNkFg6fDtIdlbxCvKkcVecV57C3fawRGpVG7OFh9sEm/RUJkAlkxWUxMm0hmTKavhpFhzyAyPDKAP4EQ3Z80bQWz8gNGM9a2f0BCDnzv9zDg/ECXKqAOVx9m49GNfHXkKzYVbeLbsm99gWEz24ymJ3tf+sb29TVJZcZkygKCQviBNG3RjWokbid8/hdjdV7thu/eC9/5CZhD6y9nj/aw+9huIziOfsXGIxs5WG1MKYoKj2JE8gjmj5jPmb3OZGD8QBIjE6XfQogAkBpJsMn/FP71M2NdrAHT4HtLICE70KXqEvXuerYWbzVC4+hGNh7dSKXD6PxOsiYxutdoRvcezaheoxgYP1D6L4ToYlIjCXaVR4x1sTb/HWIzYdZLMOjiHr0u1rG6Y2wq2uSrbWwt2YrT4wQgJzaHC/te6AuO9Oh0qW0IEaQkSALN7YINy2H1b8FZC5Puhkk/63GjsbTWHKg60KSZanf5bsAYZjs0cSizB89mVK9RjOw1kvjI+ACXWAhxqkIqSIKuj2T/F/Cvu+BwnjEr/eJHe8z+IC6Pi2/LvvU1UW08spGjtca+2XaznZG9RnJpv0sZ1WsUQxOHysgpIbqxkAqSoFkipboEPvyVsdGUvY+xP8iQmd26GavGWUNecZ6vtvF10dfUuGoASLWlMjZlLKN7jWZU71H0j+svy34I0YOEVJAEnMcDX+Uaa2M5quA7txvb3EbYA12ydtFac6j6EFuKt7CpaBMbj2xke+l23NqNQjEgfgCX9bvMCI5eo0iNTg10kYUQnUiCpKsc3Ahv32WsjdX3HLjkD9BrcKBLdUpKakvYWrKVLcVb2FK8ha0lWymtKwUgwhTB8KTh3DDsBkb3Hs2I5BEyZ0OIEBNSQRKQPpLaMmOTqQ0rwJYMlz8DI74ftM1YVY4qtpVsY0uJNzSKt/rmbigUObE5TEybyLCkYQxLHMYZCWdgNpkDXGohRCDJPJLO4vHA1y8be6XXlsJZN8N3F4I1rnOftx3q3fV8U/qNLzC2lGwhvzwf7d2QMi06zRcYQ5OGMiRxCDazLcClFkIEiswj6UqHt8C/7ob9n0P6WXDJ65B6ZkCL5PK42H1st9E8VWIEx86ynbi0sbtekjWJYYnDuDj7YoYlDWNo4lAZgiuEOCUSJP5UVwFrFxtb3UbGwvQ/w8jZENa1I5S01uyr3EdecZ5R0yjewjel3/j27rZb7AxNHMqcYXN8tY3eUb1lwp8Q4rRIkPiD1rDlNXjvXqg6AmPmwNRfQVRCFzy15kjNEV/TVENneMPSIpGmSM5IOIMrB17J0KShDE8aToY9Q4bfCiH8RoKko4p2GM1Y+f+B1JHG0ibpYzrlqWpdtRyoPMD+yv3sKNvhC4/i2mIAwlU4A+IHMC1rGsMShzEsaRj94vrJmlRCiE4VUp8wfh21VV8Fn/wOPnsSLDZjOO+YuR3aqdCjPRTXFlNYWUhhVSGFlYXsr9zvu90QGGCMoMqKzeLs1LMZmjSUYUnDGBQ/SGaICyG6nIzaai+tYfub8O5CqCiEkdfC+Q9CdPIpPbzOVceBqgPHQ8IbGA1h0XhfcIUixZZCuj2d9Oh00u3pZNgzSI9OJys2C7ule01kFEJ0bzJqyx9Kdhs7Fe7+CHoPg/9bBn3PbnKK1pqSupIWtYmG2423eQVjX410uxEME9MmGqHhDY4+0X2wmCxd+RMKIUS7SZC0x0e/hv1fUH/Bbzgw+CIKaw6zf/uLTWoUhZWFvtFRYNQqett6kx6dzjlp55Ae7a1VeAMjPiJeRksJIbo1CZJ2eDxjAG/qfRzd9SzsetZ33BpuJd2eTqY9k+/0+U6Tpqi06DSpVQghejQJknboHd+fs911x/sqvIGREJkgtQohRMiSIGmHWWfMCnQRhBAi6HT7IFFKhQG/AWKADVrrvwW4SEIIEVICOr1ZKbVCKXVUKbWl2fGLlFI7lFK7lFL3nOQyM4A0wAkUdlZZhRBCtC7QNZJc4Ang+YYDSikT8CRwAUYwrFdKvQmYgMXNHn8DMAj4TGv9V6XUq8BHXVBuIYQQXgENEq31J0qprGaHxwG7tNZ7AJRSrwAztNaLgUubX0MpVQg4vDfdnVdaIYQQrQnGlfvSgP2Nbhd6j7XldWCaUurPwCdtnaSUmqeU2qCU2lBUVNTWaUIIIdop0E1brWltHG2b67horWuAG092Ua31M0qpQ8BlFoulc1ZVFEKIEBSMNZJCIKPR7XTgoD8urLV+S2s9LzY21h+XE0IIQXAGyXpggFIqWyllAWYBb/rjwkqpy5RSz5SXl/vjckIIIQjw6r9KqZeByUAScARYpLVerpS6GFiKMVJrhdb6IT8/bxFQ4M9rBkASUHzSs0KHvB7HyWvRlLwex3X0teirtW6x1HlILiPfEyilNrS2nHOoktfjOHktmpLX47jOei2CsWlLCCFENyJBIoQQokMkSLqvZwJdgCAjr8dx8lo0Ja/HcZ3yWkgfiRBCiA6RGokQQogOkSAJckqpDKXUGqXUdqXUVqXUHd7jCUqpD5RSO73/xge6rF1JKWVSSm1USr3tvZ2tlFrnfT3+7p2DFBKUUnFKqVeVUt943ydnh+r7Qyn1U+/vyRal1MtKqchQem+0tqJ6W+8FZfiTd5X1zUqp0af7vBIkwc8F3K21HgxMAG5VSg0B7gE+0loPwFjx+GTL7fc0dwDbG91eAjzmfT3KOIVlc3qQx4F3tdZnAGdivC4h9/5QSqUBPwHGaq2HYcxDm0VovTdygYuaHWvrvfA9YID3ax7wl9N9UgmSIKe1PqS1/sr7fSXGh0Qaxj4sDZt4/Q2YGZgSdj2lVDpwCbDMe1sBU4BXvaeEzOuhlIoBzgWWA2itHVrrY4Tu+yMcsCqlwoEo4BAhzFPleAAABvdJREFU9N7QWn8ClDY73NZ7YQbwvDZ8DsQppVJP53klSLoR75L7o4B1QG+t9SEwwgboFbiSdbmlwALA472dCBzTWru8t0+2YnRPkgMUAc95m/qWKaVshOD7Q2t9AHgU2IcRIOXAl4Tue6NBW++F9q603iYJkm5CKRUNvAbcqbWuCHR5AkUpdSlwVGv9ZePDrZwaKsMRw4HRwF+01qOAakKgGas13rb/GUA20AewYTTfNBcq742T8dvvjQRJN6CUMmOEyIta69e9h480VEO9/x4NVPm62DnAdKVUPvAKRrPFUoxqecO2CH5bMbobKAQKtdbrvLdfxQiWUHx/nA/s1VoXaa2dGHsVfYfQfW80aOu94LeV1iVIgpy3/X85sF1r/cdGd70JXO/9/nrgn11dtkDQWv9Sa52utc7C6EhdrbW+FlgDXOk9LZRej8PAfqXUIO+hqcA2QvP9sQ+YoJSK8v7eNLwWIfneaKSt98KbwHXe0VsTgPKGJrD2kgmJQU4pNRH4D5DH8T6BhRj9JKuATIxfoKu01s072Xo0pdRk4Gda60uVUjkYNZQEYCMwW2tdH8jydRWl1EiMgQcWYA8wF+OPxJB7fyilHgSuxhjtuBG4CaPdPyTeG62tqA78g1beC96wfQJjlFcNMFdrveG0nleCRAghREdI05YQQogOkSARQgjRIRIkQgghOkSCRAghRIdIkAghhOgQCRIhvJRSc5RS2jusuEdRSuUrpdYGuhyiZ5IgET2SUmqyNxQavtxKqTLv8uJ/U0pd5B1HL4ToIJlHInokb61iDfAy8A7GukJ2YBDG6qeZwIcYk7OOeR9jAsyAQ2vtaeWy3ZZSKgLQWmtHoMsiep7wk58iRLf2ldb6hcYHlFJ3Ab8D7sIImv/f3rmFxlVFYfj7UUrFNkRRKvrgBSnRig/Fgg9CxVovFbT1QVCKFApKhYKiCcXirQoWvKCV6JNVtCgo0kuoGCsSbUOrgqhVMGhiBWlE1DS0qQqmy4e9R47HmXSSManM/B9s9jlr1tl7nQnMyj5rnb2uB4iIcWB8xi2cAZr1TW7z/8CPtkzLERHjEXEvsAe4Lm9DUzVGUpAtkfSgpO8l/ZYr7l2edRZL2iNpTNKwpAeqzSvpMklbJf0s6Q9JA5LWFzYUrOj15ZjG2bnK30geu1fS/JLubEkP57GOSjokab+kJ0p6VWMkkpZL6pd0JLd+STdV0TuQ7eqQtFPSYUmjSpUZz6r3uzfNiR2JaWVezP0NdehuJD0SexZ4hFQHpFfSctIus7uB+4CvgQ2SVhYvlrQM6AfmA0+RKvntBTaQVkVlTgU+JK2Q7ge6SXsobc+P4Cp0k/ZT2kdaYa0nVcG76ng3JOkuYCtpD6rHgEfz8TZJd1S55Bygj7RfUyfwGnAz8Mrx5jJNTkS4uTVdI/3oBmlTx1o6C7POW/l8VT6/sqBTkX0KzCrIb8zyP4FFBfksUlGlvQXZbOBHkmM4uWTDPVXm7MuyrpJuZ5ZfW5D9Crxdx/dxAOgrnJ8GHAG+BdoK8jZgEDgMtJeuD+CW0rjdWd5xov/mbieueUViWplKgbC2OnRfiH8Gqnfnfl9EfFIRZp2PSXWwKywF5gEvkWpjnFFppEQAgGtK8x0DNpVk7+e+OPYosEDSJXXcQ5GlpFXPpigUSsvHzwFzSPU9ihyMiDdq2HThJOc3TYSD7aaVqTiQeipODhVPImIkZw9/V0V3hFT+t8JFud88wfjzSucHI+L3kuyX3BfHvht4FdgvaYiUqdYD9MTEmWfn5/6rKp99mfsLSvKhsmINm0yLYUdiWplLcz9Qh26tbK56srwq76t0Ap/V0ClXppto3L/ff4mI7ZLOA5YBi0mriNXAbklXR+1036m8Q1OXTab1sCMxrczq3O+c5nm+yf1YRLz3Xw8eqWDVFmBLfslyI9BFql/+Zo3LBnO/gBScL3Jx7qutQIz5F46RmJZD0kmSngSuIAWq+6d5yl5Snex1kk6vYs8pkuZOdtB8H+1FWUQEqQogpAysWuwCxoC1xbnz8VpSIH7XZG0yrYlXJKbZWVhIxS2+2X4u8C5w23QbEBFjkm4nlTwdkLSZlC3VDnSQUmhXkLK1JsNcYFjSDpLz+IkU+1hDitP0TGDTIUldpKyrjyS9nD9aRQqc3xkRo5O0x7QodiSm2bk1t2Ok/7J/AD4AXo+Id2bKiIjolbQIWAesBM4k/dgPAk8DX0xh2KPAM8ASUmxkDin1eAfweESU4y5lm56XNEyK3TyUxZ8DKyJi2xTsMS2K99oyxhjTEI6RGGOMaQg7EmOMMQ1hR2KMMaYh7EiMMcY0hB2JMcaYhrAjMcYY0xB2JMYYYxrCjsQYY0xD2JEYY4xpCDsSY4wxDfEXGikhchKz5PwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(dim_range, time_range_matmul, label=\"Matmul\")\n",
    "plt.plot(dim_range, time_range_numba_matmul, label=\"Matmul Numba\")\n",
    "plt.plot(dim_range, time_range_np, label=\"Numpy\")\n",
    "plt.legend(fontsize=18)\n",
    "plt.xlabel(\"Dimension\", fontsize=18)\n",
    "plt.ylabel(\"Time\", fontsize=18)\n",
    "plt.yscale(\"log\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Why is na&iuml;ve implementation slow?\n",
    "It is slow due to two issues:\n",
    "\n",
    "- It does not use the benefits of fast memory (cache) and in general memory architecture\n",
    "- It does not use available parallelization ability (especially important for GPU) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Memory architecture\n",
    "<img width=80% src=\"Memory-Hierarchy.jpg\">\n",
    "\n",
    "- Fast memory is small\n",
    "- Bigger memory is slow"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Making algorithms more computationally intensive \n",
    "\n",
    "<font color='red'>**Implementation in NLA**</font>: use block version of algorithms. <br>\n",
    "\n",
    "This approach is a core of **[BLAS (Basic Linear Algebra Subroutines)](http://www.netlib.org/blas/)**, written in Fortran many years ago, and still rules the computational world."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Split the matrix into blocks! For illustration consider splitting in $2 \\times 2$ block matrix:\n",
    "\n",
    "$$\n",
    "   A = \\begin{bmatrix}\n",
    "         A_{11} & A_{12} \\\\\n",
    "         A_{21} & A_{22}\n",
    "        \\end{bmatrix}, \\quad B = \\begin{bmatrix}\n",
    "         B_{11} & B_{12} \\\\\n",
    "         B_{21} & B_{22}\n",
    "        \\end{bmatrix}$$\n",
    "\n",
    "Then,  \n",
    "\n",
    "$$AB = \\begin{bmatrix}A_{11} B_{11} + A_{12} B_{21} & A_{11} B_{12} + A_{12} B_{22} \\\\\n",
    "            A_{21} B_{11} + A_{22} B_{21} & A_{21} B_{12} + A_{22} B_{22}\\end{bmatrix}.$$  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "If $A_{11}, B_{11}$ and their product fit into the cache memory (which is 12 Mb (L3) for the [recent Intel Chip](https://en.wikipedia.org/wiki/List_of_Intel_microprocessors#Desktop)), then we load them only once into the memory.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## BLAS\n",
    "BLAS has three levels:\n",
    "1. BLAS-1, operations like $c = a + b$\n",
    "2. BLAS-2, operations like matrix-by-vector product\n",
    "3. BLAS-3, matrix-by-matrix product\n",
    "\n",
    "What is the principal differences between them?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The main difference is the number of operations vs. the number of input data!\n",
    "\n",
    "1. BLAS-1: $\\mathcal{O}(n)$ data, $\\mathcal{O}(n)$ operations\n",
    "2. BLAS-2: $\\mathcal{O}(n^2)$ data, $\\mathcal{O}(n^2)$ operations\n",
    "3. BLAS-3: $\\mathcal{O}(n^2)$ data, $\\mathcal{O}(n^3)$ operations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Why BLAS is so important and actual?\n",
    "\n",
    "1. The state-of-the-art implementation of the basic linear algebra operations\n",
    "2. Provides standard names for operations in any new implementations (e.g. [ATLAS](https://www.netlib.org/atlas/), [OpenBLAS](https://www.openblas.net/), [MKL](https://software.intel.com/en-us/mkl)). You can call matrix-by-matrix multiplication function (GEMM), link your code with any BLAS implementation and it will work correctly\n",
    "3. Formulate new algorithms in terms of BLAS operations\n",
    "4. There are wrappers for the most popular languages"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Packages related to BLAS\n",
    "\n",
    "1. [ATLAS](http://math-atlas.sourceforge.net) - Automatic Tuned Linear Algebra Software. It automatically adapts to a particular system architechture. \n",
    "2. [LAPACK](http://www.netlib.org/lapack/) - Linear Algebra Package. It provides high-level linear algebra operations (e.g. matrix factorizations), which are based on calls of BLAS subroutines.\n",
    "3. [Intel MKL](https://software.intel.com/en-us/intel-mkl) - Math Kernel Library. It provides re-implementation of BLAS and LAPACK, optimized for Intel processors. Available in Anaconda Python distribution: \n",
    "```\n",
    "conda install mkl\n",
    "```\n",
    "MATLAB uses Intel MKL by default.\n",
    "\n",
    "4. OpenBLAS is an optimized BLAS library based on [GotoBLAS](https://en.wikipedia.org/wiki/GotoBLAS). \n",
    "\n",
    "5. PyTorch [supports](https://pytorch.org/docs/stable/torch.html#blas-and-lapack-operations) some calls from BLAS and LAPACK\n",
    "\n",
    "6. For GPU it was implemented special [cuBLAS](https://docs.nvidia.com/cuda/cublas/index.html).\n",
    "\n",
    "\n",
    "For comparison of OpenBLAS and Intel MKL, see this [review](https://software.intel.com/en-us/articles/performance-comparison-of-openblas-and-intel-math-kernel-library-in-r)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Faster algorithms for matrix multiplication\n",
    "\n",
    "Recall that matrix-matrix multiplication costs $\\mathcal{O}(n^3)$ operations.\n",
    "However, storage is $\\mathcal{O}(n^2)$.\n",
    "\n",
    "**Question:** is it possible to reduce number operations down to $\\mathcal{O}(n^2)$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Answer**: a quest for $\\mathcal{O}(n^2)$ matrix-by-matrix multiplication algorithm is not yet done."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "* Strassen gives $\\mathcal{O}(n^{2.807\\dots})$ –– sometimes used in practice \n",
    "\n",
    "* [Current world record](http://arxiv.org/pdf/1401.7714v1.pdf) $\\mathcal{O}(n^{2.37\\dots})$ –– big constant, not practical, based on [Coppersmith-Winograd_algorithm](https://en.wikipedia.org/wiki/Coppersmith%E2%80%93Winograd_algorithm).\n",
    "- It improved the previous record (Williams 2012) by $3\\cdot 10^{-7}$\n",
    "- The papers still study multiplication of $3 \\times 3$ matrices and interpret it from different sides ([Heule, et. al. 2019](https://arxiv.org/pdf/1905.10192.pdf)) \n",
    "\n",
    "Consider Strassen in more details. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Na&iuml;ve multiplication\n",
    "\n",
    "Let $A$ and $B$ be two $2\\times 2$ matrices. Na&iuml;ve multiplication $C = AB$\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix} c_{11} & c_{12} \\\\ c_{21} & c_{22}  \\end{bmatrix}  =\n",
    "\\begin{bmatrix} a_{11} & a_{12} \\\\ a_{21} & a_{22}  \\end{bmatrix}\n",
    "\\begin{bmatrix} b_{11} & b_{12} \\\\ b_{21} & b_{22}  \\end{bmatrix} =\n",
    "\\begin{bmatrix} \n",
    "a_{11}b_{11} + a_{12}b_{21} & a_{11}b_{21} + a_{12}b_{22} \\\\ \n",
    "a_{21}b_{11} + a_{22}b_{21} & a_{21}b_{21} + a_{22}b_{22} \n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "contains $8$ multiplications and $4$ additions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Strassen algorithm\n",
    "\n",
    "In the work [Gaussian elimination is not optimal](http://link.springer.com/article/10.1007%2FBF02165411?LI=true) (1969) Strassen found that one can calculate $C$ using 18 additions and only 7 multiplications:\n",
    "$$\n",
    "\\begin{split}\n",
    "c_{11} &= f_1 + f_4 - f_5 + f_7, \\\\\n",
    "c_{12} &= f_3 + f_5, \\\\\n",
    "c_{21} &= f_2 + f_4, \\\\\n",
    "c_{22} &= f_1 - f_2 + f_3 + f_6,\n",
    "\\end{split}\n",
    "$$\n",
    "where\n",
    "$$\n",
    "\\begin{split}\n",
    "f_1 &= (a_{11} + a_{22}) (b_{11} + b_{22}), \\\\\n",
    "f_2 &= (a_{21} + a_{22}) b_{11}, \\\\\n",
    "f_3 &= a_{11} (b_{12} - b_{22}), \\\\\n",
    "f_4 &= a_{22} (b_{21} - b_{11}), \\\\\n",
    "f_5 &= (a_{11} + a_{12}) b_{22}, \\\\\n",
    "f_6 &= (a_{21} - a_{11}) (b_{11} + b_{12}), \\\\\n",
    "f_7 &= (a_{12} - a_{22}) (b_{21} + b_{22}).\n",
    "\\end{split}\n",
    "$$\n",
    "\n",
    "Fortunately, these formulas hold even if $a_{ij}$ and $b_{ij}$, $i,j=1,2$ are block matrices.\n",
    "\n",
    "Thus, Strassen algorithm looks as follows. \n",
    "- First of all we <font color='red'>split</font> matrices $A$ and $B$ of sizes $n\\times n$, $n=2^d$ <font color='red'> into 4 blocks</font> of size $\\frac{n}{2}\\times \\frac{n}{2}$\n",
    "- Then we <font color='red'>calculate multiplications</font> in the described formulas <font color='red'>recursively</font>\n",
    "\n",
    "This leads us again to the **divide and conquer** idea."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Complexity of the Strassen algorithm\n",
    "\n",
    "#### Number of multiplications\n",
    "\n",
    "Calculation of number of multiplications is a trivial task. Let us denote by $M(n)$ number of multiplications used to multiply 2 matrices of sizes $n\\times n$ using the divide and conquer concept.\n",
    "Then for na&iuml;ve algorithm we have number of multiplications\n",
    "\n",
    "$$ M_\\text{naive}(n) = 8 M_\\text{naive}\\left(\\frac{n}{2} \\right) = 8^2 M_\\text{naive}\\left(\\frac{n}{4} \\right)\n",
    "= \\dots = 8^{d-1} M(1) = 8^{d} = 8^{\\log_2 n} = n^{\\log_2 8} = n^3 $$\n",
    "\n",
    "So, even when using divide and coquer idea we can not be better than $n^3$.\n",
    "\n",
    "Let us calculate number of multiplications for the Strassen algorithm:\n",
    "\n",
    "$$ M_\\text{strassen}(n) = 7 M_\\text{strassen}\\left(\\frac{n}{2} \\right) = 7^2 M_\\text{strassen}\\left(\\frac{n}{4} \\right)\n",
    "= \\dots = 7^{d-1} M(1) = 7^{d} = 7^{\\log_2 n} = n^{\\log_2 7} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Number of additions\n",
    "\n",
    "There is no point to estimate number of addtitions $A(n)$ for naive algorithm, as we already got $n^3$ multiplications.  \n",
    "For the Strassen algorithm we have:\n",
    "\n",
    "$$ A_\\text{strassen}(n) = 7 A_\\text{strassen}\\left( \\frac{n}{2} \\right) + 18 \\left( \\frac{n}{2} \\right)^2 $$\n",
    "\n",
    "since on the first level we have to add $\\frac{n}{2}\\times \\frac{n}{2}$ matrices 18 times and then go deeper for each of the 7 multiplications. Thus,\n",
    "\n",
    "<font size=2.0>\n",
    "\n",
    "$$ \\begin{split}\n",
    "A_\\text{strassen}(n) =& 7 A_\\text{strassen}\\left( \\frac{n}{2} \\right) + 18 \\left( \\frac{n}{2} \\right)^2 = 7 \\left(7 A_\\text{strassen}\\left( \\frac{n}{4} \\right) + 18 \\left( \\frac{n}{4} \\right)^2 \\right) + 18 \\left( \\frac{n}{2} \\right)^2 =\n",
    "7^2 A_\\text{strassen}\\left( \\frac{n}{4} \\right) + 7\\cdot 18 \\left( \\frac{n}{4} \\right)^2 +  18 \\left( \\frac{n}{2} \\right)^2 = \\\\\n",
    "=& \\dots = 18 \\sum_{k=1}^d 7^{k-1} \\left( \\frac{n}{2^k} \\right)^2 = \\frac{18}{4} n^2 \\sum_{k=1}^d \\left(\\frac{7}{4} \\right)^{k-1} = \\frac{18}{4} n^2 \\frac{\\left(\\frac{7}{4} \\right)^d - 1}{\\frac{7}{4} - 1} = 6 n^2 \\left( \\left(\\frac{7}{4} \\right)^d - 1\\right) \\leqslant 6 n^2 \\left(\\frac{7}{4} \\right)^d = 6 n^{\\log_2 7}\n",
    "\\end{split}\n",
    "$$\n",
    "</font>\n",
    "\n",
    "(since $4^d = n^2$ and $7^d = n^{\\log_2 7}$).\n",
    "\n",
    "\n",
    "Asymptotic behavior of $A(n)$ could be also found from the [master theorem](https://en.wikipedia.org/wiki/Master_theorem)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Total complexity\n",
    "\n",
    "Total complexity is $M_\\text{strassen}(n) + A_\\text{strassen}(n)=$ <font color='red'>$7 n^{\\log_2 7}$</font>. Strassen algorithm becomes faster\n",
    "when\n",
    "\n",
    "\\begin{align*}\n",
    "2n^3 &> 7 n^{\\log_2 7}, \\\\\n",
    "n &> 667,\n",
    "\\end{align*}\n",
    "\n",
    "so it is not a good idea to get to the bottom level of recursion."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Strassen algorithm reloaded\n",
    "\n",
    "- Recent paper [Strassen algorithm reloaded](http://jianyuhuang.com/papers/sc16.pdf)\n",
    "claim to **break conventional wisdom** that Strassen algorithm is not very practical."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- Conventional wisdom: it is only\n",
    " practical for very large matrices. The proposed implementation is practical\n",
    " for small matrices. \n",
    "- Conventional wisdom: the matrices being\n",
    " multiplied should be relatively square. The proposed implementation is\n",
    " practical for rank-$k$ updates, where $k$ is relatively small (a shape\n",
    " of importance for libraries like LAPACK). \n",
    "- Conventional wisdom:\n",
    "  it inherently requires substantial workspace. The proposed implementation\n",
    "  requires no workspace beyond buffers already incorporated\n",
    "  into conventional high-performance DGEMM implementations.\n",
    "- Conventional wisdom: a Strassen DGEMM interface must pass\n",
    " in workspace. The proposed implementation requires no such workspace\n",
    " and can be plug-compatible with the standard DGEMM interface.\n",
    "- Conventional wisdom: it is hard to demonstrate speedup\n",
    " on multi-core architectures. The proposed implementation demonstrates\n",
    " speedup over conventional DGEMM even on an IntelR Xeon\n",
    " PhiTM coprocessor utilizing 240 threads. It is shown how a distributed\n",
    " memory matrix-matrix multiplication also benefits from\n",
    " these advances."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Strassen algorithm and tensor rank (advanced topic)\n",
    "\n",
    "- It is not clear how Strassen found these formulas. \n",
    "- However, now we can see that they are not artificial.\n",
    "- There is a general approach based on the so-called tensor decomposition technique. \n",
    "- Here by tensor we imply a multidimensional array - generalization of the matrix concept to many dimensions.\n",
    "\n",
    "Let us enumerate elements in the $2\\times 2$ matrices as follows\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix} c_{1} & c_{3} \\\\ c_{2} & c_{4}  \\end{bmatrix} =\n",
    "\\begin{bmatrix} a_{1} & a_{3} \\\\ a_{2} & a_{4}  \\end{bmatrix}\n",
    "\\begin{bmatrix} b_{1} & b_{3} \\\\ b_{2} & b_{4}  \\end{bmatrix}=\n",
    "\\begin{bmatrix} \n",
    "a_{1}b_{1} + a_{3}b_{2} & a_{1}b_{3} + a_{3}b_{4} \\\\ \n",
    "a_{2}b_{1} + a_{4}b_{2} & a_{2}b_{3} + a_{4}b_{4} \n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "This can be written as\n",
    "\n",
    "$$ c_k = \\sum_{i=1}^4 \\sum_{j=1}^4 x_{ijk} a_i b_j, \\quad k=1,2,3,4 $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "$x_{ijk}$ is a 3-dimensional array, that consists of zeros and ones:\n",
    "\n",
    "$$\n",
    "\\begin{split}\n",
    "x_{\\ :,\\ :,\\ 1} = \n",
    "\\begin{pmatrix}\n",
    "1 & 0 & 0 & 0 \\\\\n",
    "0 & 0 & 0 & 0 \\\\\n",
    "0 & 1 & 0 & 0 \\\\\n",
    "0 & 0 & 0 & 0 \\\\\n",
    "\\end{pmatrix}\n",
    "\\quad\n",
    "x_{\\ :,\\ :,\\ 2} = \n",
    "\\begin{pmatrix}\n",
    "0 & 0 & 0 & 0 \\\\\n",
    "1 & 0 & 0 & 0 \\\\\n",
    "0 & 0 & 0 & 0 \\\\\n",
    "0 & 1 & 0 & 0 \\\\\n",
    "\\end{pmatrix} \\\\\n",
    "x_{\\ :,\\ :,\\ 3} = \n",
    "\\begin{pmatrix}\n",
    "0 & 0 & 1 & 0 \\\\\n",
    "0 & 0 & 0 & 0 \\\\\n",
    "0 & 0 & 0 & 1 \\\\\n",
    "0 & 0 & 0 & 0 \\\\\n",
    "\\end{pmatrix}\n",
    "\\quad\n",
    "x_{\\ :,\\ :,\\ 4} = \n",
    "\\begin{pmatrix}\n",
    "0 & 0 & 1 & 0 \\\\\n",
    "0 & 0 & 0 & 0 \\\\\n",
    "0 & 0 & 0 & 0 \\\\\n",
    "0 & 0 & 0 & 1 \\\\\n",
    "\\end{pmatrix}\n",
    "\\end{split}\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Trilinear decomposition\n",
    "\n",
    "To get Strassen algorithm we should do the following trick –– decompose $x_{ijk}$ in the following way\n",
    "\n",
    "$$ x_{ijk} = \\sum_{\\alpha=1}^r u_{i\\alpha} v_{j\\alpha} w_{k\\alpha}. $$\n",
    "\n",
    "This decomposition is called **trilinear tensor decomposition** and has a meaning of separation of variables: we have a sum of $r$ (called rank) summands with separated $i$, $j$ and $k$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Strassen via trilinear\n",
    "\n",
    "Now we have\n",
    "\n",
    "$$ c_k = \\sum_{\\alpha=1}^r w_{k\\alpha} \\left(\\sum_{i=1}^4  u_{i\\alpha} a_i \\right) \\left( \\sum_{j=1}^4 v_{j\\alpha} b_j\\right), \\quad k=1,2,3,4. $$\n",
    "\n",
    "Multiplications by $u_{i\\alpha}$ or $v_{j\\alpha}$ or $w_{k\\alpha}$ do not require recursion since $u, v$ and $w$ are known precomputed matrices. Therefore, we have only $r$ multiplications of $\\left(\\sum_{i=1}^4  u_{i\\alpha} a_i \\right)$ $\\left( \\sum_{j=1}^4 v_{j\\alpha} b_j\\right)$ where both factors depend on the input data.\n",
    " \n",
    "As you might guess array $x_{ijk}$ has rank $r=7$, which leads us to $7$ multiplications and to the Strassen algorithm!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Summary of MM part\n",
    "- MM is the core of NLA. You have to think in block terms, if you want high efficiency\n",
    "- This is all about computer memory hierarchy\n",
    "- Concept of block algorithms\n",
    "- (Advanced topic) Strassen and trilinear form"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>\n",
       "<style>\n",
       "    @font-face {\n",
       "        font-family: \"Computer Modern\";\n",
       "        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "    }\n",
       "    div.cell{\n",
       "        /*width:80%;*/\n",
       "        /*margin-left:auto !important;\n",
       "        margin-right:auto;*/\n",
       "    }\n",
       "    h1 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\n",
       "    h2 {\n",
       "        font-family: 'Fenix', serif;\n",
       "    }\n",
       "    h3{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "        margin-top:12px;\n",
       "        margin-bottom: 3px;\n",
       "       }\n",
       "\th4{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "       }\n",
       "    h5 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\t   \n",
       "    div.text_cell_render{\n",
       "        font-family: 'Alegreya Sans',Computer Modern, \"Helvetica Neue\", Arial, Helvetica, Geneva, sans-serif;\n",
       "        line-height: 1.2;\n",
       "        font-size: 120%;\n",
       "        /*width:70%;*/\n",
       "        /*margin-left:auto;*/\n",
       "        margin-right:auto;\n",
       "    }\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "    .text_cell_render h1 {\n",
       "        font-weight: 200;\n",
       "        font-size: 50pt;\n",
       "\t\tline-height: 110%;\n",
       "        color:#CD2305;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\t\n",
       "    .text_cell_render h5 {\n",
       "        font-weight: 300;\n",
       "        font-size: 16pt;\n",
       "        color: #CD2305;\n",
       "        font-style: italic;\n",
       "        margin-bottom: .5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\n",
       "    \n",
       "    li {\n",
       "        line-height: 110%;\n",
       "    }\n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "\n",
       "</style>\n",
       "\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"]\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"../styles/custom.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python (pytorch)",
   "language": "python",
   "name": "pytorch"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.10"
  },
  "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
}