16_performance.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "<table>\n",
    " <tr align=left><td><img align=left src=\"https://i.creativecommons.org/l/by/4.0/88x31.png\">\n",
    " <td>Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Kyle T. Mandli</td>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Performance\n",
    "\n",
    "We have not really paid much attention to performance of code beyond basic instruction counting up until now.  These lectures are dedicated to a deeper dive into these issues and what are some important things to keep in mind."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Example: Matrix-Matrix Multiplication\n",
    "\n",
    "To start our discussion let us consider the algorithm for matrix-matrix multiplication which algorithmically looks like\n",
    "```\n",
    "do i=1:N\n",
    "    do j=1:N\n",
    "        do k=1:N\n",
    "            C[i, j] = C[i, j] + A[i, k] * B[k, j]\n",
    "        end do\n",
    "    end do\n",
    "end do\n",
    "```\n",
    "\n",
    "Consider the follow approaches to this problem:\n",
    "1. Matrix multiplication via a GCC (Fortran) intrinsic.\n",
    "1. Straight forward three-loop multiplication\n",
    "1. Parallelized double-loop using BLAS intrinsic\n",
    "1. BLAS intrinsic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "```fortran\n",
    "real function matrix_multiply_test(N,method)\n",
    "\n",
    "    use mod_rand\n",
    "    implicit none\n",
    "    \n",
    "    external DGEMM,DDOT\n",
    "    \n",
    "    double precision :: DDOT\n",
    "    integer, intent(in) :: N,method\n",
    "    integer :: start,finish,count_rate\n",
    "    double precision, dimension(:,:), allocatable :: A,B,C\n",
    "    \n",
    "    ! Local\n",
    "    integer :: i,j,k\n",
    "    \n",
    "    ! Create the random arrays\n",
    "    call init_random_seed()\n",
    "    allocate(A(N,N),B(N,N),C(N,N))\n",
    "    call random_number(A)\n",
    "    call random_number(B)\n",
    "    \n",
    "    ! Start the timer and start multiplying\n",
    "    call system_clock(start,count_rate)\n",
    "    select case(method)\n",
    "        case(1) ! Default method provided as an intrinsic method\n",
    "            C = matmul(A,B)\n",
    "        case(2) ! Simple three loop multiplication\n",
    "            !$OMP parallel do private(j,k)\n",
    "            do i=1,N\n",
    "                do j=1,N\n",
    "                    do k=1,N\n",
    "                        C(i,j) = C(i,j) + A(i,k) * B(k,j)\n",
    "                    enddo\n",
    "                enddo\n",
    "            enddo\n",
    "        case(3) ! OpenMP parallelized double loop\n",
    "            !$OMP parallel do private(j)\n",
    "            do i=1,N\n",
    "                do j=1,N\n",
    "                    C(i,j) = DDOT(N, A(i,:), 1, B(:,j), 1)\n",
    "                enddo\n",
    "            enddo\n",
    "        case(4) ! BLAS Routine call\n",
    "            ! call DGEMM(transa,transb,l,n,m,alpha,a,lda,b,ldb,beta,c,ldc)\n",
    "            call DGEMM('N', 'N', N, N, N, 1.d0, A, N, B, N, 0.d0, C, N)\n",
    "        case default\n",
    "            print *, \"***ERROR*** Invalid multiplication method!\"\n",
    "            matrix_multiply_test = -1\n",
    "            return\n",
    "    end select\n",
    "    call system_clock(finish,count_rate)\n",
    "    \n",
    "    matrix_multiply_test = float(finish - start) / float(count_rate)\n",
    "    \n",
    "end function matrix_multiply_test\n",
    "    \n",
    "program matrix_multiply\n",
    "    \n",
    "    use omp_lib\n",
    "\n",
    "    implicit none\n",
    "    \n",
    "    integer :: N, method, threads\n",
    "    character(len=10) :: input_N, input_method, input_threads\n",
    "    real :: matrix_multiply_test, time\n",
    "    \n",
    "    select case(iargc())\n",
    "        case(0)\n",
    "            N = 1000\n",
    "            method = 1\n",
    "            threads = 1\n",
    "        case(1)\n",
    "            call getarg(1,input_N)\n",
    "            read(input_N,'(I10)') N\n",
    "            method = 1\n",
    "        case(2)\n",
    "            call getarg(1,input_N)\n",
    "            call getarg(2,input_method)\n",
    "            read(input_N,'(I10)') N\n",
    "            read(input_method,'(I10)') method\n",
    "        case(3)\n",
    "            call getarg(1,input_N)\n",
    "            call getarg(2,input_method)\n",
    "            call getarg(3,input_threads)\n",
    "            read(input_N,'(I10)') N\n",
    "            read(input_method,'(I10)') method\n",
    "            read(input_threads,'(I10)') threads\n",
    "        case default\n",
    "            print *, \"***ERROR*** Too many arguments!\"\n",
    "            stop\n",
    "    end select\n",
    "    \n",
    "    !$ call omp_set_num_threads(threads)\n",
    "\n",
    "    time = matrix_multiply_test(N, method)\n",
    "    \n",
    "    print '(\"Timing for \",i5,\"x\",i5,\" matrices: \",f10.5,\" s\")',N,N,time\n",
    "    \n",
    "end program matrix_multiply\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Results\n",
    "\n",
    "Based on $1000 \\times 1000$ matrix-matrix multiply compiled with `gfortran` version 6.3.0 with the compile time flags `-O3 -funroll-loops -finline-functions -fdefault-real-8 -fopenmp`.\n",
    "\n",
    "\n",
    "Method                           | No-Threads            | Threaded\n",
    "---------------------------------|-----------------------|---------------------------\n",
    "Default mat_mult function        |            35.79600 s |                36.19100 s\n",
    "3 loop multiplication            |            39.24700 s |                10.04000 s                   \n",
    "Double loop (internal BLAS)      |             6.80500 s |                 1.76600 s                   \n",
    "BLAS Routine call                |             0.00300 s |                 0.00300 s                   \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Machine Architecture Considerations\n",
    "\n",
    "Before we can talk about performance we need to understand a bit about modern computing architectures.  Note that this glossing over a lot of important details as we will focus on only the details that we will specifically deal with."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Von Neumann Architecture\n",
    "![image](./images/vonneumann_architecture.png)\n",
    "*Kapooht - Wikipedia Commons*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Instruction Pipeline\n",
    "![image](./images/pipeline_1.png)\n",
    "*Cburnett - Wikipedia Commons*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "![image](./images/memory_architecture.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Moore's Law\n",
    "\n",
    "In 1965, Gordon Moore (co-founder of Intel) predicted that the transistor density (and hence the speed) of chips would double every 18 months for the forseeable future. This is know as Moore’s law This proved remarkably accurate for more than 40 years, see the graphs at. Note that doubling every 18 months means an increase by a factor of 4096 every 14 years."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Is Moore's Law Doomed?\n",
    "![image](./images/moores_law.png)\n",
    "*Steve Jurvetson - Wikipedia Commons - https://www.flickr.com/photos/jurvetson/31409423572/*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Current Performance Bottlenecks\n",
    "\n",
    " - Transistors can no longer be packed more densely in a single core\n",
    " - Memory is really the bottleneck\n",
    " - Hard limit due to the speed of light\n",
    " - Power consumption and therefore heat dissipation\n",
    " \n",
    "### Solutions?\n",
    " - Many-core technologies\n",
    " - Memory hierarchies\n",
    " - Algorithms that take into account these limitations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "![image](./images/memory_single_core.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "![image](./images/pipeline_2.png)\n",
    "*Cburnett - Wikipedia Commons*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Roof-line Model\n",
    "![image](./images/roofline.png)\n",
    "*Giu.natale - Wikipedia Commons*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Many-Core Architectures\n",
    "\n",
    "![image](./images/kepler_arch.png)\n",
    "![image](./images/kepler_smx.png)\n",
    "*NVIDIA - Kepler GK110/210 White Paper - http://images.nvidia.com/content/pdf/tesla/NVIDIA-Kepler-GK110-GK210-Architecture-Whitepaper.pdf*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Parallelization\n",
    "\n",
    "Parallelization is one of the primary ways we can increase performance on today's computing architectures.  There are 2+1 major types of parallelization paradigms:\n",
    " - Shared memory - each pipeline can access the memory for the entire problem\n",
    " - Distributed memory - each pipeline can only access part of the memory for the entire problem\n",
    " - Hybrid parallelism - use both..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Shared Memory\n",
    "\n",
    " - Basic construct is a *thread* - each thread has a pipeline and in the simplest case each core runs one thread\n",
    " - OpenMP, CUDA, OpenCL, OpenACC\n",
    " - Single nodes on a cluster, GPU, Xeon Phi, etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Distributed Memory\n",
    "\n",
    " - Basic contruct is a *process* \n",
    " - Each process is memory local but can communicate to other processes either on the same CPU or across a network\n",
    " - Each process can have multiple threads (hybrid parallelism)\n",
    " - MPI is most common\n",
    " - Clusters, super-computers, etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Scalability\n",
    "\n",
    "Measures of parallel performance:\n",
    "\n",
    " - Strong Scaling:  Execution time decreases inversely proportional to the number of processes\n",
    "   - Fixed size problem\n",
    " - Weak Scaling: Execution time remains constant as problem size and processes number are increased proportionally\n",
    "   - Variable size problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### OpenMP\n",
    "\n",
    "OpenMP is defined by a set of *directives* that are put into code which on compilation a compiler can turn into multi-threaded code."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Simple hello world from OpenMP.  Here we fetch the number of threads and print out the unique ID given to each.\n",
    "```fortran\n",
    "program hello_world_omp\n",
    "    \n",
    "    use omp_lib\n",
    "\n",
    "    implicit none\n",
    "    integer :: num_threads, thread_id\n",
    "\n",
    "    !$OMP parallel private(num_threads, thread_id)\n",
    "    !$ num_threads = omp_get_num_threads()\n",
    "    !$ thread_id = omp_get_thread_num()\n",
    "    print *, 'Hello from thread number', thread_id + 1, &\n",
    "             ' of ', num_threads, ' processes'\n",
    "\n",
    "    !$OMP end parallel\n",
    "\n",
    "end program hello_world_omp\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "```fortran\n",
    "program yeval\n",
    "   \n",
    "   use omp_lib\n",
    "\n",
    "   implicit none\n",
    "\n",
    "   integer(kind=8), parameter :: n = 2**16\n",
    "   integer(kind=4) :: i, nthreads\n",
    "   real(kind=8), dimension(n) :: y\n",
    "   real(kind=8) :: dx, x\n",
    "\n",
    "   ! Specify number of threads to use:\n",
    "   !$ print *, \"How many threads to use? \"\n",
    "   !$ read *, nthreads\n",
    "   !$ call omp_set_num_threads(nthreads)\n",
    "   !$ print \"('Using OpenMP with ',i3,' threads')\", nthreads\n",
    "\n",
    "   dx = 1.d0 / (n+1.d0)\n",
    "\n",
    "   !$omp parallel do private(x) \n",
    "   do i=1, n\n",
    "      x = i * dx\n",
    "      y(i) = exp(x) * cos(x) * sin(x) * sqrt(5.d0 * x + 6.d0)\n",
    "   enddo\n",
    "   !$omp end parallel do\n",
    "\n",
    "   print *, \"Filled vector y of length\", n\n",
    "\n",
    "end program yeval\n",
    "```\n",
    "*Modified from amath 583 - R.J. LeVeque - http://faculty.washington.edu/rjl/classes/am583s2014/notes/openmp.html*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Fine-Grain vs. Coarse-Grain Parallelism\n",
    "\n",
    "Consider the problem of normalizing a vector which requires two steps:\n",
    "1. Compute the norm of the vector, and\n",
    "1. Divide each entry of the vector by the norm.\n",
    "\n",
    "Unfortunately we need to loop over every entry in the vector to compute the norm **before** we can perform the division of each entry.  There are two ways to tackle this problem,\n",
    " - Let the compiler decide what thread takes what entries (fine grain) - large number of small tasks\n",
    " - Let the programmer explicitly control which entries are handled by each thread (coarse grain) - small number of large tasks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "```fortran\n",
    "program fine_grain\n",
    "   \n",
    "    use omp_lib\n",
    "    implicit none\n",
    "    integer :: i, thread_num\n",
    "    integer, parameter :: n = 1000\n",
    " \n",
    "    real(kind=8), dimension(n) :: x, y\n",
    "    real(kind=8) :: norm,ynorm\n",
    " \n",
    "    integer :: nthreads\n",
    "    \n",
    "    ! Specify number of threads to use:\n",
    "    nthreads = 1       ! need this value in serial mode\n",
    "    !$ nthreads = 4    \n",
    "    !$ call omp_set_num_threads(nthreads)\n",
    "    !$ print \"('Using OpenMP with ',i3,' threads')\", nthreads\n",
    "\n",
    "    ! Specify number of threads to use:\n",
    "    !$ call omp_set_num_threads(4)\n",
    " \n",
    "    ! initialize x:\n",
    "    !$omp parallel do \n",
    "    do i=1,n\n",
    "        x(i) = real(i, kind=8)  ! convert to double float\n",
    "    enddo\n",
    "\n",
    "    norm = 0.d0\n",
    "    ynorm = 0.d0\n",
    "\n",
    "    !$omp parallel private(i)\n",
    "\n",
    "    !$omp do reduction(+ : norm)\n",
    "    do i=1,n\n",
    "        norm = norm + abs(x(i))\n",
    "    enddo\n",
    "\n",
    "     !$omp barrier   ! not needed (implicit)\n",
    "\n",
    "    !$omp do reduction(+ : ynorm)\n",
    "    do i=1,n\n",
    "        y(i) = x(i) / norm\n",
    "        ynorm = ynorm + abs(y(i))\n",
    "    enddo\n",
    "    \n",
    "    !$omp end parallel\n",
    "\n",
    "    print *, \"norm of x = \",norm, \"  n(n+1)/2 = \",n*(n+1)/2\n",
    "    print *, 'ynorm should be 1.0:   ynorm = ', ynorm\n",
    "\n",
    "end program fine_grain\n",
    "\n",
    "```\n",
    "*Modified from amath 583 - R.J. LeVeque - http://faculty.washington.edu/rjl/classes/am583s2014/notes/openmp.html*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "```fortran\n",
    "program coarse_grain\n",
    "    \n",
    "    use omp_lib\n",
    "    implicit none\n",
    "    integer, parameter :: n = 1000\n",
    "    real(kind=8), dimension(n) :: x,y\n",
    "    real(kind=8) :: norm,norm_thread,ynorm,ynorm_thread\n",
    "    integer :: nthreads, points_per_thread,thread_num\n",
    "    integer :: i,istart,iend\n",
    "\n",
    "    ! Specify number of threads to use:\n",
    "    nthreads = 1       ! need this value in serial mode\n",
    "    !$ nthreads = 4    \n",
    "    !$ call omp_set_num_threads(nthreads)\n",
    "    !$ print \"('Using OpenMP with ',i3,' threads')\", nthreads\n",
    "\n",
    "    ! Determine how many points to handle with each thread.\n",
    "    ! Note that dividing two integers and assigning to an integer will\n",
    "    ! round down if the result is not an integer.  \n",
    "    ! This, together with the min(...) in the definition of iend below,\n",
    "    ! insures that all points will get distributed to some thread.\n",
    "    points_per_thread = (n + nthreads - 1) / nthreads\n",
    "    print *, \"points_per_thread = \",points_per_thread\n",
    "\n",
    "    ! initialize x:\n",
    "    do i=1,n\n",
    "        x(i) = dble(i)  ! convert to double float\n",
    "        enddo\n",
    "\n",
    "    norm = 0.d0\n",
    "    ynorm = 0.d0\n",
    "\n",
    "    !$omp parallel private(i,norm_thread, &\n",
    "    !$omp                  istart,iend,thread_num,ynorm_thread) \n",
    "\n",
    "    thread_num = 0     ! needed in serial mode\n",
    "    !$ thread_num = omp_get_thread_num()    ! unique for each thread\n",
    "\n",
    "    ! Determine start and end index for the set of points to be \n",
    "    ! handled by this thread:\n",
    "    istart = thread_num * points_per_thread + 1\n",
    "    iend = min((thread_num+1) * points_per_thread, n)\n",
    "\n",
    "    !$omp critical\n",
    "    print \"(\"Thread \",i2,\" will take i = \",i6,\" through i = \",i6)\", thread_num, istart, iend\n",
    "    !$omp end critical\n",
    "\n",
    "    norm_thread = 0.d0\n",
    "    do i=istart,iend\n",
    "        norm_thread = norm_thread + abs(x(i))\n",
    "        enddo\n",
    "\n",
    "    ! update global norm with value from each thread:\n",
    "    !$omp critical\n",
    "      norm = norm + norm_thread\n",
    "      print *, \"norm updated to: \",norm\n",
    "    !$omp end critical\n",
    "\n",
    "    ! make sure all have updated norm before proceeding:\n",
    "    !$omp barrier\n",
    "\n",
    "    ynorm_thread = 0.d0\n",
    "    do i=istart,iend\n",
    "        y(i) = x(i) / norm\n",
    "        ynorm_thread = ynorm_thread + abs(y(i))\n",
    "        enddo\n",
    "\n",
    "    ! update global ynorm with value from each thread:\n",
    "    !$omp critical\n",
    "      ynorm = ynorm + ynorm_thread\n",
    "      print *, \"ynorm updated to: \",ynorm\n",
    "    !$omp end critical\n",
    "    !$omp barrier\n",
    "\n",
    "    !$omp end parallel \n",
    "\n",
    "    print *, \"norm of x = \",norm, \"  n(n+1)/2 = \",n*(n+1)/2\n",
    "    print *, 'ynorm should be 1.0:   ynorm = ', ynorm\n",
    "\n",
    "end program coarse_grain\n",
    "```\n",
    "*Modified from amath 583 - R.J. LeVeque - http://faculty.washington.edu/rjl/classes/am583s2014/notes/openmp.html*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### MPI\n",
    "\n",
    "The Message Passing Interface (MPI) standard specifies a programming interface for communicating between processes possibly across a heterogeneous network.  Unlike OpenMP that uses directives MPI uses a series of function calls to communicate between parallel units.  This requires a lot more coordination from the programmer and lends itself to coarse-grain parallelism."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Additional Considerations for MPI\n",
    "\n",
    " - Topology of network\n",
    " - One-sided communication vs. two-sided\n",
    " - Halting operations vs. asynchronous communication"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example:  Hello World (Fortran)\n",
    "\n",
    "Hello world from MPI.  Each process will print out it's process number `proc_num` and the number of processes `num_procs`.  Note we need to run this using `mpiexec -n 4 hello_world_mpi` for this to work."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Basic MPI calls being used here:\n",
    " - `mpi_init(integer error)` - Needs to be called at the beginning of all MPI programs.\n",
    " - `mpi_comm_size(MPI_Comm comm, integer num_processes, integer error)` - Needs the MPI_Comm `comm` and returns the number of processes `num_processes`.\n",
    " - `mpi_comm_rank(MPI_Comm comm, integer process_id, integer error)` - Needs the MPI_Comm `comm` and returns this process's ID.\n",
    " - `mpi_finalize(integer error)` - Stop MPI, should be called at the end of an MPI program.\n",
    "\n",
    "Note that all of the Fortran definitions of MPI calls have `integer error` at the end which have values dependent on what might have gone wrong (call specific).\n",
    "\n",
    "```fortran\n",
    "program hello_world_mpi\n",
    "    \n",
    "    use mpi\n",
    "\n",
    "    implicit none\n",
    "    integer :: error, num_procs, proc_num\n",
    "\n",
    "    call mpi_init(error)\n",
    "    call mpi_comm_size(MPI_COMM_WORLD, num_procs, error)\n",
    "    call mpi_comm_rank(MPI_COMM_WORLD, proc_num, error)\n",
    "\n",
    "    print *, 'Hello from Process number', proc_num + 1, &\n",
    "             ' of ', num_procs, ' processes'\n",
    "\n",
    "    call mpi_finalize(error)\n",
    "\n",
    "end program hello_world_mpi\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example:  Hello World (Python)\n",
    "\n",
    "We can use MPI in Python via the package `mpi4py` which contains nearly a 1-to-1 correspondance to the C interface.\n",
    "\n",
    "Note that this cannot as is be run from a notebook.  We need to run this from the command line with `mpiexec -n 4 python hello_world.py`.\n",
    "```python\n",
    "#!/usr/bin/env python\n",
    "# encoding: utf-8\n",
    "\n",
    "\"\"\"MPI Hello World from Python\"\"\"\n",
    "\n",
    "from __future__ import absolute_import\n",
    "from __future__ import print_function\n",
    "\n",
    "from mpi4py import MPI\n",
    "\n",
    "comm = MPI.COMM_WORLD\n",
    "rank = comm.Get_rank()\n",
    "size = comm.Get_size()\n",
    "\n",
    "print('Hello from Process number %s of %s processes.' % (rank + 1, size))\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example: Computing $\\pi$ (Fortran)\n",
    "\n",
    "Compute $\\pi$ in parallel by approximating the integral \n",
    "$$\n",
    "    \\int^1_0 \\frac{1}{1 + x^2} dx = \\frac{\\pi}{4}\n",
    "$$\n",
    "with the midpoint rule.  In this case we will divide up the interval into a user specified number of points and let each processor take the same number of points (give or take a point).  This process has the following steps:\n",
    "\n",
    "1. Process 0 asks the user how many points they want to use.\n",
    "1. This number $N$ is then broadcast to every other process.\n",
    "1. Each process then computes the points it is responsible for and computes its portion of the sum.\n",
    "1. The sums from each process are then reduced (added) together and sent to process 0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The following MPI functions are used below:\n",
    "\n",
    " - `mpi_bcast( void *buffer, integer count, integer datatype, int root, MPI_Comm comm, integer error)` where `count` is the number of entries in the `buffer`, `datatype` is the type of data in the buffer, `root` is the process that is broadcasting the buffer, and `comm` is the MPI_communicator (`MPI_COMM_WORLD`).\n",
    " - `mpi_Reduce(const void *send_buffer, void *receive_buffer, integer count, integer datatype, MPI_Op operation, integer root, MPI_Comm comm, integer error)` where `send_buffer` is the buffer being of data being sent, `receive_buffer` is the buffer used to receive the data (on the process chosen to do so), `count` is the size of the `buffer`, `datatype` the type of data, `operation` the type of operation to perform, `root` the process that will receive the buffer, and `comm` the MPI_comm world.\n",
    "\n",
    "```fortran\n",
    "! Example program that computes pi in parallel\n",
    "program compute_pi\n",
    "\n",
    "    use mpi\n",
    "    \n",
    "    implicit none\n",
    "\n",
    "    integer :: error, num_procs, proc_id, points_per_proc, n, i, start, end\n",
    "    real(kind=8) :: x, dx, pi_sum, pi_sum_proc\n",
    "\n",
    "    real(kind=8), parameter :: pi = 3.1415926535897932384626433832795d0\n",
    "\n",
    "    call mpi_init(error)\n",
    "    call mpi_comm_size(MPI_COMM_WORLD, num_procs, error)\n",
    "    call mpi_comm_rank(MPI_COMM_WORLD, proc_id, error)\n",
    "\n",
    "    ! Proc 0 will ask the user for the number of points\n",
    "    if (proc_id == 0) then\n",
    "        print *, \"Using \",num_procs,\" processors\"\n",
    "        print *, \"Input n ... \"\n",
    "        read *, n\n",
    "    end if\n",
    "\n",
    "    ! Broadcast to all procs; everybody gets the value of n from proc 0\n",
    "    call mpi_bcast(n, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, error)\n",
    "\n",
    "    dx = 1.d0 / n\n",
    "\n",
    "    ! Determine how many points to handle with each proc\n",
    "    points_per_proc = (n + num_procs - 1) / num_procs\n",
    "    ! Only print out the number of points per proc by process 0\n",
    "    if (proc_id == 0) then   \n",
    "        print *, \"points_per_proc = \", points_per_proc\n",
    "    end if\n",
    "\n",
    "    ! Determine start and end index for this proc's points\n",
    "    start = proc_id * points_per_proc + 1\n",
    "    end = min((proc_id + 1) * points_per_proc, n)\n",
    "\n",
    "    ! Diagnostic: tell the user which points will be handled by which proc\n",
    "    print '(\"Process \",i2,\" will take i = \",i8,\" through i = \",i8)', &\n",
    "          proc_id, start, end\n",
    "\n",
    "    pi_sum_proc = 0.d0\n",
    "    do i=start,end\n",
    "        x = (i - 0.5d0) * dx\n",
    "        pi_sum_proc = pi_sum_proc + 1.d0 / (1.d0 + x**2)\n",
    "    enddo\n",
    "\n",
    "    call MPI_REDUCE(pi_sum_proc, pi_sum, 1, MPI_DOUBLE_PRECISION, MPI_SUM, 0, &\n",
    "                        MPI_COMM_WORLD, error)\n",
    "\n",
    "    if (proc_id == 0) then\n",
    "        print *, \"The approximation to pi is \", 4.d0 * dx * pi_sum\n",
    "        print *, \"Difference = \", abs(pi - 4.d0 * dx * pi_sum)\n",
    "    endif\n",
    "\n",
    "    call mpi_finalize(error)\n",
    "\n",
    "end program compute_pi\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example: Computing $\\pi$ (Python)\n",
    "```python\n",
    "#!/usr/bin/env python\n",
    "# encoding: utf-8\n",
    "\n",
    "\"\"\"Note passing from Python\"\"\"\n",
    "\n",
    "from __future__ import absolute_import\n",
    "from __future__ import print_function\n",
    "\n",
    "from mpi4py import MPI\n",
    "import numpy\n",
    "\n",
    "comm = MPI.COMM_WORLD\n",
    "rank = comm.Get_rank()\n",
    "size = comm.Get_size()\n",
    "\n",
    "# This does not seem to work well in Python\n",
    "# if rank == 0:\n",
    "#     print(\"Using %s processes\" % size)\n",
    "#     N = input(\"  Input N...\")\n",
    "\n",
    "N = 100\n",
    "N = comm.bcast(N, root=0)\n",
    "\n",
    "dx = 1.0 / N\n",
    "\n",
    "points_per_proc = (N + size - 1) / size\n",
    "\n",
    "if rank == 0:\n",
    "    print(\"Points/process = %s\" % points_per_proc)\n",
    "\n",
    "start = rank * points_per_proc + 1\n",
    "end = min((rank + 1) * points_per_proc, N)\n",
    "\n",
    "print(\"Process %s will take i = %s through i = %s\" % (rank, start, end))\n",
    "\n",
    "pi_sum_local = 0.0\n",
    "for i in range(start, end + 1):\n",
    "    x = (i - 0.5) * dx\n",
    "    pi_sum_local += 1.0 / (1.0 + x**2)\n",
    "pi_sum_local = numpy.array(pi_sum_local, dtype='d')\n",
    "\n",
    "pi_sum = numpy.zeros(1)\n",
    "comm.Reduce([pi_sum_local, MPI.DOUBLE], [pi_sum, MPI.DOUBLE],\n",
    "            op=MPI.SUM, root=0)\n",
    "\n",
    "if rank == 0:\n",
    "    pi = 4.0 * dx * pi_sum\n",
    "    print(\"The approximation to pi is %s\" % (pi))\n",
    "    print(\"Difference = %s\" % (numpy.abs(numpy.pi - pi)))\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example:  Passing Notes (Fortran)\n",
    "\n",
    "This example shows how one process can communicate to one other process until the message has been passed to the last process.  In this case we are only passing a single double precision number (not so interesting).  In this case each process that has not yet received the note will wait until it does.  This is called *blocking* and is not a desired behavior."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The only new function in this case is\n",
    " - `MPI_Recv(void *buffer, integer count, integer datatype, integer source, integer tag, MPI_Comm comm, MPI_Status status, integer error)` - Here we again have the data beening received in `buffer`, the number of values `count`, the type of data `datatype`, the source process is `source`, an identifying `tag`, and the MPI communicator `comm`.  The values `status` says the status of the transfer.  Note that here once we call `mpi_recv` the process will *block* until it receives the data (not great).\n",
    " \n",
    "```fortran\n",
    "! Demonstration for message passing between MPI processes\n",
    "program note_passing\n",
    "\n",
    "    use mpi\n",
    "\n",
    "    implicit none\n",
    "\n",
    "    integer :: proc_id, num_procs, error, tag\n",
    "    real(kind=8) :: important_note\n",
    "    integer, dimension(MPI_STATUS_SIZE) :: status\n",
    "\n",
    "    call mpi_init(error)\n",
    "    call mpi_comm_size(MPI_COMM_WORLD, num_procs, error)\n",
    "    call mpi_comm_rank(MPI_COMM_WORLD, proc_id, error)\n",
    "\n",
    "    if (num_procs == 1) then\n",
    "        print *, \"Only one process, cannot do anything.\"\n",
    "        call MPI_FINALIZE(error)\n",
    "        stop\n",
    "    endif\n",
    "\n",
    "    ! Not really important in this case but important to sort through messages\n",
    "    tag = 42\n",
    "\n",
    "    ! If we are process 0 then set the value to be passed\n",
    "    if (proc_id == 0) then\n",
    "        important_note = 2.718281828d0\n",
    "        print '(\"Process \",i3,\" note = \",e18.8)', proc_id, important_note\n",
    "\n",
    "        call mpi_send(important_note, 1, MPI_DOUBLE_PRECISION, 1, tag, &\n",
    "                      MPI_COMM_WORLD, error)\n",
    "\n",
    "    ! If we are one of the processes in between pass it on to the next process\n",
    "    else if (proc_id < num_procs - 1) then\n",
    "\n",
    "        call MPI_RECV(important_note, 1, MPI_DOUBLE_PRECISION, proc_id-1, tag, &\n",
    "                      MPI_COMM_WORLD, status, error)\n",
    "\n",
    "        print '(\"Process \",i3,\" received note = \",e18.8)', proc_id, important_note\n",
    "\n",
    "        call mpi_send(important_note, 1, MPI_DOUBLE_PRECISION, proc_id + 1, &\n",
    "                      tag, MPI_COMM_WORLD, error)\n",
    "\n",
    "        print '(\"Process \",i3,\"    sent note = \",e18.8)', proc_id, important_note\n",
    "\n",
    "    ! If we are the last process in the class to find out, well...\n",
    "    else if (proc_id == num_procs - 1) then\n",
    "\n",
    "        call mpi_recv(important_note, 1, MPI_DOUBLE_PRECISION, proc_id - 1, &\n",
    "                      tag, MPI_COMM_WORLD, status, error)\n",
    "\n",
    "        print '(\"Process \",i3,\" ends up with note = \",e18.8)', proc_id, important_note\n",
    "      endif\n",
    "\n",
    "    call MPI_FINALIZE(error)\n",
    "\n",
    "end program note_passing\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example:  Passing Notes (Python)\n",
    "\n",
    "```python\n",
    "#!/usr/bin/env python\n",
    "# encoding: utf-8\n",
    "\n",
    "\"\"\"Note passing from Python\"\"\"\n",
    "\n",
    "from __future__ import absolute_import\n",
    "from __future__ import print_function\n",
    "\n",
    "from mpi4py import MPI\n",
    "import numpy\n",
    "\n",
    "comm = MPI.COMM_WORLD\n",
    "rank = comm.Get_rank()\n",
    "size = comm.Get_size()\n",
    "\n",
    "tag = 42\n",
    "N = 10\n",
    "\n",
    "if rank == 0:\n",
    "    data = numpy.arange(N, dtype=numpy.float64)\n",
    "    print(\"Process %s note = %s\" % (rank, data))\n",
    "    # Note here that MPI datatype discovery is automatic\n",
    "    comm.Send(data, dest=rank + 1, tag=tag)\n",
    "\n",
    "elif rank < size - 1:\n",
    "    data = numpy.empty(N, dtype=numpy.float64)\n",
    "    comm.Recv(data, source=rank - 1, tag=tag)\n",
    "\n",
    "    print(\"Process %s note = %s\" % (rank, data))\n",
    "\n",
    "    comm.Send(data, dest=rank + 1, tag=tag)\n",
    "\n",
    "elif rank == size - 1:\n",
    "    data = numpy.empty(N, dtype=numpy.float64)\n",
    "    comm.Recv(data, source=rank - 1, tag=tag)\n",
    "\n",
    "    print(\"Process %s note = %s\" % (rank, data))\n",
    "\n",
    "else:\n",
    "    raise Exception(\"Invalid rank.\")\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example - Poisson Equation Solution Using Jacobi\n",
    "\n",
    "Finally lets look at the parallelization of the solution of the Poisson equation using Jacobi iteration.  The way we parallelize Jacobi is straight forward, we divide up the row evaluations between the available threads/processes.  There are three different ways to do this that we will go over:\n",
    "1. Fine grained where we put a OpenMP directive around individual loops and let the threads take work when available.\n",
    "1. Coarse grained parallelism and OpenMP where we explicitly decide which threads get what work.\n",
    "1. Coarse grained parallelism almost identical to the above but with MPI."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### OpenMP - Fine-Grained Parallelism\n",
    "```fortran\n",
    "! Solve Poisson equation\n",
    "!  u_{xx} = f(x)    x \\in [a, b]\n",
    "! with \n",
    "!  u(a) = alpha, u(b) = beta\n",
    "! using Jacobi iteration and a fine-grain parallelism approach using OpenMP.\n",
    "program jacobi_omp1\n",
    "\n",
    "    use omp_lib\n",
    "    \n",
    "    implicit none\n",
    "    \n",
    "    ! Problem specification and storage\n",
    "    real(kind=8) :: a, b, alpha, beta\n",
    "    real(kind=8), dimension(:), allocatable :: x, u, u_new, f\n",
    "    real(kind=8) :: dx, tolerance, du_max\n",
    "\n",
    "    integer(kind=4) :: n, num_threads\n",
    "    integer(kind=8) :: i, num_iterations\n",
    "    integer(kind=8), parameter :: MAX_ITERATIONS = 100000\n",
    "    real(kind=8) :: time(2)\n",
    "\n",
    "    ! Boundaries\n",
    "    a = 0.d0\n",
    "    b = 1.d0\n",
    "    alpha = 0.d0\n",
    "    beta = 3.d0\n",
    "\n",
    "    ! Specify number of threads to use:\n",
    "    num_threads = 4\n",
    "    !$ call omp_set_num_threads(num_threads)\n",
    "    !$ print \"('Using OpenMP with ',i3,' threads')\", num_threads\n",
    "\n",
    "    N = 100\n",
    "\n",
    "    ! Allocate storage for boundary points too\n",
    "    allocate(x(0:N + 1), u(0:N + 1), u_new(0:N + 1), f(0:N + 1))\n",
    "\n",
    "    call cpu_time(time(1))\n",
    "\n",
    "    ! grid spacing:\n",
    "    dx = (b - a) / (N + 1.d0)\n",
    "\n",
    "    ! Set iniital guess and construct the grid\n",
    "    ! Note that here we are breaking up the problem already to the threads\n",
    "    !$omp parallel do\n",
    "    do i=0, N + 1\n",
    "        ! grid points:\n",
    "        x(i) = i * dx\n",
    "        ! source term:\n",
    "        f(i) = exp(x(i))\n",
    "        ! initial guess (satisfies boundary conditions and sets them)\n",
    "        u(i) = alpha + x(i) * (beta - alpha)\n",
    "    enddo\n",
    "\n",
    "    ! Tolerance\n",
    "    tolerance = 0.1d0 * dx**2\n",
    "\n",
    "    ! Main Jacobi iteration loop\n",
    "    ! Copy old values to new\n",
    "    num_iterations = 0\n",
    "    du_max = 1d99\n",
    "    do while (du_max >= tolerance .and. num_iterations <= MAX_ITERATIONS)\n",
    "        du_max = 0.d0\n",
    "        !$omp parallel do reduction(max : du_max)\n",
    "        do i=1, N\n",
    "            u_new(i) = 0.5d0 * (u(i-1) + u(i+1) - dx**2 * f(i))\n",
    "            ! print *, abs(u(i) - u_old(i))\n",
    "            du_max = max(du_max, abs(u_new(i) - u(i)))\n",
    "        end do\n",
    "\n",
    "        if (mod(num_iterations, 1000) == 0) then\n",
    "            print *, \"du_max, iteration = \", du_max, num_iterations\n",
    "        end if\n",
    "\n",
    "        ! Copy old data into new\n",
    "        !$omp parallel do \n",
    "        do i=1, N\n",
    "            u(i) = u_new(i)\n",
    "        end do\n",
    "        num_iterations = num_iterations + 1\n",
    "    end do\n",
    "\n",
    "    call cpu_time(time(2))\n",
    "    if (num_iterations > MAX_ITERATIONS) then\n",
    "        print *, \"Iteration failed!\"\n",
    "        stop\n",
    "    end if\n",
    "    print '(\"CPU time = \",f12.8, \" seconds\")', time(2) - time(1)\n",
    "    print *, \"Total number of iterations: \", num_iterations\n",
    "\n",
    "    ! Write out solution for checking\n",
    "    open(unit=20, file=\"jacobi_omp1.txt\", status=\"unknown\")\n",
    "    do i=0,n+1\n",
    "        write(20,'(2e20.10)') x(i), u(i)\n",
    "    enddo\n",
    "    close(20)\n",
    "\n",
    "end program jacobi_omp1\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### OpenMP - Coarse Grained Parallelism\n",
    "\n",
    "Here we coordinate which threads get what but also have to be very careful about what code gets executed so as to not introduce race-conditions.\n",
    "\n",
    "```fortran\n",
    "program jacobi_omp2\n",
    "\n",
    "    use omp_lib\n",
    "    \n",
    "    implicit none\n",
    "    \n",
    "    ! Problem specification and storage\n",
    "    real(kind=8) :: a, b, alpha, beta\n",
    "    real(kind=8), dimension(:), allocatable :: x, u, u_new, f\n",
    "    real(kind=8) :: dx, tolerance, du_max, du_max_thread\n",
    "\n",
    "    integer(kind=4) :: n, num_threads, points_per_thread, thread_id, start, end\n",
    "    integer(kind=8) :: i, num_iterations\n",
    "    integer(kind=8), parameter :: MAX_ITERATIONS = 100000\n",
    "    real(kind=8) :: time(2)\n",
    "\n",
    "    ! Boundaries\n",
    "    a = 0.d0\n",
    "    b = 1.d0\n",
    "    alpha = 0.d0\n",
    "    beta = 3.d0\n",
    "\n",
    "    ! Specify number of threads to use:\n",
    "    num_threads = 2\n",
    "    !$ call omp_set_num_threads(num_threads)\n",
    "    !$ print \"('Using OpenMP with ',i3,' threads')\", num_threads\n",
    "\n",
    "    N = 100\n",
    "\n",
    "    ! Allocate storage for boundary points too\n",
    "    allocate(x(0:N + 1), u(0:N + 1), u_new(0:N + 1), f(0:N + 1))\n",
    "\n",
    "    call cpu_time(time(1))\n",
    "\n",
    "    ! grid spacing:\n",
    "    dx = (b - a) / (N + 1.d0)\n",
    "\n",
    "    ! Tolerance\n",
    "    tolerance = 0.1d0 * dx**2\n",
    "\n",
    "    ! Determine how many points to handle with each thread.\n",
    "    ! Note that dividing two integers and assigning to an integer will\n",
    "    ! round down if the result is not an integer.  \n",
    "    ! This, together with the min(...) in the definition of iend below,\n",
    "    ! insures that all points will get distributed to some thread.\n",
    "    points_per_thread = (n + num_threads - 1) / num_threads\n",
    "    print *, \"points_per_thread = \", points_per_thread\n",
    "\n",
    "    ! Start of the parallel block... \n",
    "    ! ------------------------------\n",
    "\n",
    "    ! This is the only time threads are forked in this program:\n",
    "    !$omp parallel private(thread_id, num_iterations, start, end, i,   &\n",
    "    !$OMP                  du_max_thread)\n",
    "\n",
    "    ! Set thread is, default to 0 if in serial\n",
    "    thread_id = 0\n",
    "    !$ thread_id = omp_get_thread_num()\n",
    "\n",
    "    ! Determine start and end index\n",
    "    start = thread_id * points_per_thread + 1\n",
    "    end = min((thread_id + 1) * points_per_thread, N)\n",
    "\n",
    "    ! Output some thread information and indexing\n",
    "    !$omp critical\n",
    "    print '(\"Thread \",i2,\" will take i = \",i6,\" through i = \",i6)', &\n",
    "          thread_id, start, end\n",
    "    !$omp end critical\n",
    "\n",
    "    ! Set iniital guess and construct the grid\n",
    "    do i=start, end\n",
    "        ! grid points:\n",
    "        x(i) = i * dx\n",
    "        ! source term:\n",
    "        f(i) = exp(x(i))\n",
    "        ! initial guess (satisfies boundary conditions and sets them)\n",
    "        u(i) = alpha + x(i) * (beta - alpha)\n",
    "    enddo\n",
    "    ! Note that the above does not set the boundaries, do this in a single thread\n",
    "    !$omp single\n",
    "    x(0) = a\n",
    "    x(N + 1) = b\n",
    "    u(0) = alpha\n",
    "    u(N + 1) = beta\n",
    "    !$omp end single nowait\n",
    "\n",
    "    ! Main Jacobi iteration loop\n",
    "    do num_iterations=1, MAX_ITERATIONS\n",
    "        \n",
    "        ! Make one thread reset the global du_max\n",
    "        !$omp single\n",
    "        du_max = 0.d0\n",
    "        !$omp end single\n",
    "\n",
    "        ! Private to each thread\n",
    "        du_max_thread = 0.d0\n",
    "        do i=start, end\n",
    "            u_new(i) = 0.5d0 * (u(i-1) + u(i+1) - dx**2 * f(i))\n",
    "            du_max_thread = max(du_max_thread, abs(u_new(i) - u(i)))\n",
    "        end do\n",
    "\n",
    "        ! Compute global du_max\n",
    "        !$omp critical\n",
    "        du_max = max(du_max, du_max_thread)\n",
    "        !$omp end critical\n",
    "\n",
    "        ! Make sure all threads are done contributing to du_max\n",
    "        !$omp barrier\n",
    "\n",
    "        ! Have one thread print out the convergence info\n",
    "        !$omp single\n",
    "        if (mod(num_iterations, 1000) == 0) then\n",
    "            print '(\"After \",i8,\" iterations, dumax = \",d16.6,/)', num_iterations, du_max\n",
    "        end if\n",
    "        !$omp end single nowait\n",
    "\n",
    "        ! Copy new data into u\n",
    "        do i=start, end\n",
    "            u(i) = u_new(i)\n",
    "        end do\n",
    "\n",
    "        ! Check exit criteria\n",
    "        if (du_max < tolerance) then\n",
    "            exit\n",
    "        end if\n",
    "\n",
    "        ! Make sure all threads are caught up to this point before starting\n",
    "        ! another iteration\n",
    "        !$omp barrier\n",
    "    end do\n",
    "\n",
    "    !$omp end parallel\n",
    "\n",
    "    call cpu_time(time(2))\n",
    "    if (num_iterations > MAX_ITERATIONS) then\n",
    "        print *, \"Iteration failed!\"\n",
    "        stop\n",
    "    end if\n",
    "    print '(\"CPU time = \",f12.8, \" seconds\")', time(2) - time(1)\n",
    "    print *, \"Total number of iterations: \", num_iterations\n",
    "\n",
    "    ! Write out solution for checking\n",
    "    open(unit=20, file=\"jacobi_omp2.txt\", status=\"unknown\")\n",
    "    do i=0,n+1\n",
    "        write(20,'(2e20.10)') x(i), u(i)\n",
    "    enddo\n",
    "    close(20)\n",
    "\n",
    "end program jacobi_omp2\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### MPI\n",
    "```fortran\n",
    "program jacobi1d_mpi\n",
    "    use mpi\n",
    "\n",
    "    implicit none\n",
    "\n",
    "    integer, parameter :: maxiter = 100000, nprint = 5000\n",
    "    real (kind=8), parameter :: alpha = 20.d0, beta = 60.d0\n",
    "\n",
    "    integer :: i, iter, istart, iend, points_per_task, itask, n\n",
    "    integer :: ierr, ntasks, me, req1, req2\n",
    "    integer, dimension(MPI_STATUS_SIZE) :: mpistatus\n",
    "    real (kind = 8), dimension(:), allocatable :: f, u, uold\n",
    "    real (kind = 8) :: x, dumax_task, dumax_global, dx, tol\n",
    "\n",
    "    ! Initialize MPI; get total number of tasks and ID of this task\n",
    "    call mpi_init(ierr)\n",
    "    call mpi_comm_size(MPI_COMM_WORLD, ntasks, ierr)\n",
    "    call mpi_comm_rank(MPI_COMM_WORLD, me, ierr)\n",
    "\n",
    "    ! Ask the user for the number of points\n",
    "    if (me == 0) then\n",
    "        print *, \"Input n ... \"\n",
    "        read *, n\n",
    "    end if\n",
    "    ! Broadcast to all tasks; everybody gets the value of n from task 0\n",
    "    call mpi_bcast(n, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)\n",
    "\n",
    "    dx = 1.d0/real(n+1, kind=8)\n",
    "    tol = 0.1d0*dx**2\n",
    "\n",
    "    ! Determine how many points to handle with each task\n",
    "    points_per_task = (n + ntasks - 1)/ntasks\n",
    "    if (me == 0) then   ! Only one task should print to avoid clutter\n",
    "        print *, \"points_per_task = \", points_per_task\n",
    "    end if\n",
    "\n",
    "    ! Determine start and end index for this task's points\n",
    "    istart = me * points_per_task + 1\n",
    "    iend = min((me + 1)*points_per_task, n)\n",
    "\n",
    "    ! Diagnostic: tell the user which points will be handled by which task\n",
    "    print '(\"Task \",i2,\" will take i = \",i6,\" through i = \",i6)', &\n",
    "        me, istart, iend\n",
    "\n",
    "\n",
    "    ! Initialize:\n",
    "    ! -----------\n",
    "\n",
    "    ! This makes the indices run from istart-1 to iend+1\n",
    "    ! This is more or less cosmetic, but makes things easier to think about\n",
    "    allocate(f(istart-1:iend+1), u(istart-1:iend+1), uold(istart-1:iend+1))\n",
    "\n",
    "    ! Each task sets its own, independent array\n",
    "    do i = istart, iend\n",
    "        ! Each task is a single thread with all its variables private\n",
    "        ! so re-using the scalar variable x from one loop iteration to\n",
    "        ! the next does not produce a race condition.\n",
    "        x = dx*real(i, kind=8)\n",
    "        f(i) = 100.d0*exp(x)               ! Source term\n",
    "        u(i) = alpha + x*(beta - alpha)    ! Initial guess\n",
    "    end do\n",
    "    \n",
    "    ! Set boundary conditions if this task is keeping track of a boundary\n",
    "    ! point\n",
    "    if (me == 0)        u(istart-1) = alpha\n",
    "    if (me == ntasks-1) u(iend+1)   = beta\n",
    "\n",
    "\n",
    "    ! Jacobi iteratation:\n",
    "    ! -------------------\n",
    "\n",
    "    do iter = 1, maxiter\n",
    "        uold = u\n",
    "\n",
    "        ! Send endpoint values to tasks handling neighboring sections\n",
    "        ! of the array.  Note that non-blocking sends are used; note\n",
    "        ! also that this sends from uold, so the buffer we're sending\n",
    "        ! from won't be modified while it's being sent.\n",
    "        !\n",
    "        ! tag=1 is used for messages sent to the left\n",
    "        ! tag=2 is used for messages sent to the right\n",
    "\n",
    "        if (me > 0) then\n",
    "            ! Send left endpoint value to process to the \"left\"\n",
    "            call mpi_isend(uold(istart), 1, MPI_DOUBLE_PRECISION, me - 1, &\n",
    "                1, MPI_COMM_WORLD, req1, ierr)\n",
    "        end if\n",
    "        if (me < ntasks-1) then\n",
    "            ! Send right endpoint value to process on the \"right\"\n",
    "            call mpi_isend(uold(iend), 1, MPI_DOUBLE_PRECISION, me + 1, &\n",
    "                2, MPI_COMM_WORLD, req2, ierr)\n",
    "        end if\n",
    "\n",
    "        ! Accept incoming endpoint values from other tasks.  Note that\n",
    "        ! these are blocking receives, because we can't run the next step\n",
    "        ! of the Jacobi iteration until we've received all the\n",
    "        ! incoming data.\n",
    "\n",
    "        if (me < ntasks-1) then\n",
    "            ! Receive right endpoint value\n",
    "            call mpi_recv(uold(iend+1), 1, MPI_DOUBLE_PRECISION, me + 1, &\n",
    "                1, MPI_COMM_WORLD, mpistatus, ierr)\n",
    "        end if\n",
    "        if (me > 0) then\n",
    "            ! Receive left endpoint value\n",
    "            call mpi_recv(uold(istart-1), 1, MPI_DOUBLE_PRECISION, me - 1, &\n",
    "                2, MPI_COMM_WORLD, mpistatus, ierr)\n",
    "        end if\n",
    "\n",
    "        dumax_task = 0.d0   ! Max seen by this task\n",
    "\n",
    "        ! Apply Jacobi iteration on this task's section of the array\n",
    "        do i = istart, iend\n",
    "            u(i) = 0.5d0*(uold(i-1) + uold(i+1) + dx**2*f(i))\n",
    "            dumax_task = max(dumax_task, abs(u(i) - uold(i)))\n",
    "        end do\n",
    "\n",
    "        ! Take global maximum of dumax values\n",
    "        call mpi_allreduce(dumax_task, dumax_global, 1, MPI_DOUBLE_PRECISION, &\n",
    "            MPI_MAX, MPI_COMM_WORLD, ierr)\n",
    "        ! Note that this MPI_ALLREDUCE call acts as an implicit barrier,\n",
    "        ! since no process can return from it until all processes\n",
    "        ! have called it.  Because of this, after this call we know\n",
    "        ! that all the send and receive operations initiated at the\n",
    "        ! top of the loop have finished -- all the MPI_RECV calls have\n",
    "        ! finished in order for each process to get here, and if the\n",
    "        ! MPI_RECV calls have finished, the corresponding MPI_ISEND\n",
    "        ! calls have also finished.  Thus we can safely modify uold\n",
    "        ! again.\n",
    "\n",
    "        ! Also periodically report progress to the user\n",
    "        if (me == 0) then\n",
    "            if (mod(iter, nprint)==0) then\n",
    "                print '(\"After \",i8,\" iterations, dumax = \",d16.6,/)', &\n",
    "                    iter, dumax_global\n",
    "            end if\n",
    "        end if\n",
    "\n",
    "        ! All tasks now have dumax_global, and can check for convergence\n",
    "        if (dumax_global < tol) exit\n",
    "    end do\n",
    "\n",
    "    print '(\"Task number \",i2,\" finished after \",i9,\" iterations, dumax = \",&\n",
    "            e16.6)', me, iter, dumax_global\n",
    "\n",
    "\n",
    "    ! Output result:\n",
    "    ! --------------\n",
    "\n",
    "    ! Note: this only works if all processes share a file system\n",
    "    ! and can all open and write to the same file!\n",
    "\n",
    "    ! Synchronize to keep the next part from being non-deterministic\n",
    "    call mpi_barrier(MPI_COMM_WORLD, ierr)\n",
    "\n",
    "    ! Have each task output to a file in sequence, using messages to\n",
    "    ! coordinate\n",
    "\n",
    "    if (me == 0) then    ! Task 0 goes first\n",
    "        ! Open file for writing, replacing any previous version:\n",
    "        open(unit=20, file=\"heatsoln.txt\", status=\"replace\")\n",
    "        write(20,*) \"          x                  u\"\n",
    "        write(20, '(2e20.10)') 0.d0, u(0)    ! Boundary value at left end\n",
    "\n",
    "        do i = istart, iend\n",
    "            write(20, '(2e20.10)') i*dx, u(i)\n",
    "        end do\n",
    "\n",
    "        close(unit=20)\n",
    "        ! Closing the file should guarantee that all the ouput \n",
    "        ! will be written to disk.\n",
    "        ! If the file isn't closed before the next process starts writing,\n",
    "        ! output may be jumbled or missing.\n",
    "\n",
    "        ! Send go-ahead message to next task\n",
    "        ! Only the fact that the message was sent is important, not its contents\n",
    "        ! so we send the special address MPI_BOTTOM and length 0.\n",
    "        ! tag=4 is used for this message.\n",
    "\n",
    "        if (ntasks > 1) then\n",
    "            call mpi_send(MPI_BOTTOM, 0, MPI_INTEGER, 1, 4, &\n",
    "                          MPI_COMM_WORLD, ierr)\n",
    "            endif\n",
    "\n",
    "    else\n",
    "        ! Wait for go-ahead message from previous task\n",
    "        call mpi_recv(MPI_BOTTOM, 0, MPI_INTEGER, me - 1, 4, &\n",
    "                          MPI_COMM_WORLD, mpistatus, ierr)\n",
    "        ! Open file for appending; do not destroy previous contents\n",
    "        open(unit=20, file=\"heatsoln.txt\", status=\"old\", access=\"append\")\n",
    "        do i = istart, iend\n",
    "            write(20, '(2e20.10)') i*dx, u(i)\n",
    "        end do\n",
    "\n",
    "        ! Boundary value at right end:\n",
    "        if (me == ntasks - 1) write(20, '(2e20.10)') 1.d0, u(iend+1)    \n",
    "\n",
    "        ! Flush all pending writes to disk\n",
    "        close(unit=20)\n",
    "\n",
    "        if (me < ntasks - 1) then\n",
    "            ! Send go-ahead message to next task\n",
    "            call mpi_send(MPI_BOTTOM, 0, MPI_INTEGER, me + 1, 4, &\n",
    "                          MPI_COMM_WORLD, ierr)\n",
    "        end if\n",
    "    end if\n",
    "\n",
    "    ! Notify the user when all tasks have finished writing\n",
    "    if (me == ntasks - 1) print *, \"Solution is in heatsoln.txt\"\n",
    "\n",
    "    ! Close out MPI\n",
    "    call mpi_finalize(ierr)\n",
    "\n",
    "end program jacobi1d_mpi\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABCMAAAEWCAYAAABcwF1xAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xu8p2O5+PHPZVKaTYQphzJ2u9r7lw6qSdJJZxXZlUqi\ndJJKsqNUOiqdY3cQiRhmJaWU8zHnQqOQQyoxYSuDCIuYcf/+uJ5neebpu+b4fdb6rrU+79drXmbW\n957vetbCNddc93Vfd5RSkCRJkiRJGisrjfcDSJIkSZKkqcVihCRJkiRJGlMWIyRJkiRJ0piyGCFJ\nkiRJksaUxQhJkiRJkjSmLEZIkiRJkqQxZTFCU1pEHBgRnxzv51haEbFBRNwVEdM6eO/PRMScfr+v\nJE01EfGWiDi1o/c+KyLe1cV7S5I0lixGTEERsWNE/C4ihiPirxFxQESsMUaf+7qIuC8i1m59/LcR\nUSJiw+rXh1Xr7oqI2yLitIj4r+X8nDtGxMLqveof3wYopexcSvncin5di/nc20bEhRFxd0TcXP38\nfRERy/N+pZS/lFJWLaUs7PezSuqfqRZnq/dbNyIOiYibIuLOiPh9RHw2Iv5txb6ibkXEQyPiUxFx\ndRWrb4yIkyLi5cv7nqWUoVLKcv9+SeqnpflzYXF/Jrhhpa5YjJhiImJ34MvAh4HVgU2BmcBpEfHQ\nMXqMa4E3N57pKcD0Huu+UkpZFXgMcDNw2Ap8zl9Vf4mvf+yyAu+1VKrv9TeArwLrAI8GdgaeC/T8\nXnfR8SBpbE3FOBsRawK/Ah4OPKeUshrwMmAN4D+W5z0X87ke0s/3A44GtgbeCjwS+Hcydr96jD6/\nJI2FpflzoZ+5t7REFiOmkIh4BPBZ4AOllJNLKfeXUq4D3ghsCGxfrftMRBwdEUdVu1u/iYinNd5n\nvYj4SUTMj4hrI2LXxmufiYgfRcTh1e+9IiJmtR7lCDLpq70NOHy05y6lDAM/AJ68Qt+AHqoq8Oer\nn28eETdExO5VF8NNEfH2xtqHRcTXIuIvEfG3yCMeDx/lfVcH9gbeV0o5upRyZ0m/LaW8pZTyz8bn\nPyAiToyIu4EXRcSrq0r1PyLi+oj4TON9N6wq2A+pfn1WRHwuIs6vvt+nNqveEbFpRPwyIm6PiEsj\nYvPGa/8eEWdXv+80YJFquaRlN4Xj7IeAO4Htq6+XUsr1pZQPllIuq557s4j4dUTcUf1zs8bX9PaI\nuKr6ev4cEe9pvFbH5j0j4q/AoRGxdkQcX8W22yLi3IhYqVo/6veuLSJeShZNti6lXFhKua/6cXIp\n5YONdddVn/8y4O6IeEhEfDQirqme+cqIeG1j/Y4RcV7j1yUido6IP1bPvH/Egx1yEfGO6uv/e0Sc\nEhEzG6+9LLLL5I7Irr7l6qyTNOUt9Z8LXebeUpPFiKllM2AV4KfND5ZS7gJOJBOy2tbAj4E1yWD0\ns4hYuUr2jgMuBdYHXgLsFhGvaPze1wA/JHfEjgW+3XqOC4BHRMT/i+wE2BYYtfUrIlYF3gL8dpm+\n2uWzDrmTuT7wTmD/iHhk9dqXgCcCGwOPr9Z8apT3eQ7wMODnS/E5twP2AVYDzgPuJv+wWIPcmXtv\nRPz3En7/24FHkR0XewBExPrACcDnyX+PewA/iYgZ1e/7AXAxWYT4HPmHkqQVM1Xj7EuBn5ZSHhjl\n/dck49E3gbWAfYETImKtasnNwJbAI8h4tl9EPKPxFuuQ36eZwE7A7sANwAyy6+zjQFnK7137uS8s\npdywFF/jm8mYvEYpZQFwDfB88s+MzwJzImLdxfz+LYFnAU8li1OvAIiIravnf1319ZwLHFm9tjb5\n39InyFh9DdldJ0nLaqn/XBjj3FtTmMWIqWVt4JYqiWq7iUV3xi+udvTvJ5PGVchW42cBM0ope1e7\nR38GvkcGtNp5pZQTq7kGRwBP41/V1dmXAVcBN/ZYs0dE3A78CVgV2HHpv9R/sWm1G1X/2HSUdfcD\ne1e7mScCdwH/We1g7QT8TynltlLKncAXWPTrbvqX73WjQ+GeiHhBY+3PSynnl1IeKKXcW0o5q5Ty\nu+rXl5FJ6QsX87UdWkr5QynlHuBHZLEEcgf2xOrfxQOllNOAucCrImID8t/lJ0sp/yylnEMm8JJW\nzFSNs2tVX99oXg38sZRyRCllQSnlSOD3wFYApZQTSinXVB1kZwOnkn/Rrz0AfLqKV/eQsXpdYGYV\nr88tpRSW7nvXtDbw1/oXEbFmFafviIh7W2u/WXV73FM9849LKf9XxdejgD8Cmyzme/ClUsrtpZS/\nAGfyYKzeGfhiKeWq6r+bLwAbV90RrwKuaPx38r/N55WkZbSkPxf6mXtLS+S5x6nlFmDtiHhIj0R5\n3er12vX1T0opD0TEDcB6QAHWqwJVbRq5k1NrJkrDwCo9PucRwDnk2dzRWoe/Vkr5xOK+oOov1Vc2\nnnXVUZZeUEp53uLeq3Jr6zmHyWA8gzxXd3Gzs5b82omIk3gwcX4PcCut73UpZbNq7Q0sWgi8vvFz\nIuLZZBfGk8lOh4eRu6ejaX+/6+/BTOANEbFV4/WVySR4PeDvpZS7G6/NAx67mM8jacmmapy9lfz6\nRrMeGWOa5pHdC0TEK4FPk91nK5Hx9neNtfNLKc3iwFeBzwCnVjH5oFLKl8i4N+r3LiLuanz8SdVz\nP6Hxtd0GrBERjyeLC03tWP1W8njKhtWHVmXxx90WF6u/ERFfb749+b1Zj0X/OykRschzSNIyWNKf\nC0v8M0HqJzsjppZfAf8kW0FHVK1YrwTOaHz4sY3XVyIH2fwfmRRdW0pZo/FjtVLKq5blQUop88hB\nOq+i1c68jO9T3y6x6mIKEf1wC3APsFHj6169/pyllFc2nmOIB7/XWy/Fe5fWr39Atl0/tpSyOnAg\ny3dG+HrgiNa/q3+rEvabgEfGolPuN1iOzyFpUVM1zp4OvLae29DD/5F/6W7aALgxIh4G/AT4GvDo\nUsoa5JGWZtxbJE6WnMOzeynlceSRlQ9FxEtYwveuLDrI+C/kv49nRcRjluZbUf+k6lr4HrALsFb1\nzJez/LH6Pa1nfngp5ZdkrG7+dxJYNJa0nPr154LULxYjppBSyh3kudZvRcQW1dnkDcnW/hvIamnt\nmRHxushBibuRyfUFwEXAndUgr4dHxLSIeHJEPGs5HumdwItbu/MDqeQ56O+R55gfBTmTYbRzyKWU\n28nv9XciYpuIWC0iVoqIjYElXXO3GnBbKeXeiNiEnAmxPOYAW0XEK6p/T6tEDoJ7TPWH0Vzgs5HX\n2j2Pql1a0vKbwnF2X3Lew+x6+GIVI/eNiKeSxYUnRsR2kcMf30R2JhzPgx1g84EFVZfEYq/FjIgt\nI+Lx1V/O7wAWkkc5lul7V0o5lewW+1lEPLuKhyuTx2UW59/I4sT86nnezvIPejsQ+FhEbFS91+oR\n8YbqtROAjRr/nexKzs+QpOU1YfJvTX4WI6aYUspXyEFZXwP+AVxI7sq8pFQ3PFR+DrwJ+DuwA/C6\n6lzuQnII18ZkZfUW4GBygNeyPss1pZS5K/DljLU9yTN0F0TEP8idwP8cbXH1vf4Q8BHgb9WP71bv\n88vFfJ73AXtHxJ3kgMwfLc/DllKuJzszPk4mzNeTVw3W/99vBzwbuI1sjx510r6kpTcV42x1vGEz\ncpbDhVX8OoMsFPyplHIr+TXtTh6N+AiwZSnllpIzeHYlY93fydh07BI+5RPIGHwX2Y3ynVLKmcv5\nvXstWRSZA9xe/b63UA2YHOXrvRL4evW5/wY8BTh/Cc882nsdQ14F+8Pqz5bLyS4aSim3AG8gj+7V\nR0qW6/NIEkzI/FuTWJTS7hDXVBd5leTjSynbj/ezSNJkZJyVJElTnZ0RkiRJkiRpTFmMkCRJkiRJ\nY8pjGpIkSZIkaUzZGSFJkiRJksbUQ8b7AZbV2muvXTbccMPxfgxJ+hcXX3zxLaWUGeP9HGPBWCxp\nEBmHJWn8LW0snnDFiA033JC5c72NRtLgiYh54/0MY8VYLGkQGYclafwtbSz2mIYkSZIkSRpTFiMk\nSZIkSdKYshghSZIkSZLGlMUISZIkSZI0pixGSJIkSZKkMdVZMSIiVomIiyLi0oi4IiI+22PNwyLi\nqIj4U0RcGBEbdvU8krREQ0Ow4Yaw0kr5z6Gh8X6iFWYsljThTLJYbByWNOGMURzusjPin8CLSylP\nAzYGtoiITVtr3gn8vZTyeGA/4MsdPo8kjW5oCHbaCebNg1LynzvtNOGTYIzFkiaSyRmLjcOSJo4x\njMOdFSNKuqv65crVj9JatjUwu/r50cBLIiK6eiZJGtVee8HwMHOBL9UfGx7Oj09gxmJJE8pee7Fw\neJiPANfVH5vgsdg4LGlCqXLii4Cv1B/rKA53OjMiIqZFxCXAzcBppZQLW0vWB64HKKUsAO4A1urx\nPjtFxNyImDt//vwuH1nSVDQ0BPPmcQHwEuC7wO31a3/5y7g9Vr8YiyUNvKoleMG8ebwV+CpwfPP1\nCR6LjcOSJoQqJ/4l8FIyJ/5H/VoHcbjTYkQpZWEpZWPgMcAmEfHk5Xyfg0ops0ops2bMmNHfh5Q0\ntVWtaOcDLwdmAOcAa9Svb7DBeD1Z3xiLJQ20Kg7fP28e2wE/IDvUdmmumeCx2DgsaeBVsfhc4BXA\nOsDZwCPq1zuIw2Nym0Yp5XbgTGCL1ks3Ao8FiIiHAKsDt47FM0ma4urBPNtvz9nDw7wCWJcMuo+t\n10yfDvvsM15P2HfGYkkDpRGH/zk8zBuBHwNfB/ZsrptEsdg4LGngNGLxmcPDbEFWTc+q/gl0Foe7\nvE1jRkSsUf384cDLgN+3lh0LvK36+TbAL0op7TN0ktRfjcE8pwGvBGaSQXf95rqDDoK3vGUcHrB/\njMWSBlIjDt8LvA74GfAt4EPNdTNnTvhYbByWNLAasfgU4FXAv5M58XrNdR3F4Yf0/R0ftC4wOyKm\nkUWPH5VSjo+IvYG5pZRjgUOAIyLiT8BtwLYdPo8kpWowz4lkAvyfwOnkEY0RM2dO6OS3wVgsafBU\ncXiYnNx4Bnk2eafmmpkz4brrxuHh+s44LGkwVbH4OLIK+iTgNGDt5poOc+LOihGllMuAp/f4+Kca\nP78XeENXzyBJixgayqA7bx7HAG8CngqcQmtK2ORqCTYWSxos1YC0O4EtgfOAQ3mwLQAwDktSlxo5\n8dHAm8kgdTKwZnNdx7F4TGZGSNK4a7ShHUlmfM8kOyIWKURMgpZgSRpYVSy+nRwafD4wh1Yhwjgs\nSd1p5MRD5ObcJmRHxCKFiDGIxV0e05CkwTA0BG97GyxcyGHAO4Dnk9fGrVavmT7d5FeSutLYhbuV\nnNR+GTmw8rX1GuOwJHWrkRMfArwb2JwcWrNqvWYMY7GdEZImt7r6u3Ah3wHeTt6bfBKNQoS7cJLU\nncYu3F/JxPdy4BgahQgwDktSlxo58beAd5EdaifQKESMcU5sZ4Skya0azPM14MPAa4CjgFXq1yfP\ngDRJGkxVHL4BeAl5h+WJwIubaybP0GBJGkxVLP4y8FHgv4EfAg+rXx+HnNjOCEmTU3Vncpk3j8+S\nhYg3AUfTKERMogFpkjRw6rvr583jGvJ43F+BU2kVIozFktSdRk78KbIQ8WbgRzQKEeMUh+2MkDT5\nVG1oZXiYDwNfJ49nfA+YVq+ZNs2WYEnqSt0OPDzMVeTxuH8CvyCHB4+YOTMTYGOxJPVfIyfeHdgP\neCd5lfIg5MR2RkiaPOpduO2354HhYd5HFiJ2AQ6mEXSnT4fZs01+JakL9YC04WF+C7wAeAA4i0Yh\nYvp0mDMnW4KNxZLUX42ceOHwMO8hCxEfpLU5N845scUISZNDY0Da/cBbgQOBPYFv0gh2DquUpO40\nBqSdD7wIeDhwDvDkeo1xWJK608qJtycLEB8nCxJRrxuAWOwxDUmTQzWU55/kbIifA18APtZc47BK\nSepWFYtPB7YG1gdOBzaoXzcOS1K3qjh8L/BG4DjgS+QG3YgBicV2Rkia+IaGYN487ga2JAsR36JV\niHBAmiR1pzGs8ufAq4H/IDsiRgoRxmFJ6laVE99FxuHjgP1pFSIGKBZbjJA0sVWtaH8HXkYORzuU\nnBMxYgDa0CRp0mq0BB8BvB7YmJwRsU69xqHBktStKhbfRg4NPhs4HHhfc82A5cQWIyRNTI3BPH8b\nHuZFwMXAj4Ed6zUOSJOk7jTiMMPDfIec1/NC8mjGmvU6hwZLUncasfivw8NsDvyWvM5+h3rNgObE\nzoyQNPE0roybR3ZE3AgcX/18xABVfiVpUmnE4ULO6PkEsBV5d/0q9Tqv7pSk7jRi8XVkHnwTcCLw\nkua6Ac2JLUZImniqwTy/J4PuXcBpwGbNNTNnDmTQlaRJoYrDBfgweY3y9sD3gZXrNQMyIE2SJq0q\nFl9J5sT3kJ1pmzbXDHBO7DENSRNHY0DaxcDzgfvJM3GLFCIGaDCPJE061YC0BcC7yULE+4HZNAoR\nxmFJ6k4jJ/418ALgATInXqQQMeCx2GKEpImhMSDtLPLu+n8DzgWe2lw3YIN5JGlSqWLxveQ1yoeQ\nxzO+RSOpNA5LUncaOfEZwIuB1YDzgKc0102AWOwxDUmDb2gI3vY2WLiQn5MJ8H8Ap5J32ANZ+R3w\ngCtJE9bQULYDz5vHncBrgTOA/YDd6jXGYUnqViMnPgbYFngicAqwXr1mAsViOyMkDba6+rtwIYfx\n4JVx59AoREyAyq8kTViNXbhbyKFoZ5HHMnZrrjMOS1J3Gjnx94FtgGeSRzNGChETLCe2M0LSYGrs\nwkGeSd6DHM7zU2DVep0D0iSpG604fD3wcuA6Mg6/prl2gAekSdKE1orFXwH2BF4B/IQ8tgxMyJzY\nYoSkwdO6Mu7jwJeANwKHAw+r1w34UB5JmrAacRjg92Qh4g6yHfgFzbXGYknqRisn3hP4Knk8Yzbw\n0HrdBI3DHtOQNHiqa4oWAO8iCxE7Az+gUYiYNm1CtaFJ0oRSxWGAi4DnAfeRxzMWKURMsJZgSZpQ\nGjnxO8hCxPuAIRqFiAmcE1uMkDQ4GtcU3Qu8gbyz/tPAd4Bp9brp02H27AkZdCVp4FVXdwKcRk5q\nXx04H3h6vWb6dJgzJ1uCjcWS1F+NnPge4HXAYcBngW/T+Ev8BM+JLUZIGgyNAWm3k+fgfk5eF/cZ\nIOp17sJJUnfqWAwcCbwaeDx5Zdx/1GuMw5LUnUZO/HfyiNzx5Mbcp5hcObEzIySNv8Y1Rf8HbEGe\nTz6SvMYTmFDXFEnShNMakPYN8qaMF5KF4dXBOCxJXWvkxDeSOfEfgKPIjmFgUsXizjojIuKxEXFm\nRFwZEVdExAd7rNk8Iu6IiEuqH5/q6nkkDajGNUV/AJ4L/Bk4gUYhYhJUfseDcVjSUmnswtVDg3cD\nXgucTFWIAOPwcjIWS1oqjZz498Bm5O1FJ9EoREyynLjLzogFwO6llN9ExGrAxRFxWinlyta6c0sp\nW3b4HJIGWTWY5yKyHTiAM4Fn1a9PwGuKBohxWNKSVXH4fmAn8lzyTrRm9Xh154owFktasioWXwBs\nScbfs4Bn1q9Pwpy4s86IUspNpZTfVD+/E7gKWL+rzydpgmkM5jkZeBGwGjkgbaQQMUGvKRoUxmFJ\ni9WIw3eTnRCHkUODD6Q1NNhYvNyMxZIWqxGLTySHBq8B/JJGIWKSxuExGWAZERuSA5gv7PHycyLi\n0og4KSI2GuX37xQRcyNi7vz58zt8UkljotESPBvYCngiGXSfUK+ZwNcUDaIVjcPVexiLpcmiEYdv\nAV4CnAgcgEODu2ROLGkRjVh8KPAa4P+Rm3MjQ4MncU7ceTEiIlYFfgLsVkr5R+vl3wAzSylPI4fm\n/6zXe5RSDiqlzCqlzJoxY0a3DyypO3Xld/vtKcPDfAnYkRyQdjawTr1ugl9TNGj6EYfBWCxNGvWA\ntOFhriVn9VxKBomd6zVe3dl35sSSRrRy4n2Ad5BdEWcBj67XTfKcuNNiRESsTAbdoVLKT9uvl1L+\nUUq5q/r5icDKEbF2l88kaZw0Kr8LgV2BjwHbkbtxj6jXuQvXV8ZhSYtoDEi7hByQNh84nTymARiH\nO2AsljSilRPvAnwC2J68wnO1et0UiMWdDbCMiAAOAa4qpew7ypp1gL+VUkpEbEIWR27t6pkkjZPG\nNUX3kMH2p8DuwFdoVEUn4WCe8WQcljSidXXnacDryXPJZwBPqtcZh/vOWCxpRCsn3o5sg/oI8EWm\nXk7c5W0azwV2AH4XEZdUH/s4sAFAKeVAYBvgvRGxALgH2LaUUjp8JkljrbELdxuwNXkObj/y6rgR\nk3QwzzgzDkt6MA4PDwMwB3g7WYA4kcYkReNwV4zFkhbJiW8l50P8Cvgm8IHmuikUi2OixblZs2aV\nuXPnjvdjSFpa1XTgecArgWuAI4A3NtfMnJlBd4K3oUXExaWUWeP9HGPBWCxNIFUcLsCXySNyLwKO\nAVav10ybNinOJRuHJQ2sKhZfS+bE15HF4W2aa6ZYTjwmt2lImoIa1xT9FtgUuAk4lUYhwgFpktSd\nRhxeCLyfLES8GTiJRiFikg9Ik6Rx1YjFFwPPAW4mj8uNFCKmaE7c5TENSVNVoyX4FDLQPpIckDZy\nV9kkvqZIksZdIw4PkwWIYxnlXPIk2IWTpIHUiMUnAW8A1gLOJK/wBKZ0TmxnhKT+aVxTxPAwhwKv\nBh5HnokbKUS4CydJ3Wlc3XkzeSTjOPK+yC9TJX9TdBdOksZEKyc+GNgKeAKZE48UIqZ4TmwxQlJ/\nNK4pKsBnefC+5HNpDEibAtcUSdK4aQxI+wN5dedl5A1Gu9RrjMOS1J1WTvwp4N3AS4FzgPXqdcZi\nj2lI6oPGNUX3Ae8BDgPeBnwPWLleN0WuKZKkMde6uvOX5KT2INuBN63XGYclqTutnPhd5OD2dwAH\nYk7cZmeEpBXT2IW7gzyWcRjwGeBQGkF3Cl1TJEljqrELB3A02ZW2JtkOPFKIMA5LUncaOfHt5I0Z\nRwCfAw7GnLgXixGSlk/rLNxfgOcBZ5HFiE+TO3KAbWiS1IVWHC7A18gBac8guyMeX6+dwgPSJKlT\nrVg8D3gueUz5cOATmBOPxmMakpZdYzIwwG+ALYG7gZOBl9Trpk834EpSF1pxeAGwK3AAWYyYDTy8\nXmsslqRutGLxXHJQ5T1kTvziep1xuCc7IyQtu732Ggm6xwMvIFvPzqdRiHAXTpK604jDdwJbk4WI\nDwM/pFGIcBdOkrrTiMXHAi8EHkZ2po0UIsyJR2UxQtLSq9vQqnPJ+5MJ8H8CFwBPrtdN8WuKJKlT\nQ0MjcfhGsiB8MvAd4Ct4dackda6VE38T+G/gSWRO/KR6nTnxYlmMkLR0GgPSFgIfIq+JezV5TdG6\n9Tp34SSpO3UsBi4Fng38iexSe2+9xjgsSd1p5cQfrH68hpydtk69zli8RM6MkLRkjWuK7ga2I1vR\ndgX2BaaBZ+EkqUutqztPBN4ErA6cBzwNjMOS1LVGTnwX8GayGPw/wFcxJ15WdkZIWrzGNUX/R7YD\nHw98C/gGVdC18itJ3Wld3fltckDaE4ALqQoRYByWpC41cuL6iNyJ5BG5kc05c+JlYmeEpN5au3CX\nkDdm3EF2Rby6XjdzZp5JliT1X2MXrj4i902yGPEDYNV63cyZJr+S1IVWTvxbMif+B7lB98p6nTnx\nMrMYIelfta4pOh7YFngkjXZgyDa0ffYZl0eUpEmvsQt3J9kOfAKtdmAwFktSV1o58bFkLF6LvEXu\nqfU64/By8ZiGpAfVk4G33x6GhynkUYytgf8CLqJRiPCaIknqRisWXw88nwdvzBhpBwZbgiWpCz1y\n4n3JGzM2InPikUKEOfFyszNCUmpVfu8nB1QeCLwWOAL4t3qtg3kkqRutWHwRWRAeJrsiXlGvMw5L\nUjd65MS7AAcBrwcOB6bXa43FK8TOCEkPnkmugu7t5EyIA4E9gaNpFCLchZOk/mvtwgH8GHghsArw\nSxqFCHfhJKkbrZz47+RMiIOAjwE/olGIMCdeYXZGSFNd40wy5H31W1X/PAR4R73Oyq8kdaO1C1eA\nfYBPApsBPwNm1GuNxZLUjVZO/EdyUOW1wKHAjvU643Df2BkhTXV77TWSAJ8NPBu4GTidRiHCyq8k\ndacRh+8FdiALEdsDZ9AoRBiLJak7jVh8JpkT30rG4R3rNcbhvrIYIU1VdUtwdU3RIcBLyaT3QrI1\nmOnTYc6cvKbIoCtJ/Tc0NBKH/wpsDgyRnRGHk0c0jMWS1KFWTnwQ8HJgHTInfj4YhzviMQ1pKmq0\nBC8EPkJOCH45cBSwBngmWZK6Vsdi4FLyiNytwE+A19VrZs7M6+KMxZLUf42ceAGwB3mT3BbAD4HV\nwZy4Q3ZGSFNJa0DaHcBryELEB8hJ7WtAVn9nzzboSlIXWrH4Z8BzyVkR51EVItyFk6Tu9MiJtyIL\nEbsBx1EVIsyJO2VnhDRVtAakXUMG3T8CBwA71+vchZOk7jRicQG+COwFbEIOqly3XucunCR1o5UT\nN4e3fxfYqV5nTty5zjojIuKxEXFmRFwZEVdExAd7rImI+GZE/CkiLouIZ3T1PNKU1rqm6CxyKM9f\ngVNpFSLchZtUjMXSgGjtwt1DDqjcC3gLOUB4pBAxc6ZxeBIxDksDpJUT/4IsBs8nh7cvUogwJ+5c\nl8c0FgC7l1KeBGwKvD8intRa80rgCdWPncgNWkn91Lqm6ADgZcCjgIuAF9Xrpk/P6q8mG2OxNN7q\nOFwNR7uRHBJ8JPAF4AiqQZVgLJ6cjMPSIGjkxAXYn5yXti6N4e1gHB5DnRUjSik3lVJ+U/38TuAq\nYP3Wsq2Bw0u6AFgjItZF0opr7cLdD7wfeB8ZeH8FPL5e6zVFk5axWBpnrV24i4Bnkf8j/gz4GBD1\nWmPxpGQclsZZKye+D3gvsAtZBfwV8B/1WuPwmBqTmRERsSHwdLLo1LQ+cH3j1zdUH7up9ft3ouqa\n2WCDDbp6TGnyaJ2FuwV4A3k8Yw/gS8A0yMqvAXfKMBZLY6zVmTYHeBe5C3cK8JR6nbF4yjAOS2Os\nlRPPB7YwdrZZAAAgAElEQVQBzgH2JK9RNiceP53fphERq5K3VO1WSvnH8rxHKeWgUsqsUsqsGTNm\n9PcBpcmkVfkFuJw8C/cr8s76r1IFXa8pmlKMxdIYasXi+grlHcge/V/TKES4CzdlGIelMdQjJ76M\n7Ey7iCwOj2zOmROPm047IyJiZTLoDpVSftpjyY3AYxu/fkz1MUnLqlX5BTiGTH4fQQ5He3b9gtXf\nKcVYLI2hViy+HdgOOIkcFvxNYGUwDk8xxmFpDPXIiX8CvJW8wv4csigBGIvHWZe3aQRwCHBVKWXf\nUZYdC7y1miC8KXBHKeWmUdZKGk3rTPIDwN7kXfUbAXNpFCLchZtSjMXSGOmxC3c12QlxGjmN8ACq\nQoS7cFOKcVgaQz1y4k+TRzOeSubEI4UIc+Jx12VnxHPJTdnfRcQl1cc+DmwAUEo5EDgReBV5resw\n8PYOn0eanFpnku8E3kZ2RbyVvC95FbDyO3UZi6Wu9diFOxF4M/Aw4AzgBfULxuKpyDgsjYUeOfEO\nwM+BHcmCsDnxYOmsGFFKOY/GgOhR1hRywL+k5bXXXiMJ8DXkOO6rgH2B3aj+J5w5M68oMuhOOcZi\naQw04nAhzyHvBWxMFoZn1uuMxVOScVgaI41Y/CcyJ74a+F9gV8yJB1HnAywldaRuCa7urT+VbDu7\niZzS/j9ATJ8Oc+bAddcZdCWpC0NDI3H4buBN5Jb3tsB5VIUIY7EkdaeVE59M5sR/I/PjD2JOPKgs\nRkgTUd2GNm8eBfgaeU/yY8gp7S8FzyRLUtfqWAz8GXgOOSTtK8AQMB08kyxJXWrlxF8mzzvNJHPi\nF4M58QCzGCFNJK0BacPAW4APA68nr+98HOQu3OzZBl1J6kIrFp9O7sLdQN6a8WHchZOkTrXi8N3k\nnJ6PAm8Ezgf+HcyJB5zFCGmiaFR+Aa4FNgN+CHwBOAr4N3AXTpK61KMz7RXAeuTd9S+v1xmHJakb\nrZz4z2RO/GNyZs+RmBNPFF3epiGpX+priqrpwKeT55IfAE4gj2gAGXSvu248nlCSJrehoRyO1pgP\n8S6yILwNcCiwar125kyTX0nqQisnPpWc0QPZmTZSEDYnnhDsjJAGXeOaokKeRX4FsC55Fm6kEDF9\nek4HliT1V2sX7hpyPsRRwBeBH9EoRBiLJakbrZz4S8AWPDgzbaQQYRyeMCxGSIOqdRbuLrLyuyc5\nH+IC4PH1WtvQJKkb9S5cdV3cycAscj7EyeT55JE7G43FktR/rZz4TuANwMfITuFfAf9RrzUOTyge\n05AGUV35rZLfPwKvBa4ipwR/mCr5nT7dgCtJXWnswj1AdkF8EngKcAzVwGAwFktSV1o58dVkTnw1\n8FVgd8yJJzI7I6RB0qr8AhxH7sL9FTgF+AhV0PWaIknqRisW3wG8DvgEOa195OYicBdOkrrQIyf+\nObAJMB84DdgDc+KJzs4IaVC0Kr8Lgc8AnweeQd5dv2G91uqvJHWjFYuvIAsR1wDfAD6Au3CS1Kke\nOfGnyNvjZpE58Qb1WmPxhGZnhDQIWmeSbwO2JAsROwLn0ShEuAsnSf3XYxfuR8CzgTuAM4FdcRdO\nkjrVyolvBV5FFiLeCZxLoxBhTjzh2RkhjbfGmWSA35K7cDcCBwDvwV04SepUaxfufnIw5b7krRlH\nA+vVa43FktSNVk58MTm0/SbgIODd9Trj8KRhZ4Q0Xnrswh0GbAYsICu/O1MVIqz8SlI3WrtwfwVe\nShYiPgCcRaMQYSyWpP7rkRMfAjwXeIDMiUcKEcbhScXOCGk8tHbh7iXbf78HvBj4ITADrPxKUpda\nu3DnAW8kj2UMAdvV64zFktSNHjnxLmQx4mXAD4C1wTg8SdkZIY2lHpXf64DnkYWIj5I3ZswAzyRL\nUldasbgA+wGbA6sCF9AoRLgLJ0n91yMnvpbshjgE+DhwElUhwpx40rIzQhorrcovZJDdnpwS/DNg\n6/oFq7+S1I1WLL6THIr2Y/Lu+kOB1cE4LEld6ZETnwDsABTyWvst6xeMxZOanRHSWGidSV4IfBp4\nNfAYYC6NQoS7cJLUfz124a4AngX8FPgKeV3c6uAunCR1pUdO/Emy+DCTzIlHChHmxJOenRFS11pn\nkueT3RCnktd27g9MByu/ktSVHrtwPyAHoq0GnAG8sH7BWCxJ3WjlxDcDbwFOJzvUvgU8HIzDU4id\nEVJXeuzC/Qp4BnA2eUXR96kKEVZ+JakbrV24fwLvIxPgZwK/oVGIMBZLUv/1yInPJ3Pi88gZEQdT\nFSKMw1OKnRFSF1q7cAX4JrAHsAFZlHg6WPmVpC61duGuJW/LmAt8GPgCVSJkLJakbvTIifcD9iSP\nZfwK2BiMw1PUMnVGRMQjI2KjiHhcRNhVIfXS2oW7A3gDsBvwKuBiqkKEZ5K1HIzD0lLosQt3HLkL\n90fgGHJGxEPAXTgtF2OxtAQ94vDtwOuB3YGtyJx4YzAnnsKW2BkREasD7wfeDDyUPPK+CvDoiLgA\n+E4p5cxOn1KaKFq7cJeQhYhrga8BHwICrP5qmRiHpWXQ2oW7H9gL+CpZjPgx8DgwDmuZGYulpdRj\nTs9vyJz4L8C+5CadObGW5pjG0cDhwPNLKbc3X4iIZwI7RMTjSimHdPGA0oQwNAR77QXz5gHZgvY9\nYFfyfuSzgOfVa2fOhH32MehqWRiHpSVpxWGAG4BtybPJO5OtwauAu3BaXsZiaUnqDuFqY64A3wU+\nCDyKnJu2Wb3WnHjKW2IxopTyssW8djHZYSNNXa3q711k0jsEvByYA8wAK79absZhaQl67MKdQt5c\ndA95c8ab6xeMxVpOxmJpCVodwncC7wGOBLYAjiA36YzDqi31GbeIeGfr19Mi4tOLWf/9iLg5Ii4f\n5fXNI+KOiLik+vGppX9saQD0OAt3OXln/ZHA54CTqAoR7sKpD5Y1DldrjMWa3FpzehYAnwBeCaxD\nDqscKUQ4H0J9YE4stfTIiS8DZgFHAfsAJ1AVIozDaliWgTsviYgTI2LdiNgIuIC8nns0h5FFsMU5\nt5SycfVj72V4Fml81ZXfxrGMQ4FNyOE8p5PJ8EqQ1d/Zsw266odljcNgLNZk1tqF+z/gpWTi+w7g\nQuC/IOPwnDlw3XXGYvWDObFU65ETHwI8m+yM+AXwcWAl47B6WOqrPUsp20XEm4DfAXcD25VSzl/M\n+nMiYsMVfkJp0LTOwt1F3ll/BPAish14nXqtZ+HUR8sah6vfYyzW5NNjPsSp5LGMu4HZwFvrF4zD\n6jNzYqnSyonvBN5LHlV+SfXPR4MdwhrVshzTeAI5e+QnwDxySM/0Ffz8z4mISyPipKqyPNrn3iki\n5kbE3Pnz56/gp5RWQGsX7nfksYw5wGeA06gKEVZ/1YGO4jAYizWRtHbh6mMZW5DD0X5NVYgwDqsj\n5sQS/5ITX0oeyzgS2Juc2/NosENYi7UsxzSOAz5ZSnkP8ELyqu5fr8Dn/g0ws5TyNOBbwM9GW1hK\nOaiUMquUMmvGjBkr8Cml5dQ6C1fflrEJ8HfyWMangWngWTh1qd9xGIzFmih6nEm+AXgxeSzj7cBF\nwJPAXTh1zZxYU1ePnPi7PHgs4wzgk5gTa+ks9TENYJNSyj8ASikF+HpEHLe8n7h+r+rnJ0bEdyJi\n7VLKLcv7nlInWlPa/0HelnEkeTZ5Do3KrwFX3eprHK7ex1iswdfjtowTyQ6Ie8ljctvXLxiL1T1z\nYk1NrVh8B3lbxlHkDXJHkB1qxmEtrSV2RkTE82DRQFkrpfwhIh4REU9e1k8cEetERFQ/36R6lluX\n9X2kzvTYhbsYeAYZdD9PowXNXTh1qKs4XL23sViDrXVbxn3AHsCrgfXJuDxSiHAXTh0yJ9aU1SMn\nnkvmxEcDXyBvkHsUmBNrmSxNZ8TrI+IrwMnkn/nzgVWAx5Pz+mYCu7d/U0QcCWwOrB0RN5Bd7CsD\nlFIOBLYB3hsRC8hrwLetqsvS+GtVfgvZN7kHWXw4G3hevdbqr7q3XHEYjMWa4Fpnkq8FtiWPY7wP\n+Dr5P4JxWGPEnFhTT4+c+BvAR8g5aWcDz63XGou1jGJpYl1ErAm8nvxvbR0yUF4FnLCkSe79NmvW\nrDJ37tyx/JSaalqTgW8lzyIfB2xFXuG5Vr3WKe1qiIiLSymzOnrvgYnDYCxWx3rclvFj4F1AAAeT\nf3sDjMNaRJdxuHr/gYnFxmF1rpUT30LmxMcDryFz4jXrtcZiNSxtLF6qmRGllNsiYj1yU+La+sPk\nrS1jngRLnWntwp0LbAf8DfhfYFcyEbbyq7FmHNaU0dqFuwfYDTiIHJD2Q2BDMA5rXBiLNWW0cuJz\nyJx4PvBNYBfMibXiluU2jbsaPxYAr6TKB6QJr3UWbiF5LdHmwMOAX5F3eAV4JlnjyTisyavHmeQr\nyOuTDwI+TBaINwTPJGu8GYs1ebVi8QLy+voXAdPJnPgDmBOrP5b6No1Sytebv46Ir5Hz+6SJrbUL\ndwM5DO3s6p/fAVYDK78ad8ZhTVo9ziQfRHZEPII8oP+Keq2xWOPMWKxJqxWLrydz4XPI24u+jTmx\n+mtZOiPapgOP6deDSGOuxy7cz4GnkROCDyOvKFoN3IXToDIOa+Jr3ZZxGzkPYmfg+cClNAoR7sJp\nMBmLNbH1yImPIXPi3wCHA7MxJ1b/LXVnRET8jtysAJgGzCA72aWJp8eZ5A8D+5PXFB0JPLFea/VX\nA8I4rEmndSb5PPJM8k3AV8hrCVYC47AGirFYk0qPnPhDwIHALDInfny91lisPlvqYgSwZePnC4C/\nlVIW9Pl5pO61JgNfQV4VdzkZfL9AzokAnAysQWMc1uTQui1jAfB54HPAvwO/JGdFAMZhDSJjsSaH\nVk78O+DNZG68B7AP8NB6rbFYHViWmRHzlrxKGlDNxDcCSqGQVd8PkWeSTwK2qNdb+dUAMg5rQusR\nhwHmkWeSzwN2IDvUPJOsQWYs1oQ2Sk68P1mAWIMcgPLyer2xWB1als4IaWJqtZ9RCrcA7wSOJc8i\nzwYeXa+38itJ/dUjDgMcBbwHeICc0bN9vd4zyZLUfz1i8XzgHcDx5LUwhwGPqtebE6tjKzLAUhps\nPYbxAJwOPJWczr4fcCJVIWL6dJgzB667zqArSf3SGlAJcCfwdvKI3H8Bl9AoREyfDrNnG4clqV9G\nyYlPI3PiU4FvACdQFSLMiTVGLEZocqorv/Me7KS8jxxS+TJgdeBC8tq4lcBdOEnqQmtAJcBFwNPJ\n6eyfAM4FHle/6G0ZktRfPXLif5IDgl8OrAn8GtgVCDAn1piyGKHJp8cu3O+BTYGvkdfFXQxsXL/o\nLpwk9VePXbiF5IDg5wL3A2eRAytXBnfhJKkLPXLiK8mceF/gfeR19k+tXzQn1hizGKHJpbULV4Dv\nktd1/gX4GXAAMD0i17sLJ0n91WMX7i/Ai4C9gNcDlwLPNw5LUnd65MTfAZ4J3EDOTdsfeLixWOPI\nYoQmhx67cPOB/yY7IZ5PXle0NWSwPeKIHKDmLpwk9ccoZ5J/SO66XUIezTgSWGPaNOOwJHWhRyy+\nGXgN8H7ghWROvBWYE2vceZuGJr72ZGDyms63A38n29A+CKzk1USS1I0ecfgOMvEdAp4DzKGaDWEs\nlqRu9IjFJ5C3ZdxBDqncBXNiDQ47IzRx9aj8DgMfAF4FzCDPwf0PsJLDeCSpGz3OJJ8DPI3sivhs\n9evHgW3AktSFUXLi9wNbAuuQOfGumBNrsNgZoYmpR+X3t8BbgKvIWzK+CKwC7sJJUldaZ5LvAz4N\nfJksPpxHDkozDktSR3rkxBeTOfHVwIfI4cEPA2OxBo6dEZp4WrtwC8nCw7PJFrTTgP2oChHuwklS\n//XYhasntH+JbAm+pPq1cViSOtIjJ96HjL13AacDX6cqRBiLNYDsjNDE0tqFuxZ4K7n7tg1wILAW\nWPmVpK60duEeICeyfwRYFfgp8FowDktSl1o58Z+BHYBfAm8ib85YE4zFGmh2RmhiaO3CFeBQckL7\nZeSE9h9RFSKs/EpS//XohrgR2II8h/xickL7awE8kyxJ3eiREx9Czum5ghwafCRVIcKcWAPOzggN\nvtYu3HzgPcAx5PVEs4GZYOVXkrrS40zyUcB7gX8CB5BxOcBYLEldacXim4GdgJ8DLwIOAzYA47Am\nDDsjNLh67MIdDzyFvKboa8AvqAoR7sJJUjdaZ5L/DmwHbAs8kZwNsTNVIcJdOEnqvx458bFkTnwy\neY396VSFCHNiTSAWIzR4hoZg7bUz4M6bB8CdZOV3Kx68nmh3qv+Ap0+H2bMNupLUL3XiGwE77DBy\nJvl0Mvn9MbA3Oa/nCZBxeM4cuO46Y7Ek9UuPnPgfwLuArYH1yJsz/gdzYk1MFiM0WOr2s1tvHfnQ\nueQ5uIOBPYELyWQYcBdOkvqtjsNV4kspDAMfAF4GrAZcAHyS6qyncViS+q9HTnwOmRMfCnyMzIk3\nql80FmsCshihwdCj/execjr7C8n233PIK+NG7kl2F06S+qdHHIZMdp8OfBvYDfgN8EwwDktSF0bJ\nifcANgemkRt1XwAeCsZiTWidFSMi4vsRcXNEXD7K6xER34yIP0XEZRHxjK6eRQOuvQsH/BaYBXwV\neDdwKfC8+kUrv9JSMxZrqfSIw/cBewGbkYnwL4D9gIeDZ5KlZWAc1lLrEYsvJgvAXycHBV9CxmXA\nnFgTXpedEYeRN36N5pXkUdMnkOMADujwWTSIelR+7wc+B2wC3AacCHyXvLveyq+0XA7DWKzFaQ2o\nhLyi89nkztvbyCuUX1S/6JlkaVkdhnFYizNKTvxZYFPgdnJQ5QGYE2ty6awYUUo5h/z75Gi2Bg4v\n6QJgjYhYt6vn0YDpUfm9kqz0fgp4A3A5+aczAGutZeVXWg7GYi1WHYurAZULgC+Su3A3kdfFfR9Y\nPSLXuwsnLTPjsBarR058BVmE+AzwJjInfkX9ojmxJpHxnBmxPnB949c3VB/7FxGxU0TMjYi58+fP\nH5OHU4dau3ALydazZwDXkVPafwCsCZn4zpkDt9xi0JW6YSyeinrswl1NHof7OPBaMvl9DWQcPuII\nKMVdOKkbxuGpqkdO/BUyJ74e+AkwB3gkmBNrUpoQAyxLKQeVUmaVUmbNmDFjvB9HK6K1C/dHckDl\nHmT/4uXANmD7mTSAjMWTRGsX7gFyFsTGZEw+EjgKWNs4LA0c4/Ak0sqJ/wA8n7w57tVkTvw6MCfW\npDaexYgbgcc2fv2Y6mOajFq7cA8A3yKvJ7oCOAI4Bng02AYsjS1j8VTS2oW7hpzO/iHy2s4rgG3B\nAZXS2DIOTyU9cuJvkAXh3wNDZEfEo8CcWJPeeBYjjgXeWk0Q3hS4o5Ry0zg+j/qtDrYRsMMOI7tw\n1wIvAXYlk+DLge2BsPIrjQdj8WTXjsULF/IAsD/wVHI45aHkfIh1wAGV0tgzDk92o+TEfyaHA+8G\nvJjMibfDnFhTx0O6euOIOJL8u+baEXED8GlgZYBSyoHkRQmvAv4EDANv7+pZNA7q1rN6OnspPEDe\njPFhsgp2CPkvPcBdOKkjxuIprkcsvhZ4B3AWORDtYHIbFshduH32MRZLfWQcnuJGyYkPAD5C/mXs\nUPLmInNiTTWdFSNKKW9ewusFeH9Xn1/jqG4Drs7AQQ6mfBdwBtkKfDCwQf3i9OkGXakjxuIpamgI\n9tprkensDwAHkTN6VgK+B7yTKvk1DkudMQ5PYT1y4mvJ2HsmWRD+Ho0zOsZiTTETYoClJpDWMJ4C\nHAg8BbiQ7Iw4hUYhwrNwktRfPa6Ju44sBL+XvEL5crJA7C6cJHWklRM/AHyHzInnkkWIk2gUIsyJ\nNQV11hmhKabHLtx1LNoN8T1gZv2ilV9J6r/WLlz7eNxBNIoQYCyWpH7rkRP/mYy9dTfEQdghLIGd\nEeqHHtfE7Q88GbiIDLinADOjSn+t/EpS/7V24f5MDgt+H/Bcshvi3UAYiyWpGz1y4m/xYDfEwWQ3\nxAbGYQmwM0Irokfl9xryHNzZwMvJbogNINuAnc4uSf3XisV1Qfij5B/yB5MDKwMcUClJXeiRE/+J\nzInPAbYgN+ceC+bEUoPFCC27oSH44Afh1ltHPrSQrPx+nBwP7WA0SepQM/GNgFIA+AMZe8+jlfwa\nhyWp/0bJib8BfAJ4KK3b44zF0iIsRmjZtK8nAn5P7rr9iryX6rt4TZwkdabHNXELgP2ATwGrAIcB\nb8UBlZLUmR458ZVkQfgCYEtyiPv69YvmxNK/cGaEls7QEGy4IWy//UjQXQB8CdgYuBo4HDieqhAx\nfTrMmQPXXWfQlaR+qQdUNpLfK8gbMj5CDka7ksZ99dOn2w4sSf3UIye+H/gC8HTgj8Ac4FiqQoQ5\nsTQqixFash7XxF0KPBv4GLAVmfzuQONMsrtwktRfrQGV9wF7k8nvtcAPgWOAdev1xmJJ6q8eOfFv\nyZx4L+C/yZz4LZgTS0vDYoRG16Py+0/gk8As4EbgaODHwKPByq8kdaFHLP41GYc/DWxDJr9votEN\nYSyWpP7pEYfvJQsQzwJuAn4CHAU8CozD0lJyZoR663EO7pfkHclXkWeR9wPWrF9cay34xjcMuJLU\nT61YPEwWIPYF1iHbgLeCB4dYeiZZkvqrR058Pjkb4mpgRzImP7J+0ZxYWmp2Ruhftc4k3wV8EHge\ncDd5P/JsqkLEzJlZ+b3lFoOuJPVLj124M4GnAl8jC8NXUhUipk2DI47IYoS7cJLUP62c+E7gA8Dz\nyc6IU4BDqQoR5sTSMrMYoVQnvhGwww4jZ5JPAZ5MXtv5PuBy8ro4288kqQNDQ7D22lmEqM4k3w68\nG3gxUIBfkLcWrQ4OqJSkfhslJz4J2AjYH9iFzIlfDubE0gqwGKF/HcZTCreSRzG2AB4OnAt8G1gN\nHMYjSV2oY3HjvvpjgCcB3wf2AH4HvKh+0VgsSf3VIye+hRzS/ipgVfKIxjernxuHpRVjMWIq69EG\nXIAjgf9X/XMvckrwc8HKryR1oUcsvgl4PfA6chjahcBXgelgLJakfhslJx4ic+KjyAHuvwWeA8Zh\nqU8cYDlV9RjGMw94L9mGtglwBvCU+sVp06z8SlK/tWLxA8DBwEfI88hfBHYHVq7XO6BSkvqrR058\nHbAzeVz52WRcfnL9ojmx1Dd2RkxFrWE8C4FvkOfgziFvyfgljUKEZ5Ilqb967MJdTR7BeA/wdOAy\n4KNUhQh34SSp/3rkxPuROfF5ZH58Po1ChDmx1FcWI6aauvpbDeO5lGw32w14AXBF9fNpEbnes3CS\n1F+tM8n3AZ8jb8q4jNyB+wXwxHr9WmsZhyWp31o58SVkF8SHyMLwlcCumBNLXbIYMVW0duHuAT4G\nPJNsRfsBcAIwEzLYek2cJPVfaxful8AzgE8B/w1cRd5dH+A1cZLUhVZOPAzsCcwCrgd+CBwHbADm\nxFLHnBkxmQ0NwV575e5bRAZS4HTyHNw1wI7knfVrQbaeWfGVpP7qEYvvIAvCBwKPAY4FtqrXG4sl\nqb9GyYlPJXPia4F3kIOC1wTjsDRGLEZMVu1hPKUwnxyEdgTwBLINeOSKOIfxSFL/tWJxKYWfkq2/\nN1X//BzVtcnggEpJ6rceOfHN5HGMIfJI3FnAC+v15sTSmPGYxmTUagMuwKHAf5GtZ58gzyWPFCIc\nxiNJ/dVjQOVfgK2BbXjwus7/pSpEOKBSkvqvlRM/ABxC5sQ/Io/IXUqjEGFOLI0pixGTTWsYz+/J\nosM7yHuSLyF34Vap1zuMR5L6qzWgcgE5nf1J5JXJXwV+DTyrXu8unCT1XysnvgrYHHgXeTvGpcBn\nMSeWxpPFiMmitQt3L/Bp4GlksD2IvLbzSfV6d+Ekqf9au3C/BjYh24FfSN5YtAeNM5LuwklSf/UY\n2v5JMie+nLyx6Cxykw4wJ5bGkcWIiawOthGwww4ju3BnkFfE7Q28nuyOeDewklcTSVL/tWPxwoX8\nA/gAeU3cX4EfA8cDG0KuA2OxJPXLKDnx6WRO/HngjWRO/E7MiaVB0WkxIiK2iIirI+JPEfHRHq/v\nGBHzI+KS6se7unyeSaXVBlwP49kBeCl5Ju4U8srOR0O2AXs1kTTlGIc71orFpRSOJnfc9gfeT7YG\nb0Pjuk5jsTTlGIs71CMn/hvwFuBl1ZLTgDnkvB5zYmlwdHabRkRMI3OxlwE3AL+OiGNLKVe2lh5V\nStmlq+eYlOo24OoM3ANky9mewN3AXtWPh9frvZ5ImpKMwx1qXhNXuZYsPpwEbAwcQx7RAIzD0hRm\nLO5Qj5z4IOCjwDB5POPjNOZCGIulgdJlZ8QmwJ9KKX8updxHXuSwdYefb2poDeO5DHge8B4enA/x\neRqFCNvPpKnMONyF1i7cfcAXgY3I2Tz78uCsCMABlZKMxV1o5cSXAM8F3gs8ncyR98YBldIg67IY\nsT5wfePXN1Qfa3t9RFwWEUdHxGN7vVFE7BQRcyNi7vz587t41sHXGsZzFzkE7RnAH4HDgDNxGI+k\nRfQtDoOxGPiXAZVnk10QHwdeSR7J+B8cUClpEebE/dTKie8khwQ/E7gGOBz4BXl9J2BOLA2w8R5g\neRywYSnlqeRxrtm9FpVSDiqlzCqlzJoxY8aYPuC46jGMpwA/JYsOXyev7LwaeBsQDuORtOyWKg6D\nsbg5oPJmMu5uDtxDDqf8CfBYcEClpOVhTrw4o+TE9Yye/yWv7LyanJ9mTixNDJ3NjABupMrLKo+p\nPjailHJr45cHA1/p8Hkmlrr1rNp9oxT+DOwKnEAeyfgR8Jx6/bRp7r5JajMOr6hWLF5YCgcDHwPu\nIjsi9gKm1+tnzoR99jEWS2oyFq+IHjnxNeSNRfWMnqOBTev15sTShNFlZ8SvgSdExL9HxEOBbYFj\nmwsiYt3GL19Ddriq1Qb8T3IOxEZkS/C+wFwahQjbgCX1ZhxeXq02YICLgc2Ancnk91JgH6pChG3A\nktm8D8gAABG1SURBVEZnLF5erZz4XnIOxEbAeWRHxK9pFCLMiaUJpbPOiFLKgojYhbxhchrw/VLK\nFRGxNzC3lHIssGtEvAZYANwG7NjV80wYrWE8pwK7kHMh3gDsR+uQobtwkkZhHF5OrV2424FPAN8h\nr4UbAt5MdVUnOKBS0mIZi5dTKyc+mcyJrwHeRG7Orddcb04sTThRShnvZ1gms2bNKnPnzh3vx+i/\n1jVxNwC7k0cxngB8G3h5c71XE0kDJyIuLqXMGu/nGAuTOhZX18QV8l76PYBbyGs79wbWaK43FksD\nxTj8/9u7+2C7qvKO49+niRFiKToJOFTITVBoDTBtNTLYjlIFHKBKyhQpCMXaiC9UpNBqGSKMIBEY\n5UUqUwiVqiG8WLV4QSi1BIZKC4WR9yA0QAIBRhJQoE2BcFn9Y53AZnMg597ss/c9+3w/M2e4L/tc\n1so9+WXNOs9+Vgt0WRMfTb4VY0fymniv4vXmsDTp9JrFTTew1NKlMHNmLgVetYr1wNfJHYBHyQvf\nO+lsRNiMR5Kq16VB5Z3A7sBhwBzyrXFn09mIMIslqXqlNfHz5MYZv03ul3Yy+bjOvcAcllqinw0s\ntTGlUuDryO+8LQc+AnyDvAgGbMYjSf1QyuGnU+JEcv6+GTiffGrRSzv3lgFLUvVKWbyMfEvGPcB8\ncm+I2RuudU0stYaVEU0oNUZ7DDgE+ACwDvgRuSripY0Im/FIUrVKOZyAi8nvwJ3Jy8cmf5LOP5Q2\nqJSk6pWy+FFyT549yA3cLwcuo7AR4ZpYahU3I+rS5Xzk9cDp5PvffgCcQK6K2K/4PMvPJKk6pTJg\ngLvIm8EfIzdDuxFYDMzY8BwbVEpSdV5jTfx14LeAfwa+TM7mDxef55pYah1v06hDl/ORryWXny0H\n/ohcEvz24nNsxiNJ1SrfkkFe8J4NbAmcS66EmFJ8jlksSdXpsia+BjiSfEvGR8i3ZGxffI45LLWW\nlRH9VjofeTX5OKIPkm/JGAWuoLQRMWOGoStJVelyS8YSclXaWeRbMu4DPk1nI8LGaJJUvdKa+GHg\nQGBPXr4lY5TSRoRrYqnV3Izop8L5yM8Bp5LLz0bJ78ZtaFT5kpGRfE/y2rWGriRVYUMOd27JuA14\nH/mUjBHgJkq3ZIyMwJIlkJL9ISSpKqU18VfJPXouJ58cdzddbslwTSy1nrdp9EPpfOSrgKOA/wb+\nGDiDQnNKsPxMkqpWyuEngePJt2LMAP4B+ASFHXlzWJKqV8riHwN/BawA9ieviWcXrzeLpaFiZURV\nujTjuZ/cjHLfziVXkZvyzAHLgCWpH0oNKseA88i3ZJwLHEE+JWMBhX8ALQOWpOp0WROvIFc+fJh8\nO9zVwA/pbES4JpaGlpURVSg14/nflDiF3BX4DcBp5F3gaRuu93xkSapeKYtvAD4P/AzYHfg7YJfi\n9SMjsGiRWSxJVSnl8P8U1sTTgK+Rc9k1sSRwM2LTbWjGMzZGAi4FvkBuVHkoeSPiN4vXW34mSdUq\nlQE/CvwtcCGwLXAJuUlabLjeHJak6pXWxJeQ18SPAH9GXhNvU7zeLJaGnrdpTES5/GxsjNvJ77wd\nDGwF/JTcrf0VGxGWn0lSNbqUAT9HXuzuCHwPWAj8nHyC0UsbEeawJFWny5r4VuD9wMeAt5Kr1L5L\naSPCLJaElRHjVyo/W5sSx5O7sb+FfG/yAjynXpL6ppTDKSWuAI4hN0WbD5xO6chkc1iSqtVlTfwl\nXj6haDH56GTXxJJei5URvSqdU/8C8E3yO3DnA39JPi3jU3hOvST1Temc+p8D+5CbBU8F/gW4jNJG\nhA0qJak6pTXxeuBsYAfySUWfJ6+JD8c1saTXZ2VEL0o7v9eQj+q8G9gTOAvYqXi9zXgkqXqFc+p/\nBZxI3hR+E3AmeVP4DcXrbVApSdUqrYl/Qm7SvhzYi7wmnlu83jWxpNdhZcTGFN6Fu598JvKewP+R\njyT6V0obEdOnG7qSVKXCu3Bj69axmPwO3DfIJcD3kRfDL21ETJ8OF14IK1eaxZJUlcKaeMMtcR8C\nniNXpF1NaSPCNbGkjXAzoptSM55nxsY4lhywPwFOIVdF7E+nKZrlZ5JUrS4NKq8D3g18Gngn+cjO\n84CtwRyWpH4oZfHTY2N8kbwmXgacSl4Tz8c1saTx8zaNskL52YvAt1PiOOAXwMeBr9LlhAzLgCWp\nOqUy4AdT4gvAD4BZ5JMyDqBwQoZlwJJUvUIWj5HXxAvJa+JPAIvockKGa2JJ4+BmxAalc+p/Su4L\n8TPgvcAosGvxersBS1L1CufUP0PeAD6D/I/VScDfAJsXrzeLJalapTXx9eRb4W4F/gC4HHhP8Xpz\nWNIEeZvG0qUwc2buCLxqFSuBA4H3AY8DF5HPR37FRsSUKYauJFWlVAb84tgYF5D7QpwK/Cm5L8Tx\nlDYiLAOWpOqU1sQPAh8FdgfWAhcD/05pI8I1saRNMNyVEYXys2fIvSDOIB9DdCL5Hbjp5ee4+ytJ\n1SndknF9ShzN61SlgTksSVUrZPHT5Kq0s8hr4q8Af01pMxjMYkmbbDgrI0qd2b9FfgfuFPIO8L3A\nCRQ2ImzGI0nVKp1T/wC5D8TuwBq6VKWZw5JUvdKa+Hzymvg0Xq5K+xKFjQizWFKFhqcyonj/WwSk\nxDLgGOB2XucdOJvxSFJ1li6Fo46CJ54A4CngZOBs8tGcXd+Bs0GlJFWny5r4GvKa+A5yX4grKN2O\nAa6JJVVuODYjSmXA96XEF4EfASPAJeQ+EVF8jqVnklStQha/AJxPrkJ7gnxa0SJKpxWBWSxJVSqt\nie/tnFZ0OTAbuJRcJeyaWFId+nqbRkTsHRH3RsSKiDi2y/ffGBGXdr5/U0TMrnwQGzqzr1vHk8DR\nwE7ANeSF7z3kMrRXhO6MGYaupNZoPItLt2RcDfwucAT5rPqbgX+ksBFhGbCklmk8h+FVa+KjgJ2B\n68i3Kt9DlzfnXBNL6qO+bUZExBTgHGAf8nrz4IiYW7psAfDLlNI7gDPJt6hVp7P7+/zYGGcB7yCX\nAv8FsAI4ji6d2S+8ENauNXQltULjWbzhXbhVq7gL2LvzeBb4IXkR/O7i9SMjsGQJpAQrV5rFkgZe\n4zkMr1gTnwm8Hfgm8EnymvhYYLPi9a6JJdWgn5URuwIrUkoPpJSeJ98NMb90zXzgO52Pvw/sERFB\nVRYu5N5169iJXBExD7gNOA94a/G66dNz4LrwldQ+zWbxwoWwbh3HAb8D3EQ+tWg5sD+Fd+DMYUnt\nNSnWxPesW8dccm+I3cg90/4e2Lp4nVksqUb93Ix4G/Bw4fPVna91vSal9AK5l9mM8g+KiE9FxC0R\nccuaNWt6H8FDDzFCroi4Erga2OXlH5r/axmwpHZrNosfegiAOcCR5HfgjgamFa+xDFhSu02KNfFs\n8kkZV3UeO7/8Q/N/XRNLqtlAHO2ZUlqcUpqXUpq31VZb9f7EWbPYjBy4+1B4B27KFMuAJWmcJpTF\ns2YBcDj5zPpXrKwtA5akcdmUNfHm5DXx3sWvuyaW1KB+bkY8AmxX+Hzbzte6XhMRU4EtyY3Vq7Fo\nUS43K5o+3SPiJA2TZrP4tXLYMmBJw8M1sSR10c/NiJuBHSJiTkRMAw4CRkvXjJJPdAM4AFiWUkqV\njeCQQ3K52chILkGz/EzS8Gk2i81hSXJNLEldTO3XD04pvRARnyO3apgCXJBSujsiTgJuSSmNAt8C\nlkTECuBJcjhX65BDDFpJQ2tSZLE5LGmITYocBrNY0qTTt80IgJTSleTekcWvnVD4+Fngo/0cgyQN\nO7NYkpplDkvSqw1EA0tJkiRJktQebkZIkiRJkqRauRkhSZIkSZJq5WaEJEmSJEmqVVR5alAdImIN\nsGoCT50JrK14OE1r45zAeQ2SNs4JJj6vkZTSVlUPZjKaYBb7ehksbZxXG+cEzqvIHN64Nr5e2jgn\ncF6DpI1zgj6viQduM2KiIuKWlNK8psdRpTbOCZzXIGnjnKC982paW/9cndfgaOOcwHlpfNr459rG\nOYHzGiRtnBP0f17epiFJkiRJkmrlZoQkSZIkSarVMG1GLG56AH3QxjmB8xokbZwTtHdeTWvrn6vz\nGhxtnBM4L41PG/9c2zgncF6DpI1zgj7Pa2h6RkiSJEmSpMlhmCojJEmSJEnSJOBmhCRJkiRJqlWr\nNiMiYu+IuDciVkTEsV2+/8aIuLTz/ZsiYnb9oxy/HuZ1TEQsj4g7IuKaiBhpYpzjtbF5Fa77k4hI\nETHpj8vpZU4RcWDn93V3RFxU9xgnoofX4KyIuDYibu28DvdtYpzjEREXRMTjEXHXa3w/IuLszpzv\niIh31T3GQWUWD04WtzGHwSw2i2UOD04OQzuz2Bw2h3uSUmrFA5gC3A9sD0wDbgfmlq45Aji38/FB\nwKVNj7uieX0AmN75+LNtmVfnui2A64EbgXlNj7uC39UOwK3AWzqfb930uCua12Lgs52P5wIrmx53\nD/N6P/Au4K7X+P6+wFVAALsBNzU95kF4mMWDk8VtzOFx/K7M4knyMIsbe62Yw5Pk0cYsNofN4V4f\nbaqM2BVYkVJ6IKX0PHAJML90zXzgO52Pvw/sERFR4xgnYqPzSildm1Ja1/n0RmDbmsc4Eb38vgC+\nApwGPFvn4CaolzkdDpyTUvolQErp8ZrHOBG9zCsBv9H5eEvg0RrHNyEppeuBJ1/nkvnAd1N2I/Dm\niNimntENNLM4G4QsbmMOg1lsFssczgYhh6GdWWwOZ+bwRrRpM+JtwMOFz1d3vtb1mpTSC8BTwIxa\nRjdxvcyraAF552qy2+i8OiVA26WUflznwDZBL7+rHYEdI+KGiLgxIvaubXQT18u8vgwcGhGrgSuB\nI+sZWl+N9++eMrM4G4QsbmMOg1lsFssczgYhh6GdWWwOm8M9mVrFD9HkEBGHAvOA3Zsey6aKiF8D\nzgD+vOGhVG0quSztD8m79ddHxC4ppV81OqpNdzDw7ZTS6RHxXmBJROycUnqx6YFJdWtLFrc4h8Es\nllqtLTkMrc5ic1itqox4BNiu8Pm2na91vSYippJLZ56oZXQT18u8iIg9gYXAfiml52oa26bY2Ly2\nAHYGrouIleT7k0YnecOeXn5Xq4HRlNL6lNKDwH3kIJ7MepnXAuB7ACml/wQ2A2bWMrr+6envnl7F\nLB6cLG5jDoNZbBbLHB6cHIZ2ZrE5jDncizZtRtwM7BARcyJiGrkZz2jpmlHg452PDwCWpU5Xjkls\no/OKiN8DziOH7iDcbwUbmVdK6amU0syU0uyU0mzyfX/7pZRuaWa4PenlNXgZeQeYiJhJLlF7oM5B\nTkAv83oI2AMgIt5JDt41tY6yeqPAYZ0OwrsBT6WUHmt6UAPALB6cLG5jDoNZbBbLHB6cHIZ2ZrE5\njDnck/F0u5zsD3Knz/vIXU4Xdr52EvkvLOQXwz8BK4D/ArZveswVzevfgF8At3Ueo02PuYp5la69\njkneObjH31WQS+2WA3cCBzU95ormNRe4gdxV+DbgQ02PuYc5XQw8Bqwn784vAD4DfKbwuzqnM+c7\nB+H1N1keZvHgZHEbc7jH35VZPEkeZnFjrxVzeBI92pjF5rA53MsjOv8DSZIkSZKkWrTpNg1JkiRJ\nkjQA3IyQJEmSJEm1cjNCkiRJkiTVys0ISZIkSZJUKzcjJEmSJElSrdyMkCRJkiRJtXIzQpIkSZIk\n1crNCA2diHhPRNwREZtFxJsi4u6I2LnpcUnSMDGLJalZ5rCaFimlpscg1S4iTgY2AzYHVqeUTml4\nSJI0dMxiSWqWOawmuRmhoRQR04CbgWeB308pjTU8JEkaOmaxJDXLHFaTvE1Dw2oG8OvAFuTdYElS\n/cxiSWqWOazGWBmhoRQRo8AlwBxgm5TS5xoekiQNHbNYkpplDqtJU5segFS3iDgMWJ9SuigipgD/\nEREfTCkta3pskjQszGJJapY5rKZZGSFJkiRJkmplzwhJkiRJklQrNyMkSZIkSVKt3IyQJEmSJEm1\ncjNCkiRJkiTVys0ISZIkSZJUKzcjJEmSJElSrdyMkCRJkiRJtfp/P4b30O5WfegAAAAASUVORK5C\nYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x111b65d10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plotting code for each of the methods above\n",
    "import os\n",
    "\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "u_true = lambda x: (4.0 - numpy.exp(1.0)) * x - 1.0 + numpy.exp(x)\n",
    "x_fine = numpy.linspace(0.0, 1.0, 100)\n",
    "\n",
    "base_path = \"./src\"\n",
    "file_names = [\"jacobi_omp1.txt\", \"jacobi_omp2.txt\", \"jacobi_mpi.txt\"]\n",
    "titles = [\"OpenMP - Fine-Grained\", \"OpenMP - Coarse-Grained\", \"MPI\"]\n",
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 3)\n",
    "for (i, title) in enumerate(titles):\n",
    "    path = os.path.join(base_path, file_names[i])\n",
    "    data = numpy.loadtxt(path)\n",
    "    x = data[:, 0]\n",
    "    U = data[:, 1]\n",
    "    \n",
    "    axes = fig.add_subplot(1, 3, i + 1)\n",
    "    axes.plot(x, U, 'ro', label='computed')\n",
    "    axes.plot(x_fine, u_true(x_fine), 'k', label=\"exact\")\n",
    "    axes.set_title(title)\n",
    "    axes.set_xlabel('x')\n",
    "    axes.set_ylabel('u(x)')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Software Packages\n",
    "\n",
    "### Linear Algebra\n",
    " - BLAS - Basic Lineag Algebra Subroutines\n",
    "   - Level 1: Scalar and vector operations\n",
    "   - Level 2: Matrix-vector operations\n",
    "   - Level 3: Matrix-matrix operations\n",
    " - LAPACK - Linear Algebra Package\n",
    " - ScaLAPACK - Parallel LAPACK\n",
    "\n",
    "### Solving ODEs\n",
    "ODEPACK\n",
    "\n",
    "### Solving PDEs"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}