\setcounter{equation}{0}
\setcounter{figure}{0}

\chapter{Case Studies of Multi-Paradigm Programming in the Software
Integration Framework}
\markright{APPENDIX B. CASE STUDIES IN CHAPTER 3}

This appendix details three application case studies used in
Chapter 3 to demonstrate the integration of multiple-programming
models in a unified dataflow-based programming framework.
The three studies are chosen from atmosphere science, 
computational electromagnetics, and financial modeling. For
each application, we describe a prototype system that illustrates 
system configuration, software integration,
performance, functionality, and other issues. Some background information
about the CM5 and PVM are also given to help explain the 
approaches we adopted when working in a specific context of
machine and programming tools.

\section{Parallel Programming on Connection Machine CM5}

A detailed description of the CM5 architecture and software system
can be found in \cite{cm5} and related CM5 documents from TMC. Here, 
we briefly outline the existing parallel programming environment
on the CM5 for the purpose of our discussion in later sections.

The CM5 system combines the features of SIMD and MIMD designs, integrating 
them into a single parallel architecture. The unique feature of the CM5 provides 
 alternatives with which to exploit both paradigms in a parallel 
programming environment. 

The CM5 architecture is designed to support data-parallel architecture
 especially well,
which implies relatively synchronous parallel programming 
featuring a single conceptual thread of control and a global 
approach to data layout and I/O
in which arrays of data are laid out across all processors with all
elements computed upon simultaneously by their respective
processors. Although the conceptual model itself is completely synchronous,
in practice processors may asynchronously execute any operations that are
independent of other processors (such as conditional operations
on locally stored data). The CM5 data-parallel compilers
take full responsibility and advantage of this
freedom to exploit MIMD hardware capabilities of the CM5. Currently,
the CM5 supports data-parallel Fortran (CMfortran), data-parallel C (C*), and
data-parallel Lisp (*Lisp) programming languages.

\begin{figure}
\centerline{\epsfxsize=420pt \epsfbox[50 600 550 720]{chapter3-multiparadigm/fig3.eps}}
%BoundingBox: 50 600 550 720
\caption{Components of a CM5 system}
\label{cm5}
\end{figure}

An explicit message-passing model is extended from the data-parallel model on the CM5
by using a package of macros and runtime primitives called CMMD supporting
MIMD style low-level communications operations. Programs that use
CMMD typically operate in an SPMD(single program multiple data) style in
that each processing node of the CM5 runs an independent copy of
a single program and manages its own computations and data layout with
communication between nodes handled by calls to CMMD message-passing primitives. The message-passing
programs can be written in Fortran 77, CMFortran, C, C++, C* or an 
assembly language DPEAC.

As shown in Figure \ref{cm5}, each processing node (PN) is a general-purpose computer 
that can fetch and interpret its own instruction stream, execute arithmetic
and logical instructions, calculate memory address, and perform
interprocessor communication. All PNs can perform independent tasks or
collaborate on a single problem. A control processor (CP) is similar to 
a standard
high-end workstation acting as a partition manager to communicate with the 
rest of CM5 system through the Control Network and Data Network.
CMOST, a parallel timesharing operating
system enhanced from UNIX,
runs on the CP to make all CM5 resource allocation and swapping
decisions, as well as most system calls for process execution, memory
management, I/O and access to local networks. 
The Control Network provides tightly coupled communication services. 
Optimized for fast response and low latency, its functions include
synchronizing the PNs, broadcasting data to every node, combining
a value from every node to produce a single result, and computing
certain parallel prefix operations. The Data Network provides
loosely coupled communication services. Optimized for high bandwidth,
it provides point-to-point data delivery.
A data-parallel or CMMD message-passing process running on the CP, plus 
the nodes, is a single process time-sharing the CM5 system 
with other processes.

The current CM programming environment provides some restricted ways to 
support the full integration of multi-paradigm programming in a single
language, mostly through specialized library routines
that encapsulate optimized implementation details of the functions provided. 
As an example, the CM Fortran Utility Library provides convenient 
access from CM Fortran 
to the capabilities of lower-level CM software. The purpose is 
typically to 
achieve functionality or performance beyond what is currently 
available from the compiler.
CM Fortran programmers can use the library procedures in situations where
one is normally tempted to make explicit calls to lower-level software.
The utility procedures take CM Fortran array names and other CM Fortran data 
objects as arguments so that there is no need to convert CM Fortran 
objects into data
types by lower-level software. The CM Fortran Utility Library provides 
procedures for
inquires, random number generation, dynamic array allocation, interprocessor and 
local data motion, and parallel I/O\cite{cmf}. 

In order to exploit the special high-performance hardware arithmetic 
accelerator (called vector unit)
in each PN of the CM5 system, current CM5 programming 
software provides some limited ways to mix data-parallel
paradigms in a message-passing program. Explicit message-passing CMMD 
programs can use certain language constructs from CMFortran or C*,
in which case the programs use both the microprocessor and 
the vector units of a PN, with
serial arrays stored in the microprocessor and parallel arrays in the VU memory.
It can thus take advantage of the VUs for high-speed floating-point
computations. Parallel arrays can be passed as arguments to CMMD message-passing
primitives in much the same way as serial arrays. Within this integration model, data-parallel
operations are confined in a ``node level'' while 
the overall program is in a explicit message-passing
paradigm.

If adopting this embedded integration approach on a CM5, 
an ideal integrated program would be as follows:
Executed as a single process on the control processor (CP) of the CM5, 
the main program in Fortran77 (or C) would
first act as a host sequential program in the 
CMMD host/node programming model and start a CMMD node program by
enabling and broadcasting some data to the node program, which is written in
Fortran77, C, or even CMfortran (node-level data-parallel). Then 
this host program would invoke another subroutine written
in CMFortran (or C*). After the control is back to the
main program from the CMFortran subroutine, it collects the results sent 
from the node CMMD program and finally ends. Before and after enabling the CMMD node
program or invoking the CMFortran subroutine, 
there can be any other non-parallel computing
subroutines in the sequential host program.
Technically, however, this approach would pose significant difficulties 
on the complier (linker), as it is now the compiler's responsibility, as 
opposed to that of a third-party daemon or the OS
in the process-based approach described below, to manage the time-sharing
of the CM5's PNs and internal networks between the CMFortran subroutine and 
the CMMD node program.


\section{Parallel Virtual Machine (PVM)}

The PVM \cite{pvm,pvm2} is a public-domain software package that 
supports utilization of heterogeneous
networks of parallel and serial computers as a single computational resource.
The PVM provides a general-purpose, message-passing interface to TCP/IP networking protocol,
portable across multiple platforms.

The PVM consists of two parts: a daemon process on each host machine of the virtual machine,
and a user library that contains primitive routines for initiating processes on other hosts to allow
communication and data transfer between processes and for configuration
of the virtual machine. The PVM assumes that each host in a virtual machine is
directly connected to each other host via IP. 
The daemons communicate with one another using UDP sockets.
TCP sockets are used between a daemon and the local PVM tasks, and also directly
between tasks on the same host or on different hosts if necessary. The message-passing
system is implemented on top of the UDP and TCP sockets.

The PVM software package has no built-in visualization capability.
To develop a graphical user interface (GUI) for a networking application program 
written in PVM, one has to link to other graphics packages such as Xlib
or X Toolkit. Most previous work in developing a GUI component in a PVM system
was focused on performance monitoring, graphical debugging, and
visual development tools for aiding the programmer in the tasks of
designing, compiling, executing, debugging, and 
analyzing a parallel program \cite{hence,xpvm}.

\section {Comparison of Numerical Advection Models}

In this section we use a case study to demonstrate the feasibility of 
integrating CMFortran and CMMD
modules on a CM5 in the proposed environment.

\begin{figure}[tb]
\centerline{\epsfxsize=395pt \epsfbox[50 330 550 720]{chapter3-multiparadigm/fig6.eps}}
%BoundingBox: 50 330 550 720
\caption{System configuration of the prototype system}
\label{config1}
\end{figure}

\begin{figure}[tb]
\centerline{\epsfxsize=360pt \epsfbox[36 188 576 612]{chapter3-multiparadigm/fig7.ps}}
%BoundingBox: 36 188 576 612
\caption{A graphical user interface of the prototype system}
\label{gui1}
\end{figure}

In fluid dynamics, the simple one-dimensional advection equation (1) 
with constant positive
velocity serves as a basic model for
the shape-conserving movement of an initial distribution of fluid volume
toward positive $x$. 
Since the analytic solution is known in this
simple case, numerical approximations to the advective process 
can be critically evaluated on fundamental properties such as 
accuracies, stabilities,
coordinate systems, and computer time and memory requirements.

\begin{equation}
\frac{\partial \mu}{\partial t} + u \frac{\partial \mu}{\partial x} = 0
\end{equation}

where $\mu$ is the fluid velocity and
$u$ is the mixing ratio $\mu$ $\equiv$ {\em C}/$\rho$ in which {\em C} is the 
constituent density and $\rho$ the density of the fluid. 

Our case study is to create a prototype interactive computing/visualization
environment in which a couple of 1D advection models are purposely implemented
in different programming paradigms on the CM5 and the results
graphically compared through a graphical user interface.

We choose from \cite{rood} five numerical advection algorithms with different
diffusion, dispersion and monotonicity properties. They are 
{\em Donor-cell}, {\em Partial donor-cell}, {\em van Leer}, 
{\em Hain} and {\em Lax-Wendroff},
named after the authors of the algorithms.
{\em Donor-cell} is then implemented in CMFortran and {\em van Leer}
in CMMD, while the 
rest of the models are implemented in Fortran77. All 
the {\em AVS/CM interface part} codes are
written in C.

To demonstrate the system's integration capability on a heterogeneous
architectures, we purposely run the AVS kernel on a IBM RS/6000, 
the {\em Hain} module on a DEC5000, and the {\em Lax-Wendroff} module on a SUN4, 
while the data-parallel {\em Donor-cell} module, the message-passing {\em van Leer} module
and the sequential module {\em Partial donor-cell} are running
on a 32-node CM5. The system configuration is shown in Figure \ref{config1}.

Figure \ref{gui1} is a screen-dumped picture showing the system parameter control panel (left),
model output windows (top) and the flow network (bottom). The flow-chart-like
diagram is the process configuration built by the AVS visual programming interface
Network Editor, in which {\em Advection Interface} is a system control module for
model input parameters steering and the {\em graph viewer} is an AVS system module 
for model graphical output. Both run on the same machine
as the AVS kernel. The other modules correspond to the 
respective advection models.
 
This prototype environment has the following features:

\begin{enumerate}
\item 
{\em User interactive control of the modeling}. 
Before or during a modeling run, the user can graphically choose 
advection models, set model parameters, display time steps and
the total number of grid points, and time steps.
 
\item
{\em Graphical output of the model results}.
The user can decide the displaying pattern, such as line, area, bar,
and other graphics parameters such as title, color, etc.

\item
{\em Flexible enough to include other advection models}.
It can be easily extended to add other old/new
models with varying velocity fields, multidimensions
and non-rectangular coordinated systems, or sequential, data-parallel, or 
message-passing implementations. All the models can be transparently configured to run
in a distributed computing environment.
\end{enumerate}
 
\section {Option Price Modeling on an IBM SP2 
and a DEC Alpha Cluster}
  
Stock option pricing models are used to calculate a price for an option
contract based on a set of market variables (e.g., exercise price, 
risk-free rate, time to maturity) and a set of model 
parameters.  Model price estimates are highly sensitive to parameter
values for volatility of stock price, variance of the volatility, and 
correlation between volatility and stock price.  These model parameters 
are not directly observable and must be estimated from market data.

We use a set of four option-pricing models in this study.  Simple 
models treat stock price volatility as a constant, and 
price only European (option exercised only at maturity of contract) options.  
More sophisticated
models incorporate stochastic volatility processes, and price American
contracts (option exercised at any time in life of contract).
These models are computationally intensive and have significant communication 
requirements.  The four pricing models are: BS --- the Black-Scholes constant
volatility, European model; AMC --- the American binomial, constant volatility
model; EUS --- the European binomial, stochastic volatility model;
and AMS --- the American binomial, stochastic volatility model. 

\begin{figure}
\centerline{\epsfxsize=430pt \epsfbox[60 75 550 705]{chapter3-multiparadigm/pvm/configure.eps}}
%BoundingBox: 0 0 612 725
\label{config2}
\caption{System configuration of the interactive distributed
financial modeling prototype}
\end{figure}

Figure \ref{config2} is the system configuration of a prototype 
interactive simulation-on-demand system for the 
option-price modeling application, using the AVS/PVM framework we
proposed in Section 3.3 and utilizing the network infrastructure
and distributed computing facility at NPAC.

The AVS kernel runs on a SUN10 workstation which acts as both an AVS server
to coordinate data-flow and top-level concurrent control among
remote modules and as a network gateway that links the NPAC in-house
host machines locally networked by an Ethernet to the regional end-user 
through the ATM-based WAN (in the actual demonstration, a NYNET link to 
Rome Laboratory, Griffiss Airbase, NY). The ATM-based link is built around
two Fore switches that operate at 155 Mbps (OC3c), while the wide area
network portion of the network operates at OC48(2400 Mbps) speed.

The two parallel pricing models (EUS model and AMS model) are 
implemented in PVM and run,
respectively, on an 8-node IBM SP2, networked by an Ethernet at the time of evaluation,
and an 8-node DEC Alpha cluster interconnected by an FDDI-based GIGAswitch.
They are coupled under the proposed AVS environment with the other two
sequential simple models (BS model and AMC model)
running on a SUN4 and a DEC5000 workstation, respectively. 
The nodal processor of SP2 is an IBM RISC/6000 processor running at 62.5 MHz
and is one of the most powerful processors available. The DEC Alpha
farm consists of eight Alpha
model 4000 workstations that are supported by a high-performance networking
backbone of dedicated, switched FDDI segments. The GIGAswitch
provides full FDDI bandwidth and low-latency switching to every workstation
in the farm. 

While displayed on the end-user's home machine, a user interface actually 
runs on a remote SUN4 that combines user runtime input (model parameters, network 
configuration) with historical market databases stored on disk, and 
broadcasts this data to remote compute nodes. Top-level system synchronization 
occurs with each broadcast.  

An IBM RS/6000 is used as a file server for non graphical
output of model data.  In this application, model prices calculated at
remote compute nodes and corresponding market data are written to databases 
for later analysis.

All models' output are graphically displayed on the end-user's home machine (a SUN10)
in AVS graph viewers. 
Figure \ref{gui2} gives the user interface
showing the simulation control panel (left), model output windows (top) and
the flow network (bottom).

Notice that with the use of PVM software, only one AVS license, i.e., the one on the Sun10,
is needed for this heterogeneous distributed system. 

\begin{figure}
\centerline{\epsfxsize=450pt \epsfbox[95 384 517 773]{chapter3-multiparadigm/pvm/screen.ps}}
%BoundingBox:  95 370 517 773
\caption{The Graphical User Interface}
\label{gui2}
\end{figure}

By using portable PVM software, we expect to build similar 
integrated parallel computing/interactive visualization 
applications on other MPP architectures where PVM is supported, while
AVS is not available currently on host processors such as the Intel 
Paragon and Cray T3D.

\section{Distributed Computational Electromagnetics Systems}

In this section we describe our work of developing a real-world
electromagnetic application on distributed computing systems. A
computational electromagnetics (CEM) simulation for radar
cross-section (RCS) modeling of full-scale airborne is ported to three
networked workstation cluster systems: an IBM RS/6000 cluster with Ethernet
connection, a DEC Alpha farm connected by a FDDI-based Gigswitch and an
ATM-connected SUN IPXs testbed, and two MPP machines: an Intel iPSC/860
and a CM-5. In the cluster implementation we use a ScaLAPACK LU solver from
Oak Ridge National Laboratory/University of Tennessee in the parallel 
implementation for solving a dense matrix
that forms the computationally intensive kernel part of this
application, and adopt BLACS as the message-passing interface
in overall code development to achieve high portability across the three
configurations. In the MPP implementation on the CM5, we use the
mixtured programming models in which the
matrix solver is from the vendor-supplied CMSSL, which is written 
in data-parallel CMFortran. The rest of the application is written
in explict message passing on the CMMD.
Their performance data are reported, as well as
(for comparison purposes) the timing data from all the platforms implemented.
This work shows that in order to achieve optimized performance
on MPPs coupled with clusters, a processes-based programming framework
is needed in which components and modules written in different
programming models can be seamlessly combined.

\subsection{Introduction}

Traditional electromagnetic engineering simulations are largely limited by
memory requirements, as well as by sequential processing time. 
Most electromagnetic applications on high-performance computing
systems so far experimented with or implemented have been on either massively 
parallel processors or traditional vector-based supercomputers \cite{tom1,tom2}. 
The advanced development in workstation cluster technology, represented by
clusterable IBM RS/6000 systems and Digital's Alpha farm with
both high-performance node processors and large memory capacity, provides new
opportunities and solutions to computational electromagnetic applications
that are both CPU and memory intensive. Together with the emergence of 
portable distributed linear algebra packages such as BLACS and ScaLAPACK and 
portable message-passing interface such as PVM and MPI, cluster-based computational
electromagnetics systems give an attractive solution to numerical
simulations in engineering electromagnetic analysis and design with
modest problem sizes. They achieve better cost/performance than
the MPP solutions, which is especially important to real-world
applications where workstation clusters are more accessible to
engineers. On the other hand, MPP is the only solution for solving
big CEM problems when the memory requirement cannot be fulfilled by
distributed systems.

A major motivation in this related work 
was to port the CEM application
to multiple MPP and distributed computing systems and evaluate/demonstrate
their capabilities for solving this problem. 
Portability, as well as performance scalability, was among the top objectives
of our effort \cite{project}.
With the distributed linear algebra 
package Scalapack \cite{scalapack1,scalapack2}
and BLACS \cite{blacs2} in our 
parallel implementation, we were able 
to achieve true code portability 
across Intel machines (the Intel Delta, an Intel Paragon, and an Intel iPSC/860),
an IBM SP-2 with either Ethernet or switch connection, 
and the three distributed configurations.
We focused our discussion on benchmarking and performance requirement aspects
of this CEM application on the experimented cluster and MPP systems.

\subsection{The RCS Application Problem}

The CEM application used in this work is a well-established CEM
package from Syracuse Research Corporation (SRC) called ParaMoM, which
stands for Parametric Patch Method of Moments (MoM). The most distinctive
feature of the ParaMoM is its use of basis functions that conform to
parametric surfaces with curvature and thus generate accurate and stable
simulation results.  

ParaMoM utilizes a general surface formulation. All surface integrals are 
performed on the curved surface rather than on a flat facet approximation,
as has been traditionally done. The basis (expansion) functions are
a parametric surface generalization of the well-known Rao-Witon-Glisson 
functions \cite{Rao}. The domains of these functions are 
curved three-sided surface patches that conform exactly to 
the underlying global surface representation.
With the parametric surface-patch approach as a foundation, 
a comprehensive set of capabilities was developed and placed in
ParaMoM 1.0\cite{src}, which has the following features:
\begin{itemize}
\item
Accepts IGES 114 and IGES 128 parametric surfaces, flat facets,
and IGES 110 wires.
\item 
Wire antenna modeling includes user-defined loads, feed points, and
connections to surfaces. Wires may radiate or scatter.
\item
E, H, and combined-field formulations.
\item
Up to three symmetry planes may be used to reduce memory and 
computation requirements for problems with geometric symmetry.
\item
Impedance-sheet formulation for treatment of resistive cards (R-cards)
\item
Graphics-based preprocessor to perform gridding and other model
preparation functions not included in a general-purpose CAD package.
\end{itemize}

\subsection{Parallel MoM Implementation Strategy}
The ParaMoM code consists of five major consequent processing phases that have
different computational requirements in terms of CPU time and memory
consumption, as shown in Table B.1.
\begin{table}[hbtp]
\caption {Processing Phases of Sequential  ParaMoM and 
Their Computational Requirements}
\begin{center}
\begin{tabular}{|c|c|c|c|} 
 \hline
Phase  &      Component              &  CPU time &  Memory \\ \hline\hline
 1     &       Setup                 &  $O(N)$     &  $O(N)$   \\ \hline
 2     &     Matrix fill             &  {\boldmath $O(N^{2})$}   & {\boldmath $O(N^{2})$} \\ \hline
 3     &    RHS vector fill          &  $O(N)$     &  $O(N)$   \\ \hline
 4     &  Matrix factor/solve        &  {\boldmath $O(N^{3})$}   &  {\boldmath $O(N^{2})$} \\ \hline
 5     & Far-field computation &  $O(N)$     &  $O(N)$   \\ \hline
\end{tabular}
\end{center}
\end{table}

Where $N$ is the number of unknowns, i.e., the number of basis functions,
it is proportional to the surface area of the target. 

The role of the setup phase is to read in the geometric target description
(i.e., the parametric surface representation along with the triangulation,
material, and wire connection information) and to precompute specific
geometry information required by the fill algorithm (including
positions and tangent vectors for the surface points used in matrix 
fill expansion and testing integrals). 
Because these geometry results will not be modified during the subsequent phases, 
there is no data coherency concern. At this phase we used
a block-cyclic
partition so that each node precomputes evenly apiece of the 
geometry arrays and
broadcasts its results to the rest of the nodes. Each node then maintains 
a private copy of these arrays of reasonable size ($O(N)$) in its
local memory.

The MoM matrix filling is well suited to distributed memory systems.
To better understand its parallel implementation, Program \ref{code} 
gives the sequential 
filling algorithm in pseudo-Fortran code.
The sequential algorithm has the following features, which are
well utilized and highlighted in our SPMD parallel algorithm:
\begin{enumerate}
\item
Calculation of any Z-matrix element is completely independent
of other elements in the matrix and is only dependent on
its memory allocation. This property makes it
possible to perform the Z-matrix calculations locally ``in-place''
of the processors on which the Z-matrix elements are stored, and the 
performance of the parallel algorithm depends solely on the
locality of the Z-matrix partition, not on the networking topology
of physical processors.

\noindent \mbox{}\hrulefill\mbox{} 
{\tt
\begin{tabbing}
c loop o\= ver source patch \\
        \>Do isource = 1, nfaces \\
c loop over field  patch \\
        \>\hspace{0.25in}\= Do ifield = 1, nfaces\\
c calcul\>ate integration points for patch pair (ifield,isource)\\
        \>   \>\hspace{0.25in}\= call patch-interaction(ifield,isource)\\
c loop o\>ver basis functions of the source patch\\
        \>   \>   \>Do isrcnode = 1, 3\\
c look u\>p the index in pre-defined table\\
        \>   \>   \>\hspace{0.25in}\=j = lookup(isource,isrcnode)\\
c loop o\>ver integration points of the field patch\\
        \>   \>   \>   \>Do ipnt = 1, 4\\
c calcul\>ate potential integrals for the patch pair -- main routine\\
        \>   \>   \>   \>\hspace{0.25in}\=call pots(isource,ifield,isrcnode)\\
c loop o\>ver basis functions of the field patch\\
        \>   \>   \>   \>   \>Do ifldnode = 1, 3\\
c look u\>p the index in pre-defined table\\
        \>   \>   \>   \>   \>\hspace{0.25in}\=i = lookup(ifield,ifldnode)\\
c calcul\>ate intermediate result              \\  
        \>   \>   \>   \>   \>   \>tmp = exp(ifield,ifldnode,ipnt)\\
c accumu\>late the Z-matrix element where cz is the N by N Z-matrix\\
        \>   \>   \>   \>   \>   \>cz(i,j) = cz(i,j) + tmp\\
        \>   \>   \>   \>   \>Enddo ifldnode\\
        \>   \>   \>   \>Enddo ipnt\\
        \>   \>   \>Enddo isrcnode \\
        \>   \>Enddo ifield \\
        \> Enddo isource 
\end{tabbing}}
\noindent \mbox{}\hrulefill\mbox{}\\
\label{code}
Program \ref{code} The Sequential Matrix Filling Algorithm\\

\begin{figure}
\centerline{\epsfxsize=450pt \epsfbox[60 450 550 710]{appendix/rcs/paramom.eps}}
%%BoundingBox: 84 590 461 767
\label{paramom}
\caption{Program structure of the ParaMom implementation}
\end{figure}

\begin{figure}
\centerline{\epsfxsize=450pt \epsfbox[50 570 550 720]{appendix/rcs/block.eps}}
%%BoundingBox: 84 590 461 767
\label{block}
\caption{Patch-patch interaction and the Z-matrix}
\end{figure}

\item 

Determined by the characteristics of the numbering of 
basis functions among
patches, calculation for a patch pair $(isource,ifield)$ 
is based on the 9 associated basis functions $(isrcnode=1:3,ifldnode=1:3)$ 
and is eventually filled into the Z-matrix
as shown in Figure B.2, where $j, l$ are column-index related to
the source path ($isource$) and $i, k$ are row-index related to
the field patch ($ifield$).
From Figure B.2 it is clear that a
column-major blocked partition of the Z-matrix achieves the best locality 
for a patch pair. It guarantees that at least 
6 Z-matrix elements are mapped to a processor and that the maximum
number of processors used to calculate the total 9 elements
is 2. In an average case when a general SSB partition is used,
the maximum number of processors for a patch pair is 4.
This is to say that, compared to the ideal 
partition in which all the 9 elements are mapped to a single
processor, the only overhead for the column-major block partition
is one redundant computation of the routine ``patch-interaction''
which takes less than 20\% of the total computation for a patch pair.

\item The locality issue examined in (2) is for a single 
patch pair. We must also investigate the load-balance issue for all patch pairs
to see if the overall Z-matrix computation is 
evenly distributed to all processors.
Fortunately, the numerical model statistically assures that the Z-matrix 
computation of all patch pairs is evenly distributed. Actually,
for all Z-matrix elements except those related to boundary patches,
an element's final value is accumulated from four patch pairs and
from (2) we know a patch pair contributes to nine different Z-matrix elements
in most cases. Our measurements in testing cases also proved that this is
statistically true.

\item 
It becomes clear that in (2) a larger Z-matrix size gives better
data locality, and in (3) a larger machine size makes a relatively better 
balanced load of Z-matrix among processors. Therefore,
the parallel filling algorithm will guarantee a scalable performance 
in terms of machine size and problem size.

\end{enumerate}

Our complete parallel implementation includes:
\begin {enumerate}

\item 
In the setup phase, a parallel file input for target geometry 
data and a parallel algorithm for 
precomputation in which only global broadcast operations are 
used.
 
\item 
An embarrassingly parallel algorithm for the Z-matrix filling that
   requires some redundant calculation of patch-patch interactions but
   no inter-node communication.

\item 
A similar embarrassingly parallel algorithm for the RHS excitation vectors 
   filling that has no inter-node communication.

\item
A COMPLEX parallel LU solver. Our approach is to utilize parallel
matrix equation solution software developed by others. Such 
software includes vendor-specific commercial software
such as the CMSSL matrix equation solver on the CM5 that
is written in data-parallel and is highly optimized for vector-processing, 
as well as the portable, publicly available software ScaLAPACK from
Oak Ridge National Lab.

\item 
A parallel algorithm for far-field calculation that requires 
only global summation operations.

\end{enumerate}

Our algorithm design is focused on Phase 2 of filling the impedance
matrix (Z-matrix).
The major task in our implementation is to partition the computation 
of Z-matrix elements among
processors by distributing its memory among processors.  

The Basic Linear Algebra Communication Subprograms (BLACS)\cite{blacs1} comprise a 
library that provides ease of use and portability for message-passing,
matrix-based operations in
parallel linear algebra programs on distributed memory machines. It supports
both point-point operations between processors and collective communications
on a virtual two-dimensional processor grid.
ScaLAPACK, the two-dimensional distributed memory version of common
LAPACK\cite{LAPACK} routines, relies on calls to the BLAS for local
computation and on calls to the PB-BLAS\cite{pbblas} for global computation.
For portability reasons, communications take place inside
the PB-BLAS through calls to the BLACS.

We were fortunate in gaining early access to the ScaLAPACK library
and adopted an early release of the REAL variable arithmetic version 
in writing our COMPLEX routines.
It is important to parallelize the {\it full}
applications code so as to optimize and balance the overall performance of this
application \cite{hpcn}. The BLACS interface also provided a suitable level of
abstraction for us to implement the ``Setup'' and ``Far-field''
components. The full schematic of our 
parallelized ParaMoM code in cluster implementation
is shown in Figure B.2.

ScaLAPACK is currently not an efficient matrix solution option on the
CM5, because a BLAS library that is optimized for the CM5 vector unit does not
exist. Therefore, the platform-specific CMSSL matrix equation
solver is being used. 

The CMSSL library uses a data-parallel programming model and cannot be directly
interfaced with the message-passing, matrix-filling algorithm.
Therefore, in our implementation matrix filling and factoring
are two distinct program modules. The high-bandwidth parallel file
system on the CM5 (DataVault) is used as a buffer between the two stages.
The message-passing MoM matrix-filling program fills the matrix and writes it 
to the file system in the format required for the factoring stage in CMSSL.
The matrix is subsequently read in by the data-parallel matrix solver stage.

There is very little performance penalty for using separate program
modules for the filling and factoring, because the DataVault has very high
data rates. For example, Table B.2 gives the timings for a 
complex matrix size = 2000*2000.

\begin{table}[tb]
\caption {Timings of parallel I/O and LU solver for $N = 2000$ (in seconds)}
\begin{center}
\begin{tabular}{|l|c|} \hline
Operation  & Time\\ \hline\hline
Write matrix & 2.3 \\ \hline
Read matrix & 1.5\\ \hline
LU factorization & 24.5 \\ \hline
Solve & 3.2 \\ \hline
\end{tabular}
\end{center}
\end{table}

For this moderate-sized problem the I/O time is small compared to the factor time.
Since the I/O time scales as N*N and the LU time scales as N**3, it is evident that for large problems
the I/O time will not be a limiting factor.

The initial implementation of the matrix-filling algorithm on the CM5 is being
written in F77 with calls to the CMMD message-passing library. The resulting code
is adequate to investigate the efficiency of the parallel algorithm as measured by the
speed-up
(ratio of elapsed time for one node versus elapsed time for M nodes). However, the per-node
performance is poor because the F77 code cannot utilize the CM5 vector units efficiently,
and the matrix-filling part is best optimized in explicit message-passing programs.
As indicated in our benchmarking timing data, to optimize the overall performance
of this application, a coupled CM5 and cluster solution is required so that
the matrix solver achieves the best performance on a CM5 and the filling is optimized
on a clustered IBM SP2. To integrate the two components written in different programming
models interactively, a unified integration framework is required so that results of the
matrix filling part on the SP2 can be directly fed into the LU solver running on the CM5,
bypassing the buffering on disk files as currently used in the stand-alone CM5 approach.

\subsection{Performance Comparison and Discussion}

We report in this section timings of the Z-matrix filling and LU 
factorization parts which take more than 95\% of the total execution time.
The filling part is almost 50\% dominated by MFLOPS and 50\% by memory
access, while the LU part is more than 90\% dominated by MFLOP. 
Elapsed time for filling and both elapse time and MFLOPS
for LU are listed.

\begin{figure}[thb]
\centerline{\epsfxsize=430pt \epsfbox[29 285 583 510]{appendix/rcs/lu.ps}}
%\centerline{\epsfxsize=430pt \epsfbox[50 470 550 720]{appendix/rcs/lu.ps}}
\label{time}
\caption{Fill Time v.s Processors (Matrix Size = 988)}
\end{figure}

The DEC Alpha farm consists of eight Alpha model 4000 workstations running at
133 MHz and with 64 MB memory per processor,
interconnected by a dedicated Gigaswitch that provides full FDDI
bandwidth and low-latency switching to every processor in the farm.
The Ethernet-based IBM cluster has eight IBM RS/6000 model 370 processors
running at 62.5 MHz each and with 64 MB memory per processor.
The ATM link between two SUN IPXs operates at a peak bandwidth of 155 Mbps(OC3c).
Each SUN IPX has 32 MB memory and runs at 40 MHz.

Different, modest problem sizes are used to show the
performance comparison and scalability.
For comparison
purposes, timings of this application on an Intel iPSC/860 (40 MHz, 16 MB/node) 
and a CM5 (25 MHz, 16 MB/node) are also listed. The Intel 
implementation, with BLACS linked 
to its NX version, has
exactly the same source code as that on the clustered systems,
while the CM5 implementation
used a completely different approach in which the fill part was
written in a message-passing CMMD program and the LU part was a
data-parallel CMFortran program using a vendor-supplied CMSSL LU
solver to utilize CM5's high-performance vector units.   
The CM-5 implementation used out-of-core DataVault (or SDA) disk memory
to hold and transfer the Z-matrix between the two different components.
For some problem sizes, more processors have to be used on the 
iPSC/860 and CM-5 due to the memory requirement. Notice that a block
size $N_b$ = $N/N_p$, $P=1$ and $Q=N_p$ are used for all the tests reported 
here, where $N_p$ is the total number of processors used, $N$ the matrix
size and  $P$ $\times$ $Q$ the virtual processor grid. As discussed in
Section 3, this partition is optimized for the matrix fill algorithm
on all the platform and good for the LU part on the three distributed
systems, but it is not optimized for the LU solver on the iPSC/860. 

Table B.2 gives an overall breakup of timings that indicates the two
dominated components and relatively balanced results.
In Tables B.3 and B.4, the clustered systems show compatible or even
better performances
than the MPP platforms with a comparable number of processors,
although the LU part is poorer, limited by the networking bandwidth
and high latency of message-passing protocols on the cluster systems.
The LU performance on the ATM-based configuration (Tables B.5 and B.6) is extremely
poor and non scalable with the problem size. There may be two reasons:
(1) TCP/IP instead of ATM's own protocol is used, whose high latency and
that of PVM outweigh the overall bandwidth; (2) no optimized BLAS
library is utilized.

ScaLAPACK makes use of
any system-supplied BLAS routines that will help performance
optimization at the node level. This is important, since many of the
recent parallel architectures employ vector operations at the node
level, and processors such as the i860 employed by the Intel iPSC/860
are well endowed with vector registers and can carry out BLAS operations
very efficiently. Furthermore, even on the super-scalar nodes used by
machines like the IBM RS/6000 cluster, the BLAS operations provide a
valuable level of abstraction for node level optimization.

The scalability of fill and LU parts is shown in Figure \ref{time}.

\begin{table}
\begin{center}
\caption {Timing Comparisons for $N = 4889$ (in seconds)}
\begin{tabular}{|l|r|r|r|c|r|c|r|} \hline
Platform  & $N_{p}$ &  Fill  &  LU   & LU (mflop)  & Setup & RHS+Field & Total\\ \hline\hline
Alpha(FDDI) & 8 & 1420.0  & 1119.9 & 295.5  & 12.3 & 18.0 & 2570.2 \\ \hline
IBM RS(Ethernet) & 8 & 1500.7 & 1804.7 & 183.4 & 51.4 & 28.2 & 3385.0\\ \hline
Intel iPSC/860 & 64 & 525.6 & 280.9 & 1178.8 & 45.4 & 53.0 & 904.9 \\ \hline
CM-5 & 32 & 3295.2 & 170.7 & 1938.8 & 21.1 & 4.3 & 3491.3 \\ \hline
%SUN(ATM) & 2 & \multicolumn{3}{c|}{not enough memory}  \\ \hline
\end{tabular}
\bigskip 

\caption {Timing Comparisons for $N = 3060$ (in seconds)}
%\begin{center}
\begin{tabular}{|l|r|r|r|c|} \hline
Platform  & $N_{p}$ &  Fill  &  LU   & LU (mflop)\\ \hline\hline
Alpha(FDDI) & 8 &  450.4  & 713.5 & 91.6    \\ \hline
IBM RS(Ethernet) & 8 & 495.1 & 1210.3 & 54.0 \\ \hline
Intel iPSC/860 & 16 & 830.9 & 175.6 & 372.1 \\ \hline
CM-5 & 32 & 1207.0 & 79.5 & 821.9   \\ \hline
%SUN(ATM) & 2 & \multicolumn{3}{c|}{not enough memory}  \\ \hline
\end{tabular}
\bigskip 

\caption {Timing Comparisons for $N = 1504$ (in seconds)}
\begin{tabular}{|l|r|r|r|c|} \hline
Platform  & $N_{p}$ &  Fill  &  LU   & LU (mflop) \\ \hline\hline
Alpha(FDDI) & 8 &  114.9  & 76.2 & 110.5 \\ \hline
IBM RS(Ethernet) & 8 & 127.1 & 131.8 & 63.9   \\ \hline
Intel iPSC/860 & 8 & 478.6 & 61.5 & 136.9 \\ \hline
CM-5 & 32 & 358.9 & 11.0 & 824.3 \\ \hline \hline
SUN(ATM) & 2 & 2087.4 & 1941.4 & 4.3 \\ \hline
Alpha(FDDI) & 2 &  411.0  & 278.6 & 30.2    \\ \hline
IBM RS(Ethernet) & 2 & 475.6 & 262.0 & 32.1 \\ \hline
Intel iPSC/860 & 2 & 1950.7 & 93.9 & 89.7 \\ \hline
\end{tabular}
\end{center}
\end{table}
\bigskip 

\begin{table}
\begin{center}
\caption {Timings of Different $N$ on 2 SUN IPXs Interconnected by an ATM Link}
\begin{tabular}{|c|r|r|c|} \hline
$N$ & Fill & LU & LU (mflop)\\  \hline\hline  
1504 & 2087.4 & 1941.4 & 4.3 \\ \hline
988 & 981.1 & 560.6 & 4.3 \\ \hline
672 & 449.3 & 176.4 & 4.5  \\ \hline
468 & 229.2 & 60.7 & 4.3 \\ \hline
368 & 147.8 & 29.9 & 4.1 \\  \hline
187 &  45.4 & 4.0 & 4.1 \\  \hline
\end{tabular}
\end{center}
\end{table}

\subsection{A Scheduling and
Visualization Framework for the RCS Simulation on a Meta-Computer}

The high portability achieved in the parallel code development
across multiple parallel platforms, along with the
significant importance of an interactive
visualization requirement, provides an excellent opportunity
for us to develop a real-time scheduling and visualization interface
for this application. The objectives of this prototype are three-folded,
as follows:

\begin {enumerate}
\item
An interactive visualization interface to allow real-time
simulation parameter steering and graphical display of
simulation results of far-field scattering, along with the
three-dimensional geometry of the RCS target,

\item
A high-level scheduling mechanism to facilitate runtime allocation
of simulation job among multiple parallel systems available 
for this application.

\item
A single integrated framework to couple all the components required
by this prototype system that are executed under a 
locally-networked environment, including a GUI interface for system control
and visualization, a runtime scheduling
interface for parallel job allocation, and
the parallel simulation job to be executed on one or more of the scheduled
parallel or distributed computers.

\end {enumerate}

\begin{figure}
\centerline{\epsfxsize=430pt \epsfbox[50 110 550 720]{appendix/rcs/schedule.eps}}
%%BoundingBox: 11 3 601 790
\caption{A Integrated Scheduling, 
Visualization and Simulation Framework on a Meta-Computer}
\label{schedule}
\end{figure}

We envision using the proposed dataflow integration framework
as an integrator to facilitate all the components.

Figure \ref{schedule} gives a skeleton of a flow diagram for the framework, in which
the scheduling factor is platform-dependent and is calculated
based on the following common performance parameters:

\begin {enumerate}
\item
Simulation-related parameters: (for current determined RCS target)
in-core memory, number of flop counts (or CPU, memory measurements)
\item
Platform-related parameters: (for current meta machine in place)
clock rate, theoretical peak performance, number of processors,
in-core memory size/node, communication bandwidth,
(parallel) I/O bandwidth and memory bandwidth
(in general, CPU, memory, I/O, and networking)
\item
Current runtime system load parameters: (determined at runtime)
system status (up/down), number of users, free memory available,
load average (uptime, the average number of jobs in the run
queue over the last 1, 5 and  15  minutes)

\item
User runtime defined parameter: by default all values are 1.
User can set it via a GUI to define the priority machine.
Simple or sophisticated scheduling functions
can be enforced. For example, here is a simple scheduling scheme:

\[ f_{i}(a,b,c,d) = \left\{ \begin{array}{ll}
			10\%*a + 20\%*b + 30\%*c + 40\%*d & \mbox{if the machine is up}\\
			0 & \mbox{otherwise}
		      \end{array}
	      \right. \]

where {\it a,b,c} is a parameter value calculated from the above 1,2,3, respectively.
{\it d} is a user-defined number that can be graphically steered from each
module's control parameter panel associated with the AVS remote module,
and $0 <= a,b,c,d <= 1$.

For the filter module, the decision is based on
$Max(f_{1},f_{2},f_{3},f_{4},f_{5})$. 
\end {enumerate}

Figure \ref{vfy} is the GUI of this prototype showing an EM scattering field from 
the full-scale airborne system, which is generated by the ParaMOM MPP code
on an Intel Paragon, while Figure \ref{sphere} also shows the runtime parameters
control panel on the GUI.

\begin{figure}[tb]
\centerline{\epsfxsize=360pt \epsfbox[59 448 560 780]{appendix/rcs/vfy.ps}}
%BoundingBox: 36 188 576 612
\caption{A graphical output of the prototype system (a VFY)}
\label{vfy}
\end{figure}

\begin{figure}
\centerline{\epsfxsize=250pt \epsfbox[106 456 513 787]{appendix/rcs/rcs.ps}}
%BoundingBox: 36 188 576 612
\caption{The graphical user interface of the prototype system (a sphere}
\label{sphere}
\end{figure}
