Lec06.LoopOptimizations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loop Optimization in OpenMP\n",
    "\n",
    "This lecture is almost all by example in file [stencil.c](./openmp/omp_c/stencil.c). This produces the timing results that show the performance benefits. Run this code without compiler optimization\n",
    "```\n",
    "gcc -fopenmp -O0 stencil.c \n",
    "clang -fopenmp -O0 stencil.c\n",
    "```\n",
    "\n",
    "This makes the results more understandable because it prevents the compiler from vectorizing the code and optimizing loops. The compiler will automatically make some of the optimizations that we want to demonstrate.\n",
    "\n",
    "### Loop iteration order\n",
    "\n",
    "For 2-d dense data, the array must be serialized to memory, i.e. in a linear order.\n",
    "The serialization strategies are named by which dimension (row versus column) \n",
    "occurs sequentially in memory.\n",
    "\n",
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/4/4d/Row_and_column_major_order.svg\" width=256 title=\"Row versus column major order.\" />\n",
    "\n",
    "Choosing a memory efficient order for loops has a big impact on performance.\n",
    "  * Successive loop iterations access adjacent elements or\n",
    "  * Successive loop iterations access strided elements\n",
    "  \n",
    "The different orders are also associated with programming languages that use these conventions.\n",
    "  \n",
    "<img src=\"https://images.slideplayer.com/23/6540072/slides/slide_3.jpg\" width=512 title=\"from Edgar Gabriel at UH\" />\n",
    "\n",
    "There are many conventions about loop ordering and they get confusing.  Reason carefully about how the loops variables are enumerated and the data layout.  For example, images are almost always in Fortran order so that programming then in C looks weird."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sequential access in `stencil.c`\n",
    "\n",
    "We provide two routines that show the difference between sequential and strided access in C.\n",
    "\n",
    "Which of the following performs sequential access?\n",
    "\n",
    "```c\n",
    "void initializeyx ( double* array )\n",
    "{\n",
    "    /* Initialize the array to random values */\n",
    "    for (int y=0; y<DIM; y++) {\n",
    "        for (int x=0; x<DIM; x++) {\n",
    "            array[x*DIM+y] = (double)rand()/RAND_MAX;\n",
    "        }        \n",
    "    }\n",
    "}\n",
    "\n",
    "void initializexy ( double* array )\n",
    "{\n",
    "    /* Initialize the array to random values */\n",
    "    for (int x=0; x<DIM; x++) {\n",
    "        for (int y=0; y<DIM; y++) {\n",
    "            array[x*DIM+y] = (double)rand()/RAND_MAX;\n",
    "        }        \n",
    "    }\n",
    "}\n",
    "```\n",
    "\n",
    "The performance difference reflects the latency difference between sequential and strided acceess. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A Parallel Stencil\n",
    "\n",
    "A common pattern in numerical computing is to compute a [compact stencil](https://en.wikipedia.org/wiki/Compact_stencil). \n",
    "\n",
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/thumb/c/c5/CompactStencil.svg/300px-CompactStencil.svg.png\" width=256 title=\"Compact stencil.\" />\n",
    "\n",
    "The following function computes an average over a compact stencil at each (well defined) cell in a 2-d grid.  This computation pattern is used frequently in convolutional neural networks, graphics, spatial data processing, etc.\n",
    "\n",
    "```c\n",
    "void stencil_average ( double* input_ar, double* output_ar )\n",
    "{\n",
    "    double partial = 0.0;\n",
    "\n",
    "    for (int x=HWIDTH; x<DIM-HWIDTH; x++) {\n",
    "        for (int y=HWIDTH; y<DIM-HWIDTH; y++) {\n",
    "            for (int xs=-1*HWIDTH; xs<=HWIDTH; xs++) {\n",
    "                for (int ys=-1*HWIDTH; ys<=HWIDTH; ys++) {\n",
    "                    partial += input_ar[DIM*(x+xs)+(y+ys)];\n",
    "                }   \n",
    "            }   \n",
    "            output_ar[DIM*x+y] = partial/((2*HWIDTH+1)*(2*HWIDTH+1));\n",
    "            partial=0.0;\n",
    "        }       \n",
    "    }\n",
    "}\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parallelizing `stencil_average`\n",
    "\n",
    "```c\n",
    "void stencil_average_omp ( double* input_ar, double* output_ar )\n",
    "{\n",
    "    omp_set_num_threads(4);\n",
    "    #pragma omp parallel for \n",
    "    for (int x=HWIDTH; x<DIM-HWIDTH; x++) {\n",
    "        for (int y=HWIDTH; y<DIM-HWIDTH; y++) {\n",
    "            double partial = 0.0;\n",
    "```\n",
    "\n",
    "To parallelize this computation, we need to do two things:\n",
    "  (1) add parallel directives around the outer loop.\n",
    "  (2) move the variable into the inner scope.\n",
    "\n",
    "_Why the outer loop?_ This creates a groups of threads that divide the iterations of the outer most loop. The parallel context exists for the entire computation.  \n",
    "\n",
    "_What happens if we parallelize the inner loop?_\n",
    "\n",
    "```c\n",
    "    // This is wrong\n",
    "    for (int x=HWIDTH; x<DIM-HWIDTH; x++) {\n",
    "        #pragma omp parallel for \n",
    "        for (int y=HWIDTH; y<DIM-HWIDTH; y++) {\n",
    "            double partial = 0.0;\n",
    "```\n",
    "\n",
    "Threads are created and destroyed for each iteration of the outer loop. The overhead of thread creation will slow down the program and loose speedup.  Refer to the function `stencil_avg_omp_inner()`.\n",
    "\n",
    "### loop independence \n",
    "\n",
    "For parallelism, we require that the iterations of the loops are _independent_, i.e. the computation of that iteration does not depend on prior iterations and does not affect subsequent iterations. This allows a loop to be run in parallel and produce the same result as the serial program.  We say that _loop independence_ is neccessary to achieve _serial equivalance_.\n",
    "\n",
    "The code also moves the declaration and initialization of variable `partial` inside the loop. This is needed for correctness. The variable `partial` sums the contribution of the stencil for each iteration. By declaring partial inside the loop, each iteration has its own private variable. This also means that each thread has it's own copy of partial.  We say that this variable is _thread private_.\n",
    "\n",
    "If we leave partial defined in the outer scope, it is a \"shared variable\" among threads (see `stencil_average_omp_bad`). The parallel threads update a single copy of the variable.  In addition to being incorrect, this has a negative performance affect, because the shared variable leads to _interference_ in the form _of memory contention_.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loop Unrolling\n",
    "\n",
    "Loop unrolling is a time-space tradeoff typically made by compilers\n",
    "  * time savings: eliminate branching instructions in evaluating loop conditional\n",
    "  * space increase: make a bigger program with more statements\n",
    "\n",
    "This example unrolls the entire stencil (5x5) eliminating the two inner loops.\n",
    "    \n",
    "```c\n",
    "            partial = input_ar[DIM*(x-2)+(y-2)];\n",
    "            partial += input_ar[DIM*(x-2)+(y-1)];\n",
    "            partial += input_ar[DIM*(x-2)+(y)];\n",
    "            partial += input_ar[DIM*(x-2)+(y+1)];\n",
    "            partial += input_ar[DIM*(x-2)+(y+2)];\n",
    "\n",
    "            partial += input_ar[DIM*(x-1)+(y-2)];\n",
    "            partial += input_ar[DIM*(x-1)+(y-1)];\n",
    "            partial += input_ar[DIM*(x-1)+(y)];\n",
    "            partial += input_ar[DIM*(x-1)+(y+1)];\n",
    "            partial += input_ar[DIM*(x-1)+(y+2)];\n",
    "\n",
    "            partial += input_ar[DIM*(x)+(y-2)];\n",
    "            partial += input_ar[DIM*(x)+(y-1)];\n",
    "            partial += input_ar[DIM*(x)+(y)];\n",
    "            partial += input_ar[DIM*(x)+(y+1)];\n",
    "            partial += input_ar[DIM*(x)+(y+2)];\n",
    "\n",
    "            partial += input_ar[DIM*(x+1)+(y-2)];\n",
    "            partial += input_ar[DIM*(x+1)+(y-1)];\n",
    "            partial += input_ar[DIM*(x+1)+(y)];\n",
    "            partial += input_ar[DIM*(x+1)+(y+1)];\n",
    "            partial += input_ar[DIM*(x+1)+(y+2)];\n",
    "\n",
    "            partial += input_ar[DIM*(x+2)+(y-2)];\n",
    "            partial += input_ar[DIM*(x+2)+(y-1)];\n",
    "            partial += input_ar[DIM*(x+2)+(y)];\n",
    "            partial += input_ar[DIM*(x+2)+(y+1)];\n",
    "            partial += input_ar[DIM*(x+2)+(y+2)];\n",
    "\n",
    "            output_ar[DIM*x+y] = partial/((2*HWIDTH+1)*(2*HWIDTH+1));\n",
    "            partial = 0.0;\n",
    "```\n",
    "\n",
    "The example code shows:\n",
    "  * unrolling improves serial performance (`stencil_average_unrolled()`)\n",
    "  * unrolling improves parallel performance (`stencil_average_omp_unrolled()`)\n",
    "  \n",
    "In both cases the benefit comes from reducing the number of instructions.\n",
    "  * _What instructions are eliminated?_\n",
    "  * _Approximately what fraction of instructions are eliminated_?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We say that this loop is _fully unrolled_ in that we have written all of its iterations sequentially. Loops can be partially unrolled. It's common to refer to a loop being unrolled _X times_ which means that X iterations of the original loop happen in each iteration of the unrolled loop. \n",
    "\n",
    "A loop that has been unrolled 4 times looks like.\n",
    "\n",
    "```c\n",
    "// original loop\n",
    "for (i=0; i<n; i++)\n",
    "{\n",
    "    do_stuff(i)\n",
    "}\n",
    "\n",
    "// unrolled four times\n",
    "for (i=0; i<n; i=i+4)\n",
    "{\n",
    "    do_stuff(i)\n",
    "    do_stuff(i+1)\n",
    "    do_stuff(i+2)\n",
    "    do_stuff(i+3)    \n",
    "}\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loop Fusion\n",
    "\n",
    "Another effective optimization. The concept is to do the work of multiple loops in a single loop. For serial code, this has the benefit:\n",
    "  * evaluate loop conditional variables once for all fused loops\n",
    "  \n",
    "In OpenMP, this has the additional benefit:\n",
    "  * create and destroy threads once for the fused loops, rather than in each loop\n",
    "\n",
    "The code example implements a function that sums two arrays in parallel:\n",
    "\n",
    "```c\n",
    "    #pragma omp parallel for \n",
    "    for (int x=0; x<DIM; x++) {\n",
    "        for (int y=0; y<DIM; y++) {\n",
    "            output_ar[x*DIM+y] = input_ar1[x*DIM+y] + input_ar2[x*DIM+y];\n",
    "        }        \n",
    "    }\n",
    "```\n",
    "\n",
    "Now consider that we want to compute a stencil average on two arrays and then add the result.  This can be done with three separate function calls each that has its own loop:\n",
    "```c\n",
    "    stencil_average_omp(rand_ar1, avg_ar1);\n",
    "    stencil_average_omp(rand_ar2, avg_ar2);\n",
    "    array_sum_omp(avg_ar1, avg_ar2, sum_ar);\n",
    "```\n",
    "\n",
    "The results for `seperate loops` shows the performance of doing each loop independently in parallel.  We can _fuse_ these loops and compute the average and sum in a single loop with one function call\n",
    "\n",
    "```c\n",
    "void fused_stencil_sum_omp ( double* input_ar1, double* input_ar2, double* output_ar )\n",
    "{\n",
    "    omp_set_num_threads(4);\n",
    "    #pragma omp parallel for \n",
    "    for (int x=HWIDTH; x<DIM-HWIDTH; x++) {\n",
    "        for (int y=HWIDTH; y<DIM-HWIDTH; y++) {\n",
    "            double partial1 = 0.0;\n",
    "            double partial2 = 0.0;\n",
    "            for (int xs=-1*HWIDTH; xs<=HWIDTH; xs++) {\n",
    "                for (int ys=-1*HWIDTH; ys<=HWIDTH; ys++) {\n",
    "                    partial1 += input_ar1[DIM*(x+xs)+(y+ys)];\n",
    "                    partial2 += input_ar2[DIM*(x+xs)+(y+ys)];\n",
    "                }   \n",
    "            }   \n",
    "            output_ar[DIM*x+y] = partial1/((2*HWIDTH+1)*(2*HWIDTH+1)) + partial1/((2*HWIDTH+1)*(2*HWIDTH+1));\n",
    "            partial1=0.0;\n",
    "            partial2=0.0;\n",
    "        }       \n",
    "    }\n",
    "}\n",
    "\n",
    "fused_stencil_sum_omp(rand_ar1, rand_ar2, sum_ar);\n",
    "```\n",
    "\n",
    "This again has performance benefits. The benefits come from at least two sources.\n",
    "What are they?\n",
    "<!--\n",
    "  * Better memory locality. We only iterate over each array once.\n",
    "  * Reduce branching statements (eliminate for loop)s\n",
    "-->\n",
    "\n",
    "\n",
    "### Loop Fission\n",
    "\n",
    "This optimization divides the work of a single loop into multiple loops. I have not built an example, because I can't find one that is natural and effective. Fission can be used to make the data references in a loop smaller so that it fits into a smaller cache."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Seperable Dependencies and Reduction\n",
    "\n",
    "It is a common pattern to compute an aggregate quantity (mean, sum, maximum) in a loop. This is known as a _reduction_ because you are reduce a larger amount of data into a single quantity. The natural implementation of this uses a shared variable. This is **inefficient**.\n",
    "\n",
    "```c\n",
    "void max_el_shared ( double* input_ar )\n",
    "{\n",
    "    double max_el = 0;\n",
    "    omp_set_num_threads(4);\n",
    "    \n",
    "    #pragma omp parallel for\n",
    "    for (int x=0; x<DIM; x++) {\n",
    "        for (int y=0; y<DIM; y++) {\n",
    "            max_el = max_el > input_ar[x*DIM+y] ? max_el : input_ar[x*DIM+y]; \n",
    "        }        \n",
    "    }\n",
    "}\n",
    "```\n",
    "\n",
    "All parallel threads are reading and writing a single variable so that memory location must be shared between all threads either through L3 (multicore, single processor) or memory (SMP or NUMA).\n",
    "\n",
    "One might observe that the dependency can be _seperated_ through the following process:\n",
    "  * give each thread a private variable `thread_max_el`\n",
    "  * compute a thread local maximum in each thread\n",
    "  * after all threads complete take the maximum of all the thread local maximums\n",
    "\n",
    "This can be done manually and requires care. Naive implementations often result in false sharing. For example, if one creates an array to hold thead local variables, the array will reside in one cache line and performance will be poor due to _false sharing_.\n",
    "\n",
    "OpenMP provides a directive for _reduction_ that takes care of all the details.\n",
    "\n",
    "```c\n",
    "   #pragma omp parallel for reduction ( max: max_el )\n",
    "   for (int x=0; x<DIM; x++) {\n",
    "        for (int y=0; y<DIM; y++) {\n",
    "            max_el = max_el > input_ar[x*DIM+y] ? max_el : input_ar[x*DIM+y]; \n",
    "        }        \n",
    "    }\n",
    "```\n",
    "`reduction` describe the pattern and the clause specifies a reduction operator and the variable name. There are:\n",
    "* Arithmetic reductions: +,*,-,max,min .\n",
    "* Logical operator reductions in C: & && | || ^\n",
    "\n",
    "All of these are seperable dependencies and OpenMP will compute the partial result in each thread and accumulate the final result after threads complete.  Again, we can see the performance benefit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### OK what is false sharing?\n",
    "\n",
    "The following image is credited to Intel and appears on https://d2l.ai/chapter_computational-performance/hardware.html\n",
    "\n",
    "<img src=\"https://d2l.ai/_images/falsesharing.svg\" width=512 />\n",
    "\n",
    "Almost all figures that describe false sharing use this format, including the one I learned in graduate school. False sharing arises when two cores/processors update _different_ addresses in the same cache line. The hardware that manages cache coherency invalidates the cache line on the other processor.  This should take latency equal to the level in the hierarchy at which the sharing occurs. Using the i7 Nehalem numbers for Lecture 05\n",
    "  * 35 clock cycles for two cores on same processor (L3)\n",
    "  * 130 clock cycles for two cores on different processors (memory)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Turn Optimization On\n",
    "\n",
    "Rebuilding and running the code with optimization reveals that compiler optimizations are likely separating dependencies and unrolling loops.\n",
    "```\n",
    "gcc -fopenmp -O3 stencil.c \n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loop Scheduling\n",
    "\n",
    "This is really an aside. I just want you to know that it exists.\n",
    "\n",
    "The full looping directive includes the specification of a scheduling directive and a chunk size\n",
    "```c\n",
    "#pragma omp parallel for schedule(kind [,chunk size])\n",
    "```\n",
    "in which schedule can be one of:\n",
    "* Static -- divide loop into equal sized chunks\n",
    "* Dynamic -- build internal work queue and dispatch blocksize at a time\n",
    "* Guided -- dynamic scheduling with decreasing block size for load balance\n",
    "* Auto -- compiler chooses from above\n",
    "* Runtime -- runtime configuration chooses from above\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}