\documentclass{mylncse}
\input{psfig.sty}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Peter's macros %%%
\newcommand{\nref}[1]{(\ref{#1})}
\newcommand\dst{\displaystyle}
\newcommand\C{\makebox[1em]{\rm\makebox[.25em][l]{C}\rule{.05em}{1.5ex}}}
\newcommand\R{\hbox{\rm I\kern-2pt R}}
\newcommand\N{\hbox{\rm I\kern-2pt N}}
\newcommand{\Rn}{\R^{n}}
\newcommand{\Rtwonn}{\R^{2n\times n}}
\newcommand{\Rnn}{\R^{n\times n}}
\newcommand{\Rpn}{\R^{p\times n}}
\newcommand{\Rpm}{\R^{p\times m}}
\newcommand{\Rnm}{\R^{n\times m}}
\newcommand{\Rmm}{\R^{m\times m}}
\newcommand{\Rpp}{\R^{p\times p}}
\newcommand{\Rrn}{\R^{r\times n}}
%%%% common double bar norms 
\newcommand {\nrm}	[1] {\Vert #1 \Vert}
\newcommand {\nrmF}	[1] {\nrm{#1}_F}
\newcommand {\nrmone}	[1] {\nrm{#1}_1}
\newcommand {\nrmtwo}	[1] {\nrm{#1}_2}
\newcommand {\nrminf}	[1] {\nrm{#1}_\infty}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% some greek letters
\newcommand{\ve}{\varepsilon}
\newcommand{\la}{\lambda}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% equations
\newcommand{\be}[1]{\begin{equation}\label{#1}}
\newcommand{\ee}{\end{equation}}
\newcommand{\bea}{\begin{eqnarray}}
\newcommand{\eea}{\end{eqnarray}}
\newcommand{\bean}{\begin{eqnarray*}}
\newcommand{\eean}{\end{eqnarray*}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\cond}[1]{{\mathrm cond}\left(#1\right)}
\newcommand{\condF}[1]{{\mathrm cond}_F\left(#1\right)}
\newcommand{\condtwo}[1]{{\mathrm cond}_2\left(#1\right)}
\newcommand{\subsp}[2]{{\mathrm span}\left\{#1,\ldots,#2\right\}}
\newcommand{\diagmat}[1]{{\mathrm diag}\left(#1\right)}
\newcommand{\sign}[1]{{\mathrm sign}\left(#1\right)}
\newcommand{\vecop}[1]{{\mathrm vec}\left(#1\right)}
\newcommand{\trace}[1]{{\mathrm trace}\left(#1\right)}
\newcommand{\eig}[1]{\sigma\,(#1)}
\newcommand{\Matrix}[4]{\left[ \begin{array}{cc} #1 & #2 \\
                         #3 & #4 \end{array} \right] }
\newcommand{\Vector}[2]{\left[ \begin{array}{c} #1 \\ #2 \end{array} \right] }
\newcommand{\sep}{\mathrm{sep}}
\newcommand{\lyap}[1]{A^T #1 E\; +\; E^T #1 A}
\newcommand{\cP}{\mathcal{P}}
\newcommand{\rank}[1]{{\mathrm rank}\left( #1 \right)}
\newcommand{\mlg}{\; \mathrm{ log} \,}
\newcommand{\tA}{\tilde{A}}
\newcommand{\tE}{\tilde{E}}
\newcommand{\tG}{\tilde{G}}
\newcommand{\tQ}{\tilde{Q}}
\newcommand{\tX}{\tilde{X}}
\newcommand{\ti}{{\mathrm iter}}
\newcommand{\tidc}{\mathrm{iter_{dc}}}
\newcommand{\tinw}{\mathrm{iter_{nw}}}
\newcommand{\tism}{\mathrm{iter_{sm}}}
\newcommand{\ric}[1]{{\mathcal{R}_d}\left( #1 \right)}
%
%\newtheorem{theorem}{Theorem}[section]
%\newtheorem{remark}[theorem]{Remark}
%%%%% Marginal Notes %%%%%
% to get a clean copy, comment out the next line and use the 
% next one .
%\newcommand {\mpar} [1] {\marginpar{\fussy \footnotesize{#1}}}
\newcommand {\mpar} [1] {\marginpar{}}
%

\begin{document}
\title{Parallel solvers for \\
       discrete-time algebraic Riccati equations}

\author{{\sc rafael mayo\inst{1}
       \thanks{Correspondence to: Enrique S. Quintana-Ort\'{\i},
               Depto. de Inform\'atica,
               Universidad Jaume I,
               12080--Castell\'on, Spain. \newline
	       Contract/grant sponsor:
               Conseller\'{\i}a de Cultura, Educaci\'on y Ciencia de la 
	       Generalidad Valenciana GV99-59-1-14 and
               the Fundaci\'o Caixa-Castell\'o Bancaixa.}
           \and
        enrique s. quintana-ort\'{\i}\inst{1}
           \and
        gregorio quintana-ort\'{\i}\inst{1}
           \and
        vicente hern\'andez\inst{2}}
}
\institute{Depto. de Inform\'atica,
           Universidad Jaume I,
           12080--Castell\'on, Spain;
           {\tt \{mayo,quintana,gquintan\}@inf.uji.es}.
           Tel.: +34-964-728000. Fax: +34-964-728435.
              \and
           Depto. de Sistemas Inform\'aticos y Computaci\'on,
           Universidad Polit\'ecnica de Valencia,
           46071--Valencia, Spain;
           {\tt vhernand@dsic.upv.es}.
           Tel.: +34-96-3877350. Fax: +34-96-3877359.
}
\maketitle

\vspace*{0.5cm}
\hrule


\vspace*{1cm}
\noindent
{\bf SUMMARY} \\
{\bf
We investigate the numerical solution of discrete-time 
algebraic Riccati equations on a parallel
distributed architecture. Our solvers obtain
an initial solution of the Riccati equation via the disc function method, 
and then refine this solution using Newton's method.
The Smith iteration is employed to solve the Stein equation that
arises at each step of Newton's method.

The numerical experiments on an Intel Pentium-II cluster, connected
via a Myrinet switch, 
report the performance and scalability of the new algorithms.
}

\vspace*{1cm}
\noindent
{\bf 1. INTRODUCTION}
\vspace*{1ex}

\noindent
Consider the discrete-time algebraic Riccati equation (DARE)
\be{dare}
\begin{array}{l@{\quad}c@{\quad}l}
0 = \ric{X} &:=& Q + A^T X A - X - A^T X B (R+B^TXB)^{-1} B^T X A,
\end{array}
\ee
where $A,Q,X \in\Rnn$, $Q=Q^T$, $B \in\Rnm$, $R \in \Rmm$, $R=R^T$, 
and $X=X^T$ is the
sought-after solution.  

Algebraic Riccati equations arise in many
analysis and synthesis algorithms in control theory. 
In particular, in the last years
model reduction for large-scale control systems has become one of the most
important problems in systems and control theory. Many of the 
algorithms recently proposed for this problem
require the solution  
of equations of the form \nref{dare}; see, e.g., \cite{SafC88,Sch96,ZhoDG96}.

In general,
the desired solution of the DARE (hereafter $X_d$) is stabilizing
in the sense that all the eigenvalues
of $A_d = A - B (R+B^TX_dB)^{-1} B^T X_d A$ 
are inside the unit circle (i.e., $A_d$ is d-stable, denoted as $\rho(A_d)<1$).
Throughout this paper we will assume that such a stabilizing solution
exists. Note that if it exists, it is unique \cite{LanR95}. 

The need for parallel computing in this area is clear
as already for control systems with state-space
dimension $n=1000$, \nref{dare} represents a set of nonlinear equations
with more than half a million unknowns (exploiting the symmetry of
$X$). Systems of such a dimension are 
common in chemical engineering applications, are standard for
second order systems, and represent rather coarse grids when derived from the
discretization of a partial differential equation; 
see, e.g., \cite{GarL88b,LauG88,RosW95}. 

The inverse-free iteration for the matrix disc function~\cite{Mal93} 
provides an approximate numerical solution for
the DARE~\cite{Ben97a}. The iterative procedure for the matrix disc function
is composed of QR factorizations and matrix products which can be efficiently
parallelized on a parallel distributed computer.

This approximate solution can then be refined using, e.g.,
Newton's me\-thod~\cite{Hew71,Kle68}, which converges 
to the solution of the DARE~\cite{Hew71}.
At each iteration, the method requires a few basic matrix computations
and the solution of a Stein equation (or discrete-time Lyapunov equation)
of the form
\be{stein}
0 = \tA^T \tX \tA - \tX + \tQ,
\ee
where $\tA,\tQ,\tX \in\Rnn$, 
$\tA$ is d-stable, and $\tQ = \tQ^T$.

Numerical methods for the Stein equation are studied in~\cite{Bar77,Sim96}. 
In the initial stage of these methods, $\tA$ is reduced to real Schur form
by means of the QR algorithm~\cite{GolV89}.
This is followed by a quite less expensive back substitution process. 
The QR algorithm has been reported as difficult to parallelize
on distributed memory architectures due to two factors: 
the algorithm is composed of fine-grain computations with a 
low ratio of computation/communica\-tion and the use of traditional
block scattered data layouts leads to an unbalanced distribution
of the computational load~\cite{HenV97,HenWD97}.
Different approaches have been proposed to avoid these drawbacks
in~\cite{HenV97,HenWD97,WatE91}.
The experimental studies however report parallelism and scalability
results which are still far from those of the kernels
arising in linear systems (see, e.g.,~\cite{Scalapack,DonSS86}).

In this paper we study the numerical solution of
the DARE on a parallel distributed computer. Our Stein solver is based on
the Smith iteration and only
requires matrix products,
which are well-known to be highly efficient 
on these architectures~\cite{Scalapack}.
Thus, we avoid the parallelization difficulties in the QR algorithm.


In Sections~2 and~3 we summarize, respectively,
the disc function and Newton's method for solving the DARE.
In Section~4 we review our Stein solver based on the Smith iteration.
A brief study of the computational and communication cost 
of the parallel implementations
is presented in Section~5. In Section~6 we analyze 
the experimental performance and scalability of our solvers 
on a parallel distributed architecture, an 
Intel Pentium-II cluster, connected via a Myrinet switch. 
Concluding remarks are given in Section~7.

\vspace*{3ex}
\noindent
{\bf 2. THE DISC FUNCTION METHOD FOR THE DARE} 
\vspace*{1ex}

\noindent
In~\cite{Mal93}, Malyshev proposed an inverse-free iteration for
the computation of a deflating subspace of a matrix pencil.
Specifically, the iteration divides the spectrum of the matrix pencil
along the unit circle.
The iteration was made truly inverse free by Bai, Demmel,
and Gu in~\cite{BaiDG97}, and the subspace extraction stage
was improved and the cost further reduced by Sun and Quintana 
in~\cite{SunQ96}. 

The inverse-free iteration for the solving the DARE can be formulated as
follows~\cite{Ben97a}.
We denote the identity matrix of order $n$ as $I_n$.

\vspace*{3ex}
\begin{quote}
\label{alg-disc} 
\underline{\rm The disc function method for the DARE.}  \\[.5ex]
\setlength{\parindent}{0pt}
\hspace*{3ex} \textbf{Input:} 
	       \parbox[t]{10cm}{$A,Q \in \Rnn$,
               $Q=Q^T$, $B \in\Rnm$, $R \in \Rmm$, $R=R^T$, \\ 
               and a user-defined tolerance \textit{tol}.}\\ 
\hspace*{3ex} \textbf{Output:} \parbox[t]{10cm}{
	       Approximate solution $X_{\ast} \in \Rnn$ 
               of \nref{dare}. (Starting guess \\
	       for Newton's method.)} 

\renewcommand{\labelenumii}{\arabic{enumi}.\arabic{enumii}}
%\medskip
\upshape
\begin{enumerate} 
\item Set
\(
    L_0 ~:=~ \left[
        \begin{array}{cc}
            A         & 0
\\          Q     & ~I_n~
        \end{array}
    \right]
\), 
\(
    M_0 ~:=~ \left[
        \begin{array}{cc}
          ~I_n~ & -B R^{-1} B^T
\\         0    & A^T
        \end{array}
    \right]
\), $j ~:=~ 0$.
\item \texttt{WHILE} $\nrmF{R_{j+1} - R_j} ~>~ tol \cdot \nrmF{R_j}$
\begin{enumerate}
\item Compute the QR factorization
\[
    \left[ 
        \begin{array}{c}
          ~ L_j
\\        -M_j
        \end{array}
    \right] = 
    \left[ 
        \begin{array}{cc}
          U_{11} & U_{12}
\\        U_{21} & U_{22}
        \end{array}
    \right]
    \left[ 
        \begin{array}{c}
          R_j
\\        0
        \end{array}
    \right].
\]
\item $L_{j+1} ~:=~ U_{22}^T L_j$.
\item $M_{j+1} ~:=~ U_{12}^T M_j$.
\item $j ~:=~ j+1$.
\end{enumerate}
\texttt{END WHILE}
\item Let $M_j = \left[ M_{1},  M_{2} \right]$,
with $M_{1}$, $M_{2} \in \Rtwonn$,
then  $M_1 X_* = M_2$.
\end{enumerate} 
\end{quote}

The convergence of the previous iteration can be shown to
be globally quadratic~\cite{BaiDG97}. The matrix disc function is related
with the iteration defined by 2.1--2.3 as
\[
{\rm disc}(L_0,M_0) = \lim_{j \rightarrow \infty} (L_j+M_j)^{-1}L_j;
\]
see~\cite{Ben97a}.

\vspace*{3ex}
\noindent
{\bf 3. NEWTON'S METHOD FOR THE DARE}
\vspace*{1ex}

\noindent
Newton's method was first employed by Kleinman, in~\cite{Kle68}, 
to solve the continuous-time version of the algebraic Riccati equation.
Based on this formulation, Hewer extended Newton's method for the
solution of the DARE in~\cite{Hew71}. A discussion of its convergence
properties can be found in~\cite{LanR95,Meh91}.

Specifically, the iteration can be stated as follows.

\vspace*{3ex}
\begin{quote}
\label{alg-newt} 
\underline{\rm Newton's method for the DARE.}  \\[.5ex]
\setlength{\parindent}{0pt}
\hspace*{3ex} \textbf{Input:} 
	       \parbox[t]{10cm}{$A,Q \in \Rnn$,
               $Q=Q^T$, $B \in\Rnm$, $R \in \Rmm$, $R=R^T$, a \\ 
               user-defined tolerance \textit{tol}, and
	       an initial starting guess $X_0 = X^T_0$.} \\
\hspace*{3ex} \textbf{Output:} \parbox[t]{10cm}{
	       Approximate solution $X_{j+1} \in \Rnn$ 
               of \nref{dare}.} 

\renewcommand{\labelenumii}{\arabic{enumi}.\arabic{enumii}}
%\medskip
\upshape
\begin{enumerate} 
\item $j ~:=~ 0$.
\item
\texttt{WHILE} $\nrmF{\ric{X_j}} ~>~ tol\cdot \nrmF{X_j}$
\begin{enumerate}
\item 
$A_j ~:=~ A - B (R+B^TX_jB)^{-1} B^T X_j A$.
\item Solve for $N_j$ the Stein equation
\[
  0 = \ric{X_j} + A_j^T N_j A_j - N_j .
\]
\item 
$X_{j+1} ~:=~ X_j + N_j$.
\item 
$j ~:=~ j+1$.
\end{enumerate}
\texttt{END WHILE}
\end{enumerate} 
\end{quote}

Under mild conditions, all the iterates $X_j$ are stabilizing
(i.e., $\rho(A_j)<1$),
$\lim_{j \rightarrow \infty} X_j = X_d$, and the convergence is
globally quadratic; see~\cite{LanR95}.

\vspace*{3ex}
\noindent
{\bf 4. THE SMITH ITERATION FOR THE STEIN EQUATION}
\vspace*{1ex}

\noindent
We can rewrite equation \nref{stein} in fixed point form, 
$\tX = \tA^T \tX \tA + \tQ$,
and form the fixed point iteration 
\[
X_0 := \tQ, \qquad X_{j+1} := \tQ + \tA^T X_j \tA , \quad j = 0,1,2,\ldots.
\]
Then this iteration converges to $\tX$ if $\rho(\tA) < 1$, 
and the convergence rate is linear. A quadratically 
convergent version of this fixed point iteration is suggested in
\cite{DavM68,Smi68}.
Specifically, the iteration can be written as follows.

\begin{quote}
\label{alg-smith} 
\underline{\rm The squared Smith iteration for the Stein equation.}  \\[.5ex]
\setlength{\parindent}{0pt}
\hspace*{3ex} \textbf{Input:} 
	       \parbox[t]{10cm}{$\tA,\tQ \in \Rnn$,
               $\tQ=\tQ^T$, and 
               a user-defined tolerance \textit{tol}.} \\
\hspace*{3ex} \textbf{Output:} \parbox[t]{10cm}{
	       Approximate solution $X_{j+1} \in \Rnn$ 
               of \nref{stein}.} 

\renewcommand{\labelenumii}{\arabic{enumi}.\arabic{enumii}}
%\medskip
\upshape
\begin{enumerate} 
\item Set $A_0 := \tA$, $X_0 := \tQ$, $j ~:=~ 0$.
\item \texttt{WHILE} $\nrmF{A_j} ~>~ tol$
\begin{enumerate}
\item $X_{j+1} ~:=~  A_j^T X_j A_j + X_j$.
\item $A_{j+1} ~:=~  A_j^2$.
\item $j ~:=~ j+1$.
\end{enumerate}
\texttt{END WHILE}
\end{enumerate} 
\end{quote}

The convergence theory of the Smith iteration derived in \cite{Smi68}
yields that for $\rho(\tA) < 1$ there exist real constants $0 < \mu$ and
$0 < \rho <1$ such that
\be{abserr}
\nrmtwo{\tX - X_j} \leq \mu \nrmtwo{\tQ} (1 - \rho)^{-1}\rho^{2^j}.
\ee
This shows that the method converges
for all equations with $\rho(\tA)<1$, i.e.,
convergence is guaranteed if we employ the Smith iteration for
solving the Stein equations arising at each step of Newton's method. 
Nevertheless, if $\tA$ is highly non-normal such that
$\nrmtwo{\tA} > 1$ 
then overflow may occur in the early stages of the iteration due to increasing
$\nrmtwo{A_j}$ 
although eventually, $\lim_{j\to\infty} A_j = 0$ 
if $\rho(\tA) <1$.

The parallelization of this approach was analyzed in~\cite{BenQQ99c}.

\vspace*{3ex}
\noindent
{\bf 5. PARALLELIZATION AND IMPLEMENTATION DETAILS}
\vspace*{1ex}

\noindent
The most appealing feature of all the numerical methods exposed in
the three previous sections  regarding
its parallel implementation is that all the computational cost comes
from matrix products, linear system solvers, and QR factorizations. 
These kernels are known to be highly parallelizable;
see, e.g., \cite{Scalapack}.

Our parallel algorithms are implemented using the libraries
PB-BLAS (parallel block BLAS) and
ScaLAPACK (scalable linear algebra package)~\cite{Scalapack}.
These are public-domain parallel libraries for MIMD computers which
provide scalable parallel distributed
routines for many matrix algebra kernels.
These libraries employ BLAS 
and LAPACK~\cite{Lapack} for serial computations, and
BLACS (basic linear algebra communication subprograms) for communication.
BLACS can be used on any machine
which supports either PVM~\cite{PVM} or MPI~\cite{MPI},
thus providing a highly portable environment.

For scalability purposes, we employ 
square-like topologies with one process per processor.
The following problem and machine parameters are used:
\begin{list}{-}{}
\item $n$: Size of the problem. (For simplicity, we assume $m=n$.)
\item $p = \sqrt{p} \times \sqrt{p}$: Dimension of the square
       grid of processors. 
\item $\tau$: Time of a double-precision floating-point arithmetic operation.
\item $\alpha$ (latency): Time required to communicate 
      a zero-length message between two processors.
\item $\beta$ (reciprocal of bandwidth): 
      Time required to communicate 
      a double-precision floating-point number between two processors.
\end{list}
For simplicity, we assume the time required to communicate a
message of length $l$ is $\alpha + \beta l$.

We have employed the following parallel (double-precision)
routines from ScaLAPACK~\cite{Scalapack} in our DARE and Stein solvers:
\begin{list}{-}{}
\item {\tt PDGETRF}: LU factorization with partial pivoting.
\item {\tt PDTRSM}:  Solution of a triangular linear system. 
\item {\tt PDGEMM}:  Matrix product.
\item {\tt PDGEQRF}:  QR factorization.
\item {\tt PDORMQR}: Application of Householder reflectors to a matrix.
\end{list}
Following the performance model in~\cite{BaiDDPRS97,BenBQQ98},
we present in Table~\ref{table_1} approximate 
computation and communication costs for these routines
(lower order expression have been neglected).

\begin{table}[hbt]
\begin{center}
\renewcommand{\arraystretch}{1.4}
\begin{tabular}{lccc}\hline  \hline
Parallel      & Computation 
              & \multicolumn{2}{c}{Communication cost} \\
routine       & cost  
	      & Latency
	      & Bandwidth$^{-1}$ \\\hline 
{\tt PDGETRF} & $\frac{2}{3}\times \frac{n^3}{p} \; \tau$ 
	      & $(6+ \mlg p)n \times \; \alpha$ 
	      & $(3+ \frac{1}{4} \mlg p) \times \frac{n^2}{\sqrt{p}}\;\beta$\\ 
{\tt PDTRSM}  & $\frac{n^3}{p} \; \tau$           
	      & $n \times \; \alpha            $ 
	      & $(1+ \frac{3}{4} \mlg p) \times \frac{n^2}{\sqrt{p}}\;\beta$ \\ 
{\tt PDGEMM}  & $2\times \frac{n^3}{p} \; \tau$           
	      & $(1+\frac{1}{2}\mlg p) \sqrt{p} \times \; \alpha $ 
	      & $(1+\frac{1}{2} \mlg p) \times \frac{n^2}{\sqrt{p}}\;\beta$
	      \\
{\tt PDGEQRF} & $\frac{4}{3}\times \frac{n^3}{p} \; \tau$  
	      & $3n \mlg p \times \; \alpha$
	      & $(\frac{3}{4} \mlg p) \times \frac{n^2}{\sqrt{p}}\;\beta$ \\ 
{\tt PDORMQR} & $2 \times \frac{n^3}{p} \; \tau$  
	      & $3n \mlg p \times \; \alpha$                       
	      & $2 \mlg p \times \frac{n^2}{\sqrt{p}}\;\beta$ \\ 
\hline \hline
\end{tabular}
\vspace*{2ex}
\caption{Approximate computation and communication costs for several parallel
routines in PB-BLAS and ScaLAPACK.
\label{table_1}}
\end{center}
\end{table}

A model can be easily constructed from the data in the above table
if we consider the number of calls to each of these basic kernels
in our parallel DARE and Stein solvers:
\begin{list}{--}{}
\item {\tt PDGEDRDC}: the disc function method for the DARE.
\item {\tt PDGEDRNW}: Newton's method for the DARE.
\item {\tt PDGEDLSM}: the Smith iteration for the Stein equation.
\end{list}
Table~\ref{table_2} reports the numbers of calls for
each of the DARE and Stein solvers.
In the table, $\tidc$, $\tinw$, and $\tism$ stand
for the number of iterations of the corresponding iterative method.

\begin{table}[hbt]
\begin{center}
\begin{tabular}{lrrr}\hline  \hline
              & {\tt PDGEDRDC}   & {\tt PDGEDRNW} & {\tt PDGEDLSM}  \\
Kernel        &                  &                &                 \\ \hline
{\tt PDGETRF} &                   & $\tinw$~        &       \\ 
{\tt PDTRSM}  & 1~                & ~$2\times\tinw$~ &       \\ 
{\tt PDGEMM}  & $2\times\tidc$~   & ~$7\times\tinw$~ & ~$3\times\tism$~ \\ 
{\tt PDGEQRF} & ~$1+\tidc$~ &                &       \\ 
{\tt PDORMQR} & ~$1+\tidc$~ &                 &       \\ 
{\tt PDGEDLSM}&                   & $\tinw$~       &       \\ 
\hline\hline
\end{tabular}
\end{center}
\caption{Number of calls to different parallel kernels from 
	 PB-BLAS ScaLAPACK required by the parallel DARE and 
	 Stein solvers.}
\label{table_2}
\end{table}

We observed moderate differences between our theoretical performance
model and the experimental results; these are mostly due to simplifications and
approximations in the performance model (constant cost for any
type of floating-point arithmetic operation -i.e., not taking into account
the memory hierarchy-, linear model for communications, 
ideal implementation of collective communications, etc.);
e.g., for large problems the theoretical execution times are within 13--17\% of
the real execution times for routines {\tt PDGEDLSM} and {\tt PDGEDRNW}.
A slightly larger difference is obtained for routine {\tt PDGEDRDC}.

Efficient convergence criteria for the previous iterative algorithms
have been proposed in~\cite{BenQ99,BenQQ99c,Qui96}.

\vspace*{3ex}
\noindent
{\bf 6. EXPERIMENTAL RESULTS}

\noindent
All our experiments were performed on a Beowulf cluster 
with 32 nodes connected via a {\em Myrinet} switch.
Each node is an Intel Pentium-II processor
at 300MHz with 128MBytes RAM. We employ 
double-precision floating-point {\sc ieee} arithmetic 
($\varepsilon \approx 2.2 \times 10^{-16}$),
and the libraries BLACS (on top of MPI), BLAS, LAPACK, and
ScaLAPACK~\cite{Lapack,Scalapack}.
Routine {\tt DGEMM} achieved on a single node a performance
around 180 Mflops (millions of
floating-point arithmetic operation per second).
A simple ping-pong test reported, for an implementation
of MPI specially tuned for this architecture,
a latency $\alpha \approx$ 18 ns., 
and a bandwidth $\beta^{-1} \approx$ 32 Mb/s. 

Our first experiment analyzes the performance of
the parallel algorithms. Fig.~\ref{fig:time} compares the
execution time of one iteration of serial specially-tuned
implementations of the DARE and Stein solvers 
(using the kernels in the BLAS and LAPACK libraries) 
with those of the corresponding parallel algorithms.
The problem size is fixed at $n$ = 700, and the number of
processors is $p$ = 2,\ldots,12.
When solving the Stein equation in the iteration of Newton's method
(routine {\tt PDGEDRNW}) the number of Smith iterations ($\tism$) was fixed
at 9 in both the serial and the parallel implementations.

\begin{figure}[hbt]
\centerline{\psfig{figure=Plots/extime2.eps,height=7cm,width=10cm}}
\caption{Execution time of one iteration of the serial and parallel
	 algorithms involved in the solution of the DARE.
	 Legend: symbols ``$- \!\!-\! \times \!-\!\!-$'' for {\tt PDGEDRDC},
	                 ``$- \cdot + \cdot -$'' for {\tt PDGEDRNW}, and
	                 ``$\cdots {\rm o} \cdots$'' for {\tt PDGEDLSM}.}
\label{fig:time}
\end{figure}

The results of this experiment also show the speed-up
of the parallel algorithms. Table~\ref{tab:speed-up}
reports a speed-up of 7--9.5 for the largest number of processors.
The higher speed-up of routine {\tt PDGEDRDC} is due to
the larger size of the matrices used in this case
($2n \times n$).

\begin{table}[thb]
\begin{center}
\begin{tabular}{rccc} \hline\hline
    $p$ & ~{\tt PDGEDRDC}~ & ~{\tt PDGEDRNW}~ & ~{\tt PDGEDRSM}~ 
\\ \hline 
      2  &   1.85           &   1.55           &    1.57
\\    4  &   3.45           &   2.84           &    2.91
\\    6  &   5.19           &   3.73           &    4.05
\\    8  &   6.56           &   4.68           &    4.81
\\   10  &   7.94           &   6.09           &    6.34
\\   12  &   9.47           &   7.22           &    7.57
\\ \hline\hline
\end{tabular}
\vspace*{2ex}
\caption{Speed-up of the parallel algorithms involved in
	 the solution of the DARE.}
\label{tab:speed-up}
\end{center}
\end{table}

Our next experiment reports the scalability of the parallel algorithms.
Figure~\ref{fig:scal} shows the Mflops 
per node 
of the serial and the parallel
implementations.
We evaluate the parallel algorithms on $p$=4, 9, 15, 25, and 30 processors,
with $n/\sqrt{p} = 1000$.
The results in the figure show a high scalability of the parallel routines
as the performance remains almost constant as $p$ is increased.

\begin{figure}[hbt]
\centerline{\psfig{figure=Plots/scalab2.eps,height=7cm,width=10cm}}
\caption{Mflops per node for the serial and parallel
	 algorithms involved in the solution of the DARE.
	 Legend: symbols ``$- \!\!-\! \times \!-\!\!-$'' for {\tt PDGEDRDC},
	                 ``$- \cdot + \cdot -$'' for {\tt PDGEDRNW}, and
	                 ``$\cdots {\rm o} \cdots$'' for {\tt PDGEDLSM}.}
\vspace*{-2ex}
\label{fig:scal}
\end{figure}

\vspace*{3ex}
\noindent
{\bf 7. CONCLUDING REMARKS}
\vspace*{1ex}
 
\noindent
We have analyzed the parallel solution of the discrete-time algebraic
Riccati equation by means of the inverse-free iteration for
the matrix disc function and Newton's method.
Our implementation of Newton's method employs the Smith
iteration to solve the corresponding Stein equations.

Our parallel solvers are composed of highly parallel 
medium-grain computational kernels like, e.g.,
QR factorizations, linear system solvers, matrix products, etc. 
The experimental results
on an cluster of Intel Pentium-II, connected with a Myrinet switch,
show the performance and scalability of our parallel solvers.
 
\bibliographystyle{plain}
%\bibliography{control,nla,parallel}
\input{dare.bbl}

\end{document}
