New Parallel SOR Method by Domain Partitioning.
In this paper we propose and analyze a new parallel SOR method, the PSOR method, formulated by using domain partitioning and interprocessor data communication techniques. We prove that the PSOR method has the same asymptotic rate of convergence as the Red/Black (R/B) SOR method for the five-point stencil on both strip and block partitions, and as the four-color (R/B/G/O) SOR method for the nine-point stencil on strip partitions. We also demonstrate the parallel performance of the PSOR method on four different MIMD multiprocessors (a KSR1, an Intel Delta, a Paragon, and an IBM SP2). Finally, we compare the parallel performance of PSOR, R/B SOR, and R/B/G/O SOR. Numerical results on the Paragon indicate that PSOR is more efficient than  R/B SOR and R/B/G/O SOR in both computation and interprocessor data communication.
Introduction
. The successive over-relaxation (SOR) iterative method is an
important solver for large linear systems [22]. It is also a robust smoother as well as an
efficient solver of the coarsest grid equations in the multigrid method [20]. However, the
SOR method is essentially sequential in its original form. With the increasing use of
parallel computers, several parallel versions of the SOR method have been studied by a
number of authors. Adams and Ortega [1], Adams and Jordan [2], and Adams, Leveque
and Young [3] have written about the multicolor SOR method; Block, Frommer and
Mayer [5] wrote about a general block multicolor SOR method; and White [19] wrote
about the multisplitting SOR method. We also note the papers by Evans [6] and Patel
and Jordan [12], which presented two parallel SOR methods for particular parallel
computers.
The motivation for us to develop a new parallel version of SOR is to provide parallel
multigrid methods with an efficient parallel solver of the coarsest equations [21]. For
some scientific computing problems, the size of the coarsest equations of the multigrid
method is required to be large enough. On the other hand, the convergence rate of
multigrid methods is dependent of the number of grid levels, especially for singular
perturbed problems. So, with the size of the coarsest equations being properly large,
we can improve the performance of the multigrid method, along with avoiding idle
processors on parallel machines [16]. Since domain decomposition is a widely-used
approach in the implementation of parallel multigrid methods on MIMD computers, it
is attractive to develop a parallel SOR method based on the domain decomposition.
Using a domain decomposition technique, we recently proposed a simple parallel
version of the SOR method, the JSOR method [15]. Let a grid mesh domain be
partitioned into p disjoint subgrids. The JSOR method is obtained by concurrently
applying one sweep of the SOR algorithm to each subgrid using the previous iterates
as the "boundary values". By mapping each subgrid into one processor, JSOR can
be easily implemented on MIMD computers with only one step of interprocessor-date-
communication per iteration. However, due to the slow convergence speed, the JSOR
Submitted to SIAM J. Sci. Comput. and 1996 COPPER MOUNTAIN CONFERENCE ON ITERATIVE
METHODS. This work was supported in part by the National Science Foundation through
award number DMS-9105437 and ASC-9318159.
y Courant Institute of Mathematical Sciences, New York University, 251 Mercer Street, New York,
NY 10012, (xie@monod.biomath.nyu.edu).
rarely used as a parallel solver for linear systems. Instead, it can be an
efficient smoother for parallel multigrid methods as shown in [20] and [21].
To improve the convergence properties of JSOR, we modify the pattern of inter-processor
data communication of the JSOR method. As a result, a new type of parallel
SOR method, which is called the PSOR method, is generated from JSOR. In this paper
we prove that PSOR can have the same asymptotic rate of convergence as the sequential
method. Therefore, an efficient parallel SOR method based on domain decomposition
is obtained in this paper for solving the linear systems arising from finite deference
or finite element approximations of partial differential equations.
A model analysis of the PSOR method is presented in Section 2. For the 5-point
stencil model problem, it can be easily proved that PSOR and the SOR using the natural
row-wise ordering have the same asymptotic rate of convergence. In fact, PSOR can
be equivalent to a SOR method using a new ordering. We show that the new ordering
based on either a strip partition as shown in Fig.5 or a block partition as shown in
Fig.6 leads to a consistently ordered matrix. Hence, the SOR theory in [22] follows
that PSOR has the same asymptotic rate of convergence as the SOR using the natural
row-wise ordering.
However, for the general linear system arising from finite deference or finite element
approximations of elliptic boundary value problems, it is difficult to show a global
ordering, with which the SOR method is equivalent to the PSOR method, leads to
a consistently ordered matrix. Instead, we want to give the PSOR method a direct
analysis in this paper.
We present a general description of PSOR in Section 3. Here the PSOR method is
defined on a general domain partition. Then, in Section 4, we prove a basic convergence
theorem for the linear systems in which the corresponding matrix is symmetric positive
definite. The main result of the paper is also presented in Section 4, which shows that the
spectral radius of the PSOR iteration matrix can be the same as that of the SOR iteration
matrix for a wide class of linear systems in which the corresponding sub-matrix on each
subgrid is "consistently ordered". In Section 5, we first confirm the PSOR method
based on either a strip partition or a block partition has the same asymptotic rate of
convergence. Then, we demonstrate the parallel performance of the PSOR method on
a shared memory MIMD computer (a KSR1) and three distributed memory MIMD
computers (the Intel Delta, an Intel Paragon L38 and an IBM POWERparallel System
9076 SP2). The numerical results show that PSOR is very efficient on these distinct
multiprocessor machines. Since the multi-color SOR method is a widely-used parallel
version of the SOR method, we give a comparison between the PSOR method and the
the Red/Black SOR method for the 5-point stencil as well as the four-color SOR method
for the 9-point stencil [3] on parallel computers in Section 6. The numerical results show
that PSOR can have much better performance than both Red/Black and four-color SOR
methods in either floating operations or interprocessor date-communication.
2. The Model Problem Analysis. We consider the 5-point approximation to
the Poisson's equation on a unit square with zero boundary data
(1)
and
@\Omega h . Here u ij denote the approximation of u(x
,\Omega h and
@\Omega h are the sets of the interior and
boundary mesh points, respectively.
Fig. 1. Natural row-wise
ordering.
Fig. 2. Red-Black ordering.
Under some ordering of unknowns, (1) can be written in a matrix form
A being a (m \Gamma matrices. Obviously,
there are many ways to order the unknowns, but the natural row-wise ordering as shown
in Fig.1 and the Red-Black ordering as shown in Fig.2 are two widely-used orderings
in practice.
The SOR method using the natural row-wise ordering generates a sequence of iterates
from a given initial guess u (0)
ij and a real number ! 2 (0; 2) by the form
(2)
which is completely sequential.
The SOR method with the Red-Black ordering, which is usually called the Red-Black
method, takes the form
which can be entirely implemented in parallel on the same colors.
For the model problem, it has been shown that the SOR using the natural ordering
and the Red-Black SOR method have the same convergence rate [22].
JSOR is a simple parallel version of SOR by domain decomposition, which has
been analyzed in [15] recently. For simplicity, we suppose that the grid
mesh\Omega h is
partitioned into p
strips\Omega h;- for and each of them contains at least two
grid lines. We denote
h;- the set of the mesh points of the first grid line
h;-
=\Omega h;-
h;- . Then the JSOR method takes the form
p. Note that (3) is the same as (2).
Clearly, the above JSOR scheme can be implemented in parallel on p processors
by mapping the equations
on\Omega h;- , (3) and (4), into Processor - for
pseudo-code of JSOR on p processors is given by
JSOR Algorithm: For in parallel
Compute
on\Omega h;- by using (3) and (4).
Communicate u (k+1)
ij to other processors as needed.
However, the numerical experiments in [15] showed that the convergence rate of
JSOR is slow down almost linearly with respect to the number p of strips. Hence, it
is rarely used as a parallel solver for linear systems instead of an efficient smoother for
parallel multigrid methods [20] [21].
Notice that the interprocessor data-communication of JSOR on the strip partition
case is as follows:
a) Send u (k+1)
h;- from Processor - to Processor
ij on the last grid line
of\Omega h;- from Processor - to Processor -
If we carry out Step a) before the computation of (4)
h;- , then the update
can be available when we compute the updates u (k+1)
on the last grid line
of\Omega h;- , such that they are defined by
Consequently, a new type of parallel SOR, the PSOR method, is generated from the
method. A pseudo-code of PSOR on p processors is given by
PSOR Algorithm: For in parallel
Compute
h;- by using (3).
h;- from Processor - to Processor
Compute
h;- by using (4) and (5).
ij on the last grid line
of\Omega h;- from Processor - to Processor -
Remarkably, it can be easily shown that PSOR can have the same asymptotic rate
of convergence as the SOR iteration (2) for the model problem.
In fact, the PSOR with a strip partition is equivalent to the SOR scheme using a
new ordering as shown in Fig.5. Let be the matrix associated to the new
ordering. We then can show that A is a consistently ordered matrix as follows.
For the model problem, we have that a ij 6= 0 if and only if mesh node i is adjacent
to node j. According to the definition of a consistently ordered matrix [22], we construct
a disjoint partition of the index set
Wh,1 h,2 h,3
Wh,7 h,8 h,9
Fig. 3. A block partition of the grid
mesh
domain\Omega h .
Fig. 4. The local ordering number
of each
block\Omega h;- , which can lead to a
consistently ordered sub-matrix.
26 27 28 29
28
12192Fig. 5. Global ordering of PSOR on a strip
partition, which can lead to a consistently ordered
matrix.
Fig. 6. Global ordering of PSOR on a block
partition, which also can lead to a consistently ordered
matrix.
such that if a i;j 6= 0 and i This shows
that A is a consistently ordered matrix. Therefore, the SOR theory in [22] follows that
the PSOR with strip partition has the same asymptotic rate of convergence as the SOR
iteration (2).
Similarly, we also can define the PSOR method on a block partition as shown in
Fig.3. At each
block\Omega h;- , PSOR uses a particular local ordering of the mesh points
as shown in Fig.4 in order to communicate the data between processors efficiently. The
vector U - associated
with\Omega h;- has the following partitioning:
where the superscript t denotes a vector transpose,
We can show that such a local ordering within each
sub-grid\Omega h;- can lead to a
consistently ordered matrix. For example, we construct the following subsets for the
ordering shown in Fig.4
such that they satisfy the definition of a consistently ordered matrix.
The implementation of the PSOR method on a block partition is the same as the
strip partition case except that the first component u 1 of each U 1
- should be sent to
other processors as needed as soon as it is updated.
Like the strip partition case, PSOR on a block partition is also equivalent to the
SOR method which is with a new global ordering as shown in Fig.6. We also can show
that the new ordering can lead to a consistently ordered matrix. In fact, we have the
following subsets for the ordering shown in Fig.6:
which satisfy the definition of a consistently ordered matrix. Therefore, the PSOR on
the block partition has the same convergence rate as the SOR using the natural ordering.
3. PSOR Method. In this section, we shall give the PSOR method a general
description. We consider the solution of linear system
which is supposed to arise from a finite element or a finite different discretization of
an elliptic boundary value problem. Here n\Thetan is a sparse matrix, f and u are
real column vectors of order n.
Let\Omega h denote the set of mesh points on which the unknown vector u is defined.
We partition it into p disjoint
subgrids\Omega h;- such that
Based on this mesh partition, the linear system (7) can be written into a block form
Here U - and F - comprise, respectively, the components of u and f associated with the
mesh points on the
subgrid\Omega h;- , and A - is formed from A by deleting all rows except
those corresponding
to\Omega h;- and all columns except those corresponding
to\Omega h;- .
If the entry a ij of matrix A is not zero, then the corresponding mesh points i and j
are said to be coupled. According to this, we denote
h;- the set that comprises the
mesh points
of\Omega h;- that are coupled to the mesh points
of\Omega h;- with -. We then
h;-
=\Omega h;-
h;- , and assume
h;- is nonempty. Thus, each subgrid mesh
\Omega h;- can be partitioned into two disjoint nonempty subsets
\Omega h;-
h;-
Further, we assume that the mesh points of each
subgrid\Omega h;- are ordered in such
a way that the matrix A - related to the
subgrid\Omega h;- is consistently ordered, and the
sub-vector U - satisfies
Here U 1
- and U 2
- comprise the components of u associated with the mesh points on the
h;-
h;- , respectively. One of such examples has been illustrated in (6).
Associated with partition (9), each sub-matrix A - can be written into the following
block form
A 11
A 21
- A 22
Moreover, from the definition
h;- it follows that a nonzero matrix A - with - 6= -
can be reduced to
for -, and A
A 21
for -:
Here A 12
- and A 21
are nonzero.
In fact, let a ij be a nonzero entry of A - with -. The definition of A - follows
that i
2\Omega h;- and j
2\Omega h;- . In other words, mesh point i
of\Omega h;- is said to be coupled
to mesh point j
in\Omega h;- . According to the definition
h;- , we have that mesh point i
belongs
h;- . Similarly, noting that the symmetry of A gives a 0, we can
show mesh point j is
h;- . Hence, a ij must be an entry of A 12
- . This shows A 12
- is
nonzero. On the other hand, other sub-matrices of A - must be zero.
Consequently,
#/
A 21
A 12
A 21
Hence, (8) is written into an equivalent form
A 11
A 21
- A 22
#/
A 12
A 21
Let the k-th PSOR iterate be denoted by u
Suppose that u (k) is given. We can construct the following linear system
h;-
A 11
A 12
and define U 1;(k+1)
- as the value of one SOR iterate for solving (11).
Clearly, the p linear systems in (11) are independent. So, we can compute U 1;(k+1)
in parallel on p processors. When they are available, we communicate U 1;(k+1)
- for
to other processors as needed, such that we get a linear system
h;-
as follows:
A 22
A 21
We then define U 2;(k+1)
- as the value of one SOR iterate for solving (12). Obviously, the
computation of U 2;(k+1)
- can also be implemented in parallel on p processors. Therefore,
one step of PSOR iteration is defined by the following algorithm.
PSOR Algorithm: For
- as the value of one SOR iterate for solving (11)
h;- .
Communicate U 1;(k+1)
- to other processors as needed.
- as the value of one SOR iterate for solving (12 )
h;- .
Communicate U 2;(k+1)
- to other processors as needed.
Let c -;1 and c -;2 denote the right hand sides of (11) and (12), respectively. With
the matrix expression of SOR [22], we express U 1;(k+1)
- and U 2;(k+1)
- as follows:
U
where D -;j is the diagonal matrix of A jj
I is an identity matrix, and L -;j and U -;j
are respectively strictly lower and upper triangular matrices such that
-;j A jj
2:
We then write the above two equations into a block form
where D - is the diagonal matrix of A - , and L - and U - are respectively strictly lower
and upper triangular matrices such that
Here we also used the following equalities:
A 12
A 21
Since the p block equations in (13) are independent, we can number the subgrids
f\Omega h;- g in an arbitrary ordering to form a global matrix . For the simplicity,
in the following, we still denote
Further, we introduce the following notation:
A
A 31 A
A p1 A p2
and
Using them, we can write (13) into a block form
or, equivalently,
Since the determinant of I \Gamma !(B+M) is 1, matrix I \Gamma !(B +M) is nonsingular. So we
can solve (15) for u (k+1) obtaining a matrix expression of the PSOR iterate as follows
e
and F We refer to M PSOR (!) as the PSOR iteration
matrix.
the SOR method reduces to the Gauss-Seidel method. Similarly, we refer
to the PSOR with as the PGS (parallel Gauss-Seidel) method. From (17) it
follows the PGS iteration matrix
4. PSOR Analysis. In this section we study the convergence of the PSOR
method. We denote by ae(A) the spectral radius of matrix A, which is defined as
the maximum of the moduli of the eigenvalues of A. The determinant of a matrix A
is denoted as det(A). It is well known that a necessary and sufficient condition for the
convergence of a linear stationary iteration
is that
Here G is the iteration matrix. Hence, we only need to study ae(M PSOR (!)) for the
convergence of the PSOR method.
We first establish a necessary condition for the convergence of the PSOR method.
theorem 1. Let M PSOR be the PSOR iteration matrix. If ae(M PSOR
the PSOR method is convergent, then 2:
Proof. By a well-known theorem of linear algebra, we know that det(A) is equal to
the product of the eigenvalues of matrix A. With the definitions of B; C; M , and N , we
have
Y
and
Y
Thus,
Hence, if ae(M PSOR 2. This completes
the proof of Theorem 1.
A sufficient convergence condition for the PSOR method is given in the following
theorem.
theorem 2. Let M PSOR be the PSOR iteration matrix. Suppose that A is a symmetric
positive definite matrix.
Proof. Let x be an eigenvector of M PSOR (!), and - be its et x be an eigenvector
of M PSOR (!), and - be its corresponding eigenvalue such that
By (18), the above equality is written as
Multiplying x t D to the both sides of (21), we find
r 2 are real numbers, and
\Gamma1: Clearly, the symmetry of matrix A gives that
Thus,
By the positive definiteness of A, we have that for any x 6= 0,
Hence, j-j ! 1. This follows that ae(M PSOR (!)) ! 1. This completes the proof of
Theorem 2.
We now turn to show that both PSOR and SOR can have the same asymptotic
convergent rate. Essentially, we only need to assume that each sub-matrix A - with
consistently ordered in the following theorems. However, due to
the independence of the p block equations in (13), we can chose a ordering of the
subgrids
f\Omega h;- g as shown in Fig.7 so that the corresponding global matrix
consistently ordered too. Hence, in the following we can assume that A is consistently
ordered.
Fig. 7. A Global ordering which results a matrix partition A = (A- ), such that both A and all
of the sub-matrices A- are consistently ordered for the 5-point stencil.
Similar to Theorem 3.3 (p.147) of the SOR theory in [22], we have the following
theorem for PSOR.
theorem 3. Let A have a block partition
. Let ff and k be real numbers. If A and A - for all are
consistently ordered, then for ff 6= 0 and for all k,
is independent of ff.
Proof. Let be a permutation defined on the integers
is an entry of D \Gamma1 Bg; is an entry of D \Gamma1 Cg;
is an entry of D \Gamma1 Ng; and is an entry of D \Gamma1 Mg:
The general term of \Delta is
Y
a ioe(i) ff nB +nM \Gamman C \Gamman N k n\Gamma(n B +nM +nC +nN
are, respectively, the number of values of i such that (i; oe(i)) 2
, such that (i; oe(i)) 2 such that (i; oe(i)) 2 TM and such that (i; oe(i)) 2 TN . Since
only if there exists one a only need to consider the terms of
a
Let nL and nU be , respectively, the number of values of i such that i ? oe(i) and
such that i ! oe(i). Obviously, we have
where l k and - k are, respectively, the number of values of i such that i ? oe(i) and such
that well as such that a ioe(i) is an entry of A kk .
Since A and A - for all are consistently ordered, from the proof of
Theorem 3.3 (p.147) in [22] we have
It follows that
(l
and then
Therefore,
This shows that t(oe) is independent of ff and the proof of Theorem 3 is completed.
Let D be a diagonal matrix of A, L and U are respectively strictly lower and upper
triangular matrices such that
Obviously, we have
From [22] we know that the SOR iteration matrix
and the Jacobi iteration matrix
The following theorem shows that PSOR and SOR have the same asymptotic convergent
rate.
theorem 4. Let M J ; M SOR and M PSOR be the iteration matrices of the Jacobi,
SOR and PSOR methods, respectively. We assume that matrix A has a block partition
such that both A and the sub-matrix A - for all are
consistently ordered. If M J has real eigenvalues and ae(M J
is the optimal relaxation parameter, which has the following expression:
Proof. Let - 6= 0. With (18) and (19), we have
det
det
det
using Theorem 3 in the last step.
If - is is an eigenvalue of M PSOR , then
det
which implies that
is an eigenvalue of M J . Conversely, if - ? 0 is an eigenvalue of M J , then there
always exists an eigenvalue of M PSOR (!) which satisfies (25), and the value of - can
be computed from the equation
From [22] we know that the eigenvalues - of the SOR iteration matrix M SOR also satisfy
equation (26). Hence, the sets of the eigenvalues of both M PSOR and M SOR are the
same. Therefore, we have
Noting that ae(M SOR (!)) has been shown to have expression (24) in [22], we conclude
that (24) holds for the PSOR method also. This completes the proof of Theorem 4.
5. Numerical Examples. In this section, we first confirm that the PSOR method
on either a strip partition or a block partition has the same asymptotic convergence rate
as the SOR method. We then highlight the parallel performance of the PSOR method
on four large parallel MIMD machines: a KSR1, the Intel Delta, an Intel Paragon L38
and an IBM POWERparallel System 9076 SP2. The numerical results demonstrate the
high efficiency of the PSOR method on parallel MIMD computers.
Two PSOR programs were written in Pfortran [4] and MPI (a Message-Passing
Interface Standard) [17], respectively. The PSOR programs were compiled with optimization
level \GammaO2 on the KSR1 and the SP2 and \GammaO4 on the Intel Delta and the
Paragon, respectively. In the Pfortran program, the CPU time was computed by using
the dclock() system routine on the Intel Delta and the Paragon, the user seconds()
system routine on the KSR1, and the mclock() on the SP2. In the MPI program, we
used MP I WTIME() function of MPI.
In the figures and tables, we have used the following annotations.
Linear Time for p processor case is defined by T (1)
is the CPU time
needed to solve the given problem on one processor. This stands for the ideal case.
T otal Time represents the CPU time spent from the beginning of the iteration until
either (28) is satisfied or the limited number of iterations is reached. It does not include
the CPU time spent on the calculation of f , the initial guess, and the input/output
of data. Comm: Time represents the CPU time spent on the interprocessor data
communication. Comp: Time represents the CPU time spent on the computation of
the iteration, including the L 2 -norm of the residual. Other
Comm: Time \Gamma Comp: Time, which includes the time spent on the global summations
in the computation of (28). This also indicates the accuracy of our time measurements.
We first considered the model problem (1) with f(x; sin -x sin -y. Clearly,
sin -y is the exact solution of the model problem. For simplicity, we fixed the
1=513, the relaxation parameter and the initial guess u
for all of the numerical experiments in this section. In this experiment, we also fixed
the number of PSOR iterations as 1000. With the Pfortran program, we implemented
PSOR on 1, 4, processors of Paragon, respectively. In the strip case, the
grid mesh was partitioned into p strips with equal sizes when PSOR was implemented
on p processors. In the block case, the grid mesh was divided into 2 \Theta 2, 4 \Theta 4, 8 \Theta 8
and blocks with equal sizes when PSOR was implemented on 4, 16 , 64 and 256
processors, respectively.

Table

1 shows that the PSOR on either a strip partition or a block partition has
the same asymptotic rate of convergence as the corresponding sequential SOR method.
Due to the different computing formula of PSOR at different groups of processors, the
relative residual
and the relative error
would be a little different. Here k \Delta k 2 is the L 2 norm, u (1000) is the 1000th iterate of
PSOR, and u is the exact solution on the grid mesh. From the table we also see that
the performance of PSOR on the strip partition is better than that of PSOR on the
block partition. Hence, we only considered the PSOR method on a strip partition in
the reminder of the paper.
We next considered the performance of PSOR on different parallel machines. We
implemented the PSOR method using
which leaded to the following stopping criterion for the PSOR iteration:

Table
The parallel performance of PSOR on either a block partition or a strip partition for solving the
5-point stencil of Poisson equation with
Processor Total Time Comm. Time
strips blocks strips blocks strips blocks strips blocks
Number of Processors (on KSR1)
Time
in
Seconds
Total Time : solid line
Linear Time dotted line
Comp. Time
Other Time
Comm. Time
Number of Processors (on DELTA)
Time
in
Seconds
Total Time the solid line
Linear Time : the dotted line
Comm. Time
Other Time
Comp. Time
Fig. 8. The parallel performance of the
PSOR method for solving (1) with
a KSR1. The floating point performance of the
on the KSR1 is 6.58
Mflops.
Fig. 9. A parallel performance of the PSOR
method for solving (1) with 1on the Intel
Delta. The floating point performance of the SOR
program (i.e. on the Intel Delta is 6.09
Mflops.
Number of Processors (on Paragon L38)
Time
in
Seconds
Total Time the solid line
Linear Time : the dotted line
Comm. Time
Other Time
Comp. Time
Number of Processors (on SP2)
Time
in
Seconds
Total Time the solid line
Linear Time the dotted line
Comm. Time
Other Time
Comp. Time :
Fig. 10. A parallel performance of the PSOR
method for solving (1) with 1on the
Paragon L38. The floating point performance of
the SOR program (i.e. on the Paragon is
7.678 Mflops.
Fig. 11. A parallel performance of the PSOR
method for solving (1) with 1on the SP2.
The floating point performance of the SOR program
(i.e. on the SP2 is 40.63 Mflops.
Number of Processors
Time
in
Seconds
Total Time the solid line
Comm. Time the dotted line
Number of Processors
Time
in
Seconds
Total Time : the solid line
Comm. Time the dotted line
Fig. 12. A comparison of the performances
of the PSOR method on the KSR1 and the Intel
Delta.
Fig. 13. A comparison of the performances
of the PSOR method on the Paragon L38 and the
SP2.
As a result, a lot of computer time in checking the convergence of PSOR could be saved.
The test problem was the model problem (1) with f(x; 1. The same codes in
Pfortran were run on the KSR1 and the Intel Delta, while the same codes in MPI were
run on the Paragon and the SP2. The results were reported in Fig. 8 to Fig.11.
Fig. 8 to Fig. 11 provide the CPU time for the PSOR method on the KSR1, the Intel
Delta, the Paragon and the SP2, respectively, as a function of the number of processors.
The total numbers of PSOR iterations, satisfying the convergence criterion (28), for 1, 2,
4, 8, 16, 32, 64 processors are 1027, 1026, 1025, 1023, 1018, 1008 and 942, respectively.
From the figures we see that Comp: Time is almost the same as Linear Time, and
both Comm: Time and Other Time are very small. These results demonstrated that
the PSOR method is an efficient parallel version of the SOR method using the optimal
relaxation parameter.
Fig. 12 presents a comparison of the performance of PSOR on the KSR1 and the
Intel Delta. From this we see that the PSOR method has largely the same performance
on these two different architecture machines. We also compare the performance of
PSOR on the Paragon and the SP2 in Fig. 13. The SP2 is the latest distributed
memory machine from IBM. From the figure we see that it is more powerful in both
floating operations and interprocessor message passing than the Paragon.
6. Comparison of PSOR with Multicolor SOR. In this section we present a
comparison of the parallel performance of the PSOR method with the multicolor SOR
method on the Paragon. The multicolor SOR method for solving the 5-point stencil (1)
is the Red-Black SOR method, which has been described in Section 2.
We also considered the 9-point approximation to the Poisson equation as follows:
and
@\Omega h .
Red Black Glue Orange
Fig. 14. Four color ordering for the 9-point stencil
Number of Processors (on Paragon)
Time
in
Seconds
Total Time: Solid line
Comm. Time: Dotted line
Red/Black
Number of Processors (on Paragon)
Time
in
Seconds
Total Time: Solid line
Comm. Time: Dotted line
Fig. 15. A Comparison of the parallel perfor-
mation of PSOR with the Red/Black SOR method
Fig. 16. A Comparison of the parallel perfor-
mation of PSOR with the four-color SOR method
for the 9-point stencil
The multi-color SOR method involves four colors for the 9-point stencil. There
are many distinct four-color ordering, but we only considered one of them, which is
illustrated in Fig.14 with ordering R=B=G=O. This ordering has been shown to be
equivalent to the natural ordering in [2].
For the simplicity, we fixed the total number of iterations as 1000 for both PSOR
and the multi-color SOR. We also set for the numerical
experiments. The results were reported in Table 2 and 3. We also draw figures in
Fig.15 and Fig.16 according to the data in Table 2 and 3. Residual term in Table 2 and
3 was computed by using kf \Gamma Au (1000) k 2 . Note that the PSOR method with
reduces to the sequential SOR method. Since the PSOR iteration formula is different
from the different groups of processors, different residuals were obtained. In contrast,
the multicolor SOR had the same residual for all the groups of processors.
Fig.15 and Fig.16 clearly show that PSOR has a much better parallel performance
than the multi-color SOR. Here the Red-Black and four-color SOR methods were programmed
based on a strip partition, which is one of the efficient ways to program them
on a MIMD computer. Due to this, they required two and four steps of interprocessor

Table
A comparison of the parallel performance of PSOR with the Red-Black SOR method for solving
the 5-point stencil of Poisson equation with Paragon). Here the L2-norm of the
residual of the Red-Black SOR was 2:57 \Theta 10 \Gamma5
Processor Total Time Comm. Time Residual (PSOR)

Table
A comparison of the parallel performance of PSOR with the Four-Color SOR method for solving
the 9-point stencil of Poisson equation with Paragon). Here the L2-norm of the
residual of the Four-Color SOR was 4:88 \Theta 10 \Gamma6
Processor Total Time Comm. Time Residual (PSOR)
8 43.2 78.5 1.76 6.5 2:61 \Theta 10 \Gamma6
128 4.4 11.2 1.82 6.52 5:65 \Theta 10 \Gamma6
data communication per iteration, respectively. In contrast, each PSOR iteration only
exchanges interprocessor data once for both 5-point and 9-point stencils. Hence, PSOR
took much fewer CPU time on interprocessor data-communication than the multi-color
SOR. Moreover, our experiments indicated that PSOR also took much few CPU time
in floating operations than the multi-color SOR in either the sequential case (i.e. one
processor case) or the parallel case.

Acknowledgments

. The author would like to thank his advisor Professor L.
Ridgway Scott for valuable discussions and his continuous support. He is also grateful
to Professor Tamar Schlick for her comments and support. Access to the the Delta and
a Paragon L38 from Intel Corporation and an IBM POWERparallel System 9076 SP2
has been provided by the center for Advanced Computing Research at Caltech and the
Theory Center.



--R

A multi-color SOR method for parallel computation


Analysis of the SOR iteration
for the 9-point Laplacian
A Parallel Dialect of Fortran
coloring schemes for the SOR method on local
memory parallel computers.
Parallel SOR Iterative Methods
Farhat: Thesis in Civil Engineering at Berkeley

Determination of stripe structures for finite element matrices

The SOR method on parallel computers.
Solution of partial differential equations on vector and parallel


method on a multiprocessor
Designing Efficient Algorithms for Parallel Computers
Purdue University
Parallel Linear Stationary Iterative Methods

Parallel U-Cycle Multigrid Method
University of Tennessee
Matrix Iterative Analysis
Multisplittings and Parallel Iterative Methods
Mechanics and Engineering
Parallel Multiplicative Smoother and Analysis

Iterative Solution of Large Linear System
--TR
SOR;JSOR;nonmigratory permutation;convergence analysis;PSOR;multicolor SOR;parallel computing