.. title: PAR Class 9, Wed 2018-03-21
.. slug: class09
.. date: 2018-03-21
.. tags: class
.. category: 
.. link: 
.. description: 
.. type: text

.. raw:: html

   <style> .red {color:red} </style>
   <style> .blue {color:blue} </style>

.. role:: red
.. role:: blue

.. sectnum::
.. contents:: Table of contents

parallel.ecse hardware details
------------------------------

I put the invoice on parallel.ecse in /parallel-class/ .   It gives the hardware specifics.

	      

Nvidia GPU summary
------------------

Here's a summary of the Nvidia Pascal GP104 GPU architecture as I understand it.  It's more
compact than I've found elsewhere.  I'll add to it from time to time.  Some numbers are probably wrong.

#. The **host** is the CPU.

#. The **device** is the GPU.

#. The device contains 20 **streaming multiprocessors** (SMs).

   Different GPU generations have used the terms SMX or SMM.

#. A **thread** is a sequential program with private and shared memory, program counter, etc.

#. Threads are grouped, 32 at a time, into **warps**.

#. Warps of threads are grouped into **blocks**.  

   Often the warps are only implicit, and we consider that the threads are grouped directly into blocks.

   That abstract hides details that may be important; see below.

#. Blocks of threads are grouped into a **grid**, which is all the threads in the kernel.

#. A **kernel** is a parallel program executing on the device.

   a. The kernel runs potentially thousands of **threads**.

   #. A kernel can create other kernels and wait for their completion.

   #. There may be a limit, e.g., 5 seconds, on a kernel's run time.

#. Thread-level resources:

   a. Each thread can use up to 255 fast **registers**.  Registers are private to the thread.

      All the threads in one block have their registers allocated from a fixed pool of 65536 registers.  The more registers that each thread uses, the fewer warps in the block  can run simultaneously.

   #. Each thread has 512KB slow **local memory**, allocated from the global memory.

   #. Local memory is used when not enough registers are available, and to
      store thread-local arrays. 

#. Warp-level resources:

   a. Threads are grouped, 32 at a time, into **warps**.

   #. Each warp executes as a SIMD, with one instruction register.  At each cycle,
      every thread in a warp is either executing the same instruction, or is disabled.
      If the 32 threads want to execute 32 different instructions, then they will
      execute one after the other, sequentially.

      If you read in some NVidia doc that threads in a warp run independently, then
      continue reading the next page to get the info mentioned in the previous paragraph.

   #. If successive instructions in a warp do not depend on each other, then,
      if there are enough warp schedulers available, they may be executed in
      parallel.   This is called **Instruction Level Parallelism (ILP)**.

   #. For an array in local memory, which means that each thread will have
      its private copy, the elements for all the threads in a warp are
      **interleaved** to potentially increase the I/O rate.

      Therefore your program should try to have successive threads read successive
      words of arrays.

   #. A thread can read variables from other threads in the same warp, with the
      **shuffle** instruction.  Typical operation are to read from the K-th next
      thread, to do a butterfly permutation, or to do an indexed read.  This happens in
      parallel for the whole warp, and does not use shared memory.

   #. A **warp vote** combines a bit computed by each thread to report
      results like *all* or *any*.

#. Block-level resources:

   a. A block may contain up to 1024 threads.

   #. Each block has access to 65536 fast 32-bit **registers**,
      for the use of its threads.

   #. Each block can use up to 49152 bytes of the SM's fast **shared**
      **memory**.  The block's shared memory is shared by all the threads in
      the block, but is hidden from other blocks.

      Shared memory is basically a user-controllable cache of some global
      data.  The saving comes from reusing that shared data several times
      after you loaded it from global memory once.

      Shared memory is interleaved in banks so that some access patterns are faster than others.

   #. Warps in a block run asynchronously and run different instructions.  They
      are scheduled and executed as resources are available.

   #. The threads in a block can be synchonized with **__syncthreads()**.

      Because of how warps are scheduled, that can be slow.

   #. The threads in a block can be arranged into a 3D array, up to
      1024x1024x64.

      That is for convenience, and does not increase performance (I think).

   #. I'll talk about **textures** later.


#. Streaming Multiprocessor (SM) - level resources:

   a. Each SM has 128 single-precision CUDA cores, 64
      double-precision units, 32 special function units, and
      32 load/store units.    

   #. In total, the GPU has 2560 CUDA cores.

   #. A **CUDA core** is akin to an ALU.  The cores, and all the units, are
      pipelined.

   #. A CUDA core is much less powerful than one core of an Intel Xeon.  My
      guess is 1/20th.

   #. Beware that, in the CUDA C Programming Guide, NVidia sometimes calls an
      SM a core.

   #. The limited number of, e.g., double precision units means that an DP
      instruction will need to be scheduled several times for all the threads
      to execute it.  That's why DP is slower.

   #. Each SM has 4 warp schedulers and 8 instruction dispatch units.

   #. 64 warps can simultaneously reside in an SM.

   #. Therefore up to 32x64=2048 threads can be executed in parallel by an
      SM.

   #. Up to 16 blocks that can simultaneously be resident in an SM.
  
      However, if each block uses too many resources, like shared memory,
      then this number is reduced.

      Each block sits on only one SM; no block is split.  However a block's
      warps are executed asynchronously (until synced).

   #. Each SM has 64KiB (?) fast memory to be divided between **shared** memory and an **L1 cache**.  Typically, 48KiB (96?) is used for the shared memory, to be divided among its resident blocks, but that can be changed.

   #. The 48KB L1 cache can cache local or global memory.

   #. Each SM has a read-only data cache of 48KB to cache the
      global constant memory.

   #. Each SM has 8 texture units, and many other graphics capabilities.

   #. Each SM has 256KB of L2 cacha.

#. Grid-level resources:

   a. The blocks in a grid can be arranged into a 3D array.  
      up to :math:`(2^{31}-1, 2^{16}-1, 2^{16}-1)`.

   #. Blocks in a grid might run on different SMs.
      
   #. Blocks in a grid are queued and executed as resources are
      available, in an unpredictable parallel or serial order.
      Therefore they should be independent of each other.

   #. The number of instructions in a kernel is limited.

   #. Any thread can stop the kernel by calling **assert**.

#. Device-level resources:

   a. There is a large and slow 8GB **global memory**, which
      persists from kernel to kernel.

      Transactions to global memory are 128 bytes.

      Host memory can also be memory-mapped into global memory, although the
      I/O rate will be lower.

      Reading from global memory can take hundreds of cycles.  A warp that
      does this will be paused and another warp started.  Such context
      switching is very efficient.  Therefore device throughput stays high,
      although there is a latency.  This is called **Thread Level
      Parallelism (TLP)** and is a major reason for GPU performance.

      That assumes that an SM has enough active warps that there is always
      another warp available for execution.  That is a reason for having
      warps that do not use all the resources (registers etc) that they're
      allowed to.

   #. There is a 2MB L2 cache, for sharing data between SMs.
      
   #. There is a 64KiB Small and fast global **constant memory**, ,
      which also persists from kernel to kernel.  It is implemented as a
      piece of the global memory, made fast with caches.

      (Again, I'm still resolving this apparent contradiction).

   #. **Grid Management Unit (GMU)** schedules (pauses, executes, etc) grids on
      the device.  This is more important because grids can start other
      grids **(Dynamic Parallelism)**.

   #. **Hyper-Q**: 32 simultaneous CPU tasks can launch kernels into the
      queue; they don't block each other.  If one kernel is waiting, another runs.

   #. **CUDA Work Distributor (CWD)** dispatches 32 active grids at
      a time to the SMs.  There may be 1000s of grids queued and waiting.

   #. **GPU Direct**: Other devices can DMA the GPU memory.

   #. The base clock is 1607MHz.

   #. GFLOPS: 8873.

   #. Memory bandwidth: 320GB/s

#. GPU-level resources:

   a. Being a Geforce product, there are many graphics facilities that we're not using.

   #. There are 4 **Graphics processing clusters** (GPCs) to do graphics stuff.

   #. Several perspective projections can be computed in parallel, for systems with several displays.

   #. There's HW for texture processing.

#. Generational changes:

   a. With each new version, Nvidia tweaks the numbers.   Some get higher, others get lower.

      i. E.g., Maxwell had little HW for double precision, and so that was slow.

      #. Pascal's clock speed is much higher.      
      
#. Refs:

   a. The CUDA program deviceDrv.

   #. http://developer.download.nvidia.com/compute/cuda/compute-docs/cuda-performance-report.pdf

   #. http://international.download.nvidia.com/geforce-com/international/pdfs/GeForce_GTX_1080_Whitepaper_FINAL.pdf

   #. `Better Performance at Lower Occupancy <http://www.cs.berkeley.edu/~volkov/volkov10-GTC.pdf>`_, 
      Vasily Volkov, UC Berkeley, 2010.

   #. https://www.pgroup.com/lit/articles/insider/v2n1a5.htm - well written but old.

   *(I'll keep adding to this. Suggestions are welcome.)*

More CUDA
---------

#. CUDA function qualifiers:

   a. *__global__*   device function called from host, starting a kernel.

   #. *__device__* device function called from device function.

   #. *__host__* (default)  host function called from host function.

#. CUDA variable qualifiers:

   a. *__shared__*
   #. *__device__* global
   #. *__device__ __managed__* automatically paged between host and device.
   #. *__constant__*
   #. (nothing) register if scalar, or local if array or if no more registers
      available.

#. If installing CUDA on your machine, this repository seems best:

   http://developer.download.nvidia.com/compute/cuda/repos/ubuntu1604/x86_64

   That includes the Thrust headers but not example programs.

   
Thrust
------

#. Thrust is an API that looks like STL. Its backend can be CUDA,
   OpenMP, or sequential host-based code.

#. The online Thrust directory structure is a mess.  Three main
   sites appear to be these:

   a. https://github.com/thrust -

      i. The best way to install it is to clone from here.

      #. The latest version of the examples is also here.

      #. The wiki has a lot of doc.
      
   #. https://thrust.github.io/

      This points to the above site.

   #. https://developer.nvidia.com/thrust

      This has links to other Nvidia docs, some of which are obsolete.

   #. http://docs.nvidia.com/cuda/thrust/index.html

      easy-to-read, thorough, obsolete, doc

   #. https://code.google.com/  - no longer exists.

#. The latest version is 1.8.3.

#. Functional-programming philosophy.

#. Many possible backends:  host, GPU, OpenMP, TBB...

#. Easier programming, once you get used to it.

#. Code is efficient.

#. Uses some unusual C++ techniques, like overloading **operator()**.

#. Since the Stanford slides were created, Thrust has adopted
   unified addressing, so that pointers know whether they are
   host or device.

#. On parallel in /parallel-class/thrust/ are many little demo programs from the thrust distribution, with my additions.

#. CUDACast videos on Thrust:

   `CUDACast #.15 - Introduction to Thrust <https://www.youtube.com/watch?v=mZJEbO9Eros>`_

   `CUDACast #.16 - Thrust Algorithms and Custom Operators <https://www.youtube.com/watch?v=xtWJCL7LMqU>`_


#. **Thrust is fast because** the functions that look like they
   would need linear time really take only log time in parallel.

#. In functions like **reduce** and **transform**, you often see an argument like **thrust::multiplies<float>()**.  The syntax is as follows:

   a. **thrust::multiplies<float>** is a class.

   #. It overloads **operator()**.

   #. However, in the call to reduce, **thrust::multiplies<float>()** is calling the default
      constructor to construct a variable of class
      **thrust::multiplies<float>**, and passing it to **reduce**.

   #. **reduce** will treat its argument as a function name and call it with an argument, triggering **operator()**.

   #. You may also create your own variable of that class, e.g., **thrust::multiplies<float> foo**.   Then you just say **foo** in the argument list, not **foo()**.

   #. The optimizing compiler will replace the **operator()** function call
      with the defining expression and then continue optimizing.  So, there
      is no overhead, unlike if you passed in a pointer to a function.

#. Sometimes, e.g., in **saxpy.cu**, you see **saxpy_functor(A)**.    

   a. The class **saxpy_functor** has a constructor taking one argument.   

   #. **saxpy_functor(A)** constructs and returns a variable of class **saxpy_functor** and stores **A** in the variable.

   #. The class also overloads **operator()**.

   #. (Let's call the new variable **foo**).  **foo()** calls **operator()** for
      **foo**; its execution uses the stored **A**.

   #. Effectively, we did a **closure** of **saxpy_functor**; this is, we
      bound a property and returned a new, more restricted, variable or
      class.

#. The Thrust `examples <https://github.com/thrust/thrust/tree/master/>`_ teach several non-intuitive paradigms.  As I figure them out, I'll describe a few.  My descriptions are modified and expanded versions of the comments in the programs.  This is not a list of all the useful programs, but only of some where I am adding to their comments.

   a. **arbitrary_transformation.cu** and **dot_products_with_zip.cu**. show
      the very useful zip_iterator.  Using it is a 2-step process.

      i. Combine the separate iterators into a tuple.

      #. Construct a zip iterator from the tuple.

      Note that operator() is now a template.

   #. **boundingbox.cu** finds the bounding box around a set of 2D points.

      The main idea is to do a reduce.  However, the combining operation, instead of addition, is to combine two bounding boxes to find the box around them.

   #. **bucket_sort2d.cu** overlays a grid on a set of 2D points and finds
      the points in each grid cell (bucket).

      i. The tuple is an efficient class for a short vector of fixed length.

      #. Note how random numbers are generated.  You combine an engine that produces random output with a distribution.

         However you might need more complicated coding to make the numbers good when executing in parallel.  See **monte_carlo_disjoint_sequences.cu**.

      #. The problem is that the number of points in each cell is unpredictable.

      #. The cell containing each point is computed and that and the points are sorted to bring together the points in each cell.

      #. Then **lower_bound** and **upper_bound** are used to find each bucket in that sorted vector of points.

      #. See this `lower_bound description <http://thrust.github.io/doc/group__vectorized__binary__search.html>`_.

   #. **mode.cu** shows:

      a. Counting the number of unique keys in a vector.

	 i. Sort the vector.

	 #. Do an **inner_product**.  However, instead of the operators being **times** and **plus**, they are **not equal to the next element** and **plus**.

      #. Counting their multiplicity.

	 i. Construct vectors, sized at the number of unique keys, to hold the unique keys and counts.

	 #. Do a **reduce_by_keys** on a constant_iterator using the sorted vector as the keys.  For each range of identical keys, it sums the constant_iterator.  That is, it counts the number of identical keys.

	 #. Write a vector of unique keys and a vector of the counts.

      #. Finding the most used key (the mode).

	 i. Do **max_element** on the counts vector.

   #. **repeated_range.cu** repeats each element of an N-vector K times:
      repeated_range([0, 1, 2, 3], 2) -> [0, 0, 1, 1, 2, 2, 3, 3].  It's a lite
      version of **expand.cu**, but uses a different technique.

      a. Here, N=4 and K=2.

      #. The idea is to construct a new iterator, **repeated_range**, that, when read and
	 incremented, will return the proper output elements.

      #. The construction stores the relevant info in structure components of the variable.

      #. Treating its value like a subscript in the range [0,N*K), it divides
	 that value by K and returns that element of its input.

      See also **strided_range.cu** and **tiled_range.cu**.


Unionfs: Linux trick of the day
-------------------------------

#. aka overlay FS, translucent FS.

#. If a, b are directories, and m is an empty directory, then

   unionfs -o cow a=RW:b m

   makes m to be a combo of a and b, with a being higher priority

#. Writing a file into m writes it in a.

#. Changing a file in b writes the new version into a

#. Deleting a file in b causes a white-out note to be stored in a.

#. Unmount it thus:

   fusermount -u m

#. None of this requires superuser.  

#. Application: making a read-only directory into a read-write directory.

#. Note: IBM had a commercial version of this idea in its CP/CMS OS in the 1960s.