\documentstyle[12pt,algorithm]{cmu-art}
\title{Mapping Extended Kalman Filters onto Linear Arrays}
\author{R. S. Baheti, D. R. O'Hallaron, and H. R. Itzkowitz\thanks{Manuscript received August 10, 1989; revised May 16, 1990. Paper recommended by Associate Editor at Large A. Laub. This work was done while the authors were at General Electric, Schenectady, NY. Currently, R. Baheti is with the National Science Foundation, Washington, DC; D. O'Hallaron is with Carnegie Mellon University, Pittsburgh, PA; H. Itzkowitz is with AT\&T Bell Laboratories, Naperville, IL.}}
\date{}
\begin{document}
\bibliographystyle{ieeetr}
\maketitle
\begin{abstract} 
The paper describes techniques for mapping extended Kalman filters
onto linear arrays of programmable cells designed for real--time
applications. First, we introduce a general method for mapping a
standard (non--square root) Kalman filter, where the columns of the
covariance matrix are updated in parallel. Next, we introduce a
general method for mapping a factorized (square root) filter, where
fast Givens rotations are used to triangularize the pre--matrix, and
where rotations of the rows of the pre--matrix are performed in
parallel. We use these mappings to implement an extended Kalman filter
commonly used in target tracking applications on the Warp computer.
The Warp is a commercially--available linear array of 10 or more
programmable cells connected to an MC68020--based workstation. The
Warp implementation of the standard Kalman filter running on 8 Warp
cells achieves a measured speedup of 7 over the same filter running on
a single cell. The Warp implementation of the factorized filter
running on 10 Warp cells achieves a measured of speedup of 2.  The
Warp implementations of both the standard and factorized filters
achieve roughly two orders of magnitude speedup over the same filters
implemented on an MC68020--based workstation. With the commercial
availability of integrated Warp chips, these results suggest radically
new solutions to on--board target tracking and other Kalman filtering
applications.
\end{abstract}
\section{Introduction}
The extended Kalman filter is widely used in real--time estimation
applications. The filter mechanism is expensive, in general requiring
$O(n^3)$ floating--point computations, where $n$ is the number of state
variables.  This time--complexity has motivated the design of a number
of proposed parallel algorithms
\cite{andrews81,chen86,jover86,shaffer87b,sung87}, which are
designed for special--purpose VLSI architectures such as systolic
arrays and wavefront arrays.
For example, Jover and Kailath \cite{jover86} developed a parallel
algorithm for the Kalman filter measurement update based on stable
orthogonal transformations. Their paper provides many useful insights
into the measurement update, including the generalization of Bierman's
U-D factorization to handle vector updates. The filter realization
requires a systolic--type architecture with $2n + 4$ elementary
processors and a bank of delay lines. Based on these results, Shaffer
\cite{shaffer87b} designed a target tracking filter than can be
realized using a bit--serial silicon compiler. The investigation
uncovered some limitations of the systolic--type architectures,
including the use of a large number of custom VLSI chips, low
utilization of processing power, and inflexibility of the
architecture.

This paper reports on the results of an investigation into the
feasibility of implementing the extended Kalman filter on Warp, a
programmable commercially--available linear array.  A linear array is
the simplest configuration for a distributed--memory parallel machine,
and as such, can be embedded in many existing parallel configurations,
including hypercubes \cite{seitz85}, and two--dimensional meshes
\cite{borkar88,borkar90}. Thus, results obtained for linear arrays are
applicable to a wide variety of parallel machines.

In our particular investigation, we devised general methods for
mapping the standard (non--square root) Kalman filter onto a linear
array of $n+2$ cells, and for mapping a factorized (square root)
filter onto a linear array of $n+m+1$ cells, where $n$ is the
dimension of the state vector, and $m$ is the dimension of the
measurement vector.  Using these general mapping techniques, we
implemented a target tracking filter on the Warp parallel computer.
The Warp is a linear array of 10 or more identical and programmable
cells connected to an MC68020--based workstation.  Each cell has local
storage for programs and data, a pipelined floating point multiplier
and adder, and achieves a maximum performance of 10 MFLOPS. Each cell
can send and receive up to 20 million 32--bit data words to and from
its neighbors each second.  The details of the Warp architecture are
described in \cite{annaratone87e}.

The results of the investigation were mixed. On the one hand, the Warp
implementation of the standard Kalman filter (using 8 cells) achieved
a remarkable measured speedup of 7 over the single--cell filter, with
a corresponding efficiency of 87\%. On the other hand, the Warp
implementation of the factorized filter (using 10 cells) achieved a
somewhat disappointing measured speedup of 2 over the single--cell
filter, with a corresponding efficiency of 20\%. Thus, while we have
clearly demonstrated the feasibility of the standard Kalman filter, we
have also uncovered the need for further work in the design of
parallel factorized filters on general purpose parallel machines.

The paper is organized as follows. Sections \ref{sec:ekf}, 
\ref{sec:givens}, and \ref{sec:trackingfilter} give background
information on the extended Kalman filter, fast Givens rotations, and
a target tracking filter.  In Section \ref{sec:mapping} we describe
general techniques for mapping both the standard Kalman filter and a
factorized filter onto a simple linear array. In Section
\ref{sec:implementation} we use these techniques to
implement a 9--state tracking filter on Warp, and report on measured
performance. Finally, Section
\ref{sec:conclusions} is devoted to concluding remarks.
 
\section{Extended Kalman Filter }
\label{sec:ekf}
The extended Kalman filter requires on--line updates of the covariance
matrix that can be implemented in several ways.  In the paper, we
consider two methods of updating the covariance equations. The first
algorithm, which we call the {\em standard Kalman filter}, updates the
covariance matrix using scalar measurements \cite{bierman77}.  The
second algorithm, which we call the {\em factorized Kalman filter},
updates a factorized covariance matrix using fast Givens rotations
\cite{jover86}.  The system model and filter equations are presented in the
following.

A discrete--time state space model for an extended Kalman filter is given by 
\cite{anderson79}
\begin{eqnarray}
{\bf x}_{t+1} & = & {\bf f}_t({\bf x}_t)+ {\bf g}_t({\bf x}_t){\bf
w}_t \nonumber \\ {\bf z}_t & = & {\bf h}_t({\bf x}_t)+{\bf v}_t
\label{eqnnon}
\end{eqnarray}
where subscript $t$ is the discrete-time, ${\bf x}_t \in {\cal R}^n$ is a 
state vector, and
${\bf y}_t \in {\cal R}^m$ is a measurement vector. 
${\bf f}_t(\cdot),{\bf h}_t(\cdot) $ are in general nonlinear and 
$ {\bf g}_t(\cdot) $ is a nonconstant.  
The process noise $\{{\bf w}_t\}$
and the measurement noise $\{{\bf v}_t\}$ are zero mean, white gaussian random 
sequences.  The initial state ${\bf x}_0$ is a gaussian random vector with
mean ${\bf \bar{x}}_0$ and covariance $P_0$.
It is assumed that $\{{\bf w}_t\}$
,$\{{\bf v}_t\},\{{\bf x}_t\},{\bf x}_0$
are mutually independent, 
$ E({\bf w}_t {\bf w}_ {\tau} ^ T )  =  Q_t \delta _{t \tau}$ , and
$E({\bf v}_t {\bf v}_ {\tau} ^ T )  =  R_t \delta _{t \tau}$.
Let $\Phi _t$ and $H_t$ denote the Jacobian matrices evaluated at the
conditional means 
${\bf \hat{x}}_{t|t}$ and ${\bf \hat{x}}_{t|t-1}$, respectively, defined as
\begin{eqnarray}
\Phi _t  = \left. \frac { \partial {\bf f}_t({\bf x})}  { \partial {\bf x}} 
\right|  {{\bf{x}}  = {{\bf \hat{x}}_{t|t}}} 
,~~~~H = \left. \frac { \partial {\bf h(x)}}  { \partial {\bf x}} 
\right|  {{\bf{x}}  = {{\bf \hat{x}}_{t|t-1}}}
,~~~~G_t={\bf g}_t({\bf {\hat x}}_{t|t}) \label{eqnjac}
\end{eqnarray}
where ${\bf \hat{x}}_{t|t}$ and ${\bf \hat{x}}_{t|t-1}$ denote the
conditional mean estimates of ${\bf {x}}_t$ based on the measurements
${\bf y}_{t}$ and ${\bf y}_{t-1}$ respectively. If the nonlinear
functions in Equation \ref{eqnnon} are sufficiently smooth, then using
Taylor series expansion about $ {\bf \hat{x}}_{t|t}$ and ${\bf
\hat{x}}_{t|t-1}$, the extended-Kalman filter equations are given as

{\it Time Update:}
\begin{eqnarray}
{\bf {\hat{x}}}_{t+1|t} & = & {\bf f}_t( {\bf \hat{x}}_{t|t}) \label{eqntux} \\
P_{t+1|t} & = & \Phi _t P_{t|t} {\Phi _t} ^ T + G_tQ_tG_t^T \label{eqntup}
\end{eqnarray}
{\it Measurement Update:}
\begin{eqnarray}
{{\bf \hat{x}}}_{t|t}& = & {\bf \hat{x}}_{t|t-1}+K_t[{\bf y}_t-{\bf h}_t({\bf \hat{x}}_{t|t-1})] \label{eqnmux} \\
R_{e,t} &= & H_tP_{t}H_t ^ T + R_t \label{eqninnovcov} \\ 	
K_t & = &  P_t H_t^T  { R_{e,t}} ^ {-1} \label{eqngain} \\
P_{t|t} & = & P_{t|t-1}-K_tH_tP_{t|t-1} \label{eqnmup}
\end{eqnarray}
where $ R_{e,t} $ denotes  covariance of the innovations,  $K_t$
denotes the Kalman gain matrix, $P_{t|t}$ and $P_{t|t-1}$ denote the
appropriate covariance matrices.
The initial conditions for the recursion are given by ${\bf {\hat x}}_{0|-1}={\bf x}_0$ and $P_{0|-1}=P_0$.
 

\subsection{Scalar Measurement Update}

When the measurement noise covariance matrix $R_t$ is diagonal, it can
be computationally efficient to process one measurement at a time
\cite{bierman77}. In fact, it is possible to transform the $R_t$ matrix to
diagonal form using a Cholesky decomposition. Let
$R_t=diag[r_1,\ldots,r_m]$. Let $y_i$ and $h_i(\cdot)$ denote the $ith$
element of ${\bf y}_t$ and ${\bf h}_t(\cdot)$, respectively. The
subscript $t$ is omitted to simplify the notations. Let $H_i$ denote
the $ith$ row of the measurement partial derivative matrix $H_t$.
Then state estimates in Equation \ref{eqnmux} for $i=1,\ldots,m$ can be
written as
\begin{eqnarray}
{\bf \hat{x}}^{(0)} & = & {\bf \hat{x}}_{t|t-1} \nonumber \\
{\bf \hat{x}}^{(i)} & = & {\bf \hat{x}}^{(i-1)} + K^{(i)}[ y_{i}(t) - h_{i}
({\bf \hat{x}}^{(i-1)})] \label{eqnmuxk}\\
{\bf \hat{x}}_{t|t} & = & {\bf \hat{x}}^{(m)} \nonumber
\end{eqnarray}
A scalar measurement update of the covariance matrix using
Equations \ref{eqninnovcov}, \ref{eqngain}, and \ref{eqnmup}
can be written as 
\begin{eqnarray}
P^{(0)} & = & P_{t|t-1} \nonumber \\ 
R_e^{(i)}& =& H_iP^{(i-1)}H_i^T+r_i \label{eqninnovcovk} \\
K^{(i)} &=& P^{(i-1)}H_i^T (R_e^{(i)})^{-1} \label{eqngaink} \\
P^{(i)} & = & P^{(i-1)} - K^{(i)}H_{i}P^{(i-1)}, \label{eqnmupk} \\
P_{t|t} & = & P^{(m)} \nonumber
\end{eqnarray}

\subsection{Factorized Covariance Update}

In the following we consider a factorized filter formulation based on
stable orthogonal transformations.  A factorization of the form
$LDL^T$ rather than the Cholesky decomposition $LL^T$ is used to avoid
computations of the arithmetic square roots. The algorithm combines
the time and measurement update equations to estimate the state and
the covariance matrix. Jover and Kailath \cite{jover86} and Kailath
\cite{kailath87} have shown that Bierman's U-D factorized update is a
special case of the factorized filter described in this section. Even
though we consider the combined update to develop the mapping
strategy, in practice it may be preferable to use Bierman's U-D scalar
measurement update because of the possibility of missing
observations, non--uniform sampling rates, and data outliers.
 
The time update and the measurement update of the estimated state and the
covariance matrix can be combined as 
\begin{eqnarray}
{\bf \hat{x}}_{t+1|t} & = & \Phi _t +{\bf \hat{x}}_{t|t-1}+
\Phi _tK_t[{\bf y}_t-{\bf h}_t({\bf \hat{x}}_{t|t-1})] \label{time-x} \\
P_{t+1|t}&  =& \Phi _t P_{t|t-1} \Phi _t ^ T + G_tQ_tG_t^{T} - \Phi _t K_tR_{e,t}
K_t^T \Phi _t^T		\label{time-P}
\end{eqnarray}
Consider the matrix
  \[
	S = \left [ \begin{array}{cc}
			R_{e,t}	&  H_t P_{t|t-1} \Phi_t^T \\
			\Phi_t P_{t|t-1}H_t^{T}	& \Phi_t P_{t|t-1} \Phi_t^T+G_tQ_tG_t^T
		    \end{array}
	    \right ]
  \]
Define the following factorizations:
  \begin{eqnarray*}
	P_{t|t-1}	 &=& L \overline{D} L^{T}			\\
	P_{t|t}	 &=& L_{+}D_{+}L_{+}^{T}	\\
	R_{e,t} &=& L_{e}D_{e}L_{e}^{T}.
  \end{eqnarray*}
 The matrix $S$ can be factored as
  \begin{eqnarray*}
    S &=& ADA^{T} \\
      &=&
	\left [ \begin{array}{ccc}
		I & H_t L & 0		\\
		0 & \Phi_t L & G_t
		\end{array}
	\right ]
	\left [ \begin{array}{ccc}
		R_t & 0 & 0	\\
		0 & \overline{D} & 0	\\
		0 & 0 & Q_t
		\end{array}
	\right ]
	\left [ \begin{array}{cc}
		    I		&       0		\\
		L^{T}H_t^{T}	& L^{T} \Phi_t^{T}	\\
		    0		&	G_t^{T}		
		\end{array}
	\right ],
  \end{eqnarray*}
where $A$ is an $(n+m) \times (2n+m)$ matrix
and $D$ is a $(2n+m) \times (2n+m)$ diagonal matrix.
\footnote{The matrices $R_t$ and $Q_t$ are assumed to be diagonal.
If these matrices are not diagonal, the factorization can be changed
so that $D$ remains diagonal.}
Similarly, it is easily verified that $S$ also can be factored as

  \begin{eqnarray*}
    S &=& A'D'(A')^{T} \\
      &=&
	\left [ \begin{array}{ccc}
		L_{e}		&   0	& 0	\\
	     \Phi_t K_tL_{e}		& L_{+} & 0
		\end{array}
	\right ] 
	\left [ \begin{array}{ccc}
		D_{e}	& 0 	& 0	\\
		0 	& D_{+} & 0	\\
		0 	& 0 	& D_{a}
		\end{array}
	\right ] 
	\left [ \begin{array}{cc}
		L_{e}^{T}	& L_{e}^{T}K_t^{T}\Phi_t^T \\
		     0		& 	L_{+}^{T}	\\
		     0		&	    0	
		\end{array}
	\right ]
  \end{eqnarray*}
where $A'$ is an $(n + m) \times (2n + m)$ unit lower triangular matrix
and $D'$ is a $(2n + m) \times (2n + m)$ diagonal matrix.
Because of the nature of this factorization, the matrix $D_{a}$ is
arbitrary.

From these two block factorizations of $S$, it is known that there
exists an orthogonal transformation matrix $\Theta$ such that
\begin{eqnarray} \lefteqn{ 	\left [ \begin{array}{ccc}
I & H_t L & 0		\\ 		0 & \Phi_t L & G_t
\end{array} 	\right ] 	\left [ \begin{array}{ccc}
R_t & 0 & 0	\\ 		0 & \overline{D} & 0	\\
0 & 0 & Q_t 		\end{array} 	\right ]^{\frac{1}{2}} \Theta
= 	 }					\nonumber \\ & &
 	\left [ \begin{array}{ccc}
L_{e}		& 0	& 0	\\ 		\Phi_t
K_tL_{e}		& L_{+} & 0 		\end{array} 	\right
] 	\left [ \begin{array}{ccc} 		D_{e}	& 0 	&
0	\\ 		0 	& D_{+} & 0	\\ 		0
& 0 	& D_{a} 		\end{array} 	\right ]^{\frac{1}{2}}
\label{pre-ADThetaA'D'} \end{eqnarray} or, for a suitably defined
$\Psi$, \begin{eqnarray} 	\left [ \begin{array}{ccc}
I & H_t L & 0		\\ 		0 & \Phi_t L & G_t
\end{array} 	\right ] \Psi &=& 	\left [ \begin{array}{ccc}
L_{e}		& 0	& 0	\\ 		\Phi_t
K_tL_{e}		& L_{+} & 0 		\end{array} 	\right
].  \label{pre-APsiA'} 
\end{eqnarray} 
Note that $\Psi$ can be any orthogonal matrix that triangularizes the
pre--matrix on the left--hand side of Equation~\ref{pre-APsiA'}.  Using
fast Givens rotations \cite{gentleman73}, we can find $\Psi$ as a
product of elementary matrices, none of which contain any scalar
square roots.

Since the matrix $A'$ contains $L_{+}$ and the matrix $D'$ contains $D_{+}$, the 
matrix $P_{t|t}$ can be explicitly computed as soon as the triangularization
is complete. 
The state update given by Equation~\ref{time-x} requires $\Phi _tK_t$.
Since the matrix $A'$ contains both $L_{e}$ and
$\Phi_tK_t L_{e}$, the required matrix can be computed by 
  \[ \Phi_t K_t = [ \Phi_tK_t L_{e}][L_{e}]^{-1}. \]
Note that $L_e$ is a unit lower triangular matrix.   

\section{Fast Givens Rotations}
\label{sec:givens}
In this section we show how to use fast Givens rotations to find
the appropriate matrix $\Psi$ in Equation~\ref{pre-APsiA'}.
The fast Givens rotation, due to Gentleman\cite{gentleman73}, eliminates the
necessity for computing square roots.  Therefore, it may be computationally
efficient in certain implementations of the factorized filter.

Consider a matrix factorization
  \begin{equation}
	  S = ADA^{T},		\label{M-ADA}
  \end{equation}
where $D = diag[d_1,d_2,\ldots,d_n]$ is a diagonal matrix.
The goal is to transform the factorization in
Equation~\ref{M-ADA} to the factorization
  \begin{equation}
	  S = A^{'}D^{'}(A^{'})^{T},		\label{M-A'D'A'}
  \end{equation}
where $D^{'}= diag[d_1^{'}, d_2^{'} ,\ldots, d_n^{'}]$ is a diagonal matrix and $A^{'}$ is unit lower triangular.
As a simple case first consider the
use of a fast Givens algorithm to rotate
a $1 \times n$ vector 
$A = [a_{1},a_2, \ldots ,a_{n}]$, resulting in a $1 \times n$
vector $A^{'} = [1 \;\; 0 \; \ldots \; 0]$.
This can be accomplished by zeroing out elements
2 through $n$ of $A$ one at a time.  


Consider the $n \times n$ matrix $\Psi _{i,j}$ as follows:
  \begin{equation}
     \Psi _{i,j} =
	\left [ \begin{array}{ccccrcc}
		1 & \cdots & 0 & \cdots & & \cdots & 0 \\
		\vdots & \ddots &&&&& \vdots \\
		0 &  & a_{i}d_{i}/d^{'}_{i} & & - a_{j} & &  \\
		\vdots &&& \ddots &&& \vdots \\
		&& a_{j}d_{j}/d^{'}_{i} & & a_{i} && \\
		\vdots &&&&& \ddots & \vdots \\
		0 & \cdots && \cdots && \cdots & 1
		\end{array}
	\right ],
     \label{Psiij}
  \end{equation}
where
  \begin{equation}
	d^{'}_{i} = a_{i}^{2}d_{i} + a_{j}^{2}d_{j}.
     \label{d'i}
  \end{equation}
Matrix $\Psi _{i,j}$ differs from the identity matrix in four 
places --- the intersections of rows $i$ and $j$ with
columns $i$ and $j$.
It is seen that
  \begin{equation} 
	[a_{1} \ldots a_{i} \ldots a_{j} \ldots a_{n}] \Psi _{i,j}  =
	[a_{1} \ldots 1 \ldots 0 \ldots a_{n}],
     \label{APsiA}
  \end{equation}
where the $jth$ element becomes zero,
the $i$th element becomes unity, and all other elements of $A$ remain
unchanged.  Let
  \begin{equation}
    d^{'}_{j} = \frac{d_{i}d_{j}}{d^{'}_{i}}		\label{d'j}
  \end{equation}
and 
  \begin{equation}
	d^{'}_{k} = d_{k}				\label{d'k}
  \end{equation}
for integers $1 \leq k \leq n$ such that $k \neq i$ and $k \neq j$.
It can be verified that 
 \[ \Theta _{i,j} = D^{- \frac{1}{2}} \Psi _{i,j} (D^{'})^{\frac{1}{2}} \]
results in an orthogonal matrix.  Therefore, by computing the matrix 
$D^{'}$ as specified above, the $j$th element of row vector $A$
is set to zero.

The above procedure to zero out 
elements in a $1 \times n$ vector can be extended to triangularize 
an  $m \times n$ matrix with $m \leq n$.
The strategy is to zero out the rows of $A$ one at a time to obtain
a unit lower triangular matrix.

Let $A_{i}$ and $A'_{i}$ represent the $i$th row of the $m \times n$ matrices 
$A$ and $A'$, respectively.
The first step is to zero out entries in $A_{1}$, where
the matrix $A$ is multiplied by $\Psi _{1,j}$ as
  \[
	\left [ \begin{array}{c}
		A_{1} \\
		A_{2} \\
		\vdots  \\
		A_{m} 
		\end{array}
	\right ]
	\Psi _{1,j} =
	\left [ \begin{array}{c}
		A^{'}_{1} \\
		A^{'}_{2} \\
		\vdots  \\
		A^{'}_{m} 
		\end{array}
	\right ].
  \]
The equation $A_{1} \Psi _{1,j} = A^{'}_{1}$ is Equation~\ref{APsiA}
with $i=1$.
The remaining row transformations can be represented by
  \begin{equation}
	[a_{k1} \ldots a_{kj} \ldots a_{kn}] \Psi _{1,j} =
	[a^{'}_{k1} \ldots a^{'}_{kj} \ldots a_{kn}],
     \label{AkPsiAk}
  \end{equation}
where $a_{kj}$ and $a'_{kj}$ represent the $jth$ elements of $A_{k}$
and $A'_{k}$, respectively.  Only elements
1 and $j$ of $A_{k}$ change and are given by
  \begin{equation}
	a^{'}_{k1} = \frac{a_{11}d_{1}}{d^{'}_{1}}a_{k1} +
			\frac{a_{1j}d_{j}}{d^{'}_{1}}a_{kj}
     \label{pre-a'k1}
  \end{equation}
  \begin{equation}
	a^{'}_{kj} = - a_{1j}a_{k1} + a_{11}a_{kj}.
     \label{pre-a'kj}
  \end{equation}

For $A_{i}$, we need to compute $\Psi _{i,j}$ for all $j > i$ 
in order to complete the matrix triangularization.  
When zeroing $A_{i}$, Equations~\ref{pre-a'k1}~and~\ref{pre-a'kj}
become
  \begin{equation}
	a^{'}_{ki} = \frac{a_{ii}d_{i}}{d^{'}_{i}}a_{ki} +
			\frac{a_{ij}d_{j}}{d^{'}_{i}}a_{kj}
     \label{a'ki}
  \end{equation}
  \begin{equation}
	a^{'}_{kj} = - a_{ij}a_{ki} + a_{ii}a_{kj}.
     \label{a'kj}
  \end{equation}
Note that all rotations of $A_{i}$ affect all other rows.  
However, each row above $A_{i}$
contains zeroes as $i$th and $j$th elements, and as evidenced by
Equations~\ref{a'ki}~and~\ref{a'kj}, these elements will remain as zeroes
for any rotation applied to them.  

A procedure to triangularize the $m \times n$ matrix $A$ using fast
Givens rotations is summarized in
Algorithm~\ref{givens-mxn}.  Row $A_{1}$ is first zeroed out, then $A_{2}$,
and so on until a unit lower triangular matrix is achieved by zeroing
out all rows of $A$.
The order in which the elements of $A$ are zeroed out can affect the 
efficiency of the algorithm.

\subsection{\bf  Efficiency Improvement to 
		 Algorithm~\protect\ref{givens-mxn}  } 
%\subsection{Efficiency Improvement to Algorithm}

The computation of $a^{'}_{ki}$ and $a^{'}_{kj}$ using 
Equations~\ref{a'ki}~and~\ref{a'kj} requires four multiplications and
two additions for each row $A_{k}$ being rotated.  However, after the
first rotation zeroing an element of row $A_{i}$, notice that 
$a_{ii} = 1$.  Hence, for all subsequent rotations zeroing elements in
$A_{i}$, Equation~\ref{a'kj} can be simplified to
  \begin{equation}
	a^{'}_{kj} = - a_{ij}a_{ki} + a_{kj}.
     \label{a'kj-opt}
  \end{equation}
Using $a_{ii} = 1$ and the results of Equation~\ref{a'kj-opt}, 
Equation~\ref{a'ki} can be written as follows \cite{gentleman73}:
  \begin{eqnarray}
	a^{'}_{ki} &=& \frac{d_{i}}{d^{'}_{i}}a_{ki} +
			\frac{a_{ij}d_{j}}{d^{'}_{i}}a_{kj} \nonumber \\
%		   &=& \frac{d_{i}}{d^{'}_{i}}a_{ki} +
%			\frac{a_{ij}d_{j}}{d^{'}_{i}}(a^{'}_{kj}+a_{ij}a_{ki})
%							 \nonumber \\
		   &=& \frac{d_{i}+a_{ij}^{2}d_{j}}{d^{'}_{i}}a_{ki} +
			  a_{ij} \frac{d_{j}}{d^{'}_{i}} a^{'}_{kj}
							 \nonumber \\
		   &=& a_{ki} + (a_{ij} \frac{d_{j}}{d^{'}_{i}}) a^{'}_{kj}.
     \label{a'ki-opt}
  \end{eqnarray}
Since $a_{ij}d_{j}/d^{'}_{i}$ is an element of the
transformation matrix $\Psi _{i,j}$, this quantity need not be
recalculated in Equation~\ref{a'ki-opt}.  Hence, 
Equations~\ref{a'ki-opt}~and~\ref{a'kj-opt} compute the values of
$a^{'}_{ki}$ and $a^{'}_{kj}$, respectively, 
using only two multiplications and two additions.

A triangularization procedure that employs these more efficient 
equations is presented in \mbox{Algorithm~\ref{givens-mxn-opt}}.
The 
improved algorithm is more complicated than \mbox{Algorithm~\ref{givens-mxn}}
due to the separation of the first zeroing rotation of a row from
the rest of the rotations, but the improvement in performance can
be significant.  In fact, Algorithm~\ref{givens-mxn-opt}
triangularizes an $n \times n$ matrix ($m=n$) in order $\frac{2}{3} n^{3}$
operations\cite{golub83}, compared with the order $\frac{4}{3} n^{3}$
operations required by Algorithm~\ref{givens-mxn}. 

\section{Target Tracking Application}
\label{sec:trackingfilter}
In this section we describe a Kalman filter for target tracking
applications.  It will serve as a motivating example in the sections
to follow. The filter is used to estimate the position, velocity and
acceleration of a maneuvering target from noisy measurements of
range, azimuth, and elevation from a single sensor.
The tracking filter is based on the Newtonian mechanics and a
stochastic acceleration model.  In this formulation, it is assumed
that the target moves at a constant velocity if there is no maneuver
or atmospheric turbulence.  Evasive maneuvers and other disturbances
are modeled as perturbations upon a constant-velocity trajectory.

The target acceleration is assumed to be correlated in time.  For
example, a lazy turn will use a highly correlated acceleration.  On
the other hand, evasive maneuvers will be less correlated in time.
The correlation function $ \gamma ( \tau ) $ for the target
acceleration $ a(t) $ can be represented by
\begin{displaymath}
\gamma ( \tau ) = E\{ a(t)a(t+ \tau ) \} = \sigma _ m  ^ 2 e ^{ - \alpha \mid \tau \mid }
\end{displaymath}
where $E$ denotes the expected value, $ \sigma _ m  ^ 2 $ denotes the variance
of the target acceleration, and $ 1 / \alpha $ is the maneuver time constant.
The target dynamic model in one dimension  has the form\cite{singer70}:
\begin{displaymath}
d/dt  \left [
       \begin{array}{c}
        x \\
        \dot x \\
        \ddot x
       \end{array} \right ] =
\left [ \begin{array}{c c c}
      0 & 1 & 0 \\
      0 & 0 & 1 \\
      0 & 0 & - \alpha
     \end{array}
      \right ]  \left [ \begin{array}{c}
         x \\
         \dot x \\
         \ddot x \end{array}
    \right ]
 + \left [ \begin{array}{c}
    0 \\
    0 \\
    1 
 \end{array}
 \right ]   
 w(t)    
\end{displaymath}
where $ x,  \dot x, $ and $ \ddot x $ denote the target position, velocity,
and acceleration at time $ t $, and  $ \{ w(t) \} $ is a continuous
white gaussian
noise process with zero mean and variance density $ \sigma _ m ^ 2 $.  Note that
the target acceleration is represented by a first-order
Gauss-Markov process. The model is clearly a simple approximation of the vehicle acceleration 
dynamics.  However, in active tracking modes, the measurements may be
available at high data rates, and a simple acceleration model is useful.

The filtering may be performed in polar coordinates with range, azimuth, 
and elevation angles, and their derivatives as the state variables.  In 
polar coordinates the measurement system model is linear and uncoupled.
However, even for a target moving in a straight line at constant speed,
the dynamic model is nonlinear.
In Cartesian coordinates, a constant-speed target does obey
a linear law, but the measurement equations are nonlinear.  Both 
coordinate systems have been successfully used in many target tracking
applications.
A particular choice of the coordinate system depends on many other factors
and specific situations.

A discrete--time equation of target motion in Cartesian coordinates
can be expressed by a state variable model given by \cite{baheti86} 
\begin{eqnarray}
{\bf x}_{t+1} & = & \Phi ( \Delta t, {\bf \Lambda })  {\bf x}_{t} + {\bf w}(t) \nonumber \\
{\bf y}_t   & = & {\bf h}_t({\bf x}_t) + {\bf v}_t, ~~t=0, 1, \ldots \label{model}
\end{eqnarray}
where $t$ is the normalized discrete-time, $ \Delta t $ denotes the
sampling interval, $ {\bf x}_t $ is a nine-dimensional state vector with position,
velocity, and acceleration in each of the Cartesian coordinates axes $x,y,z$.
The first--order Markov parameters of the discretized stochastic acceleration models in
$ x,y,$ and $ z $ axes are denoted by the $  3 \times 3  $ matrix
$ {\Lambda } $ = diag$[ \lambda _ x , \lambda  _ y , \lambda _ z ]$,
where $ \lambda  = (1 - e ^ {- \alpha \Delta t}).  $
The matrix $ \Phi ( \Delta t, {\Lambda } )$ is given by
\begin{equation}
\Phi ( \Delta t, {\Lambda } ) = \left [
\begin{array}{c c c}
    \Phi _x & 0 & 0  \\
    0 &  \Phi _y       & 0 \\
    0 &  0          &  {\Phi _z }
\end{array}
\right ]
\label{eqnphi}
\end{equation}
where $ \Phi _x $ is a $ 3 \times 3 $ matrix given by
\begin{displaymath}
\Phi _x = \left [
\begin{array} {c c c}
1  & \Delta t  & 1/2 (\Delta t) ^ 2 \\
0  &  1        & \Delta t \\
0  &  0        & \lambda _x
\end{array}
\right ]
\end{displaymath}
The symbols $ \Phi _x,\Phi _y, \Phi _z $ denote the state transition matrices
for $x,y$ and $z$ axes respectively.

The state vector
$ {\bf x}(t) $ is represented by

\begin{displaymath}
{\bf x}_t = [ x, \dot x,\ddot x,y,\dot y,\ddot y ,z,\dot z,
\ddot z ] ^ T
\end{displaymath}
The vector $ {\bf h}_t({\bf x}_t) $
is the transformation from the Cartesian coordinates to the polar
coordinates defined by   
\begin{eqnarray}
 r & = & \sqrt{x ^ 2 + y ^ 2 + z ^ 2} \nonumber \\
 \theta _ T & = & \arctan (y/x) \nonumber \\
 \theta _ E & = & \arcsin (z/r) \nonumber
\end{eqnarray}
where $ r, \theta _ T, $ and $ \theta _ E $ denote the target range,
azimuth angle, and elevation angle, respectively.  The noisy 
measurements of range, azimuth, and elevation angles are denoted by
$ {\bf y}_t $.
\begin{displaymath}
{\bf y}_t = [r, \theta _ T, \theta _ E ] ^ T.
\end{displaymath}
It is assumed that the state and measurement noise sequences $ \{ {\bf w}(t) \}$ and $\{ {\bf v}(t) \} $ are zero mean gaussian random vectors
with covariance $ Q_t $ and $ R_t $, respectively.  The 
covariance $Q_t$ is given by
\begin{equation}
Q_t= \int_{t-1}^{t} \Phi (\tau) Q_c(t-\tau ) \Phi ^T (\tau) d \tau
\end{equation}
where $Q_c$ denotes the continuous--time process noise covariance density
in the stochastic acceleration model.
In the derivation of the discrete--time model from the continuous--time
a constant sampling rate is assumed for simplicity.
It is straightforward to include 
variable data rates, missing or bad quality observations, and variable
dimensions of the state and the measurement vectors.

The  measurement errors in the target range, azimuth angle, and
elevation angle are assumed uncorrelated and may depend on the target
range.  Let $ \sigma ^ 2 _ R , \sigma ^ 2 _ {\theta T }, $ and
$ \sigma ^ 2 _ {\theta E} $ denote the variance of the measurement
noise in range, azimuth angle, and elevation angle, respectively.  The
measurement noise covariance matrix $ R_t $ is given by
\begin{displaymath}
R_t = diag [ \sigma ^ 2 _ R, \sigma ^ 2 _ {\theta T}, \sigma ^ 2 _ {\theta E} ].
\end{displaymath}
A selection of the initial state ${\bf x}_0$ and covariance $P_0$ is
based on the first few observations of the target position and 
the measurement noise covariance $R_0$.  When in doubt, a large value of
$P_0$ is normally chosen, since after a few iterations the estimates 
are independent of the initial $P_0$.

The dynamics and measurement models presented above are commonly used in
many tracking and navigation applications.  An extended Kalman filter
can be used to recursively estimate the state vector $ {\bf x}_t $ 
using the dynamics and the measurement models in Equation \ref{model}.
The nonlinear measurement model  is 
linearized  at the state estimate $ {\bf \hat{x}}_{t|t-1} $ 
so that it can be used at time $t$ in a measurement update. 
The linearized observation matrix $H_t$ is given by
\begin{displaymath}
H = 
\left[
\begin{array}{ccccccccc}
{x}/{r} &0&0&{y}/{r} &0&0&{z}/{r} &0&0 \\
{-y}/{d^2} &0&0&{x}/{d^2} &0 &0  &0 &0  &0 \\
{-xz}/{r^{2}d} &0&0 &{-yz}/{r^{2}d} &0&0 &{d}/{r^2} &0  &0\\
\end{array}
\right]
\label{eqnjacstruct}
\end{displaymath}
where $d=(x^2 + y^2)^{1/2}$ and $r=(x^2 + y^2 + z^2)^{1/2}$.

\section{Mapping on Linear Array}
\label{sec:mapping}
In this section, we describe general methods for mapping the standard
Kalman filter onto a linear array of $n+2$ cells, and for mapping the
factorized Kalman filter onto a linear array of $n+m+1$ cells, where
$n$ is the size of the state vector, and $m$ is the size of the
measurement vector. First, we describe a simple model of a linear
array.

\subsection{Machine Model}
\label{sec:machinemodel}
The machine model is the linear array depicted in
Figure~\ref{fig:array}.  A square denotes one of $p$ logical processors,
which we call {\em logical cells}.  Arrows denote unidirectional
communications channels. Logical cells are numbered from left to
right, starting with unity.  Logical cells are separately programmable
and asynchronous.  Each logical cell has its own local code and data
storage.  Each logical cell runs at most one process, which is
assigned to the cell at compile time. A communications channel is a
bounded FIFO buffer. Each logical cell can send and receive data from
its immediate right and left neighbors.  Logical cell $1$ is connected
on the left to an external host.  Logical cell $k$ sends a word of
data to logical cell $k+1$ using a statement of the form
\begin{displaymath}
\mbox{\tt {\bf send}(R, val)}
\end{displaymath}
which blocks until the buffer between cells $k$ and $k+1$ is nonfull,
and then places the local variable {\tt val} at the tail of the
buffer.  Similarly, logical cell $k+1$ receives a word of data from
logical cell $k$ using a statement of the form
\begin{displaymath}
\mbox{\tt {\bf receive}(L, val)}
\end{displaymath}
which blocks until the buffer between cells $k+1$ and $k$ is nonempty,
and then retrieves the first item from the buffer and assigns it to
variable {\tt val}. 

\subsection{Standard Kalman Filter}
In this section, we describe a general method for mapping a standard
Kalman filter onto a linear array of $n+2$ logical cells, where
$n$ is the size of the state vector.  The result was first reported in
\cite{baheti88e}. An earlier, less efficient mapping was reported in
\cite{ohallaron88d}. First we discuss some general design guidelines.
Then we give a detailed description of the mapping.

The first design guideline is to devote the bulk of the logical cells
in the array to the covariance updates.  It is well known that the
covariance updates are the dominant computations, so it would seem
obvious to allocate the majority of the computational resources to
these computations.  The second design guideline is to distribute the
columns of the covariance matrix across the logical cells of the
linear array and to update these columns in place. The motivation here
is to minimize the data transfers involving the largest data
structure.  The third design guideline is to minimize the number of
data transfers between the linear array and the external host. The
motivation here is that on real systems, the transfer rate between the
external host and the linear array is typically an order of magnitude
less than the transfer rate between adjacent logical cells.

The mapping of the standard Kalman filter onto a linear array is shown
in Figure~\ref{fig:skfmapping}.  The mapping assumes a linear array of
$n+2$ logical cells: a collection of $n$ {\em p-cells}, a single {\em
k-cell}, and a single {\em x-cell}.
The $kth$ p-cell contains the $kth$ column of the covariance
matrix and is responsible for performing the time and measurement
updates of that column (Equations \ref{eqntup} and \ref{eqnmupk}). 
The k-cell computes Kalman gain vectors (Equation \ref{eqngaink}),
which are distributed to the x-cell and the array of p-cells, and
Jacobian transformations, which are distributed to the array of
p-cells.  The x-cell contains the state vector, ${\bf \hat{x}}$, and
is responsible for updating the state vector (Equations \ref{eqntux}
and \ref{eqnmuxk}), and computing the coordinate transformations.
Details on the computations performed by the x-cell and the k-cell can
be found in \cite{ohallaron88d}. In this paper, we will focus on the
covariance updates performed by the p-cells.

{\em Covariance Time Update: } Let the covariance matrix $P_{t|t}$
be written as
\begin{equation}
P_{t|t} = \left[ P_1 ~ P_2 ~\cdots~ P_n \right]
\label{eqnmap}
\end{equation}
where $P_k$ denotes the $kth$ column of $P_{t|t}$. The time
subscripts for $P_k$ are dropped to simplify the notation.

For the target tracking filter, we can exploit the sparse structure of
$\Phi_t$ to develop a simple parallel algorithm for the covariance
time update.  The algorithm consists of two steps: shifting and
computing.  Initially, the $kth$ p-cell contains the $kth$
column of $P_{t|t}$:
\begin{center}
\begin{tabular}{rccccccccc}
Cell:&1&2&3&4&5&6&7&8&9\\ 
Contents:&$P_1$&$P_2$&$P_3$&$P_4$&$P_5$&$P_6$&$P_7$&$P_8$&$P_9$\\ 
\end{tabular}
\end{center}
First, columns $P_2$, $P_3$, $P_5$, $P_6$, $P_8$, $P_9$ are shifted
(in parallel) one position to the left.  After the shift the array is
in the following state:
\begin{center}
\begin{tabular}{rccccccccc}
Cell:&1&2&3&4&5&6&7&8&9\\ 
Contents:&$P_1$&$P_2$&$P_3$&$P_4$&$P_5$&$P_6$&$P_7$&$P_8$&$P_9$ \\
&$P_2$&$P_3$&     &$P_5$&$P_6$&     &$P_8$&$P_9$&      \\
\end{tabular}
\end{center}
Next, columns $P_3$, $P_6$, and $P_9$ are again shifted one position
to the left. After the shift the array is in the following state:
\begin{center}
\begin{tabular}{rccccccccc}
Cell:&1&2&3&4&5&6&7&8&9\\ 
Contents:&$P_1$&$P_2$&$P_3$&$P_4$&$P_5$&$P_6$&$P_7$&$P_8$&$P_9$ \\
&$P_2$&$P_3$&     &$P_5$&$P_6$&     &$P_8$&$P_9$&      \\
&$P_3$&     &     &$P_6$ &    &    &$P_9$  &    &      \\
\end{tabular}
\end{center}

No further shift operations are required due to the block-diagonal
structure of $\Phi_t$.  With the appropriate columns in place, each cell
computes the columns of $P_{t+1|t}$ in parallel according to Equation
\ref{eqntup}.

Note that while this mapping ($kth$ column of the covariance matrix
$P_{t|t}$ on the $kth$ p-cell) is general, the shift pattern may vary
depending on the specific structure of $\Phi_t$. For example, the
structure of $\Phi_t$ depends on the ordering of elements in the state
vector.  Different orderings will result in different shift patterns.
We purposely do not exploit the symmetry of $P$.  This choice allows a
very regular and balanced algorithm for the measurement update presented
in the following.

{\em Covariance Scalar Measurement Update:}
The measurement update consists of the innovations covariance, 
Kalman gain, and the estimation error covariance matrix given
by Equations \ref{eqninnovcovk}, \ref{eqngaink}, \ref{eqnmupk}
respectively. Again, the strategy is to update the columns of the
covariance matrix $P_{t|t}$ in parallel. Note that the $kth$
column of $P_{t|t-1}$ is in the $kth$ p-cell after the time update. 

Define 
\begin{equation}
P_{t|t-1} = \left[ P_1^{(0)} ~ P_2^{(0)} ~\cdots~ P_n^{(0)} \right]
\label{eqnblockp}
\end{equation}
where $P_k^{(0)}$ denotes the $kth$ column of $P_{t|t-1}$.
The time subscripts for $P_k^{(0)}$ are dropped to simplify
the notation. Note that $H_i$ denotes the $ith$ row of the 
measurement partial derivative matrix $H_t$. It is easy to
see that
\begin{equation}
H_{i}P_{t|t-1} = \left[ H_{i}P_{1}^{(0)} ~ H_{i}P_{2}^{(0)} ~\cdots~ H_{i}P_{n}^{(0)} \right]. 
\label{eqnblockhp}
\end{equation}

Initially, before the $ith$ $(1 \leq i \leq m)$ scalar measurement
update of $P$, the $kth$ p-cell contains $P_k^{(i-1)}$. The superscript
$(i-1)$ is now dropped to simplify the notation.
\begin{center}
\begin{tabular}{rcccc}
Cell:&1&2&$\cdots$&$n$\\
Contents:&$P_1$&$P_2$&$\cdots$&$P_n$  \\
\end{tabular}
\end{center}
The scalar elements of $H_i$ are sent by the k-cell to the array of
p-cells.
\begin{displaymath}
\cdots ~H_{i2} ~H_{i1} \rightarrow \framebox{1} \rightarrow \framebox{2} 
\rightarrow \cdots \rightarrow \framebox{n}
\end{displaymath}
As the elements of $H_i$ flow through the $kth$ p-cell, the cell
computes $H_i P_k$ and passes the result back through the array
to the k-cell.
\begin{displaymath}
H_{i}P_1 ~H_{i}P_2 ~\cdots \leftarrow \framebox{1} \leftarrow
\framebox{2}
\leftarrow \cdots \leftarrow \framebox{n}
\end{displaymath}
After this computation, the array of p-cells is in the following
state.
\begin{center}
\begin{tabular}{rcccc}
Cell:&1&2&$\cdots$& $n$ \\
Contents:&$P_1$&$P_2$&$\cdots$&$P_n$ \\
&$H_i P_1$&$H_i P_2$&$\cdots$&$H_i P_n$ \\
\end{tabular}
\end{center}
Next, the k-cell uses $H_i P$ to compute the gain vector
$K^{(i)}$. As each element of $K^{(i)}$ is computed it is sent 
by the k-cell to the array of p-cells.
\begin{displaymath}
\cdots ~K_2^{(i)} ~K_1^{(i)} \rightarrow \framebox{1} \rightarrow \framebox{2} 
\rightarrow \cdots \rightarrow \framebox{n}
\end{displaymath}
As $K^{(i)}$ flows through the array, the $kth$ p-cell uses
$K^{(i)}$ and $H_{i} P_k^{(i-1)}$ to perform the scalar measurement update of
$P_k$ according to the following equation.
\begin{equation}
P_k^{(i)}  =  P_k^{(i-1)} - K^{(i)}H_{i}P_k^{(i-1)}, ~~k = 1,...,n .
\label{eqnmupcol}
\end{equation}

\subsection{Factorized Kalman Filter}
In this section, we describe a general method for mapping a factorized
Kalman filter onto a linear array of $n+m+1$ logical cells, where
$n$ is the dimension of the state vector, and $m$ is the dimension of
the measurement vector.  The result was first reported in
\cite{baheti89}.  First, we discuss some general design guidelines,
and then give a detailed description of the mapping.

The design guidelines for the factorized filter are similar to
those for the standard Kalman filter. The first guideline is to allocate the
bulk of the computational resources to the combined covariance update,
which is the dominant computation.  The second guideline is
to distribute the rows of the pre-matrix across the logical cells of
the linear array and to update these rows in place. The third
guideline is to minimize the number of data transfers between the
linear array and the external host.

The mapping of the factorized Kalman filter onto a linear array is
shown in Figure~\ref{figsqr-array}.  The mapping assumes a linear
array of $n+m$ {\em a-cells} and a single {\em x-cell}.
The $kth$ p-cell contains the $kth$ row of the pre--matrix
and is responsible for performing the combined covariance update on
that row (Equation \ref{pre-APsiA'}).  The x-cell performs the combined state
update in a straightforward manner using Equation~\ref{time-x}. In the
following we focus on the combined covariance update performed by 
the a-cells.

{\em Combined Covariance Update:} 
A combined update of the covariance matrix $P_{t|t-1}$ is performed using
fast Givens rotations of the matrices
\begin{equation}
A = 
	\left [ \begin{array}{ccc}
		I & H_t L & 0		\\
		0 & \Phi_t L & G_t
		\end{array}
	\right ]
\label{eqncombined1}
\end{equation}
and
\begin{equation}
 D = 
	\left [ \begin{array}{ccc}
		R_t & 0 & 0	\\
		0 & \overline{D} & 0	\\
		0 & 0 & Q_t
		\end{array}
	\right ]
\label{eqncombined2}
\end{equation}
where $A$ is an $(n+m) \times (2n+m)$ {\em pre-matrix} and $D$ is a
$(2n+m) \times (2n+m)$ diagonal matrix.
For notational convenience, define $A$ as an $M \times N$ matrix.
In order to facilitate the development of the mapping for
fast Givens rotations, it is useful to represent 
Algorithm~\ref{givens-mxn-opt} 
in a condensed manner, shown in Algorithm~\ref{alggivens-short}.
Below is a textual description of what each procedure listed in 
Algorithm~\ref{alggivens-short} represents:
\begin{description}
\item[$W_{i,N}$:]  
	Rotate $A_{i}$ to zero out the $N$th element
	and compute resulting $A'_{i}$ and $D'$.
\item[$X_{i,N,k}$:]  
	Rotate row $A_{k}$ of matrix $A$,
	changing only elements $i$ and $N$.
\item[$Y_{i,j}$:]  
	Rotate $A_{i}$ to zero out the $j$th element (when $a_{ii}=1$)
	and compute resulting $A'_{i}$ and $D'$.
\item[$Z_{i,j,k}$:]  
	Rotate row $A_{k}$ of matrix $A$ (when $a_{ii}=1$),
	changing only elements $i$ and $j$.
\end{description}



A summary of the
data used in the computation and the data changed by each procedure
is presented in Table~\ref{tabdata-used}.
The primes have been dropped from the notation because the
primed variables are moved back into the unprimed variables as soon as all
the computations are completed.

From the table, it is apparent that procedure $W_{i,N}$ can be
executed by a cell which stores only $D$ and row $A_{i}$ of matrix $A$.
This cell computes $\Psi _{i,N}$ and communicates it to other cells.

It is also apparent that procedure $X_{i,N,k}$ can be executed by a cell
which can access the rotation matrix $\Psi _{i,N}$ and which stores
only row $A_{k}$ of matrix $A$.  Hence, the procedures $X_{i,N,k}$ 
($i < k \leq M$) can be executed by separate cells in parallel.

The data details for procedures $Y_{i,j}$ and $Z_{i,j,k}$ are similar to
the details for procedures $W_{i,N}$ and $X_{i,N,k}$, respectively.
The only difference is notational, where the rotation matrix
is represented by $\Psi _{i,j}$ rather than $\Psi _{i,N}$.

The above analysis leads to a natural mapping of the algorithm onto a
linear array.  Consider the problem of zeroing row $A_{i}$ of matrix $A$.
What is desired is an array where rows $A_{i}$ through $A_{M}$ of the
pre--matrix are 
stored on a separate cells.  
%\footnote{Recall that, in the triangularization algorithm, zeroing
%row $A_{i}$ does not affect the rows above it in matrix $A$.  Therefore,
%this storage is not considered here.}
Cell $i$ contains
$A_{i}$ and $D$ and performs procedure $W_{i,N}$, computing
$\Psi _{i,N}$ in the process.  Then $\Psi _{i,N}$ is passed to
cells $i+1$ through $M$ so that cell $k$ can appropriately rotate $A_{k}$
using procedure $X_{i,N,k}$.
Subsequently, for each remaining element of $A_{i}$ to be zeroed, 
\begin{enumerate}
\item Cell $i$ performs procedure $Y_{i,j}$, 
\item Cell $i$ sends $\Psi _{i,j}$ to all cells $k$ ($k=i+1,\ldots,M$), and
\item Each cell $k$ performs procedure $Z_{i,j,k}$.
\end{enumerate}

The next step is to extend the above reasoning from the case of zeroing
one row to the case of zeroing all rows (Table~\ref{tabzero-matrix}).  
The first step of triangularizing
a matrix is zeroing the first row, so the above reasoning applies directly
$(i=1)$.  Hence, what is required is an array of
$n+m$ a-cells, where a-cell
$k$ contains row $A_{k}$ of pre--matrix $A$.
Initially, the matrix $D$ is stored on a-cell 1.  
As described in the previous paragraph, a-cell 1
performs procedures $W$ and $Y$; a-cells 2 through $n+m$ perform
procedures $X$ and $Z$.  
The subscripts have been dropped for convenience.

Once $A_{1}$ is zeroed, the next step is to zero $A_{2}$.  As shown in 
Table~\ref{tabzero-matrix}, this requires that a-cell 2 store matrix $D$.
Therefore, $D$ must be passed to a-cell 2 by a-cell 1 
before row $A_{2}$
can be zeroed.  Once this passage is complete, the zeroing process
follows the single row reasoning outlined above ($i$=2).  Hence, 
a-cell 2 now
performs procedures $W$ and $Y$, and a-cells 3 through $n$+$m$ perform
procedures $X$ and $Z$.  Since row $A_{1}$
will not be changed by subsequent rotations in the algorithm, a-cell 1
will be dormant.

This process continues for the zeroing of all rows, as shown in
Table~\ref{tabzero-matrix}.  Notice that zeroing $A_{n+m}$ does not require
the execution of procedures $X$ and $Z$ (there are no rows in $A$ below
row $A_{n+m}$).  When the triangularization process is complete, the
resulting matrix $A'$ (stored in $A$) is partitioned just as the original
matrix: row $k$ resides on a-cell $k$.  The new diagonal matrix $D'$ 
(stored in $D$) resides on a-cell $n$+$m$, 
whereas the initial diagonal matrix resided on a-cell 1.  
Thus, an algorithm is defined on an $n$+$m$-cell
linear array to compute a unit lower triangular $A'$ and a diagonal $D'$ 
such that Equation~\ref{M-A'D'A'}
is satisfied.

To set the $A$ matrix in Equation \ref{eqncombined1} for Kalman
filtering, the a-cells contain the following matrices:
\begin{center}
\begin{tabular}{rccccccc}
Cell:	& 1 & 2 & $\cdots$ & $m$ & $m$+1 & $\cdots$ & $m$+$n$		\\
Contents:& $Q_t,R_t$&&$\cdots$& & $\Phi _{1},G_{1}$ &$\cdots$& $\Phi _{n},G_{n}$\\
\end{tabular}
\end{center}
where $\Phi _{k}$ and $G_{k}$ represent the $k$th rows of $\Phi_t$ and $G_t$,
respectively.  The first a-cell initializes $D$, given by Equation
\ref{eqncombined2}, therefor it contains 
$Q_t$ and $R_t$.  In addition, the a-cells contain the rows of matrix $L$:
\begin{center}
\begin{tabular}{rccccccc}
Cell:	& 1 & 2 & $\cdots$ & $m$ & $m$+1 & $\cdots$ & $m$+$n$		\\
Contents: & & & $\cdots$ & & $L_{1}$ & $\cdots$ & $L_{n}$,		\\
\end{tabular}
\end{center}
where $L_{k}$ represents the $k$th row of $L$.

To begin the covariance update computation, matrix $H_t$ is received from
the x-cell one row at a time such that $H_{k}$, the $kth$ row of $H_t$,
is stored on the $kth$ a-cell.
\begin{displaymath}
\cdots ~H_{2} ~H_{1} \rightarrow \framebox{1} \rightarrow \framebox{2} 
\rightarrow \cdots \rightarrow \framebox{m}
\end{displaymath}
Then the rows of matrix $L$ are distributed to all of the a-cells so
that the contents of the cells are:
\begin{center}
\begin{tabular}{rccccccc}
Cell:	& 1 & 2 & $\cdots$ & $m$ & $m$+1 & $\cdots$ & $m$+$n$		\\
Contents: & $H_{1},L$ & $H_{2},L$ & $\cdots$ & 
				$H_{m},L$ & $L$ & $\cdots$ & $L$	\\
\end{tabular}
\end{center}
At this point, each cell has all the information it needs to compute its
row of matrix $A$.

Once pre--matrix $A$ is formed, we can triangularize the matrix by performing 
the appropriate fast Givens rotations, and the covariance update is 
complete.  However, it is useful to examine the triangularization process
in two stages:  zeroing the first $m$ rows of $A$ and zeroing the last
$n$ rows of $A$.  After completion of the first stage, the first $m$ rows
and columns of matrix $A$ will not be affected by subsequent rotations.
Since the matrices $L_{e}$ and $K_{t}L_{e}$ 
reside within the first $m$ columns of
the triangularized matrix $A'$, they are completely determined upon completion
of the first stage of the fast Givens process.  Hence, these matrices
may be sent to the x-cell before the second stage of the 
triangularization is performed.
%\begin{displaymath}
%L_{e} ~K_{t}L_{e} \leftarrow \framebox{1} \leftarrow
%\framebox{2}
%\leftarrow \cdots \leftarrow \framebox{n+m}
%\end{displaymath}

\section{Warp Implementations}
\label{sec:implementation}
In this section, we describe implementations of the standard and
factorized filters for target tracking applications on the Warp and
report the measured performance.

\subsection{The Warp Computer}
The Warp computer\cite{annaratone87e} (or simply Warp) is a linear
array of  identical and programmable cells connected to a
general-purpose host. Figure \ref{warp} shows the architecture of
a 10-cell Warp.  An MC6802--based workstation called the {\em master}
is connected via a VME repeater to a pair of {\em clusters}. Each
cluster consists of an MC68020 processor(P) with an MC68881 floating
point coprocessor, 3 Mbytes of data memory(M), 1 Mbyte of program
memory, and a switch(S), all connected by a local VSB bus. The
clusters are known collectively as the external host, so called
because they are external to the Master.

Each Warp cell contains a pipelined floating-point adder and
multiplier, 8K instruction words, and 32K data words. The clock cycle
of each cell is 200ns. Each functional unit can emit one result per
cycle, for a maximum performance of 5 MFLOPS per functional unit, 10
MFLOPS per cell, and 100 MFLOPS for the 10-cell Warp array. Each cell
is connected to its nearest neighbor by two 32-bit buffered data
channels (X and Y).  Data flows from left to right along the X
channel; the Y channel can be statically reconfigured to pass data in
either direction. Each cell can send and receive up to 20 million
32-bit data words to and from its neighboring cells each second.

Warp is programmed in a high-level Pascal-like language with {\em
send} and {\em receive} statements for cell--cell and host--array data
transfers.  The compiler inputs a source file, generates
microcode for the Warp array, and generates C code for the
clusters.  Warp is accessed with a programming environment
which provides a programmable shell for interactively running and
debugging source programs.

\subsection{Implementation of Tracking Filter}
The Warp implementation of the 9-state tracking filter using the
standard Kalman equations is shown in Figure~\ref{fig:skfwarp}.
Since the current Warp has only 10 physical cells, and since the
mapping assumes there are 11 logical cells (9 p-cells, 1 k-cell, and 1
x-cell), more than one logical cell had to be assigned to some of the
physical Warp cells. Through experimentation we found that the best
assignment is the one shown in Figure~\ref{fig:skfwarp}. Note that
this assignment uses only 8 of the 10 Warp cells. With this
assignment, p-cells 1, 4, and 7 are each assigned to a single Warp cell,
p-cells 2 and 3 are assigned to one Warp cell, p-cells 5 and 6
are assigned to one Warp cell, and p-cells 8 and 9 are assigned
to one Warp cell.

Since data may enter the Warp array only from the left hand
side and may exit the array only from the right hand side, the
implementation reverses the left-right orientation of the logical
cells, using cell 0 (corresponding to p-cells 8 and 9) to input
the measurement vector ${\bf y}_t$ and
to stream it along the X channel through the array to cell 9
(corresponding to the x-cell).  The updated state vector ${\bf
\hat{x}_{t|t}}$ exits the array from the right along the X channel.

The implementation of the 9-state tracking filter using the
factorized update is shown
in Figure~\ref{fig:fkfwarp}. An approach similar to the standard
filter is used, with the left--right orientation of the logical cells
reversed, and with multiple a-cells assigned to some Warp cells to
compensate for the undersized array.

\subsection{Performance}
The key performance metric for a Kalman filter is {\em throughput},
the number of measurements that can be processed each second.
We let $T_p$ denote the throughput achieved using $p$ Warp cells.
Two other metrics are of interest, {\em speedup} and {\em efficiency}.
The speedup achieved using $p$ Warp cells is given by
\begin{displaymath}
S_p = \frac{T_p}{T_1}.
\end{displaymath}
The efficiency achieved using $p$ Warp cells is given by
\begin{displaymath}
E_p = \frac{S_p}{p}.
\end{displaymath}
Measured values for throughput, speedup, and efficiency are given in
Table \ref{tab:performance}.  Throughput was measured using a built-in
timing facility on the Warp external host. We ran the filter using 300
measurements from a target trajectory. We noted the time
that the input cluster started executing before processing the first
sample and the time that the output cluster finished
executing after processing the last sample.  Dividing 300 by the
difference between the start and finish times gave the throughput measure.
The throughput measures include the time
required to transfer data between the external host and the Warp
array, as well as startup times for the external host and Warp array.
Note that we use 8 cells to implement the standard filter and 10
cells to implement the factorized filter.

The performance of the standard Kalman filter on Warp is encouraging.
The filter is an eclectic mix of matrix operations, ranging in
complexity from linear to cubic in the number of states. Coordinating
all of these various operations and still achieving a measured
efficiency of 87\% is very encouraging. There are a number of reasons
for the good performance of the standard Kalman filter on Warp. First,
we use a specific ordering of the state variables and exploit the
structure of $\Phi_t$ in the time update.  Second, using Bierman's
scalar measurement updates avoids computing a matrix
inverse in the gain computation.  Third, we keep
the overhead associated with data transfers small by updating the columns
of the covariance matrix in--place and by overlapping data transfers with
computation.

The performance of the factorized Kalman filter on Warp is less
encouraging. If we consider that the throughput for the factorized
Kalman filter on an MC68020--based workstation is roughly 20 Hz, then
the throughput of 1250 Hz for the 10--cell Warp is quite impressive.
Unfortunately, most of that performance is due to the fast cycle--time
of the Warp cell. Only a factor of two is attributable to parallelism
exploited by the implementation. The problem stems from the fact that
a cell is idle once it zeros its assigned row. This imbalance could be
mitigated somewhat by choosing either a much bigger problem (more
states) or a much smaller array (fewer cells). In that case, rows of
the pre--matrix could be interleaved across the cells of the array to
create a more balanced load.

The performance data in Table~\ref{tab:performance} illustrates an
interesting tension.  We see that the standard Kalman filter can be
implemented very efficiently. However, because of truncation errors,
the standard filter may not provide the necessary numerical
robustness, and we may want to use the factorized filter. But given
the mappings discussed here, the factorized filter is not nearly as
efficient as the standard filter. Thus, we have the interesting
trade--off between numerical robustness and efficiency.  If a parallel
machine has sufficient word--length, then there is a great incentive
to use the standard filter rather than the factorized filter. Note
that this tension does not arise to this degree in the implementation
of Kalman filters on sequential computers.

\section{Conclusions}
\label{sec:conclusions}
We have reported on the results of an investigation into the efficient
implementation of extended Kalman filters on a programmable
commercially--available parallel computer designed for real--time
applications.  We have described general methods for mapping the
standard Kalman filter and the factorized filter on linear
arrays, and we have used these methods to implement standard and
factorized versions of a 9--state tracking filter on the Warp
computer.

The performance results demonstrate for the first time the feasibility
of implementing the standard filter on a real parallel machine, but
also point to the need for further work in the design and
implementation of factorized filters. There are several areas that
need further attention.  We need to develop better mappings for
factorized filters running on real linear arrays, or we need to show
that such mappings are not possible. Further, with the advent of
programmable general--purpose two-dimensional arrays like iWarp
\cite{borkar88,borkar90}, we
need to investigate efficient mappings for both the standard and
factorized filters on such architectures.

\section*{Acknowledgments}
Thanks to Prof. H. T. Kung for graciously allowing us unlimited access
to the Warp computer at Carnegie Mellon, to Chuck Fiduccia for a
discussion that ultimately led to the algorithm for the standard
filter's covariance time update, and to the anonymous reviewers
for their thoughtful suggestions.

\bibliography{/usr/droh//bib/defs,/usr/droh/bib/refs,/usr/droh/bib/me}

\newpage

\section*{Figures}
\begin{figure}[h]
\setlength{\unitlength}{0.175in}
\begin{center}
\begin{picture}(20,5)(1,0)
\scriptsize
\put(5,1){\framebox(2,2){\shortstack{{cell}\\{\em 1}}}}
\put(8,1){\framebox(2,2){\shortstack{{cell}\\{\em 2}}}}
\put(11,1){\framebox(2,2){\shortstack{{cell}\\{\em 3}}}}
\put(14,1){\makebox(2,2){$\cdots$}}
\put(17,1){\framebox(2,2){\shortstack{{cell}\\{\em p}}}}
\put(3,2.5){\vector(1,0){2}}
\put(5,1.5){\vector(-1,0){2}}
%\put(4,3.1){\makebox(0,0){$y$}}
%\put(4,1.1){\makebox(0,0){$\hat{x}$}}
\multiput(7,2.5)(3,0){4}{\vector(1,0){1}}
\multiput(8,1.5)(3,0){4}{\vector(-1,0){1}}
\end{picture}
\end{center}
\caption{Linear array model}
\label{fig:array}
\end{figure}

\newpage

\begin{figure}[h]
\setlength{\unitlength}{0.175in}
\begin{center}
\begin{picture}(26,5)(0,0)
\scriptsize
\put(5,1){\framebox(2,2){x-cell}}
\put(8,1){\framebox(2,2){k-cell}}
\put(11,1){\framebox(2,2){\shortstack{{p-cell}\\{1}}}}
\put(14,1){\framebox(2,2){\shortstack{{p-cell}\\{2}}}}
\put(17,1){\makebox(2,2){$\cdots$}}
\put(20,1){\framebox(2,2){\shortstack{{p-cell}\\{$n$}}}}
\put(3,2.5){\vector(1,0){2}}
\put(5,1.5){\vector(-1,0){2}}
\put(4,3.1){\makebox(0,0){$y$}}
\put(4,1.1){\makebox(0,0){$\hat{x}$}}
\multiput(7,2.5)(3,0){5}{\vector(1,0){1}}
\multiput(8,1.5)(3,0){5}{\vector(-1,0){1}}
\end{picture}
\end{center}
\caption{Mapping of the standard Kalman filter}
\label{fig:skfmapping}
\end{figure}

\newpage

\begin{figure}[h]
\setlength{\unitlength}{0.175in}
\begin{center}
\begin{picture}(22,5)(0,0)
\scriptsize
\put(5,1){\framebox(2,2){x-cell}}
\put(8,1){\framebox(2,2){\shortstack{{a-cell}\\{1}}}}
\put(11,1){\framebox(2,2){\shortstack{{a-cell}\\{2}}}}
\put(14,1){\makebox(2,2){$\cdots$}}
\put(17,1){\framebox(2,2){\shortstack{{a-cell}\\{$n+m$}}}}
\put(3,2.5){\vector(1,0){2}}
\put(5,1.5){\vector(-1,0){2}}
\put(4,3.1){\makebox(0,0){${y}$}}
\put(4,1.1){\makebox(0,0){${\hat{x}}$}}
\multiput(7,2.5)(3,0){4}{\vector(1,0){1}}
\multiput(8,1.5)(3,0){4}{\vector(-1,0){1}}
\end{picture}
\end{center}
\caption{Mapping of the factorized filter}
\label{figsqr-array}
\end{figure}

\newpage

\begin{figure}[h]
\setlength{\unitlength}{0.2in}
\begin{center}
\begin{picture}(20,13)(0,0)
\scriptsize
%host
\put(9,11.5){\framebox(3,2){\shortstack{Master}}}
\put(10.5,9.5){\line(0,1){2}}
%external host
%vme global bus
\put(7,9.5){\line(1,0){7}}
%components
\put(7,8){\framebox(.5,1){P}}
\put(8,8){\framebox(.5,1){M}}
\put(9,7){\framebox(.5,2){S}}
\put(11.5,7){\framebox(.5,2){S}}
\put(12.5,8){\framebox(.5,1){M}}
\put(13.5,8){\framebox(.5,1){P}}
%connections to vme
\multiput(7.25,9)(1,0){3}{\line(0,1){.5}}
\multiput(11.75,9)(1,0){3}{\line(0,1){.5}}
%connections to vsb
\multiput(7.25,7.5)(1,0){2}{\line(0,1){.5}}
\multiput(12.75,7.5)(1,0){2}{\line(0,1){.5}}
%vsb local busses
\put(7,7.5){\line(1,0){2}}
\put(12,7.5){\line(1,0){2}}
%array
\put(1,2){\framebox(1,1){\shortstack{0}}}
\put(3,2){\framebox(1,1){\shortstack{1}}}
\put(5,2){\framebox(1,1){\shortstack{2}}}
\put(7,2){\framebox(1,1){\shortstack{3}}}
\put(9,2){\framebox(1,1){\shortstack{4}}}
\put(11,2){\framebox(1,1){\shortstack{5}}}
\put(13,2){\framebox(1,1){\shortstack{6}}}
\put(15,2){\framebox(1,1){\shortstack{7}}}
\put(17,2){\framebox(1,1){\shortstack{8}}}
\put(19,2){\framebox(1,1){\shortstack{9}}}
%from input cluster to array
\put(9.25,5){\line(0,1){2}}
\put(.5,5){\line(1,0){8.75}}
\put(.5,2.75){\line(0,1){2.25}}
\put(.5,2.75){\vector(1,0){.5}}
%from array to input cluster
\put(11.75,5){\vector(0,1){2}}
\put(20.5,5){\line(-1,0){8.75}}
\put(20.5,2.75){\line(0,1){2.25}}
\put(20,2.75){\line(1,0){.5}}
% X and Y channel
\multiput(2,2.75)(2,0){9}{\vector(1,0){1}}
\multiput(19,2.25)(-2,0){9}{\vector(-1,0){1}}
\put(10.35,2.85){X}
\put(10.35,1.75){Y}
%dashed boxes labeling clusters and array
\put(.75,1.5){\dashbox{.1}(19.5,2)}
\put(6.5,6.5){\dashbox{.1}(3.5,3.5)}
\put(11,6.5){\dashbox{.1}(3.5,3.5)}
\put(8.3,10.5){\makebox(0,0){Input Cluster}}
\put(12.8,10.5){\makebox(0,0){Output Cluster}}
\put(10.5,1){\makebox(0,0){Warp Array}}
\end{picture}
\end{center}
\caption{The Warp computer}
\label{warp}
\end{figure}

\newpage

\begin{figure}[h]
\setlength{\unitlength}{0.175in}
\begin{center}
\begin{picture}(30,5)(1,0)
\scriptsize
\put(5,1){\framebox(2,2){\shortstack{Cell\\0}}}
\put(6,4.25){\makebox(0,0){$\overbrace{}^{\shortstack{\mbox{p-cells}\\9,8}}$}}

\put(8,1){\framebox(2,2){\shortstack{Cell\\1}}}
\put(9,4.25){\makebox(0,0){$\overbrace{}^{\shortstack{\mbox{p-cell}\\7}}$}}

\put(11,1){\framebox(2,2){\shortstack{Cell\\2}}}
\put(12,4.25){\makebox(0,0){$\overbrace{}^{\shortstack{\mbox{p-cell}\\6,5}}$}}

\put(14,1){\framebox(2,2){\shortstack{Cell\\3}}}
\put(15,4.25){\makebox(0,0){$\overbrace{}^{\shortstack{\mbox{p-cell}\\4}}$}}

\put(17,1){\framebox(2,2){\shortstack{Cell\\4}}}
\put(18,4.25){\makebox(0,0){$\overbrace{}^{\shortstack{\mbox{p-cells}\\3,2}}$}}

\put(20,1){\framebox(2,2){\shortstack{Cell\\5}}}
\put(21,4.25){\makebox(0,0){$\overbrace{}^{\shortstack{\mbox{p-cell}\\1}}$}}

\put(23,1){\framebox(2,2){\shortstack{Cell\\6}}}
\put(24,4.25){\makebox(0,0){$\overbrace{}^{\mbox{k-cell}}$}}

\put(26,1){\framebox(2,2){\shortstack{Cell\\7}}}
\put(27,4.25){\makebox(0,0){$\overbrace{}^{\mbox{x-cell}}$}}

\put(3,2.5){\vector(1,0){2}}
\put(4,3.1){\makebox(0,0){$y_t$}}

\put(28,2.5){\vector(1,0){2}}
\put(29,3.1){\makebox(0,0){$\hat{x}_{t|t}$}}

\multiput(7,2.5)(3,0){7}{\vector(1,0){1}}
\multiput(8,1.5)(3,0){7}{\vector(-1,0){1}}
\put(7.3,2.8){X}
\put(7.3,0.9){Y}
\end{picture}
\end{center}
\caption{Implementation of standard Kalman filter}
\label{fig:skfwarp}
\end{figure}

\newpage

\begin{figure}[h]
\setlength{\unitlength}{0.175in}
\begin{center}
\begin{picture}(35,5)(1,0)
\scriptsize
\put(5,1){\framebox(2,2){\shortstack{Cell\\0}}}
\put(6,4.25){\makebox(0,0){$\overbrace{}^{\shortstack{\mbox{a-cell}\\12}}$}}

\put(8,1){\framebox(2,2){\shortstack{Cell\\1}}}
\put(9,4.25){\makebox(0,0){$\overbrace{}^{\shortstack{\mbox{a-cell}\\11}}$}}

\put(11,1){\framebox(2,2){\shortstack{Cell\\2}}}
\put(12,4.25){\makebox(0,0){$\overbrace{}^{\shortstack{\mbox{a-cell}\\10}}$}}

\put(14,1){\framebox(2,2){\shortstack{Cell\\3}}}
\put(15,4.25){\makebox(0,0){$\overbrace{}^{\shortstack{\mbox{a-cell}\\9}}$}}

\put(17,1){\framebox(2,2){\shortstack{Cell\\4}}}
\put(18,4.25){\makebox(0,0){$\overbrace{}^{\shortstack{\mbox{a-cell}\\8}}$}}

\put(20,1){\framebox(2,2){\shortstack{Cell\\5}}}
\put(21,4.25){\makebox(0,0){$\overbrace{}^{\shortstack{\mbox{a-cell}\\7}}$}}

\put(23,1){\framebox(2,2){\shortstack{Cell\\6}}}
\put(24,4.25){\makebox(0,0){$\overbrace{}^{\shortstack{\mbox{a-cells}\\6,3}}$}}

\put(26,1){\framebox(2,2){\shortstack{Cell\\7}}}
\put(27,4.25){\makebox(0,0){$\overbrace{}^{\shortstack{\mbox{a-cells}\\5,2}}$}}

\put(29,1){\framebox(2,2){\shortstack{Cell\\8}}}
\put(30,4.25){\makebox(0,0){$\overbrace{}^{\shortstack{\mbox{a-cells}\\4,1}}$}}

\put(32,1){\framebox(2,2){\shortstack{Cell\\9}}}
\put(33,4.25){\makebox(0,0){$\overbrace{}^{\mbox{x-cell}}$}}

\put(3,2.5){\vector(1,0){2}}
\put(4,3.1){\makebox(0,0){$y_t$}}
\put(34,2.5){\vector(1,0){2}}
\put(35,3.1){\makebox(0,0){$\hat{x}_{t|t}$}}

\multiput(7,2.5)(3,0){9}{\vector(1,0){1}}
\multiput(8,1.5)(3,0){9}{\vector(-1,0){1}}
\put(7.3,2.8){X}
\put(7.3,0.9){Y}
\end{picture}
\end{center}
\caption{Implementation of factorized filter}
\label{fig:fkfwarp}
\end{figure}

\newpage

\section*{Tables}
\begin{table}[h]
\begin{center}
\begin{tabular}{||c|c|c||}    \hline
Procedure & Data Used & Data Changed \\  \hline
$W_{i,N}$   & $A_{i}$, \, $D$ &  $A_{i}$, \, $D$, \, $\Psi _{i,N}$	\\
$X_{i,N,k}$ & $A_{k}$, \, $\Psi _{i,N}$	&       $A_{k}$		\\
$Y_{i,j}$     & $A_{i}$, \, $D$	&  $A_{i}$, \, $D$, \, $\Psi _{i,j}$	\\
$Z_{i,j,k}$   & $A_{k}$, \, $\Psi _{i,j}$	&       $A_{k}$		\\
	\hline
\end{tabular}
\end{center}
\caption{Data Used and Changed by Procedures in 
			Algorithm~\protect\ref{alggivens-short} }
\label{tabdata-used}
\end{table}

\newpage

\begin{table}[h]
%\scriptsize
%\small
\begin{center}
\begin{tabular}{|c|ccccc|}
	\hline
Row   & a-cell 1 & a-cell 2 & a-cell 3 & $\cdots$ & 
	 a-cell $M$ \\ \hline
 1 & $A_{1}$,$D$ & $A_{2}$ & $A_{3}$ & $\cdots$ & 
		 $A_{M}$ \\
      & $W$,$Y$     & $X$,$Z$ & $X$,$Z$ & $\cdots$ & $X$,$Z$ \\
	&&&&& \\
 2 & $A_{1}$ & $A_{2}$,$D$ & $A_{3}$ & $\cdots$ 
		 & $A_{M}$ \\
      &             & $W$,$Y$ & $X$,$Z$ & $\cdots$ &  $X$,$Z$ \\
	&&&&& \\
 3 & $A_{1}$ & $A_{2}$ & $A_{3}$,$D$ & $\cdots$ & 
		 $A_{M}$ \\
      &             &         & $W$,$Y$ & $\cdots$ & $X$,$Z$ \\
  $\vdots$&$\vdots$&$\vdots$&$\vdots$&&$\vdots$ \\
 $M$-1 & $A_{1}$ & $A_{2}$ & $A_{3}$ & $\cdots$ & 
		 $A_{M}$ \\
      &             &         &  &                 &      $X$,$Z$ \\
        &&&&& \\
 $M$ & $A_{1}$ & $A_{2}$ & $A_{3}$ & $\cdots$ & 
		 $A_{M}$,$D$ \\
      &             &         &  &                 &          $W$,$Y$ \\
		\hline
\end{tabular}
\end{center}
\caption{Data and execution pattern of linear array throughout the
	 fast Givens algorithm.  The first column lists the row currently
	 being zeroed, and the remaining columns show the data stored
	 (first row) and the procedures executed (second row) 
	 on the a-cells of the array.}
\label{tabzero-matrix}
\end{table}

\newpage

\begin{table}[h]
\begin{center}
\begin{tabular}{|c|c|c|}
\hline
\hline
&Standard Filter&Factorized Filter\\
\hline
\hline
Single Warp Cell (sequential)&1140 Hz& 625 Hz\\
\hline
Multiple Warp Cells (parallel)&8000 Hz&1250 Hz\\
\hline
Speedup&7&2\\
\hline
Efficiency&87\%&20\%\\
\hline
\end{tabular}
\end{center}
\caption{Tracking filter performance}
\label{tab:performance}
\end{table}

\newpage

\section*{Algorithms}
\begin{algorithm}[h]
\begin{tabbing}
For \= For \= For \= each $A_{k}$ below $A_{i}$ ($i < k \leq m$) \= \kill
For each row $A_{i}$ of $A$ ($i=1,\ldots,m$) \\
\> For each element $j$ of $A_{i}$ ($j=n,\ldots,i+1$) \\
\> \> Calculate $d^{'}_{i}$ \> \> [Equation \ref{d'i}] \\
\> \> Calculate $d^{'}_{j}$ \> \> [Equation \ref{d'j}] \\
\> \> Calculate matrix $\Psi _{i,j}$  \> \> [Equation \ref{Psiij}] \\
\> \> Calculate $A^{'}_{i}$, the rotation of $A_{i}$
					 \> \> [Equation \ref{APsiA}] \\
\> \> \\
\> \> For each $A_{k}$ below $A_{i}$ ($k=i+1,\ldots,m$) \\
\> \> \> Calculate $a^{'}_{ki}$ \> [Equation \ref{a'ki}] \\
\> \> \> Calculate $a^{'}_{kj}$ \> [Equation \ref{a'kj}] \\
\> \> \> Reset $A_{k}$ to $A^{'}_{k}$ \> \\
\> \> EndFor \\
\> \> \\
\> \> Reset $A_{i}$ to $A^{'}_{i}$ \> \\
\> \> Reset $D$ to $D^{'}$ \> \\
\> EndFor \\
EndFor 
\end{tabbing}
\caption{Fast Givens rotations: $m \times n$ matrix}
\label{givens-mxn}
\end{algorithm}

\newpage

\begin{algorithm}[h]
\begin{tabbing}
For \= For \= For \= each $A_{k}$ below $A_{i}$ ($i < k \leq m$) xx \= \kill
For each row $A_{i}$ of $A$ ($i=1,\ldots,m$) \\
\> Calculate $d^{'}_{i}$ \> \> \> [Equation \ref{d'i}] \\
\> Calculate $d^{'}_{n}$ \> \> \> [Equation \ref{d'j}] \\
\> Calculate matrix $\Psi _{i,n}$ \> \> \> [Equation \ref{Psiij}] \\
\> Calculate $A^{'}_{i}$, the rotation of $A_{i}$
					 \> \> \> [Equation \ref{APsiA}] \\
\> \\
\> For each $A_{k}$ below $A_{i}$ ($k=i+1,\ldots,m$) \\
\> \> Calculate $a^{'}_{ki}$ \> \> [Equation \ref{a'ki}] \\
\> \> Calculate $a^{'}_{kn}$ \> \> [Equation \ref{a'kj}] \\
\> \> Reset $A_{k}$ to $A^{'}_{k}$ \> \\
\> EndFor \\
\> \\
\> Reset $A_{i}$ to $A^{'}_{i}$ \> \\
\> Reset $D$ to $D^{'}$ \> \\
\> \\
\> For each element $j$ of $A_{i}$ ($j=n-1,\ldots,i+1$) \\
\> \> Calculate $d^{'}_{i}$ \> \> [Equation \ref{d'i}] \\
\> \> Calculate $d^{'}_{j}$ \> \> [Equation \ref{d'j}] \\
\> \> Calculate matrix $\Psi _{i,j}$  \> \> [Equation \ref{Psiij}] \\
\> \> Calculate $A^{'}_{i}$, the rotation of $A_{i}$
					 \> \> [Equation \ref{APsiA}] \\
\> \> \\
\> \> For each $A_{k}$ below $A_{i}$ ($k=i+1,\ldots,m$) \\
\> \> \> Calculate $a^{'}_{kj}$ \> [Equation \ref{a'kj-opt}] \\
\> \> \> Calculate $a^{'}_{ki}$ \> [Equation \ref{a'ki-opt}] \\
\> \> \> Reset $A_{k}$ to $A^{'}_{k}$ \> \\
\> \> EndFor \\
\> \> \\
\> \> Reset $A_{i}$ to $A^{'}_{i}$ \> \\
\> \> Reset $D$ to $D^{'}$ \> \\
\> EndFor \\
EndFor 
\end{tabbing}
\caption{Improved fast Givens rotations \protect\cite{gentleman73} : $m \times n$ matrix}
\label{givens-mxn-opt}
\end{algorithm}

\newpage

\begin{algorithm}[h]
\begin{tabbing}
For \= For \= For \= each $A_{k}$ below $A_{i}$ xx \= \kill
For each row $A_{i}$ of $A$ ($i=1,\ldots,M$) \\
\> Procedure $W_{i,N}$ \> \> \> 
	[Equations \ref{d'i},\ref{d'j},\ref{Psiij},\ref{APsiA}] \\
\> \\
\> For each $A_{k}$ below $A_{i}$ ($k=i+1,\ldots,M$) \\
\> \> Procedure $X_{i,N,k}$ \> \>
	[Equations \ref{a'ki},\ref{a'kj}] \\
\> EndFor \\
\> \\
\> For each element $j$ of $A_{i}$ ($j=N-1,\ldots,i+1$) \\
\> \> Procedure $Y_{i,j}$ \> \>
	[Equations \ref{d'i},\ref{d'j},\ref{Psiij},\ref{APsiA}] \\
\> \> \\
\> \> For each $A_{k}$ below $A_{i}$ ($k=i+1,\ldots,M$) \\
\> \> \> Procedure $Z_{i,j,k}$ \>
	[Equations \ref{a'kj-opt},\ref{a'ki-opt}] \\
\> \> EndFor \\
\> \> \\
\> EndFor \\
EndFor 
\end{tabbing}
\caption{Condensed Description of Fast Givens Rotations Algorithm}
\label{alggivens-short}
\end{algorithm}

\newpage

\section*{Biography}

David R. O'Hallaron received the Ph.D. degree in Computer Science from
the University of Virginia, Charlottesville, in 1986.

From 1986 to 1989 he was with the GE Research and Development Center,
Schenectady, NY, as a Computer Scientist. In 1989, he joined the
faculty in the School of Computer Science, Carnegie Mellon University,
Pittsburgh, PA, as a Systems Scientist. His main research interests
are in the areas of programming models, algorithms, and applications
for parallel machines.

\end{document}


