\documentstyle[twocolumn,epsf]{article} \pagestyle{empty} \setlength{\textwidth}{6.9in} \setlength{\textheight}{8.75in} \setlength{\columnsep}{2.5pc} \setlength{\topmargin}{-0.3in} \setlength{\oddsidemargin}{-.3in} \setlength{\parindent}{1pc} \makeatletter \def\@normalsize{\@setsize\normalsize{12pt}\xpt\@xpt \abovedisplayskip 11pt plus2pt minus5pt\belowdisplayskip \abovedisplayskip \abovedisplayshortskip \z@ plus3pt\belowdisplayshortskip 6pt plus3pt minus3pt\let\@listi\@listI} %the following line was changed \def\subsize{\@setsize\subsize{12pt}\xipt\@xipt} \def\section{\@startsection{section}{1}{\z@}{24pt plus 2 pt minus 2 pt} {12pt plus 2pt minus 2pt}{\large\bf}} \def\subsection{\@startsection {subsection}{2}{\z@}{12pt plus 2pt minus 2pt}{12pt plus 2pt minus 2pt}{\subsize\bf}} \makeatother \input{psfig} \begin{document} \bibliographystyle{plain} \date{} \title{\LARGE \bf Parallel Inversion of Polynomial Matrices} \author{ \begin{tabular}[t]{ccc} Alina Solovyova-Vincent & Frederick C. Harris, Jr. & M. Sami Fadali \\ Computer Science & Computer Science & Electrical Engineering \\ University of Nevada, Reno & University of Nevada, Reno & University of Nevada, Reno\\ Reno, NV 89557 & Reno, NV 89557 & Reno, NV 89557 \\ alina@cs.unr.edu & fredh@cs.unr.edu & fadali@ee.unr.edu \end{tabular} } \maketitle \thispagestyle{empty} \subsection*{\centering Abstract} \vspace*{-3mm} %\hspace{1em} This paper presents an overview of different methods proposed in the last several decades for computing the inverse of a polynomial matrix, concentrating on Bus\l owicz's algorithm. A detailed description of Bus\l owicz's algorithm and its sequential implementation is followed by the presentation of a new parallel algorithm, based on Bus\l owicz's. The distributed and shared memory versions of this parallel algorithm are discussed, and the resulting computation times are analyzed and compared. \noindent {\bf keywords:} parallel algorithm, polynomial matrix inversion \section{Introduction} The problem of inverting polynomial matrices (or, more generally, rational matrices) has been under investigation for over half of a century. This research is well motivated because the computation of such inverses is needed in many fields. For instance, in multivariable control systems, a system is often described by a matrix of rational transfer functions. The problem of finding the inverse of a rational matrix arises in analysis and design using the inverse Nyquist array method~\cite{Mac,Rosenbrock-74}, in parameterization design of linear decoupling controllers~\cite{Lin2,Pavel}, in robust stability analysis~\cite{Fadali}, and in design using the QFT method~\cite{Horowitz,Mac}. The inversion of polynomial matrices is also required in various fields of control system synthesis~\cite{Kail,Wol}. Furthermore, the inversion of rational matrices is required in the analysis and synthesis of passive and active RLC networks for inversion of admittance or impedance matrices~\cite{Keon} and in the analysis of power systems using the method of diakoptics~\cite{Bram}. When a rational matrix is expressed as a ratio of a numerator polynomial matrix and a denominator scalar polynomial, the computation of the inverse essentially reduces to the computation of the inverse of a polynomial matrix~\cite{Lin}. Thus, in many cases, the problem of finding the inverse of a rational matrix can be solved by inverting the corresponding polynomial matrix. The rest of this paper is outlined as follows: Section 2 introduces definitions and notations and provides the overview of other existing inversion methods along with their advantages and disadvantages. It also introduces Bus\l owicz's algorithm and outlines reasons for selecting this algorithm as the basis for a parallel implementation. Section 3 describes the details of the sequential implementation of the algorithm as well as the changes necessary to parallelize it. Discussions of the shared memory and distributed memory parallel implementations complete this section. Section 4 presents and analyzes the results of the sequential and parallel versions of the program. Conclusions and directions for future work are provided in Section 5. \section{Notation, Literature Review, and Bus\l owicz's Algorithm} % \hspace{1em} % This section presents the problem of inverting a polynomial % matrix. Section 2.1 introduces the notation and defines some of % the terms commonly used in this area of research. Section 2.2 % provides an overview of different methods and algorithms % appearing in the literature for inversion of rational and % polynomial matrices. Section 2.3 focuses on Bus\l owicz's % algorithm~\cite{Bus} and gives its detailed description. % We concentrate on this algorithm because it is the basis for % our parallel implementation discussed in Section 3.2. \subsection{Introduction of Notation} A polynomial matrix is a matrix which has polynomials in all of its entries. Consider a polynomial matrix $H(s)$ of degree $n$ \[ H(s)=H_{n}s^{n}+H_{n-1}s^{n-1}+H_{n-2}s^{n-2}+...+H_{0}, \] where $H_{i}$ are $r \times r$ constant square matrices, $i=0,...,n$. An example of such a matrix is % \[ H(s)=\left[ \begin{array}{ll} s+2 & s^{3}+3s^{2}+s \\ s^{3} & s^{2}+1 \end{array} \right] . \] In this case, the degree of the polynomial matrix is $n=3,$ and the size of the matrix $H_{i}$ is $r$ $=$ $2$. For this example, \[ H_{0}=\left[ \begin{array}{ll} 2 & 0 \\ 0 & 1 \end{array} \right] , H_{1}=\left[ \begin{array}{ll} 1 & 1 \\ 0 & 0 \end{array} , \right] \] \[ H_{2}=\left[ \begin{array}{ll} 0 & 3 \\ 0 & 1 \end{array} \right] , H_{3}=\left[ \begin{array}{ll} 0 & 1 \\ 1 & 0 \end{array} \right] . \] % % % \[ % H(s)=\left[ % \begin{array}{ll} % s^{3}+2s+1 & s^{3}+s^{2} \\ % s & s^{2}+s % \end{array} % \right] . % \] % In this case, the degree of the polynomial matrix is $n=3,$ and % the size of the matrix $H_{i}$ is $r$ $=$ $2$. % \newpage % For this example, % % \[ % H_{0}=\left[ % \begin{array}{ll} % 1 & 0 \\ % 0 & 0 % \end{array} % \right] ,H_{1}=\left[ % \begin{array}{ll} % 2 & 0 \\ % 1 & 1 % \end{array} % \right] ,H_{2}=\left[ % \begin{array}{ll} % 0 & 1 \\ % 0 & 1 % \end{array} % \right] ,H_{3}=\left[ % \begin{array}{ll} % 1 & 1 \\ % 0 & 0 % \end{array} % \right] . % \] $H(s)$ is considered column proper if its highest degree coefficient matrix $H_{n}$ is non-singular~\cite{Wol}. $H(s)$ is row proper if its transpose, $H^{T}(s)$, is column proper. The notation used to denote the inverse of a matrix is $H^{-1}(s).$ Only unimodular matrices ($i.e.$, polynomial matrices with a non-zero determinant that is independent of $s$) have inverses that are themselves polynomial matrices~\cite{Mun}. Rational matrices are matrices whose entries are rational functions in $s$, which are non-singular at $s=0$. A rational function can be expressed as $Hd^{-1}$, where $H$ is a polynomial matrix and $d$ is a scalar polynomial. Thus the problem of inverting a rational matrix can be reduced to inverting a polynomial matrix. \subsection{Literature Review} \label{sec:litreview} % \hspace{1em} Over the last several decades there have been % several methods proposed for the inversion of polynomial or % rational matrices. This subsection provides an overview of those % methods and considers some of the advantages and disadvantages % of each approach. We begin with a review of the special case of inverting the resolvent matrix $\left[ sI_{r}-H\right]$, where $I$ is the unit matrix and $H$ is an $r\times r$ matrix of constants. A process for finding $\left[ sI_{r}-H\right]^{-1}$ is well documented and is known as Leverrier's algorithm~\cite{LEVERRIER}. Leverrier's algorithm as well as multiple extensions of this method ($i.e.$,~Leverrier-Faddeev algorithm~\cite{FADDEEV}, Souriau-Frame-Faddeev algorithm~\cite{ROSENBROCK}, {\it etc.}) serve as a basis for several matrix inversion techniques that follow. A number of different approaches for the inversion of polynomial matrices have been proposed over the past years. The assumptions made by different authors vary, and the results do not always have the same form. One of the first papers on this topic by Kosut~\cite{Kos} gives a direct algorithm based on a generalization of the Leverrier method. His method contains many polynomial operations and is not very general. Munro and Zakian~\cite{Mun} used the approach suggested by Kosut for the inversion of rational polynomial matrices by the Souriau-Frame-Faddeev algorithm. They considered two distinct methods in their paper: one based on the Gaussian elimination algorithm and the other one based on the Faddeev algorithm. Both methods involve performing direct computation of the adjoint matrix obtained by polynomial operations. However, their methods have downfalls in that such operations are lengthy, require a ``large degree of involved bookkeeping''~\cite{Mun}, and are known to cause numerical problems. In addition, operations in the field of rational functions utilized in both methods are not suitable for computer programming. Downs~\cite{D3,D4} presented another approach, based on exact Gaussian elimination for matrices with integer coefficients. His method still contained many polynomial operations. Almost at the same time, Emre $et $ $al.$~\cite{Emre1} proposed a method of inversion of rational matrices based on Cramer's rule. The primary motivation for introducing this new method was to avoid polynomial arithmetic and to establish an algorithm systematically dealing with constant matrices. This approach required only simple arithmetic. Their method originally required restrictive assumptions that the polynomial matrix $H(s)$ is non-singular at $s=0$ and that the determinant is known at the outset. The problem of polynomial cancellation was not taken into account by Emre $et$ $al$. Downs was the one to point out the many restrictions and problems of their approach. In a series of papers that followed~\cite{D5,Emre3,Emre2}, most of these problems were resolved. Another point worth mentioning is that the inversion presented by Emre $et$ $al.$ was carried out by computing the determinants recursively. This inversion method requires that the determinants of $(n+1)r$ constant matrices be evaluated in order to compute the determinant of a polynomial matrix $H(s)$ of order $r$ whose degree is $n$. Computation time for this method is large for large $r$ and $n$. Inouye~\cite{Ino} approached the problem of inverting polynomial matrices by generalizing Fadeev's recursive formula. His method is an extension of the Souriau-Frame-Faddeev algorithm. It does not require prerequisite determinants and requires operations with only constant matrices. It simultaneously determines the adjoint matrices and the coefficients of the determinants. The author showed that his algorithm is ``faster than the existing ones.'' One of the disadvantages of his method is that it works only for row- or column-proper polynomial matrices. It also gives the inverse in the minimal degree form only if the polynomial matrix to be inverted is not a special form, but it cannot ensure that the denominator and inversion numerator matrix obtained will be irreducible for a general case. Num~\cite{Num}, and much later Schuter and Hippe~\cite{Sch}, proposed finding the inverse by generalizing known polynomial interpolation approaches. Both techniques require a careful choice of base points in order to avoid ill-conditioned equations. Both methods require complex computations. Another problem with interpolation methods is that only upper bounds for the degrees of the determinant and the adjoint are usually available. The interpolation thus involves redundant equations and polynomials with unnecessarily high degrees. In 1980 Bus\l owicz~\cite{Bus} published a paper with a method that is similar to the method proposed by Inouye~\cite{Ino} but more general in that it works for any non-singular polynomial matrix (as opposed to only row- or column-proper matrices). Bus\l owicz's recursive algorithm computes the inverse by Cramer's rule, explicitly calculating the adjoint matrix and the determinant starting from the coefficient matrices. It requires operations with only constant matrices. The drawback of Bus\l owicz's algorithm is that the irreducible form cannot be ensured in general. Still another approach was developed independently by Zhang~\cite{Zhang} and Chang $et$ $al$.~\cite{Chang}. They both used a division algorithm for polynomial matrices to compute the inverse in irreducible form; however, their algorithms had increased computational complexity. Fragulis $et$ $al$.~\cite{Fra} proposed an algorithm that is a generalization of the Leverrier-type algorithm. The inverse is calculated using the recursive formula. Their approach does not seem to be significantly different from the one proposed by Bus\l owicz and does not provide any clear advantage over it. The method of finding the inverse of a polynomial matrix based on state space realizations is used by Lin and Hsieh~\cite{Lin}. They compute neither the determinant nor the adjoint matrix. This method does not yield the exact solution, but experiments show that the algorithm gives accurate results for rational matrices that arise in the analysis and design of linear multivariable control systems. \subsection{Bus\l owicz's Algorithm} % \hspace{1em}This subsection explains the reasons for choosing % the algorithm developed by Bus\l owicz~\cite{Bus} for finding the % inverse of a polynomial matrix. It also describes the algorithm % itself. \subsubsection{Advantages of Bus\l owicz's Algorithm} \hspace{1pc} As mentioned in Section~\ref{sec:litreview}, Bus\l owicz based his approach for finding the inverse of a polynomial matrix on a generalization of Fadeev's recursive formula. A similar method was proposed by Inouye~\cite{Ino} in 1979, and it was the fastest and most general method at that time. %One major drawback of the latter algorithm was that fact %that it was suitable only for row- or column-proper polynomial %matrices. Bus\l owicz's algorithm is even more general, does not require knowledge of the determinant at the outset, and works for any non-singular polynomial matrix. The only operations required are those on constant matrices. There were several reasons for choosing this method for implementation. First, methods proposed before Bus\l owicz published his paper were obviously less general. Second, we wanted to implement an exact method, thus eliminating available algorithms that use approximations or interpolations such as \cite{Num}, \cite{Sch} and \cite{Lin}. Two other newer methods \cite{Chang,Zhang}, provide only slight improvement of the results in that they yield the inverse in the already irreducible form. However, both authors agree that their algorithms require additional ``complex computations''. Finally, Bus\l owicz claimed that his algorithm was suitable for computer programming \cite{Bus}. We agreed with this assessment and also saw a potential for great speedup in the parallel implementation. \subsubsection{The Algorithm} \hspace{1pc} One of the general ways to compute the inverse of a matrix $H(s)$ is to evaluate the expression given by \begin{equation} H^{-1}(s)=\frac{{\rm adj}\, H(s)}{\det H(s)}, \label{H^-1} \end{equation} where ${\rm adj}\,H(s)$ denotes the adjacent matrix $H(s)$, which is found as \begin{equation} {\rm adj}\,H(s)=\sum\limits_{k=0}^{n(r-1)}Q_{k}s^{k},\ Q_{k}\in R^{r\times r} \label{adj} \end{equation} and \begin{equation} \det H(s)=\sum\limits_{k=0}^{nr}a_{k}s^{k},\ a\in R. \label{detH} \end{equation} The problem of finding the inverse of a polynomial matrix comes down to finding an efficient method for calculating matrices $Q_{k}$, $ k=0,1,...,n(r-1)$, and the coefficients $a_{k},$ $k=0,1,...,rn$, from the given matrices $H_{i}$, $i=0,1,...,n.$ Bus\l owicz showed in his paper that the matrices $Q_{k}$ of ${\rm adj}\,H(s)$ can be computed as \begin{equation} Q_{k}=(-1)^{r+1}R_{r-1,k}, \ k=0,1,...,n(r-1), \label{Q} \end{equation} and the coefficients $a_{k}$ of $\det H(s)$ can be found using the formula \begin{equation} a_{k}=\frac{(-1)^{r+1}}{r} {\rm tr}G_{r,k}, \ k=0,1,...,nr \label{a} \end{equation} where tr denotes the trace of a matrix. The matrices $R_{r-1,k}$ and $G_{r,k}$ appearing in the above expressions are computed from the following iterative formulae: {\small \begin{eqnarray} G_{i,k} &=&H_{0}R_{i-1,k}+H_{1}R_{i-1,k-1}+...+H_{n}R_{i-1,k-n}, \label{G} \end{eqnarray} } \begin{equation} a_{i,k}=-\frac{1}{i} {\rm tr}G_{i,k},\ i=1,2,...,r, \ {\rm and} \label{a(i,k)} \end{equation} \begin{equation} R_{i,k}=G_{i,k}+I_{r}a_{i,k}, \ i=1,2,...,r-1 \ {\rm and }\ k=0,1,...,in, \label{R} \end{equation} \noindent where \begin{equation} R_{0,k}=\left\{ \begin{array}{ll} I_{r} & {\rm for} \ k=0 \\ R_{0,k}= 0 & {\rm for } \ k\neq 0 \end{array} \right. \end{equation} \noindent and \begin{equation} R_{i,k}=0 \ {\rm for} \ j<0 \ {\rm or} \ k<0 \ {\rm or} \ k>jn. \label{Rconditions} \end{equation} \noindent In addition, \begin{equation} R_{r,k}=G_{r,k}+I_{r}a_{r,k}=0, \ k=0,1,...,rn. \end{equation} The algorithm proposed by Bus\l owicz for inversion of the polynomial matrices then consists of the following steps: \begin{enumerate} \item Using formulae (\ref{G})-(\ref{Rconditions}), calculate $G_{i,k},a_{i,k}$ and $R_{i,k}$ for $i=1,2,...,r-1$ and $k=0,1,...,in$, and calculate from formula~(\ref{Q}) the matrices $Q_{k} \ for \ k=0,1,...,n(r-1)$. \item Using formulae (\ref{Q}) and (\ref{adj}), calculate the matrix $adjH(s).$ \item Calculate the matrices $G_{r,k}$ for $k=0,1,...,rn$ from the formulae (\ref {G})-(\ref{R}) and the coefficients $a_{i,k}$ of the polynomial $\det H(s)$ from formula (\ref{a}). \item From formula (\ref{detH}) calculate the polynomial $\det H(s).$ \item Calculate the matrix $H^{-1}(s)$ from formula (\ref{H^-1}). \item The computations could be checked using the following equation: \[ G_{r,k}+(-1)^{r}a_{k}I_{r}=0, \ k=0,1,...,rn. \] Note: In the case where the polynomial matrix has no inverse, the coefficients $a_{k}$, $k=0,1,...,rn$, calculated by formula (\ref{a}) will be equal to zero. \end{enumerate} \section{Sequential and Parallel Algorithms} \subsection{Sequential Algorithm} %Figures \ref{fig:sequential} and \ref{fig:check} show %important parts of our implementation of Bus\l owicz's %sequential algorithm. Before the sequential code is discussed in detail, we introduce the variables that appear in the program. \vspace{-6pt} %\begin{itemize} \begin{list}{$\bullet$}{\parsep=0pt \itemsep=0pt} \item Each $H_{i}$ is an $r\times r$ matrix. $n+1$ of them form the matrix $H(s)$, the inverse of which is to be computed (see the first formula in 2.1). Each $H_{i}$ is represented by a 3-dimensional array $H[i][x][y]$, where $i=0,1,...,n$, $x=0,1,...,r-1$ and $ y=0,1,...,r-1$. \item Each $G_{i,k}$ is an $r\times r$ matrix (see formula (6) in 2.3.2), that is represented by a 4-dimensional array $G[i][k][x][y]$, where $i=0,2,...,r$, $k=0,1,...,rn$, $ x=0,1,...,r-1$ and $y=0,1,...,r-1.$ \item Each $R_{i,k}$ is an $r\times r$ matrix (see formulae (8)-(11)), that is represented by a 4-dimensional array $R[i][k][x][y]$, where $i=0,2,...,r$, $k=0,1,...,rn$, $ x=0,1,...,r-1$ and $y=0,1,...,r-1.$ %i=1,..r-1; k=0,..in \item Each $a_{i,k}$ is a coefficient (see formula (7)), that is represented by a 2-dimensional array $a[i][k]$, where $i=0,1,...,r$ and $k=0,1,...,rn$. \item Each $a_{k}$ is a coefficient of det$H(s)$ (see formula (3)), that is represented by a 1-dimensional array $alpha[k]$, where $k=0,1,...,rn$. \item Each $Q_{k}$ is one of $r\times r$ matrices that compose adj$H(s)$ (see formula (\ref{adj})). They are represented by a 3-dimensional array $Q[k][x][y]$, where $i=0,1,...,n$, $x=0,1,...,r-1$ and $ y=0,1,...,r-1$. \item $Ident[x][y]$ is an $r\times r$ unit matrix. \end{list} %\end{itemize} \vspace{-6pt} % \begin{figure} % \input{s3-sequential} % \caption{\label{fig:sequential} The sequential algorithm.} % \end{figure} % Figure \ref{fig:check} presents the final steps of the algorithm % along with the condition that checks for the correctness of the % computations (formula (11)). At this point, checking % the value of $a_{k}$ can be used to verify if a given polynomial % matrix is nonsingular (nonzero determinant). The % coefficients $a_{k}$ for $k=0,1,....,rn$ will be equal to zero % if the matrix is singular. Several changes had to be made to the algorithm outlined in Bus\l owicz's paper. First of all, step 1 of the algorithm (see Section 2.3.2) specifies the calculation of $G_{i,k},$ $a_{i,k}$ and $R_{i,k}$ for $i=1,2,...,r-1$ and $k=0,1,...,in$. However, our algorithm separates this step into two steps because different approaches are required for calculating the variables for $i=1$ and $i>1$. Thus, the algorithm first calculates $G_{i,k},$ $a_{i,k}$ and $R_{i,k}$ for $i=1$ and $k=0,1,...in$ and then continues with the rest of the calculations for $i>1$. Second, as mentioned above, steps 1 and 3 of the algorithm are combined. Computations of $Q_{k}$ are delayed until everything else in steps 1 and 3 is calculated. Hence, step 2 is also performed later in the program. % (see Figure~\ref{fig:check}). % \begin{figure}[h] % \input{s3-check} % \caption{\label{fig:check} Calculate and check the % coefficients of the inverse matrix.} % \end{figure} \subsection{Parallel Algorithm} Armed with a working sequential version of Bus\l owicz's algorithm, we began analyzing program dependencies in order to decide on parallelization techniques. This section presents the details of the parallel implementation of Bus\l owicz's algorithm and outlines the changes made and challenges encountered in the process of parallelizing the sequential version of the program. \begin{figure} \input{s3-parallel} \caption{\label{fig:parallel} The parallel algorithm.} \end{figure} The parallel algorithm is given in Figure~\ref{fig:parallel}. Two new variables appear in this parallel code segment: {\it NUMPROC} is the number of processors used for calculations, and {\it p} is the distinct number associated with each processor {\it p=}0{\it ,...,NUMPROC--}1. Because a SPMD (single program multiple data) programming structure was used, each processor executed the code shown on its portion of data. This algorithm was implemented for both distributed memory and shared memory machines. Implementing the program in a shared memory environment allowed the creation of variables that could be accessed directly by every process. In the shared memory environment, the shared memory segments are created using the {\it shmget()} system calls. Because there is a limit on the number of shared memory segments that can be created, 2-, 3- and 4-dimensional matrices are represented as 1-dimensional arrays. Shared memory segments are attached to the data segments of the calling process before performing calculations using {\it shmat()} and then are detached after the computations are completed. In a distributed memory environment, variables computed by one process that are required by another have to be passed explicitly by the program. MPI was chosen to provide the functionality required for programming in a distributed memory environment. Because most parallelism occurs in the loops, a first attempt to parallelize any code typically requires looking for independent loops that can be split across multiple processors. Independent loops can be executed in any order without affecting the semantics of the program. There are several $for$ loops in the sequential program, but, unfortunately, not all are independent. Clearly, the large outer $i$-loop (line 15 in Figure~\ref{fig:parallel}) is not independent. Calculations in the $i^{th}$ iteration depend on the results of the previous $(i-1)^{st}$ iteration because the $i^{th}$ iteration involves operations on $R_{i-1,k}$ (line 31). However, the $k$-loops are independent and can be parallelized (lines 1 and 17 in Figure~\ref{fig:parallel}). This parallelization was accomplished by performing striped partitioning of the matrices across the processors. The presence of the dependent loops in the program created another challenge: synchronization of the processes and data. Looking at lines 15 and 17, one can notice that there are two nested loops, with an independent loop inside the dependent one. To make matters worse, the number of inner iterations ($ll$-loop on line 22 varies. The dependence on $k$ can be seen in lines 19-22 Thus there are so-called partially parallel loops, $i.e.$, loops whose parallelization requires synchronization to ensure that iterations are executed in the correct order and produce the correct output. Specifically, no process can go on with execution of the $i^{th}$ iteration until every other process had completed its $(i-1)^{st}$ iteration. In a shared memory environment, this synchronization is accomplished by placing a barrier before starting the next iteration of the $i$-loop (line 42 in Figure~\ref{fig:parallel}). Another barrier is placed on line 14 (Figure~\ref{fig:parallel}) to synchronize the processes, making sure that $R_{1,k}$ is calculated for all values of $k$ before continuing with calculations for $i>1.$ % In a distributed memory environment, synchronization is as important, but the data calculated by the processes must also be explicitly exchanged so it can be used in the next iteration by other processors. This explicit exchange using MPI communication calls implicitly accomplishes the process synchronization required. The rest of the program is left unchanged from the sequential version. \section{Results} The distributed memory code using MPI was tested on three different platforms. The first was a network of SGI O2 workstations. These machines have 180MHz MIPS R5000 processors with 320MB ram. The second platform was a network of Pentium IV workstations, each with a 1.8 GHz processor and 256MB ram. Both of these first two platforms have a standard 100 megabit network. The third platform is a cluster of Pentium IV Xeon processors, each at 2.2 GHz with 2GB of ram. The communication network is Myrinet 2000. The shared memory code was tested on two different platforms. The first platform was an SGI Power Challenge 10000. This machine is a shared memory multiprocessor, consisting of 8 MIPS R10000 and 1 GB of ram. The second platform was an SGI Origin 2000. This machine is a shared memory multiprocessor, consisting of 16 MIPS R12000 300MHz processors and 2 GB of ram. The complexity of the implemented sequential algorithm is ${\cal O}(n^{2}r^5)$. Thus the run times increase rapidly as the problem size increases. The problem size can be increased either by scaling the degree of the polynomial matrix $n$, the size of the matrix $r$, or both. We considered only real-life cases in the field of control theory, where neither the size of the matrix nor the degree of the polynomial typically exceeds $25$. The figures in this section illustrate the computation times of a sequential program under various conditions as well as computation times obtained on the distributed and shared memory platforms with various numbers of processors. For comparison of the platforms, the sequential run times for the largest problem size are provided in Table~\ref{tab:table-seq}. \begin{table} \begin{center} \begin{tabular}{|c|crc|}\hline Platform & \multicolumn{3}{c|}{Sequential Time (sec)} \\ \hline SGI O2 NOW & \hspace{0.3in} &2645.30& \\ \hline P IV NOW & &29.99& \\ \hline P IV Cluster & &18.75& \\ \hline SGI Power Challenge & &913.99& \\ \hline SGI Origin 2000 & &552.95& \\ \hline \end{tabular} \end{center} \caption{Sequential run times for different platforms ($n=25$, $r=25$).} \label{tab:table-seq} \end{table} \subsection{Distributed Memory Implementation} The results obtained on the distributed memory platforms were not as expected. On the first two platforms, the algorithm provided some speedup on two processors for all problem sizes. However, when more processors were added, speedup was obtained only on the larger problem sizes, and the efficiency decreased drastically. On the third platform we obtained speedup across all processors, but the efficiency was poor for more than two processors. The efficiency of the parallel algorithm on the third distributed memory platform is shown in Table~\ref{tab:cortex-eff}. \begin{table}[h] \begin{center} \begin{tabular}{|c|c|c|c|} \hline Processors & 2 & 4 & 8 \\ \hline Efficiency & 89.6\% & 68.8\% & 49.7\% \\ \hline \end{tabular} \end{center} \caption{Efficiency, P IV cluster ($n=25$, $r=25$).} \label{tab:cortex-eff} \end{table} %\newpage \subsection{Shared Memory Implementation} The results obtained on both shared memory platforms were outstanding. Figure~\ref{fig:fig1} represents the average computation times (in seconds) on the first shared memory platform for the case when the degree of the polynomial matrix was fixed ($n=10$) and the matrix size was varied from $r=2$ to $r=25$. As we saw earlier, the algorithm is ${\cal O}(n^2r^5)$. Figure~\ref{fig:8procmesh} provides a computation-time surface showing how the changes of problem size ($n$ and $r$) on 8 processors cause the computation time to increase drastically. \begin{figure}[h] \epsfxsize=4in \centerline{\epsffile{figs/nconst.ps}} \caption{Run times, SGI Power Challenge ($n=10$, $r$ is varied).} \label{fig:fig1} \end{figure} \begin{figure} \epsfxsize=4in \centerline{\epsffile{figs/8procmesh.ps}} \caption{Run times, SGI Power Challenge, 8 processors, varied $n$ and $r$.} \label{fig:8procmesh} \end{figure} Table~\ref{tab:table1} presents the average efficiency of our algorithm for $n=25$ and $r=25$ on both shared memory platforms. \begin{table}[h] \begin{center} \begin{tabular}{|c|c|c|c|c|c|c|c|c|}\hline Processors & 2 & 3 & 4 \\ \hline SGI Power Challenge & 99.65\% & 98.4\% & 98.2\\ \hline SGI Origin 2000 & 99.9\% & 101.0\% & 98.7\% \\ \hline \end{tabular} \begin{tabular}{|c|c|c|c|c|c|c|c|c|}\hline 5 & 6 & 7 & 8 & 16\\ \hline 98.2\% & 97.9\% & 97.9\% & 95.8\% & n/a \\ \hline 100.5\% & 99.0\% & 98.7\% & 98.2\% & 93.8\% \\ \hline \end{tabular} \end{center} \caption{Efficiency, shared memory platforms ($n=25$, $r=25$).} \label{tab:table1} \end{table} \section{Conclusions and Future Work} % \hspace{1em} % The results obtained reflect the fact that this problem is highly % data intensive. Therefore, the need to transfer large amounts % of data between processors after each iteration overemphasized % the weaknesses of a distributed memory environment, making the % speedup minimal. At the same time the fast shared memory % implementation on the SGI machines provided excellent performance % since the hardware handled the communication and exchange of data. The results obtained reflect the major difference between shared and distributed memory environments. Excellent performance in the shared memory environment shows that an efficient parallel algorithm can be designed for the highly data intensive problem of polynomial matrix inversion. These excellent results are due to the fact that communication and exchange of data were handled in the hardware in the fast shared memory implementation on the SGIs. In the case of a distributed memory environment, however, the dependencies of the original algorithm required transfer of large amounts of data between processors after each iteration, thus emphasizing the weaknesses of that environment and leading to a minimal speedup. In the distributed memory environment, three platforms were used. %The first had slow CPUs and a slow network interface. The second %had fast CPUs and the same slow network interface. The third %platform had the fastest CPUs and a much faster network interface. Even with faster CPUs and faster networks, the efficiency did not improve much. These poor results were due to the fact that the communication costs far outweighed the performance gain of multiple processors. Because the problem sizes were limited to real-life applications, data sets of much larger sizes were not considered. %However, if they were, performance would be expected %to increase. So when the need for such problems arises, this %distributed memory implementation may be more appropriate. In the shared memory environment, near linear speedup was achieved on both platforms as can be seen in Table~\ref{tab:table1}. This speedup means that the algorithm can take full advantage of the distributed computing power in the shared memory environment as the size of the problem increases. This great speedup can be attributed to the high degree of parallelism we were able to extract from the original algorithm as well as to the elimination of the need for the program to perform the communication explicitly. The efficiency over 100\% on the Origin 2000 can be attributed to the architecture design where the processors are 2 CPUs to a card and the cache is 8 MB per CPU. This design allows the cache for each CPU to be used by either processor on the card, and allows one CPU to use the cache of both when the other CPU is idle. The NUMA memory architecture, which is tightly coupled to the CPU cards, also contributes to the behavior. We have presented a parallel algorithm for computing the inverses of polynomial matrices. We have performed an exhaustive search of all available algorithms for polynomial matrix inversion and based our parallel algorithm on the method proposed in~\cite{Bus}. We have implemented the sequential version as well as two parallel versions. Based on the shared memory implementation results, we conclude that this new parallel algorithm is very efficient but should not be used on a distributed memory environment for small problem sizes. We see this work continuing in a variety of different ways. First, there is an algorithm for inverting multivariable polynomial matrices~\cite{ozguler-80} that has never been parallelized. Second, we anticipate evaluating the distributed memory implementation in order to minimize message passing, thus improving performance. Third, larger problem sizes may also be considered in order to determine when the computation time overtakes the communication overhead. %\baselineskip=12pt \bibliography{paper} \end{document}