sparse-direct.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra\n",
    "using SparseArrays\n",
    "using PyPlot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8×8 SparseMatrixCSC{Float64,Int64} with 22 stored entries:\n",
       "  [1, 1]  =  1.0\n",
       "  [2, 1]  =  0.01\n",
       "  [3, 1]  =  0.01\n",
       "  [4, 1]  =  0.01\n",
       "  [5, 1]  =  0.01\n",
       "  [6, 1]  =  0.01\n",
       "  [7, 1]  =  0.01\n",
       "  [8, 1]  =  0.01\n",
       "  [1, 2]  =  0.01\n",
       "  [2, 2]  =  1.0\n",
       "  [1, 3]  =  0.01\n",
       "  [3, 3]  =  1.0\n",
       "  [1, 4]  =  0.01\n",
       "  [4, 4]  =  1.0\n",
       "  [1, 5]  =  0.01\n",
       "  [5, 5]  =  1.0\n",
       "  [1, 6]  =  0.01\n",
       "  [6, 6]  =  1.0\n",
       "  [1, 7]  =  0.01\n",
       "  [7, 7]  =  1.0\n",
       "  [1, 8]  =  0.01\n",
       "  [8, 8]  =  1.0"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# A has diagonal + nonzero first row and column\n",
    "n = 8\n",
    "Iv, Jv, Vv = Int64[], Int64[], Float64[]\n",
    "for ind = 2:n\n",
    "    push!(Iv, 1, ind)\n",
    "    push!(Jv, ind, 1)\n",
    "    push!(Vv, 0.01, 0.01)\n",
    "end\n",
    "for ind = 1:n\n",
    "    push!(Iv, ind)\n",
    "    push!(Jv, ind)\n",
    "    push!(Vv, 1)\n",
    "end\n",
    "A = sparse(Iv, Jv, Vv, n, n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAGiCAYAAADa2tCeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAYQElEQVR4nO3df2wWhR3H8c9B7YNin0dBijQ8QMPYAPkhUsfKjzkFuzSEyJYxXZSVsS1pUvlhY+bQP8TF8LAs+2OLs1lx6UYIqVkUxG2AJZOiYWwFJWOdQRjEVgUZRp4r/eOI7e0vGrtC+9zTfp/rPbxfyWU+tzvvm0d83t7d057j+74vAACG2IiwBwAA5CcCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMJE3gXnxxRdVWlqqUaNGaf78+XrrrbfCHmlAhw4d0ooVK1RSUiLHcbR79+6wR8pIKpXSvffeq6KiIhUXF2vlypU6efJk2GMNqK6uTnPmzFE8Hlc8Hld5ebn27t0b9liBpVIpOY6jjRs3hj1KvzZv3izHcXotd955Z9hjZeSjjz7SY489prFjx+qWW27R3XffrWPHjoU91oCmTJnS5z13HEc1NTWhzJMXgXn55Ze1ceNGPfPMM3r33Xe1ZMkSVVZWqq2tLezR+tXZ2am5c+fqhRdeCHuUQJqbm1VTU6MjR46oqalJn3/+uSoqKtTZ2Rn2aP2aOHGitm7dqqNHj+ro0aN64IEH9NBDD6m1tTXs0TLW0tKi+vp6zZkzJ+xRMnLXXXfp3LlzPcuJEyfCHmlAn332mRYtWqSbbrpJe/fu1b///W/98pe/1G233Rb2aANqaWnp9X43NTVJklatWhXOQH4e+OpXv+pXV1f3Wjd9+nT/pz/9aUgTBSfJ37VrV9hjZOXChQu+JL+5uTnsUQK7/fbb/ZdeeinsMTLS0dHhT5s2zW9qavLvu+8+f8OGDWGP1K9nn33Wnzt3bthjBPbUU0/5ixcvDnuMIbFhwwZ/6tSpfnd3dyjHj/wZzJUrV3Ts2DFVVFT0Wl9RUaHDhw+HNNWNJZ1OS5LGjBkT8iSZ6+rqUmNjozo7O1VeXh72OBmpqanR8uXLtWzZsrBHydipU6dUUlKi0tJSPfLIIzpz5kzYIw1oz549Kisr06pVq1RcXKx58+Zp27ZtYY8V2JUrV7Rjxw6tXbtWjuOEMkPkA3Px4kV1dXVp/PjxvdaPHz9e58+fD2mqG4fv+6qtrdXixYs1a9assMcZ0IkTJ3TrrbcqFoupurpau3bt0syZM8Mea0CNjY165513lEqlwh4lYwsWLND27du1f/9+bdu2TefPn9fChQv16aefhj1av86cOaO6ujpNmzZN+/fvV3V1tdavX6/t27eHPVogu3fv1qVLl7RmzZrQZigI7chD7P8L7ft+aNW+kTz++OP65z//qbfffjvsUTLyla98RcePH9elS5f0yiuvqKqqSs3NzcM6Mu3t7dqwYYPeeOMNjRo1KuxxMlZZWdnz17Nnz1Z5ebmmTp2qP/zhD6qtrQ1xsv51d3errKxMW7ZskSTNmzdPra2tqqur0/e///2Qp8vc7373O1VWVqqkpCS0GSJ/BnPHHXdo5MiRfc5WLly40OesBkNr3bp12rNnj958801NnDgx7HEyUlhYqC996UsqKytTKpXS3Llz9atf/Srssfp17NgxXbhwQfPnz1dBQYEKCgrU3NysX//61yooKFBXV1fYI2Zk9OjRmj17tk6dOhX2KP2aMGFCn//gmDFjxrD/0tAXffDBBzpw4IB+9KMfhTpH5ANTWFio+fPn93xb4qqmpiYtXLgwpKnym+/7evzxx/Xqq6/qr3/9q0pLS8MeKWu+78vzvLDH6NfSpUt14sQJHT9+vGcpKyvTo48+quPHj2vkyJFhj5gRz/P03nvvacKECWGP0q9Fixb1+dr9+++/r8mTJ4c0UXANDQ0qLi7W8uXLQ50jLy6R1dbWavXq1SorK1N5ebnq6+vV1tam6urqsEfr1+XLl3X69Ome12fPntXx48c1ZswYTZo0KcTJ+ldTU6OdO3fqtddeU1FRUc/ZYyKR0M033xzydNf39NNPq7KyUslkUh0dHWpsbNTBgwe1b9++sEfrV1FRUZ/7W6NHj9bYsWOH9X2vJ598UitWrNCkSZN04cIFPf/883JdV1VVVWGP1q8nnnhCCxcu1JYtW/Td735X//jHP1RfX6/6+vqwR8tId3e3GhoaVFVVpYKCkD/iQ/numoHf/OY3/uTJk/3CwkL/nnvuicRXZt98801fUp+lqqoq7NH6da2ZJfkNDQ1hj9avtWvX9vwZGTdunL906VL/jTfeCHusrETha8oPP/ywP2HCBP+mm27yS0pK/G9/+9t+a2tr2GNl5PXXX/dnzZrlx2Ixf/r06X59fX3YI2Vs//79viT/5MmTYY/iO77v++GkDQCQzyJ/DwYAMDwRGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATeRUYz/O0efPmYf+rP/4fc+cWc+deVGdn7sHJqx+0dF1XiURC6XRa8Xg87HEyxty5xdy5F9XZmXtw8uoMBgAwfBAYAICJnP+qze7ubn388ccqKioa8geCua7b63+jgrlzi7lzL6qzM3dfvu+ro6NDJSUlGjGi/3OUnN+D+fDDD5VMJnN5SADAEGtvbx/wQYM5v0RWVFSU60MCAIZYJp/lOQ/MUF8WAwDkXiaf5dzkBwCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACayCsyLL76o0tJSjRo1SvPnz9dbb7011HMBACIucGBefvllbdy4Uc8884zeffddLVmyRJWVlWpra7OYDwAQUY7v+36QHRYsWKB77rlHdXV1PetmzJihlStXKpVKDbi/67pKJBLBJwUADBvpdFrxeLzfbQKdwVy5ckXHjh1TRUVFr/UVFRU6fPhw8AkBAHmrIMjGFy9eVFdXl8aPH99r/fjx43X+/Plr7uN5njzP63ntum4WYwIAoiarm/yO4/R67ft+n3VXpVIpJRKJniWZTGZzSABAxAQKzB133KGRI0f2OVu5cOFCn7OaqzZt2qR0Ot2ztLe3Zz8tACAyAgWmsLBQ8+fPV1NTU6/1TU1NWrhw4TX3icViisfjvRYAQP4LdA9Gkmpra7V69WqVlZWpvLxc9fX1amtrU3V1tcV8AICIChyYhx9+WJ9++ql+9rOf6dy5c5o1a5b+8pe/aPLkyRbzAQAiKvDPwQwWPwcDANE35D8HAwBApggMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMBH4iZaInhw/Uw5AHgvy0EjOYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYCByYQ4cOacWKFSopKZHjONq9e7fBWACAqAscmM7OTs2dO1cvvPCCxTwAgDxREHSHyspKVVZWWswCAMgjgQMTlOd58jyv57XrutaHBAAMA+Y3+VOplBKJRM+STCatDwkAGAbMA7Np0yal0+mepb293fqQAIBhwPwSWSwWUywWsz4MAGCY4edgAAAmAp/BXL58WadPn+55ffbsWR0/flxjxozRpEmThnQ4AEB0Ob7v+0F2OHjwoO6///4+66uqqvT73/9+wP1d11UikQhySAxSwH/EAHBdVz/D0+m04vF4v9sGPoP5xje+wQcWAGBA3IMBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJgI/D+ZGFtXn4DiOE/YIWYvqew6AMxgAgBECAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJgIFJpVK6d5771VRUZGKi4u1cuVKnTx50mo2AECEBQpMc3OzampqdOTIETU1Nenzzz9XRUWFOjs7reYDAESU4/u+n+3O//3vf1VcXKzm5mZ9/etfz2gf13WVSCSyPWSoBvFWhcpxnLBHyFpU33MgX139DE+n04rH4/1uWzCYA6XTaUnSmDFjrruN53nyPK/XcACA/Jf1TX7f91VbW6vFixdr1qxZ190ulUopkUj0LMlkMttDAgAiJOtLZDU1Nfrzn/+st99+WxMnTrzudtc6g4lqZKJ6uYZLZACGivklsnXr1mnPnj06dOhQv3GRpFgsplgsls1hAAARFigwvu9r3bp12rVrlw4ePKjS0lKruQAAERcoMDU1Ndq5c6dee+01FRUV6fz585KkRCKhm2++2WRAAEA0BboHc71r+Q0NDVqzZk1Gfw++ppx73IMBMFTM7sHwLzsAIFP8LjIAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwEeuAYoinKD4qL6tM4o/yeA0OFMxgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADARKDB1dXWaM2eO4vG44vG4ysvLtXfvXqvZAAARFigwEydO1NatW3X06FEdPXpUDzzwgB566CG1trZazQcAiCjH931/MH+DMWPG6Be/+IV++MMfZrS967pKJBKDOWRoBvlWIQuO44Q9Qlb4s4J8dfUzPJ1OKx6P97ttQbYH6erq0h//+Ed1dnaqvLz8utt5nifP83oNBwDIf4Fv8p84cUK33nqrYrGYqqurtWvXLs2cOfO626dSKSUSiZ4lmUwOamAAQDQEvkR25coVtbW16dKlS3rllVf00ksvqbm5+bqRudYZTFQjw2WP3OMSGTC8BLlENuh7MMuWLdPUqVP129/+NtBwUcSHRu4RGGB4CRKYQf8cjO/7vc5QAACQAt7kf/rpp1VZWalkMqmOjg41Njbq4MGD2rdvn9V8AICIChSYTz75RKtXr9a5c+eUSCQ0Z84c7du3Tw8++KDVfACAiBr0PZiguAeDILgHAwwvOb0HAwDAtRAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYCLQEy2BXIvqg7t4UBrAGQwAwAiBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJgYVGBSqZQcx9HGjRuHaBwAQL7IOjAtLS2qr6/XnDlzhnIeAECeyCowly9f1qOPPqpt27bp9ttvH+qZAAB5IKvA1NTUaPny5Vq2bNmA23qeJ9d1ey0AgPxXEHSHxsZGvfPOO2ppaclo+1Qqpeeeey7wYACAaAt0BtPe3q4NGzZox44dGjVqVEb7bNq0Sel0umdpb2/PalAAQLQ4vu/7mW68e/dufetb39LIkSN71nV1dclxHI0YMUKe5/X6/67FdV0lEonsJw5RgLcKNzjHccIeISv8GcdArn6Gp9NpxePxfrcNdIls6dKlOnHiRK91P/jBDzR9+nQ99dRTA8YFAHDjCBSYoqIizZo1q9e60aNHa+zYsX3WAwBubPwkPwDAROBvkf2/gwcPDsEYAIB8wxkMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmBv3AMQB9+b4f9ghZcRwn7BGyEtX3O99xBgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYCBWbz5s1yHKfXcuedd1rNBgCIsIKgO9x11106cOBAz+uRI0cO6UAAgPwQODAFBQWctQAABhT4HsypU6dUUlKi0tJSPfLIIzpz5ozFXACAiAt0BrNgwQJt375dX/7yl/XJJ5/o+eef18KFC9Xa2qqxY8decx/P8+R5Xs9r13UHNzEAIBIc3/f9bHfu7OzU1KlT9ZOf/ES1tbXX3Gbz5s167rnnsh5wOBnEWwVEguM4YY+QFf7dzB3XdZVIJJROpxWPx/vddlBfUx49erRmz56tU6dOXXebTZs2KZ1O9yzt7e2DOSQAICIC3+T/Is/z9N5772nJkiXX3SYWiykWiw3mMACACAp0BvPkk0+qublZZ8+e1d///nd95zvfkeu6qqqqspoPABBRgc5gPvzwQ33ve9/TxYsXNW7cOH3ta1/TkSNHNHnyZKv5AAARFSgwjY2NVnMAAPIMv4sMAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATgR44BiC/+b4f9ghZcRwn7BGyFtX3PBOcwQAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwETgwH330kR577DGNHTtWt9xyi+6++24dO3bMYjYAQIQVBNn4s88+06JFi3T//fdr7969Ki4u1n/+8x/ddtttRuMBAKIqUGB+/vOfK5lMqqGhoWfdlClThnomAEAeCHSJbM+ePSorK9OqVatUXFysefPmadu2bf3u43meXNfttQAA8l+gwJw5c0Z1dXWaNm2a9u/fr+rqaq1fv17bt2+/7j6pVEqJRKJnSSaTgx4aADD8Ob7v+5luXFhYqLKyMh0+fLhn3fr169XS0qK//e1v19zH8zx5ntfz2nXdyEYmwFsFIIccxwl7hKxF7XPFdV0lEgml02nF4/F+tw10BjNhwgTNnDmz17oZM2aora3tuvvEYjHF4/FeCwAg/wUKzKJFi3Ty5Mle695//31Nnjx5SIcCAERfoMA88cQTOnLkiLZs2aLTp09r586dqq+vV01NjdV8AICICnQPRpL+9Kc/adOmTTp16pRKS0tVW1urH//4xxnvf/X6XRRF7VopcKPgHkzuBLkHEzgwg0VgAAw1ApM7Zjf5AQDIFIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJgrCHgAABitqD+36oig/LG0gnMEAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAICJQIGZMmWKHMfps9TU1FjNBwCIqIIgG7e0tKirq6vn9b/+9S89+OCDWrVq1ZAPBgCItkCBGTduXK/XW7du1dSpU3XfffcN6VAAgOgLFJgvunLlinbs2KHa2lo5jnPd7TzPk+d5Pa9d1832kACACMn6Jv/u3bt16dIlrVmzpt/tUqmUEolEz5JMJrM9JAAgQhzf9/1sdvzmN7+pwsJCvf766/1ud60zmKhGJsu3CgCuq78rQMNZOp1WPB7vd5usLpF98MEHOnDggF599dUBt43FYorFYtkcBgAQYVldImtoaFBxcbGWL18+1PMAAPJE4MB0d3eroaFBVVVVKijI+jsCAIA8FzgwBw4cUFtbm9auXWsxDwAgT2R9kz9brusqkUjk8pBDhpv8AIZaPt/k53eRAQBMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABM5fyRllJ+p4rpu2CMAwLCQyWd5zgPT0dGR60MOmag+KA0AhlpHR8eAn4k5f6Jld3e3Pv74YxUVFQ35k9xc11UymVR7e/uAT1obTpg7t5g796I6O3P35fu+Ojo6VFJSohEj+r/LkvMzmBEjRmjixImmx4jH45H6w3AVc+cWc+deVGdn7t4yvZrDTX4AgAkCAwAwkVeBicVievbZZxWLxcIeJRDmzi3mzr2ozs7cg5Pzm/wAgBtDXp3BAACGDwIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABM/A/p3eutasyBTwAAAABJRU5ErkJggg==",
      "text/plain": [
       "Figure(PyObject <Figure size 640x480 with 1 Axes>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "PyObject <matplotlib.image.AxesImage object at 0x7ff070dc6340>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "spy(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# From https://docs.julialang.org/en/v1/stdlib/SparseArrays/\n",
    "#=\n",
    "struct SparseMatrixCSC{Tv,Ti<:Integer} <: AbstractSparseMatrix{Tv,Ti}\n",
    "    m::Int                  # Number of rows\n",
    "    n::Int                  # Number of columns\n",
    "    colptr::Vector{Ti}      # Column j is in colptr[j]:(colptr[j+1]-1)\n",
    "    rowval::Vector{Ti}      # Row indices of stored values\n",
    "    nzval::Vector{Tv}       # Stored values, typically nonzeros\n",
    "end\n",
    "=#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([1, 2, 3, 4, 5, 6, 7, 8], [1.0, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cp = A.colptr\n",
    "col1_inds = cp[1]:(cp[2]-1)\n",
    "A.rowval[col1_inds], A.nzval[col1_inds]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([1, 2], [0.01, 1.0])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "col2_inds = cp[2]:(cp[3]-1)\n",
    "A.rowval[col2_inds], A.nzval[col2_inds]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8×8 Array{Float64,2}:\n",
       " 1.0   0.01  0.01  0.01  0.01  0.01  0.01  0.01\n",
       " 0.01  1.0   0.0   0.0   0.0   0.0   0.0   0.0\n",
       " 0.01  0.0   1.0   0.0   0.0   0.0   0.0   0.0\n",
       " 0.01  0.0   0.0   1.0   0.0   0.0   0.0   0.0\n",
       " 0.01  0.0   0.0   0.0   1.0   0.0   0.0   0.0\n",
       " 0.01  0.0   0.0   0.0   0.0   1.0   0.0   0.0\n",
       " 0.01  0.0   0.0   0.0   0.0   0.0   1.0   0.0\n",
       " 0.01  0.0   0.0   0.0   0.0   0.0   0.0   1.0"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_full = Matrix(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8×8 Array{Float64,2}:\n",
       "  1.0007    -0.010007    -0.010007    …  -0.010007    -0.010007\n",
       " -0.010007   1.0001       0.00010007      0.00010007   0.00010007\n",
       " -0.010007   0.00010007   1.0001          0.00010007   0.00010007\n",
       " -0.010007   0.00010007   0.00010007      0.00010007   0.00010007\n",
       " -0.010007   0.00010007   0.00010007      0.00010007   0.00010007\n",
       " -0.010007   0.00010007   0.00010007  …   0.00010007   0.00010007\n",
       " -0.010007   0.00010007   0.00010007      1.0001       0.00010007\n",
       " -0.010007   0.00010007   0.00010007      0.00010007   1.0001"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inv(A_full)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LU{Float64,Array{Float64,2}}\n",
       "L factor:\n",
       "8×8 Array{Float64,2}:\n",
       " 1.0    0.0          0.0         …   0.0          0.0         0.0\n",
       " 0.01   1.0          0.0             0.0          0.0         0.0\n",
       " 0.01  -0.00010001   1.0             0.0          0.0         0.0\n",
       " 0.01  -0.00010001  -0.00010002      0.0          0.0         0.0\n",
       " 0.01  -0.00010001  -0.00010002      0.0          0.0         0.0\n",
       " 0.01  -0.00010001  -0.00010002  …   1.0          0.0         0.0\n",
       " 0.01  -0.00010001  -0.00010002     -0.00010005   1.0         0.0\n",
       " 0.01  -0.00010001  -0.00010002     -0.00010005  -0.00010006  1.0\n",
       "U factor:\n",
       "8×8 Array{Float64,2}:\n",
       " 1.0  0.01     0.01     0.01         0.01        …   0.01         0.01\n",
       " 0.0  0.9999  -0.0001  -0.0001      -0.0001         -0.0001      -0.0001\n",
       " 0.0  0.0      0.9999  -0.00010001  -0.00010001     -0.00010001  -0.00010001\n",
       " 0.0  0.0      0.0      0.9999      -0.00010002     -0.00010002  -0.00010002\n",
       " 0.0  0.0      0.0      0.0          0.9999         -0.00010003  -0.00010003\n",
       " 0.0  0.0      0.0      0.0          0.0         …  -0.00010004  -0.00010004\n",
       " 0.0  0.0      0.0      0.0          0.0             0.9999      -0.00010005\n",
       " 0.0  0.0      0.0      0.0          0.0             0.0          0.9999"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L, U, p = lu(A_full)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAGiCAYAAADa2tCeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAYD0lEQVR4nO3db2yV9f3/8dcFtQfFnqMgRRoO0DA2QP6I1LHyZ34V7NIQIlvGdFFWxrakSeWPjZlDb4iL4bAsu7HF2ay4dCOE1CwK4jbAkknRMLaCkjFmEAaxVUGGkXOV3riI7fW78YvNukJ7rtPzPlevw/ORXNnOtevq9Y4x57nPuU57Ob7v+wIAIMdGhD0AAKAwERgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAICJggnMiy++qPLyco0aNUrz58/XW2+9FfZIgzp06JBWrFihsrIyOY6j3bt3hz1SRlKplO69916VlJSotLRUK1eu1KlTp8Iea1ANDQ2aM2eO4vG44vG4KisrtXfv3rDHCiyVSslxHG3cuDHsUQa0efNmOY7TZ7vzzjvDHisjH330kR577DGNHTtWt9xyi+6++24dO3Ys7LEGNWXKlH7/zB3HUV1dXSjzFERgXn75ZW3cuFHPPPOM3n33XS1ZskTV1dVqb28Pe7QBdXV1ae7cuXrhhRfCHiWQ1tZW1dXV6ciRI2ppadHnn3+uqqoqdXV1hT3agCZOnKitW7fq6NGjOnr0qB544AE99NBDOnnyZNijZaytrU2NjY2aM2dO2KNk5K677tL58+d7txMnToQ90qA+++wzLVq0SDfddJP27t2rf/3rX/rFL36h2267LezRBtXW1tbnn3dLS4skadWqVeEM5BeAr371q35tbW2ffdOnT/d/8pOfhDRRcJL8Xbt2hT1GVi5evOhL8ltbW8MeJbDbb7/df+mll8IeIyOdnZ3+tGnT/JaWFv++++7zN2zYEPZIA3r22Wf9uXPnhj1GYE899ZS/ePHisMfIiQ0bNvhTp071e3p6Qrl+5FcwV69e1bFjx1RVVdVnf1VVlQ4fPhzSVDeWdDotSRozZkzIk2Suu7tbzc3N6urqUmVlZdjjZKSurk7Lly/XsmXLwh4lY6dPn1ZZWZnKy8v1yCOP6OzZs2GPNKg9e/aooqJCq1atUmlpqebNm6dt27aFPVZgV69e1Y4dO7R27Vo5jhPKDJEPzKVLl9Td3a3x48f32T9+/HhduHAhpKluHL7vq76+XosXL9asWbPCHmdQJ06c0K233qpYLKba2lrt2rVLM2fODHusQTU3N+udd95RKpUKe5SMLViwQNu3b9f+/fu1bds2XbhwQQsXLtSnn34a9mgDOnv2rBoaGjRt2jTt379ftbW1Wr9+vbZv3x72aIHs3r1bly9f1po1a0KboSi0K+fY/xba9/3Qqn0jefzxx/WPf/xDb7/9dtijZOQrX/mKjh8/rsuXL+uVV15RTU2NWltbh3VkOjo6tGHDBr3xxhsaNWpU2ONkrLq6uve/z549W5WVlZo6dap+//vfq76+PsTJBtbT06OKigpt2bJFkjRv3jydPHlSDQ0N+t73vhfydJn77W9/q+rqapWVlYU2Q+RXMHfccYdGjhzZb7Vy8eLFfqsa5Na6deu0Z88evfnmm5o4cWLY42SkuLhYX/rSl1RRUaFUKqW5c+fql7/8ZdhjDejYsWO6ePGi5s+fr6KiIhUVFam1tVW/+tWvVFRUpO7u7rBHzMjo0aM1e/ZsnT59OuxRBjRhwoR+/4djxowZw/5LQ//tgw8+0IEDB/TDH/4w1DkiH5ji4mLNnz+/99sSX2hpadHChQtDmqqw+b6vxx9/XK+++qr+8pe/qLy8POyRsub7vjzPC3uMAS1dulQnTpzQ8ePHe7eKigo9+uijOn78uEaOHBn2iBnxPE/vvfeeJkyYEPYoA1q0aFG/r92///77mjx5ckgTBdfU1KTS0lItX7481DkK4iOy+vp6rV69WhUVFaqsrFRjY6Pa29tVW1sb9mgDunLlis6cOdP7+ty5czp+/LjGjBmjSZMmhTjZwOrq6rRz50699tprKikp6V09JhIJ3XzzzSFPd31PP/20qqurlUwm1dnZqebmZh08eFD79u0Le7QBlZSU9Lu/NXr0aI0dO3ZY3/d68skntWLFCk2aNEkXL17U888/L9d1VVNTE/ZoA3riiSe0cOFCbdmyRd/5znf097//XY2NjWpsbAx7tIz09PSoqalJNTU1KioK+S0+lO+uGfj1r3/tT5482S8uLvbvueeeSHxl9s033/Ql9dtqamrCHm1A15pZkt/U1BT2aANau3Zt778j48aN85cuXeq/8cYbYY+VlSh8Tfnhhx/2J0yY4N90001+WVmZ/61vfcs/efJk2GNl5PXXX/dnzZrlx2Ixf/r06X5jY2PYI2Vs//79viT/1KlTYY/iO77v++GkDQBQyCJ/DwYAMDwRGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBRUYz/O0efPmYf+nP/4Xc+cXc+dfVGdn7qEpqF+0dF1XiURC6XRa8Xg87HEyxtz5xdz5F9XZmXtoCmoFAwAYPggMAMBE3v/UZk9Pjz7++GOVlJTk/IFgruv2+c+oYO78Yu78i+rszN2f7/vq7OxUWVmZRowYeI2S93swH374oZLJZD4vCQDIsY6OjkEfNJj3j8hKSkryfUkAQI5l8l6e98Dk+mMxAED+ZfJezk1+AIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYCKrwLz44osqLy/XqFGjNH/+fL311lu5ngsAEHGBA/Pyyy9r48aNeuaZZ/Tuu+9qyZIlqq6uVnt7u8V8AICIcnzf94OcsGDBAt1zzz1qaGjo3TdjxgytXLlSqVRq0PNd11UikQg+KQBg2Ein04rH4wMeE2gFc/XqVR07dkxVVVV99ldVVenw4cPBJwQAFKyiIAdfunRJ3d3dGj9+fJ/948eP14ULF655jud58jyv97XrulmMCQCImqxu8juO0+e17/v99n0hlUopkUj0bslkMptLAgAiJlBg7rjjDo0cObLfauXixYv9VjVf2LRpk9LpdO/W0dGR/bQAgMgIFJji4mLNnz9fLS0tffa3tLRo4cKF1zwnFospHo/32QAAhS/QPRhJqq+v1+rVq1VRUaHKyko1Njaqvb1dtbW1FvMBACIqcGAefvhhffrpp/rpT3+q8+fPa9asWfrzn/+syZMnW8wHAIiowL8HM1T8HgwARF/Ofw8GAIBMERgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBgIvATLXMlk4fVDDeO44Q9AgBEBisYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYCB+bQoUNasWKFysrK5DiOdu/ebTAWACDqAgemq6tLc+fO1QsvvGAxDwCgQBQFPaG6ulrV1dUWswAACkjgwATleZ48z+t97bqu9SUBAMOA+U3+VCqlRCLRuyWTSetLAgCGAfPAbNq0Sel0unfr6OiwviQAYBgw/4gsFospFotZXwYAMMzwezAAABOBVzBXrlzRmTNnel+fO3dOx48f15gxYzRp0qScDgcAiC7H930/yAkHDx7U/fff329/TU2Nfve73w16vuu6SiQSSqfTisfjQS4dOsdxwh4BAIaFTN7DA69g/u///k8BmwQAuAFxDwYAYILAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYCPw8mBtZVJ+Dw4PSAISBFQwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAE4ECk0qldO+996qkpESlpaVauXKlTp06ZTUbACDCAgWmtbVVdXV1OnLkiFpaWvT555+rqqpKXV1dVvMBACLK8X3fz/bk//znPyotLVVra6u+/vWvZ3SO67pKJBJKp9OKx+PZXhoBOI4T9ggACkwm7+FFQ72AJI0ZM+a6x3ieJ8/zel+7rjuUSwIAIiLrm/y+76u+vl6LFy/WrFmzrntcKpVSIpHo3ZLJZLaXBABESNYfkdXV1elPf/qT3n77bU2cOPG6x11rBZNMJvmILI/4iAxArpl9RLZu3Trt2bNHhw4dGjAukhSLxRSLxbK5DAAgwgIFxvd9rVu3Trt27dLBgwdVXl5uNRcAIOICBaaurk47d+7Ua6+9ppKSEl24cEGSlEgkdPPNN5sMCACIpkD3YK73WX5TU5PWrFmT0c/ga8r5xz0YALmW83swQ/iVGQDADYa/RQYAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgIlADxxDNEX5QXE8jROILlYwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBgIlBgGhoaNGfOHMXjccXjcVVWVmrv3r1WswEAIixQYCZOnKitW7fq6NGjOnr0qB544AE99NBDOnnypNV8AICIcnzf94fyA8aMGaOf//zn+sEPfpDR8a7rKpFIKJ1OKx6PD+XSuAE4jhP2CACuIZP38KJsf3h3d7f+8Ic/qKurS5WVldc9zvM8eZ7X+9p13WwvCQCIkMA3+U+cOKFbb71VsVhMtbW12rVrl2bOnHnd41OplBKJRO+WTCaHNDAAIBoCf0R29epVtbe36/Lly3rllVf00ksvqbW19bqRudYKJplM8hEZMsJHZMDwlMl7+JDvwSxbtkxTp07Vb37zm4yO5x4MgiAwwPCUyXv4kH8Pxvf9PisUAACkgDf5n376aVVXVyuZTKqzs1PNzc06ePCg9u3bZzUfACCiAgXmk08+0erVq3X+/HklEgnNmTNH+/bt04MPPmg1HwAgooZ8DyYo7sEgCO7BAMNTXu7BAABwLQQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmAj0REsg3/L8PLyc4UFpACsYAIARAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwMaTApFIpOY6jjRs35mgcAEChyDowbW1tamxs1Jw5c3I5DwCgQGQVmCtXrujRRx/Vtm3bdPvtt+d6JgBAAcgqMHV1dVq+fLmWLVs26LGe58l13T4bAKDwFQU9obm5We+8847a2toyOj6VSum5554LPBgAINoCrWA6Ojq0YcMG7dixQ6NGjcronE2bNimdTvduHR0dWQ0KAIgWx/d9P9ODd+/erW9+85saOXJk777u7m45jqMRI0bI87w+/9u1uK6rRCKhdDqteDye/eTAMOY4TtgjAKYyeQ8P9BHZ0qVLdeLEiT77vv/972v69Ol66qmnBo0LAODGESgwJSUlmjVrVp99o0eP1tixY/vtBwDc2PhNfgCAicDfIvtfBw8ezMEYAIBCwwoGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATQ37gGID+fN8Pe4SsOI4T9ggoIKxgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYIDADARKDAbN68WY7j9NnuvPNOq9kAABFWFPSEu+66SwcOHOh9PXLkyJwOBAAoDIEDU1RUxKoFADCowPdgTp8+rbKyMpWXl+uRRx7R2bNnLeYCAERcoBXMggULtH37dn35y1/WJ598oueff14LFy7UyZMnNXbs2Gue43mePM/rfe267tAmBgBEguP7vp/tyV1dXZo6dap+/OMfq76+/prHbN68Wc8991y//el0WvF4PNtLAzDgOE7YIyAiMnkPH9LXlEePHq3Zs2fr9OnT1z1m06ZNSqfTvVtHR8dQLgkAiIjAN/n/m+d5eu+997RkyZLrHhOLxRSLxYZyGQBABAVawTz55JNqbW3VuXPn9Le//U3f/va35bquampqrOYDAERUoBXMhx9+qO9+97u6dOmSxo0bp6997Ws6cuSIJk+ebDUfACCiAgWmubnZag4AQIHhb5EBAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGAi0APHABQ23/fDHiErjuOEPQKugRUMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABOBA/PRRx/pscce09ixY3XLLbfo7rvv1rFjxyxmAwBEWFGQgz/77DMtWrRI999/v/bu3avS0lL9+9//1m233WY0HgAgqgIF5mc/+5mSyaSampp6902ZMiXXMwEACkCgj8j27NmjiooKrVq1SqWlpZo3b562bds24Dme58l13T4bAKDwBQrM2bNn1dDQoGnTpmn//v2qra3V+vXrtX379uuek0qllEgkerdkMjnkoQEAw5/j+76f6cHFxcWqqKjQ4cOHe/etX79ebW1t+utf/3rNczzPk+d5va9d11UymVQ6nVY8Hh/C6ADw/zmOE/YIN5xM3sMDrWAmTJigmTNn9tk3Y8YMtbe3X/ecWCymeDzeZwMAFL5AgVm0aJFOnTrVZ9/777+vyZMn53QoAED0BQrME088oSNHjmjLli06c+aMdu7cqcbGRtXV1VnNBwCIqED3YCTpj3/8ozZt2qTTp0+rvLxc9fX1+tGPfpTx+a7rKpFIcA8GQM5wDyb/MnkPDxyYoSIwAHKNwORfzm/yAwCQKQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBSFPQAADFWen5uYU4X8sDRWMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYCJQYKZMmSLHcfptdXV1VvMBACKqKMjBbW1t6u7u7n39z3/+Uw8++KBWrVqV88EAANEWKDDjxo3r83rr1q2aOnWq7rvvvpwOBQCIvkCB+W9Xr17Vjh07VF9fL8dxrnuc53nyPK/3teu62V4SABAhWd/k3717ty5fvqw1a9YMeFwqlVIikejdkslktpcEAESI4/u+n82J3/jGN1RcXKzXX399wOOutYJJJpNKp9OKx+PZXBoACsZAnwANZ5m8h2f1EdkHH3ygAwcO6NVXXx302Fgsplgsls1lAAARltVHZE1NTSotLdXy5ctzPQ8AoEAEDkxPT4+amppUU1OjoqKsvyMAAChwgQNz4MABtbe3a+3atRbzAAAKROAlSFVVlbL8XgAA4AbC3yIDAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJvL+SMovniXjum6+Lw0AyJFMnguW98B0dnZKkpLJZL4vDQDIkc7OTiUSiQGPcfw8P56yp6dHH3/8sUpKSuQ4Tk5/tuu6SiaT6ujoUDwez+nPtsTc+cXc+RfV2Zm7P9/31dnZqbKyMo0YMfBdlryvYEaMGKGJEyeaXiMej0fqX4YvMHd+MXf+RXV25u5rsJXLF7jJDwAwQWAAACYKKjCxWEzPPvusYrFY2KMEwtz5xdz5F9XZmXto8n6THwBwYyioFQwAYPggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwMT/A6FUhbwu1I1PAAAAAElFTkSuQmCC",
      "text/plain": [
       "Figure(PyObject <Figure size 640x480 with 1 Axes>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "PyObject <matplotlib.image.AxesImage object at 0x7ff04045f040>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "spy(U)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAGiCAYAAADa2tCeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAYdElEQVR4nO3df2wXhf3H8ddB7QfFfj4KUqThAzSMDZAfInWs/JhTsEtDiGwZ00VZGduSJpUfNmYO/UNcDB+WZX9scTYrLt0IITWLgrgNsGRSNIytoGSMGYRBbFWQYeRzpX8csb3vH9+sWVdoP/fp5/253ofnI7nMz2d33Du1vad39ynn+L7vCwCAHBsR9gAAgMJEYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYKJjAvvviiysvLNWrUKM2fP19vvfVW2CMN6tChQ1qxYoXKysrkOI52794d9kgZSaVSuvfee1VSUqLS0lKtXLlSp06dCnusQTU0NGjOnDmKx+OKx+OqrKzU3r17wx4rsFQqJcdxtHHjxrBHGdDmzZvlOE6f5c477wx7rIx89NFHeuyxxzR27Fjdcsstuvvuu3Xs2LGwxxrUlClT+n3NHcdRXV1dKPMURGBefvllbdy4Uc8884zeffddLVmyRNXV1Wpvbw97tAF1dXVp7ty5euGFF8IeJZDW1lbV1dXpyJEjamlp0eeff66qqip1dXWFPdqAJk6cqK1bt+ro0aM6evSoHnjgAT300EM6efJk2KNlrK2tTY2NjZozZ07Yo2Tkrrvu0vnz53uXEydOhD3SoD777DMtWrRIN910k/bu3at//vOf+vnPf67bbrst7NEG1dbW1ufr3dLSIklatWpVOAP5BeDLX/6yX1tb2+e96dOn+z/+8Y9Dmig4Sf6uXbvCHiMrFy9e9CX5ra2tYY8S2O233+6/9NJLYY+Rkc7OTn/atGl+S0uLf9999/kbNmwIe6QBPfvss/7cuXPDHiOwp556yl+8eHHYY+TEhg0b/KlTp/o9PT2h7D/yZzBXr17VsWPHVFVV1ef9qqoqHT58OKSpbizpdFqSNGbMmJAnyVx3d7eam5vV1dWlysrKsMfJSF1dnZYvX65ly5aFPUrGTp8+rbKyMpWXl+uRRx7R2bNnwx5pUHv27FFFRYVWrVql0tJSzZs3T9u2bQt7rMCuXr2qHTt2aO3atXIcJ5QZIh+YS5cuqbu7W+PHj+/z/vjx43XhwoWQprpx+L6v+vp6LV68WLNmzQp7nEGdOHFCt956q2KxmGpra7Vr1y7NnDkz7LEG1dzcrHfeeUepVCrsUTK2YMECbd++Xfv379e2bdt04cIFLVy4UJ9++mnYow3o7Nmzamho0LRp07R//37V1tZq/fr12r59e9ijBbJ7925dvnxZa9asCW2GotD2nGP/W2jf90Or9o3k8ccf19///ne9/fbbYY+SkS996Us6fvy4Ll++rFdeeUU1NTVqbW0d1pHp6OjQhg0b9MYbb2jUqFFhj5Ox6urq3n+ePXu2KisrNXXqVP3ud79TfX19iJMNrKenRxUVFdqyZYskad68eTp58qQaGhr03e9+N+TpMveb3/xG1dXVKisrC22GyJ/B3HHHHRo5cmS/s5WLFy/2O6tBbq1bt0579uzRm2++qYkTJ4Y9TkaKi4v1hS98QRUVFUqlUpo7d65+8YtfhD3WgI4dO6aLFy9q/vz5KioqUlFRkVpbW/XLX/5SRUVF6u7uDnvEjIwePVqzZ8/W6dOnwx5lQBMmTOj3HxwzZswY9h8a+m8ffPCBDhw4oB/84AehzhH5wBQXF2v+/Pm9n5b4j5aWFi1cuDCkqQqb7/t6/PHH9eqrr+rPf/6zysvLwx4pa77vy/O8sMcY0NKlS3XixAkdP368d6moqNCjjz6q48ePa+TIkWGPmBHP8/Tee+9pwoQJYY8yoEWLFvX72P3777+vyZMnhzRRcE1NTSotLdXy5ctDnaMgLpHV19dr9erVqqioUGVlpRobG9Xe3q7a2tqwRxvQlStXdObMmd7X586d0/HjxzVmzBhNmjQpxMkGVldXp507d+q1115TSUlJ79ljIpHQzTffHPJ01/f000+rurpayWRSnZ2dam5u1sGDB7Vv376wRxtQSUlJv/tbo0eP1tixY4f1fa8nn3xSK1as0KRJk3Tx4kU9//zzcl1XNTU1YY82oCeeeEILFy7Uli1b9O1vf1t/+9vf1NjYqMbGxrBHy0hPT4+amppUU1OjoqKQD/GhfHbNwK9+9St/8uTJfnFxsX/PPfdE4iOzb775pi+p31JTUxP2aAO61syS/KamprBHG9DatWt7v0fGjRvnL1261H/jjTfCHisrUfiY8sMPP+xPmDDBv+mmm/yysjL/m9/8pn/y5Mmwx8rI66+/7s+aNcuPxWL+9OnT/cbGxrBHytj+/ft9Sf6pU6fCHsV3fN/3w0kbAKCQRf4eDABgeCIwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYKKjCe52nz5s3D/q/++F/MnV/MnX9RnZ25h6agftHSdV0lEgml02nF4/Gwx8kYc+cXc+dfVGdn7qEpqDMYAMDwQWAAACby/ldt9vT06OOPP1ZJSUnOHwjmum6f/40K5s4v5s6/qM7O3P35vq/Ozk6VlZVpxIiBz1Hyfg/mww8/VDKZzOcuAQA51tHRMeiDBvN+iaykpCTfuwQA5Fgmx/K8BybXl8UAAPmXybGcm/wAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYIDADARFaBefHFF1VeXq5Ro0Zp/vz5euutt3I9FwAg4gIH5uWXX9bGjRv1zDPP6N1339WSJUtUXV2t9vZ2i/kAABHl+L7vB9lgwYIFuueee9TQ0ND73owZM7Ry5UqlUqlBt3ddV4lEIvikAIBhI51OKx6PD7hOoDOYq1ev6tixY6qqqurzflVVlQ4fPhx8QgBAwSoKsvKlS5fU3d2t8ePH93l//PjxunDhwjW38TxPnuf1vnZdN4sxAQBRk9VNfsdx+rz2fb/fe/+RSqWUSCR6l2Qymc0uAQAREygwd9xxh0aOHNnvbOXixYv9zmr+Y9OmTUqn071LR0dH9tMCACIjUGCKi4s1f/58tbS09Hm/paVFCxcuvOY2sVhM8Xi8zwIAKHyB7sFIUn19vVavXq2KigpVVlaqsbFR7e3tqq2ttZgPABBRgQPz8MMP69NPP9VPfvITnT9/XrNmzdKf/vQnTZ482WI+AEBEBf49mKHi92AAIPpy/nswAABkisAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAE4GfaJkrmTysZrhxHCfsEbKS52fKAQggqseVTHAGAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMBE4MAcOnRIK1asUFlZmRzH0e7duw3GAgBEXeDAdHV1ae7cuXrhhRcs5gEAFIiioBtUV1erurraYhYAQAEJHJigPM+T53m9r13Xtd4lAGAYML/Jn0qllEgkepdkMmm9SwDAMGAemE2bNimdTvcuHR0d1rsEAAwD5pfIYrGYYrGY9W4AAMMMvwcDADAR+AzmypUrOnPmTO/rc+fO6fjx4xozZowmTZqU0+EAANHl+L7vB9ng4MGDuv/++/u9X1NTo9/+9reDbu+6rhKJhNLptOLxeJBdh85xnLBHyErAf8UA8iiqx5VMjuGBz2C+9rWvccACAAyKezAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYIDADARODnwdzIovocnKg+0EiK7tccAGcwAAAjBAYAYILAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwECkwqldK9996rkpISlZaWauXKlTp16pTVbACACAsUmNbWVtXV1enIkSNqaWnR559/rqqqKnV1dVnNBwCIKMf3fT/bjf/973+rtLRUra2t+upXv5rRNq7rKpFIKJ1OKx6PZ7trBOA4TtgjZG0I355AJET15zOTY3jRUHcgSWPGjLnuOp7nyfO83teu6w5llwCAiMj6Jr/v+6qvr9fixYs1a9as666XSqWUSCR6l2Qyme0uAQARkvUlsrq6Ov3xj3/U22+/rYkTJ153vWudwSSTSS6R5VFUT8ElLpGh8EX159PsEtm6deu0Z88eHTp0aMC4SFIsFlMsFstmNwCACAsUGN/3tW7dOu3atUsHDx5UeXm51VwAgIgLFJi6ujrt3LlTr732mkpKSnThwgVJUiKR0M0332wyIAAgmgLdg7netcKmpiatWbMmoz+DjynnX1Sv8Urcg0Hhi+rPZ87vwfDDDgDIFH8XGQDABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJgI9cAzRFOUHxUX1aX9R/poDucIZDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmAgUmIaGBs2ZM0fxeFzxeFyVlZXau3ev1WwAgAgLFJiJEydq69atOnr0qI4ePaoHHnhADz30kE6ePGk1HwAgohzf9/2h/AFjxozRz372M33/+9/PaH3XdZVIJJROpxWPx4eya9wAHMcJe4SsDPHHCjeQqH6PZ3IML8r2D+/u7tbvf/97dXV1qbKy8rrreZ4nz/N6X7uum+0uAQAREvgm/4kTJ3TrrbcqFouptrZWu3bt0syZM6+7fiqVUiKR6F2SyeSQBgYAREPgS2RXr15Ve3u7Ll++rFdeeUUvvfSSWltbrxuZa53BJJNJLpEhI1G9fMAlMmQqqt/jmRzDh3wPZtmyZZo6dap+/etfZ7Q+92AQRFR/+AgMMhXV7/FMjuFD/j0Y3/f7nKEAACAFvMn/9NNPq7q6WslkUp2dnWpubtbBgwe1b98+q/kAABEVKDCffPKJVq9erfPnzyuRSGjOnDnat2+fHnzwQav5AAARNeR7MEFxDwZBRPX6NPdgkKmofo/n5R4MAADXQmAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACAiUBPtATyLaoP7orqQ6Si+vXG8MQZDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmBhSYFKplBzH0caNG3M0DgCgUGQdmLa2NjU2NmrOnDm5nAcAUCCyCsyVK1f06KOPatu2bbr99ttzPRMAoABkFZi6ujotX75cy5YtG3Rdz/Pkum6fBQBQ+IqCbtDc3Kx33nlHbW1tGa2fSqX03HPPBR4MABBtgc5gOjo6tGHDBu3YsUOjRo3KaJtNmzYpnU73Lh0dHVkNCgCIFsf3fT/TlXfv3q1vfOMbGjlyZO973d3dchxHI0aMkOd5ff6/a3FdV4lEQul0WvF4PPvJgWHMcZywR8hKgMMBciSq3yuZHMMDXSJbunSpTpw40ee9733ve5o+fbqeeuqpQeMCALhxBApMSUmJZs2a1ee90aNHa+zYsf3eBwDc2PhNfgCAicCfIvtfBw8ezMEYAIBCwxkMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmhvzAMQD9+b4f9ghZcRwn7BGyEtWvd6HjDAYAYILAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwECszmzZvlOE6f5c4777SaDQAQYUVBN7jrrrt04MCB3tcjR47M6UAAgMIQODBFRUWctQAABhX4Hszp06dVVlam8vJyPfLIIzp79qzFXACAiAt0BrNgwQJt375dX/ziF/XJJ5/o+eef18KFC3Xy5EmNHTv2mtt4nifP83pfu647tIkBAJHg+L7vZ7txV1eXpk6dqh/96Eeqr6+/5jqbN2/Wc8891+/9dDqteDye7a4BGHAcJ+wRsjKEw1joovo1z+QYPqSPKY8ePVqzZ8/W6dOnr7vOpk2blE6ne5eOjo6h7BIAEBGBb/L/N8/z9N5772nJkiXXXScWiykWiw1lNwCACAp0BvPkk0+qtbVV586d01//+ld961vfkuu6qqmpsZoPABBRgc5gPvzwQ33nO9/RpUuXNG7cOH3lK1/RkSNHNHnyZKv5AAARFSgwzc3NVnMAAAoMfxcZAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmAj1wDEBh830/7BGy4jhO2CPgGjiDAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGAicGA++ugjPfbYYxo7dqxuueUW3X333Tp27JjFbACACCsKsvJnn32mRYsW6f7779fevXtVWlqqf/3rX7rtttuMxgMARFWgwPz0pz9VMplUU1NT73tTpkzJ9UwAgAIQ6BLZnj17VFFRoVWrVqm0tFTz5s3Ttm3bBtzG8zy5rttnAQAUvkCBOXv2rBoaGjRt2jTt379ftbW1Wr9+vbZv337dbVKplBKJRO+STCaHPDQAYPhzfN/3M125uLhYFRUVOnz4cO9769evV1tbm/7yl79ccxvP8+R5Xu9r13WVTCaVTqcVj8eHMDoA/D/HccIe4YaTyTE80BnMhAkTNHPmzD7vzZgxQ+3t7dfdJhaLKR6P91kAAIUvUGAWLVqkU6dO9Xnv/fff1+TJk3M6FAAg+gIF5oknntCRI0e0ZcsWnTlzRjt37lRjY6Pq6uqs5gMARFSgezCS9Ic//EGbNm3S6dOnVV5ervr6ev3whz/MeHvXdZVIJLgHAyBnuAeTf5kcwwMHZqgIDIBcIzD5l/Ob/AAAZIrAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABNFYe04kUiEteus5fnZbAAyFOWfzUJ+WBpnMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYCJQYKZMmSLHcfotdXV1VvMBACKqKMjKbW1t6u7u7n39j3/8Qw8++KBWrVqV88EAANEWKDDjxo3r83rr1q2aOnWq7rvvvpwOBQCIvkCB+W9Xr17Vjh07VF9fL8dxrrue53nyPK/3teu62e4SABAhWd/k3717ty5fvqw1a9YMuF4qlVIikehdkslktrsEAESI4/u+n82GX//611VcXKzXX399wPWudQYT1chk+aUCgOsa6ArQcJZOpxWPxwdcJ6tLZB988IEOHDigV199ddB1Y7GYYrFYNrsBAERYVpfImpqaVFpaquXLl+d6HgBAgQgcmJ6eHjU1NammpkZFRVl/RgAAUOACB+bAgQNqb2/X2rVrLeYBABSIrG/yZ8t1XSUSiXzuMme4yQ8g1wr5Jj9/FxkAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwkfdHUkb5mSqu64Y9AgAMC5kcy/MemM7OznzvMmei+qA0AMi1zs7OQY+JeX+iZU9Pjz7++GOVlJTk/ElurusqmUyqo6Nj0CetDSfMnV/MnX9RnZ25+/N9X52dnSorK9OIEQPfZcn7GcyIESM0ceJE033E4/FIfTP8B3PnF3PnX1RnZ+6+Mr2aw01+AIAJAgMAMFFQgYnFYnr22WcVi8XCHiUQ5s4v5s6/qM7O3EOT95v8AIAbQ0GdwQAAhg8CAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATPwfRSK+2pt88fIAAAAASUVORK5CYII=",
      "text/plain": [
       "Figure(PyObject <Figure size 640x480 with 1 Axes>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "PyObject <matplotlib.image.AxesImage object at 0x7ff0804b3b50>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q = collect(1:n)\n",
    "q[1] = n\n",
    "q[n] = 1\n",
    "AQ = A[:, q]\n",
    "spy(AQ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LU{Float64,Array{Float64,2}}\n",
       "L factor:\n",
       "8×8 Array{Float64,2}:\n",
       " 1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0\n",
       " 0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0\n",
       " 0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0\n",
       " 0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0\n",
       " 0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0\n",
       " 0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0\n",
       " 0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0\n",
       " 0.01  0.01  0.01  0.01  0.01  0.01  0.01  1.0\n",
       "U factor:\n",
       "8×8 Array{Float64,2}:\n",
       " 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.01\n",
       " 0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.01\n",
       " 0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.01\n",
       " 0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.01\n",
       " 0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.01\n",
       " 0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.01\n",
       " 0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.01\n",
       " 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.9993"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "AQ_full = Matrix(AQ)\n",
    "L, U, p = lu(AQ_full)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAGiCAYAAADa2tCeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAYWklEQVR4nO3df2xVhf3/8deB2oti71WQIg0XaBgbID9E6lj5Madgl4YQ2TKmi7IytiVNKj9szBz6h7gYLsuyP7Y4mxWXboSQmkVB3AZYMikaxlZQMsYMwiC2Ksgwck/pH4fYns8f36zfdYX2ntv7vqfn8nwkJ/PencN5p7bn6TnnluP4vu8LAIAcGxH2AACAwkRgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJgomMC+++KLKy8s1atQozZ8/X2+99VbYIw3q0KFDWrFihcrKyuQ4jnbv3h32SBlJpVK69957VVJSotLSUq1cuVKnTp0Ke6xBNTQ0aM6cOYrH44rH46qsrNTevXvDHiuwVColx3G0cePGsEcZ0ObNm+U4Tp/lzjvvDHusjHz00Ud67LHHNHbsWN1yyy26++67dezYsbDHGtSUKVP6fc0dx1FdXV0o8xREYF5++WVt3LhRzzzzjN59910tWbJE1dXVam9vD3u0AXV1dWnu3Ll64YUXwh4lkNbWVtXV1enIkSNqaWnR559/rqqqKnV1dYU92oAmTpyorVu36ujRozp69KgeeOABPfTQQzp58mTYo2Wsra1NjY2NmjNnTtijZOSuu+7S+fPne5cTJ06EPdKgPvvsMy1atEg33XST9u7dq3/+85/6+c9/rttuuy3s0QbV1tbW5+vd0tIiSVq1alU4A/kF4Mtf/rJfW1vb573p06f7P/7xj0OaKDhJ/q5du8IeIysXL170Jfmtra1hjxLY7bff7r/00kthj5GRzs5Of9q0aX5LS4t/3333+Rs2bAh7pAE9++yz/ty5c8MeI7CnnnrKX7x4cdhj5MSGDRv8qVOn+j09PaHsP/JnMFevXtWxY8dUVVXV5/2qqiodPnw4pKluLOl0WpI0ZsyYkCfJXHd3t5qbm9XV1aXKysqwx8lIXV2dli9frmXLloU9SsZOnz6tsrIylZeX65FHHtHZs2fDHmlQe/bsUUVFhVatWqXS0lLNmzdP27ZtC3uswK5evaodO3Zo7dq1chwnlBkiH5hLly6pu7tb48eP7/P++PHjdeHChZCmunH4vq/6+notXrxYs2bNCnucQZ04cUK33nqrYrGYamtrtWvXLs2cOTPssQbV3Nysd955R6lUKuxRMrZgwQJt375d+/fv17Zt23ThwgUtXLhQn376adijDejs2bNqaGjQtGnTtH//ftXW1mr9+vXavn172KMFsnv3bl2+fFlr1qwJbYai0PacY/9baN/3Q6v2jeTxxx/X3//+d7399tthj5KRL33pSzp+/LguX76sV155RTU1NWptbR3Wkeno6NCGDRv0xhtvaNSoUWGPk7Hq6uref549e7YqKys1depU/e53v1N9fX2Ikw2sp6dHFRUV2rJliyRp3rx5OnnypBoaGvTd73435Oky95vf/EbV1dUqKysLbYbIn8HccccdGjlyZL+zlYsXL/Y7q0FurVu3Tnv27NGbb76piRMnhj1ORoqLi/WFL3xBFRUVSqVSmjt3rn7xi1+EPdaAjh07posXL2r+/PkqKipSUVGRWltb9ctf/lJFRUXq7u4Oe8SMjB49WrNnz9bp06fDHmVAEyZM6PcfHDNmzBj2Hxr6bx988IEOHDigH/zgB6HOEfnAFBcXa/78+b2flviPlpYWLVy4MKSpCpvv+3r88cf16quv6s9//rPKy8vDHilrvu/L87ywxxjQ0qVLdeLECR0/frx3qaio0KOPPqrjx49r5MiRYY+YEc/z9N5772nChAlhjzKgRYsW9fvY/fvvv6/JkyeHNFFwTU1NKi0t1fLly0OdoyAukdXX12v16tWqqKhQZWWlGhsb1d7ertra2rBHG9CVK1d05syZ3tfnzp3T8ePHNWbMGE2aNCnEyQZWV1ennTt36rXXXlNJSUnv2WMikdDNN98c8nTX9/TTT6u6ulrJZFKdnZ1qbm7WwYMHtW/fvrBHG1BJSUm/+1ujR4/W2LFjh/V9ryeffFIrVqzQpEmTdPHiRT3//PNyXVc1NTVhjzagJ554QgsXLtSWLVv07W9/W3/729/U2NioxsbGsEfLSE9Pj5qamlRTU6OiopAP8aF8ds3Ar371K3/y5Ml+cXGxf88990TiI7NvvvmmL6nfUlNTE/ZoA7rWzJL8pqamsEcb0Nq1a3u/R8aNG+cvXbrUf+ONN8IeKytR+Jjyww8/7E+YMMG/6aab/LKyMv+b3/ymf/LkybDHysjrr7/uz5o1y4/FYv706dP9xsbGsEfK2P79+31J/qlTp8IexXd83/fDSRsAoJBF/h4MAGB4IjAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJgoqMJ7nafPmzcP+r/74X8ydX8ydf1GdnbmHpqB+0dJ1XSUSCaXTacXj8bDHyRhz5xdz519UZ2fuoSmoMxgAwPBBYAAAJvL+V2329PTo448/VklJSc4fCOa6bp//jQrmzi/mzr+ozs7c/fm+r87OTpWVlWnEiIHPUfJ+D+bDDz9UMpnM5y4BADnW0dEx6IMG834GU1JSku9d5kw6nQ57BAAFJpFIhD1CVjI5luc9MLm+LJZPUfoUCQBYyuRYzk1+AIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYCKrwLz44osqLy/XqFGjNH/+fL311lu5ngsAEHGBA/Pyyy9r48aNeuaZZ/Tuu+9qyZIlqq6uVnt7u8V8AICIcnzf94NssGDBAt1zzz1qaGjofW/GjBlauXKlUqnUoNu7rhvZR4QG/FIBwKCi+pTfdDo96FN+A53BXL16VceOHVNVVVWf96uqqnT48OHgEwIAClZRkJUvXbqk7u5ujR8/vs/748eP14ULF665jed58jyv97XrulmMCQCImqxu8v/vKZ3v+9c9zUulUkokEr1LMpnMZpcAgIgJFJg77rhDI0eO7He2cvHixX5nNf+xadMmpdPp3qWjoyP7aQEAkREoMMXFxZo/f75aWlr6vN/S0qKFCxdec5tYLKZ4PN5nAQAUvkD3YCSpvr5eq1evVkVFhSorK9XY2Kj29nbV1tZazAcAiKjAgXn44Yf16aef6ic/+YnOnz+vWbNm6U9/+pMmT55sMR8AIKIC/x7MUPF7MADw//F7MAAABERgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgInAT7TMlUweVjPcRPXBQDwoDUAYOIMBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYCJwYA4dOqQVK1aorKxMjuNo9+7dBmMBAKIucGC6uro0d+5cvfDCCxbzAAAKRFHQDaqrq1VdXW0xCwCggAQOTFCe58nzvN7Xruta7xIAMAyY3+RPpVJKJBK9SzKZtN4lAGAYMA/Mpk2blE6ne5eOjg7rXQIAhgHzS2SxWEyxWMx6NwCAYYbfgwEAmAh8BnPlyhWdOXOm9/W5c+d0/PhxjRkzRpMmTcrpcACA6HJ83/eDbHDw4EHdf//9/d6vqanRb3/720G3d11XiURC6XRa8Xg8yK5D5zhO2CNkJeC/YgB5FNXjSibH8MBnMF/72tc4YAEABsU9GACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGAi8PNgbmRRfQ5OVB9oJEX3aw6AMxgAgBECAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJgIFJpVK6d5771VJSYlKS0u1cuVKnTp1ymo2AECEBQpMa2ur6urqdOTIEbW0tOjzzz9XVVWVurq6rOYDAESU4/u+n+3G//73v1VaWqrW1lZ99atfzWgb13WVSCSUTqcVj8ez3TUCcBwn7BGyNoRvTyASovrzmckxvGioO5CkMWPGXHcdz/PkeV7va9d1h7JLAEBEZH2T3/d91dfXa/HixZo1a9Z110ulUkokEr1LMpnMdpcAgAjJ+hJZXV2d/vjHP+rtt9/WxIkTr7vetc5gkskkl8jyKKqn4BKXyFD4ovrzaXaJbN26ddqzZ48OHTo0YFwkKRaLKRaLZbMbAECEBQqM7/tat26ddu3apYMHD6q8vNxqLgBAxAUKTF1dnXbu3KnXXntNJSUlunDhgiQpkUjo5ptvNhkQABBNge7BXO9aYVNTk9asWZPRn8HHlPMvqtd4Je7BoPBF9ecz5/dg+GEHAGSKv4sMAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATgR44hmiK8oPiovq0vyh/zYFc4QwGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMBApMQ0OD5syZo3g8rng8rsrKSu3du9dqNgBAhAUKzMSJE7V161YdPXpUR48e1QMPPKCHHnpIJ0+etJoPABBRju/7/lD+gDFjxuhnP/uZvv/972e0vuu6SiQSSqfTisfjQ9k1bgCO44Q9QlaG+GOFG0hUv8czOYYXZfuHd3d36/e//726urpUWVl53fU8z5Pneb2vXdfNdpcAgAgJfJP/xIkTuvXWWxWLxVRbW6tdu3Zp5syZ110/lUopkUj0LslkckgDAwCiIfAlsqtXr6q9vV2XL1/WK6+8opdeekmtra3Xjcy1zmCSySSXyJCRqF4+4BIZMhXV7/FMjuFDvgezbNkyTZ06Vb/+9a8zWp97MAgiqj98BAaZiur3eCbH8CH/Hozv+33OUAAAkALe5H/66adVXV2tZDKpzs5ONTc36+DBg9q3b5/VfACAiAoUmE8++USrV6/W+fPnlUgkNGfOHO3bt08PPvig1XwAgIga8j2YoLgHgyCien2aezDIVFS/x/NyDwYAgGshMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMBEoCdaAvkW1Qd3RfUhUlH9emN44gwGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMDCkwqVRKjuNo48aNORoHAFAosg5MW1ubGhsbNWfOnFzOAwAoEFkF5sqVK3r00Ue1bds23X777bmeCQBQALIKTF1dnZYvX65ly5YNuq7neXJdt88CACh8RUE3aG5u1jvvvKO2traM1k+lUnruuecCDwYAiLZAZzAdHR3asGGDduzYoVGjRmW0zaZNm5ROp3uXjo6OrAYFAESL4/u+n+nKu3fv1je+8Q2NHDmy973u7m45jqMRI0bI87w+/9+1uK6rRCKhdDqteDye/eTAMOY4TtgjZCXA4QA5EtXvlUyO4YEukS1dulQnTpzo8973vvc9TZ8+XU899dSgcQEA3DgCBaakpESzZs3q897o0aM1duzYfu8DAG5s/CY/AMBE4E+R/a+DBw/mYAwAQKHhDAYAYILAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABNDfuAYgP583w97hKw4jhP2CFmJ6te70HEGAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJgIFZvPmzXIcp89y5513Ws0GAIiwoqAb3HXXXTpw4EDv65EjR+Z0IABAYQgcmKKiIs5aAACDCnwP5vTp0yorK1N5ebkeeeQRnT171mIuAEDEBTqDWbBggbZv364vfvGL+uSTT/T8889r4cKFOnnypMaOHXvNbTzPk+d5va9d1x3axACASHB83/ez3birq0tTp07Vj370I9XX119znc2bN+u5557r9346nVY8Hs921wAMOI4T9ghZGcJhLHRR/Zpncgwf0seUR48erdmzZ+v06dPXXWfTpk1Kp9O9S0dHx1B2CQCIiMA3+f+b53l67733tGTJkuuuE4vFFIvFhrIbAEAEBTqDefLJJ9Xa2qpz587pr3/9q771rW/JdV3V1NRYzQcAiKhAZzAffvihvvOd7+jSpUsaN26cvvKVr+jIkSOaPHmy1XwAgIgKFJjm5marOQAABYa/iwwAYILAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABOBHjgGoLD5vh/2CFlxHCfsEXANnMEAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMBE4MB999JEee+wxjR07VrfccovuvvtuHTt2zGI2AECEFQVZ+bPPPtOiRYt0//33a+/evSotLdW//vUv3XbbbUbjAQCiKlBgfvrTnyqZTKqpqan3vSlTpuR6JgBAAQh0iWzPnj2qqKjQqlWrVFpaqnnz5mnbtm0DbuN5nlzX7bMAAApfoMCcPXtWDQ0NmjZtmvbv36/a2lqtX79e27dvv+42qVRKiUSid0kmk0MeGgAw/Dm+7/uZrlxcXKyKigodPny4973169erra1Nf/nLX665jed58jyv97Xrukomk0qn04rH40MYHQD+H8dxwh7hhpPJMTzQGcyECRM0c+bMPu/NmDFD7e3t190mFospHo/3WQAAhS9QYBYtWqRTp071ee/999/X5MmTczoUACD6AgXmiSee0JEjR7RlyxadOXNGO3fuVGNjo+rq6qzmAwBEVKB7MJL0hz/8QZs2bdLp06dVXl6u+vp6/fCHP8x4e9d1lUgkuAcDIGe4B5N/mRzDAwdmqAgMgFwjMPmX85v8AABkisAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAE0Vh7TiRSIS1awBAHnAGAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJgIFZsqUKXIcp99SV1dnNR8AIKKKgqzc1tam7u7u3tf/+Mc/9OCDD2rVqlU5HwwAEG2BAjNu3Lg+r7du3aqpU6fqvvvuy+lQAIDoCxSY/3b16lXt2LFD9fX1chznuut5nifP83pfu66b7S4BABGS9U3+3bt36/Lly1qzZs2A66VSKSUSid4lmUxmu0sAQIQ4vu/72Wz49a9/XcXFxXr99dcHXO9aZzBEBgCiLZ1OKx6PD7hOVpfIPvjgAx04cECvvvrqoOvGYjHFYrFsdgMAiLCsLpE1NTWptLRUy5cvz/U8AIACETgwPT09ampqUk1NjYqKsv6MAACgwAUOzIEDB9Te3q61a9dazAMAKBBZ3+TPluu6SiQS+dwlACDHMrnJz99FBgAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEzkPTB5fvwMAMBAJsfyvAems7Mz37sEAORYJsfyvD/RsqenRx9//LFKSkrkOE5O/2zXdZVMJtXR0THok9aGE+bOL+bOv6jOztz9+b6vzs5OlZWVacSIgc9RinK65wyMGDFCEydONN1HPB6P1DfDfzB3fjF3/kV1dubuK9PH3nOTHwBggsAAAEwUVGBisZieffZZxWKxsEcJhLnzi7nzL6qzM/fQ5P0mPwDgxlBQZzAAgOGDwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABP/B5hPuNS0SVcsAAAAAElFTkSuQmCC",
      "text/plain": [
       "Figure(PyObject <Figure size 640x480 with 1 Axes>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "PyObject <matplotlib.image.AxesImage object at 0x7ff080627a60>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "spy(AQ[p, :])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.5.1",
   "language": "julia",
   "name": "julia-1.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}