%% LATEX TWICE TO GET THE REFERENCES RIGHT
%%\documentstyle[12pt]{article}
\documentstyle[12pt,epsf,psfig]{article}

\newcommand{\postscript}[2]
 {\setlength{\epsfxsize}{#2\hsize}
  \centerline{\epsfbox{#1}}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newtheorem{theorem}{Theorem}
\newtheorem{itlemma}{Lemma}[section] %number by section (set in \em by default)
%%%same numbering counter as itlemmas:
\newtheorem{itproposition}[itlemma]{Proposition}
\newtheorem{itcorollary}[itlemma]{Corollary}
\newtheorem{itremark}[itlemma]{Remark}
\newtheorem{itnote}[itlemma]{Note}
\newtheorem{itobservation}[itlemma]{Observation}
\newtheorem{itdefinition}[itlemma]{Definition}
\newtheorem{itexercise}[itlemma]{Exercise}
\newtheorem{itexample}[itlemma]{Example}

\newenvironment{lemma}{\begin{itlemma}\rm}{\end{itlemma}} %no-italics
\newenvironment{remark}{\begin{itremark}\rm}{\end{itremark}} %no-italics
\newenvironment{observation}{\begin{itobservation}\rm}{\end{itobservation}} %no-italics
\newenvironment{corollary}{\begin{itcorollary}\rm}{\end{itcorollary}}
\newenvironment{proposition}{\begin{itproposition}\rm}{\end{itproposition}}
\newenvironment{definition}{\begin{itdefinition}\rm}{\end{itdefinition}}
\newenvironment{exercise}{\begin{itexercise}\rm}{\end{itexercise}}
\newenvironment{example}{\begin{itexample}\rm}{\end{itexample}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This are for choices.  They take 4 or 6 args
%% (carriage returns allowed in between}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\twoif}[4]{
\left\{ \begin{array}{ll}#1&#2\\#3&#4\end{array}\right.
}
\newcommand{\threeif}[6]{
\left\{ \begin{array}{ll}#1&#2\\#3&#4\\#5&#6\end{array}\right.
}

\def\partmlab#1#2#3#4#5#6#7#8{\bordermatrix{&#1&#2\cr#3&#5&#6\cr

\newcommand{\bigframebox}[1]{\framebox{\parbox{\smallerwidth}{#1}}}

\def \beq {\begin{eqnarray}}
\def \eeq {\end{eqnarray}}
\def \beqn {\begin{eqnarray*}}
\def \eeqn {\end{eqnarray*}}

\newcommand{\rf}[1]{[\ref{#1}]}            %bibliographical references

\newcommand{\text}[1]{\hbox{\rm \ #1\ \/}}
\newcommand{\be}[1]{\begin{equation}\label{#1}}
\newcommand{\ee}{\end{equation}}
\newcommand{\bl}[1]{\begin{lemma}\label{#1}}
\newcommand{\ble}[1]{\begin{lemmaex}\label{#1}}
\newcommand{\br}[1]{\begin{remark}\label{#1}}
\newcommand{\bn}[1]{\begin{note}\label{#1}}
\newcommand{\bo}[1]{\begin{observation}\label{#1}}
\newcommand{\bt}[1]{\begin{theorem}\label{#1}}
\newcommand{\bd}[1]{\begin{definition}\label{#1}}
\newcommand{\bp}[1]{\begin{proposition}\label{#1}}
\newcommand{\bc}[1]{\begin{corollary}\label{#1}}
\newcommand{\bfact}[1]{\begin{fact}\label{#1}}
\newcommand{\ber}[1]{\begin{exercise}\label{#1}}
\newcommand{\bex}[1]{\begin{example}\label{#1}}
\newcommand{\bem}[1]{\begin{example}\label{#1}}    %Yes, 2 different ones...
\newcommand{\ec}{\mybox\end{corollary}}
\newcommand{\efact}{\mybox\end{fact}}
\newcommand{\eer}{\mybox\end{exercise}}
\newcommand{\eex}{\mybox\end{example}}
\newcommand{\eem}{\mybox\end{example}}
\newcommand{\el}{\mybox\end{lemma}}
\newcommand{\ele}{\mybox\end{lemmaex}}
\newcommand{\er}{\mybox\end{remark}}
\newcommand{\en}{\mybox\end{note}}
\newcommand{\eo}{\mybox\end{observation}}
\newcommand{\et}{\qed\end{theorem}}
\newcommand{\ed}{\mybox\end{definition}}
\newcommand{\ep}{\mybox\end{proposition}}
\newcommand{\epr}{\end{proof}}
\newcommand{\bpr}{\begin{proof}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  These are just for statements, no "qed" box is put in.
%%  (Also for when I have to finish something after a displayed equation
%%        or itemized.)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\ecs}{\end{corollary}}
\newcommand{\eers}{\end{exercise}}
\newcommand{\eexs}{\end{example}}
\newcommand{\eems}{\end{example}}
\newcommand{\els}{\end{lemma}}
\newcommand{\eles}{\end{lemmaex}}
\newcommand{\ers}{\end{remark}}
\newcommand{\ens}{\end{note}}
\newcommand{\eos}{\end{observation}}
\newcommand{\ets}{\end{theorem}}
\newcommand{\eds}{\end{definition}}
\newcommand{\eps}{\end{proposition}}

\newcommand{\qed}{\hfill \halmos} %put \qed at right margin
\newcommand{\mybox}{\hfill $\Box$} %put \qed at right margin (white square)
 
\def \R {\bb R} % reals
\def \C {\bb C} % complexes
\def \Z {\bb Z} % integers
\def \X {\bb X} % for this paper only!
\def \inter {{\rm int \,}}
\def \rank {{\rm rank \,}}
%%%\def \inf {{\rm inf \,}}  %better in Latex -- subscripts OK!
%%%\def \min {{\rm min \,}}
\def \im {{\rm im \,}}
\def \ker {{\rm ker \,}}
\def \norm #1{\Vert #1 \Vert}
\def \ip #1\by #2{\langle #1 , #2 \rangle}
\def \Sone {{\bf S}^1} % circle
\def \wigglearrow {\sim \mkern -6mu >}
        
\newcommand{\comment}[1]{}
\newcommand{\halmos}{\rule{1ex}{1.4ex}}
%%%alternative to \qed (see below)
\newenvironment{proof}{\noindent {\em Proof}.\ }{\hspace*{\fill}$\halmos$\medskip}
 
\def\bi{\begin{itemize}}
\def\ei{\end{itemize}}
\def\ben{\begin{enumerate}}
\def\een{\end{enumerate}}
\def\st{\, | \,}

\newcommand{\myboxineq}{\eqno\Box} %put \qed at margin (white square) in eqn
\newcommand{\qedineq}{\eqno\halmos} %put \qed at margin in eqn

%%%% FROM HERE, macros special for paper
\def\mm{m}
\def\pp{p}
\def\proj{\pi _1} % projection on first group of coordinates
\def\x{u}  % vector argument
\def\SEC{(SEC)}
\def\SECo{(SEC$^0$)}
\def\INV{(INV)}
\def\Fp{{\cal F}_\pp}
\def\Fpm{{\cal F}_\pp^\mm}  %there is a 3,3 and a 2,1 in text

\def\Rm{\R^m}
\def\Rn{\R^n}
\newcommand{\ball}[2]{{\bb B}(#1,#2)}
\newcommand{\cball}[2]{\overline {\bb B} (#1,#2)}
\def\Fnm{{\cal F}_n^\mm}
\def\clos{{\rm clos\,}}
\def\inter{{\rm int\,}}

\def\CCF{compatible class of functions}
\begin{document}


{\sl Example with MPI}


{\bf 
Computing $\pi$ via Simpson's Rule and Message Passing Interface (MPI)
}


\begin{center}
Saleh Elmohamed\\
Northeast Parallel Architectures Center,\\
Syracuse University\\
Syracuse, NY 13244-4100\\
Saleh@npac.syr.edu\\
\end{center}

\tableofcontents

\section{Computation of $\pi$ - Equations}
\setcounter{page}{1}

This example computes $\pi$ using Simpson's rule and MPI
on the SP2 and AlphaCluster. 

\be[pi-eq1]
\pi = {\displaystyle \int_{0}^1} {4 \over {1+x^2}} dx
\ee

and this is equal to:

\be{pi-eq2}
4 \ln(2) + \text{constant}
\ee

This integration can be approximated numerically by:

\be{pi-eq3}
\pi \approx {\displaystyle \sum_{i=1}^n} {{4n} \over {n^2 + {(i + 0.5)}^2}} 
\ee
Of course, the larger the value of n, the better the approximation.

When using Simpson's rule and integrating on a multiprocessor, the
logical approach would be to partition the interval to subintervals
and make each processor run on one of these subintervals as follows:

\begin{verbatim}

  x1 = ((double) my_rank)/((double) Procs);
  x2 = ((double) (my_rank + 1))/((double) Procs);

\end{verbatim}

In turn, each processor applies the Simpson's rule to its own interval,
again, partitioning it into subintervals of equal number to the
above. On each of these subinterval $(p_j,p_{j+1})$ where 
$l = p_{j+1} - p_j$, we have:

\be{pi-eq4}
{\displaystyle \int_{p_j}^{p_{j+1}}} f(x) dx \approx 
{l \over 6}(f(p_j) + f({{p_j + p_{j+1}} \over 2}) + f(p_{j+1}))
\ee

Double-precision arithmetic is used throughout the C code, and if
needed with {\it casting}. Also, each point in the interval $1\ldots
N$ is  computed as follows:

\begin{verbatim}
  for(i = 1; i < N ; i++)
    {
      y = x1 + (x2 - x1)*((double) i)/((double) N);
                               /* call Simpson's rule  */
      sum = (double) simpson(i, y, h, sum);
 
      ... etc. ...
    }
\end{verbatim}

Among the processor, one node (processor) takes charge (usually the
node running the program) receiving and adding up the results. The
final result is according to the above \ref{pi-eq1} and \ref{pi-eq4}
equations. One other important point is that my knowledge of calculus
says that equation \ref{pi-eq1} is equal to $4$ multiplied by the
inverse trigonometric function of $\tan$ which is:
\be{pi-eq5}
4 \arctan(1) 
\ee

In C, we use the double-precision function $atan()$ to take care of
that. Other function we experimented with:
\be{pi-eq6}
4 \log(value)
\ee
where $value$ is the constant 2.19328059, but the results are not as
accurate as those obtained from $\arctan$. Take a look at the 3rd
section.

\section{Computation of $\pi$ - The Program}

\begin{verbatim}

#include  <stdio.h>
#include  <math.h>
#include  <mpi.h>              /* MPI lib routines */

#define approx_val  2.19328059

#define N 32                   /* Number of intervals in each processor */
double integrate_f();          /* Integral function */

double simpson();

double integrate_f(double x)
{
  return 4.0/(1.0 + x * x);     /* compute and return value */
}

double simpson(int i, double y, double h, double sum)
{
                               /* store result in sum */
  sum += integrate_f(y);
  sum += 2.0 * integrate_f(y - h);
  if(i == (N - 1))
    sum += 2.0 * integrate_f(y + h);
  return sum;
} /* simpson */

void main(int argc, char *argv[])
{
  int Procs;                   /* Number of processors */
  int my_rank;                 /* Processor number */
  double total; 
  double exact_val_of_Pi, pi, y, processor_output_share[8], x1, x2, h, sum;
  int i;
  MPI_Status status;
                               /* Let the system do what it needs to start up MPI */
  MPI_Init(&argc, &argv);
                               /* Get my process rank */
  MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
                               /* Find out how many processes are being used. */
  MPI_Comm_size(MPI_COMM_WORLD, &Procs);
                               /* Each processor computes its interval */
  x1 = ((double) my_rank)/((double) Procs);
  x2 = ((double) (my_rank + 1))/((double) Procs);
                               /* l is the same for all processes. */
  h = 1.0/((double) (2 * N * Procs));
  sum = 0.0;
  for(i = 1; i < N ; i++)
    {
      y = x1 + (x2 - x1)*((double) i)/((double) N);
                               /*--------------- call Simpson's rule  */
      sum = (double) simpson(i, y, h, sum);
     }
                               /* Include the endpoints of the intervals */
  sum += (integrate_f(x1) + integrate_f(x2))/2.0;
  total = sum;
                               /* Add up the integrals calculated by each process. */
  if(my_rank == 0)
    {
      processor_output_share[0] = total;
                               /* source = i, tag = 0 */
      for(i = 1; i < Procs; i++)
        MPI_Recv(&(processor_output_share[i]), 1, MPI_DOUBLE, i, 0, MPI_COMM_WORLD, &status);
    }
  else
    {
                               /* dest = 0, tag = 0 */
      MPI_Send(&total, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
    }
                               /* Add up the value of Pi and print the result.  */
  if(my_rank == 0)
    {
      pi = 0.0;
      for(i = 0; i < Procs; i++)
        pi += processor_output_share[i];
      pi *= 2.0 * h/ 3.0;
      printf("-------------------------------------------------\n");
      printf("The computed Pi of the integral for %d grid points is  %30.18e \n",
             (N * Procs), pi);

/* This is directly derived from the integration of the formula. See
   section 1 above. */

      exact_val_of_Pi = 4.0 * atan(1.0);

/*---------
     exact_val_of_Pi = 4.0 * log(approx_val); 
----------*/
/*=====================
   Since we're dealing with double values, fabs() is the function we need
   rather than abs().
=======================*/

      printf("The discrepancy between exact and computed values of Pi : %30.18e\n
             fabs(pi - exact_val_of_Pi));
    }
  MPI_Finalize();
}
\end{verbatim}


\section{Computation of $\pi$ - The Output}

\subsection{On the SP2 - NPAC}

For $N = 16$, $Proc = 4$, and the use of $\arctan()$, we get:

\begin{verbatim}
saleh:merlin2[44] gcc simpson.c -o simpson -L/usr/local/lib/rs6000/ch_p4 -lmpi -lm
saleh:merlin2[45] mpirun -mr_machinefile machinelist -mr_np 4 simpson

The computed Pi of the integral for 64 grid points is
3.1415926535897833e+00 
The discrepancy between exact and computed value of Pi :
9.7699626167013776e-15
\end{verbatim}

For $N = 32$, $Proc = 4$, and the use of $\arctan()$, we get:

\begin{verbatim}
The computed Pi of the integral for 128 grid points is
3.1415926535897936e+00 
The discrepancy between exact and computed value of Pi :
4.4408920985006262e-16
\end{verbatim}

For $N = 32$, $Proc = 4$, and the use of $\log()$, we get:

\begin{verbatim}
The computed Pi of the integral for 128 grid points is
3.1415926535897936e+00 
The discrepancy between exact and computed value of Pi :
9.8348027721684161e-07
\end{verbatim}

For $N = 32$, $Proc = 5$, and the use of $\arctan()$, we get:

\begin{verbatim}
The computed Pi of the integral for 160 grid points is
3.1415926535897931e+00 
The discrepancy between exact and computed value of Pi :
0.0000000000000000e+00
\end{verbatim}

For $N = 32$, $Proc = 6$, and the use of $\arctan()$, we get:

\begin{verbatim}
The computed Pi of the integral for 192 grid points is
3.1415926535897931e+00 
The discrepancy between exact and computed value of Pi :
0.0000000000000000e+00
\end{verbatim}

For $N = 32$, $Proc = 7$, and the use of $\arctan()$, we get:

\begin{verbatim}
The computed Pi of the integral for 224 grid points is
3.1415926535897931e+00 
The discrepancy between exact and computed value of Pi :
0.0000000000000000e+00
\end{verbatim}

For $N = 32$, $Proc = 8$, and the use of $\arctan()$, we get:

\begin{verbatim}
The computed Pi of the integral for 256 grid points is
3.1415926535897936e+00 
The discrepancy between exact and computed value of Pi :
4.4408920985006262e-16
\end{verbatim}

For $N = 64$, $Proc = 4$, and the use of $\arctan()$, we get:

\begin{verbatim}
The computed Pi of the integral for 256 grid points is
3.1415926535897922e+00 
The discrepancy between exact and computed value of Pi :
8.8817841970012523e-16
\end{verbatim}

\subsection{On the Alphafarm - NPAC}

For $N = 32$, $Proc = 7$, and the use of $\arctan()$, we get:

\begin{verbatim}
The computed Pi of the integral for 224 grid points is
3.1415926535897931e+00 
The discrepancy between exact and computed value of Pi :
0.0000000000000000e+00
\end{verbatim}

For $N = 64$, $Proc = 7$, and the use of $\arctan()$, we get:

\begin{verbatim}
The computed Pi of the integral for 448 grid points is
3.1415926535897931e+00 
The discrepancy between exact and computed value of Pi :
0.0000000000000000e+00
\end{verbatim}

For $N = 32$, $Proc = 8$, and the use of $\arctan()$, we get:

\begin{verbatim}
The computed Pi of the integral for 256 grid points is
3.1415926535897936e+00
The discrepancy between exact and computed value of Pi :
4.4408920985006262e-16
\end{verbatim}

For $N = 64$, $Proc = 8$, and the use of $\arctan()$, we get:

\begin{verbatim}
The computed Pi of the integral for 512 grid points is 
3.1415926535897931e+00
The discrepancy between exact and computed value of Pi :
0.0000000000000000e+00
\end{verbatim}

As we could see that the error is gone up when going to 8 processors
with the same number of interval points. Still, doubling the interval
points and going to only 4 processors didn't improve the error.

\end{document}
