\documentclass[12pt]{article}
\usepackage[dvips]{epsfig}

\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}}

%\tolerance=1600 %allow some tolerance in extending after line.
\parskip=3pt plus 1pt minus 1pt
%\overfullrule=0pt %no dark lines if overfull

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% These 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.
}
\newcommand{\fourif}[8]{
\left\{ \begin{array}{ll}#1&#2\\#3&#4\\#5&#6\\#7&#8\end{array}\right.
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\partmlab#1#2#3#4#5#6#7#8{\bordermatrix{&#1&#2\cr#3&#5&#6\cr
                                        \noalign{\medskip}#4&#7&#8\cr}}
% For 2by2 partitioned matrices; arguments are: two labels of top, side,
% 4 entries (row by row).  If entries are themselves matrix or pmatrix,
% better for readability of source to define first: \def\matrix1{...} etc.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\newcommand{\bigframebox}[1]{\framebox{\parbox{\smallerwidth}{#1}}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Warning: use \noindent before!
% AND IN TEXT, after begin doc, do:
% \newlength{\smallerwidth}  %for bigboxes
% \setlength{\smallerwidth}{\textwidth}
% \addtolength{\smallerwidth}{-8pt}        %MUST BE DONE AFTER begin doc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\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)

%\newcommand{\pmat}[2]{\left( \begin{array}{#1} #2 \end{array} \right) }
%use as e.g. \pmat{cc}{ a&b\\ a &c }  Note 'cc' for 2 Centered cols
%can also use \pmatrix{ a&b\cr a &c \cr} for plainTeX standard construct
%\setcounter{page}{27} should set the CURRENT page to 27

\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 >}

\def \today {\ifcase \month \or January\or February\or  
  March\or April\or May\or June\or July\or August\or 
  September\or October\or November\or December\fi\   
  \number \day, \number \year}                       
\newcommand{\comment}[1]{}
\newcommand{\halmos}{\rule{1ex}{1.4ex}}
%%%alternative to \qed (see below)
\newenvironment{proof}{\noindent {\em Proof}.\ }{\hspace*{\fill}$\halmos$\medskip}

%\def\qed{\kern 1em \vrule height .5em width .5em depth 0em}

\def\vbar{\mathchoice{\vrule height6.3ptdepth-.5ptwidth.8pt\kern-.8pt}
   {\vrule height6.3ptdepth-.5ptwidth.8pt\kern-.8pt}
   {\vrule height4.1ptdepth-.35ptwidth.6pt\kern-.6pt}
   {\vrule height3.1ptdepth-.25ptwidth.5pt\kern-.5pt}}
\def\fudge{\mathchoice{}{}{\mkern.5mu}{\mkern.8mu}}
\def\bbc#1#2{{\rm \mkern#2mu\vbar\mkern-#2mu#1}}
\def\bbb#1{{\rm I\mkern-3.5mu #1}}
\def\bba#1#2{{\rm #1\mkern-#2mu\fudge #1}}
\def\bb#1{{\count4=`#1 \advance\count4by-64 \ifcase\count4\or\bba A{11.5}\or
   \bbb B\or\bbc C{5}\or\bbb D\or\bbb E\or\bbb F \or\bbc G{5}\or\bbb H\or
   \bbb I\or\bbc J{3}\or\bbb K\or\bbb L \or\bbb M\or\bbb N\or\bbc O{5} \or
   \bbb P\or\bbc Q{5}\or\bbb R\or\bbc S{4.2}\or\bba T{10.5}\or\bbc U{5}\or
%%%   \bbb P\or\bbc Q{5}\or\bbb R\or\bba S{8}\or\bba T{10.5}\or\bbc U{5}\or
   \bba V{12}\or\bba W{16.5}\or\bba X{11}\or\bba Y{11.7}\or\bba Z{7.5}\fi}}
%
%
\renewcommand{\thefootnote}{\fnsymbol{footnote}}
%%% read as: theFOOTNOTE (i.e. the command  executed when doing a footnote) is
%%% now to be set as a FOOTNOTESYMBOL (type of number) value of counter
%%% FOOTNOTE (Instead of Roman).  See LaTeX manual pp.91-92 as well as 175.
\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

\topmargin -0.5in %default is 0.0, i.e. 1 inch from top.  Now it is 0.5
\textheight 9in
\oddsidemargin 0.25in %more reasonable margins!  (This is 1.25; can use <0 too)
\evensidemargin 0.25in %(same as above)
\textwidth 6.25in %(could use also -.5,-.5,7.5 for instance)
%\tolerance=1600 %allow some tolerance in extending after line.
\parskip=3pt plus 1pt minus 1pt
%\overfullrule=0pt %no dark lines if overfull

%%%% 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}

\bibliographystyle{alpha}
\newpage
\begin{flushright}
{\sl NPAC Discussion/NHSE Notes}  % Remove for submission to journal
\end{flushright}
\vspace{1cm}

\begin{center} 

{\LARGE\bf 
Notes on Parallel and Distributed 
\vspace{.5cm}

Computing}
\vspace{5cm}

{\large
Saleh Elmohamed\\
Geoffrey Fox\\}
\vspace{.5cm}
Northeast Parallel Architectures Center\\
Syracuse University, Syracuse, NY 13244-4100\\
\vspace{.25cm}
\{saleh, gcf\}@npac.syr.edu\\
\end{center}

%%\vspace{1cm}
%%\centerline{March 11, 1997}

%%\vspace{1.5cm}
%%\centerline{ABSTRACT}
%%A set of notes on programming models as well as various theoretical
%%and practical aspects of parallel and distributed computing.

\newpage
\pagestyle{empty}

\tableofcontents
\newpage
%% \listoffigures
%% \listoftables

\newpage
\pagestyle{plain}
\pagestyle{myheadings}
%% \markright{\sl SCCS report }

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{On Parallelism}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\setcounter{page}{1}

\subsection{Principles}

Principles of parallel computing are as follows:
\begin{itemize}
\item Parallelism and Amdahl's Law. 
\item Granularity. 
\item Locality. 
\item Load balance. 
\item Coordination and synchronization. 
\item Performance modeling. 
\end{itemize}

\subsection{Motivations}

Here are few of the main reasons for looking into parallelism: 
\begin{itemize}
\item It provides an interesting alternative computing paradigm. 
\item It applies at all levels of system design. 
\item It is an interesting perspective from which to view
  architecture.
\item It is becoming increasingly central in information processing.
\end{itemize}

Parallelism is exploited at several levels: 
\begin{enumerate}
\item Instruction-level parallelism.
\item Multiprocessor servers.
\item Large-scale multiprocessors (MPPs).
\end{enumerate}

Our focus in these notes is on the multiprocessor level of
parallelism.

\subsection{Limits on Performance}

It is well-known that the performance of any parallel program is
limited by the following factors:
\begin{itemize}
\item Available parallelism. 
\item Load balance: means that some processors work
  while others wait due to insufficient parallelism or unequal size 
  tasks. Few examples of unequal size tasks are:
\begin{enumerate}
\item The problems are fundamentally unstructured. 
\item Adapting to only the "interesting parts of the domain". 
\item Dealing with trees-structured computations. 
\end{enumerate}
\item Extra work of managing parallelism and dealing with redundant
  computation.
\item Communication.
\end{itemize}

\newpage

\subsection{Parallel Machines and Programming Models}

The key three models are:
\begin{enumerate}
\item {\bf Hardware-based}: Shared Memory Architecture. {\bf
    Software-based}: Shared Address Space.
\item {\bf Hardware-based}: Distributed Memory Architecture. {\bf
    Software-based}: Message Passing.
\item {\bf Hardware-based}: SIMD (Single Instruction Multiple Data)
  Architecture. {\bf Software-Based}: Data Parallel.
\end{enumerate}

\begin{figure}[tbp]
\centerline{\psfig{figure=arch-models.ps,height=355pt,width=470pt}}
\vspace{0.20cm}
\caption{
Illustrates the high-level mapping of these three programming models
to the hardware. (Mem: memory, CA: communication assistant, \$: cache,
and P: processor).}
\label{models}
\end{figure}

\newpage
\section{Shared-Memory Multiprocessing}

There are primarily two ways to share memory in a multiprocessing
environment: through centralized shared memory or through distributed
shared memory. In the centralized approach, which used in symmetric
multiprocessor (SMP) systems and is usually hardware-based, each
processor is connected to a shared bus. This shared bus handles all
accesses to main memory and I/O. In the distributed approach, each
processor has its own local memory and I/O system, and the processors
communicate by messages passed through a switched interconnection
network. Because of this architecture, distributed shared memory
systems are usually software-based.

The centralized approach is often simply called {\it shared memory},
while the distributed approach is often called {\it message
  passing}. In both approaches (as in any modern architecture), each
processor has a local cache, which as we will see, introduces
interesting problems with coherence in a shared memory environment. 
See Figure~\ref{centralized} and Figure~\ref{distributed} for
schematics of the two approaches.

\begin{figure}[tbp]
\centerline{\psfig{figure=shared-memory.ps,height=285pt,width=385pt}}
\vspace{0.20cm}
\caption{
Illustrates the approach of centralized shared memory in a symmetric
multiprocessor system (SMP), where each processor accesses memory
through a shared bus. (Processors are denoted by P, caches by \$, and
memory by M).}
\label{centralized}
\end{figure}
%%\vspace{1.0cm}
\begin{figure}[tbp]
\centerline{\psfig{figure=message-passing.ps,height=285pt,width=390pt}}
\vspace{0.20cm}
\caption{
Illustrates the approach of distributed shared memory, where each
processor has local memory and processors communicate through
interconnection network. (processors are denoted by P, caches by \$,
and memory by M).}
\label{distributed}
\end{figure}

In the future, it is likely that multiprocessing architecture will
consist of clusters of SMPs, that is, a message passing system in
which each ``processor'' is actually an SMP.

There are some key tradeoffs between shared memory and message passing
systems:
\begin{itemize}
\item Interface to Communication 
\begin{itemize}
\item Shared Memory Systems: Communication between processors is
  implicit and transparent. Processors access memory through the
  shared bus.  
\item Message Passing Systems: Processors must explicitly communicate
  with each other through messages.  
\end{itemize}

\item Complexity of the Architecture
\begin{itemize}
\item Shared Memory Systems: Leverages conventional architecture
  better since existing processors can be added to the shared bus
  system easily.
\item Message Passing Systems: Since there are fewer assumptions on
  the model, it leads to a simpler multiprocessing architecture
  overall. The catch is that this requires code to be rewritten for
  new platforms due to the explicit interface to communication.
\end{itemize}

\item Convenience
\begin{itemize}
\item Shared Memory Systems: Serial code runs without modification.
\item Message Passing Systems: Message passing libraries (i.e. mpi,
  pvm,...) are available for a wide variety of platforms. 
\end{itemize}

\item Protocols 
\begin{itemize}
\item Shared Memory Systems: Processors do not explicitly communicate
  with each other so communication protocols are hidden within the
  system. This means communication can be close to the hardware (or
  built into the processor), with the shared bus system deciding on
  how most efficiently to manage communication.
  Since communication occurs as part of the memory system, a smart
  shared memory architecture can make communication faster by taking
  advantage of the memory hierarchy. 
\item Message Passing Systems: Communication protocols are fully under
  user control. These protocols are complex to the programmer causing
  communication to be treated as an I/O call for portability
  reasons. This can be expensive and slow.
\end{itemize}
\end{itemize}

\subsection{Cache coherence}

Recall that the purpose of a cache is to keep the freshest data in
fast memory in order to exploit {\bf temporal locality} in memory
access. Thus the overall goal of caching is to reduce communication
between a processor and the main memory, since this communication is
typically slow compared to memory accesses to a cache. When shared
memory architecture is introduced, however, new issues relating to
cache coherence emerge.

Suppose we have two processors, A and B, and a variable {\bf x} stored
in shared memory which initially has the value 17 and is present in
neither processor's cache. Consider the following sequence of events:
\begin{enumerate}
\item Processor A reads {\bf x}, causing a cache fault, so {\bf x} =
  17 is loaded into its cache.
\item Processor B reads {\bf x}, causing a cache fault, so {\bf x} =
  17 is loaded into its cache.
\item Processor A writes a value of 42 to its cache copy of {\bf x},
  without updating main memory.
\item Processor B reads {\bf x}, causing a cache hit, so it still sees
  a value of 17.
\item Processor A flushes {\bf x} from its cache, so it writes back
  the value of 42 to main memory.
\item Processor B reads {\bf x}, causing a cache hit, so it still sees
  a value of 17.
\end{enumerate}

The above sequence of events elucidates a couple of cache coherence
issues. At step4, processor B does not see the new value of {\bf x}
written at step 3 by processor A, because processor A only updates its
own cache copy. But even after processor A writes back the new value
to main memory in step 5, processor B still sees the old value in step
6 because it reads from its own cache. Figure~\ref{issues} represents
this sequence of events.

\begin{figure}[tbp]
\centerline{\psfig{figure=cache-coherence-prob.ps,height=400pt,width=400pt}}
\vspace{0.20cm}
\caption{Cache coherence issues: The problem is that processor B still
  see the old value of {\bf x} in step 6 even though processor A has
  written back the new value to main memory in step 5. (caches is
  denoted by \$)}
\label{issues}
\end{figure}

The problem highlighted above is one of cache coherence. Even after a
processor writes back and updates memory location to main memory,
another processor may not see the new value. One simple solution to
this problem is to discard caches altogether, but this leads to
unacceptable loss in performance. A more reasonable solution is to use
a {\bf coherence protocol} in which each line is marked with a state
and these states are used to ensure coherence.

For example, a simple protocol might include three cache states:
\begin{itemize}
\item INVALID. There is no usable data in the cache line.
\item SHARED. The processor only has permission to read the data in
  the cache line.
\item EXCLUSIVE. The processor has read and write permissions for the
  data in the cache line.
\end{itemize}

The protocol rquires that before a processor writes a cache line, it
must first obtain an exclusive copy. This is accomplished by
invalidating all other copies that may exist in the caches of other
processors. 

There are other protocol requirements, but without specifying these or
explaining yet how the protocol is implemented, let us reconsider the
previous example with processors A and B sharing a variable {\bf
  x}. This time, however, we will apply this simple coherence
protocol. Recall that {\bf x} is initially 17 and is not present in
either processor's cache initially.

\begin{enumerate}
\item Processor A reads {\bf x}, causing a cache fault, so 17 is
  loaded into its cache and the line is marked EXCLUSIVE since it is
  the only processor that has a cache copy of {\bf x} at this point.
\item Processor B reads {\bf x}, causing a cache fault, so 17 is
  loaded into its cache. Since processor A already has a copy of {\bf
    x}, both processors must mark their copies of {\bf x} as SHARED.
\item Processor A writes a value of 42 to its cache copy of {\bf x}
  without updating main memory. Processor B's cache copy of {\bf x} is
  invalidated and processor A's copy is marked as EXCLUSIVE.
\item Processor B reads {\bf x}, causing a cache fault (since its
  former copy of {\bf x} has been marked INVALID). It retrieves the
  new value of 42 from processor A and both copies are marked as
  SHARED.
\item Processor A flushes {\bf x} from its cache, so it writes back
  the value of 42 to main memory. Since processor B is now the only
  processor with a cache copy of {\bf x}, it marks its copy as
  EXCLUSIVE.
\item Processor B reads {\bf x}, causing a cache hit.
\end{enumerate}

Due to the invalidation of outdated cache copies in step 3 (and the
other mechanism of the protocol), processor B never reads an outdated
copy of x. Figure~\ref{solution} represents this sequence of events,
with coherence protocol in use.

\begin{figure}[tbp]
\centerline{\psfig{figure=cache-coherence-sol.ps,height=410pt,width=410pt}}
\vspace{0.20cm}
\caption{Cache coherence issues: The solution is to use a coherence
  protocol, which marks each line with a state (E for EXCLUSIVE, S for
  SHARED, and a dash(-) for INVALID) and invalidates processor B's
  copy of {\bf x} when processor A writes to it in step 3. (caches is
  denoted by \$)}
\label{solution}
\end{figure}

In general, cache coherence protocols are very complex. In real
implementations, a cache line can be in one of 1000's of different
states. Thus, it is no surprise that algorithm verification programs
are routinely used to evaluate their correctness.

\subsubsection{Snooping}

How would the above protocol be implemented in an SMP? In particular,
when processor B reads {\bf x} in step 4, how does it know to get the
new value from processor A's cache? After all, in shared memory
architecture, processors access memory through a shared bus, not by
message passing.

One general technique, called {\it snooping}, take advantage of the
fact that all memory accesses are serialized by the shared bus. With
snooping, all the caches of each processor listen to the bus and
handle memory requests if appropriate. In step 4 above, processor A's
cache notices that processor B is reading x. Since it has an exclusive
copy of x and therefore knows that it has the latest value, processor
A supplies the value x to the bus, preempting the main memory (which
contains stale data). Furthermore, both processor A and B know to mark
their copies of x as SHARED. See Figure~\ref{snooping} for an
illustration of this example.

\begin{figure}[tbp]
\centerline{\psfig{figure=snooping.ps,height=400pt,width=400pt}}
\vspace{0.20cm}
\caption{An example of snooping. Processor B issues a read to a
  location owned by processor A. Processor A ``snoops'' the read
  request from the shared bus and supplies the correct value,
  preempting the outdated value from main memory. Both cache copies
  are marked SHARED.}
\label{snooping}
\end{figure}

With snooping, the shared bus is treated as a broadcast communication
channel. Of course, this communication is handled at the architecture
level, so the user need not be aware of it besides relying on the fact
that there is in fact cache coherence.

\subsubsection{Directory-based coherence}

Another way to implement a coherence protocol is to use a directory,
where a central {\it home} unit keeps track of which memory locations
are being shared in multiple caches and which are held exclusively in
a single processor's cache. All memory accesses are handled by the
home, which knows based on its directory which caches need to be
updated or invalidated on any given memory access.

For instance, suppose a variable is being shared by multiple
processors, and one of the processors wants to write a new value to
the variable. In a directory-based coherence scheme, the writing
processor would issue a write request to home. Then the home, knowing
which other processors shared the variable, would send invalidate
messages to these processors in order to invalidate their copies of
the variable, and grant the exclusive copy to the requesting
processor. See Figure~\ref{directory} for an illustration of this
example. Figure~\ref{directories} shows a scalable multiprocessor with
directories.


\begin{figure}[tbp]
\centerline{\psfig{figure=directory.ps, height=390pt,width=385pt,angle=-90}}
\vspace{0.20cm}
\caption{Directory-based coherence. In this example processor 1 issues
  a write request, causing home to invalidate the processors which
  have a cache copy of the memory location. The home unit updates its
  directory, noting the fact the now processor 1 has an exclusive
  copy, and that the other processors no longer have shared copies.}
\label{directory}
\end{figure}

\begin{figure}[tbp]
\centerline{\psfig{figure=directories.ps, height=385pt,width=425pt}}
\vspace{0.20cm}
\caption{Scalable multiprocessor with directories. Every block of main
  memory, the size of a cache block, has a directory entry that keeps
  track of its cached copies and their state.}
\label{directories}
\end{figure}

Directory-based coherence is used by all software-based
implementations of shared memory. One problem associated with
directory-based coherence, however, is that it is particularly
sensitive to the interconnection network (IN) speed. Furthermore,
because all memory operations are handled centrally, latency is
increased for these operations.
 
\subsubsection{False sharing}

Suppose we have two processors, A and B, and that processor A is
accessing variable {\bf x} and processor B is accessing variable{\bf
  y}. Since the two processors are accessing different variables, they
should each obtain and maintain exclusive copies of these variables,
correct? Not necessarily. Typical caches contain more than a single
word per cache line. So if {\bf x} and {\bf y} reside on the same
cache line, then the two processors will be sharing a single cache
line though they are not actually sharing the same variable. 

Since the coherence protocols are based on cache lines, a write by
processor A to {\bf x} will cause processor B's copy of the shared
cache line to be invalidated. And likewise a write by processor B to
{\bf y} will cause processor A's copy of the shared cache line to be
invalidated. This unfortunate side effect of using coherence protocols
with multiple-word cache lines, which leads to unneeded invalidations,
is called {\bf false sharing}. See Figure~\ref{sharing} for an
illustration of this phenomenon.

\begin{figure}[tbp]
\centerline{\psfig{figure=sharing.ps,height=410pt,width=410pt,angle=-90}}
\vspace{0.20cm}
\caption{Flase sharing: (a) When {\bf x} and {\bf y} do not reside on
  the same cache line, the two processors acquire and maintain
  exclusive copies. (b) When {\bf x} and {\bf y} share a cache line,
  however, false sharing can occur since writes to one variable cause
  the other variable to be invalidated since they lie on the same
  cache line.}
\label{sharing}
\end{figure}

Depending on how a program accesses memory with respect to cache lines
and time, false sharing can depreciate the performance of a parallel
program significantly. Currently, it is up to the programmer to detect
false sharing and to take steps to alleviate the problem. The simple
solution to false sharing is to split the data across cache lines,
which can be done manually by a knowledgeable programmer.

\subsection{Memory consistency}

The {\it memory consistency model} of a shared memory system specifies
how memory will behave in the context of multiple processors reading
and writing shared locations. The most basic consistency model, based
on the intuition of a single-processor system, is 
{\it sequential consistency}. A system is sequentially consistent if
the following two conditions are met:
\begin{enumerate}
\item The result of any execution is the same as if the operations of
  all the processors were executed in {\it some} sequential order.

One can think of sequential order as a master controller picking one
of the processors to run next while halting all of the others. This
way, threads can be interleaved, yet only one processor is active at a
time.

\item The operations of each individual processor appear in this
  sequence in the order specified by the program.
\end{enumerate}

Sequential consistency gives the programmer a straightforward,
intuitive model of how multiple processors interact regarding shared
memory. Unfortunately, it is difficult to implement sequential
consistency cheaply.

There are four main types of consistency models:
\begin{enumerate}
\item R $\longrightarrow$ R
\item W $\longrightarrow$ R
\item R $\longrightarrow$ W
\item W $\longrightarrow$ W
\end{enumerate}

The operations in the first, a read followed by a read, can be
performed in any order. For the last three combinations, however,
order must always be preserved if the operations are performed to the
same location. As we shall see, there are various synchronization
methods used to ensure the proper ordering of these operations.

\subsubsection{Relaxed consistency}

Programs can often be made to run faster by allowing memory operations
to complete in a different order from that specified by the
program. For example, {\it total store order (TSO) consistency}
(see ref) allows a read to bypass an already issued write if the two
operations access different locations: given that {\bf x} and {\bf y}
are distinct memory locations, a write to {\bf x} can execute after a
read to {\bf y} even if the read is prior to the write in the program
order, since there is no semantic difference between the two orderings
when run on a single processor. This is an example of a {\it relaxed
  consistency mode}, in which the requirements of sequential
consistency are relaxed in order to realize higher performance.

\subsubsection{Relaxed consistency anomalies}

Unfortunately, relaxed consistency can lead to unexpected behavior
with multiple processors. For example, consider the following
situation:
\begin{itemize}
\item Two variables in shared memory, {\bf x} and {\bf y}, are both
  initially 0.
\item Processor A runs the following code:

\begin{verbatim}
       x = 1;
       if (y == 0)
         foo();
\end{verbatim}

\item Processor B runs the following code:

\begin{verbatim}
       y = 1;
       if (x == 0)
         foo();
\end{verbatim}
\end{itemize}

How many times can the procedure foo() be run? It depends on the
consistency model. Under sequential consistency, foo() can be run
either zero one times. It will run zero times if both x and y are
set to 1 before either processor executes the if statement. To see
that foo() will be run at most once, suppose processor A calls
foo(). Then y must still be be 0 and so processor B has not yet
assigned 1 to y, and thus, due to the execution order requirement of
sequential consistency, processor B has not yet tested the value of
x. When processor B does read x, it will see a value of 1, since
processor A must have set x to 1 before calling foo(), again due to
the ordering requirement of sequential consistency.

With total store ordering, however, foo() could be executed twice --
once by each processor! This can happen because TSO will allow the
reads (which test whether x or y equal 0) to be executed before the
writes (which set x or y to 1), despite the specified program
order. After all, x and y are different memory locations, so why
shouldn't these operations be executed out of order if it leads to
higher performance? As we see, the reason is that a reordering that is
semantically identical in a single processor situation is not
necessarily so with multiple processors. (Note that instruction
reordering may be performed by a compiler or even by the hardware
itself).

A solution to the above anomaly, where foo() can be executed twice
despite the intention of the code to use a {\it memory fence}. A
memory fence is an instruction around which memory accesses are not
allowed to be reordered. Any memory operation before a fence must be
executed before any operation after the fence. Operations between
fences can still be reordered according to TSO (or whatever other
relaxed consistency model in use).

\subsection{Locking}

A {\it lock} is a synchronization variable that can be acquired by
only one processor at a time. If a processor is unable to obtain a
lock, it {\it blocks}, waiting for the lock to be freed. Locks are
used to guarantee {\it mutual exclusion}, the condition that a
fragment of code will not be executed simultaneously by more than
processor. For example, if two processors A and B both run the
following code:

\begin{verbatim}
   lock(S);
   foo();
   unlock(S);
\end{verbatim}

where S is a shared lock, then foo() (the {\it critical section}) will
be run by at most one of the two processors at any given time.

\subsection{Atomic exchange}

A very common way to implement locks is with {\it atomic exchange}
instruction, which is provided by many machines. The semantic of
atomic exchange are simple:
\begin{itemize}
\item The values of two memory locations are exchanged.
\item The exchange occurs atomically, that is, other processors see
  either the state of the variables before or after the exchange,
  never an intermediate state where the exchange is only partially
  completed.
\end{itemize}

Given this instruction, locking can be implemented easily:

\begin{itemize}
\item The shared lock variable S is initially 1, indicating that the
  lock has not been acquired by any processor.
\item Each processor contending for the lock has a private variable s,
  initially set to 0, which indicates whether that processor owns the
  lock.
\item When locking, a processor does an atomic exchange between the
  shared lock variable S and its private indicator s. If S equals 1
  and therefore the lock is free to be acquired, the instruction will
  set s = 1 and S = 0, giving the locking processor the lock and
  indicating to other processor that the lock is no longer free. If S
  equals 0, the processor blocks by looping until S becomes 1 and gets
  exchanged into s. The code for lock(S) is

\begin{verbatim}
  While (s != 1)
    AtomicExchange(s, S);
\end{verbatim}

\item When unlocking, a processor simply does a single atomic exchange
  between S and s, which appropriately sets the shared lock variable
  back to 1 and the private indicator back to 0. The code for
  unlock(s) is the one instruction:

\begin{verbatim}
   AtomicExchange(s, S);
\end{verbatim}
\end{itemize}

\subsubsection{Load-linked / store conditional}

A more modern approach to implementing locks is based on two new
instructions: {\it load linked} (LL) and {\it store conditional}
(SC). A store conditional is a store from x to y which completes the
store only if no change has occured to y since the last time y was
loaded with the load linked instruction, with x being set to 0 or 1 to
indicate whether the store occured or not. 

For example, consider the following code which implements an atomic
exchange between the location pointed to by R1 (written as *R1 and the
value of R4):

\begin{verbatim}
 try:  MOV    R3, R4
       LL     R2, *R1
       SC     R3, *R1
       BEQZ   R3, try
       MOV    R4, R2
\end{verbatim}

The code works as follows:

\begin{enumerate}
\item The value of R4 is assigned to a temporary register R3.
\item The load linked instruction is used to load *R1 into R2.
\item The store conditional instruction is used to store R3 (which
  equals R4) into *R1, with two possible results: *R1 has not changed
  since the LL instruction, so *R1 gets R3 and R3 gets 1, indicating a
  successful store; or *R1 has changed since the LL instruction (due
  to another processor), so *R1 is unchanged and R3 gets set to 0,
  indicating the store failed.
\item If the store failed (i.e. RC == 0), then the program tries again
  by branching to the beginning.
\item Otherwise, the exchange is completed by storing the value of R2
  (which equals *R1) into R4.
\end{enumerate}

Notice that this code depends on one of the exchanged variables (R4)
being a private local register, which cannot be accessed by other
processors. 

The main advantage to using the load-linked/store
conditional scheme for implementing locks is that it exploits {\it
  speculative computing}. Speculative computing is when a processor
proceeds to execute subsequent instructions even when the work may
need to be aborted depending on later results. Improved performance is
realized when the processor does not need to abort. If it must abort,
then nothing has been lost since the processor does not have any other
useful work to do.

In the case of pipelining, for example, when a branch instruction is
reached, a processor will choose one of the two branches to follow. If
the wrong one is chosen, the instructions after the branch are flushed
from the pipeline and the processor backtracks to process the correct
branch. The only other alternative is to execute NOP's until the
branch instruction is evaluated, which never reduces processing time.

Similarly, in the case of load-linked / store conditional, the
procedure proceeds assuming that the store will be successful. If it
fails, then the code simply aborts to the original state before the
exchange attempt. If it succeeds, then it can complete the locking
process quickly and proceed with subsequent operations.

\subsubsection{Deadlock and the dining philosophers problem}

Once upon a time, before the guilded age of computer science and keen
insights into locking, there were five philosophers sitting at a round
table. In front of each philosopher was a bowl of rice, and also
each had only one chopstick. Deciding that the most enlighted solution
to this dilemma was to share chopsticks in a systematic way, each
philosopher placed his chopstick to the left of his bowl, resulting in
a single chopstick between every two adjacent bowls of rice. Then they
agreed to the following rule: when someone wanted to eat, he would
first grab the chopstick to his left, then the one to his right, and
finally replace the two chopsticks after satisfying his hunger.

For many moons the philosophers thought their grat thoughts, eating
whenever they felt hungry. (The ancient traditions do not tell us why
the bowls of rice never emptied ...) One day, the five philosophers
had each picked up the chopstick to his left at the very same moment,
and then, seeing that the chopstick to his right was not there, begain
waiting for it to be replaced. Unfortunately, no chopstick was ever
put back on the table, since each philosopher was waiting for the
philosopher on his right, thus no one could ever get two chopsticks to
eat his rice. This is what is referred to as deadlock, hence,
starvation. 

What could the philosophers have done to prevent such an unfortunate
end? They could have first numbered the chopsticks from 1 to 5, then
made the rule that when one of them wanted to eat, he would first pick
up the chopstick next to him with the lower number, and then pick up
the other chopstick. Such an algorithm would have prevented deadlock
because in any situation in which the philosophers contend for the
same set of chopsticks, the person who manages to pick up the highest
numbered chopstick will always be able to proceed with eating, and
will thus eventually replace his chopsticks, allowing others who are
waiting for his chopsticks to also proceed.

The moral of the dining philosophers story is that to prevent
deadlock, locks need to be ordered, thus breaking the symmetry in the
lock-acquiring algorithm. This allows at least one processor to
acquire all of the required locks in order to continue processing.

Relating the story to processors and shared resources, the problem of
deadlock can be prevented by numbering the shared resources (e.g the
chopsticks), and then having each processor (e.g. the philosophers)
acquire locks (e.g. pick up a chopstick) according to a specific total
ordering, with lowest first. Thus, the invariant is that the processor
which acquires the highest numbered resources (according to the total
order) will always be able to proceed.

\subsection{Nonuniform memory access}

One possible approach to build a scalable shared memory system is to
maintain the uniform memory access or what is called the ``dancehall''
and provide a scalable interconnect between the processors and
memory. Figure~\ref{dancehall} illustrates this approach.

\begin{figure}[tbp]
\centerline{\psfig{figure=dancehall.ps,height=285pt,width=395pt}}
\vspace{0.20cm}
\caption{Dancehall multiprocessor organization. Even if there is no
  communication or sharing between the processors, the bandwidth
  requirement on the interconnect increases linearly with the number
  of processors.}
\label{dancehall}
\end{figure}

The main disadvantage is that latencies incurred by accesses to memory
would require large bandwidth to be supplied to every processor in the
network. 

An alternative approach is interconnecting a set of complete
processors, each with not only its own cache but also local memory. In
other words, each node is a UMA SMP. This approach is referred to as
NonUniform Memory Access (NUMA). Figure~\ref{NUMA} illustrates a
simple case of two-node NUMA, each is a two-processor SMP.

\begin{figure}[tbp]
\centerline{\psfig{figure=CC-NUMA.ps,height=420pt,width=250pt,angle=-90}}
\vspace{0.20cm}
\caption{A simple CC-NUMA system: consisting of two nodes each is a
  UMA SMP. CA: denotes communication assist, which acts as local
  surrogates for all other nodes of the system. They redirect memory
  requests appropriately across an interconnect (NI) to other nodes.}
\label{NUMA}
\end{figure}

\begin{itemize}
\item In this approach, every processor has access to all of memory
  using normal loads and stores. However, there is a difference from
  the standard SMP model and that is there will be ``access delay''
  depending on what parts of memory are accessed; hence the
  ``nonuniform'' part.
\item From the above figure, in general, in each node 
  access to local memory is faster than access to remote memories in 
  other nodes.
\item The ability to access local memory quickly does not effect much,
  if at all, the access time of remote data.
\item I/O systems may either be part of every node or consolidated
  into special I/O nodes.
\item The ``N'' (of nonuniform) in NUMA is measured as the ratio of
  access time to remote memory to access time to local memory.
\item The NUMA approach is widely used because of its inherent
  performance advantages and it makes a better use of the mainstream
  processor memory technology.
\item Systems such as SGI Origin2000 was designed as NUMA system, with
  built-in hardware to support cache coherency. Its architecture is
  referred to as cache coherent NUMA or CC-NUMA. It has an ``N'' ratio
  of (3 to 2):1.
\end{itemize}

\subsection{Some concluding remarks}

Shared memory architectures are becoming the standard approach to
multiprocessing, with every major manufacturer offering an SMP
product. As we have seen, such architectures can be tricky for the
user, especially regarding memory consistency and issues of locking. A
great step forward in this regard would be high-level programming
support which frees the programmar from having to use user-level
protocols to manage these complexities.

\subsection{Hardware Examples}
The Intel Pentium Pro 4-processor "quad pack" is a good example 
of a shared memory (centralized) architecture.

\subsubsection{Pentium Pro 4-processor -- quad pack}

\begin{figure}[tbp]
\centerline{\psfig{figure=quad-pack.ps,height=400pt,width=400pt}}
\vspace{0.20cm}
\caption{A rough sketch of the Pentium Pro 4-processor quad pack
  organization.}
\label{quad}
\end{figure}

Figure~\ref{quad} illustrates the organization of a Pentium Pro
motherboard used in many multiprocessor servers nowadays. 

\begin{itemize}

\item It accommodates up to 4 processor modules, each containing a
  Pentium Pro processor, first-level cache, translation lookaside
  buffer, a 256-KB second-level cache, an interrupt controller, and a
  bus interface in a single chip connecting directly to a 64-bit
  memory bus.
\item The bus operates at 66 MHz, and memory transactions are
  pipelined to give a peak bandwidth of 528 MB/sec. 

\item A two-chip memory controller and 4-chip memory interleave unit
  (MIU) connect the bus to multiple banks of DRAM. Bridges connect the
  memory bus to 2 independent PCI buses, which host display, network,
  SCSI, and lower-speed I/O connections.

\item The Pentium Pro module contains all the logic necessary to
  support the multiprocessor communication architecture, including
  that required for memory and cache consistency.

\item The structure of the Pentium Pro "quad pack" is similar to a
  large number of earlier SMP designs but shows an expanded view of a
  typical Pentium Pro SMP.
\end{itemize}



\section{Message Passing}

Message passing is a parallel programming style used typically on
distributed memory machines. Processors communicate using messages,
and it is the most popular parallel programming paradigm. 
One reason for this popularity is the wide number of platforms which
can support a message-passing model. Programs written in a
message-passing style can run on distributed or shared memory
multiprocessors, and network of workstations. The point of having the
paradigm, just as in the sequential case, is that the programmer knows
that his or her algorithms should in principle be portable to any
architecture that supports a message passing model\footnote{In
  reality, just as in the sequential case, any programmer who cares
  about performance should give some thought to the target
  architecture.}. MPI and PVM are the two message passing libraries
that are  most widely used.

Some of the parallelism issues that apply to the shared memory model
also apply to message passing. For example, deadlocks, which typically
involve a processor waiting to recieve a message that never arrives,
and serialization. Other issues such as locks and global queues don't
apply. Load balancing is usually done in an explicit fashion to take
into account the slow communication speeds between nodes (each will
have at least one processor, its own memory, communication
interconnect, and possibly its own I/O devices): typically
tens or handreds of megabits per second verses Gigabytes per second 
within SMPs.

\subsection{Generic Message Passing}

Message passing libraries support basic primitives to send and receive
data through complementary operations. Also, they support a machine
model in which each pair of processors have a virtual communication
channel between them.

\subsubsection{Typical Primitives}

\begin{verbatim}
 void send (processor, void *msg, size)
 int  recv (processor, void *msg_buffer, max_buffer_size)
\end{verbatim}

The two most widely used operations are blocking send and blocking
receive. Blocking send does not return control until it is safe to
modify the message buffer, and blocking receive does not return
control until the message is received. 

These type of communication primitives are also referred to as
point-to-point communication since only the sending and the receiving
processors need to know anything about the message. There are several
variations on how the sending of a message can interact with the
execution of the sub-program. 

There are two important aspects or distinctions in the message passing
model. The first common distinction is between {\it synchronous} and {\it
  asynchronous} sends. Synchronous sends are provided with information
about the completion of the message, while asynchronous sends only
know when the message has left. The other important distinction is
between {\it blocking} and {\it non-blocking}. Blocking operations
only return from subroutine call when the operation has
completed. Non-blocking operations return straight away and allow the
sub-program to continue to perform other work. At some later time the
sub-program can test for the completion of the non-blocking operation.

\subsubsection{Pseudo-implementation}

\begin{itemize}
\item For the receive operation:

\begin{verbatim}

 int recv (processor, void *msg_buffer, max_buffer_size)
   while channel from processor empty
     do  nothing
   if max_buffer_size < msg_size
     then  error
 read message from channel into msg_buffer
   return msg_size

\end{verbatim}

\item For the send operation:

\begin{verbatim}

 send (processor, void *msg, size)
   while channel to processor busy
     do  nothing
   while recv not called by processor
     do  nothing
 write size, msg into channel
   return

\end{verbatim}

\end{itemize}

The above is a {\it synchronous} blocking send, the most common type
of send operations used in message passing libraries. It does not complete
until a recv initiated by the receiving processor. 

\subsection{Message passing programming}

Message passing programs are usually of the SPMD (Single Program,
Multiple Data) style, refering to a restriction of the message-passing
model which requires that all processes run the same executable. Some
vendors provide parallel environments which only support SPMD parallel
programs. 

Here is a simple message passing example which exchanges data between
two processors, 0 and 1:

\begin{verbatim}
   if self = 0
       then  send (1, &data, 100)
             recv (1, &data, 100)
   if self = 1
       then  recv (1, &data, 100)
             send (1, &data, 100)
\end{verbatim}

Notice that if the send and the recv in the second {\bf if} are
switched, this program will deadlock. 

Synchronous blocking communication is the most common protocol for
sending large amounts of data. It up to the programmer to explicitly
program both ends of the communication and take care to avoid any
deadlock.

\subsection{Collective Communication}

Message passing libraries also support global, collective
communication operations such as broadcast, reduce, scan, and
barriers.

\subsection{Summary}
The reasons behind the success of message passing are:
\begin{itemize}
\item Message passing programs are portable, and the model makes few
  hardware assumptions.
\item It provides a simple model and it is not difficult to get a
  simple message passing program working.
\item One possible solution to deadlocks is to increase buffer sizes.
\item It exploits hierarchical memory nicely.
\end{itemize}

\subsection{Hardware Examples}

The IBM SP-2 represent a good example of a message passing
architecture.

\subsubsection{IBM SP-2}

\begin{figure}[tbp]
\centerline{\psfig{figure=SP-2.ps,height=420pt,width=430pt}}
\vspace{0.20cm}
\caption{Diagram of a single node of the SP-2.}
\label{SP2}
\end{figure}

A summary of the machine's basic characteristics: 
\begin{itemize}

\item Constructed essentially out of almost complete RS6000
  workstations.
\item Some modifications are made to package the RS6000s into standing
  racks.
\item The rack-level topology varies from machine to machine but it is
  typically a variant of a butterfly network.
\item The physical topology is not fixed in hardware so the network
  interface inserts the route for each outgoing message via a table
  lookup.
\item Since a complete system is used for the node, the communication
  assist and network interface (NI) are part of a card that plugs into
  the system. This is a MicroChannel adapter card connecting to a
  switch in the base of the rack.
\item Figure~\ref{SP2} illustrates in some details the organization 
  of an SP-2 node. As shown, a network interface card (NIC) connected
  to the MicroChannel I/O bus. The NIC contains the drivers for the
  actual link into the network, a large amount of memory to buffer
  message data, a direct memory access (DMA) engine, and a complete
  i860 microprocessor to move data between host memory and the
  network.
\item Not shown in the figure is the network itself which is a
  butterfly-like structure, constructed by cascading 8 x 8 crossbar
  switches. The individual links operate at 40 MB/sec in each
  direction, which is the full capability of the I/O bus.
\item Between 8 to 16 nodes fill a rack, so large configurations are
  built by cabling together a number of racks, including additional
  switch racks.
\item With a complete system at every node, there is the option of
  distributed disks and other I/O devices over the entire
  machine. Disks on compute nodes are typically used only for 
   swap space and temporary storage. Most of the I/O devices are
   concentrated on dedicated I/O nodes. 
\item This style of design allows for some degree of heterogeneity
  among the nodes since all that is required is that the network
  interface card can be inserted. Therefore, SP systems support
  several models of workstations as nodes, including SMP nodes. 
\end{itemize}

\section{Data-Parallel Computing}

In this paradigm, the parallelism lies in the data, not in the
control. Also, it relies on the {\it vector model} of
parallelism. This model assumes multiple data and serial
control. Figure~\ref{VRAM} depicts the vector RAM (VRAM)
architecture. Notice that the program resides in the scalar memory and
that there is one program and one program counter. This model is also
referred to as the SIMD model (Single Instruction Stream, Multiple
Data Stream), a term that was coined by Flynn in the late
sixities. See Figure~\ref{SIMD} for a schematic view of the model.
A competing model, that was mentioned in the shared memory section, is
the MIMD model (Multiple Instruction Stream, Multiple Data Stream). 

\begin{figure}[tbp]
\centerline{\psfig{figure=VRAM.ps,height=250pt,width=400pt}}
\vspace{0.20cm}
\caption{Diagram of a vector architecture. The program resides in
  scalar memory.}
\label{VRAM}
\end{figure}

\begin{figure}[tbp]
\centerline{\psfig{figure=SIMD.ps,height=260pt,width=410pt}}
\vspace{0.20cm}
\caption{A schematic view of the SIMD architecture supporting the
  data-parallel model.}
\label{SIMD}
\end{figure}

\subsection{Vector Operations}

In this section we present the major classes of VRAM instructions. We
consider vectors A and B, where:

\begin{verbatim}
 A = [1 2 3],   B = [4 5 6]
\end{verbatim}

{\bf Elementwise Operations}: execute one operation on all vector
elements.

\begin{verbatim}
Example:  A + B --> [5 7 9].
Example:  A * B --> [4 10 18].
Example:  select(A, B, [1 0 1]) --> [1 5 3].
\end{verbatim}

{\bf Reduction Operations}: take vectors as arguments and output
scalars.

\begin{verbatim}
Example: size(A) --> 3.
Example: sum(A)  --> 6.
\end{verbatim}

{\bf Distribute Operations}: take scalars as arguments and return
vectors.

\begin{verbatim}
Example: dist(3,5) --> [5 5 5].
\end{verbatim}

{\bf Routing Operations}: move data among processors.

\begin{verbatim}
Example: B[2 1 3] --> [5 4 6].
Example: B[2 1 2] --> [5 4 5].
\end{verbatim}

Notice that executing B[2 1 2] is a qualitatively harder operation
than executing B[2 1 3], because in B[2 1 2] communication is
unbalanced. More specifically, B[2 1 2] requires a {\it concurrent
  read}. That is, several processors simultaneously read the value of
the second element of vector B. The implementation of concurrent reads
in hardware affects performance substantially.

\begin{verbatim}
Example: C[2 3 1] = B  ==> C --> [6 4 5].
Example: C[2 1 2] = B  ==> C --> ??
\end{verbatim}

The operation C[2 3 1] = B involves a write, but the operation C[2 1
2= B invloves a {\it concurrent write}. In a concurrent write,
multiple processors simultaneously attempt to write to the same memory
cell.

We can define the memory value after a concurrent write in several
ways. If concurrent writes are {\it illegal}, then on a concurrent
write, the computer overheats, catches fire, and blows up :). On
arbitrary concurrent write, the processor with the highest priority
succeeds in writing its value. On a {\it reduction} concurrent write,
all values written to any memory location are {\it combined} together
-- e.g. added, anded, or multiplied.

{\bf Scan Operations}: Scan operations are also called {\it parallel
  prefix} operations.

\begin{verbatim}
Example: +-scan(A) --> [1 3 6].
Example: *-scan(B) --> [4 20 120].
\end{verbatim}

{\bf Multidimensional Operations}: operations on multidimensional
arrays are multidimensional operations. Assume:
 
\begin{verbatim}
       | 1  2 |
       |      |
   C = | 3  4 |
       |      |
       | 5  6 | 

Example: row-sum(C) --> [3 7 11].
\end{verbatim}

{\bf Ordinary Serial Control}: includes operations such as if, while,
procedure calls, etc.

\subsection{Algorithmic Measures}

In this section we look at several measures that gauge the efficiency
of parallel algorithms. 

The {\it work} of a parallel program is the total number of operations
performed collectively by the processors. Alternatively, it is the
total time necessary to execute the parallel program on {\it one}
processor. The work is approximately the sizes of the vectors
manipulated summed over all computation steps.

The {\it critical path length} of a computation is the total number of
vector operations performed. Alternatively, it is the time necessary
to execute the parallel program on an {\it infinite} number of
processors.

A third measure is the {\it computation time on P processors}. The
following extremely important theorem, called Brent's theorem, relates
these three measures.

\bt{Brent-Graham}
(Due to Brent and Graham) Consider an algorithm having work $W$ and
critical path length $C$. On $P$ processors, one can achieve execution
time on $P$ processors $T_P$, where

$\;\;\;\;\;\; T_p \le {W \over P} + C.$
\ets

The proof is quite simple and will be left as an exercise for the
reader. 

{\bf Falsehoods in Simplification}. Several falsehoods derive from the
VRAM abstractions. Some of the operations listed in the previous
section are more work-consuming than others. many of the operations,
such as routing operations, have ``hidden logs.'' That is, routing
requires a time logarithmic in the number of processors.

\section{Active Messages}
Active Messages are a primitive communication mechnism which can be
used to implement traditional parallel programming paradigms such as
shared memory of MPI. It is a protocol layer intended to lie close to
the hardware and provide a flexible abstraction to allow programming.
The basic concept behind Active Messages is quite simple. Messages are
passed with a header indicating a user-level routine which extracts
and integrates the data from the network. At this level, the messages
are not buffered.

\subsection{Network Interface}
The network interface hardware which Active Messages relies upon
consists of two FIFOs which receive and send messages to and from the
network and a set of status bits (See
Figure~\ref{active-messages1}). The status bits represent the state of
the messages in the send and receive FIFOs. Typically, the network
interface is connected to the I/O bus, the cache bus, or integrated
with the processor registers. Some network interfaces are associated
with direct memory access (DMA) engines which process them directly to
memory.

Messages are sent and received through the network interface via a few
simple steps. To send a message, the status bits for the send FIFO are
checked to determine if it is free. Then the message is written into
the send FIFO and dispatched. To receive a message, the status bits of
the receive FIFO are checked to determine if a message has arrived. If
so, the message is read from the receive FIFO, and the FIFO is
subsequently cleared. To ensure proper transport of the messages, the
active message protocol adheres to the following contract coined the
consumption assumption.

{\bf Consumption Assumption}: The network promises to deliver all
messages eventually if processors promise to receive every message
eventually. 

This assumption translates into the requirement that message sending
routines cannot block. Otherwise, situations can arise to violate this
assumption resulting in deadlock. 

\subsection{Program Structure}

Active messages is flexible, low-kevel protocol layer. Its message
passing semantics are essentially restricted forms of remote procedure
calls (RPC). A data packet is attached with a pointer to a handler,
some remote procedure, which is invoked when the message arrives. The
remote procedure's task is to quickly remove the message data from the
receive FIFO and incorporate it into the computation with minimal
work. The handler is run on top of the user's stack. There is no
separate space for handlers or buffering of messages. To facilitate
the sender in specifying the address of the handler, programs are
implemented within the context of the SPMD environment.

In the SPMD model, identical copies of a program run on all
processors. Usually, the same memory map is utilized for both code and
data to simplify references to each other's variables and
handlers. However, in some implementations, the pointers to handlers
are indexes to the addresses of remote functions. The following
example in which a split-phase remote fetch is performed illustrates
the structure of a typical parallel program implemented with active
messages.

\begin{verbatim}

Example: y = processor 42's value of x.

Main function:
-------------
 main()
 {
  if (self != 1) {
    while(1)
       poll();
  }
  else {

   arrived = FALSE;

 /*  (a)  (b)  (c)   (d) */

 rpc(42, &foo, self, &bar);

 /* (a): destination processor id */
 /* (b): request handler          */
 /* (c): caller's processor id.   */
 /* (d): reply handler            */

  ...
  ...
 while (!arrived) poll();
  ...
  ...
 }
}

Handlers:
--------
 foo(m, *reply_handler) {
   rpc(m, reply_handler, x);
 }

 bar(v) {
   y = v;
   arrived = TRUE;
 }

\end{verbatim}

Processor 1 sends a message tp invoke the foo() handler on processor
42, which returns the value of x. Upon receipt of x at processor 1,
the bar() handler is invoked to set the values of y and the {\it
  arrived} flag, indicating the transaction is complete. The execution
of the above code is depicted in Figure~\ref{active-messages2}.

\begin{figure}[tbp]
\centerline{\psfig{figure=active-messages1.ps,height=300pt,width=270pt,
    angle=-90}}
\vspace{0.20cm}
\caption{Active messages: network interface hardware.}
\label{active-messages1}
\end{figure}

\begin{figure}[tbp]
\centerline{\psfig{figure=active-messages2.ps,height=300pt,width=270pt,
    angle=-90}}
\vspace{0.20cm}
\caption{Active messages: execution of example program.}
\label{active-messages2}
\end{figure}

\subsection{Sending and Receiving Messages}
...

...


\section{Programming Examples}

\subsection{Java -- MPI}

\subsubsection{The problem}
Computational Fluid Dynamics (CFD) 2-D Inviscid flow simulation. 

\subsubsection{Issues of performance}

\begin{itemize}
\item Parallelism was achieved by adding
message-passing (MPI) to the sequential version with very little
modifications. Message-passing was initiated between methods
computing artificial dissipation, where 2nd and 4th order derivatives
are computed. Those are the methods that are used continuously
throughout the simulation. Because of that, we thought that those are
the only methods in the code that could benefit from parallelism. We
confined our message-passing to simple send/receive and bcast. 

\item The parallel version benefited greatly from 
the efficient and well-designed sequential version, in which inlining
and other tune ups were already done. So the result was a robust
parallel version. 

\item  The inlining helped great deal to speed up
execution without creating any degradation to the quality of the
solution. Depending on the network and computation nodes, we have
noticed that the Java/MPI version was running close to even faster than
the CFD HPF (high performance fortran) version using the same number
of processors. In addition, we noticed that it definitely ran much
faster on a dedicated network of nodes. But the architecture of these
nodes are different than those used for the HPF tests.

\item There is still more room for improvements. In
addition to inlining, JARing can greatly help improve the speed of
class downloading. Another thing to do is to use thread pooling. This
is basically to create a ready supply of sleeping threads at the 
beginning of execution. Because of the thread startup process is
expensive in terms of the system resources, thread pooling makes the
startup process a little slower, but improves runtime performance
because sleeping (or suspended) threads are awakened only when they
are needed to perform new tasks.

\item Method inlining enhances performance by reducing the number of
  method calls the program makes. Other things that can be done is
  perhaps to try streamlined synchronization, JIT, and other
  performance enhancing approaches. These are all quite useful tools
  to try specifically for large scale applications such as
  computational fluid dynamics.
\end{itemize}

Here is the full listing of the source code:


\begin{verbatim}

/**********************************************************************/
/* CFD.java                                                           */
/* The Perallel version.                                              */
/* Sang Lim and Saleh Elmohamed   (1998)                              */
/**********************************************************************/
/* This version is a modification to the sequential version of CFD
 * which was written by Saleh Elmohamed. This modified version uses
 * mpiJava and was tested on upto 4 processors. However, it was designd
 * to run any number of processors. 
 * Initially, processor 0 (acting as the master) reads in the data
 * from a datafile. Then computation is done by all processors
 * invloved including processor 0. In addition , processor 0 handles
 * all work invloved in the display part. In a nutshell, processor 0
 * takes care of the I/O as well as doing its share of computation.
*/
/**********************************************************************/
/* CFD.java                                                           */
/* A 2-D Inviscid Flow Simulation.                                    */
/* Saleh Elmohamed.            (1997)                                 */ 
/**********************************************************************/
/* a note to experienced object oriented programmers:                 */
/* Some parts of this code will look sloppy because I used inline     */
/* code for some things which could have been done more elegantly     */
/* using public methods.                                              */
/* This was done because inlining by hand increased the execution     */
/* speed by about 20%.                                                */
/**********************************************************************/
/*The simulation:
 *-----------
 *  The code simulates a 2-D inviscid flow through an axisymmetric 
 *  nozzle. The simulation yields countour plots of all flow variables, 
 *  including velocity components, pressure, mach number, density and 
 *  entropy, and temperature. The plots show the location of any shock
 *  wave that would reside in the nozzle. Also, the code finds the 
 *  steady state solution to the  2-D Euler equations 
 *  (see report and reference). 
 *
 *  It uses a 4-stage Runge-Kutta time-stepping algorithm and a finite
 *  volume centeral-difference technique to find the solution. 
 *  At each stage, the residual is multiplied by a different value. 
 *  Generally, this gives a bit more accurate results since you go
 *  over more than one stage.
 *  A numerical dissipation model is emplyed to dampen any spurious 
 *  oscillations and prevent the soultion from blowing up in the 
 *  presence of shock waves (Jameson et al, 81).  The dampening order 
 *  is set by:  (1) compute the even derivatives (i.e. 2nd and 4th), 
 *  (2) set the order proportional to this.
 *
 *  If the inlet and exit conditions are supersonic, the code requires
 *  about 1800 iterations to complete execution. However, if subsonic 
 *  boundary conditions are employed, the number of iterations 
 *  required jumps to about 6000, as new errors are introduced into
 *  the domain by these conditions (see Dang, 1997). 
 *  The flow is going from left to right across the nozzle. To see the
 *  flow through, a "tuft" key is added as well as a "grid" key. 
 *  The flow can be either "subsonic", "transonic", or "supersonic".
 *  The subsonic flow requires < 1, the transonic flow is partly sub 
 *  and partly super, and the supersonic requires a flow of greater 
 *  than and equal 1.
 *
 * References:
 ************
 *  1. "Numerical Solutions of the Euler Equations by Finite Volume 
 *     Methods Using Runge-Kutta Time-Stepping Schemes", A. Jameson,
 *     W. Schmidt, and E. Turkel, paper no: AIAA-81-1259, AIAA 14th
 *     Fluid and Plasma Dynamics Conference, June 23-35, 1981, 
 *     Palo Alto, Calif.
 *
 *  2. "Monograph in CFD", T. Q. Dang, Dept. of Mechnical, Aerospace 
 *     and Manufacturing Engineering. Syracuse Univ. 1997.
 *      
 ********************************************************************/
import java.awt.*;
import java.awt.event.*;
import java.io.InputStream;         /* Lets you read files */
import java.io.StreamTokenizer;     /* Lets you interprit what you read */
import java.net.*;                  /* Let you open a URL */
import java.io.IOException;         /* Deals with io exceptions */
import java.io.*;
import java.util.*;
import mpi.*;

public class CFD extends Frame implements Runnable,
                                          ItemListener,
                                          ActionListener{

          /* Declare "globals" */
  Checkbox drawGridBox, pauseBox, tuftsBox;
  Choice plotChoice, timeChoice, viewChoice;
  Button quit;
  Color ctable[];
  DataPanel dataPanel;
  DampingPanel sdampingPanel, fdampingPanel;
  MachffPanel machffPanel;
  String inname = null;
  StringBuffer buffer = new StringBuffer();
  Thread drawThread = null;

  double machffmin;
  double machffmax; 
  final double dampingMin = 0.1;
  final double dampingMax = 5.0;
  int me;
  int[] terminate = new int[1];
  double sTime,eTime;   

  public CFD() {
    me = MPI.COMM_WORLD.Rank();
    InputStream instream = null;
    Panel controlPanelSouth = new Panel();
    Panel controlPanelNorth = new Panel();
    Panel controlPanelEast = new Panel();

    inname = "Tran.dat";
    try {
        instream = new DataInputStream(new FileInputStream(inname));
        } catch (IOException e){
            System.out.println("File not opened properly\n"+e.toString() );
            System.exit(1);
          }  

    dataPanel = new DataPanel(instream);

    /* Make the control panel */
    controlPanelNorth.setLayout(new FlowLayout(FlowLayout.LEFT, 5, 5));
    controlPanelSouth.setLayout(new FlowLayout(FlowLayout.LEFT, 5, 5));
    controlPanelEast.setLayout(new FlowLayout(FlowLayout.LEFT, 5, 5));

    controlPanelNorth.setBackground(Color.lightGray);
    controlPanelSouth.setBackground(Color.lightGray);
    controlPanelEast.setBackground(Color.lightGray);

    quit = new Button("QUIT");
    quit.addActionListener(this);
    
    plotChoice = new Choice();
    plotChoice.addItem("Contours");
    plotChoice.addItem("Line Plot");

    plotChoice.addItemListener(this);
    
    viewChoice = new Choice();
    viewChoice.addItem("Density");
    viewChoice.addItem("Pressure");
    viewChoice.addItem("Temperature");
    viewChoice.addItem("Mach Number");
    viewChoice.addItem("Axial Velocity");
    viewChoice.addItem("Normal Velocity");
    viewChoice.addItem("Total Pressure");

    viewChoice.addItemListener(this);

    drawGridBox = new Checkbox("Grid", null, false); 
    tuftsBox    = new Checkbox("Tufts", null, false);

    timeChoice  = new Choice();
    timeChoice.addItem("Local Timestep");
    timeChoice.addItem("Time Accurate");

    timeChoice.addItemListener(this);
    
    pauseBox = new Checkbox("Pause", null, false);
     
    pauseBox.setFont(new Font("Helvetica", Font.PLAIN, 15));
    drawGridBox.setFont(new Font("Helvetica", Font.PLAIN, 15));
    tuftsBox.setFont(new Font("Helvetica", Font.PLAIN, 15));
 
    drawGridBox.addItemListener(this);
    tuftsBox.addItemListener(this);
    pauseBox.addItemListener(this);
    /* Set limits on mach slider.          */
    /* These shifting limits are necessary */
    /* because there is a singularity      */
    /* in the boundary conditions when the */
    /* mach number == 1.0                  */
    if (dataPanel.machff < 0.9) {
      machffmin = 0.2;
      machffmax = 0.9;
    }
    else if (dataPanel.machff > 1.3) {
      machffmin = 1.3;
       machffmax = 2.5;
    }
    else {
      machffmin = 0.2;
      machffmax = 0.9;
    }

    machffPanel = new MachffPanel("Inlet Mach Number", machffmin, machffmax,
                                  dataPanel.machff);
    sdampingPanel = new DampingPanel("2nd Order Damping",
                                     dampingMin, dampingMax,
                                     dataPanel.secondOrderDamping);
    fdampingPanel = new DampingPanel("4th Order Damping", dampingMin,
                                     dampingMax,
                                     dataPanel.fourthOrderDamping);
    
    /* Add the button panel */
    controlPanelSouth.setFont(new Font("Helvetica", Font.PLAIN, 15));
    dataPanel.setFont(new Font("Helvetica", Font.PLAIN, 15));
    Label label2 = new Label (" of");
    controlPanelNorth.add(plotChoice);
    controlPanelNorth.add(label2);
    controlPanelNorth.add(viewChoice);
    controlPanelNorth.add(drawGridBox);
    controlPanelNorth.add(tuftsBox);
    controlPanelNorth.add(pauseBox);
    controlPanelNorth.add(timeChoice);
    controlPanelSouth.add(machffPanel);
    controlPanelSouth.add(fdampingPanel);
    controlPanelSouth.add(sdampingPanel);
    controlPanelEast.add(quit);

    add("East", controlPanelEast);
    add("South", controlPanelSouth);
    add("North", controlPanelNorth);
    /* Add contour plot to main screen */
    add("Center", dataPanel);
//System.out.println("2 : "+controlPanelSouth.getBackground());
  }

  public void start() {
    terminate[0] = 1;
    if (drawThread == null) {
      drawThread = new Thread(this);
      /* Priority has to be turned down so numerics  */
      /* don't hose graphics */
      drawThread.setPriority(Thread.MIN_PRIORITY+2);
      drawThread.start();
    }
  }

  public void run(){
    while (drawThread != null) {
      MPI.COMM_WORLD.Bcast(terminate,0,1,MPI.INT,0);
      if(terminate[0] == 0){
          MPI.Finalize();
        System.exit(0);
      }
      if (me == 0){
        dataPanel.machff = machffPanel.sliderValue;
        dataPanel.secondOrderDamping = sdampingPanel.sliderValue;
        dataPanel.fourthOrderDamping = fdampingPanel.sliderValue;
      }
      dataPanel.doIteration();
      if (me == 0)
        try {
          drawThread.sleep(1000);
        } catch (InterruptedException e) {
          System.err.println("Execution interrupted by another thread.");
          System.err.println(e.getMessage());
        }
    }
  }
  
  public void stop() {
    if (me == 0){
      if (drawThread != null && drawThread.isAlive()) {
               drawThread.stop();
      }
      drawThread = null;
    }
  }

  public void actionPerformed(ActionEvent evt){
    String arg = evt.getActionCommand();
    if (arg.equals("QUIT")) terminate[0] = 0;
 //   System.exit(0);
  }
  public void itemStateChanged(ItemEvent e) { 
    if (e.getItemSelectable() instanceof Choice){
      dataPanel.nplot = plotChoice.getSelectedIndex();
      dataPanel.nview = viewChoice.getSelectedIndex();
      dataPanel.ntime = timeChoice.getSelectedIndex();
    }
    else if (e.getItemSelectable() instanceof Checkbox){
      dataPanel.setGridSwitch(drawGridBox.getState());
      dataPanel.setTuftsSwitch(tuftsBox.getState());
      if(pauseBox.getState() == true) {
        drawThread.suspend();
      }
      else
        drawThread.resume();
    }    
  }
  public static void main(String[] args)
  {
    MPI.Init(args);
    int root = MPI.COMM_WORLD.Rank();
    
    CFD f = new CFD();

    f.setSize(1024,850);
    if(root == 0)
      f.show();

    f.start();
  }
}

/* The fundamental class is the data          */
/* This executes when you create a new object */
class datafield {
  private double a[][];               /* Grid cell area */
  private double local_a[][];
  int n_bar; 

  double deltat[][];                  /* Timestep */
  double local_deltat[][];

  private double machff;              /* Inflow mach number */
  public double machfftarget;         /* Target inflow mach number */
  public double secondOrderDamping;
  public double fourthOrderDamping;
  public int ntime;

  private double opg[][], pg[][], pg1[][];   /* Pressure */
  private double local_opg[][], local_pg[][], local_pg1[][];

  private double sxi[][], seta[][];
  private double local_sxi[][], local_seta[][];

  private double tg[][], tg1[][];           /* Temperature */
  private double local_tg[][], local_tg1[][];

  private double xnode[][], ynode[][];       /* Storage of node coordinates */
  private double local_xnode[][], local_ynode[][];

  double xmin=Double.MAX_VALUE, xmax=Double.MIN_VALUE,
         local_xmin=Double.MAX_VALUE, local_xmax=Double.MIN_VALUE,
         ymin=Double.MAX_VALUE, ymax=Double.MIN_VALUE,
         local_ymin=Double.MAX_VALUE, local_ymax=Double.MIN_VALUE;
                                    /* Max and Min of node coordinates */
  double hold_MinMax[];
  double cff, uff, vff, pff, rhoff, tff, jplusff, jminusff;
                                    /* Far field values */
  double datamax, datamin;
  int iter,my_rank,np;
  int imax, jmax;                   /* Number of nodes in x and y direction */
  private Statevector d[][];        /* Damping coefficients */
  private Statevector local_d[][];

  private Statevector f[][], g[][]; /* Flux Vectors */
  private Statevector local_f[][], local_g[][];

  private Statevector r[][], ug1[][];
  private Statevector local_r[][];
  private Statevector local_ug1[][];
  private Statevector ug[][];      /* Storage of data */
  private Statevector local_ug[][];

  final double Cp = 1004.5;         /* specific heat, const pres. */
  final double Cv=717.5;            /* specific heat, const vol. */
  final double gamma=1.4;           /* Ratio of specific heats */
  final double rgas=287.0;          /* Gas Constant */
  final int keytextwidth=50;        /* Width of text field in color table */
  final int ncontour = 16;          /* Number of contours between min and max */
  final double fourthOrderNormalizer = 0.02; /* Damping coefficients */
  final double secondOrderNormalizer = 0.02;

  datafield(InputStream instream) throws IOException {
    int i, j, n;                    /* Dummy counters */
    double scrap, scrap2;           /* Temporary storage */

                                    /* Convert the stream into 
                                       tokens (which helps you parse it) 
                                     *-----------------*/
      Reader reader = new BufferedReader(new InputStreamReader(instream));
      StreamTokenizer intokens = new StreamTokenizer(reader);

    my_rank = MPI.COMM_WORLD.Rank();
    np = MPI.COMM_WORLD.Size();
    if(my_rank == 0){
      /* Read header */
      if (intokens.nextToken() == StreamTokenizer.TT_NUMBER)
        imax = (int) intokens.nval;
      else
        throw new IOException();

      if (intokens.nextToken() == StreamTokenizer.TT_NUMBER)
        jmax = (int) intokens.nval;
      else
        throw new IOException();

      if (intokens.nextToken() == StreamTokenizer.TT_NUMBER) {
        machfftarget = intokens.nval;
        machff = machfftarget;
      }
      else
        throw new IOException();

      if (intokens.nextToken() == StreamTokenizer.TT_NUMBER)
        iter = (int) intokens.nval;
      else
        throw new IOException();
    }
    
    double basic_ele[] = new double[5];
    if(my_rank == 0){
      basic_ele[0] = imax;
      basic_ele[1] = jmax;
      basic_ele[2] = machfftarget;
      basic_ele[3] = machff;
      basic_ele[4] = iter;
    }
   /* Send initial value to other processors */
    MPI.COMM_WORLD.Bcast(basic_ele,0,5,MPI.DOUBLE,0);

    imax = (int)basic_ele[0];
    jmax = (int)basic_ele[1];
    machfftarget = basic_ele[2];
    machff = basic_ele[3];
    iter = (int)basic_ele[4];

    n_bar = (imax -1)/np;
  
    deltat = new double[imax+1][jmax+2];
    local_deltat = new double[n_bar+2][jmax+2];

    opg = new double[imax+2][jmax+2];
    local_opg = new double[n_bar+2][jmax+2];

    pg = new double[imax+2][jmax+2];
    local_pg = new double[n_bar+2][jmax+2];

    pg1 = new double[imax+2][jmax+2];
    local_pg1 = new double[n_bar+2][jmax+2];

    sxi =new double[imax+2][jmax+2];;
    local_sxi = new double[n_bar+2][jmax+2];

    seta = new double[imax+2][jmax+2];;
    local_seta = new double[n_bar+2][jmax+2];

    tg = new double[imax+2][jmax+2];
    local_tg = new double[n_bar+2][jmax+2];

    tg1 = new double[imax+2][jmax+2];
    local_tg1 = new double[n_bar+2][jmax+2];

    ug = new Statevector[imax+2][jmax+2];
    local_ug = new Statevector[n_bar+2][jmax+2];

    a = new double[imax][jmax];
    local_a = new double[n_bar+2][jmax];

    d =  new Statevector[imax+2][jmax+2];
    local_d = new Statevector[n_bar+2][jmax+2];

    f =  new Statevector[imax+2][jmax+2];
    local_f = new Statevector[n_bar+2][jmax+2];

    g =  new Statevector[imax+2][jmax+2];
    local_g = new Statevector[n_bar+2][jmax+2];

    r =  new Statevector[imax+2][jmax+2];
    local_r = new Statevector[n_bar+2][jmax+2];

    ug1 =  new Statevector[imax+2][jmax+2];
    local_ug1 = new Statevector[n_bar+2][jmax+2];

    xnode = new double[imax][jmax];
    local_xnode = new double[n_bar+2][jmax];

    ynode = new double[imax][jmax];
    local_ynode = new double[n_bar+2][jmax];
    
    for (i = 0; i < n_bar+2 ; i++)
       for (j = 0; j < jmax+2 ; j++){
          local_ug1[i][j] = new Statevector();
          local_ug[i][j] = new Statevector();
          local_d[i][j] = new Statevector();
          local_f[i][j] = new Statevector();
          local_g[i][j] = new Statevector();
          local_r[i][j] = new Statevector();
       }

    int initUg = 0;
    for (i = 0; i < imax+2; ++i)
      for (j = 0; j < jmax+2; ++j) {
        d[i][j] =  new Statevector();
        f[i][j] =  new Statevector();
        g[i][j] =  new Statevector();
        r[i][j] = new Statevector();
        ug[i][j] = new Statevector();
        ug1[i][j] = new Statevector();
        initUg ++; 
      }

    /* Set farfield values (we use normalized units for everything */    
    cff = 1.0;
    vff = 0.0;
    pff = 1.0 / gamma;
    rhoff = 1.0;
    tff = pff / (rhoff * rgas);

    /* Now load the data */
    hold_MinMax = new double[4];
    if(my_rank == 0){
      for (i = 1; i < imax+1; i++)
        for (j = 1; j < jmax+1; ++j) {
          switch(intokens.nextToken()) {
          case StreamTokenizer.TT_NUMBER:
            xnode[i-1][j-1] = (double) intokens.nval;
            xmax = Math.max(xmax, xnode[i-1][j-1]);
            xmin = Math.min(xmin, xnode[i-1][j-1]);
            if (intokens.nextToken() == StreamTokenizer.TT_NUMBER) {
              ynode[i-1][j-1] = (double) intokens.nval;
              ymax = Math.max(ymax, ynode[i-1][j-1]);
              ymin = Math.min(ymin, ynode[i-1][j-1]);
              if(intokens.nextToken() == StreamTokenizer.TT_NUMBER) {
                ug[i][j].a =(double) intokens.nval;
                if(intokens.nextToken() == StreamTokenizer.TT_NUMBER) {
                  ug[i][j].b =(double) intokens.nval;
                  if(intokens.nextToken() == StreamTokenizer.TT_NUMBER) {
                    ug[i][j].c =(double) intokens.nval;
                    if(intokens.nextToken() == StreamTokenizer.TT_NUMBER) {
                      ug[i][j].d =(double) intokens.nval;
                      scrap = ug[i][j].c/ug[i][j].a;
                      scrap2 = ug[i][j].b/ug[i][j].a;
                      tg[i][j] = ug[i][j].d/ug[i][j].a 
                        - (0.5 * (scrap*scrap + scrap2*scrap2));
                      tg[i][j] = tg[i][j] / Cv;
                      pg[i][j] = rgas * ug[i][j].a * tg[i][j];
                    }
                  }
                }
              }
            }
            break;
          case StreamTokenizer.TT_EOF:
            throw new IOException();
          default:
            throw new IOException();
          }
        }
      }

/*
------------- send boundary elements to other processors ----------------
*/ 
    int initPg = 0;
    int initXY = 0;
    double[] global_array = new double[(imax+2)*(jmax+2)*16];
    if(my_rank == 0){
      for (i = 1; i <= imax; ++i) 
        for (j = 1; j <= jmax; ++j){ 
          global_array[initPg++] = pg[i][j];
          global_array[initPg++] = tg[i][j];
          global_array[initPg++] = ug[i][j].a;
          global_array[initPg++] = ug[i][j].b;
          global_array[initPg++] = ug[i][j].c;
          global_array[initPg++] = ug[i][j].d;
        }    
      global_array[initPg++] = xmin;
      global_array[initPg++] = xmax;
      global_array[initPg++] = ymin;
      global_array[initPg++] = ymax;

      for (i=0; i < imax;i++)
        for(j=0;j < jmax;j++){
          global_array[initPg++] = xnode[i][j];
          global_array[initPg++] = ynode[i][j];
        }
    }

    MPI.COMM_WORLD.Bcast(global_array,0,(imax+2)*(jmax+2)*16,MPI.DOUBLE,0);

    if (my_rank != 0){
      for (i = 1; i <= imax; ++i) 
        for (j = 1; j <= jmax; ++j){ 
          pg[i][j] = global_array[initPg++];  
          tg[i][j] = global_array[initPg++];
          ug[i][j].a = global_array[initPg++];
          ug[i][j].b = global_array[initPg++];
          ug[i][j].c = global_array[initPg++];
          ug[i][j].d = global_array[initPg++];
        }
      xmin = global_array[initPg++];
      xmax = global_array[initPg++];
      ymin = global_array[initPg++];
      ymax = global_array[initPg++];

      for (i=0; i < imax;i++)
        for(j=0;j < jmax;j++){
          xnode[i][j] = global_array[initPg++];
          ynode[i][j] = global_array[initPg++];
        }
    }

    int global = 0;
    for (i = 1; i <= n_bar; ++i) 
      for (j = 1; j <= jmax; ++j){
        global = my_rank * n_bar + i;
        local_pg[i][j] = pg[global][j];
        local_tg[i][j] = tg[global][j];
        local_ug[i][j].a = ug[global][j].a;
        local_ug[i][j].b = ug[global][j].b;
        local_ug[i][j].c = ug[global][j].c;
        local_ug[i][j].d = ug[global][j].d;
      }

    int lastPro = np - 1;
    int position = 0;
    if (my_rank == lastPro)
      position = imax;
    else 
      position = (n_bar+1)*(my_rank+1);

    int pgPosition = (n_bar+1)*my_rank-1;

    for (j = 1; j<=jmax;++j){
      local_pg[n_bar+1][j]   = pg[position][j];
      local_tg[n_bar+1][j]   = tg[position][j];
      local_ug[n_bar+1][j].a = ug[position][j].a;
      local_ug[n_bar+1][j].b = ug[position][j].b;
      local_ug[n_bar+1][j].c = ug[position][j].c;
      local_ug[n_bar+1][j].d = ug[position][j].d;
      if(my_rank != 0)
        local_pg[0][j] = pg[pgPosition][j]; 
    }

    global = 0;
    for (i = 0; i < n_bar; ++i) 
      for (j = 0; j < jmax; ++j){
        global = my_rank * n_bar + i;
        local_xnode[i][j] = xnode[global][j];
        local_ynode[i][j] = ynode[global][j];
      }

    if (my_rank == lastPro)
       position = imax - 1;
    else 
       position = n_bar * (my_rank + 1);

    for (i = 0;i<jmax ; i++){
      local_xnode[n_bar][i] = xnode[position][i];
      local_ynode[n_bar][i] = ynode[position][i];
    }    

    /* Calculate grid cell areas */

    for (i = 1; i <= n_bar; ++i)
      for (j = 1; j < jmax; ++j){
        local_a[i][j] = 0.5 * ((local_xnode[i][j] - local_xnode[i-1][j-1]) 
                            * (local_ynode[i-1][j] - local_ynode[i][j-1])-
                              (local_ynode[i][j] - local_ynode[i-1][j-1]) 
                            * (local_xnode[i-1][j] - local_xnode[i][j-1]));
      }
    double send_a[] = new double[jmax];
    double recv_a[] = new double[jmax];
    int right = my_rank + 1;
    int left  = my_rank - 1;
    if (my_rank == 0)
      left = np -1;
    else if(my_rank == (np -1))
      right = 0;
       
    for(j=1;j<jmax;j++)
      send_a[j] = local_a[n_bar][j];

    MPI.COMM_WORLD.Sendrecv(send_a,0,jmax,MPI.DOUBLE,right,1,
                            recv_a,0,jmax,MPI.DOUBLE,left ,1);
    if(my_rank != 0)
      for(j=1;j<jmax;j++)
        local_a[0][j] = recv_a[j];

    for(j=1;j<jmax;j++)
      send_a[j] = local_a[1][j];

    MPI.COMM_WORLD.Sendrecv(send_a,0,jmax,MPI.DOUBLE,left ,1,
                            recv_a,0,jmax,MPI.DOUBLE,right,1);
    if (my_rank != (np-1))
      for(j=1;j<jmax;j++)
        local_a[n_bar+1][j] = recv_a[j];  
  }

  public void drawKey(Color ctable[], Graphics g, int xkeywidth, 
                      int ykeywidth, int xcanvaswidth, int ycanvaswidth) {
    int colorbandheight, colorbandwidth;
    int n;
    int yband;
    int xoffset, yoffset;
    FontMetrics keyFontMetrics;

    xoffset = xcanvaswidth - xkeywidth;
    yoffset = ycanvaswidth - ykeywidth;
    colorbandwidth = xcanvaswidth - keytextwidth - (xoffset+5)-20;
    colorbandheight = (ykeywidth-10) / (ctable.length-1) + 1;

    g.setColor(Color.black);
    g.fillRect(xoffset, yoffset, xkeywidth+20, ykeywidth);

    /* Draw the colors themselves */
    for (n = 0; n < ctable.length; ++n) {
      g.setColor(ctable[n]);
      yband =  ycanvaswidth - 5 -
        (int) (((double)(n+1)/(double) ctable.length) 
               * (double)(ykeywidth-10));
      g.fillRect(xoffset+5, yband, colorbandwidth, colorbandheight);
    }

    /* Write in the text for the key */
    keyFontMetrics = g.getFontMetrics();
    keyFontMetrics.getHeight();
    g.setColor(Color.white);
    StringBuffer writeBuffer = new StringBuffer().append(datamin);
    writeBuffer.setLength(6);
    g.drawString(" "+writeBuffer, xoffset+5+colorbandwidth, ycanvaswidth - 5);
    StringBuffer writeBuffer2 = new StringBuffer().append(datamax);
    writeBuffer2.setLength(6);
    g.drawString(" "+writeBuffer2, xoffset+5+colorbandwidth, 
                 keyFontMetrics.getHeight() + 5);

    /* Draw outline frame */
    g.setColor(Color.white);
    g.drawRect(xoffset, yoffset, xkeywidth, ykeywidth);
  }

  void doIteration() {
    double error;
    double scrap;
    int i, j;
    /* Record the old pressure values */
    for (i = 1; i < imax; ++i)
      for (j = 1; j < jmax; ++j) {
        opg[i][j] = pg[i][j];
      }

    /* Set new far field mach number (if necessary) */
    if (machff != machfftarget) {
      if (Math.abs(machfftarget-machff) <= 0.02)
        machff = machfftarget;
      else {
        if (machfftarget > machff) {
          machff += 0.02;
        }
        else
          machff -= 0.02;
      }
    }

    calculateDummyCells(local_pg, local_tg, local_ug);


    calculateDeltaT();
    calculateDamping(local_pg, local_ug);

    /* Do the integration */
    /*...... Step 1 ......*/
    calculateF(local_pg, local_tg, local_ug);
    calculateG(local_pg, local_tg, local_ug);
    calculateR();

    for (i = 1; i <= n_bar; ++i)
     for(j = 1; j < jmax; ++j){
        local_ug1[i][j].a=
           local_ug[i][j].a-0.25*local_deltat[i][j]/local_a[i][j]
           *(local_r[i][j].a-local_d[i][j].a);

        local_ug1[i][j].b=
           local_ug[i][j].b-0.25*local_deltat[i][j]/local_a[i][j]
           *(local_r[i][j].b-local_d[i][j].b);

        local_ug1[i][j].c=
           local_ug[i][j].c-0.25*local_deltat[i][j]/local_a[i][j]
           *(local_r[i][j].c-local_d[i][j].c);

        local_ug1[i][j].d=
           local_ug[i][j].d-0.25*local_deltat[i][j]/local_a[i][j]
           *(local_r[i][j].d-local_d[i][j].d);
    }

    calculateStateVar(local_pg1, local_tg1, local_ug1);

    /*...... Step 2 ......*/
    calculateDummyCells(local_pg1,local_tg1, local_ug1);
    calculateF(local_pg1, local_tg1, local_ug1);
    calculateG(local_pg1, local_tg1, local_ug1);
    calculateR();

    for (i = 1; i <= n_bar; ++i)
      for (j = 1; j < jmax; ++j) {
        local_ug1[i][j].a=
          local_ug[i][j].a-0.33333*local_deltat[i][j]/local_a[i][j]
          *(local_r[i][j].a-local_d[i][j].a);

        local_ug1[i][j].b=
          local_ug[i][j].b-0.33333*local_deltat[i][j]/local_a[i][j]
          *(local_r[i][j].b-local_d[i][j].b);

        local_ug1[i][j].c=
          local_ug[i][j].c-0.33333*local_deltat[i][j]/local_a[i][j]
          *(local_r[i][j].c-local_d[i][j].c);

        local_ug1[i][j].d=
          local_ug[i][j].d-0.33333*local_deltat[i][j]/local_a[i][j]
          *(local_r[i][j].d-local_d[i][j].d);
      }

    calculateStateVar(local_pg1, local_tg1, local_ug1);
    
    /*...... Step 3 ......*/
    calculateDummyCells(local_pg1, local_tg1, local_ug1);
    calculateF(local_pg1, local_tg1, local_ug1);
    calculateG(local_pg1, local_tg1, local_ug1);
    calculateR();
    for (i = 1; i <= n_bar; ++i)
      for (j = 1; j < jmax; ++j) {
        local_ug1[i][j].a=
          local_ug[i][j].a-0.5*local_deltat[i][j]/local_a[i][j]
          *(local_r[i][j].a-local_d[i][j].a);

        local_ug1[i][j].b=
          local_ug[i][j].b-0.5*local_deltat[i][j]/local_a[i][j]
          *(local_r[i][j].b-local_d[i][j].b);
        local_ug1[i][j].c=
          local_ug[i][j].c-0.5*local_deltat[i][j]/local_a[i][j]
          *(local_r[i][j].c-local_d[i][j].c);
        local_ug1[i][j].d=
          local_ug[i][j].d-0.5*local_deltat[i][j]/local_a[i][j]
          *(local_r[i][j].d-local_d[i][j].d);
      }
    calculateStateVar(local_pg1, local_tg1, local_ug1);
    
    /*..... Step 4 (final step) .....*/
    calculateDummyCells(local_pg1, local_tg1, local_ug1);
    calculateF(local_pg1, local_tg1, local_ug1);
    calculateG(local_pg1, local_tg1,local_ug1);
    calculateR(); 

    int initUg = 0;
    for (i = 1; i <= n_bar; ++i)
      for (j = 1; j < jmax; ++j) {
        local_ug[i][j].a -= local_deltat[i][j]/local_a[i][j]
                            *(local_r[i][j].a-local_d[i][j].a);

        local_ug[i][j].b -= local_deltat[i][j]/local_a[i][j]
                            *(local_r[i][j].b-local_d[i][j].b);

        local_ug[i][j].c -= local_deltat[i][j]/local_a[i][j]
                            *(local_r[i][j].c-local_d[i][j].c);

        local_ug[i][j].d -= local_deltat[i][j]/local_a[i][j]
                            *(local_r[i][j].d-local_d[i][j].d);
     } 

    calculateStateVar(local_pg, local_tg, local_ug);

    double one_dim_array[] = new double[(n_bar+1) * jmax * 6];
    int initDim = 0;
    for (i = 1; i <= n_bar+1; ++i) 
       for (j = 1; j <= jmax; ++j) {
         one_dim_array[initDim++] = local_ug[i][j].a;
         one_dim_array[initDim++] = local_ug[i][j].b;
         one_dim_array[initDim++] = local_ug[i][j].c;
         one_dim_array[initDim++] = local_ug[i][j].d;
         one_dim_array[initDim++] = local_pg[i][j];
         one_dim_array[initDim++] = local_tg[i][j];
       }
    double global_dim_array[] = new double[(imax +2)*(jmax + 2)*6];
    MPI.COMM_WORLD.Gather(one_dim_array,0,initDim,MPI.DOUBLE,
             global_dim_array,0,initDim,MPI.DOUBLE,0);

    MPI.COMM_WORLD.Bcast(global_dim_array,0,(imax+2)*(jmax+2)*6,MPI.DOUBLE,0);

    initDim = 0;
    int k = 1;
    int boundary = n_bar + 1;
    for (i = 1; i <= (n_bar+1)*np;++i){
      for (j = 1; j <= jmax; ++j)
        if (i != boundary){
           ug[k][j].a = global_dim_array[initDim++];
           ug[k][j].b = global_dim_array[initDim++];
           ug[k][j].c = global_dim_array[initDim++];
           ug[k][j].d = global_dim_array[initDim++];
           pg[k][j] = global_dim_array[initDim++];
           tg[k][j] = global_dim_array[initDim++];
        }
      if (i == boundary){
         initDim += jmax*6;
         k -= 1;
         boundary += boundary; 
      } 
      k++;
    }

    /* calculate RMS Pressure Error. */
    error = 0.0;
    for (i = 1; i < imax; ++i)
      for (j = 1; j < jmax; ++j) {
        scrap = pg[i][j] - opg[i][j];
        error += scrap*scrap;
      }
    error = Math.sqrt(error / (double)((imax-1) * (jmax-1)) );
//System.out.println("error : "+error); 
  }


  private void calculateStateVar(double localpg[][], double localtg[][],
                         Statevector localug[][])
    /* Calculates the new state values for Range-Kutta */
    /* Works for default values, 8/15 at 9:45 pm */
    {
      double temp, temp2;
      int i, j;

      for (i = 1; i <= n_bar; ++i) {
        for (j = 1; j < jmax; ++j) {
          temp = localug[i][j].b;
          temp2 = localug[i][j].c;
          localtg[i][j] = localug[i][j].d/localug[i][j].a - 0.5 *
            (temp*temp + temp2*temp2)/(localug[i][j].a*localug[i][j].a);

          localtg[i][j] = localtg[i][j] / Cv;
          localpg[i][j] = localug[i][j].a * rgas * localtg[i][j];
        }
      }
    }

  private void calculateR() {
    /* Works for default values, straight channel (all 0's) 8/11, 9:15 pm */
    
    double deltax, deltay;
    double temp;
    int i,j;
    Statevector scrap;
    
    for (i = 1; i <= n_bar; ++i) {
      for (j = 1; j < jmax; ++j) {
        
        /* Start by clearing R */
        local_r[i][j].a = 0.0;
        local_r[i][j].b = 0.0;
        local_r[i][j].c = 0.0;
        local_r[i][j].d = 0.0;
        
        /* East Face */
        deltay = (local_ynode[i][j] - local_ynode[i][j-1]);
        deltax = (local_xnode[i][j] - local_xnode[i][j-1]);
        temp = 0.5 * deltay;
        local_r[i][j].a += temp*(local_f[i][j].a + local_f[i+1][j].a);
        local_r[i][j].b += temp*(local_f[i][j].b + local_f[i+1][j].b);
        local_r[i][j].c += temp*(local_f[i][j].c + local_f[i+1][j].c);
        local_r[i][j].d += temp*(local_f[i][j].d + local_f[i+1][j].d);

        temp = -0.5*deltax;
        local_r[i][j].a += temp * (local_g[i][j].a+local_g[i+1][j].a);
        local_r[i][j].b += temp * (local_g[i][j].b+local_g[i+1][j].b);
        local_r[i][j].c += temp * (local_g[i][j].c+local_g[i+1][j].c);
        local_r[i][j].d += temp * (local_g[i][j].d+local_g[i+1][j].d);
        
        /* South Face */
        deltay = (local_ynode[i][j-1] - local_ynode[i-1][j-1]);  
        deltax = (local_xnode[i][j-1] - local_xnode[i-1][j-1]);

        temp = 0.5 * deltay;
        local_r[i][j].a  += temp*(local_f[i][j].a+local_f[i][j-1].a);
        local_r[i][j].b  += temp*(local_f[i][j].b+local_f[i][j-1].b);
        local_r[i][j].c  += temp*(local_f[i][j].c+local_f[i][j-1].c);
        local_r[i][j].d  += temp*(local_f[i][j].d+local_f[i][j-1].d);

        temp = -0.5*deltax;
        local_r[i][j].a += temp * (local_g[i][j].a+local_g[i][j-1].a);
        local_r[i][j].b += temp * (local_g[i][j].b+local_g[i][j-1].b);
        local_r[i][j].c += temp * (local_g[i][j].c+local_g[i][j-1].c);
        local_r[i][j].d += temp * (local_g[i][j].d+local_g[i][j-1].d);
        
        /* West Face */
        deltay = (local_ynode[i-1][j-1] - local_ynode[i-1][j]);
        deltax = (local_xnode[i-1][j-1] - local_xnode[i-1][j]);

        temp = 0.5 * deltay;
        local_r[i][j].a  += temp*(local_f[i][j].a+local_f[i-1][j].a);
        local_r[i][j].b  += temp*(local_f[i][j].b+local_f[i-1][j].b);
        local_r[i][j].c  += temp*(local_f[i][j].c+local_f[i-1][j].c);
        local_r[i][j].d  += temp*(local_f[i][j].d+local_f[i-1][j].d);

        temp = -0.5*deltax;
        local_r[i][j].a += temp * (local_g[i][j].a+local_g[i-1][j].a);
        local_r[i][j].b += temp * (local_g[i][j].b+local_g[i-1][j].b);
        local_r[i][j].c += temp * (local_g[i][j].c+local_g[i-1][j].c);
        local_r[i][j].d += temp * (local_g[i][j].d+local_g[i-1][j].d);
                
        /* North Face */
        deltay = (local_ynode[i-1][j] - local_ynode[i][j]);
        deltax = (local_xnode[i-1][j] - local_xnode[i][j]);

        temp = 0.5 * deltay;
        local_r[i][j].a  += temp*(local_f[i][j].a+local_f[i+1][j].a);
        local_r[i][j].b  += temp*(local_f[i][j].b+local_f[i+1][j].b);
        local_r[i][j].c  += temp*(local_f[i][j].c+local_f[i+1][j].c);
        local_r[i][j].d  += temp*(local_f[i][j].d+local_f[i+1][j].d);
        
        temp = -0.5*deltax;
        local_r[i][j].a += temp * (local_g[i][j].a+local_g[i][j+1].a);
        local_r[i][j].b += temp * (local_g[i][j].b+local_g[i][j+1].b);
        local_r[i][j].c += temp * (local_g[i][j].c+local_g[i][j+1].c);
        local_r[i][j].d += temp * (local_g[i][j].d+local_g[i][j+1].d);
      }
    }
  }

  private void calculateG(double localpg[][], double localtg[][],
                  Statevector localug[][]) {
    /* Works for default values 8/15, 5:15 pm */    
    double temp, temp2, temp3;
    double v;
    int i, j;
    
    for (i = 0; i <= n_bar + 1; ++i) {
      for (j = 0; j < jmax + 1; ++j) {
        v = localug[i][j].c / localug[i][j].a;
        local_g[i][j].a = localug[i][j].c;
        local_g[i][j].b = localug[i][j].b * v;
        local_g[i][j].c = localug[i][j].c*v + localpg[i][j];
        temp = localug[i][j].b * localug[i][j].b;
        temp2 = localug[i][j].c * localug[i][j].c;
        temp3 = localug[i][j].a * localug[i][j].a;
        local_g[i][j].d = localug[i][j].c * (Cp * localtg[i][j]+ 
                 (0.5 * (temp + temp2)/(temp3)));
      }
    }

    double[] send_FG  = new double[(jmax+1)*8];
    double[] recv_FG  = new double[(jmax+1)*8];
    int initFG;
    int right = my_rank + 1;
    int left  = my_rank - 1;
    if (my_rank == 0)
      left = np -1;
    else if(my_rank == (np -1))
      right = 0;

    initFG = 0;
    for(j=0;j<jmax+1;j++){
      send_FG[initFG++] = local_g[n_bar][j].a;
      send_FG[initFG++] = local_g[n_bar][j].b;
      send_FG[initFG++] = local_g[n_bar][j].c;
      send_FG[initFG++] = local_g[n_bar][j].d;
      send_FG[initFG++] = local_f[n_bar][j].a;
      send_FG[initFG++] = local_f[n_bar][j].b;
      send_FG[initFG++] = local_f[n_bar][j].c;
      send_FG[initFG++] = local_f[n_bar][j].d;
    }
    MPI.COMM_WORLD.Sendrecv(send_FG,0,(jmax+1)*8,MPI.DOUBLE,right,1,
                            recv_FG,0,(jmax+1)*8,MPI.DOUBLE,left ,1);
    initFG = 0;
    if (my_rank != 0)
      for(j=0;j<jmax+1;j++){
        local_g[0][j].a = recv_FG[initFG++];
        local_g[0][j].b = recv_FG[initFG++];
        local_g[0][j].c = recv_FG[initFG++];
        local_g[0][j].d = recv_FG[initFG++];
        local_f[0][j].a = recv_FG[initFG++];
        local_f[0][j].b = recv_FG[initFG++];
        local_f[0][j].c = recv_FG[initFG++];
        local_f[0][j].d = recv_FG[initFG++];
      }
    initFG = 0;
    for(j=0;j<jmax+1;j++){
      send_FG[initFG++] = local_g[1][j].a;
      send_FG[initFG++] = local_g[1][j].b;
      send_FG[initFG++] = local_g[1][j].c;
      send_FG[initFG++] = local_g[1][j].d;
      send_FG[initFG++] = local_f[1][j].a;
      send_FG[initFG++] = local_f[1][j].b;
      send_FG[initFG++] = local_f[1][j].c;
      send_FG[initFG++] = local_f[1][j].d;
     }     
    MPI.COMM_WORLD.Sendrecv(send_FG,0,(jmax+1)*8,MPI.DOUBLE,left ,1,
                            recv_FG,0,(jmax+1)*8,MPI.DOUBLE,right,1);
    initFG = 0;
    if (my_rank != (np-1))
      for(j=0;j<jmax+1;j++){
         local_g[n_bar+1][j].a = recv_FG[initFG++];
         local_g[n_bar+1][j].b = recv_FG[initFG++];
         local_g[n_bar+1][j].c = recv_FG[initFG++];
         local_g[n_bar+1][j].d = recv_FG[initFG++];
         local_f[n_bar+1][j].a = recv_FG[initFG++];
         local_f[n_bar+1][j].b = recv_FG[initFG++];
         local_f[n_bar+1][j].c = recv_FG[initFG++];
         local_f[n_bar+1][j].d = recv_FG[initFG++];
       }
  }

  private void calculateF(double localpg[][], double localtg[][], 
                  Statevector localug[][]) {
     /* Works for default values 8/15, 4:50 pm */
    {
      double u;
      double temp1, temp2, temp3;
      int i, j;
      
      for (i = 0; i <= n_bar+1; ++i) {
        for (j = 0; j < jmax + 1; ++j) {          
          u = localug[i][j].b/ localug[i][j].a;
          local_f[i][j].a = localug[i][j].b;
          local_f[i][j].b = localug[i][j].b *u + localpg[i][j];
          local_f[i][j].c = localug[i][j].c * u;
          temp1 = localug[i][j].b * localug[i][j].b;
          temp2 = localug[i][j].c * localug[i][j].c;
          temp3 = localug[i][j].a * localug[i][j].a;
          local_f[i][j].d = localug[i][j].b * (Cp * localtg[i][j] + 
                 (0.5 * (temp1 + temp2)/(temp3)));
        }
      }
    }
  }

  private void calculateDamping(double localpg[][], Statevector localug[][]) {
      double adt, sbar;
      double nu2;
      double nu4;
      double tempdouble;
      int ascrap, i, j;
      Statevector temp = new Statevector();
      Statevector temp2 = new Statevector();
      Statevector scrap2 = new Statevector(), scrap4 = new Statevector();
      
      nu2 = secondOrderDamping * secondOrderNormalizer;
      nu4 = fourthOrderDamping * fourthOrderNormalizer;

      double send_pg[] = new double[jmax+1];
      double recv_pg[] = new double[jmax+1];
      int right = my_rank + 1;
      int left  = my_rank - 1;
      if (my_rank == 0)
        left = np -1;
      else if(my_rank == (np -1))
        right = 0;

      for(j=0;j<=jmax;j++)
           send_pg[j] = localpg[n_bar][j];

      MPI.COMM_WORLD.Sendrecv(send_pg,0,jmax+1,MPI.DOUBLE,right,1,
                              recv_pg,0,jmax+1,MPI.DOUBLE,left ,1);
      if(my_rank != 0)
        for(j=0;j<=jmax;j++)
          localpg[0][j] = recv_pg[j];

      for(j=0;j<=jmax;j++)
        send_pg[j] = localpg[1][j];

      MPI.COMM_WORLD.Sendrecv(send_pg,0,jmax+1,MPI.DOUBLE,left ,1,
                              recv_pg,0,jmax+1,MPI.DOUBLE,right,1);
      if (my_rank != (np-1))
        for(j=0;j<=jmax;j++)
          localpg[n_bar+1][j] = recv_pg[j];

      /* First do the pressure switches */
      /* Checked and works with defaults. */
      for (i = 1; i <= n_bar; ++i)
        for (j = 1; j < jmax; ++j){
          local_sxi[i][j] = Math.abs(localpg[i+1][j] -
                2.0 * localpg[i][j] + localpg[i-1][j])/ localpg[i][j];
          local_seta[i][j] = Math.abs(localpg[i][j+1] -
               2.0 * localpg[i][j] + localpg[i][j-1]) / localpg[i][j];
        }

      double send_sxi[] = new double[jmax+1];
      double recv_sxi[] = new double[jmax+1];

      for(j=1;j<jmax;j++)
        send_sxi[j] = local_sxi[n_bar][j];

      MPI.COMM_WORLD.Sendrecv(send_sxi,0,jmax+1,MPI.DOUBLE,right,1,
                              recv_sxi,0,jmax+1,MPI.DOUBLE,left ,1);
      if(my_rank != 0)
        for(j=1;j<jmax;j++)
          local_sxi[0][j] = recv_sxi[j];

      for(j=1;j<jmax;j++)
        send_sxi[j] = local_sxi[1][j];

      MPI.COMM_WORLD.Sendrecv(send_sxi,0,jmax+1,MPI.DOUBLE,left ,1,
                              recv_sxi,0,jmax+1,MPI.DOUBLE,right,1);
      if (my_rank != (np-1))
        for(j=1;j<jmax;j++)
          local_sxi[n_bar+1][j] = recv_sxi[j];

      double send_ug[] = new double[jmax*8];
      double recv_ug[] = new double[jmax*8];      
      int initUG = 0;
      Statevector hold_n2[][] = new Statevector[2][jmax+2];

      initUG = 0;
      for(j=1;j<jmax;j++){
        send_ug[initUG++] = localug[n_bar][j].a;
        send_ug[initUG++] = localug[n_bar][j].b;
        send_ug[initUG++] = localug[n_bar][j].c;
        send_ug[initUG++] = localug[n_bar][j].d;
        send_ug[initUG++] = localug[n_bar-1][j].a;
        send_ug[initUG++] = localug[n_bar-1][j].b;
        send_ug[initUG++] = localug[n_bar-1][j].c;
        send_ug[initUG++] = localug[n_bar-1][j].d;
      }
      MPI.COMM_WORLD.Sendrecv(send_ug,0,jmax*8,MPI.DOUBLE,right,1,
                              recv_ug,0,jmax*8,MPI.DOUBLE,left ,1);
      initUG = 0;
      if(my_rank != 0)
        for(j=1;j<jmax;j++){
          localug[0][j].a = recv_ug[initUG++];
          localug[0][j].b = recv_ug[initUG++];
          localug[0][j].c = recv_ug[initUG++];
          localug[0][j].d = recv_ug[initUG++];

          hold_n2[0][j] = new Statevector();
          hold_n2[0][j].a = recv_ug[initUG++];
          hold_n2[0][j].b = recv_ug[initUG++];
          hold_n2[0][j].c = recv_ug[initUG++];
          hold_n2[0][j].d = recv_ug[initUG++];
        }

      initUG = 0;
      for(j=1;j<jmax;j++){
        send_ug[initUG++] = localug[1][j].a;
        send_ug[initUG++] = localug[1][j].b;
        send_ug[initUG++] = localug[1][j].c;
        send_ug[initUG++] = localug[1][j].d;
 
        send_ug[initUG++] = localug[2][j].a;
        send_ug[initUG++] = localug[2][j].b;
        send_ug[initUG++] = localug[2][j].c;
        send_ug[initUG++] = localug[2][j].d;
      }
      MPI.COMM_WORLD.Sendrecv(send_ug,0,jmax*8,MPI.DOUBLE,left,1,
                              recv_ug,0,jmax*8,MPI.DOUBLE,right,1);
      initUG = 0;
      if (my_rank != (np-1))
        for(j=1;j<jmax;j++){
          localug[n_bar+1][j].a = recv_ug[initUG++];
          localug[n_bar+1][j].b = recv_ug[initUG++];
          localug[n_bar+1][j].c = recv_ug[initUG++];
          localug[n_bar+1][j].d = recv_ug[initUG++];

          hold_n2[1][j] = new Statevector();
          hold_n2[1][j].a = recv_ug[initUG++];
          hold_n2[1][j].b = recv_ug[initUG++];
          hold_n2[1][j].c = recv_ug[initUG++];
          hold_n2[1][j].d = recv_ug[initUG++];
        }

      /* Then calculate the fluxes ..... */
      for (i = 1; i <= n_bar; ++i) {
        for (j = 1; j < jmax; ++j) {
          
          /* Clear values */
          /* East Face */
          if (i > 1 && i < n_bar) {
            adt = (local_a[i][j] + local_a[i+1][j]) / 
                  (local_deltat[i][j] + local_deltat[i+1][j]);
            sbar = (local_sxi[i+1][j] + local_sxi[i][j]) * 0.5;
          }
          else {
            adt = local_a[i][j]/local_deltat[i][j];
            sbar = local_sxi[i][j];
          }

          if((my_rank < np-1 && i == n_bar)||(my_rank > 0 && i == 1)){
             adt = (local_a[i][j] + local_a[i+1][j]) / 
                   (local_deltat[i][j] + local_deltat[i+1][j]);
             sbar = (local_sxi[i+1][j] + local_sxi[i][j]) * 0.5;
          }
          tempdouble = nu2*sbar*adt;
          scrap2.a = tempdouble * (localug[i+1][j].a-localug[i][j].a);
          scrap2.b = tempdouble * (localug[i+1][j].b-localug[i][j].b);
          scrap2.c = tempdouble * (localug[i+1][j].c-localug[i][j].c);
          scrap2.d = tempdouble * (localug[i+1][j].d-localug[i][j].d);
          if (i > 1 && i < n_bar) {
            temp = localug[i+2][j].svect(localug[i-1][j]);
            temp2.a = 3.0*(localug[i][j].a-localug[i+1][j].a);
            temp2.b = 3.0*(localug[i][j].b-localug[i+1][j].b);
            temp2.c = 3.0*(localug[i][j].c-localug[i+1][j].c);
            temp2.d = 3.0*(localug[i][j].d-localug[i+1][j].d);

            tempdouble = -nu4*adt;
            scrap4.a = tempdouble*(temp.a+temp2.a);
            scrap4.b = tempdouble*(temp.a+temp2.b);
            scrap4.c = tempdouble*(temp.a+temp2.c);
            scrap4.d = tempdouble*(temp.a+temp2.d);
          }
          else {
            scrap4.a = 0.0;
            scrap4.b = 0.0;
            scrap4.c = 0.0;
            scrap4.d = 0.0;
          }

          if((my_rank < np-1 && i == n_bar)||(my_rank > 0 && i == 1)){
            if (my_rank < np-1 && i==n_bar)
               temp = hold_n2[1][j].svect(localug[i-1][j]);
            else 
               temp = localug[i+2][j].svect(localug[i-1][j]);
            temp2.a = 3.0*(localug[i][j].a-localug[i+1][j].a);
            temp2.b = 3.0*(localug[i][j].b-localug[i+1][j].b);
            temp2.c = 3.0*(localug[i][j].c-localug[i+1][j].c);
            temp2.d = 3.0*(localug[i][j].d-localug[i+1][j].d);

            tempdouble = -nu4*adt;
            scrap4.a = tempdouble*(temp.a+temp2.a);
            scrap4.b = tempdouble*(temp.a+temp2.b);
            scrap4.c = tempdouble*(temp.a+temp2.c);
            scrap4.d = tempdouble*(temp.a+temp2.d);
          }
          temp.a = scrap2.a + scrap4.a;
          temp.b = scrap2.b + scrap4.b;
          temp.c = scrap2.c + scrap4.c;
          temp.d = scrap2.d + scrap4.d;
          local_d[i][j] = temp;

          /* West Face */ 
          if(i > 1 && i < n_bar) {
            adt = (local_a[i][j] + local_a[i-1][j]) / 
                  (local_deltat[i][j] + local_deltat[i-1][j]);
            sbar = (local_sxi[i][j] + local_sxi[i-1][j]) *0.5;
          }
          else {
            adt = local_a[i][j]/local_deltat[i][j];
            sbar = local_sxi[i][j];
          }
          if((my_rank < np-1 && i == n_bar)||(my_rank > 0 && i == 1)){
            adt = (local_a[i][j] + local_a[i-1][j]) / 
                  (local_deltat[i][j] + local_deltat[i-1][j]);
            sbar = (local_sxi[i][j] + local_sxi[i-1][j]) *0.5;
          }
          tempdouble = -nu2*sbar*adt;
          scrap2.a = tempdouble * (localug[i][j].a-localug[i-1][j].a);
          scrap2.b = tempdouble * (localug[i][j].b-localug[i-1][j].b);
          scrap2.c = tempdouble * (localug[i][j].c-localug[i-1][j].c);
          scrap2.d = tempdouble * (localug[i][j].d-localug[i-1][j].d);
          if (i > 1 && i < n_bar) {
            temp = localug[i+1][j].svect(localug[i-2][j]);
            temp2.a = 3.0*(localug[i-1][j].a-localug[i][j].a);
            temp2.b = 3.0*(localug[i-1][j].b-localug[i][j].b);
            temp2.c = 3.0*(localug[i-1][j].c-localug[i][j].c);
            temp2.d = 3.0*(localug[i-1][j].d-localug[i][j].d);

            tempdouble = nu4*adt;
            scrap4.a = tempdouble*(temp.a+temp2.a);
            scrap4.b = tempdouble*(temp.a+temp2.b);
            scrap4.c = tempdouble*(temp.a+temp2.c);
            scrap4.d = tempdouble*(temp.a+temp2.d);
          }
          else {
            scrap4.a = 0.0;
            scrap4.b = 0.0;
            scrap4.c = 0.0;
            scrap4.d = 0.0;
          }
          if((my_rank < np-1 && i == n_bar)||(my_rank > 0 && i == 1)){
            if(my_rank > 0 && i==1)
               temp = localug[i+1][j].svect(hold_n2[0][j]);
            else
               temp = localug[i+1][j].svect(localug[i-2][j]);
            temp2.a = 3.0*(localug[i-1][j].a-localug[i][j].a);
            temp2.b = 3.0*(localug[i-1][j].b-localug[i][j].b);
            temp2.c = 3.0*(localug[i-1][j].c-localug[i][j].c);
            temp2.d = 3.0*(localug[i-1][j].d-localug[i][j].d);

            tempdouble = nu4*adt;
            scrap4.a = tempdouble*(temp.a+temp2.a);
            scrap4.b = tempdouble*(temp.a+temp2.b);
            scrap4.c = tempdouble*(temp.a+temp2.c);
            scrap4.d = tempdouble*(temp.a+temp2.d);
          }

          local_d[i][j].a += scrap2.a + scrap4.a;
          local_d[i][j].b += scrap2.b + scrap4.b;
          local_d[i][j].c += scrap2.c + scrap4.c;
          local_d[i][j].d += scrap2.d + scrap4.d;

          /* North Face */
          if (j > 1 && j < jmax-1) {
            adt = (local_a[i][j] + local_a[i][j+1]) / 
                  (local_deltat[i][j] + local_deltat[i][j+1]);
            sbar = (local_seta[i][j] + local_seta[i][j+1]) * 0.5;
          }
          else {
            adt = local_a[i][j]/local_deltat[i][j];
            sbar = local_seta[i][j];
          }
          tempdouble = nu2*sbar*adt;
          scrap2.a = tempdouble * (localug[i][j+1].a-localug[i][j].a);
          scrap2.b = tempdouble * (localug[i][j+1].b-localug[i][j].b);
          scrap2.c = tempdouble * (localug[i][j+1].c-localug[i][j].c);
          scrap2.d = tempdouble * (localug[i][j+1].d-localug[i][j].d);

          if (j > 1 && j < jmax-1) {
            temp = localug[i][j+2].svect(localug[i][j-1]);
            temp2.a = 3.0*(localug[i][j].a-localug[i][j+1].a);
            temp2.b = 3.0*(localug[i][j].b-localug[i][j+1].b);
            temp2.c = 3.0*(localug[i][j].c-localug[i][j+1].c);
            temp2.d = 3.0*(localug[i][j].d-localug[i][j+1].d);

            tempdouble = -nu4*adt;
            scrap4.a = tempdouble*(temp.a+temp2.a);
            scrap4.b = tempdouble*(temp.a+temp2.b);
            scrap4.c = tempdouble*(temp.a+temp2.c);
            scrap4.d = tempdouble*(temp.a+temp2.d);
          }
          else {
            scrap4.a = 0.0;
            scrap4.b = 0.0;
            scrap4.c = 0.0;
            scrap4.d = 0.0;
          }
          local_d[i][j].a += scrap2.a + scrap4.a;
          local_d[i][j].b += scrap2.b + scrap4.b;
          local_d[i][j].c += scrap2.c + scrap4.c;
          local_d[i][j].d += scrap2.d + scrap4.d;
          
          /* South Face */
          if (j > 1 && j < jmax-1) {
            adt = (local_a[i][j] + local_a[i][j-1]) / 
                  (local_deltat[i][j] + local_deltat[i][j-1]);
            sbar = (local_seta[i][j] + local_seta[i][j-1]) * 0.5;
          }
          else {
            adt = local_a[i][j]/local_deltat[i][j];
            sbar = local_seta[i][j];
          }
          tempdouble = -nu2*sbar*adt;
          scrap2.a = tempdouble * (localug[i][j].a-localug[i][j-1].a);
          scrap2.b = tempdouble * (localug[i][j].b-localug[i][j-1].b);
          scrap2.c = tempdouble * (localug[i][j].c-localug[i][j-1].c);
          scrap2.d = tempdouble * (localug[i][j].d-localug[i][j-1].d);

          if (j > 1 && j < jmax-1) {
            temp = localug[i][j+1].svect(localug[i][j-2]);
            temp2.a = 3.0*(localug[i][j-1].a-localug[i][j].a);
            temp2.b = 3.0*(localug[i][j-1].b-localug[i][j].b);
            temp2.c = 3.0*(localug[i][j-1].c-localug[i][j].c);
            temp2.d = 3.0*(localug[i][j-1].d-localug[i][j].d);

            tempdouble = nu4*adt;
            scrap4.a = tempdouble*(temp.a+temp2.a);
            scrap4.b = tempdouble*(temp.a+temp2.b);
            scrap4.c = tempdouble*(temp.a+temp2.c);
            scrap4.d = tempdouble*(temp.a+temp2.d);
          }
          else {
            scrap4.a = 0.0;
            scrap4.b = 0.0;
            scrap4.c = 0.0;
            scrap4.d = 0.0;
          }
          local_d[i][j].a += scrap2.a + scrap4.a;
          local_d[i][j].b += scrap2.b + scrap4.b;
          local_d[i][j].c += scrap2.c + scrap4.c;
          local_d[i][j].d += scrap2.d + scrap4.d;
        }
      }
  }
  
  private void calculateDeltaT() {
    double xeta, yeta, xxi, yxi;              /* Local change in x and y */
    int i, j;
    double mint;
    double c, q, r;
    double safety_factor = 0.3;
    
    for (i = 1; i <= n_bar; ++i)
      for (j = 1; j < jmax; ++j) {
        xxi = (local_xnode[i][j] - local_xnode[i-1][j] 
               + local_xnode[i][j-1] - local_xnode[i-1][j-1]) * 0.5;
        yxi = (local_ynode[i][j] - local_ynode[i-1][j] 
               + local_ynode[i][j-1] - local_ynode[i-1][j-1]) * 0.5;
        xeta = (local_xnode[i][j] - local_xnode[i][j-1] 
               + local_xnode[i-1][j] - local_xnode[i-1][j-1]) * 0.5;
        yeta = (local_ynode[i][j] - local_ynode[i][j-1] 
               + local_ynode[i-1][j] - local_ynode[i-1][j-1]) * 0.5;
        
        q = (yeta * local_ug[i][j].b - xeta * local_ug[i][j].c);
        r = (-yxi * local_ug[i][j].b + xxi * local_ug[i][j].c);
        c = Math.sqrt (gamma * rgas * local_tg[i][j]);

        local_deltat[i][j] = safety_factor * 2.8284 * local_a[i][j] /

          ( (Math.abs(q) + Math.abs(r))/local_ug[i][j].a + c * 
           Math.sqrt(xxi*xxi + yxi*yxi + xeta*xeta + yeta*yeta +
                                  2.0 * Math.abs(xeta*xxi + yeta*yxi)));
      }

    /* If that's the user's choice, make it time accurate */
    if (ntime == 1) {
      mint = 100000.0;
      for (i = 1; i <= n_bar; ++i)
        for (j = 1; j < jmax; ++j)
          if (local_deltat[i][j] < mint)
            mint = local_deltat[i][j];
      
      for (i = 1; i <= n_bar; ++i)
        for (j = 1; j < jmax; ++j)
          local_deltat[i][j] = mint;
    }

    double send_deltat[] = new double[jmax];
    double recv_deltat[] = new double[jmax];

    int right = my_rank + 1;
    int left  = my_rank - 1;
    if (my_rank == 0)
      left = np -1;
    else if(my_rank == (np -1))
      right = 0;

    for(j=1;j<jmax;j++)
      send_deltat[j] = local_deltat[n_bar][j];

    MPI.COMM_WORLD.Sendrecv(send_deltat,0,jmax,MPI.DOUBLE,right,1,
                            recv_deltat,0,jmax,MPI.DOUBLE,left ,1);
    if(my_rank != 0)
      for(j=1;j<jmax;j++)
        local_deltat[0][j] = recv_deltat[j];

    for(j=1;j<jmax;j++)
      send_deltat[j] = local_deltat[1][j];

    MPI.COMM_WORLD.Sendrecv(send_deltat,0,jmax,MPI.DOUBLE,left ,1,
                            recv_deltat,0,jmax,MPI.DOUBLE,right,1);
    if (my_rank != (np-1))
      for(j=1;j<jmax;j++)
        local_deltat[n_bar+1][j] = recv_deltat[j];
  }

  private void calculateDummyCells(double localpg[][],
                                   double localtg[][], Statevector localug[][]) {
    double c;
    double jminus;
    double jplus;
    double s;
    double rho, temp, u, v;
    double scrap, scrap2;
    double theta;
    double uprime;
    int i, j;
    Vector2 norm = new Vector2();
    Vector2 tan = new Vector2();
    Vector2 u1 = new Vector2();

    uff = machff;
    jplusff = uff + 2.0 / (gamma - 1.0) * cff;
    jminusff = uff - 2.0 / (gamma - 1.0) * cff;

    for (i = 1; i <= n_bar; ++i) {
      /* Bottom wall boundary cells */
      /* Routine checked by brute force for initial conditions. */
      /* Routine checked by brute force for random conditions. */
      /* Construct tangent vectors */
      tan.ihat = local_xnode[i][0] - local_xnode[i-1][0];
      tan.jhat = local_ynode[i][0] - local_ynode[i-1][0];
      norm.ihat = - (local_ynode[i][0] - local_ynode[i-1][0]);
      norm.jhat = local_xnode[i][0] - local_xnode[i-1][0];
      
      scrap = tan.magnitude();
      tan.ihat = tan.ihat / scrap;
      tan.jhat = tan.jhat / scrap;
      scrap = norm.magnitude();
      norm.ihat = norm.ihat / scrap;
      norm.jhat = norm.jhat / scrap;
      
      /* Now set some state variables */
      rho = localug[i][1].a;
      localtg[i][0] = localtg[i][1];
      u1.ihat = localug[i][1].b / rho;
      u1.jhat = localug[i][1].c / rho;
      
      u = u1.dot(tan) + u1.dot(norm) * tan.jhat /norm.jhat;
      u = u / (tan.ihat - (norm.ihat * tan.jhat / norm.jhat));
      
      v = - (u1.dot(norm) + u * norm.ihat) / norm.jhat;
      
      /* And construct the new state vector */
      localug[i][0].a = localug[i][1].a;
      localug[i][0].b = rho * u;
      localug[i][0].c = rho * v;
      localug[i][0].d = rho * (Cv * localtg[i][0] + 0.5 * (u*u + v*v));
      localpg[i][0] = localpg[i][1];
      
      /* Top Wall Boundary Cells */
      /* Checked numerically for default conditions. */
      /* Construct normal and tangent vectors */
      /* This part checked and works; it produces the correct vectors */
      tan.ihat = local_xnode[i][jmax-1] - local_xnode[i-1][jmax-1];
      tan.jhat = local_ynode[i][jmax-1] - local_ynode[i-1][jmax-1];
      norm.ihat = local_ynode[i][jmax-1] - local_ynode[i-1][jmax-1];
      norm.jhat = -(local_xnode[i][jmax-1] - local_xnode[i-1][jmax-1]);
      
      scrap = tan.magnitude();
      tan.ihat = tan.ihat / scrap;
      tan.jhat = tan.jhat / scrap;
      scrap = norm.magnitude();
      norm.ihat = norm.ihat / scrap;
      norm.jhat = norm.jhat / scrap;
     
      /* Now set some state variables */
      rho = localug[i][jmax-1].a;
      temp = localtg[i][jmax-1];
      u1.ihat = localug[i][jmax-1].b / rho;
      u1.jhat = localug[i][jmax-1].c / rho;
     
      u = u1.dot(tan) + u1.dot(norm) * tan.jhat /norm.jhat;
      u = u / (tan.ihat - (norm.ihat * tan.jhat / norm.jhat));
     
      v = - (u1.dot(norm) + u * norm.ihat) / norm.jhat;
      
      /* And construct the new state vector */
      localug[i][jmax].a = localug[i][jmax-1].a;
      localug[i][jmax].b = rho * u;
      localug[i][jmax].c = rho * v;
      localug[i][jmax].d = rho * (Cv * temp + 0.5 * (u*u + v*v));
      localtg[i][jmax] = temp;
      localpg[i][jmax] = localpg[i][jmax-1];
    }

   for (j = 1; j < jmax; ++j) {
     /* Inlet Boundary Cells: unchecked */
     /* Construct the normal vector; */
     norm.ihat = local_ynode[0][j-1] - local_ynode[0][j];
     norm.jhat = local_xnode[0][j] - local_xnode[0][j-1];
     scrap = norm.magnitude();
     norm.ihat = norm.ihat / scrap;
     norm.jhat = norm.jhat / scrap;
     theta = Math.acos((local_ynode[0][j-1] - local_ynode[0][j]) / 
      Math.sqrt((local_xnode[0][j] - local_xnode[0][j-1])
                *(local_xnode[0][j] - local_xnode[0][j-1]) 
               + (local_ynode[0][j-1] - local_ynode[0][j]) 
               * (local_ynode[0][j-1] - local_ynode[0][j])));
     
     u1.ihat = localug[1][j].b / localug[1][j].a;
     u1.jhat = localug[1][j].c / localug[1][j].a;
     uprime = u1.ihat * Math.cos(theta);
     c = Math.sqrt(gamma * rgas * localtg[1][j]);
     /* Supersonic inflow; works on the initial conditions */
     if (uprime < -c) {
       /* Use far field conditions */
       localug[0][j].a = rhoff;
       localug[0][j].b = rhoff * uff;
       localug[0][j].c = rhoff * vff;
       localug[0][j].d = rhoff * (Cv * tff + 0.5 * (uff*uff + vff*vff));
       localtg[0][j] = tff;
       localpg[0][j] = pff;
     }
     /* Subsonic inflow */
     /* This works on the initial conditions */
     else if(uprime < 0.0) {
       /* Calculate Riemann invarients here */
       jminus = u1.ihat - 2.0/(gamma-1.0) * c;
       s = Math.log(pff) - gamma * Math.log(rhoff);
       v = vff;
       
       u = (jplusff + jminus) / 2.0;
       scrap = (jplusff - u) * (gamma-1.0) * 0.5;
       localtg[0][j] = (1.0 / (gamma * rgas)) * scrap * scrap;
       localpg[0][j] = Math.exp(s) / Math.pow((rgas * localtg[0][j]), gamma);
       localpg[0][j] = Math.pow(localpg[0][j], 1.0 / (1.0 - gamma));
       
       /* And now: construct the new state vector */

       localug[0][j].a = localpg[0][j] / (rgas * localtg[0][j]);
       localug[0][j].b = localug[0][j].a * u;
       localug[0][j].c = localug[0][j].a * v;
       localug[0][j].d = localug[0][j].a * (Cv * tff + 0.5 * (u*u + v*v));
     }
     /* Other options: */
     /* We should throw an exception here. */
     else {
       System.err.println("You have outflow at the inlet, which is not allowed.");
     }
     /* Outlet Boundary Cells */
     /* Construct the normal vector; works. */
     norm.ihat = local_ynode[0][j] - local_ynode[0][j-1];
     norm.jhat = local_xnode[0][j-1] - local_xnode[0][j];
     scrap = norm.magnitude();
     norm.ihat = norm.ihat / scrap;
     norm.jhat = norm.jhat / scrap;
     scrap = local_xnode[0][j-1] - local_xnode[0][j];
     scrap2 = local_ynode[0][j] - local_ynode[0][j-1];
     theta = Math.acos((local_ynode[0][j] - local_ynode[0][j-1]) / 
                       Math.sqrt(scrap*scrap + scrap2*scrap2));
     u1.ihat = localug[n_bar][j].b / localug[n_bar][j].a;
     u1.jhat = localug[n_bar][j].c / localug[n_bar][j].a;
     uprime = u1.ihat * Math.cos(theta);
     c = Math.sqrt(gamma * rgas * localtg[n_bar][j]);
     /* Supersonic outflow; works for defaults conditions. */
     if (uprime > c){
     /* Use a backward difference 2nd order derivative approximation */
     /* To set values at exit */
       localug[n_bar+1][j].a = 
              2.0 * localug[n_bar][j].a - localug[n_bar-1][j].a;
       localug[n_bar+1][j].b = 
              2.0 * localug[n_bar][j].b - localug[n_bar-1][j].b;
       localug[n_bar+1][j].c = 
              2.0 * localug[n_bar][j].c - localug[n_bar-1][j].c;
       localug[n_bar+1][j].d = 
              2.0 * localug[n_bar][j].d - localug[n_bar-1][j].d;
        
       localpg[n_bar+1][j] = 2.0 * localpg[n_bar][j] - localpg[n_bar-1][j];
       localtg[n_bar+1][j] = 2.0 * localtg[n_bar][j] - localtg[n_bar-1][j];
     }
     /* Subsonic Outflow; works for defaults conditions. */
     else if (uprime < c && uprime > 0) {
       jplus = u1.ihat + 2.0/(gamma - 1) * c;
       v = localug[n_bar][j].c / localug[n_bar][j].a;
       s = Math.log(localpg[n_bar][j]) -
                    gamma * Math.log(localug[n_bar][j].a);
       
       u = (jplus + jminusff) / 2.0;
       scrap =(jplus - u)* (gamma-1.0) * 0.5;
       localtg[n_bar+1][j] = (1.0 / (gamma * rgas)) * scrap * scrap;
       localpg[n_bar+1][j] = Math.exp(s) / 
                       Math.pow((rgas * localtg[n_bar+1][j]), gamma);
       localpg[n_bar+1][j] = Math.pow(localpg[n_bar+1][j], 1.0 / 
                             (1.0-gamma));
       rho = localpg[n_bar+1][j]/ (rgas * localtg[n_bar+1][j]);
       
       /* And now, construct the new state vector. */
       localug[n_bar+1][j].a = rho;
       localug[n_bar+1][j].b = rho * u;
       localug[n_bar+1][j].c = rho * v;
       localug[n_bar+1][j].d = rho * (Cv * localtg[n_bar+1][j] 
                   + 0.5 * (u*u + v*v));
       
     }

     /* Other cases that shouldn't have to be used. */
     if (uprime < -c) {
       /* Supersonic inflow.  Use far field conditions */
       localug[0][j].a = rhoff;
       localug[0][j].b = rhoff * uff;
       localug[0][j].c = rhoff * vff;
       localug[0][j].d = rhoff * (Cv * tff + 0.5 * (uff*uff + vff*vff));
       localtg[0][j] = tff;
       localpg[0][j] = pff;
     }
     /* Subsonic inflow.  This works on the initial conditions. */
     else if(uprime < 0.0) {

       /* Debug: throw exception here? */
       /* Calculate Riemann invarients here. */

       jminus = u1.ihat - 2.0/(gamma-1.0) * c;
       s = Math.log(pff) - gamma * Math.log(rhoff);
       v = vff;
       
       u = (jplusff + jminus) / 2.0;
       scrap = (jplusff - u)* (gamma-1.0) * 0.5;
       localtg[0][j] = (1.0 / (gamma * rgas)) * scrap * scrap;
       localpg[0][j] = Math.exp(s) / Math.pow((rgas * localtg[0][j]), gamma);
       localpg[0][j] = Math.pow(localpg[0][j], 1.0 / (1.0 - gamma));
       
       /* And now: construct the new state vector. */

       localug[0][j].a = localpg[0][j] / (rgas * localtg[0][j]);
       localug[0][j].b = localug[0][j].a * u;
       localug[0][j].c = localug[0][j].a * v;
       localug[0][j].d = localug[0][j].a * (Cv * tff + 0.5 * (u*u + v*v));
     }
     else if ((uprime < c && uprime > 0) || (uprime > c)){}
     /* Other Options */
     /* Debug: throw exception here? */
     else {
       System.err.println("You have inflow at the outlet, which is not allowed.");
     }
   } 
   /* Do something with corners to avoid division by zero errors */
   /* What you do shouldn't matter */
   localug[0][0] = localug[1][0];
   localug[n_bar+1][0] = localug[n_bar+1][1];
   localug[0][jmax] = localug[1][jmax];
   localug[n_bar+1][jmax] = localug[n_bar+1][jmax-1];
  }

  void drawEdge(Graphics g, int width, int height) {
    /* Draws the outside border of the grid (as opposed to the frame) */
    int gridheight,  /* height of drawing region for grid in pixels */
        gridwidth;   /* Width of same */
    int i, j;
    int x1=0, x2, y1=0, y2;
    int xoffset, yoffset;
    double scale;

    gridheight = height;  /* Leave room at bottom for labels */
    gridwidth = width;
    scale = Math.min(gridwidth/(xmax-xmin), gridheight/(ymax-ymin));
    xoffset = (int) (xmin * scale);
    yoffset = (int) (ymin * scale);

    g.setColor(Color.white);
    for (i = 0; i < imax-1; ++i) {
      x1 = (int) (scale * xnode[i][0]) - xoffset;
      y1 = gridheight - (int) (scale * ynode[i][0]) - yoffset;
      x2 = (int) (scale * xnode[i+1][0]) - xoffset;
      y2 = gridheight - (int) (scale * ynode[i+1][0]) - yoffset;
      g.drawLine(x1, y1, x2, y2);

      x1 = (int) (scale * xnode[i][jmax-1]) - xoffset;
      y1 = gridheight - (int) (scale * ynode[i][jmax-1]) - yoffset;
      x2 = (int) (scale * xnode[i+1][jmax-1]) - xoffset;
      y2 = gridheight - (int) (scale * ynode[i+1][jmax-1]) - yoffset;
      g.drawLine(x1, y1, x2, y2);
    }
    for (j =0; j < jmax-1; ++j) {
      x1 = (int) (scale * xnode[imax-1][j]) - xoffset;
      y1 = gridheight - (int) (scale * ynode[imax-1][j]) - yoffset;
      x2 = (int) (scale * xnode[imax-1][j+1]) - xoffset;
      y2 = gridheight - (int) (scale * ynode[imax-1][j+1]) - yoffset;
      g.drawLine(x1, y1, x2, y2);

      x1 = (int) (scale * xnode[0][j]) - xoffset;
      y1 = gridheight - (int) (scale * ynode[0][j]) - yoffset;
      x2 = (int) (scale * xnode[0][j+1]) - xoffset;
      y2 = gridheight - (int) (scale * ynode[0][j+1]) - yoffset;
      g.drawLine(x1, y1, x2, y2);
    }
    int k = 1;
    String s;
    int tempX1=0;

    for (i = 0; i < imax-1; ++i)
      for (j = 0; j < jmax-1; ++j) {
        x1 = (int) (scale * xnode[i][j]) - xoffset;
        y1 = gridheight - (int) (scale * ynode[i][j]) - yoffset;
        x2 = (int) (scale * xnode[i+1][j]) - xoffset;
        y2 = gridheight - (int) (scale * ynode[i+1][j]) - yoffset;
        if (i == n_bar*k){
          g.setColor(Color.lightGray);
          g.fillRect(x2-3,0,6,height);
//        g.drawLine(x2, 0, x2, height);
          g.setColor(Color.yellow);
          s = "Proc "+k;
          tempX1 = (int)(tempX1+x1)/2;
          g.drawString(s,tempX1,(int)height/2);
          tempX1 = x1;
          k++;
        }
    }
    s = "Proc "+k;
    tempX1 = (int)(tempX1+x1)/2;
    g.drawString(s,tempX1,(int)height/2);
  }

  void drawTufts(Graphics g, int width, int height) {
    Color tuftscolor;
    double scale, theta;
    int gridwidth, gridheight;
    int i, j;
    int x1, x2, y1, y2;
    int xoffset, yoffset;

    tuftscolor = new Color(150,150,150);

    gridheight = height;  /* Leave room for labels if necessary */
    gridwidth = width;

    g.setColor(tuftscolor);
    scale = Math.min(gridwidth/(xmax-xmin), gridheight/(ymax-ymin));
    xoffset = (int) (xmin * scale);
    yoffset = (int) (ymin * scale);
    for(i = 0; i < imax-1; ++i)
      for (j = 0; j < jmax-1; ++j) {
        x1 = (int) (scale * xnode[i][j]) - xoffset;
        y1 = gridheight - (int) (scale * ynode[i][j]) - yoffset;
        if (Math.abs(ug[i][j].b) > 0.00001) {
            theta = -Math.atan(ug[i][j].c/ug[i][j].b);
          }
        else {
          if (ug[i][j].c > 0.0)
            theta = -Math.PI * 0.5;
          else
            theta = -Math.PI * 1.5;
        }
        x2 = x1 + (int) (6.0*Math.cos(theta));
        y2 = y1 + (int) (6.0*Math.sin(theta));
        g.drawLine(x1, y1, x2, y2);
      }
  }

  void drawGrid(Graphics g, int width, int height) {
    /* Draws the grid; does not draw outside edges */
    Color gridcolor;
    int i, j;               /* Dummy Counters */
    int gridheight,         /* height of drawing region for grid in pixels */
        gridwidth;          /* Width of same */
    int x1, y1, x2, y2;
    int xoffset, yoffset;
    double scale;
    java.awt.Dimension x;

    gridheight = height;    /* Leave room for labels if necessary */
    gridwidth = width;

    gridcolor = new Color(100,100,100);
    g.setColor(gridcolor);
    scale = Math.min(gridwidth/(xmax-xmin), gridheight/(ymax-ymin));
    xoffset = (int) (xmin * scale);
    yoffset = (int) (ymin * scale);
    for (i = 0; i < imax-1; ++i)
      for (j = 0; j < jmax-1; ++j) {
        x1 = (int) (scale * xnode[i][j]) - xoffset;
        y1 = gridheight - (int) (scale * ynode[i][j]) - yoffset;
        x2 = (int) (scale * xnode[i+1][j]) - xoffset;
        y2 = gridheight - (int) (scale * ynode[i+1][j]) - yoffset;
        
        g.drawLine(x1, y1, x2, y2);
        x2 = (int) (scale * xnode[i][j+1]) - xoffset;
        y2 = gridheight - (int) (scale * ynode[i][j+1]) - yoffset;
        g.drawLine(x1, y1, x2, y2);
      }
    for (i = 0; i < imax-1; ++i) {
      x1 = (int) (scale * xnode[i][jmax-1]) - xoffset;
      y1 = gridheight - (int) (scale * ynode[i][jmax-1]) - yoffset;
      x2 = (int) (scale * xnode[i+1][jmax-1]) - xoffset;
      y2 = gridheight - (int) (scale * ynode[i+1][jmax-1]) - yoffset;
      g.drawLine(x1, y1, x2, y2);
    }
    for (j =0; j < jmax-1; ++j) {
      x1 = (int) (scale * xnode[imax-1][j]) - xoffset;
      y1 = gridheight - (int) (scale * ynode[imax-1][j]) - yoffset;
      x2 = (int) (scale * xnode[imax-1][j+1]) - xoffset;
      y2 = gridheight - (int) (scale * ynode[imax-1][j+1]) - yoffset;
      g.drawLine(x1, y1, x2, y2);
    }
  }

  void drawLinePlot(Graphics g, int ndata, int width, int height) {
    int gridheight, gridwidth;
    int jmid;
    int i, j;
    int plotheight, plotwidth;
    double plotdata[][];
    double scalex, scaley;
    double mach, pres;
    double plotmin, plotmax;
    double temp, vel, c;
    double u, v;
    int xoffset, yoffset;
    int n, x, y;
    int x1, y1, x2, y2;
    int xmargin = 55;
    int ymargin = 30;

    gridheight = height;
    gridwidth = width;
    plotheight = height-ymargin-20;  /* Leave room at bottom for labels */
    plotwidth = width-xmargin-20;

    plotmin=Double.MAX_VALUE;
    plotmax = Double.MIN_VALUE;
    plotdata = new double[3][imax];
    jmid = (jmax + 1) / 2;
    for (n = 0; n < 3; ++n) {
      if (n == 0)
        j = 1;
      else if (n==1)
        j = jmid;
      else
        j = jmax;
      
      for (i=1; i < imax+1; ++i) {
        switch(ndata) {
        case 0:                /* Density */
          plotdata[n][i-1] = ug[i][j].a;
          plotmin = Math.min(plotdata[n][i-1], plotmin);
          plotmax = Math.max(plotdata[n][i-1], plotmax);
          break;
        case 1:                /* Pressure */
          u = ug[i][j].b/ug[i][j].a;
          v = ug[i][j].c/ug[i][j].a;
          plotdata[n][i-1] = ug[i][j].d/ug[i][j].a - (0.5 * (u*u + v*v));
          plotdata[n][i-1] = plotdata[n][i-1]/Cv;
          plotdata[n][i-1] = plotdata[n][i-1] * rgas * ug[i][j].a;
          plotmax = Math.max(plotmax, plotdata[n][i-1]);
          plotmin = Math.min(plotmin, plotdata[n][i-1]);
          break;
        case 2:                /* Temperature */
          u = ug[i][j].b/ug[i][j].a;
          v = ug[i][j].c/ug[i][j].a;
          plotdata[n][i-1] = ug[i][j].d/ug[i][j].a - (0.5 * (u*u + v*v));
          plotdata[n][i-1] = plotdata[n][i-1]/Cv;
          plotmax = Math.max(plotmax, plotdata[n][i-1]);
          plotmin = Math.min(plotmin, plotdata[n][i-1]);
          break;
        case 3:                /* Mach Number */
          u = ug[i][j].b/ ug[i][j].a;
          v = ug[i][j].c / ug[i][j].a;
          vel = Math.sqrt(u*u + v*v);
          temp = ug[i][j].d/ug[i][j].a - (0.5 * (u*u + v*v));
          temp = temp/Cv;
          c = Math.sqrt(gamma * rgas * temp);
          plotdata[n][i-1] = vel / c;
          plotmax = Math.max(plotmax, plotdata[n][i-1]);
          plotmin = Math.min(plotmin, plotdata[n][i-1]);
        break;
        case 4:                /* Axial Velocity */
          plotdata[n][i-1] = ug[i][j].b / ug[i][j].a;
          plotmax = Math.max(plotmax, plotdata[n][i-1]);
          plotmin = Math.min(plotmin, plotdata[n][i-1]);
          break;
        case 5:                /* Normal Velocity */
          plotdata[n][i-1] = ug[i][j].c / ug[i][j].a;
          plotmax = Math.max(plotmax, plotdata[n][i-1]);
          plotmin = Math.min(plotmin, plotdata[n][i-1]);
          break;
        case 6:                /* Stagnation (Total) Pressure */
          u = ug[i][j].b/ ug[i][j].a;
          v = ug[i][j].c / ug[i][j].a;
          vel = Math.sqrt(u*u + v*v);
          temp = ug[i][j].d/ug[i][j].a - (0.5 * (u*u + v*v));
          temp = temp/Cv;
          c = Math.sqrt(gamma * rgas * temp);
          mach = vel/c;
          mach = (1.0 + (gamma-1)*0.5*mach*mach);
          pres = ug[i][j].d/ug[i][j].a - (0.5 * (u*u + v*v));
          pres = pres/Cv;
          pres = pres * rgas * ug[i][j].a;
          plotdata[n][i-1] = pres * Math.pow(mach, gamma/(gamma-1.0));
          plotmax = Math.max(plotmax, plotdata[n][i-1]);
          plotmin = Math.min(plotmin, plotdata[n][i-1]);
          break;
        default:
          break;
        }  /* ends the switch */
      }
    }

    /* Now draw it up! */
    /* Start with Axes */
    g.setColor(Color.white);
    if (plotmin > 0 && plotmax > 0) {
      plotmin = 0;
      plotmax *=1.20;
    }
    else if (plotmin < 0 && plotmax < 0) {
      plotmax = 0;
      plotmin *= 1.20;
    }

    /* Draw x axis */
    y = gridheight + (int) ( (plotmin / (plotmax - plotmin))
                      * (double) plotheight)-ymargin;
    g.drawLine(xmargin, y, xmargin+plotwidth, y);
    /* Draw y axis */
    g.drawLine(xmargin, gridheight-ymargin,
               xmargin, gridheight-ymargin-plotheight);
    scalex = (double) plotwidth/(xmax-xmin);
    scaley = (double) plotheight/(plotmax-plotmin);
    xoffset = (int) (xmin * scalex);
    yoffset = (int) (plotmin * scaley);

    if (ndata == 3) {
      g.setColor(new Color(128, 128, 128));
      y = gridheight -ymargin - (int)(1.0*scaley)+yoffset;
      g.drawLine(xmargin, y, xmargin+plotwidth, y);
      g.setColor(Color.white);
      g.drawString("1.0", 5, y);
    }

    StringBuffer writeBuffer = new StringBuffer().append(plotmin);
    writeBuffer.setLength(6);
    g.drawString(writeBuffer.toString(), 5, gridheight-ymargin);
    g.drawString(String.valueOf(xmin), xmargin, gridheight-5);
    StringBuffer writeBuffer2 = new StringBuffer().append(plotmax);
    writeBuffer2.setLength(6);
    g.drawString(writeBuffer2.toString(), 5, gridheight-plotheight-ymargin+5);
    g.drawString(String.valueOf(xmax), plotwidth+xmargin-10, gridheight-5);
    x1 = gridwidth/2 - 10;
    g.drawString("Axial Position(x)", x1, gridheight-5);
//    g.drawString("Red: Top Wall", xmargin+15, gridheight-(ymargin+25));
    g.drawString("Red: Top Wall", xmargin+15, gridheight-(ymargin+800));
    g.drawString("Blue: Centerline", xmargin+15, gridheight-(ymargin+770));
    g.drawString("Green: Lower Wall", xmargin+15, gridheight-(ymargin+740));

    g.setColor(Color.green);
    for (i = 0; i < imax-1; ++i) {
      x1 = (int) (scalex * xnode[i][0]) - xoffset + xmargin;
      y1 = gridheight - ymargin - (int) (plotdata[0][i] * scaley) + yoffset;
      x2 = (int) (scalex * xnode[i+1][0]) - xoffset + xmargin;
      y2 = gridheight - ymargin - (int) (plotdata[0][i+1] * scaley) + yoffset;
      g.drawLine(x1, y1, x2, y2);
      g.drawOval(x1-3, y1-3, 6,6);
    }
    g.setColor(Color.blue);
    for (i = 0; i < imax-1; ++i) {
      x1 = (int) (scalex * xnode[i][jmid]) - xoffset + xmargin;
      y1 = gridheight - ymargin - (int) (plotdata[1][i] * scaley) + yoffset;
      x2 = (int) (scalex * xnode[i+1][jmid]) - xoffset + xmargin;
      y2 = gridheight - ymargin - (int) (plotdata[1][i+1] * scaley) + yoffset;
      g.drawLine(x1, y1, x2, y2);
      g.drawOval(x1-3, y1-3, 6,6);
    }
    g.setColor(Color.red);
    for (i = 0; i < imax-1; ++i) {
      x1 = (int) (scalex * xnode[i][jmax-1]) - xoffset + xmargin;
      y1 = gridheight - ymargin - (int) (plotdata[2][i] * scaley) + yoffset;
      x2 = (int) (scalex * xnode[i+1][jmax-1]) - xoffset + xmargin;
      y2 = gridheight - ymargin - (int) (plotdata[2][i+1] * scaley) + yoffset;
      g.drawLine(x1, y1, x2, y2);
      g.drawOval(x1-3, y1-3, 6,6);
    }    
  }

  void drawContour(Color ctable[], Graphics g, int ndata,
                   int width, int height) {

    /* Draws contours outlining data set specified by ndata */

    boolean drawcont=false;
        double c;  /* Speed of sound */
    double contval[];
    double data[][];
    
    double frac;
    double mach;
    double pres;
    double scale;
    double scrap;
    double u,       /* x-velocity */
          v,        /* y-velocity */
          vel,      /* total velocity */
          temp;     /* Temperature */
    double uprime, vprime;
    double xnode1, xnode2;
    double ynode1, ynode2;
    int i, j;
    int i1=0, i2=0, i3=0, i4=0;
    int j1=0, j2=0, j3=0, j4=0;
    int n;
    int gridheight,   /* height of drawing region for grid in pixels */
        gridwidth;    /* Width of same */
    int speciesindex;
    int xdraw1, xdraw2;
    int ydraw1, ydraw2;
    int xoffset, yoffset;

    /* Calculate the data you want from the state vector */

    datamin = Double.MAX_VALUE;
    datamax = Double.MIN_VALUE;
    data = new double[imax][jmax];
    for (i = 1; i < imax+1; ++i)
      for (j = 1; j < jmax+1; ++j) {
        switch(ndata){
        case 0:           /* Density */
          data[i-1][j-1] = ug[i][j].a;
          datamax = Math.max(datamax, data[i-1][j-1]);
          datamin = Math.min(datamin, data[i-1][j-1]);
        break;
        case 1:          /* Pressure */
          u = ug[i][j].b/ug[i][j].a;
          v = ug[i][j].c/ug[i][j].a;
          data[i-1][j-1] = ug[i][j].d/ug[i][j].a - (0.5 * (u*u + v*v));
          data[i-1][j-1] = data[i-1][j-1]/Cv;
          data[i-1][j-1] = data[i-1][j-1] * rgas * ug[i][j].a;
          datamax = Math.max(datamax, data[i-1][j-1]);
          datamin = Math.min(datamin, data[i-1][j-1]);
          break;
        case 2:         /* Temperature */
          u = ug[i][j].b/ug[i][j].a;
          v = ug[i][j].c/ug[i][j].a;
          data[i-1][j-1] = ug[i][j].d/ug[i][j].a - (0.5 * (u*u + v*v));
          data[i-1][j-1] = data[i-1][j-1]/Cv;
          datamax = Math.max(datamax, data[i-1][j-1]);
          datamin = Math.min(datamin, data[i-1][j-1]);
          break;
        case 3:        /* Mach Number */
          u = ug[i][j].b/ ug[i][j].a;
          v = ug[i][j].c / ug[i][j].a;
          vel = Math.sqrt(u*u + v*v);
          temp = ug[i][j].d/ug[i][j].a - (0.5 * (u*u + v*v));
          temp = temp/Cv;
          c = Math.sqrt(gamma * rgas * temp);
          data[i-1][j-1] = vel / c;
          datamax = Math.max(datamax, data[i-1][j-1]);
          datamin = Math.min(datamin, data[i-1][j-1]);
        break;
        case 4:        /* Axial Velocity */
          data[i-1][j-1] = ug[i][j].b / ug[i][j].a;
          datamax = Math.max(datamax, data[i-1][j-1]);
          datamin = Math.min(datamin, data[i-1][j-1]);
          break;
        case 5:        /* Normal Velocity */
          data[i-1][j-1] = ug[i][j].c / ug[i][j].a;
          datamax = Math.max(datamax, data[i-1][j-1]);
          datamin = Math.min(datamin, data[i-1][j-1]);
          break;
        case 6:        /* Stagnation (Total) Pressure */
          u = ug[i][j].b/ ug[i][j].a;
          v = ug[i][j].c / ug[i][j].a;
          vel = Math.sqrt(u*u + v*v);
          temp = ug[i][j].d/ug[i][j].a - (0.5 * (u*u + v*v));
          temp = temp/Cv;
          c = Math.sqrt(gamma * rgas * temp);
          mach = vel/c;
          mach = (1.0 + (gamma-1)*0.5*mach*mach);
          pres = ug[i][j].d/ug[i][j].a - (0.5 * (u*u + v*v));
          pres = pres/Cv;
          pres = pres * rgas * ug[i][j].a;
          data[i-1][j-1] = pres * Math.pow(mach, gamma/(gamma-1.0));
          datamax = Math.max(datamax, data[i-1][j-1]);
          datamin = Math.min(datamin, data[i-1][j-1]);
          break;
        default:

          /* Needs error checking here too?? */
          break;
        }
      }

    /* Determine contour values */

    contval = new double[ncontour];
    scrap = datamax-datamin;
    for (n = 0; n < ncontour; ++n)
      contval[n] = datamin + ((double)n/(double)(ncontour-1))*scrap;
    /* Set screen values */

    gridheight = height;  /* Leave room at bottom for labels */
    gridwidth = width;
    scale = Math.min(gridwidth/(xmax-xmin), gridheight/(ymax-ymin));
    xoffset = (int) (xmin * scale);
    yoffset = (int) (ymin * scale);
    /* Go through cells looking for crossings */
    for (i = 0; i < imax-1; ++i)
      for (j = 0; j < jmax-1; ++j) {
        for (n = 0; n < ncontour; ++n) {
          
          /* First check the SW triangle */

          drawcont = false;
          if ((contval[n] > data[i][j] && contval[n] < data[i+1][j]) ||
              (contval[n] < data[i][j] && contval[n] > data[i+1][j])) {
            i1 = i; j1 = j; i2 = i+1; j2  = j; drawcont = true; 
          }
          if ((contval[n] > data[i][j] &&  contval[n] < data[i][j+1]) ||
              (contval[n] < data[i][j] &&  contval[n] > data[i][j+1]))
            if(drawcont) {
              i3=i; j3=j; i4 = i; j4 = j+1;
            }
            else {
              i1 = i; j1 = j; i2 = i; j2 = j+1; drawcont = true;
            }
          if ((contval[n] > data[i][j+1] &&
               contval[n] < data[i+1][j]) ||
              (contval[n] < data[i][j+1] && 
               contval[n] > data[i+1][j])) {
            i3 = i; j3 = j+1; i4 = i+1; j4 = j;
          }
          
          if(drawcont) {  /* Go ahead, draw the contour */

            /* Sets color of contour */

            g.setColor(ctable[n]);
            frac = Math.abs((contval[n]-data[i1][j1])/
                            (data[i1][j1]-data[i2][j2]));
            xnode1 = (xnode[i2][j2]-xnode[i1][j1]) * frac + xnode[i1][j1];
            xdraw1 = (int) (scale * xnode1) - xoffset;
            ynode1 = (ynode[i2][j2]-ynode[i1][j1]) * frac + ynode[i1][j1];
            ydraw1 = gridheight - (int) (scale * ynode1) - yoffset;

            frac = Math.abs((contval[n]-data[i3][j3])/
                              (data[i3][j3]-data[i4][j4]));
            xnode2 = (xnode[i4][j4]-xnode[i3][j3]) * frac + xnode[i3][j3];
            xdraw2 = (int) (scale * xnode2) - xoffset;
            ynode2 = (ynode[i4][j4]-ynode[i3][j3]) * frac + ynode[i3][j3];
            ydraw2 = gridheight - (int) (scale * ynode2) - yoffset;
            g.drawLine(xdraw1, ydraw1, xdraw2, ydraw2);     
          }
         
          /* Then check the NE triangle */

          drawcont = false;
          if ((contval[n] > data[i][j+1] &&
               contval[n] < data[i+1][j+1]) ||
              (contval[n] < data[i][j+1] &&
               contval[n] > data[i+1][j+1])) {
            i1 = i; j1 = j+1; i2 = i+1; j2  = j+1; drawcont = true; 
          }
          if ((contval[n] > data[i+1][j] &&  contval[n] < data[i+1][j+1]) ||
              (contval[n] < data[i+1][j] && contval[n] > data[i+1][j+1]))
            if(drawcont) {
              i3=i+1; j3=j; i4 = i+1; j4 = j+1;
             }
            else {
              i1 = i+1; j1 = j; i2 = i+1; j2 = j+1; drawcont = true;
            }
          if ((contval[n] > data[i][j+1] && contval[n] < data[i+1][j]) ||
              (contval[n] < data[i][j+1] && contval[n] > data[i+1][j])) {
            i3 = i; j3 = j+1; i4 = i+1; j4 = j;
          }
          
          if(drawcont) {  /* Go ahead, draw the contour */
            g.setColor(ctable[n]);
            frac = Math.abs((contval[n]-data[i1][j1])/
                            (data[i1][j1]-data[i2][j2]));
            xnode1 = (xnode[i2][j2]-xnode[i1][j1]) * frac + xnode[i1][j1];
            xdraw1 = (int) (scale * xnode1) - xoffset;
            ynode1 = (ynode[i2][j2]-ynode[i1][j1]) * frac + ynode[i1][j1];
            ydraw1 = gridheight - (int) (scale * ynode1) - yoffset;

            frac = Math.abs((contval[n]-data[i3][j3])/
                              (data[i3][j3]-data[i4][j4]));
            xnode2 = (xnode[i4][j4]-xnode[i3][j3]) * frac + xnode[i3][j3];
            xdraw2 = (int) (scale * xnode2) - xoffset;
            ynode2 = (ynode[i4][j4]-ynode[i3][j3]) * frac + ynode[i3][j3];
            ydraw2 = gridheight - (int) (scale * ynode2) - yoffset;
            g.drawLine(xdraw1, ydraw1, xdraw2, ydraw2);     
          }
        }  /* Ends "for (n = 0; n < ncoutour...)" */
      }
  }
}

class Statevector
{
  double a;   /* Storage for Statevectors */
  double b;
  double c;
  double d;
  
  Statevector() {
    a = 0.0;
    b = 0.0;
    c = 0.0;
    d = 0.0;
  }

  /* Most of these vector manipulation routines are not          */
  /* used in this program because I inlined them for speed.      */
  /* I leave them here because they may be useful in the future. */

  public Statevector amvect(double m, Statevector that) {

    /* Adds statevectors multiplies the sum by scalar m */

     Statevector answer = new Statevector();

    answer.a = m * (this.a + that.a);
    answer.b = m * (this.b + that.b);
    answer.c = m * (this.c + that.c);
    answer.d = m * (this.d + that.d);
    return answer;
  }

  public Statevector avect(Statevector that) {
     Statevector answer = new Statevector();

    /* Adds two statevectors */

    answer.a = this.a + that.a;
    answer.b = this.b + that.b;
    answer.c = this.c + that.c;
    answer.d = this.d + that.d;
    return answer;
  }

  public Statevector mvect(double m) {
     Statevector answer = new Statevector();

    /* Multiplies statevector scalar m */    

    answer.a = m * this.a;
    answer.b = m * this.b;
    answer.c = m * this.c;
    answer.d = m * this.d;
    return answer;
  }

  public Statevector svect(Statevector that) {
     Statevector answer = new Statevector();

    /* Subtracts vector that from this */

    answer.a = this.a - that.a;
    answer.b = this.b - that.b;
    answer.c = this.c - that.c;
    answer.d = this.d - that.d;
    return answer;
  }

  public Statevector smvect(double m, Statevector that) {
     Statevector answer = new Statevector();

    /* Subtracts statevector that from this and multiplies the */
    /* result by scalar m */
    answer.a = m * (this.a - that.a);
    answer.b = m * (this.b - that.b);
    answer.c = m * (this.c - that.c);
    answer.d = m * (this.d - that.d);
    return answer;
  }
}

class Vector2
{
  double ihat;   /* Storage for 2-D vector */
  double jhat;

  Vector2() {
    ihat = 0.0;
    jhat = 0.0;
  }

  public double magnitude() {
    double mag;

    mag = Math.sqrt(this.ihat*this.ihat + this.jhat * this.jhat);
    return mag;
  }

  public double dot(Vector2 that) {
    /* Calculates dot product of two 2-d vector */
    double answer;

    answer = this.ihat * that.ihat + this.jhat * that.jhat;
    return answer;
  }
}

class DataPanel extends Panel {
  boolean gridSwitch = false;
  boolean tuftsSwitch = false;
  Color ctable[];
  datafield data2D;
  public double machff;
  public double fourthOrderDamping;
  public double secondOrderDamping;
  public int nplot;
  public int ntime;
  public int nview;
  int itercount;
  final int iterationsPerCycle = 1;  /* Number of iterations per redraw */
  final int keyWidth=100;            /* width of color key in pixels */

  /* For the screen buffering */

  private Image offScreenImage;
  private Graphics offScreenGraphics;
  private Dimension offScreenSize;

  DataPanel(InputStream instream) { /* Read in the input file */
    itercount = 0;

    try {
      data2D = new datafield(instream);
   } catch (IOException e) {
      System.err.println("  Input/output error: probably a bad filename or bad connection." + e.getMessage());
    }
    
    machff = data2D.machfftarget;
    fourthOrderDamping = 1.0;
    secondOrderDamping = 1.0;

    /* Initialize Color Table */

    ctable = new Color[16];
    ctable[0] = new Color(0,0,255);
    ctable[1] = new Color(0,64,255);
    ctable[2] = new Color(0,128,255);
    ctable[3] = new Color(0,192,255);
    ctable[4] = new Color(0,255,255);
    ctable[5] = new Color(0,255,192);
    ctable[6] = new Color(0,255,128);
    ctable[7] = new Color(0,255,64);
    ctable[8] = new Color(0,255,0);
    ctable[9] = new Color(32,255,0);
    ctable[10] = new Color(96,255,0);
    ctable[11] = new Color(160,255,0);
    ctable[12] = new Color(255,255,0);
    ctable[13] = new Color(255,160,0);
    ctable[14] = new Color(255,96,0);
    ctable[15] = new Color(255,0,0);
  }

  public void doIteration() {
    int n;
    itercount++;
/**************************** bcast here ****************************/
    double element[] = new double[3];
    element[0] = secondOrderDamping;
    element[1] = fourthOrderDamping;
    element[2] = ntime;

    MPI.COMM_WORLD.Bcast(element,0,3,MPI.DOUBLE,0);

    secondOrderDamping = element[0];
    fourthOrderDamping = element[1];
    ntime = (int)element[2];

    for (n = 0; n < iterationsPerCycle; ++n) {
      data2D.machfftarget = machff;
      data2D.secondOrderDamping = secondOrderDamping;
      data2D.fourthOrderDamping = fourthOrderDamping;
      data2D.ntime = ntime;
      data2D.doIteration();
    }
    repaint();
  }

  public void setGridSwitch(boolean nchoice) {
    gridSwitch = nchoice;
  }

  public void setTuftsSwitch(boolean nchoice) {
    tuftsSwitch = nchoice;
  }
  
  public void paint(Graphics g) {
    if (itercount == 0)
      return;
  
    g.setColor(Color.black);
    if (nplot == 0) {
      g.fillRect(0, 0, getSize().width-keyWidth, getSize().height - 1);
      if (gridSwitch)
        data2D.drawGrid(g, getSize().width-keyWidth, getSize().height-1);
      if (tuftsSwitch)
        data2D.drawTufts(g, getSize().width-keyWidth, getSize().height-1);
      data2D.drawContour(ctable, g, nview, getSize().width-keyWidth,
                         getSize().height-1);
      data2D.drawEdge(g, getSize().width-keyWidth, getSize().height-1);
      data2D.drawKey(ctable, g, keyWidth, getSize().height-1, 
                     getSize().width-1, getSize().height-1);
    }
    else {
      g.fillRect(0, 0, getSize().width-1, getSize().height - 1);
      data2D.drawLinePlot(g, nview, getSize().width-1, getSize().height-1);
    }
  }

  public final synchronized void update (Graphics theG)

    /* This section is based on  a borrowed  (Matt Gray - MIT)  */
    /* NoFlickerApplet. (impliments Double Buffering) */ 

    {
      Dimension d = getSize();
      if((offScreenImage == null) || (d.width != offScreenSize.width) ||
         (d.height != offScreenSize.height)) 
        {
          offScreenImage = createImage(d.width, d.height);
          offScreenSize = d;
          offScreenGraphics = offScreenImage.getGraphics();
          offScreenGraphics.setFont(getFont());
        }
    offScreenGraphics.setColor(Color.black);
    if (nplot == 0) {
      offScreenGraphics.fillRect(0,0,getSize().width-keyWidth, getSize().height - 1);
      if (gridSwitch){
        data2D.drawGrid(offScreenGraphics, getSize().width-keyWidth, 
                        getSize().height-1);
      }
      if (tuftsSwitch){
        data2D.drawTufts(offScreenGraphics, getSize().width-keyWidth,
                         getSize().height-1);
      }
      data2D.drawContour(ctable, offScreenGraphics, nview,
                         getSize().width-keyWidth,
                         getSize().height-1);
      data2D.drawEdge(offScreenGraphics, getSize().width-keyWidth,
                      getSize().height-1);
      data2D.drawKey(ctable, offScreenGraphics, keyWidth, getSize().height-1, 
                     getSize().width-1, getSize().height-1);
    }
    else {
      offScreenGraphics.fillRect(0, 0, getSize().width-1, getSize().height - 1);
      data2D.drawLinePlot(offScreenGraphics, nview, getSize().width-1,
                          getSize().height-1);
    }
            theG.drawImage(offScreenImage, 0, 0, null);
    }
}

class ParameterPanel extends Panel implements AdjustmentListener,ActionListener {
  final int sliderMin = 0;
  final int sliderMax = 100;
  double sliderValMin;
  double sliderValue;
  double sliderValMax;
  Scrollbar slider;
  TextField textField;

  double getValue() {
    double f;
    try {
      f = Double.valueOf(textField.getText()).doubleValue();
    } catch (java.lang.NumberFormatException e) {
      f = sliderValMin;
    }
    return f;
  }

  void makeSlidePanel(String parameterTitle) {

    /* This lays out components, sets initial appearances */

    GridBagConstraints c = new GridBagConstraints();
    GridBagLayout  gridbag = new GridBagLayout();
    setLayout(gridbag);

    c.fill = GridBagConstraints.HORIZONTAL;

    Label label = new Label (parameterTitle, Label.LEFT);
    c.weightx = 0.0;
    gridbag.setConstraints(label, c);
    add(label);

    textField = new TextField(String.valueOf(sliderValue), 6);
    c.weightx = 1.0;
    c.gridwidth = GridBagConstraints.REMAINDER;
    gridbag.setConstraints(textField, c);
    add(textField);
    textField.addActionListener(this);

    int sliderIntValue = 
      (int) ((sliderValue-sliderValMin)/(sliderValMax-sliderValMin) *
             (double) (sliderMax-sliderMin)) + sliderMin;
    slider = new Scrollbar(Scrollbar.HORIZONTAL, sliderIntValue, 10,
                           sliderMin, sliderMax);
    c.weightx = 0.0;
    c.gridheight = 1;
    c.gridwidth = GridBagConstraints.REMAINDER;
    gridbag.setConstraints(slider, c);
    add(slider);
    slider.addAdjustmentListener(this);
  }

  void setValue(int n) {
    setSliderValue(n);
    textField.setText(String.valueOf(sliderValue));
  }

  void setSliderValue(double f) {
    sliderValue = f;
    int sliderIntValue;
    sliderIntValue = 
      (int) ((f-sliderValMin) / (sliderValMax-sliderValMin) 
             * (double) (sliderMax-sliderMin)) + sliderMin;
    if (sliderIntValue > sliderMax) {
      sliderIntValue = sliderMax;
      sliderValue = sliderValMax;
    }
    if (sliderIntValue < sliderMin) {
      sliderIntValue = sliderMin;
      sliderValue = sliderValMin;
    }
    slider.setValue(sliderIntValue);
    textField.setText(String.valueOf(sliderValue));
  }

  public void adjustmentValueChanged(AdjustmentEvent e) {
      sliderValue = ((double) slider.getValue()/
        (double) (sliderMax-sliderMin)) * 
        (sliderValMax-sliderValMin) + sliderValMin;
      textField.setText(String.valueOf((float)sliderValue));
  }

  public void actionPerformed (ActionEvent e) {
     setSliderValue(getValue());
  }
}

class MachffPanel extends ParameterPanel  {

  /* Ok, here's where we handle events for this scrollbar */

  MachffPanel(String parameterTitle, double min, double max, 
              double initialMachff) {
    super();
    sliderValMin = min;
    sliderValMax = max;
    sliderValue = initialMachff;
    makeSlidePanel(parameterTitle);    
  }

  public void adjustmentValueChanged(AdjustmentEvent e) {
    double temp;

    /* Overrides ParameterPanel's default behavior */

      sliderValue = ((double) slider.getValue()/
        (double) (sliderMax-sliderMin)) * 
        (sliderValMax-sliderValMin) + sliderValMin;
      if (sliderValue > 0.9 && sliderValue < 1.3)
        sliderValue = 1.3;
      textField.setText(String.valueOf((float)sliderValue));
  }

  public void ActionPerformed (ActionEvent e){
    double temp = getValue();
      if (temp > 0.9 && temp < 1.3) {
        temp = 1.3;
        textField.setText(String.valueOf(temp));
      }
      setSliderValue(temp);  
  }
}

class DampingPanel extends ParameterPanel {

  /* Ok, here's where we handle events for this scrollbar */

  DampingPanel(String parameterTitle, double min, double max, 
              double initialDamping) {
    super();
    sliderValMin = min;
    sliderValMax = max;
    sliderValue = initialDamping;
    makeSlidePanel(parameterTitle);    
  }
}
\end{verbatim}

\subsection{C -- MPI}

\subsubsection{The problem}
Computing the value of $\pi$ using Simpson's Rule.

\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 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 usually
not as accurate as those obtained from $\arctan$.

\begin{verbatim}

/*--------------------------------------------*/
/* Saleh Elmohamed  1996.                     */
/* saleh@npac.syr.edu.                        */
/* Example for CPS-615.                       */
/*--------------------------------------------*/

#include  <stdio.h>
#include  <math.h>
#include  <mpi.h>

#define approx_val  2.19328059

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

double 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, l, 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. */
  l = 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, l, 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 * l/ 3.0;
      printf("-------------------------------------------------\n");
      printf("The computed Pi of the integral for %d grid points is  %25.16e \n",
             (N * Procs), pi);

/* This is directly derived from the integeration of the formula. See 
     the report. */
#if 1
      exact_val_of_Pi = 4.0 * atan(1.0);
#endif

#if 0
      exact_val_of_Pi = 4.0 * log(approx_val); 
#endif

      printf("The error or the discrepancy between exact and computed value of Pi : %25.16e\n",
             fabs(pi - exact_val_of_Pi));
      printf("-------------------------------------------------\n");
    }
  MPI_Finalize();
}

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

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

