
\section{\bf Introduction}

% Tom's HPF Intro.

High Performance Fortran \cite{HPFF} has been accepted as 
%the industry 
the standard
data parallel extensions to  Fortran 90 \cite{F90Book}.
 Its aim is to simplify
the process of developing data parallel applications on distributed memory
systems. HPF programs are expected to be scalable and portable, 
as the performance is preserved while moving them to different
platforms with a comparable number of processors. The relative
simplicity of developing codes in HPF comes from the fact that HPF
provides a  global name space as
sequential codes do. The tedious and error-prone task of
interprocessor communication detection and generation of appropriate
calls to the runtime system are left to the compiler. This makes HPF
codes easier to understand, develop, debug, and maintain.

HPF language features allow the user to control two basic aspects
of parallel computations: data distributions and parallelism. The user
specifies the data distributions using compiler directives. The
parallelism is expressed explicitly using Fortran 90 array syntax
augmented by the {\tt FORALL} construct, the {\tt INDEPENDENT} directive, and a
rich set of new intrinsic and library functions. The HPF is a new
technology, and only a few compilers are available at this time. We
question  in this paper  whether this technology actually
fulfills its promises. We selected a nontrivial,
well-known application, the Particle-in-Cell (PIC) plasma simulation code,
and implemented it in HPF.
Particle-in-cell plasma simulation codes are large-scale application codes
that comprise thousands of lines of code and have large memory and
high-performance
requirements. The numerical simulation of plasmas is one of the Grand Challenge
problems facing the high-performance computing community.
These characteristics of PIC codes provide a strong motivation for their
 parallelization.
%Large  memory and high computational demands of the PIC codes create a
%strong motivation for their parallelization. 
However, data accesses with multiple levels of 
indirection and complicated message-passing patterns in some parts of the
parallel PIC  codes make them hard to be implemented  in HPF. 

The PIC method is an alternative  to the direct particle-particle
methods where simulation includes the interaction of all the
particles in the simulation domain, and the complexity is O($N_p^2$) for a
system of $N_p$ particles. Instead, the  PIC method simulates the interaction
of particles with the fields resulting from the particles' charges.
The simulation domain is partitioned into cells separated by grid points. 
Each cell contains a number of particles.
PIC method is especially
preferable for simulating systems containing an order of $10^6$ or
more particles, since its computational complexity is limited to O($N_p$) for 
particles and O($N_g.log N_g$) for the grid points where $N_g$ is the number 
of grid points in the simulation domain.
In the PIC method, moving the particles to a new position
dominates the runtime of the algorithm and this phase can be parallelized
easily due to the independent behavior of particles.

The basic PIC algorithm consists of an initialization phase followed by
four processing phases that are 
repeated many times (Figure~\ref{fig:pic-phases}).
 In the {\it scatter} phase, the particle
attributes are interpolated to nearby points of a regular computational
mesh that discretizes the problem domain. The appropriate field equations
are then solved on the computational mesh ({\it field solve} phase), and 
the force on each particle due to resultant fields is found by
interpolation on the mesh ({\it gather} phase). Finally, the particles are 
repositioned under the influence of this force in the {\it particle push}
phase.

Parallel PIC methods have been previously studied by a number of
researchers on a wide variety 
of platforms using different parallel programming methodologies 
\cite{OO-PIC, PAX-PIC, CM-PIC, Conc-PIC, Cube-PIC, Decyk-PIC}.
%Native message passing mechanisms of the
%experimental parallel machines were used in many of the previous studies.

We  implemented one- and  two-dimensional
PIC simulation codes in HPF. We compiled our codes using PGI's HPF
compiler and ran them on an IBM SP-2 with up to 32 nodes. 

In this paper we talk about our data decomposition and parallelization
strategy,  and emphasize HPF features used in individual phases of the 
simulation code. We assess the power
of HPF in representing such complicated codes. 
For the sake of simplicity, the most of the following discussion involves only
the one-dimensional PIC method; the interested reader may refer to 
\cite{PIC-Book1, PIC-Book2}
for higher dimensional PIC methods.

\begin{figure}[htb]
\centerline{\psfig{figure=xfig/pic-phases.eps,width=4.00in}}
\caption{\sf Main phases of a typical PIC code.}
\label{fig:pic-phases}
\end{figure}

%===================
% End of intro.tex
%===================

\section{The Approach to Parallelization}

We applied a top-down approach when trying to identify all the high-level
parallelism afforded by the problem before moving down to other levels.
We used the most efficient parallel algorithms at each step and
analyzed the implementation tradeoffs to determine what could be done
to make the results more attractive in the scope of HPF. 

We will explore our approach to the problem's parallelization using HPF
in five levels: distribution of principal data structures, loop level
parallelization, communication optimizations,
conveying data distributions to the subroutines, and augmenting HPF
with extrinsic routines when necessary. 

