\chapter{QUADRATURE WITH PARFOR}

\label{chapter_time_work_rate}


  \begin{figure}[h]
  \begin{centering}
    \includegraphics[width=0.49\textwidth]{chap_quad_parfor/quad_picture_sequential.png}
    \includegraphics[width=0.49\textwidth]{chap_quad_parfor/quad_picture_parallel.png}
    \caption{Sequential and Parallel Integral Estimates by Quadrature}
  \end{centering}
  \end{figure}

  
\section{Introduction}

The integral of a function $f(x)$ over the interval $[a,b]$ can be interpreted
as the area between the $x$-axis and the curve $y=f(x)$.  If we can't see a
way to determine the exact value of the integral through calculus, we can
think of approximating the integral by considering the area.  This approach
is called {\it{quadrature}}.  It makes an estimate of the integral by evaluating
the function at many points.  

Quadrature is an ideal application for parallel programming, because it
involves a large number of simple, independent tasks whose properties can be
defined in advance.  Notice however, that the quadrature problem raises
a new issue that did not arise before, namely {\it{communication}}.
In the combination lock problem, there was never any interaction between 
the workers, and the client simply assigned tasks to the workers and
never asked for a result; instead, a worker that found a combination simply
printed it directly to the user.  For the quadrature problem, the workers
still won't need to communicate with each other, but once each worker has
computed its partial estimate, it must communicate this information to the
client, so that it can be gathered into a single final estimate.  

This is a logical feature of the computation, and we shall want to know that
MATLAB is able to handle it properly.

\section{The Sequential Quadrature Problem}

We suppose that we are given a function $f(x)$, and an interval $[a,b]$. 
The integral is symbolized by $I(f,[a,b]) = \int_a^b f(x) dx$.  We suppose
that it is not possible to determine the exact value of $I(f,[a,b])$,
and so instead we plan to approximate the value by a quadrature formula,
which we will write as $Q(f,[a,b])$.

One of the simplest quadrature formulas is the {\it{midpoint rule}}:
\begin{displaymath}
  Q_M(f,[a,b]) = (b-a) \cdot f(\frac{a+b}{2})
\end{displaymath}
While the accuracy of this rule is low, improved results can be obtained
by using the {\it{composite midpoint rule}}, which involves dividing the
interval $[a,b]$ into $n$ equal subintervals, applying the midpoint
rule to each, and then summing the result:
\begin{align*}
  Q_{CM}(f,[a,b],n) &= \sum_{i=1}^n(b_i-a_i) \cdot f(\frac{a_i+b_i}{2})\\
                    &= \frac{b-a}{n} \cdot \sum_{i=0}^n f(\frac{(n-i)a+ib}{n})
\end{align*}

A simple MATLAB procedure for the quadrature problem might be:
\begin{lstlisting}
function q = quad_sequential ( a, b, n, f )
%% quad_sequential.m
  q = 0.0;
  for i = 1 : n
    ai = ( ( n - i + 1 ) * a + ( i - 1 ) * b ) / n;
    bi = ( ( n - i     ) * a +   i       * b ) / n;
    xi = ( ai + bi ) / 2.0;
    q = q + f ( xi );
  end 
  q = q * ( b - a ) / n; 

  return
end
\end{lstlisting}

Our computation is written as a MATLAB M-file; we invoke it by 
setting the values of the input.  Here, our function to integrate will
be $f(x) = (4x+1) \cdot (11-4x)/8$, and for convenience, we will use
MATLAB's {\it{anonymous function}} facility to define the function
$f$ on a single command line:

\begin{lstlisting}
%% quad_sequential_call.m
  a = 0.0;
  b = 2.0;
  n = 21;
  f = @(x) (4*x+1)*(11-4*x)/8;
  q = quad_sequential ( a, b, n, f );
  fprintf ( 1, '  The integral estimate is %f\n', q );
\end{lstlisting}


\section{Quadrature Using PARFOR}

We can make a parallel version of the quadrature computation simply by
replacing {\bf{for}} with {\bf{parfor}}.  
\begin{lstlisting}
function q = quad_parfor ( a, b, n, f )
%% quad_parfor.m
  q = 0.0;
  parfor i = 1 : n
    ai = ( ( n - i + 1 ) * a + ( i - 1 ) * b ) / n;
    bi = ( ( n - i     ) * a +   i       * b ) / n;
    xi = ( ai + bi ) / 2.0;
    q = q + f ( xi );
  end 
  q = q * ( b - a ) / n; 

  return
end
\end{lstlisting}

Before we go on, though, we need to try to make a mental model of how the
variables in the parallel program are accessed by the workers while executing
the loop.  It's easy to imagine that the variables {\bf{a}}, {\bf{b}}, and 
{\bf{n}}, as well as the expression for {\bf{f}} are ``owned'' by the client 
copy of MATLAB, and that the workers either get a copy of these
values, or can refer to them at any time (and this is correct).  

However, we can't really think about the variable {\bf{q}} in the same way.
The client initializes {\bf{q}} to zero, but doesn't execute the statements 
inside the loop which update its value.  So it's the workers
that have to update the value of {\bf{q}} with the new information.  If we
suppose there's only one copy of {\bf{q}} available, and we carefully insist
that only one worker at a time can update it, then we will surely have a
considerable amount of undesirable delay while workers wait for their turn to
apply their update.  On the other hand, if we simply allow any
worker at any time to read or write the single copy of {\bf{q}}, then we
are liable to get erroneous results.  Two workers could try to update at
the same time, and only one of the updates would be correctly stored.

Luckily, MATLAB is automatically able to recognize what is going on with
the variable {\bf{q}}, and essentially, as the loop starts, it creates 
for each worker a zeroed-out private variable called {\bf{q}}, which that
worker alone is free to modify. On the completion of the loop, MATLAB quietly 
gathers the multiple copies of {\bf{q}}, and adds them
to the value of {\bf{q}} that was stored in the client before the parallel
loop began. The variable {\bf{q}} is an example of what's called a 
{\it{reduction variable}}.

Notice that the variable {\bf{x}} has a simpler behavior.  It is simply
a temporary variable, which is only used inside the loop.  We might imagine
that each worker makes a separate private copy of {\bf{x}} to be used
during its iterations.  The variable {\bf{x}} is just a temporary convenience,
and once we exit the loop, its value is of no importance.  Thus, even though
the sequential code uses a single copy of {\bf{x}}, but the parallel code
uses multiple versions of this variable, we are right in assuming that this
is not going to cause any problems.

To invoke the code, once again we have to use the {\bf{matlabpool}} command
to request the help of the workers:

\begin{lstlisting}
%% quad_parfor_call.m
  a = 0.0;
  b = 2.0;
  n = 21;
  workers = 4;
  f = @(x) (4*x+1)*(11-4*x)/8;
  matlabpool ( 'open', 'local', workers );
  q = quad_parfor ( a, b, n, f );
  matlabpool close
  fprintf ( 1, '  The integral estimate is %f\n', q );
\end{lstlisting}

\section{Comparing Sequential and Parallel Rates}

When we try the parallel program out on our sample problem, for some
small problem sizes, the results don't look very good!  We run the
sequential code (``0 workers'') versus the parallel code with 1, 2 or 4
workers, and for $n$ ranging from 1 to $2^{10}$.  The parallel
timings are consistently worse than the sequential code and even seem
to get worse as we add more workers!

However, it's important to be cautious here.  Note that even the sequential
timings are not showing a pattern of doubling as the problem size doubles.
This suggests that the computational work in the problem is insignificant
compared to the other things that MATLAB is doing for us.  We have to remember
that the benefits of parallelism come for problems with a ``significant'' amount
of work.  Certainly our quadrature function $f$ is trivial to evaluate, and 
perhaps a few hundred evaluation points of an easy function don't constitute
enough work.  

\begin{center}
\begin{tabular}{|r|r|r|r|r|}
\hline
         N/W&      0&         1&         2&         4\\
\hline
         1&    0.2923&    0.4164&    0.2306&    0.2436\\
         2&    0.0073&    0.0294&    0.0270&    0.2131\\
         4&    0.0008&    0.0234&    0.0223&    0.2150\\
         8&    0.0008&    0.0205&    0.0235&    0.0293\\
        16&    0.0008&    0.0263&    0.0263&    0.0391\\
        32&    0.0008&    0.0306&    0.0395&    0.0435\\
        64&    0.0008&    0.0207&    0.0326&    0.0489\\
       128&    0.0008&    0.0201&    0.0303&    0.0410\\
       256&    0.0009&    0.0201&    0.0332&    0.0362\\
       512&    0.0010&    0.0213&    0.0258&    0.0320\\
      1024&    0.0012&    0.0204&    0.0250&    0.0342\\
\hline
\end{tabular}
\end{center}

We still assume that there are some versions of this problem for which
parallelism provides a benefit.  How do we look for them?  We could look at 
cases where the function is more expensive to evaluate, but instead, let's 
stick with our current function, gradually increase the value of $n$, and see if
there's at least some range over which the parallel code runs faster.

  \begin{figure}[h]
  \begin{centering}
    \includegraphics[width=0.95\textwidth]{chap_quad_parfor/quad_timings.png}
    \caption{Rates for 1, 2, 4 workers, relative to sequential code}
  \end{centering}
  \end{figure}

Our comparisons will use n = 1, 2, 4, ..., $2^{23}$ points, and we will try
the sequential code ( ``0 workers''), as well as parallel runs with 1, 2, and 4
workers.  We suspect that the 1 worker code will {\it{never}} run faster than
the sequential code, but we hope that there will be a range where the 2 and
4 worker codes are about twice and four times as fast as the sequential code.

To make this trend easy to spot, we'll plot the computational rates relative
to the sequential code.  That means that, for every problem size, the rate of
the sequential code will be plotted as 1, showing up as a horizontal line.  
If a parallel code is faster for that problem size, it will appear above the 
line.

From the plot, we get some reassurance and understanding.  The sequential code
is represented by the blue line.  The 1 worker code is red, and we see that 
it is never better than the sequential code, and at the highest values of $n$
seems to be leveling off in rate.  The 2 and 4 worker codes first beat the sequential
code at $n=2^16$, and thereafter, their rates seem to be rising towards their
theoretical limits of 2 and 4.

We can see that if we had only studied values of $n$ up to $2^15$, we would never have
seen a benefit from parallelism.  The fact that we have to wait so long to see
a benefit is also because the integrand is so simple to evaluate.  This 
suggests that the small, simple problems we like to use for benchmarking will
often lead us to incorrect conclusions.  To understand whether parallelism will
benefit your problem, it's really important to examine a realistic range of problem
sizes, parameters, and difficulties.
