IRIA
laboria

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
78150 - Le Chesnay
Tél. 954.90.20

laboratoire de recherche
en informatique
et automatique

Parallel algorithms
for solving triangular
linear systems
with small parallelism

Laurent
L. Hyafil, H.T. Kung

Rapport de Recherche n°97

Février 1975.
PARALLEL ALGORITHMS FOR SOLVING TRIANGULAR
LINEAR SYSTEMS WITH SMALL PARALLELISM

L. Hyafil

IRIA/LABORIA

H.T. Kung

Department of Computer Sciences

Carnegie Mellon University

Pittsburg, Pa. 15213

Résumé :

Nous étudions dans ce rapport la résolution des systèmes linéaires triangulaires sur un multipro-
cesseur à faible degré de parallélisme. Nous montrons que le problème peut être traité en temps :

$$
\frac{n^2}{k} + O(n) \text{ avec } k \text{ processeurs si } k \le O(n)
$$

$$
O\left(n^{2-r} \log n\right) \text{ avec } O\left(n^r\right) \text{ processeurs si } 1 < r < \frac{3}{2}
$$

$$
\text{et } O\left(n^{1-\frac{r}{3}} \log^2 n\right) \text{ avec } O\left(n^r\right) \text{ processeurs si } \frac{3}{2} \le r \le 3
$$
ABSTRACT

The problem of solving triangular linear systems of size n on a parallel computer with small parallelism is considered. Assume that the time is measured by the number of parallel steps of any arithmetic operation. It is shown that the problem can be done in time

$$ \frac{n^2}{k} + O(n) \text{ with } k \text{ processors for } k \le O(n), \\ O(n^{2-r} \log n) \text{ with } O(n^r) \text{ processors for } 1 < r < \frac{3}{2}, \text{ and} \\ O(n^{1-\frac{r}{3}} \log^2 n) \text{ with } O(n^r) \text{ processors for } \frac{3}{2} \le r \le 3. $$

The results are obtained by using two principles of reducing parallel algorithms with large parallelism to parallel algorithms with small parallelism. The two principles for the reduction are expected to be useful for other problems.

1. INTRODUCTION

We consider the problem of solving triangular linear systems of size n on a SIMD (Single-Instruction Stream-Multiple-Data Stream) array computer such as ILLIAC IV. We assume that the time is measured by the number of parallel steps of any arithmetic operation. The following two results were previously proven:

a. With $n-1$ processors the problem can be done in time $O(n)$ by a direct implementation of Gaussian elimination. (It seems to us that this implementation has not appeared in the literature, so we include it in Section 3 of this paper.)

b. With $O(n^3)$ processors the problem can be done in time $O(\log^2 n)$. (This result was recently shown by a number of people, including Chen and Kuck [74], Heller [74], Munro [74], Orcutt [74], etc., by at least three different ways.)

An obvious drawback of results (a) and (b) is that the number of processors required is prohibitively large when n is large. In this paper, we shall remove the drawback by using two principles of reducing algorithms with large parallelism to algorithms with small parallelism. Our new results are summarized as follows:

a'. By applying the *algorithm-decomposition principle*, result (a) is used to show that with $k$ processors the problem can be done in time $\frac{n^2}{k} + O(n)$ for $k \le O(n)$. If $k$ is a constant or $k = o(n)$, the coefficient of the leading term in the bound is the best possible, since the minimum time sequential algorithm takes time $n^2$.

b'. By applying the *problem-decomposition principle*, result (b) is used to show that with $O(n^r)$
processors the problem can be done in time

$$ (1.1) \qquad O(n^{2-r} \log n) \text{ when } 1 < r < \frac{3}{2} $$

and in time

$$ (1.2) \qquad O(n^{1-\frac{r}{3}} \log^2 n) \text{ when } \frac{3}{2} \le r \le 3. $$

Results on some particular choices of r are given in the following table:

<table><thead><tr><th>Number of Processors</th><th>Upper Bound on Time</th></tr></thead><tbody><tr><td>1</td><td>n<sup>2</sup></td></tr><tr><td>k</td><td>(n<sup>2</sup>/k) + O(n)</td></tr><tr><td>(1/n)<sup>2</sup></td><td>(3/n)<sup>2</sup> + O(n)</td></tr><tr><td>O(n)</td><td>O(n)</td></tr><tr><td>(3/n)<sup>2</sup></td><td>(1/n<sup>2</sup>) + O(n<sup>2</sup>)</td></tr><tr><td>O(n<sup>2</sup>)</td><td>(1/3) + O(n<sup>2</sup>)</td></tr><tr><td>(5/n)<sup>2</sup></td><td>(1/6) + O(n<sup>2</sup>)</td></tr><tr><td>O(n<sup>3</sup>)</td><td>O(log<sup>2</sup>n)</td></tr></tbody></table>

Observe that in the table when the number of processors increases the upper bound on time decreases as one might expect. However, it still remains as an open problem to show whether the bounds in (1.1) and (1.2) are sharp. We believe it is very likely that these upper bounds will be improved later. Furthermore, if the upper bound for the case of O(n<sup>5</sup>) processors is improved, then that for the case of O(n<sup>3</sup>) processors with any r < s, r ≠ 1, can also be improved by the same techniques as used in the paper.

In view of the fact that our results (a') and (b') are obtained from results (a) and (b), we know it is not always a bad idea to have results requiring large or "unbounded" parallelism, as long as we have some efficient methods to reduce algorithms using large parallelism to algorithms using small parallelism. We believe the two principles of reduction used in this paper for obtaining results (a') and (b') can be applied to other problems. In fact, the algorithm-decomposition principle was first used by Brent [74] and later by Winograd [74] for the problem of parallel evaluation of arithmetic expressions. It turns out that for their problem the problem-decomposition principle can also be used to produce efficient algorithms using a fixed number of processors. This result will not be given in the paper.

In the following section we shall introduce notation used in the paper. We then prove result (a') in Section 3 and result (b') in the final section.
## 2. DEFINITIONS AND NOTATION

We consider parallel algorithms for solving the following nonsingular lower triangular linear system of size n on SIMD machines:

$$ (2.1) \quad \begin{bmatrix} a_{11} & & & \\ & a_{21} & a_{22} & \\ & & \ddots & \\ & & & a_{n1} & \\ & & & a_{n2} & \\ & & & & a_{nn} \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix} = \begin{bmatrix} b_1 \\ b_2 \\ \vdots \\ b_n \end{bmatrix}, $$

which will also be written as

$$ (2.2) \quad A_x = b. $$

We assume that the time is measured by the number of parallel steps of any arithmetic operation.

Define

$T_k(n) = \text{minimum time needed to solve the system (2.1) of size n with k processors.}$

Also define

$M_k(i,p) = \text{minimum time needed to perform i matrix-vector multiplications for pxp matrices and p-vectors, with k processors,}$

and

$A_k(i,p) = \text{minimum time needed to perform i vector additions for p-vectors, with k processors.}$

Then, it is easy to check that

$$ M_{p2}(1,p) \leq \lfloor \log k \rfloor + 1, $$

$$ A_p(1,p) = 1, $$

and hence

$$ (2.3) \quad M_k(i,p) = \left\lfloor \frac{i}{\frac{k}{2}} \right\rfloor (\log p + 1), \quad \text{for } k \geq p^2, $$

$$ (2.4) \quad A_k(i,p) = \begin{bmatrix} i \\ \lfloor \frac{k}{p} \rfloor \end{bmatrix}, \quad \text{for } k \geq p. $$

(All logarithms in this paper are taken to base 2.)
### 3. ALGORITHMS OBTAINED FROM THE ALGORITHM-DECOMPOSITION PRINCIPLE

In this section we shall first show an algorithm using n processors, which results from a direct implementation of Gaussian elimination. Then by decomposing each parallel step of the algorithm into up to $\lfloor \frac{n}{k} \rfloor$ substeps, the algorithm is reduced to an algorithm using k processors, for any k such that $k \le O(n)$.

For $i = 1, ..., n-1$, define

$$L^{(i)} = \begin{bmatrix} 1 & & & & \\ & \ddots & \frac{1}{a_{ii}} & & \\ & & \ddots & \ddots & \\ & & & \frac{1}{a_{ii}} & \\ & & & & 1 \end{bmatrix}$$

Then it is easy to check that

$$L^{(n)} L^{(n-1)} \dots L^{(2)} L^{(1)} A = I$$

Hence by (2.2) we have

$$(3.1) \quad \underline{x} = L^{(n)} L^{(n-1)} \dots L^{(2)} L^{(1)} \underline{b}.$$

Let vectors $\underline{y}^{(1)}, \dots, \underline{y}^{(n+1)}$ be generated by

$$\underline{y}^{(1)} = \underline{b}$$

$$(3.2) \quad \underline{y}^{(i+1)} = L^{(i)} \underline{y}^{(i)}, \quad i = 1, \dots, n,$$

then by (3.1), $\underline{y}^{(n+1)} = \underline{x}$. Hence (3.2) gives an algorithm for solving the system (2.1), which is in fact Gaussian elimination applied to the system. It is easy to verify that

$$(3.3) \quad L^{(i)} \underline{y}^{(i)} = \begin{bmatrix} y_1^{(i)} \\ \vdots \\ y_i^{(i)} \\ \hline \frac{y_i^{(i)}}{a_{ii}} \\ \vdots \\ y_{i+1}^{(i)} - \frac{a_{i+1,i}}{a_{ii}} y_i^{(i)} \\ \vdots \\ y_n^{(i)} - \frac{a_{ni}}{a_{ii}} y_i^{(i)} \end{bmatrix}, \quad i = 1, \dots, n.$$

If n processors are used, for each i the product $L^{(i)} \underline{y}^{(i)}$ can be obtained in three steps, which consist of a parallel division, multiplication and subtraction. Therefore, (3.2) defines a parallel algorithm in a natural way which uses n processors and takes time $3n$.
We can perform each step of the parallel algorithm by k processors instead of n processors in the following obvious way: If the step consists of l operations, l ≤ n, then it is decomposed into $\left\lfloor \frac{l}{k} \right\rfloor$ substeps, each of which consists of at most k operations and is done in a unit of time with k processors. Then the algorithm is reduced to an algorithm using k processors. We say such a reduction is based on the algorithm-decomposition principle. Observe that the new algorithm using k processors can form the product $L_{\mathbf{Y}}^{(i)}$ in (3.3) in time $2\left\lfloor \frac{n-i}{k} \right\rfloor + \left\lfloor \frac{n-i+1}{k} \right\rfloor$. Hence, since

$$ \sum_{i=1}^{n} \left( 2 \left\lfloor \frac{n-i}{k} \right\rfloor + \left\lfloor \frac{n-i+1}{k} \right\rfloor \right) \leq \frac{3}{2k} n^2 + \left( 3 - \frac{7}{2k} \right) n + \frac{2}{k} - 2. $$

We have

$$ (3.4) \quad T_k(n) \le \frac{3}{2k} n^2 + \left(3 - \frac{7}{2k}\right)n + \frac{2}{k} - 2. $$

In the rest of the section we shall improve the bound in (3.4). It is not difficult to check that

(3.2) can be replaced by

$$ (3.5) \quad \mathbf{Y}^{(i+1)} = K^{(i)} \underline{\mathbf{w}}^{(i)}, \quad i = 1, \dots, n, $$

where $\underline{\mathbf{w}}^{(i)}$ is a vector such that

$$ (3.6) \quad \left\{ \begin{array}{ll} \mathbf{w}_j^{(i)} = \mathbf{y}_j^{(i)}, & \text{for } j \neq i, \\ \mathbf{w}_i^{(i)} = \frac{\mathbf{y}_i^{(i)}}{a_{ii}}, & \end{array} \right. $$

and

$$ K^{(i)} = \begin{bmatrix} 1 & \cdots & 1 \\ & -a_{i+1,i}, & 1 \\ & \vdots & \vdots \\ & & -a_{ni} \\ & & & 1 \end{bmatrix}. $$

Thus, (3.5) defines a parallel algorithm using n-1 processors. We now use this algorithm to construct an algorithm using k processors by the algorithm-decomposition principle. Then it is easy to see that the algorithm obtained can form the product $K_{\mathbf{Y}}^{(i)}$ in (3.5) in time $2\left\lfloor \frac{n-i}{k} \right\rfloor$ rather than $2\left\lfloor \frac{n-i}{k} \right\rfloor + \left\lfloor \frac{n-i+1}{k} \right\rfloor$. Since n steps for the division in (3.6) are needed and

$$ \sum_{i=1}^{n} 2 \left\lfloor \frac{n-i}{k} \right\rfloor \leq \frac{n^2}{k} + \left(2 - \frac{3}{k}\right)n + \frac{2}{k} - 2, $$

we have the following

**Theorem 3.1**

$$ T_k(n) \le \frac{1}{k} n^2 + \left(3 - \frac{3}{k}\right)n + \frac{2}{k} - 2. $$
Hence the bound in (3.4) is improved. Note that because the minimum time sequential algorithm takes time $n^2$, the coefficient of the leading term in the theorem is the best possible when k is a constant or $k = o(n)$.

### 4. ALGORITHMS OBTAINED FROM THE PROBLEM-DECOMPOSITION PRINCIPLE

Consider parallel algorithms for solving the system (2.1) of size $n$. Let $N$ be the set of positive integers and $f, g: N \to N$ monotonically increasing functions such that $f^{-1}$ (the inverse function of $f$) exists and $g(n) = o(f(n))$. We shall first show that any parallel algorithm $F$ using $f(n)$ processors can be used to construct a parallel algorithm $G$ using $g(n)$ processors. Our approach in this section is to decompose the problem (i.e. the system (2.1)) rather than the algorithm $F$. We say this approach is based on the problem-decomposition principle. Define

$$ (4.1) \quad p = [f^{-1}(g(n))], \quad m = \lfloor \frac{n}{p} \rfloor $$

and decompose the system (2.1) into a block system of size $m$ with each block of size $p$:

$$ (4.2) \quad \begin{bmatrix} A_{11} \\ A_{21} & A_{22} \\ \vdots & \vdots \\ A_{m1} & A_{m2} & \cdots & A_{mm} \end{bmatrix} \circ \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_m \end{bmatrix} = \begin{bmatrix} b_1 \\ b_2 \\ \vdots \\ b_m \end{bmatrix} $$

(Some padding of the original system (2.1) may be necessary, since we could have $pm > n$.) Hence $A_{ij}$ are $pxp$ matrices and $x_i, b_i$ are $p$-vectors. By (4.2), we have

$$ (4.3) \quad \left\{ \begin{array}{ll} A_{11} & x_1 = b_1, \\ A_{22} & x_2 = b_2 - A_{21-1}, \\ A_{33} & x_3 = b_3 - \sum_{j=1}^{2} A_{3j} x_j, \\ \vdots & \vdots \\ A_{mm} & x_m = b_m - \sum_{j=1}^{m-1} A_{mj} x_j. \end{array} \right. $$

The equations in (4.3) and parallel algorithm $F$ lead to the following algorithm for obtaining $x_1, x_2, \dots, x_m$ by using only $g(n)$ processors:

**Algorithm G**

1. Obtain $x_1$ from the first equation in (4.3) by using algorithm $F$.

2. For $i = 1, \dots, m-1$,
2.1. compute the matrix-vector product $A_{i+1,j} x_j^i$ for $j = 1,...,i$ with $g(n)$ processors, for $i = 1,...,n$.

2.2. form $B_{i+1} = (\sum A_{i+1,j} x_j^i)$ with $g(n)$ processors,

2.3. obtain $x_{i+1}^i$ from the $(i+1)$st equation in (4.3) by using algorithm F.

Since Algorithm G is based on the equations in (4.3), it is obvious to see that the algorithm indeed finds all $x_i^i$, $i = 1,...,m$, and hence solves the original system (2.1). Steps 1 and 2.3 of Algorithm G amount to solving lower triangular linear systems of size $p$. Since from (4.1) $f(p) \le g(n)$, by using Algorithm F in these steps, at most $g(n)$ processors are required. Hence Algorithm G uses $g(n)$ processors.

By Algorithm G, we have

$$ (4.4) \quad \begin{aligned} T_{g(n)}(n) &\le T_{f(p)} + \sum_{i=1}^{m-1} (M_{g(n)}(i,p) + A_{g(n)}(i,p) + T_{f(p)}(p)) \\ &= m T_{f(p)}(p) + \sum_{i=1}^{m-1} (M_{g(n)}(i,p) + A_{g(n)}(i,p)). \end{aligned} $$

From (2.3), (2.4), (4.1) and (4.4), we can easily verify that if $g(n) \ge p^2$ then

$$ (4.5) \quad T_{g(n)}(n) \le O(mT_{f(p)}(p)) + O\left(\frac{n^2}{g(n)} \log p\right). $$

We are ready to prove the following

**Theorem 4.1**

$$ \underline{T_{\lfloor n^r \rfloor}(n)} \le \left\{ \begin{array}{ll} \frac{O(n^{2-r} \log n)}{n^{1-r} \log^2 n} & \text{for } 1 < r < \frac{3}{2}. \\ \\ \frac{O(n^{2-r} \log n)}{n^{1-r} \log^2 n} & \text{for } \frac{3}{2} \le r < 3. \end{array} \right. $$

**Proof**

It has been shown by many people, including Chen and Kuck [74], Heller [74], Munro [74], Orcutt [74], etc., that

$$ (4.6) \quad T_3(n) \le O(\log^2 n). $$

Let $f(n) = n^3$ and $g(n) = \lfloor n^r \rfloor$. Then

$$ p = (\lfloor n^r \rfloor)^{\frac{1}{3}} \quad \text{and} \quad m = O(n^{1-\frac{r}{3}}). $$

By (4.6),

$$ T_{f(p)}(p) \le O(\log^2 p) = O(\log^2 n). $$
Then by (4.5),

$$ (4.7) \quad \underset{[n^r]}{T_r(n)} \le O(n^{\frac{1-r}{3}\log^2 n}) + O(n^{2-r}\log n). $$

The theorem follows from (4.7) and the fact that

$$ 1 - \frac{r}{2} \ge 2-r \quad \text{iff} \quad r \ge \frac{3}{2}. $$

The results in Theorem 4.1 for the case $1 < r < \frac{3}{2}$ can also be obtained from the fact that

$$ T_n(n) \le O(n) $$

instead of (4.6), by a technique similar to the one used in this section.

## 5. CONCLUSIONS

In this paper we have reduced parallel algorithms with large parallelism to parallel algorithms with small parallelism. We used the following two principles for the reduction:

1. Algorithm-Decomposition Principle: We decompose each step of the original algorithm into sub-steps so that each of them can be done in a unit of time with a small number of processors.

2. Problem-Decomposition Principle: We decompose the problem into subproblems so that each of them can be done by the original algorithm with a small number of processors.

These two reduction principles are rather unsophisticated and they probably have been used without emphasis by many people. Yet the ideas work well and give us nontrivial parallel algorithms with small parallelism, as we have demonstrated for the problem of solving triangular linear systems. It is conceivable that they have many other applications to specific problems.

Since in practice we often have to solve a "large" problem with a "small" number of processors, we believe that the study of efficient reductions from large parallelism to small parallelism is one of the fundamental issues in parallel computation.

## References

Brent, R. P. [74]. "The Parallel Evaluation of General Arithmetic Expressions," J. ACM, 21, 201-206.

Chen, S. C. and D. J. Kuck [74]. Time and Parallel Processor Bounds for Linear Recurrence Systems, Department of Computer Science, University of Illinois at Urbana-Champaign.

Heller, D. [74]. On the Efficient Computation of Recurrence Relations, Report, ICASE, NASA Langley Research Center, Hampton, Virginia. Also available as a CMU Computer Science Department Report.

Munro, I. [74]. Private Communication. To appear in Borodin, A. and I. Munro, "Computational Complexity of Algebraic and Numeric Problems," American Elsevier, 1975.

Orcutt, S. E. [74]. Parallel Solution Methods for Triangular Linear Systems of Equations, Report 77, Digital Systems Laboratory, Stanford University.
Winograd, S. [74]. On the Parallel Evaluation of Certain Arithmetic Expressions, Report RC 4808, IBM Research Center, Yorktown Heights, New York.