
\subsection{Field-Solve Phase}
% IFFT -> Poisson -> FFT
 
Field solve phase consists of three steps. First, the inverse FFT of complex
charge density values are taken, then the electric field values are determined
by solving  Poisson's equations  on the complex domain. Last, the forward FFT
of the results is taken.

\subsubsection{Forward and Inverse Fast Fourier Transform (FFT)}
% \bf SUBROUTINE PFFT1R}

\begin{figure}[htbp]
\centerline{\psfig{figure=xfig/fft1.eps,width=3.00in}}
\caption{\sf Schematic description of FFT.}
\label{fig:fft}
\end{figure}

FFT has two major steps. The first step is to permute the input data. 
Element {\sc in(i)} goes to the element {\sc out(reverse(i))} of a 
temporary array. The {\sc reverse()} function performs a bit-reversal 
operation. We divide the input data among processors
evenly. The permutation step requires an all-to-all 
broadcast communication to send each element to the appropriate place. 
In the second step, the FFT of input data is calculated in $log(N)$ steps.
At each step, elements between processors are interchanged 
(Figure~\ref{fig:fft}) for an 8-element array on an 8-processor system. 
At each step $i$, the destination is found by flipping the $(i-1)^{th}$ bit of 
the original global indices. If that element is not located locally, then
it is fetched from its owner processor. 


The FFT code is too complicated to be handled efficiently by  current HPF 
compilers. We wrote a message-passing version of FFT using extrinsic functions.
Because each extrinsic function running on different processors 
sees only its own local data,whenever it needs other processors' data, it 
requests and gets needed data. Here, all processors know  when and to which 
processor it will send data and receive data from, as explained above.

A scrambling routine is executed  after the forward FFT and 
before the inverse FFT. 
Each processor exchanges its charge density array with its pair processor, 
as described in Figure~\ref{fig:scramble}, for eight processors and does some 
reduction operation on local and received data. 

\begin{figure}[htbp]
\centerline{\psfig{figure=xfig/scramble1.eps,width=4.0in}}
\caption{\sf Communication structure of scrambling.}
\label{fig:scramble}
\end{figure}

\subsubsection{Poisson Solver}
% PPOIS Routine

The one-dimensional Poisson equation in the Fourier space is solved to find the
force per unit charge {\sc fc(k)} at each FP grid point {\sc k} using periodic 
boundary conditions. 
The force per unit charge at each particle position is calculated
independently using an {\tt INDEPENDENT DO} loop (Figure~\ref{fig:POIS1}.)
A complex form factor array {\sc ffc} that is prepared in the
initialization phase is used in the computation. 

\SS
\begin{figure}[htbp]
\begin{small}
\noindent\hrulefill
{\sf
\begin{tabbing}
xxxxxxxxxx\=xxxxx\=xxxxx\=xxxxx\=xxxxx\kill
%\>\> {\bf Subroutine Poisson\_1(q,f,isign,ffc,ax,affp,wp,nx) }\\
%\>\> Dimension q(nxpmx,nblok),f(nxpmx,nblok),ffc(kxp,kblok)\\
%\> !HPF\$ Distribute *(*,block) :: q,f,ffc \\
\>\> dnx = 2*$\pi$/float(nx) \\
\>  !HPF\$ INDEPENDENT,NEW(wp,koff,j,temp)\\
\>\>   DO k = 1, Number\_of\_Fpartitions \\
\>\>\>      koff = kxp*(k-1)-1 \\
\>\>\>      DO j = 1, Number\_of\_Complex\_Grid\_Points  \\
\>\>\>\>       temp = dnx*(j+koff)*real(ffc(j,k)*aimag(ffc(j,k))\\
\>\>\>\>       fc(j,k) = temp*cmplx(aimag(q(j,k)),-real(q(j,k)))\\
\>\>\>      END DO\\
\>\>   END DO
\end{tabbing}
}
\noindent\hrulefill
\end{small}
\caption{\sf Finding the  electrical fields at FP grid points by using a 
Poisson solver.}
\label{fig:POIS1}
\end{figure}
\DS

\subsection{Gather (Interpolation) Phase}

%\subsection {\bf Transferring Data Between Complex Charge Array And 
%      Real Charge Array}
% CFPPP1 Routine 

To copy data from a complex force array {\bf fc}
to a real force array {\bf f}
we reverse the operation of the scatter copy routine above. 
Guard cells are similarly updated.
If we reverse the direction of the arrows in Figure \ref{fig:cpp2fp},
we get a picture of how this copy operation is done. 

\SS
\begin{figure}[htb]
\begin{small}
\noindent\hrulefill
{\sf
\begin{tabbing}
xxxxxxxxxx\=xxxxx\=xxxxx\=xxxxx\=xxxxx\kill 
%\>\>{\bf Subroutine Copy\_PP\_2\_FP1(q,qc,isign)}\\
\>\> FORALL(k=1:kblok,j=1:kxp )\\
\>\>\>      f(2*j-1,k) = real(fc(j,k))  \\
\>\>\>      f(2*j,k) = aimag(fc(j,k))  \\
\>\>   END FORALL\\
\\
\> C Consider contributions to guard cells \\
\>\>   FORALL k = 1 : kblok \\
\>\>\>      f(1,k) = aimag(fc(kxp,k-1+NP (mod NP)))\\
\>\>\>      f(nxpmx-1,k) = real(fc(1,k+1(mod NP)))\\ 
\>\>\>      f(nxpmx,k) = aimag(fc(1,k+1(mod NP)))\\
\>\>   END FORALL
\end{tabbing}
}
\noindent\hrulefill
\end{small}
\caption{\sf Copying data from complex charge density array into real charge 
density array.}
\label{code:CPPFP1}
\end{figure}
\DS

 The forces on the boundaries should be updated at the second step, 
as seen in the second {\tt FORALL} loop of Figure~\ref{code:CPPFP1}. 
The compiler analyzes the data access patterns and, if necessary, shifts 
the values at 
boundary grid points internally  to the neighbor partitions exchanging data 
with other processors.


\subsection{Particle Push Phase}
% PPUSH1 Routine

Using  a leap-frog scheme in time and a second-order spline interpolation 
in space, the particles' positions and velocities are updated. 
Due to the periodic boundary conditions, any particle moving out of
the left-most boundary rejoins the simulation from the right-most boundary,
and vice versa.

The electrostatic force on each particle position is calculated
by interpolation of the electric fields on the  nearest grid points to the
given position. The interaction among the electrical fields on grid points
and  nearby particles is the reverse of that shown in 
Figure~\ref{fig:pdost}.
The electrostatic force for partition {\sc k} at local grid point 
{\sc j} is kept in {\sc f(j,k)}.

The particle push phase can be performed in HPF by using an 
{\tt INDEPENDENT DO} loop over all the particle partitions. 
Since all data items used in the loop body are local, there is no need for 
interprocessor communication. As the compiler  we used 
serialized this loop portion when it detected the dynamically
computed value of {\sc nn} in the body , we wrote extrinsic routines to 
simulate the effect of the {\tt INDEPENDENT DO} loop.
computed value of {\sc nn} in the body. 
Compilers assumes that the use of  {\sc nn}
to index  the {\sc f} array will cause a wide range of elements to be
accessed that may disturb inter-iteration independence, and it therefore 
completely  serializes  the loop. 

\SS
\begin{figure}[htbp]
\begin{small}
\noindent\hrulefill
{\sf
\begin{tabbing}
xxxxxxxxxx\=xxxxx\=xxxxx\=xxxxx\=xxxxx\kill
%$\> {\bf Subroutine Ppush1(part,f,ek) }\\
%\> Dimension part(idimp,npmax,nblok),f(nxpmx,nblok)\\
%!HPF\$ Distribute *(*,*,block) :: part \\
%!HPF\$ Distribute *(*,block)   :: f \\
!HPF\$ INDEPENDENT,NEW(j,nn,dx,dv)\\
\> DO k = (1 : Number\_of\_particle\_partitions) \\
\>\> DO j = 1, number\_of\_grid\_points \\
\>\>\>   nn = $\lfloor part(1,j,k)+0.5 \rfloor$  \\
\>\>\>   dx = part(1,j,k)-nn  \\
\>\>\>   dv = part(VEL,j,k) + ((...dx...f(nn,k)...f(nn+1,k)...f(nn-1,k)...)\\
\>\>\>   part(VEL,j,k) = dv  \\
\>\>\>   part(POS,j,k) = part(POS,j,k) + dv*dt \\
\>\>  END DO\\
\>   END DO
\end{tabbing}
}
\noindent\hrulefill
\end{small}
\caption{\sf Pushing  the particles.}
\label{code: PPUSH1}
\end{figure}
\DS

\subsubsection {\bf  Particle Migration Step}
% (PMOVE1)}

Once the new positions of particles have been computed, the ones that leave
the current partition should be passed to the 
processor responsible for their new position.
As noted before, the particles are placed into partitions according to their
positions.

The main data structures used in this phase are given in the following table:

\begin{center}
{\sf 
\begin{small}
\begin{tabular}{| c | l |} \hline
Data Structure            & Definition \\ \hline \hline
sbufr(idimp,nbmax,nblok)  & Send buffer to the right partition\\
sbufl(idimp,nbmax,nblok)  & Send buffer to the left  partition\\
rbufr(idimp,nbmax,nblok)  & Receive buffer to the right partition\\
rbufl(idimp,nbmax,nblok)  & Receive buffer to the left partition\\
jsl(idimp,nblok)          & Left-going  particle counter  \\
jsr(idimp,nblok)          & Right-going particle counter \\
jss(idimp,nblok)          & Scratch array \\ \hline
\end{tabular}
\end{small}
}
\end{center}

Here {\sc nbmax} is the size of the buffer for passing particles. {\sc idimp}
enumerates particle properties, such as velocity and position. {\sc numpp}
is the number of particle partitions. 

\SS
\begin{figure}[htbp]
\begin{small}
\noindent\hrulefill
{\sf
\begin{tabbing}
xxxxxxxxxx\=xxxxxxxx\=xxxxxxxx\=xxxxxxxx\=xxxxx\kill
%\>\>{\bf Subroutine Pmove1(part,rbufl,rbufr,sbufr,sbufl,jsl,jsr,jss,edges,npp)} \\
%\>\> Real,dimension(idimp,npmax,nblok) :: part \\
%\>\> Real,dimension(idimp,nbmax,nblok) :: rbufr,rubfl,sbufr,sbufl \\
%\>\> Real,dimension(idimp,nblok)       :: jsl,jsr,jss,edges\\
%\>\> Real,dimension(ntmax,nblok)       :: ihole\\
%\>\> Real,dimension(nblok)             :: npp\\
%\> !HPF\$ Distribute *(*,*,block)     :: part,rbufl,rbufr,sbufr,sbufl \\
%\> !HPF\$ Distribute *(*,block)       :: jsl,jsr,jss,edges,ihole\\
%\> !HPF\$ Distribute *(block)         :: npp\\
\>\> npr = SUM(npp)\\
\>10 \> call Pmove1\_part1(part,sbufr,sbufl,jsl,jsr,jss,ihole,edges,npp)\\
\>20 \>rbufl     \> = CSHIFT(sbufr, dim=3, RIGHT)  \\
\>\>   jsl(2, :) \> = CSHIFT(jsr(1, :), dim=2, RIGHT) \\
\>\>   rbufr     \>  = CSHIFT(sbufl, dim=3, LEFT)  \\
\>\>   jsr(2, :) \>  = CSHIFT(jsl(1, :), dim=2, LEFT) \\
\>\>   call Pmove1\_part2(rbufl,rbufr,jsl,jsr,edges,nps)\\  
\>\>   call Pmove1\_part3(rbufl,rbufr,sbufr,sbufl,jsl,jsr,jss,edges)\\
\>\>   call Pmove1\_part4(jsl,jsr,jss,npmax,ierr)\\
\>\> call Pmove1\_part5(part,rbufl,rbufr,jsl,jsr,jss,edges,npp,ihole)\\
\>\>  IF nps $>$ 0 GOTO 20\\
\>\>  IF ibflag $>$ 0 GOTO 10\\
\>\>  npss = SUM(npp)
\end{tabbing}
}
\noindent\hrulefill
\end{small}
\caption{\sf Moving the particles into appropriate spatial regions after the
Push phase.}
\label{code:PMOVE1}
\end{figure}
\DS

The main move routine is given in Figure~\ref{code:PMOVE1}. 
We add the number of particles, {\sc npr}, across processors in one step using
HPF's {\tt SUM} intrinsic. The
{\tt Pmove1\_part1()} routine is used to determine outgoing particles and 
place the properties of particles into the send buffers. This is easily done 
using an {\tt INDEPENDENT DO} loop (Figure~\ref{code:PMOVE1-1}).

\SS
\begin{figure}[htbp]
\begin{small}
\noindent\hrulefill
{\sf
\begin{tabbing}
xxxxxxxxxx\=xxxxx\=xxxxx\=xxxxx\=xxxxx\kill
%\>{\bf Subroutine Pmove1\_part1(part,sbufr,sbufl,jsl,jsr,jss,ihole,edges,npp,ibflag)} \\
%\>\> Real,dimension(idimp,npmax,nblok) :: part \\
%\>\> Real,dimension(idimp,nbmax,nblok) :: sbufr,subfl \\
%\>\> Real,dimension(idimp,nblok)       :: jsl,jsr,jss,edges\\
%\>\> Real,dimension(ntmax,nblok)       :: ihole\\
%\>\> Real,dimension(nblok)             :: npp\\
%\>!HPF\$ Distribute *(*,*,block)     :: part,sbufl,sbufr \\
%\>!HPF\$ Distribute *(*,block)       :: jsl,jsr,jss,edges,ihole\\
%\>!HPF\$ Distribute *(block)         :: npp\\
\>\>    jsl = jsr = jss = 0 \\
\>!HPF\$ INDEPENDENT , NEW(xt,j)\\
\>\>    DO k=1, Number\_of\_Ppartitions\\
\>\>\>   DO j = 1, Number\_of\_particles\_in\_partition \\
\>\>\>        xt = part(POS, j, k)  \\
\>\>\>        IF (xt $>$ edges(RGHT, k) THEN \\
\>\>\>\>         jsr(1,k) = jsr(1,k) + 1 \\
\>\>\>\>         If k = nvp then    ! If PP k is last PP \\
\>\>\>\>            xt = xt - nx  \\
\>\>\>\>         EndIf \\
\>\>\>\>         sbufr(POS, jsr(1,k),k) = xt   \\
\>\>\>\>         sbufr(VEL, jsr(1,k),k) = part(2,j,k)  \\
\>\>\>\>         ihole(jsl(1,k)+jsr(1,k),k) = j  \\
\>\>\>        END IF \\
\>\>\>       \\
\>\>\>    ... Symmetric operations for the particles going left ... \\
\\
\>\>\>  END DO  \\ 
\>\>\>  jss(1,k)=jsl(1,k)+jsr(1,k) !total number of outgoing particles \\
\>\> END DO 
\end{tabbing}
}
\noindent\hrulefill
\end{small}
\caption{\sf Determining the outgoing particles and placing them into
send buffers (Pmove1-1).}
\label{code:PMOVE1-1}
\end{figure}
\DS

Since the particles can move only one partition at each iteration of the 
loop (Figure~\ref{code:PMOVE1-1}), doing an HPF {\tt CSHIFT} on the last 
dimension of the send buffers and the send counter arrays transfers  all the 
information about the particles moving out of their partitions to the receive 
buffers and receive buffer counter
arrays, respectively, (Figure~\ref{fig:PMOVE1}). Because of the periodic 
boundaries, the last partition in the
simulation domain exchanges data with the first partition. 

The {\tt Pmove1\_part2()} routine is responsible for checking the receive 
buffers  to determine whether any particles need to be sent further. 
If there are such  particles, then 
{\tt Pmove1\_part3()} takes over the control and places 
those particles into send buffers
again while transferring the other particles into vacant positions in the receive buffers.

The overflow condition in any specific partition is checked by the 
{\tt Pmove1\_part4()} routine.
At the last step we place the  incoming particles  into particle array using 
{\tt Pmove1\_part5()}. If any gaps remain in the particle arrays, 
then the particles
toward the end of the array are shifted forward to fill in the gaps
(Figure~\ref{code:PMOVE1-5}).

These five steps are repeated until  each partition arrives at its final 
destination.
However, since many particles travel for only short distances, there are 
usually only a few  iterations.

\begin{figure}[htbp]
\centerline{\psfig{figure=xfig/move-fig.eps,width=4.0in}}
\caption{Communication of buffers between processors.}
\label{fig:PMOVE1}
\end{figure}


\SS
\begin{figure}[htbp]
\begin{small}
\noindent\hrulefill
{\sf
\begin{tabbing}
xxxxxxxxxx\=xxxxx\=xxxxx\=xxxxx\=xxxxx\kill
%\>\>{\bf Subroutine Pmove1\_part5(part,rbufl,rbufr,jsl,jsr,jss,edges,npp,ihole)} \\
%\>\> Real,dimension(idimp,npmax,nblok) :: part \\
%\>\> Real,dimension(idimp,nbmax,nblok) :: rbufr,rubfl \\
%\>\> Real,dimension(idimp,nblok)       :: jsl,jsr,jss,edges\\
%\>\> Real,dimension(ntmax,nblok)       :: ihole\\
%\>\> Real,dimension(nblok)             :: npp\\
%\> !HPF\$ Distribute *(*,*,block)     :: part,rbufl,rbufr \\
%\> !HPF\$ Distribute *(*,block)       :: jsl,jsr,jss,edges,ihole\\
%\> !HPF\$ Distribute *(block)         :: npp\\
\> !HPF\$ INDEPENDENT, NEW(j,offset) \\
\>\> DO k=1,number of partitions(kblok)\\  
\>\>   IF (jsl(IN, k) .le. jss(OUT,k)) then\\
\>\>   ! put into holes \\
\>\>\>   DO j = 1, jsl(IN, k) \\
\>\>\>\>    part(POS, ihole(j, k), k) =  rbufl(POS, j, k) \\
\>\>\>\>    part(VEL, ihole(j, k), k) =  rbufl(VEL, j, k)  \\
\>\>\>   END DO   \\
\>\>\>   ! check right ones \\
\>\>\>   jss(IN, k) = min(jss(OUT, k) - jsl(IN, k), jsr(IN,k)) \\
\>\>\>   DO j = 1, jss(IN, k) \\
\>\>\>\>    part(POS, ihole(j, k)+jsl(IN,k), k) =  rbufr(POS, j, k) \\
\>\>\>\>    part(VEL, ihole(j, k)+jsl(IN,k), k) =  rbufr(VEL, j, k) \\
\>\>   END DO \\
\>\>   IF (jsr(IN, k) .gt. jss(OUT, k) - jsl(IN, k)) THEN \\
\>\>\>    DO j = 1, jsr(IN, k)-jss(IN, k) \\
\>\>\>\>     part(POS, npp(k)+j, k) =  rbufr(POS, j+jss(IN,k), k) \\
\>\>\>\>     part(VEL, npp(k)+j, k) =  rbufr(VEL, j+jss(IN,k), k) \\
\>\>\>    END DO \\
\>\>   END IF \\
\>\>   ELSE \\
\>\>\>   jss(IN, k) = min(jsl(IN, k), jss(OUT, k)) \\
\>\>\>   DO j = 1, jss(IN, k) \\
\>\>\>\>    part(POS, ihole(j, k), k) =  rbufl(POS, j, k) \\
\>\>\>\>    part(VEL, ihole(j, k), k) =  rbufl(VEL, j, k)  \\
\>\>\>   END DO   \\
\>\>\>   offset = jsl(IN,k) - jss(IN,k) \\
\>\>\>   DO j = 1, offset \\
\>\>\>\>    part(POS, npp(k)+j, k) =  rbufl(POS, j+jss(IN,k), k) \\
\>\>\>\>    part(VEL, npp(k)+j, k) =  rbufl(VEL, j+jss(IN,k), k) \\
\>\>\>   END DO \\   
\>\>\>   DO j = 1, jsr(IN, k)-offset \\
\>\>\>\>     part(POS, npp(k)+offset+j, k) =  rbufr(POS, j, k) \\
\>\>\>\>     part(VEL, npp(k)+offset+j, k) =  rbufr(VEL, j, k) \\
\>\>\>   END DO \\
\>\>\>   IF (jss(OUT,k) - jsl(IN,k) - jsr(IN,k)) .gt. 0) THEN\\
\>\>\>\>       offset = jss(OUT,k) - jsl(IN,k) - jsr(IN,k) \\
\>\>\>\>       DO j = 1, offset \\
\>\>\>\>        part(POS, ihole(jsl(IN,k)+jsr(IN,k)), k) = part(POS, npp(k), k)\\
\>\>\>\>        part(VEL, ihole(jsl(IN,k)+jsr(IN,k)), k) = part(VEL, npp(k), k)\\
\>\>\>\>        npp(k) = npp(k) - 1 \\
\>\>\>\>       END DO \\
\>\>\>    ENDIF \\
\>\>    ENDIF \\
\>\>  ENDDO 
\end{tabbing}
}
\noindent\hrulefill
\end{small}
\caption{\sf Distributing of incoming particles from buffers (Pmove1-5).}
\label{code:PMOVE1-5}
\end{figure}
\DS

Each of these five steps is done in parallel. The bodies of each routine 
encompass a large
{\tt INDEPENDENT DO} loop. Since the compiler  we use has a tendency to 
serialize  these loops, we simulated them by writing extrinsic routines. 

Array indices with one or more levels of indirection 
such as {\sc part(pos, ihole(j, k), k) = rbufl(pos, j, k)}, 
causes the compiler to generate a serial loop, 
and bring in all the 
possible values of the right-hand side for each processor. Since accumulation 
in the {\tt INDEPENDENT
DO} loop is not allowed, we do this by calculating partial sums in 
each processor and then 
taking a global sum using the {\tt SUM} intrinsic across all processors. 
 










