<commit_before>{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pyopencl as cl\n",
    "import pyopencl.array\n",
    "import pyopencl.clrandom\n",
    "import loopy as lp\n",
    "import loopy.transform.iname"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "tol = 1e-10\n",
    "mnorm = lambda x: np.max(np.abs(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ctx = cl.create_some_context()\n",
    "queue = cl.CommandQueue(ctx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n = 16*10**3\n",
    "# x_d = cl.clrandom.rand(queue, n, dtype=np.double)\n",
    "# A_d = cl.clrandom.rand(queue, (n, n), dtype=np.double)\n",
    "\n",
    "x_d = cl.array.arange(queue, n, dtype=np.double)/n\n",
    "A_d = cl.array.arange(queue, n*n, dtype=np.double).reshape((n,n))/(n*n)\n",
    "x = x_d.get()\n",
    "A = A_d.get()\n",
    "b = A.dot(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tiled Kernel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from loopy.transform.data import reduction_arg_to_subst_rule"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "knl = lp.make_kernel(\n",
    "        \"{ [i,j]: 0<=i,j<n }\",\n",
    "        \"out[i] = sum(j, A[i,j]*x[j])\",\n",
    "         assumptions=\"n>=0\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Andreas's transform\n",
    "\n",
    "tknl = knl\n",
    "tile_size = 16\n",
    "tknl = lp.split_iname(tknl, \"i\", tile_size, inner_tag=\"ilp.seq\")\n",
    "tknl = lp.split_iname(tknl, \"j\", tile_size)\n",
    "\n",
    "# Slight issue: It won't  let you pick an order between i_inner and j_inner.\n",
    "# https://github.com/inducer/loopy/issues/82\n",
    "\n",
    "tknl = lp.add_and_infer_dtypes(tknl, {\"A,x\": np.float64})\n",
    "print(lp.generate_code_v2(tknl).device_code())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# # Scott's old transform\n",
    "# tile_size = 16\n",
    "# knl = lp.split_iname(knl, \"j\", tile_size)\n",
    "\n",
    "# knl = lp.split_reduction_inward(knl, \"j_inner\", within=\"iname:i\")\n",
    "# #knl = lp.split_reduction_inward(knl, \"j_inner\")\n",
    "# knl = lp.split_iname(knl, \"i\", tile_size)\n",
    "# knl = reduction_arg_to_subst_rule(knl, \"j_outer\")\n",
    "# knl = lp.precompute(knl, \"red_j_outer_arg\", \"j_outer,i_inner\",\n",
    "#                     temporary_scope=lp.temp_var_scope.PRIVATE)\n",
    "# knl = lp.realize_reduction(knl)\n",
    "# knl = lp.rename_iname(knl, \"red_j_outer_arg_i_inner\", \"i_inner_\")\n",
    "# knl = lp.set_loop_priority(knl, \"j_outer_0,i_inner_\")\n",
    "\n",
    "# knl = lp.assume(knl, \"n mod %i = 0\"%(tile_size,))\n",
    "\n",
    "# print(knl)\n",
    "\n",
    "# knl_tile = knl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# cgr = lp.generate_code_v2(knl)\n",
    "# print(cgr.device_code())\n",
    "# evt, (b_d,) = knl_tile(queue, A=A_d, x=x_d); evt.wait()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#%timeit evt, _ = knl_tile(queue, A=A_d, x=x_d, n=n); evt.wait()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true<commit_msg>Timing the matvec against numpy.<commit_after>{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pyopencl as cl\n",
    "import pyopencl.array\n",
    "import pyopencl.clrandom\n",
    "import loopy as lp\n",
    "import loopy.transform.iname"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "tol = 1e-10\n",
    "mnorm = lambda x: np.max(np.abs(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ctx = cl.create_some_context()\n",
    "queue = cl.CommandQueue(ctx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n = 16*10**3\n",
    "# x_d = cl.clrandom.rand(queue, n, dtype=np.double)\n",
    "# A_d = cl.clrandom.rand(queue, (n, n), dtype=np.double)\n",
    "\n",
    "x_d = cl.array.arange(queue, n, dtype=np.double)/n\n",
    "A_d = cl.array.arange(queue, n*n, dtype=np.double).reshape((n,n))/(n*n)\n",
    "x = x_d.get()\n",
    "A = A_d.get()\n",
    "b = A.dot(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tiled Kernel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from loopy.transform.data import reduction_arg_to_subst_rule"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "knl = lp.make_kernel(\n",
    "        \"{ [i,j]: 0<=i,j<n }\",\n",
    "        \"out[i] = sum(j, A[i,j]*x[j])\",\n",
    "         assumptions=\"n>=0\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Andreas's transform\n",
    "\n",
    "tknl = knl\n",
    "tile_size = 16\n",
    "tknl = lp.split_iname(tknl, \"i\", tile_size, \n",
    "                      outer_tag=\"g.0\",\n",
    "                      inner_tag=\"ilp.seq\")\n",
    "tknl = lp.split_iname(tknl, \"j\", tile_size)\n",
    "\n",
    "# Slight issue: It won't  let you pick an order between i_inner and j_inner.\n",
    "# https://github.com/inducer/loopy/issues/82\n",
    "\n",
    "tknl = lp.add_and_infer_dtypes(tknl, {\"A,x\": np.float64})\n",
    "print(lp.generate_code_v2(tknl).device_code())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "evt, (b_d,) = tknl(queue, A=A_d, x=x_d); evt.wait()\n",
    "mnorm(b-b_d.get())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%timeit evt, _ = tknl(queue, A=A_d, x=x_d, n=n); evt.wait()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%timeit A.dot(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# # Scott's old transform\n",
    "# tile_size = 16\n",
    "# knl = lp.split_iname(knl, \"j\", tile_size)\n",
    "\n",
    "# knl = lp.split_reduction_inward(knl, \"j_inner\", within=\"iname:i\")\n",
    "# #knl = lp.split_reduction_inward(knl, \"j_inner\")\n",
    "# knl = lp.split_iname(knl, \"i\", tile_size)\n",
    "# knl = reduction_arg_to_subst_rule(knl, \"j_outer\")\n",
    "# knl = lp.precompute(knl, \"red_j_outer_arg\", \"j_outer,i_inner\",\n",
    "#                     temporary_scope=lp.temp_var_scope.PRIVATE)\n",
    "# knl = lp.realize_reduction(knl)\n",
    "# knl = lp.rename_iname(knl, \"red_j_outer_arg_i_inner\", \"i_inner_\")\n",
    "# knl = lp.set_loop_priority(knl, \"j_outer_0,i_inner_\")\n",
    "\n",
    "# knl = lp.assume(knl, \"n mod %i = 0\"%(tile_size,))\n",
    "\n",
    "# print(knl)\n",
    "\n",
    "# knl_tile = knl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# cgr = lp.generate_code_v2(knl)\n",
    "# print(cgr.device_code())\n",
    "# evt, (b_d,) = knl_tile(queue, A=A_d, x=x_d); evt.wait()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
{'fraction_non_alphanumeric': 0.23300337213096922, 'fraction_numerical': 0.004895028826280866, 'mean_word_length': 2.5064836003051107, 'pattern_counts': {'":': 162, '<': 7, '<?xml version=': 0, '>': 5, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2448668', 'n_tokens_mistral': 3922, 'n_tokens_neox': 3907, 'n_words': 604}