TB-Lecture-08-julia.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 8: Gram-Schmidt orthogonalization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modified Gram-Schmidt\n",
    "\n",
    "The classical Gram-Schmidt process is numerically unstable, as will be seen in the next lecture. A mathematically equivalent formulation, which we call MGS, turns out to be a lot better. \n",
    "\n",
    "In classical GS the key step is to project a new column of $A$ onto the orthogonal complement of the span of all the previous $q_j$s. We express this projection as $P_j=I - Q_jQ_j^*$, where $Q_j=\\bigl[ q_1 \\; q_2\\; \\cdots \\; q_j\\bigr]$. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6-element Array{Float64,1}:\n",
       " 1.0        \n",
       " 1.0        \n",
       " 1.0        \n",
       " 5.58584e-16\n",
       " 2.94898e-16\n",
       " 6.3158e-17 "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Qj = qr(rand(6,3)[1];   # get a random Q with 3 ON columns\n",
    "Pj = eye(6) - Qj*Qj';\n",
    "σ = svd(Pj)[2]       # all 1 or 0, as required for an orthogonal projector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's intuitively clear that we get the same result by subtracting off projections onto each column, one at a time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.2619003506203114e-16"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Mj = eye(6) - Qj[:,1]*Qj[:,1]' - Qj[:,2]*Qj[:,2]' - Qj[:,3]*Qj[:,3]';\n",
    "norm(Pj-Mj)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also factor this into a product of orthogonal projectors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.862094358590312e-16"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Gj = eye(6);\n",
    "for k = 1:3\n",
    "    Gj = Gj*(eye(6)-Qj[:,k]*Qj[:,k]');\n",
    "end\n",
    "norm(Gj-Pj)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This works because the \"cross-terms\" you get in the products have inner products between different $q_k$ and thus are zero. In MGS, we apply each of these factors, as soon at it is found, to \\emph{all} the columns of $A$ at once. By the time we are ready to choose a new column from $A$, it is already orthgonalized against the entire ON basis set found so far."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Operation counts\n",
    "\n",
    "It's standard in numerical analysis to compare the performance of algorithms by counting the number of floating point operations, or flops, that are performed as a function of the input size. By either form of QR factorization, the number of flops is $\\sim 2mn^2$ in asymptotic notation. Let's see if we observe the implicit proportionality with respect to $n^2$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "elapsed time: 0.026677532 seconds\n",
      "elapsed time: 0.09042511 seconds\n",
      "elapsed time: 0.110793201 seconds\n",
      "elapsed time: 0.217990293 seconds\n",
      "elapsed time: 0.347987628 seconds\n",
      "elapsed time: 0.483492849 seconds\n",
      "elapsed time: 0.687628253 seconds\n",
      "elapsed time: 0.914240892 seconds\n",
      "elapsed time: 1.125613474 seconds\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtAAAAIxCAYAAABkT4+7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XlYlXX+//HXOSyuhJm7WVC2UOoYpGmaTaaImlip1cmt8tLKzKJxcgy3iumXM5riN7LFRsPlaGmp5EIumblhA2O2qJmpFZZGKSIqIOf+/XGSNLE459xwtufjurqczn14328axRcfPvfnbTEMwxAAAACACrF6uwEAAADAnxCgAQAAABcQoAEAAAAXEKABAAAAFxCgAQAAABcQoAEAAAAXEKABAAAAFxCgAQAAABcQoAEAAAAXBHWAfvXVVxUXF6fw8HA999xz3m4HAAAAfiCoA3STJk307LPPqm/fvt5uBQAAAH4i1NsNeFNiYqIkafny5V7uBAAAAP7CL1agCwsLNWHCBHXv3l2XXHKJrFar0tPTy31vcXGxRo8eraZNm6pmzZpq166d1qxZU8UdAwAAIFD5RYDOy8vT888/r127dql169ayWCwXfO/gwYM1bdo0DRw4UNOnT1doaKh69OihzZs3V2HHAAAACFR+sYWjSZMm+vHHH9WgQQNlZ2erTZs25b5v27ZtWrhwoaZMmaKkpCRJ0sCBA9WiRQs9/fTT2rhxY1W2DQAAgADkFyvQYWFhatCgwZ++b9GiRQoNDdXQoUPLXqtWrZqGDBmiLVu2KDc3tzLbBAAAQBDwiwBdUdu3b9fVV1+t2rVrn/N627Zty66frbS0VKdOnVJpaalKSkpUVFQkh8NRZf0CAADA/wRUgP7hhx/UuHHj815v3LixDMPQwYMHz3k9JSVFNWvW1JtvvqkXXnhBNWvW1Ny5c6uqXQAAAPghv9gDXVEnT55UtWrVznu9evXqZdfPNmHCBE2YMKFCtfPy8pSZmamoqCjVqFHD82YBAABgqpMnT2r//v3q1q2b6tWrV2n3CagAXaNGDRUVFZ33+qlTp8quuyszM1MDBgxw++MBAABQNebOnav+/ftXWv2ACtCNGzc+b5uG5NzaITlP83BXVFSUJOf/ITExMW7XgWuSkpI0depUb7dRKXz5c/NWb5V938qob1ZNT+u4+/G+/PswUAXyf3Nf/tz4ulb1Nb3xdW3nzp0aMGBAWW6rLAEVoFu3bq3169fr+PHj5zxIuHXrVlksFrVu3drt2mdWr2NiYhQbG+txr6iYyMjIgP3v7cufm7d6q+z7VkZ9s2p6Wsfdj/fl34eBKpD/m/vy58bXtaqv6a2va5Jnuw4qIqAeIuzbt69Onz6t119/vey14uJizZ49W+3atVPTpk292B3cYbPZvN1CpfHlz81bvVX2fSujvlk1Pa3jy7+fcK5A/v/Klz83vq5Vfc1A/rpmMQzD8HYTFZGWlqajR48qNzdXr776qu6++27dcMMNkqSRI0cqIiJCknTvvfdqyZIlevLJJ9W8eXPNnj1b//3vf7Vu3Tp16NDB7fvn5OQoLi5O2dnZPvvdNQC4KjExUcuWLfN2GwBgiqrKa36zhWPy5Mn69ttvJUkWi0Xvvfee3nvvPUnOaYNnAvScOXM0btw4zZ07V0eOHFGrVq20fPlyj8IzAAAAcIbfBOh9+/ZV6H3h4eGaNGmSJk2aVMkdAYD/8+UfkQKAr/KbAO0rkpKSFBkZKZvNxl88APweX8cABAK73S673a78/PwquZ/f7IH2NvZAAwAA+LaqymsBdQoHAAAAUNkI0AAAAIALCNAAAACACwjQAAAAgAsI0AAAAIALCNAAAACACzgH2kWcAw0AAOBbOAfaR3EONAAAgG/jHGgAAADABxGgAQAAABcQoAEAAAAXEKABAAAAFxCgAQAAABcQoAEAAAAXEKABAAAAFzBIxUUMUgEAAPAtDFLxUQxSAQAA8G0MUgEAAAB8EAEaAAAAcAEBGgAAAHABARoAAABwAQEaAAAAcAEBGgAAAHABARoAAABwAQEaAAAAcAEBGgAAAHABo7xdxChvAAAA38Iobx/FKG8AAADfxihvAAAAwAcRoAEAAAAXEKABAAAAFxCgAQAAABcQoAEAAAAXEKABAAAAFxCgAQAAABcQoAEAAAAXEKABAAAAFxCgAQAAABcQoAEAAAAXhHq7AX+TlJSkyMhI2Ww22Ww2b7cDAAAQ9Ox2u+x2u/Lz86vkfhbDMIwquZOfy8nJUVxcnLKzsxUbG+vtdgAAAPA7VZXX2MIBAAAAuIAADQAAALiAAA0AAAC4gAANAAAAuIAADQAAALiAAA0AAAC4gAANAAAAuIAADQAAALiAAA0AAAC4gAANAAAAuIAADQAAALiAAA0AAAC4gAANAAAAuCDU2w34m6SkJEVGRspms8lms3m7HQAAgKBnt9tlt9uVn59fJfezGIZhVMmd/FxOTo7i4uKUnZ2t2NhYb7cDAACA36mqvMYWDgAAAMAFBGgAAADABQRoAAAAwAUEaAAAAMAFBGgAAADABQRoAAAAwAUEaAAAAMAFBGgAAADABQRoAAAAwAUEaAAAAMAFBGgAAADABQRoAAAAwAUEaAAAAMAFBGgAAADABQRoAAAAwAUEaAAAAMAFBGgAAADABQRoAAAAwAWh3m7A3yQlJSkyMlI2m002m83b7QAAAAQ9u90uu92u/Pz8KrmfxTAMo0ru5OdycnIUFxen7OxsxcbGersdAAAA/E5V5TW2cAAAAAAuIEADAAAALiBAAwAAAC4gQAMAAAAuIEADAAAALiBAAwAAAC4gQAMAAAAuIEADAAAALiBAAwAAAC4gQAMAAAAuIEADAAAALiBAAwAAAC4gQAMAAAAuIEADAAAALiBAAwAAAC4gQAMAAAAuIEADAAAALiBAAwAAAC4gQAMAAAAuIEADAAAALiBAAwAAAC4gQAMAAAAuIEADAAAALiBAAwAAAC4gQAMAAAAuIEADAAAALgj1dgP+JikpSZGRkbLZbLLZbN5uBwAAIOjZ7XbZ7Xbl5+dXyf0shmEYVXInP5eTk6O4uDhlZ2crNjbW2+0AAADgd6oqr7GFAwAAAHABARoAAABwAQEaAAAAcAEBGgAAAHABARoAAABwAQEaAAAAcAEBGgAAAHABARoAAABwAQEaAAAAcAEBGgAAAHABARoAAABwAQEaAAAAcAEBGgAAAHABARoAAABwAQEaAAAAcAEBGgAAAHABARoAAABwAQEaAAAAcAEBGgAAoJIYhuHtFlAJCNAAAAAmKigo0MiRExQd3UXNmt2p6OguGjlyggoKCrzdGkwS6u0GAAAAAkVBQYHat++jnTufksMxUZJFkqG0tEytW9dHW7YsVkREhJe7hKdYgQYAADBJcvLkX8NzgpzhWZIscjgStHNnksaOneLN9mASAjQAAIBJMjI2yeHoVu41hyNBy5ZtquKOUBkI0AAAACYwDEMlJbX028rz71lUUlKTBwsDAAEaAADABBaLRWFhhZIuFJANhYUVymK5UMCGvyBAAwAAmKRBgw6SMsu9ZrWuUmJix6ptCJWCUzgAAABM8PLL0rZto9SgQR/l5RlnPUhoyGpdpZiYqUpJWeztNmECVqABAAA8NH++9Pjj0lNPRWjPnsUaMSJLUVHxatq0t6Ki4jViRBZH2AUQVqABAAA8sHKlNHiw9MAD0uTJksUSodTUiUpNdT5YyJ7nwMMKNAAAgJs2bZL69JF69JDeeEP6fVYmPAcmAjQAAIAbduyQ7rhDattWWrBACuXn+kGDAA0AAOCib76RunWToqOlpUulGjW83RGqEgEaAADABT/+KHXtKkVESKtWSZGR3u4IVY0fNgAAAPyJMw8DHj3qXHk+dcq5/7lBA293Bm8gQAMAAJSjoKBAycmTlZGxSSUltRQSUqiSkg46eXKUNm6MUFSUtzuEtxCgAQAAfqegoEDt2/fRzp1PyeGYqDMDUaRMRUf30WWXLZbEmc7Bij3QAAAAv5OcPPnX8HxmmqB+/TVBBw4kaezYKV7sDt4W1AE6Ly9Pd9xxh2rXrq2YmBitW7fO2y0BAAAfkJGxSQ5Ht3KvORwJWrZsUxV3BF8S1Fs4hg8frsaNG+vnn3/WBx98oHvuuUdff/216tSp4+3WAACAlxiGoZKSWvpt5fn3LCopqcmUwSAWtCvQhYWFWrp0qZ577jlVq1ZNvXr1UqtWrbR06VJvtwYAALzIYrEoNLRQzj3P5TEUFlZIeA5ifhGgCwsLNWHCBHXv3l2XXHKJrFar0tPTy31vcXGxRo8eraZNm6pmzZpq166d1qxZc9779uzZo4iICDVu3LjstRYtWuiLL76otM8DAAD4vvx8qaSkg6TMcq9brauUmNixapuCT/GLAJ2Xl6fnn39eu3btUuvWrf/wO77Bgwdr2rRpGjhwoKZPn67Q0FD16NFDmzdvPud9x48f10UXXXTOaxdddJGOHz9eKZ8DAADwfT/+KN16q1RYOErR0S/Jal2p31aiDVmtKxUTM1UpKX/zZpvwMr/YA92kSRP9+OOPatCggbKzs9WmTZty37dt2zYtXLhQU6ZMUVJSkiRp4MCBatGihZ5++mlt3Lix7L21a9fWsWPHzvn4Y8eOqXbt2pX3iQAAAJ+1d68UH+8ckvLxxxGKilqssWOnaNmyl1RSUlNhYSeUmNhBKSmLFRHBEXbBzC8CdFhYmBpUYNTPokWLFBoaqqFDh5a9Vq1aNQ0ZMkTJycnKzc1V06ZNJUlXXXWVjh8/rh9++KFsG8dnn32mBx54oFI+BwAA4Lv+9z+pe3fnWO5Nm/TrkJQIpaZOVGqqeGAQ5/CLLRwVtX37dl199dXnrSK3bdu27PoZtWrVUu/evTVhwgSdOnVKGRkZ+vzzz9W7d+8q7RkAAHjX+vXObRvNmkkbN6rcCYOEZ5wtoAL02avJZ2vcuLEMw9DBgwfPeT0tLU25ubm65JJL9Pe//11vv/02R9gBABBE3n1X6tZNuukmad06qX59b3cEf+AXWzgq6uTJk6pWrdp5r1evXr3s+tnq1aun5cuXV0lvAADAt7z+uvToo1K/flJ6uhQe7u2O4C8CKkDXqFFDRUVF571+6tSpsuueSkpKUmRk5Dmv2Ww22Ww2j2sDAIDKZxjSP/8pjRsnjRghpaZK1oD6mXxwsNvtstvt57yWn59fJfcOqADduHHj87ZpSM6tHZLzNA9PTZ06VbGxsR7XAQAAVefMQ4AOh/TEE9LLL0vPPy8lJ0tmbm8+WXJSNcI8X7DDnytvATMnJ0dxcXGVfu+A+n6rdevW+uqrr847y3nr1q2yWCxq3bq1lzoDAABVraCgQCNHTlB0dBc1a3anoqK66NprJygtrUCvvSaNHWtueP7myDeKSYvRkl1LzCsKnxRQAbpv3746ffq0Xn/99bLXiouLNXv2bLVr167sCDsAABDYCgoK1L59H6Wltdf+/auVm7tUBw6s1p497XXppX1ksxWYer9v879V57c6q1poNbW7tJ2pteF7/GYLR1pamo4eParc3FxJ0rJly/Tdd99JkkaOHKmIiAi1bdtW/fr105gxY3To0CE1b95cs2fP1oEDBzRr1ixvtg8AAKpQcvJk7dz5lByOhLNetUhKUG6uobFjpyg1daIp98o9lqvOb3WW1WLV2kFr1ah2I1PqwndZDMMw/vxt3hcdHa1vv/223Gv79u3TZZddJsm54jxu3DjNnTtXR44cUatWrZSSkqIuXbp4dP8ze2qys7PZAw0AgI+Lju6i/ftXyxmaf89QVFS89u1b7fF9Dh0/pFtn36oTJSe04cENiqoT5XFNuK+q8prfrEDv27evQu8LDw/XpEmTNGnSpEruCAAA+CLDMFRSUkvlh2dJsqikpKbH0wXzTuSpy5wuOlZ0jPAcZPwmQPuKM8fYcXQdAAC+yWKxKDS0UJKhC61Ah4UVehSei04XKX5OvA4XHtZHD3yk5nWbu10LnjtzpB3H2PkojrEDAMC3ORxStWodJGVKSjjvutW6SomJHT26R7XQahpywxB1uryTrq13rUe14LkzC5tVdYwdARoAAAQMw5D+9jfpq69G6dJL++jgQePXBwktkgxZrasUEzNVKSmLPb7XY20f87gG/FNAHWMHAACC27/+JU2bJqWlRejLLxdrxIgsRUXFq2nT3oqKiteIEVnasmWxIiIivN0q/Bgr0AAAICDMmiX94x/S+PHS8OGSFKHU1IlKTZXHDwwCZ2MFGgAA+L2MDGnoUGnYMGnixPOvE55hJgI0AADwa5s3S/fcI/XuLb3yirnjuU87TuvU6VPmFURAIEADAAC/9cUX0h13SDfdJM2bJ4WEmFe71FGqwUsG666Fd8lP5s6hirAH2kWcAw0AgG/49lupWzfpssukpUul6tXNq+0wHBqWMUwLPl+gBX0WsAXEx1X1OdB+M8rb2xjlDQCA78jLkzp2lIqKnFs4Gjc2r7ZhGHpsxWN69b+vKv2udA1oNcC84qhUjPIGAAAoR2Ghc9vGL79ImzaZH56fynxKM/47QzN7zSQ8o1wEaAAA4DdKSqS+fZ17n9evl666yrzahmHombXPaFrWNL3c/WUNiR1iXnEEFAI0AADwCw6H9NBD0tq10ooVktkTm6dtnaYXN72oKfFTmDKIP0SABgAAfuHpp50nbdjtUpcu5tfvfW1vhYeEE57xpzjGDgAA+KwzZx38+9/SlClSaqp0772Vc68rLr6C8IwKYQUaAAD4lIKCAiUnT1ZGxiaVlNRSUVGh8vI6aNSoUXr88QhvtwewAg0AAHxHQUGB2rfvo7S09tq/f7Vyc5cqL2+1pPZaubKPCgoKvN0iQIB2VVJSkhITE2W3273dCgAAASc5ebJ27nxKDkeCpDPDSyySErRzZ5LGjp3ixe7gq+x2uxITE5WUlFQl92OQSgUxSAUAgMoXHd1F+/ev1m/h+WyGoqLitW/fao/ukX8qX5HVIz2qAd9UVXmNFWgAAOATDMNQSUktlR+eJcmikpKa8mTt7/2v3ldUapQ+/fFTt2sABGgAAOATLBaLQkMLJV0oIBsKCyuUxXKhgP3HPtj7gfq83Ue3Rd2m6+pf53afgNsB+ttvv9Ujjzyia665RnXr1tWGDRskSXl5eRo5cqT+97//mdYkAAAIfIYh1arVQVJmudet1lVKTOzoVu31+9frzgV3qusVXbWg7wKFhYR50CmCnVvH2H355Ze65ZZb5HA4dNNNN+nrr7/W6dOnJUn16tXTxo0bVVhYqDfffNPUZgEAQOCaPFn68stRatq0j374wTjrQUJDVusqxcRMVUrKYpfrbv5us+6Yf4c6XtZRi+5ZpPCQcNN7R3BxawX66aefVp06dfTVV19p7ty55+1F6tmzpz7++GNTGgQAAIFvyRJp9GhpzJgI7dy5WCNGZCkqKl5Nm/ZWVFS8RozI0pYtixUR4do50J/kfqLu87rrxiY3asl9S1Q9tHolfQYIJm6tQG/YsEHjx49X/fr19fPPP593/bLLLlNubq7HzQEAgMCXkyP17y/16SOlpEhWa4RSUycqNdX5YKG7e54PHD2gbnO76fr61yvDlqGaYTVN7hzByq0VaIfDoZo1L/yb8KefflK1atXcbgoAAASH77+XevWSrr9eeustyfq7ZOJueJakyyIv09hOY7Wy/0pFVGOCIczjVoCOjY3V8uXLy712+vRpLViwQO3atfOoMQAAENiOH3eG55AQaelS6Q/W5txisVj0VPunOPMZpnMrQI8ZM0arVq3So48+qs8//1ySdOjQIa1Zs0bx8fHauXOn/vGPf5jaKAAACBylpc5tG19/Lb3/vtS4sbc7AirOrT3Q3bt31+zZs/XEE0/o9ddflyQNGDBAhmHooosuUnp6ujp16mRqo74iKSlJkZGRstlsstls3m4HAAC/NHq0MzgvWya1auXtbuDv7Ha77Ha78vPzq+R+Ho3yLiws1AcffKCvv/5aDodDV155pbp16+byE7L+gFHeAACY4403pGHDpGnTpCee8HY3CCRVldfcWoE+o1atWrrrrrvM6gUAAAS4tWul4cOlRx+VRo40p+aRk0cUWT1SVgsDllE1PArQJSUlys3N1ZEjR8qdS89KLQAAOGPXLudRdZ07S9OnSx4csFEm70SebnvrNnW9oqte6vaS5wWBCnArQB89elSjRo3SvHnzVFxcfN71M2c2lpaWetwgAADwf3l5Us+eUtOm0ttvS6EeLeE5HTl5RPFz4nXo+CENjR3qeUGggtz67fvAAw8oIyND9913n2666SZFRnI8DAAAKF9RkXTXXVJBgbRmjWRGbDhWdEwJ8xJ0IP+APhz8oWLqx3heFKggtwL0Bx98oJEjR2rq1Klm9wMAAAKIYUhDh0qffCKtWydFR3te83jxcfWY10O783Zr3eB1atWQYzxQtdwK0JdccomaN29udi8AACDAvPCCNGeONG+edPPNntc7WXJSifZEfXroU60ZuEaxjXneClXPrcdVhw0bpgULFsjhcJjdDwAACBBvvy2NHStNmCDdf785NQcvGays3CytuH+Fbrr0JnOKAi5yawV63LhxKioq0o033qiBAwfq0ksvVUhIyHnvu/vuuz1uEAAA+J+sLGnwYMlmcwZos4y6eZQeufER3XL5LeYVBVzkVoDOzc3VunXrtH37dm3fvr3c93AKBwAAweXMKVwHDki9e0s33CD95z/mHFd3Rtumbc0rBrjJrQD90EMPKScnR2PGjOEUDgAAglhBQYGSkycrI2OTSkpqKSSkUIWFHVSr1igtWRKh6tW93SFgPrcC9MaNGzV69Gg9++yzZvcDAAD8REFBgdq376OdO5+SwzFRkkWSISlTF1/cRzVqLJYU4dUegcrg1kOEjRo1Ut26dc3uxS8kJSUpMTFRdrvd260AAOBVycmTfw3PCXKGZ/36a4K++SZJY8dO8WJ3CCZ2u12JiYlKSkqqkvtZjPJmcP+JV155RTNmzNCWLVtUu3btyujL5+Tk5CguLk7Z2dmMKAcAQFJ0dBft379av4XnsxmKiorXvn2rXa5rGIYOFx5Ww9oNPe4RwaWq8ppbWzhOnTqlsLAwNW/eXPfcc4+aNWt23ikcFoulyr4LAAAAVcswDJWU1FL54VmSLCopqVn2YKErdZPXJes///uPdo3YpTrV65jSL2AmtwL0qFGjyv73yy+/XO57CNAAAAQui8WisLBCOfc8l78CHRZW6FJ4lqSUDSn6fxv/n16Kf4nwDJ/lVoDet2+f2X0AAAA/cvq0VKNGB0mZkhLOu261rlJiYkeXav5r0780fv14vdD5BSW1ZxEOvsutAH355Zeb3QcAAPATRUXOyYK7d49S06Z99MMPxlkPEhqyWlcpJmaqUlIWV7hm6tZUjV4zWuM7jdeYW8ZUWu+AGdw6hQMAAASn48elnj2lFSukJUsitHPnYo0YkaWoqHg1bdpbUVHxGjEiS1u2LFZERMWOsHv1v6/qycwn9fTNT2viXydW7icAmKBCK9DR0dGyWq3atWuXwsLCFB0d/ad7miwWi/bu3WtKkwAAwPt++UXq0UP68ktp1Srp1lslKUKpqROVmiqXHxiUpMyvM/Xo8kc1su1IvdjlRZc/HvCGCgXoW2+9VRaLRVar9Zx/BwAAweHgQSk+Xjp0SPrwQyku7vz3uJMN/hr1V712x2saGjuUbAG/UaEAPXv27D/8dwAAELj27pW6dpVKSqSPP5auvda82tVCq2lY3DDzCgJVwK090Onp6dq/f/8Frx84cEDp6enu9gQAAHzEZ59JHTtKoaHSxo3mhmfAX7kVoB988EFt3rz5gte3bt2qBx980O2mAACA923d6tzn3KiRMzxzCBfg5FaA/rPp34WFhQoNdeuEPAAA4ANWr5Zuv126/npp/XqpQQNvdwT4jgqn3B07dmj79u1l//7xxx/r9OnT573v6NGjevXVV3X11Veb0yEAAKhSixdLNpvUpYu0aJFUs6bnNb8/9r0uvehSzwsBPqDCAfq9997Ts88+K8n5lO1rr72m1157rdz31qlThz3QAAD4of/8Rxo6VLrnHumtt6TwcM9rbsvdpi7pXTQtYZoeuuEhzwsCXlbhAD1s2DDdcccdMgxDbdu21XPPPafu3buf8x6LxaJatWrpyiuvDNgtHElJSYqMjJTNZpPNZvN2OwAAmGbKFGnUKOmRR6SXX5ZCQjyvuf3H7eo2t5taNmype66/x/OCQDnsdrvsdrvy8/Or5H4W4882NJfjo48+UkxMjBoE0YaonJwcxcXFKTs7W7Gxsd5uBwAA0xiGNHas9MIL0jPPSCkpkhlHMn9++HP9dfZfFX1xtNYMXKPI6pGeFwX+QFXlNbeWiW91jh4CAAB+zuGQRoyQZsyQ/v1v5wq0GXbn7VaX9C669KJLlTkgk/CMgBKY+ywAAMCfKimRBg+WFi6UZs6Uhgwxp+7eX/aqc3pnXVLzEq0euFp1a9Q1pzDgIwjQAAAEoRMnpH79nMfVLVwo9e1rTt3i0mIlzEtQ7fDaWjtorerXqm9OYcCHEKABAAgShmHIYrEoP1/q1UvKzpaWL3eO6TZLeEi4/q/7/6lFgxZqVLuReYUBH0KABgAggBUUFCg5ebIyMjappKSWrNZCnTzZQSUlo7RmTYTatzf/ngnNE8wvCvgQAjQAAAGqoKBA7dv30c6dT8nhmCjJIsmQlKkrr+yjFi0WS4rwao+AP6pQgH7uuedcLmyxWDRu3DiXPw4AAJgjOXnyr+H57BVhi6QE7dtnaOzYKUpNneil7gD/VaEAPXHixPNes/x6QOTvj5G2WCxle6wI0AAAeE9GxqZfV57P53AkaNmyl5SaWrU9AYHAWpE3ORyOc/757rvv1LJlS9lsNm3btk35+fnKz89XVlaW7rvvPv3lL3/Rd999V9m9AwCACzAMQyUlteRccS6PRSUlNc9bCKuowuJCHS487HZ/gD+rUID+vccee0xXXXWV5s6dqxtvvFERERGKiIhQmzZtNG/ePF155ZV67LHHzO4VAABUkMVikdVaKOee5/IYCgsrLPuJsitOlpxUL3sv3TH/DrcDOODP3ArQ69afbKlCAAAgAElEQVStU+fOnS94/fbbb9fatWvdbgoAAHjm6FHpxIkOkjLLvW61rlJiYkeX6xadLtJdC+9SVm6WpsRPcSuAA/7OrQBdvXp1bdmy5YLXN2/erOrVq7vdFAAAcN+JE1JionT69ChdeeVLslpX6reVaENW60rFxExVSsrfXKpbXFqsfu/000cHPlKGLUO3XH6L6b0D/sCtAN2/f3/NmzdPI0eO1J49e8r2Ru/Zs0ePP/645s+fr/79+5vdKwAA+BMlJdI99ziHpKxcGaH//W+xRozIUlRUvJo27a2oqHiNGJGlLVsWKyKi4kfYnXacVv93+ytzb6beu/c9dY6+8E+igUDn1jnQkyZNUl5enl5++WWlpaXJanXmcIfDIcMwZLPZNGnSJFMbBQAAf8zhkB56SPrgAykjQ78OSYlQaupEpab+NonQVaWOUg1eMlhLdi3Ron6LGJSCoOdWgA4PD9ecOXP097//XcuXL9e3334rSbr88svVvXt3/eUvfzG1SQAA8McMQ0pKkubNk+x2qVu389/j7n7lFze+qAWfL9CCPgvU+9reHnYK+D+PJhG2atVKrVq1MqsXAADgppQUafp0acYM6d57za39WNvH1LpRa/W8uqe5hQE/5VGA3rp1qz788EMdPnxYw4cP11VXXaUTJ05o165duvrqq1W7dm2z+gQAABfwyivS+PHOEP3II+bXr1O9DuEZOItbDxEWFxfr7rvvVocOHZScnKzp06eXDU6xWq2Kj49XKqONAACodHa7NGKE9OST0jPPeLsbIDi4FaDHjRun999/XzNmzNDu3bvPOUS9evXq6tevn5YuXWpak74kKSlJiYmJstvt3m4FABDkVq6UBg2SBg6UpkyROJIZwcputysxMVFJSUlVcj+3tnDY7XY9+uijGjZsmH7++efzrsfExOidd97xuDlfNHXqVMXGxnq7DQBAkNu0SerTR+reXZo5U7K6tSQGBAabzSabzaacnBzFxcVV+v3c+uN2+PBhtWzZ8oLXQ0JCdOLECbebAgAAF7Zjh3THHVKbNtLChVJYmDl1d+ftNqcQEODcCtDNmjXTrl27Lnh906ZNat68udtNAQCA8u3d6zyiLjpaWrZMqlHDnLqv/fc1XffKdcr6PsucgkAAcytA33///XrttdfOGed95mzJN954Q2+//bYGDRpkTocAAECS9MMPUny8FBEhrVolRUaaU3fW/2bpkeWP6LE2j6lt07bmFAUCmFt7oJOTk7V161Z16tRJMTExslgsSkpK0i+//KLvv/9ePXr0qLJN3AAABIMjR5wrz0VFzv3PDRqYU3f+Z/M1ZNkQDYsdptSEVLeHrQDBxK0V6PDwcK1atUqzZs3SFVdcoWuvvVZFRUVq1aqVZs+erYyMDIWEhJjdKwAAQamw0LnnOTfXOab78svNqbv4y8Ua9N4gDfrLIM24YwbhGaggtwepWCwWDRgwQAMGDDCzHwAAcJbiYqlvX+nTT6V166TrrjOnbsbuDN23+D71u76f3kx8U1YLx3gAFeXRJMKzGYahDz/8UEVFRerYsaMiIiLMKg0AQFByOKQHHpDWrpVWrJDamrQ9+bv879TvnX7qdXUvpd+ZrhArPzUGXOHWt5vJycm67bbbyv7dMAzFx8era9eu6tmzp1q2bKm9e/ea1iQAAMHGMKTHH3ceUzd/vtSli3m1m0U20zv93tGCvgsUFmLSGXhAEHErQC9evFhtz/o2eNGiRVq7dq1SUlL0/vvvq7S0VBMnTjSrRwAAgs7EidIrr0ivvurcwmG2Xtf0UnhIuPmFgSDg1haO3Nzcc855fvfdd3XddddpzJgxkqRHH31UM2bMMKdDAACCzPTp0nPPSS++KA0d6u1uAPyeWyvQoaGhKioqkuTcvrF27VolJCSUXW/YsKHy8vLM6RAAgCAyd670xBPSqFHS0097uxsA5XErQLdo0UJz587VkSNHNGvWLP3888/q2bNn2fUDBw6oXr16pjUJAEAweP9950ODDz0k/etfEqfKAb7JrS0c48ePV69evcpCcocOHc55qHD58uVq06aNOR0CABAENmyQ+vWTEhOl114zJzx/c+QbNazVULXCa3leDEAZtwJ0165dlZOTo9WrV6tOnTq69957y64dOXJEnTp1Uu/evU1rEgCAQLZ9u9Srl3Tzzc4TN0JNOGT2myPfqNOsTkponqCZiTM9LwigjNt/RK+77jpdV85p7hdffLGmTp3qUVMAAASLPXucI7qvvlpaskSqXt3zmgeOHlDntzqrVngtpXRO8bwggHN49D3u559/rhUrVmj//v2SpKioKHXv3l0tW7Y0ozcAAAKOYRhlI7Nzc6WuXaWLL3YOSjFjBlnusVzdnn67rBar1g5aq0a1G3leFMA53ArQRUVFevjhhzVnzhwZhiGr1fksosPh0JgxY9S/f3/NnDlT4eGcLwkAQEFBgZKTJysjY5NKSmopLKxQXbt20Mcfj5LDEaHVq6X69T2/z6Hjh3R7+u0qLi3Whgc36NKLLvW8KIDzuBWgR48erfT0dA0fPlyPP/64rrzySlksFn399deaPn26ZsyYobp162ratGlm9wsAgF8pKChQ+/Z9tHPnU3I4JkqySDL0xhuZCgnpo6ysxWrWzPOl57wTeeoyp4uOFR3Thgc3KKpOlMc1AZTPrWPs5s6dq4EDB+rll1/WNddco9DQUIWEhOiaa65RWlqa+vfvr7lz55rdKwAAfic5efKv4TlBzvCsX39NkGEkKT19iin3GZoxVIcLD2vd4HVqXrf5n38AALe5FaBLSkrUrl27C16/+eabdfr0abebAgAgUGRkbJLD0a3caw5HgpYt22TKfVITUrV20FpdW+9aU+oBuDC3AnS3bt2UmZl5weurVq1SfHy8200BABAIDMNQSUkt/bby/HsWlZTUlGEYHt/rssjL1KJBC4/rAPhzbu2Bfv7553XPPffo7rvv1mOPPabmzZ0/KtqzZ4/S0tJ04MABLVy4UL/88ss5H1e3bl3POwYAwE9YLBaFhRVKMlR+iDYUFlZYdioHAP/gVoCOiYmRJH322WdaunTpOdfOfBdd3hnRpaWl7twOAAC/1atXB6WlZf66B/pcVusqJSZ29EJXADzh9ihvvlsGAODPDR8+SjNm9JHDYUg68yChIat1lWJipiolZbGXOwTgKrcC9MSJE01uAwCAwHPihDRgQIQaNlys7t2naM2al1RSUlNhYSeUmNhBKSmLFeHC9JTTjtP6/PDnat2odSV2DeDPeDSJEAAAlM8wpCFDpJ07pY0bI3TDDRN/fd1w66e4pY5SDV4yWBm7M7T/yf2qW4PnigBv8ShAb9q0STk5OcrPz5fD4TjnmsVi0bhx4zxqzhclJSUpMjJSNptNNpvN2+0AAHzUCy9ICxZI77wj3XDDb6+7E54dhkPDMoZpwecLtKDPAsIz8Dt2u112u135+flVcj+L4cbZOb/88ot69uypbdu2lX0nfabMmf9tsVgC6qHBnJwcxcXFKTs7W7Gxsd5uBwDgw5Yske66S5owQfJ016NhGHpsxWN69b+vKv2udA1oNcCUHoFAVFV5za1zoP/+979rx44dmj9/vr755hsZhqHMzEx99dVXeuSRR9S6dWsdPHjQ7F4BAPB5n30mDRgg9ekjjR/vWS3DMPRU5lOa8d8ZeqPXG4RnwEe4FaBXrFihhx9+WPfee2/Zww9Wq1XNmzdXWlqaoqKi9OSTT5raKAAAvu6nn6TERKl5c+mttySrW3/LOhmGoWfWPqNpWdOU1iNNQ2KHmNcoAI+49Uf76NGjuv766yVJtWvXliQdP3687Hp8fPwfTioEACDQFBdLfftKhYXS0qVSrVqe1cvcm6kXN72ol+Jf0vA2w81pEoAp3ArQTZo00Y8//ihJqlatmho0aKBPP/207Hpubi7nRAMAgoZhSI8/Lm3ZIr33nnT55Z7X7HZlN60dtFZJ7ZM8LwbAVG6dwtGpUyetXr1aycnJkqR7771X//rXvxQSEiKHw6Fp06apW7dupjYKAICvSkuTXn9devNNqUMHc2paLBZ1ju5sTjEApnIrQD/11FNavXq1ioqKVK1aNU2cOFFffPFF2bF1nTp10v/93/+Z2igAAL5o7VrpySed/zz0kLe7AVAV3ArQLVu2VMuWLcv+/eKLL9aaNWt09OhRhYSEuDRVCQAAf7Vnj9Svn9Sli/Tvf3u7GwBVxdRJhHXq1DGzHAAAPis/X+rdW6pf3zkwJZTZvkDQqNAf9/T0dLeKDxo0yK2PAwDAl5WWSvffLx08KGVlSZ6sH2V9n6W4JnEKtZLAAX9RoT+tDzzwgMuFLRYLARoAEJDGjJFWrZJWrJCuucb9Ohm7M3T323drarepGtF2hHkNAqhUFQrQ+/btq+w+AADwC+npzv3OL70keXLg1Ad7P1Dfd/oq8ZpEPRz3sHkNAqh0FQrQl5txoCUAAH5u61Zp6FDpwQedp264a/3+9eq9oLe6XtFV9j52hYWEmdckgErn0YaroqIi5eTk6PDhw+rQoYPq1atnVl8AAPiU77+X7rxTuvFGacYMyd15YZu+3aQ75t+hWy67RYvuWaTwkHBzGwVQ6dyaRChJ06dPV+PGjdWxY0fdfffd2rFjhyQpLy9P9erV03/+8x/TmgQAwJtOnHCeuBEeLr37rlStmnt1Psn9RN3nddeNTW7UkvuWqHpodXMbBVAl3ArQs2bN0pNPPqmEhAS9+eabMgyj7Fq9evXUuXNnLViwwLQmAQDwFsNwDkjZtUtatkxq2NC9OiWlJbpv8X1q0aCF3r//fdUMq2luowCqjFtbOKZMmaLevXtr/vz5+vnnn8+7HhcXp+nTp3vcHAAA3vbCC9LChdI770itW7tfJywkTEvvW6pmFzVT7fDa5jUIoMq5tQL99ddfq3v37he8Xrdu3XKDNQAA/mTJEmnsWGnCBKlvX8/rtWjQQpHVIz0vBMCr3ArQderUUV5e3gWvf/nll2rUqJHbTQEA4G07dkgDBkh9+kjjx3u7GwC+xK0A3aNHD73++us6evToede++OILvfHGG0pMTPS4OQAAvOGnn6TEROmqq6S33pKsbj9yDyAQufUlISUlRaWlpWrRooXGjh0ri8Wit956SwMGDNCNN96oBg0aaDzfrgMA/FBxsXPV+eRJaelSqVYtb3cEwNe4FaCbNGmi7OxsJSQkaOHChTIMQ3PmzFFGRoZsNpu2bt3KmdAAAL9jGNKIEc6BKe++K112mes1Dh0/pC8Of2F+cwB8htuDVBo0aKCZM2dq5syZ+umnn+RwOFS/fn1Z+TkXAMBPpaVJb7whvfmm1KGD6x+fdyJPXeZ0kUUWbX9ku6wW/k4EApFHkwjPqF+/vhllAADwmrVrneO5n3zSee6zq46cPKL4OfE6XHhYHz3wEeEZCGCmBGgAAPzZnj1Sv35Sly7Sv//t+scfKzqmhHkJOpB/QOsHr9e19a41v0kAPoMADQAIavn5zhM36teXFiyQQl38m/F48XH1mNdDu/N2a93gdWrZsGXlNArAZxCgAQBBq7RUstmkH36QsrKkOnVc+/iTJSeVaE/Up4c+1eqBqxXbOLZyGgXgUwjQAICg9Y9/SJmZ0sqV0jXXuP7x07ZOU1Zullb1X6V2l7Yzv0EAPokADQAISunp0uTJ0tSpUny8ezVG3TxK3a/qrtaNWpvbHACfxiPCAICgs3WrNHSo9OCD0hNPuF8nLCSM8AwEIQI0ACCofP+9dOed0o03SjNmSBaLtzsC4G8I0ACAoHHihNS7txQe7pw0WK2atzsC4I/YAw0ACAqG4RyQsmuXtGmT1LChtzsC4K9YgQYABIV//lNauND58GBrF7YtG4ahlXtWyjCMymsOgF8hQAMAAt5770njxkkTJ0p9+lT84wzDUPK6ZPWY30Nbvt9Saf0B8C8EaABAQDqzYrxjhzRwoDM4jxvnWo2UDSn6fxv/n16Kf0k3N7u5EroE4I/YAw0ACBgFBQVKTp6sjIxNKimpJau1UEeOdNAVV4zSW29FyOrCstG/Nv1L49eP1wudX1BS+6TKaxqA3yFAAwACQkFBgdq376OdO5+SwzFRkkWSISlTRUV95HAslhRRoVrTs6Zr9JrRGt9pvMbcMqbymgbgl4J6C8err76quLg4hYeH67nnnvN2OwAADyQnT/41PCfIGZ71668J+vrrJI0dO6VCdV7772t6YtUTevrmpzXxrxMrqVsA/iyoA3STJk307LPPqm/fvt5uBQDgoYyMTXI4upV7zeFI0LJlm/60xsGCg3oy80mNbDtSL3Z5URamrAAoR1Bv4UhMTJQkLV++3MudAADcVVIiffKJoZ9/rqXfVp5/z6KSkpoyDOMPQ3GTiCb6ZOgnur7+9YRnABfkkyvQhYWFmjBhgrp3765LLrlEVqtV6enp5b63uLhYo0ePVtOmTVWzZk21a9dOa9asqeKOAQBVxeGQtm+XpkyRevaU6taVOnSw6PjxQjn3PJfHUFhYYYVCcYsGLQjPAP6QTwbovLw8Pf/889q1a5dat279h1/IBg8erGnTpmngwIGaPn26QkND1aNHD23evLkKOwYAVBbDcE4PfOUVqW9fqX596YYbpLFjpeJi6ZlnpKwsafjwDrJaM8utYbWuUmJixyruHECg8sktHE2aNNGPP/6oBg0aKDs7W23atCn3fdu2bdPChQs1ZcoUJSU5jxgaOHCgWrRooaefflobN24se+/8+fP18MMPy2KxaMCAAXrllVeq5HMBALjuwAFp3brf/jl4UAoNlW66SRoxQurcWWrXTqpW7bePiYkZpfXr+2jnTuOsBwkNWa2rFBMzVSkpi7316QAIMD4ZoMPCwtSgQYM/fd+iRYsUGhqqoUOHlr1WrVo1DRkyRMnJycrNzVXTpk0lSffff7/uv//+SusZAOC+Q4ekDz/8LTDv3StZLFJsrNS/vzMwd+wo1a594RoRERHasmWxxo6domXLXlJJSU2FhZ1QYmIHpaQsVkRExY6wA4A/45MBuqK2b9+uq6++WrV/9xW1bdu2ZdfPBOjylJaWqqSkpOzXoqIihYWFyerKSfsAAJcdPSp99JEzLK9dK33xhfP1666Tund3BuZbb3Xub3ZFRESEUlMnKjVVf/jAYNb3Wbr0okvV9KIL/x0BABfi1wH6hx9+UOPGjc97vXHjxjIMQwcPHvzDj09JSdGzzz5b9gX2hRde0KxZszRo0KBK6RcAglVhobRpkzMsr1sn5eQ4HwaMjnaG5WeekW67TSrnS7rbLhSet+VuU9c5XXV3zN2afeds824IIGj4dYA+efKkqp29Ae5X1atXL7v+RyZMmKAJEyZUSm8AEMyKi6WtW3/bkrF1q/O4uUaNnIH5kUecv0ZHV21f23/crm5zu6llw5Z6ucfLVXtzAAHDrwN0jRo1VFRUdN7rp06dKrtutqSkJEVGRp7zms1mk81mM/1eAOAvSkudq8pnAvPHH0snT0oXXyz99a/S1KnOwHzttc69zd7w+eHP1SW9i5rXba4V969Q7fA/2FANwOfZ7XbZ7fZzXsvPz6+Se/t1gG7cuHG52zR++OEHSc7TPMw2depUxcbGml4XAPyJYTj3LZ8JzOvXS/n5Uq1aUqdO0nPPOQPzX/4ihYR4u1tpd95udUnvoksvulSZAzIVWT3yzz8IgE8rbwEzJydHcXFxlX5vvw7QrVu31vr163X8+PFzHiTcunWrLBaLWrdu7cXuAMD3/dlkvt/eJ33zzblHyx0+LIWHSzffLP3tb87A3KaN8zVfsveXveqc3lmX1LxEqweuVt0aLj6ZCAC/49fHTfTt21enT5/W66+/XvZacXGxZs+erXbt2v3hCRwAEKwKCgo0cuQERUd3UbNmdyo6uotGjpyggoKCc96XmyvNnSs9+KAUFSU1b+7cu7xvnzRkiLR6tfM0jQ8/lMaNkzp08L3wLEnPrHtGtcNra+2gtapfq7632wEQAHx2BTotLU1Hjx5Vbm6uJGnZsmX67rvvJEkjR45URESE2rZtq379+mnMmDE6dOiQmjdvrtmzZ+vAgQOaNWuWN9sHAJ9UUFCg9u37aOfOp+RwTNSZYSNpaZlavbqPkpMXa/PmCK1bJ+3e7fyYVq2ku++Wbr9duuUWKdLPdj/M7DVTx4qOqVHtRt5uBUCAsBiGYXi7ifJER0fr22+/Lffavn37dNlll0lyrjiPGzdOc+fO1ZEjR9SqVSulpKSoS5cupvZzZk9NdnY2e6AB+K2RIycoLa39r5P6fm+lpCxdddVE3X67c0vGX//qHJ0NAP6gqvKaz65A79u3r0LvCw8P16RJkzRp0qRK7ggA/F9GxqZfV57Lk6BmzV7SV19VZUcA4H98NkD7qjPH2HF0HQB/YxiGSkpqybltozwWORw1K/xgIQD4ijNH2nGMnY/iGDsA/spiscgwCiUZKj9EGwoLKyQ8A/A7ZxY2q+oYO78+hQMAUHEffyzl5XWQlFnudat1lRITO1ZtUyY5XnxcS3ct9XYbAIIEARoAgsDSpVJ8vNS27Shde+1LslpXyrkSLUmGrNaViomZqpSUv3mzTbecLDmpRHuiBi0ZpJ8Kf/J2OwCCAAEaAALcm286j6Hr2VNavTpC27Yt1ogRWYqKilfTpr0VFRWvESOytGXLYkVERHi7XZcUnS7SXQvvUlZult63vc85zwCqBHugASBAGYb04ovSM884B6C8/LJzrHb16hFKTZ2o1NSKTyL0RcWlxer3Tj99dOAjLb9/uW65/BZvtwQgSLACDQAByOGQnnzSGZ4nTpReecUZnn/PX8Pzacdp9X+3vzL3Zuq9e99T5+jO3m4JQBBhBRoAAkxxsTR4sLRwoTM4P/qotzsyV6mjVIOXDNaSXUu0qN8iJTQvbygMAFQeArSLOAcagC8rKJD69JE++kh6+22pb19vd2S+DQc26O0v3tb8u+er97W9vd0OAB9Q1edA++wob1/DKG8Avu6nn6QePaTdu6UlS5yjuAPVN0e+0RUXX+HtNgD4mKAf5Q0AqLj9+6Vu3aT8fOfq8w03eLujykV4BuBNPEQIAH5uxw7p5pul06elTZsCPzwDgLcRoAHAj338sdSpk9SwoTM8X3mltzsCgMBHgAYAP3VmumBsrHPbRqNG3u4IAIIDARoA/NDMmb9NF1yxQrroIm93ZL63v3hbx4qOebsNADgPARoA/IhhSC+8IA0dKg0b5jzruXp1b3dlvlf/+6ruXXSv5nw6x9utAMB5CNAA4CfOTBdMTv7j6YL+btb/ZunR5Y9qZNuRGt5muLfbAYDzcIydixikAsAbzp4uOGOG9Mgj3u6ocsz/bL6GLBuih+Me1rSEaX47ahxA1WKQio9ikAoAbzl7uuC8eYE5XVCSFn+5WPcuulcD/zJQbya+KauFH5ICcA2DVAAA50wXXLVKuu02b3dUOTJ2Z+i+xfep3/X9NLPXTMIzAJ9GgAYAH7V/v/OYumPHAnu64GnHaY1ZO0a9ru6l9DvTFWINwI3dAAIKARoAfNCOHVJCglSjRuAPSAm1hmrtoLW6uMbFCgsJ83Y7APCn+BkZAPiYs6cLbt4c2OH5jIa1Gyo8JNzbbQBAhRCgAcCHLF0qde3623TBhg293REA4PcI0ADgI85MF+zVK3CnCwJAICBAA4CXnT1d8OGHpQULAnO6IAAECgI0AHiRwyE98cRv0wXT0gJzuqAkffXzV/pg7wfebgMAPMYpHADgJcEyXVCSvjnyjTq/1VkNajXQ7dG3c1QdAL9GgHYRo7wBmKGgwLnfecMG6Z13nJMGA9W3+d+q81udVSu8llb0X0F4BmA6Rnn7KEZ5AzDL2dMFly4N3OmCknSw4KA6zeokh+HQhgc36NKLLvV2SwACGKO8ASAABct0QUk6dPyQbk+/XcWlxYRnAAGFhwgBoIrs2CHdfLPzwcFNmwI7POedyFOXOV2Ufypf6wavU1SdKG+3BACmIUADQBXYsME5XbBRo8AfzS1JCz5foMOFh7Vu8Do1r9vc2+0AgKkI0ABQyZYscW7biI2V1q8PjumCj7V5TDse2aFr613r7VYAwHQEaACoRDNnOk/YCLbpghaLRQ1rB8F3CgCCEgEaACqBYUj//CfTBQEgEBGgAcBkZ6YLjh0b+NMFASAYcYwdAJgomKYLAkCwYgUaAExSUCD17Cm9+65zumCgh+fTjtNK3ZqqktISb7cCAFWKAA0AJjh8WOrcWcrKklatCuzR3JJU6ijV4CWDNWr1KH1y8BNvtwMAVYotHADgoX37pG7dgmO6oCQ5DIeGZQzTgs8XaEGfBbq52c3ebgkAqhQB2kVJSUmKjIyUzWaTzWbzdjsAvGzHDikhQapZMzgGpBiGoRErRmjW9llKvytd/a7v5+2WAEB2u112u135+flVcj+LYRhGldzJz+Xk5CguLk7Z2dmKjY31djsAfMCGDVJionTFFdLKlYE/IMUwDD2V+ZSmZU3TzF4zNSR2iLdbAoBzVFVeYw80ALjhzHTBuLjgmC5oGIaeWfuMpmVNU1qPNMIzgKBGgAYAFwXjdMFDhYc0838z9VL8SxreZri32wEAr2IP9P9v7+6Do6rvPY5/No8YWYhCbQM2gJOhOEXMEC9SKLaDQR40oWoiRI2xjUBVJlx6OzIgFNColOnlIhVUNANJI4s83CJYEY0U0UBASIH2Kog0gIUQmjaBPD+e+0ckGgmQQzb724f3a2bH4exmz+cI/PLJ4ez5AkAHWZb0/PMtA1KeeEJatixwBqR8r/v39NmTn6l3RG/TUQDAOM5AA0AHNDdLGRkt5XnhQumllwKnPF9AeQaAFpyBBoArqKtrmS64bh3TBQEAFGgAuKyKCum++1ruuLF+vf8PSAEAXBkFGgAu4exZacIE6ehRads26ac/NZ0IAOANuAYaANpRVCT9+MfSP/7RMu4T4IIAABSCSURBVF0wUMrzW4ffUmFxoekYAODVKNAA8C2HDkkjRrR8cDA/X4qNNZ3IM97+/G0lr0/Wq/teNR0FALwaBRoAvmHnTumOO6SoqMAYzX3Be8fe0/3r7lfCDxL00oSXTMcBAK9GgQaArwTadMELdhzfoYlrJ2rMTWPkut+l0OBQ05EAwKtRoAFAgTldUJLyT+brnjX3aFT0KG14YIPCgsNMRwIAr0eBBhDQLEt67jlpypSW+zuvXSuFh5tO5RmfnPpE498Yr9v63KZNkzepW0g305EAwCdwGzubZs6cqZ49eyolJUUpKSmm4wDohOZmacaMlqmCCxdK8+ZJDofpVJ6T9ZcsDb5hsN5+8G1FhEaYjgMAV83lcsnlcuncuXMe2Z/DsizLI3vycYWFhYqLi9P+/fs1dOhQ03EAdNKF6YLr10srVkjTpplO5HlNzU2qbqiWM9xpOgoAuIWn+hpnoAEEnG9PF7zvPtOJzAgOCqY8A8BVoEADCChMFwQAdBYFGkDAKCpquU1dRUXLdMFAGZACAHAv7sIBICAcPNgyXdCypF27KM8AgKtHgQbg9z78sO10wZtuMp3Ic0oqS7Q4f7H4vDgAuA8FGoBf++MfpbFjpdtuC6zpgpJUWl2q+D/Ea2nBUp2tOms6DgD4DQo0AL/12mtSUpKUmBhY0wUlqaymTHf94S6drTqr7Wnb9d3uAfSTAwB0MQo0AL9jWVJmpjR1ast0QZcrcKYLStL5uvMa98Y4nTh3QnmpeRrUe5DpSADgV7gLBwC/8s3pgs88I82dG1jTBSvrKzXhjQk6UnpE29O265bv3mI6EgD4HQo0AL9RVyc98oi0YYP0yiuBN12wpqFGia5EHSw5qPdT39fQKKamAkBXoEAD8AsVFdK990offRS40wUPlRzSoZJDeufBdzT8xuGm4wCA36JAA/B5TBdscfuNt+v4fx5X97DupqMAgF+jQAPwaUwXbIvyDABdj7twAPBZTBcEAJhAgQbgkwJ5uiAAwCwKNACfE8jTBQEA5lGgAfiUQJ4uKEmWZWnRx4t0ovyE6SgAELAo0AB8wjenCz7+eOBNF5RayvPT25/W7A9m68/H/2w6DgAELO7CYdPMmTPVs2dPpaSkKCUlxXQcICA0N0sZGdLy5YE5XfCCzJ2ZeuHjF7TkriV6NPZR03EAwGu4XC65XC6dO3fOI/tzWJZleWRPPq6wsFBxcXHav3+/hg5luhfgKd+cLvjyyy1noAPR4vzFmpU3S8+Pfl6zR802HQcAvJKn+hpnoAF4rQvTBT/+OHCnC0rSsj3LNCtvln5zx28ozwDgBSjQALyGZVlyfHVtxtmz0vjx0hdfSO++G7jTBV/d96pmvDtDT414Sgt+usB0HACAKNAADKuoqNDTT/9OW7bkq6HhWoWGVuknPxmpjz76taqrndq5U7r1VtMpzWhqbtLa/1urjGEZWhS/qPWHCwCAWRRoAMZUVFToRz+6X5999is1Ny+Q5JBk6fjxbQoNvV/79m3UkCFOwynNCQ4K1taHtio8OJzyDABehNvYATDm6ad/91V5HqeW8qyv/jtOTU0zlZX13wbTeYduId0ozwDgZSjQAIzZsiVfzc1j232uuXmcNm/O93AiAACujAINwAjLstTQcK2+PvP8bQ41NESIO20CALwNBRqAEQ6HQ6GhVZIuVZAthYZWcfkCAMDrUKABGJOQMFJBQdvafS4o6F0lJv7Yw4nMyD+Zr5c/edl0DABAB1GgARjz3HO/1s03L1FQ0FZ9fSbaUlDQVt188/8oM/O/TMbziL2n9mr8G+O17tN1amxuNB0HANABFGgAxjidTu3evVHTp+9R//53qW/fierf/y5Nn75Hu3dvlNPp37ewO3DmgMbmjtUt371FW1K2KCSIO4sCgC9gtQZglNPp1IsvLtCLL7adROjv/nb2b4rPiVfM9TF658F31D2su+lIAIAO4gw0AK8RKOX5SOkRxefE68YeN2rbw9vUs1tP05EAADZQoAHAg479+5hG54xWr4heej/1fV1/zfWmIwEAbKJAA4AHfXL6E/UI76EPHvlA37n2O6bjAACuAtdAA4AHTR48WfcOulfhIeGmowAArhJnoAHAwyjPAODbKNAAAACADRRoAAAAwAYKNAAAAGADBRoA3KyyvlLT35muspoy01EAAF2AAg0AblTTUKNEV6KyD2br72V/Nx0HANAFuI0dALhJXWOd7n3zXu05tUfvPvSu4vrEmY4EAOgCFGgAcIP6pnolr0/Whyc+1J8e/JNG9RtlOhIAoItQoAGgkxqbG/XQ/z6kbce26a3Jb2n0gNGmIwEAuhAFGgA6oam5SWmb0rTp8CZtSN6gcTHjTEcCAHQxCrRNM2fOVM+ePZWSkqKUlBTTcQAYVl5brkMlh7TmvjWaOGii6TgAEJBcLpdcLpfOnTvnkf05LMuyPLInH1dYWKi4uDjt379fQ4cONR0HgBdpbG5USBDnIwDANE/1NW5jBwCdRHkGgMBCgQYAAABsoEADAAAANlCgAQAAABso0ADQAb/f83ttOrzJdAwAgBegQAPAFbyy7xVlvJuhvaf2mo4CAPACFGgAuIxVf1mlx//0uDKGZei50c+ZjgMA8AIUaAC4hDV/XaP0zemaFjdNS8ctlcPhMB0JAOAFKNAA0I6Nn27UI398RGmxaVpx9wrKMwCgFQUaAL5ly5EtmrxxspJ/mKzXE15XkIOlEgDwNb4rAMC3HCw5qISBCcr5WY6Cg4JNxwEAeBnmzwLAt8y9Y66ampsozwCAdnEGGgDaQXkGAFwKBRoAAACwgQINAAAA2ECBBgAAAGygQAMISEdKjyj9rXTVNdaZjgIA8DHchQNAwDn272ManTNakd0iVVFfofCQcNORAAA+hDPQAALKifITGp0zWteGXqu81Dz1juhtOhIAwMdQoAEEjFPnT+nOnDsV7AjW9rTtinJGmY4EAPBBXMIBICCUVJbozpw7Vd9Ur50/36kbe9xoOhIAwEdRoAH4vdLqUsX/IV7n685r5893qn9kf9ORAAA+jEs4APi90urS1ss2Yq6PMR0HAODjOAMNwO8N6j1IhdMKFeTgnAEAoPP4bgIgIFCeAQDuwncUAAAAwAYKNAAAAGADBRoAAACwgQINwC80Njfqsc2PqeAfBaajAAD8HAUagM9ram5S2qY0ZR/M1tmqs6bjAAD8HLexA+DTmq1mTd0yVWv/tlZr71+rxB8kmo4EAPBzFGgAPsuyLE1/Z7pWHVilnHtzlPzDZNORAAABgAINwCdZlqVfbfuVXt73sl5PeF0PD3nYdCQAQIDgGmgAPseyLM35YI6W7lmq5ROWK31ouulIAIAAQoEG4HMsWSquLNaSu5boif94wnQcAECA4RIOAD4nyBGkVRNXyeFwmI4CAAhAnIEG4JMozwAAUwK2QNfX1ys9PV39+vVTZGSkRowYoYICBjAAAADg8gK2QDc2NmrAgAHatWuXysvLNWPGDCUkJKi6utp0NAAAAHixgC3QERERmjt3rvr27StJmjRpksLCwnTkyBHDyQAAAODNvLJAV1VVaf78+Ro/frx69eqloKAg5eTktPva+vp6zZo1S3379lVERISGDx+uvLw82/s8evSoysrKFBMT09n4ANxk46cblbE1Q81Ws+koAAC08soCXVpaqmeffVaHDx9WbGzsZT8slJaWpqVLlyo1NVXLli1TSEiIJkyYoF27dnV4f7W1tUpNTdWcOXPkdDrdcQgAOmnLkS2avHGy/ln9T1mWZTqO33K5XKYjAIDP8coC3adPH505c0ZFRUVavHjxJb957t27V2+++aYWLVqkRYsW6bHHHtMHH3ygfv366amnnmrz2jVr1sjpdKpHjx564omv7xvb2NiopKQkDRw4UHPnzu3S4wLQMe8de09J65OU+INE5fwsR8FBwaYj+S0KNADY55UFOjQ0VDfccMMVX7dhwwaFhIRoypQprdvCw8OVnp6u3bt369SpU63bH3zwQVVUVOj8+fNasWKFpJZpZqmpqQoJCVF2drb7DwSAbTuO79DEtRM15qYxct3vUmhwqOlIAAC04ZUFuqMOHDiggQMHqnv37m22Dxs2rPX5y5k6darOnDmjdevWcU9ZL+XPZ8e8+dhMZVvw+wW6Z809GhU9Shse2KCw4DC3vn9XHJe73rOz7+PNf57Qlj//XnnzsZnK1tX7ZV0zw6cLdHFxsaKioi7aHhUVJcuydPr06Ut+7cmTJ5WVlaW9e/eqV69erZd35Ofnd2Vk2OTNf3k6y5uPzUS2vaf2KnNFpm7rc5s2Td6kbiHd3L4PvtHAG/jz75U3HxsF2vPv6c/rmk+P8q6pqVF4ePhF27t169b6/KVER0erubnjn+y/8F6fffaZzZTojHPnzqmwsNB0jC7hzcdmItvB0wfVw+qhzJszdfivh7tkH11xXO56z86+z9V+vTf/OfRX/vz/3JuPzVS2rt4v61pbF3ra5TqgW1hebt++fZbD4bCys7Mvem7w4MFWfHz8Rds//fRTy+FwWCtXrnRbjtzcXEsSDx48ePDgwYMHDy9/5Obmuq0Dtsenz0BHRUW1e5lGcXGxpJa7ebjL2LFjlZubq/79++uaa65x2/sCAADAPWpqanT8+HGNHTu2S/fj0wU6NjZWO3bsUGVlZZsPEhYUFMjhcCg2NtZt++rdu7ceeught70fAAAA3G/kyJFdvg+f/hBhUlKSGhsbtXLlytZt9fX1Wr16tYYPH946phsAAABwF689A718+XKVl5e33st58+bN+vLLLyVJGRkZcjqdGjZsmJKTkzV79myVlJQoJiZGq1ev1okTJ7Rq1SqT8QEAAOCnHJblnTNyBwwYoJMnT7b7XFFRkaKjoyW1nHGeN2+ecnNzVVZWpiFDhigzM1Px8fGejAsAAIAA4bWXcBQVFampqandx4XyLElhYWH67W9/q1OnTqm6uloFBQVGy/Mrr7yiuLg4hYWF6ZlnnjGWAwDcob6+Xunp6erXr58iIyM1YsQIFRQUmI4FAJ0ybdo0RUVFKTIyUrfeeqvefvttW1/vtQXaV/Xp00cLFy5UUlKS6SgA0GmNjY0aMGCAdu3apfLycs2YMUMJCQmqrq42HQ0ArtrMmTNVVFSk8vJyZWVl6eGHH1ZZWVmHv54C7WaJiYm655571LNnT9NRAKDTIiIiNHfu3NYPZU+aNElhYWE6cuSI4WQAcPUGDRrUOnhPkhoaGlo/d9cRAVugq6qqNH/+fI0fP169evVSUFCQcnJy2n1tfX29Zs2apb59+yoiIkLDhw9XXl6ehxMDwOV5Yl07evSoysrKFBMT4+74AHCRrlzXnnzySUVEROj222/XnXfeqcGDB3c4V8AW6NLSUj377LM6fPiwYmNj5XA4LvnatLQ0LV26VKmpqVq2bJlCQkI0YcIE7dq1y4OJAeDyunpdq62tVWpqqubMmSOn09kVhwAAbXTlurZ8+XJVVVUpLy9PY8aMsResS+ccerH6+nqrpKTEsqzLjwvfs2eP5XA4rCVLlrRuq62ttWJiYqyRI0de8v1/+ctfWgsXLnR/cAC4hK5c1xoaGqy7777bSk1N7ZrwANCOru5rFyQkJFhbt27tcK6APQMdGhqqG2644Yqv27Bhg0JCQjRlypTWbeHh4UpPT9fu3bttXS8DAF2pq9Y1y7KUmpqqkJAQZWdnuz03AFyKp/paY2Ojvvjiiw7nCtgC3VEHDhzQwIED24wKl6Rhw4a1Pv9NTU1Nqq2tVVNTkxoaGlRXV6fm5maP5QWAK7G7rk2dOlVnzpzRunXrLvvPpwBgip117fz583K5XKqqqlJTU5PWr1+vHTt26I477ujw/ijQV1BcXKyoqKiLtkdFRcmyLJ0+fbrN9szMTEVERCgrK0vPP/+8IiIilJub66m4AHBFdta1kydPKisrS3v37lWvXr3kdDrVo0cP5efnezIyAFyWnXXN4XDotdde0/e//3317t1bixcvlsvl0pAhQzq8P68d5e0tampqFB4eftH2C7c+qampabN9/vz5mj9/vkeyAcDVsLOuRUdH869oALyenXXN6XRq+/btndofZ6Cv4JprrlFdXd1F22tra1ufBwBfwroGwN94el2jQF9BVFSUiouLL9p+YVufPn08HQkAOoV1DYC/8fS6RoG+gtjYWH3++eeqrKxss72goEAOh0OxsbGGkgHA1WFdA+BvPL2uUaCvICkpSY2NjVq5cmXrtvr6eq1evVrDhw9vHW8LAL6CdQ2Av/H0uhbQHyJcvny5ysvLW+8NuHnzZn355ZeSpIyMDDmdTg0bNkzJycmaPXu2SkpKFBMTo9WrV+vEiRNatWqVyfgAcBHWNQD+xhvXNYdlWZbb39VHDBgwQCdPnmz3uaKiIkVHR0tq+Qlm3rx5ys3NVVlZmYYMGaLMzEzFx8d7Mi4AXBHrGgB/443rWkAXaAAAAMAuroEGAAAAbKBAAwAAADZQoAEAAAAbKNAAAACADRRoAAAAwAYKNAAAAGADBRoAAACwgQINAAAA2ECBBgAAAGygQAMAAAA2UKABAAAAGyjQAAAAgA0UaAAAAMAGCjQAAABgAwUaAAAAsIECDQABZMGCBQoKCtKxY8f06KOP6rrrrlNkZKR+8YtfqLa21nQ8APAJFGgACCAOh0OS9MADD6iqqkqLFi3SpEmTlJ2drYULFxpOBwC+IcR0AACA58XFxWnlypWtvy4tLVVWVpZeeOEFg6kAwDdwBhoAAozD4dC0adPabBs1apT+9a9/qbKy0lAqAPAdFGgACEDR0dFtfn3ddddJksrKykzEAQCfQoEGgAAUHBzc7nbLsjycBAB8DwUaAAAAsIECDQAAANhAgQYAAABsoEADAAAANjgsPjECAAAAdBhnoAEAAAAbKNAAAACADRRoAAAAwAYKNAAAAGADBRoAAACwgQINAAAA2ECBBgAAAGygQAMAAAA2UKABAAAAGyjQAAAAgA0UaAAAAMAGCjQAAABgw/8DhC2uUmiWeNMAAAAASUVORK5CYII=",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x320026390>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "elapsed time: 1.404498596 seconds\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(PyObject <matplotlib.text.Text object at 0x31e6b7390>,PyObject <matplotlib.text.Text object at 0x31feff2d0>)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n_ = collect(50:50:500);\n",
    "time_ = zeros(size(n_));\n",
    "for k = 1:length(n_)\n",
    "    n = n_[k];\n",
    "    A = rand(1200,n);\n",
    "    Q = zeros(1200,n);  R = zeros(600,600); \n",
    "    \n",
    "    tic();\n",
    "    R[1,1] = norm(A[:,1]);\n",
    "    Q[:,1] = A[:,1]/R[1,1];\n",
    "    for j = 2:n\n",
    "        R[1:j-1,j] = Q[:,1:j-1]'*A[:,j];\n",
    "        v = A[:,j] - Q[:,1:j-1]*R[1:j-1,j];\n",
    "        R[j,j] = norm(v);\n",
    "        Q[:,j] = v/R[j,j];\n",
    "    end\n",
    "    time_[k] = toc();\n",
    "end\n",
    "\n",
    "using PyPlot\n",
    "loglog(n_,time_,\"-o\",n_,(n_/500).^2,\"--\")\n",
    "xlabel(\"n\"), ylabel(\"elapsed time\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Flops aren't everything. On massively parallel computers, for example, the time needed for memory access and communication often competes with or dwarfs the floating point arithmetic. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.4.3",
   "language": "julia",
   "name": "julia-0.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}

TB-Lecture-08-matlab.ipynb