lecture1/the_dream.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nprocs()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "1+1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Distributed, LinearAlgebra\n",
    "#rmprocs(workers())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A Brief Feel for the Dream"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[32m\u001b[1m    Updating\u001b[22m\u001b[39m registry at `~/data/.julia/registries/JuliaComputingRegistry.toml`\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[32m\u001b[1m   Resolving\u001b[22m\u001b[39m package versions...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[32m\u001b[1m  No Changes\u001b[22m\u001b[39m to `/mnt/data/.julia/environments/v1.7/Project.toml`\n",
      "\u001b[32m\u001b[1m  No Changes\u001b[22m\u001b[39m to `/mnt/data/.julia/environments/v1.7/Manifest.toml`\n"
     ]
    }
   ],
   "source": [
    "import Pkg; Pkg.add(\"CUDA\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Distributed\n",
    "\n",
    "if nprocs() == 1\n",
    "  addprocs(4) \n",
    "end\n",
    "\n",
    "# @everywhere begin\n",
    "#     import Pkg\n",
    "#     Pkg.activate($(Base.current_project()))\n",
    "#   end\n",
    "\n",
    "@everywhere using LinearAlgebra, DistributedArrays, CUDA\n",
    "@everywhere import CUDA: @time\n",
    "@everywhere CUDA.allowscalar(true)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## So you like linear algebra? <br>\n",
    "<table style=\"width:100%\">\n",
    "  <tr>\n",
    "    <th><img src=\"http://freecomputerbooks.com/covers/LAPACK-Users-Guide-Third-Edition.jpg\" width=\"200\"></th>\n",
    "    <th><img src=\"http://www.netlib.org/scalapack/slug/slug_cover.gif\"  width=\"200\" style=\"background-color:white\"></th>\n",
    "      <th><img src=\"https://images-na.ssl-images-amazon.com/images/I/41uT31fSxBL._SX352_BO1,204,203,200_.jpg\" width=\"200\"></th>\n",
    "    \n",
    "  </tr>\n",
    "    \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## but you are innovating beyond the traditional libraries?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  For demonstration let's make a matrix where you know more than the libraries know."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's create a custom  matrix structure that depends on n parameters, not $n^2$: <br> <br>\n",
    "Diagonal(v) + v*v'\n",
    "<br>\n",
    "\n",
    "For example (n=3):\n",
    "$$ \\begin{pmatrix}\n",
    "v_1 & & \\\\\n",
    "& v_2 & \\\\\n",
    "& & v_3\n",
    "\\end{pmatrix} + \n",
    "\\begin{pmatrix}\n",
    "v_1^2 & v_1 v_2 & v_1 v_3 \\\\\n",
    "v_2 v_1 & v_2^2 & v_2 v_3 \\\\\n",
    "v_3 v_1 & v_3 v_2 & v_3^2\n",
    "\\end{pmatrix} .\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Build a Custom Type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "struct MyMatrix{T, V<:AbstractVector{T}} <: AbstractMatrix{T}\n",
    "   v::V\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "Base.size(A::MyMatrix)         =   length(A.v), length(A.v)\n",
    "Base.getindex(A::MyMatrix,i,j) =   A.v[i]*(i==j)  + A.v[i]*A.v[j]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 MyMatrix{Int64, Vector{Int64}}:\n",
       "   2    10    100\n",
       "  10   110   1000\n",
       " 100  1000  10100"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "A = MyMatrix([1,10,100])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MyMatrix{Int64, Vector{Int64}}\n",
      "  v: Array{"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Int64}((3,)) [1, 10, 100]\n"
     ]
    }
   ],
   "source": [
    "dump( A )  # n storage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Array{Int64}((3, 3)) "
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2 10 100; 10 110 1000; 100 1000 10100]\n"
     ]
    }
   ],
   "source": [
    "dump(Matrix(A)) # n² storage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## My very own largest eigensolver for my very own matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "f(A::MyMatrix)  = λ ->  1 + mapreduce((v) -> v^2 / (v - λ)  , +, A.v)\n",
    "f′(A::MyMatrix) = λ ->      mapreduce((v) -> v^2 / (v - λ)^2, +, A.v)\n",
    "\n",
    "function LinearAlgebra.eigmax(A::MyMatrix; tol = eps(2.0), debug = false)\n",
    "    x0 = maximum(A.v) + maximum(A.v)^2\n",
    "    δ = f(A)(x0)/f′(A)(x0)\n",
    "    while abs(δ) > x0 * tol               \n",
    "        x0 -= δ\n",
    "        δ = f(A)(x0)/f′(A)(x0) \n",
    "        debug && println(\"x = $x0, δ = $δ\") # Debugging\n",
    "    end\n",
    "    x0\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10200.107083707298"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "eigmax(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10200.107083707298"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "eigmax(Matrix(A))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Go Heterogeneous"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 MyMatrix{Int64, CuArray{Int64, 1, CUDA.Mem.DeviceBuffer}}:\n",
       " 2  2   3\n",
       " 2  6   6\n",
       " 3  6  12"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "gpuA = MyMatrix(CuArray([1,2,3]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 MyMatrix{Int64, DArray{Int64, 1, Vector{Int64}}}:\n",
       " 2  2   3\n",
       " 2  6   6\n",
       " 3  6  12"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "distA = MyMatrix(distribute([1,2,3]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Int64}:\n",
       " 2  2   3\n",
       " 2  6   6\n",
       " 3  6  12"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Matrix(gpuA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Int64}:\n",
       " 2  2   3\n",
       " 2  6   6\n",
       " 3  6  12"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Matrix(distA)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compare Timings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4000000"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = 4_000_000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4000000-element Vector{Float64}:\n",
       " -0.053793817582657326\n",
       " -0.003592096969330304\n",
       " -0.18867427852644192\n",
       " -0.01582632633108473\n",
       " -0.1832277013847907\n",
       "  0.005077106926582318\n",
       "  0.09429267789201608\n",
       " -0.09386769921413878\n",
       " -0.031091414723704663\n",
       "  0.12403022622093754\n",
       "  ⋮\n",
       "  0.08202188137552724\n",
       "  0.2054528132742123\n",
       "  0.1863330435843987\n",
       " -0.04871046896917302\n",
       " -0.23073206053276693\n",
       " -0.10850741820814912\n",
       " -0.04602062235592608\n",
       " -0.08199493813278634\n",
       " -0.21094130703715375"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "v = randn(N)*.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MyMatrix{Float64, Vector{Float64}}\n",
      "  v: Array{"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Float64}((4000000,)) [-0.053793817582657326, -0.003592096969330304, -0.18867427852644192, -0.01582632633108473, -0.1832277013847907, 0.005077106926582318, 0.09429267789201608, -0.09386769921413878, -0.031091414723704663, 0.12403022622093754  …  -0.052266126744157315, 0.08202188137552724, 0.2054528132742123, 0.1863330435843987, -0.04871046896917302, -0.23073206053276693, -0.10850741820814912, -0.04602062235592608, -0.08199493813278634, -0.21094130703715375]\n"
     ]
    }
   ],
   "source": [
    "dump(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4000000×4000000 MyMatrix{Float64, Vector{Float64}}:\n",
       " -0.0509        0.000193233   0.0101495    …   0.00441082    0.0113473\n",
       "  0.000193233  -0.00357919    0.000677736      0.000294534   0.000757722\n",
       "  0.0101495     0.000677736  -0.153076         0.0154703     0.0397992\n",
       "  0.000851359   5.68497e-5    0.00298602       0.00129768    0.00333843\n",
       "  0.00985652    0.000658172   0.0345704        0.0150237     0.0386503\n",
       " -0.000273117  -1.82375e-5   -0.000957919  …  -0.000416297  -0.00107097\n",
       " -0.00507236   -0.000338708  -0.0177906       -0.00773152   -0.0198902\n",
       "  0.0050495     0.000337182   0.0177104        0.00769668    0.0198006\n",
       "  0.00167253    0.000111683   0.00586615       0.00254934    0.00655846\n",
       " -0.00667206   -0.000445529  -0.0234013       -0.0101699    -0.0261631\n",
       "  ⋮                                        ⋱                \n",
       " -0.00441227   -0.000294631  -0.0154754       -0.00672538   -0.0173018\n",
       " -0.0110521    -0.000738006  -0.0387637       -0.0168461    -0.0433385\n",
       " -0.0100236    -0.000669326  -0.0351563       -0.0152784    -0.0393053\n",
       "  0.00262032    0.000174973   0.00919041       0.00399401    0.010275\n",
       "  0.012412      0.000828812   0.0435332    …   0.0189189     0.0486709\n",
       "  0.00583703    0.000389769   0.0204726        0.00889706    0.0228887\n",
       "  0.00247562    0.000165311   0.00868291       0.00377346    0.00970765\n",
       "  0.00441082    0.000294534   0.0154703       -0.0752718     0.0172961\n",
       "  0.0113473     0.000757722   0.0397992        0.0172961    -0.166445"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "A = MyMatrix(v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "LoadError: UndefVarError: @folded_threads not defined\nin expression starting at /home/jrun/data/code/sc.ipynb:1",
     "output_type": "error",
     "traceback": [
      "LoadError: UndefVarError: @folded_threads not defined\n",
      "in expression starting at /home/jrun/data/code/sc.ipynb:1\n"
     ]
    }
   ],
   "source": [
    "@folded_threads eigmax(A) # using a macro\n",
    "\n",
    "folded_context(true)  # context variable\n",
    "eigmax(A)\n",
    "\n",
    "eigmax( Threaded(A) ) # similar to CuArray and distribute (using the type mechanism to figure out to execute)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "distA = MyMatrix(distribute(v));  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "gpuA = MyMatrix(CuArray(v));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "1+1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  3.472926 seconds (5.08 M CPU allocations: "
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "272.595 MiB, 2.29% gc time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "39991.435006823696"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@time eigmax(distA) # run twice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(value = 39991.4350068237, time = 0.004795719, cpu_bytes = 203792, cpu_gctime = 0.0, cpu_gcstats = Base.GC_Diff(203792, 0, 0, 4143, 0, 0, 0, 0, 0), gpu_bytes = 59248, gpu_memtime = 0.0005168400000000739, gpu_memstats = CUDA.AllocStats(92, 59248, 0, 0, 0.0005168400000000739))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "CUDA.@timed(eigmax(gpuA)) # run twice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.004720 seconds (4.14 k CPU allocations: 199.016 KiB) (92 GPU allocations: 57.859 KiB, 11.30% memmgmt time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "39991.4350068237"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@time eigmax(gpuA)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Abstraction\n",
    "\n",
    "1. A data structure is a mathematical matrix! <br>\n",
    "2. A matrix can be serial, distributed, gpu <br>\n",
    "BOTH ARE ABSTRACTIONS, BOTH use underlying similar mechanisms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Why Julia?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Well designed abstractions <br>\n",
    "2. Multiple dispatch <br>\n",
    "3. Careful balance between static and dynamic <br>\n",
    "4. Compiles to GPU at multiple levels, not just lowest CUDA level <br>\n",
    "5. Metaprogramming across the stack <br>\n",
    "6. Interfaces with LLVM <br>\n",
    "7. Plays nicely with Python, legacy codes, ... <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> <a href=\"https://arxiv.org/abs/1411.1607\">\n",
    "Julia: A Fresh Approach to Numerical Computing \n",
    "</a> </center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " <center> <img src=\"siam.png\" width=400> </a> </center>"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.7.3",
   "language": "julia",
   "name": "julia-1.7"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.7.3"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}