Parallel dichotomy algorithm for solving tridiagonal SLAEs

0.36 MB·

by  Andrew Terekhov

Parallel dichotomy algorithm for solving tridiagonal SLAEs
A.V. Terekhov
Institute of Computational Mathematicsand Mathematical Geophysics, 630090,Novosibirsk,Russia
Budker Institute of Nuclear Physics, 630090, Novosibirsk, Russia
9
Abstract
0
0
Aparallelalgorithmforsolvingaseriesofmatrixequationswithaconstanttridiagonalmatrixanddifferent
2
right-hand sides is proposed and studied. The process of solving the problem is represented in two steps.

n The first preliminary step is fixing some rows of the inverse matrix of SLAEs. The second step consists in
a
calculating solutions for all right-hand sides. For reducing the communication interactions, based on the
J
formulated and provedmain parallelsweep theorem, we propose an originalalgorithmfor calculating share
0
components of the solution vector. Theoretical estimates validating the efficiency of the approach for both
2
thecommon-anddistributed-memorysupercomputersareobtained. Directanditerativemethodsofsolving

] a 2D Poissonequation, which include procedures of tridiagonalmatrix inversion,are realizedusing the mpi
A
technology. Results of computational experiments on a multicomputer demonstrate a high efficiency and
N
scalability of the parallel sweep algorithm.
.
h
Key words: Parallelalgorithm, Tridiagonalmatrix algorithm (TDMA), Thomas algorithm, Sweep
t
a Method, Poissonequation, Alternating Direction Method, Fourier Method
m PACS: 02.60.Dc, 02.60.Cb, 02.70.Bf, 02.70.Hm
[

2
v
9 1. Introduction.
5
8 The progress in numerical methods of solving ”complex problems” is impossible without applying pow-
2 erful parallel computer systems. Thus, it is necessary to investigate numerical algorithms that allow for
.
1 efficient parallel implementation.
0 The problem of solving tridiagonal systems of linear algebraic equations [1, 2, 3] is one of the most
9 frequently solved problems in computational mathematics. The tridiagonal SLAEs arise in three-point
0
approximation of problems for ordinary differential equations of second order with constant and variable
:
v coefficients, and also in realization of difference schemes for equations in partial derivatives [4, 5, 6]. As a
i rule, tridiagonal SLAEs are solved using various versions of the direct difference equation method, that is,
X
a sweep method: monotonic, nonmonotonic, flux and orthogonal[1, 2, 7, 8, 9].
r
a Development and improvement of parallel sweep algorithms is of great interest, which is confirmed by
numerous publications [3, 10, 11, 12, 13, 14, 15] concerned with this difficult problem. Analyzing papers
dealing with this topic, we can conclude that presently available parallelsweep algorithms are insufficiently
efficient and, what is more important, they are insufficiently scalable. The primary cause is that efficient,
in a theoretical aspect, parallel algorithms realized on different multiprocessor computer systems become
disadvantageous due to the presence of such operations as communications and synchronizations.
Solving problems by finite-difference methods frequently requires to solve not one, but a series of tridi-
agonal SLAEs with different right-hand sides, the number of problems in the series can reach thousands.
Thus,the problemof designinganefficient parallelsweepalgorithmfor solvingseriesof tridiagonalsystems
of equations deserves consideration.
Email addresses: [email protected] (A.V.Terekhov)
Preprint submitted toElsevier January 20, 2009
In this paper, we propose a new approach to designing a parallel sweep algorithm for solving a series
of tridiagonal SLAEs with a constant matrix and different right-hand sides. The process of solving the
problem is subdivided into two steps. The first, preliminary step consists in fixing some rows of the SLAE
inverse matrix by means of a sequential procedure. Then follows calculation of solutions for all right-hand
sides; doing so, for increasing the algorithm efficiency using the formulated and proved main parallel sweep
theorem,we proposedanoriginalalgorithmfor calculating individual components from the solutionvector.
2. Statement of the problem.
Theseriesofsystemsofalgebraiclinearequationswithasymmetricalconstanttridiagonalmatrixmeans
AXn =Fn, n=1,...,N. (1)
b a 0
1 1
(cid:13) a1 b2 a2 (cid:13)
(cid:13) (cid:13)
(cid:13) a2 b3 a3 (cid:13)
A=(cid:13)(cid:13)(cid:13) ... ... ... (cid:13)(cid:13)(cid:13)
(cid:13) (cid:13)
(cid:13) a b a (cid:13)
(cid:13) n−2 n−1 n−1 (cid:13)
(cid:13) 0 a b (cid:13)
(cid:13) n−1 n (cid:13)
(cid:13) (cid:13)
, where N is the number of problems(cid:13)in the series. (cid:13)
Assuming that system (1) is nondegenerate, we aim for designing a parallel algorithm for solving the
problem and subsequent realizing on a multicomputer [16, 17].
Datadecomposition.Thecomputationalandcommunicationcomplexityofaparallelalgorithm,hence,
the execution time depend drastically on the way of decomposition of problem data. Let us dwell on the
problem of mapping the data of problem (1) onto a set of processor elements (PEs).
Designing algorithms for distributed-memory supercomputers, it is necessary to take into account the
fact that local data (data in the local memory of the same PE) are accessed much faster than data on a
distant PE. Thus, even during designing the parallel sweep algorithm, we perform computing such as to
minimize communication interactions by means of increasing the portion of local calculations (calculations
performed with local data).
We ground the proposed parallel sweep algorithm on the following specification of data distribution
between PEs:
1. Assuming thatthe number ofPEsis p,divide vectorsFn andXn into subvectorsQn , Un as follows1:
i i
Fn =(cid:0)Qn1,Qn2,...,Qnp(cid:1)T =(cid:16)f1n,f2n,...,fsnize{Fn}−1,fsnize{Fn}(cid:17)T, (2)
T
Xn =(cid:0)Un1,Un2,...,Unp(cid:1)T =(cid:16)xn1,xn2,...,xnsize{Xn}−1,xnsize{Xn}(cid:17) . (3)
2. The sizes of Qn and Un are chosen under the conditions
i i
size Qn =size Un 2 i=1,...,p
{ i} { i}≥
pi=1size{Qni}= pi=1 size{Uni}=size{Fn}=size{Xn}
P P
3. Demand that the pair of subvectors (Qn ,Un ) belong to PE number i .
i i
4. The row of A number j is on the same PE as the pair of elements xn,fn from (2),(3).
j j
(cid:0) (cid:1)
WeshouldnotethatthespecificationofdecompositionofXn,Fn,andArulesoutabsolutelyduplication
of the problem data. Exactly for this distribution we will design the parallel sweep algorithm for solving of
problem (1).
1ThenumberofelementsofavectorV,isdenotedassizesize{V}.
2
3. Parallel sweep algorithm
3.1. Basic algorithm
Lemma 1. Let the tridiagonal system of linear equations AXn =Fn be divided into subsystems of the form
AiUni =Qni, i=1,...,p, (4)
1 0
(cid:13)(cid:13) ali+1 bli+2 ali+2 (cid:13)(cid:13)
Ai =(cid:13)(cid:13)(cid:13)(cid:13) ali+2 bl.i.+.3 a.li.+.3 ... (cid:13)(cid:13)(cid:13)(cid:13),
(cid:13) (cid:13)
(cid:13) a b a (cid:13)
(cid:13) li+ti−2 li+ti−1 li+ti−1 (cid:13)
(cid:13) 0 1 (cid:13)
(cid:13) (cid:13)
(cid:13) (cid:13)
(cid:13) i−1 (cid:13)
ti =size{Uni}, li = ti
kX=1
according to the proposed approach, and let we know the values of elements 2 first Un last Un ,
{ i} { i}
then systems (4) can be solved independently, equality (3) will be fulfilled.
The proof of the lemma follows from the tridiagonal matrix of the structure itself.
Algorithm 1. Based on Lemma 1, we can propose the following parallel sweep algorithm
1. Decompose initial system (1) into subsystems of form (4).
2. Find solutions in boundary elements first Un , last Un .
{ i} { i}
3. Compute Xn , by solving independently subsystems.
Thus, the parallel sweep algorithm enables to reduce the solution of problem (1) to p independent
subproblems of form (4), if values of first Un and last Un are known. However, we have still to solve
{ i} { i}
the issue of the efficient way of computing ”boundary” elements, i.e., it is necessary to design a parallel
algorithm for computing an individual component of the solution vector Xn.
3.2. Computing arbitrary solution component.
Lemma 2. Let B be the symmetrical tridiagonal matrix. Then the value of the kth component of solution
vector (Y) of the equation BY =F can be found as
k
(Y) =GkTF, (5)
k
where the vector Gk is the solution of the following system of equations
BGk =ek, (6)
ek is the unit vector.
Proof. Let B.k be designated by the kth column of B, and Bk. by the kth row, respectively
B1.
B =(B.1,B.2,...,B.n)= B..2.. 
 
 Bn. 
 
2Thefirstandlastelements ofsomevector V,aredenotedasfirst{V}andlast{V}.
3
By virtue of the definition of the inverse matrix BB−1 =I, the solution of system (6) is the kth column of
B−1 and the kth row (under the condition of matrix symmetry).
From this follows
Gk =B−1ek =B−.k1 = B−k.1 T, (Y)k =B−k.1F=GTkF. (7)
(cid:0) (cid:1)
Thus, the lemma has been proved.
Letusrepresentthealgorithmforcomputingthearbitrarysolutioncomponentfortheseriesoftridiagonal
equations of form (1).
Algorithm 2. For calculating M different components of the solution vector (Xn) , n = 1,...,N;
km
m=1,...,M from series of equations (1), :
1. Find Gk , m=1,...,M, as the solution of AGk =ek .
m m m
2. Define (Xn) for the whole series as
km
(Xn)km =GTkmFn, n=1,...,N, m=1,...,M.
Thus,Algorithm2makesitpossibleto findaseparatecomponentofthe solutionvector. Itisimportant
that for different km the values of (Xn)km can be calculatedindependently. We shouldnote that the vector
Gk does not depend on the right-hand side of (1) , hence, it may be determined once for all Fn.
m
Let us consider Algorithm 2 in application to the parallel sweep algorithm. According to Algorithm 1,
at the second step it is required to calculate values of elements first Un , last Un for all PEs. The
{ i} { i}
specified data decomposition assumes that only one subvector Un and one subvector Qn are placed on a
i i
single PE, therefore, each PE will compute two elements (first Un ,last Un ) from the solution vector.
{ i} { i}
Parallel realization of Algorithm 2 entails a great difficulty, namely, according to (5), each PE will have to
perform about O(size Xn ) operations regardless of the number of involved computational resources.
{ }
This causes the problem of modifying Algorithm 2 in such a manner that the number of operations per
PE is about O(size Un ).
{ i}
3.3. The main parallel sweep theorem.
We willstartdesigning anefficient algorithmforcomputing ”boundary”elements; for illustration,letus
consider the following boundary-value problem
d2ϕ
= ρ(x); ϕ(x )=0, ϕ(x )=0 (8)
dx2 − 0 1
Asisknown[18],thesolutionofproblem(8)mayberepresentedintheintegralformviathecorresponding
Green function 3
x1
ϕ(x)= G(x,s)ρ(s)ds. (9)
Z
x0
Let us partition the interval (x ,x ) by three points x ,x ,x and define the right-hand side of
0 1 1/4 1/2 3/4
{ }
(8) as
0, x x x
0 1/4
≤ ≤
 k(x), x <x<x
ρ(x)= 1/4 1/2 (10)
 0, x1/2 x x3/4
≤ ≤
0, x x x
3/4 1
According to (9), the solution at the points of partition≤ing ≤may be defined as
3Inthiscase,werestrictourselvesonlytothefactofitsexistence.
4
x x x x x
0 1/4 1/2 3/4 1
Figure1:
x1/2 x1/2
ϕ(x )= G(x ,s)k(s)ds, ϕ(x )= G(x ,s)k(s)ds,
1/4 1/4 1/2 1/2
Z Z
xx1/14/2 x1/4 (11)
ϕ(x )= G(x ,s)k(s)ds.
3/4 3/4
Z
x1/4
Another wayforfinding the solutionofequation(8)atthe pointx without calculatingthe integralof
3/4
form (9), is as follows: since the point x belongs to the interval (x ,x ) he solution at it should satisfy
3/4 1/2 1
the equation
d2ϕ˜
=0, (12)
dx2
with the boundary conditions
x1/2
ϕ˜(x1/2)= G(x1/2,s)k(s)ds, ϕ˜(x1)=0. (13)
Z
x1/4
It is extremely important (from the viewpoint of computation) that the solution of problem (12),(13) is
represented analytically [19]
x x
1
ϕ˜(x)=ϕ(x ) − . (14)
1/2 x x
1 1/2
−
Thus, the solution of (8) with the right-hand side (10) at the points x ,x ,x is as follows
1/4 1/2 3/4
{ }
x1/2 x1/2
ϕ(x )= G(x ,s)k(s)ds, ϕ(x )= G(x ,s)k(s)ds,
1/4 1/4 1/2 1/2
Z Z
x1/4 x1/4 (15)
x x
3/4 1
ϕ(x )=ϕ(x ) − .
3/4 1/2 x x
1 1/2
−
Comparison of (11) with (15) in their computation complexity shows evident advantage of the latter
because it is required to compute less integrals of form (9).
For the arbitrary function ρ(x) we summarize the obtained result as a theorem.
Theorem 1. It is required to find the solution of boundary problem (8) at points with the coordinates
x x <x <x , x <x , i=1,...,N 1 . Then the following identity takes place
i 0 i N i i+1
{ | − }
i N
x x x x
ϕ(x )= αR i− N + αL i− 0, i=1,...,N 1, (16)
i j x x j x x −
Xj=1 j − N j=Xi+1 j − 0
xi
αR = G(x ,s)ρ(s)ds, i=1,...,N 1,
i Z i −
xi−1
(17)
xi
αL = G(x ,s)ρ(s)ds, i=2,...,N.
i Z i−1
xi−1
5
Let us formulate and prove the main parallel sweep theorem.
Theorem 2. Let we have the nondegenerate system of linear algebraic equations with the tridiagonal matrix
AX=F of dimension n. Then for each solution vector component from the set
Ω= (X) 1<n <n, n <n , i=1,...,p n (18)
ni| i i i+1 ≤
(cid:8) (cid:9)
the following identity holds true
i p+1
(X) = βR ZR + βL ZL , (19)
ni Xj=1 nj(cid:16) nj(cid:17)ni j=Xi+1 nj(cid:16) nj(cid:17)ni
ni−1
βL = (F) A−1 , (n =n+1)
ni j ni−1,j p+1
j=Xni−1
ni−1 (20)
βR = (F) A−1 , (n =1)
ni j ni,j 0
j=Xni−1
BLZL =eL, BRZR =eR, (21)
k k k k
where
b c 0 1 0 0
1 1
(cid:13) a1 b2 c2 (cid:13) (cid:13) ak bk+1 ck+1 (cid:13)
(cid:13) (cid:13) (cid:13) (cid:13)
(cid:13) a2 b3 c3 (cid:13) (cid:13) ak+1 bk+2 ck+2 (cid:13)
BkL =(cid:13)(cid:13)(cid:13) ... ... ... (cid:13)(cid:13)(cid:13), BkR =(cid:13)(cid:13)(cid:13) ... ... ... (cid:13)(cid:13)(cid:13) (22)
(cid:13) (cid:13) (cid:13) (cid:13)
(cid:13) a b c (cid:13) (cid:13) a b c (cid:13)
(cid:13) k−2 k−1 k−1 (cid:13) (cid:13) n−2 n−1 n−1 (cid:13)
(cid:13) 0 0 1 (cid:13) (cid:13) 0 a b (cid:13)
(cid:13) (cid:13) (cid:13) n−1 n (cid:13)
(cid:13) (cid:13) (cid:13) (cid:13)
(cid:13) (cid:13) (cid:13) (cid:13)
eR =(1,0,0,...,0)T, eL =(0,...,0,0,1)T
Proof. Let
p p
X= Xj, F= Fj (23)
Xj=1 Xj=1
AXj =Fj
Define elements of F as follows
j
0, i<n
j−1

Fj i = (F)i, nj−1 ≤i<nj (24)
(cid:0) (cid:1)
0, i n
Then the solutions (X ) , (X ) areunambiguou≥slyjdefined by the following expressions
j nj−1 j nj
6
nj−1
(X ) =A−1 F = A−1 (F ) =βL
j nj−1 nj−1. j nj−1,k j k ni
k=Xnj−1
(25)
nj−1
(X ) =A−1F = A−1 (F ) =βR
j nj nj. j nj,k j k ni
k=Xnj−1
According to Lemma 1, define (X ) for 1 k j 2 as
j nk ≤ ≤ −
BL XL =βL eL,
nj j nj
T (26)
XL = xL ,...,xL ,...,xL ,...,xL ,
j (cid:16) jn1 jn2 jnj−2 jnj−1(cid:17)
and for j+1 k p as
≤ ≤
BRXR =βReR,
nj j nj
T (27)
XR = xR ,...,xR ,...,xR ,...,xR ,
j (cid:16) jnj jnj+1 jnp−1 jnp(cid:17)
Denoting
ZR = BR −1eR, ZL = BL −1eL, (28)
ni (cid:0) ni(cid:1) ni−1 (cid:16) ni−1(cid:17)
we obtain the general formula for computing (X ) , i=1,...,p:
j ni
βL ZL , i>j
 nj(cid:16) nj(cid:17)ni
(X ) = (29)
j ni 
βR ZR , i j
 nj(cid:16) nj(cid:17)ni ≤
Substituting (29) into (23), yields (19).
The theorem has been proved.
Remark1. Since the vectors ZR,L do not depend on the right-hand side of (1), therefore, they may be
ni
defined once for the whole series of problems.
Remark2. If A=AT, then according to Lemma 2 the quantities βR,L may be defined as follows
ni
ni−1
βL = (F) G ,
ni j ni−1 j
j=Xni−1 (cid:0) (cid:1) (30)
ni−1
βR = (F) (G ) ,
ni j ni j
j=Xni−1
AGk =ek . (31)
m m
Based on Theorem 2, we formulate the algorithm for computing several components from the solu-
tion vector for the series of tridiagonal equations (1).
Algorithm 3.For computing M different components of the solution vector (Xn) , n = 1,...,N;
km
m=1,...,M from the series of equations (1), follow:
7
1. Preliminary step. (Performed once for the whole series of problems).
1.1 Find Gk , m=1,...,M, from the solution of equation (31).
m
1.2 Find ZR,L, m=1,...,M from (21).
k
m
2. Step of obtaining solutions. (Performed for each right-hand side F F , n=1,...,N.)
n
2.1 Determine βR,L, m=1,...,M according to (30) .
km
2.2 Determine (Xn) , m=1,...,M according to (19).
km
Elementary counting of arithmetic operations at the preliminary step of Algorithm 3 shows that its
realization by formulas (21) and (31) requires 24NP operations. For calculating components (18) by
≈
formulas (19),(30) requires 12N +M2 operations.
An important property of Algorithm 3 is that at each of four steps of the algorithm, calculations for
different k are independent. Therefore, the number of arithmetic operations per PE is 24NM/p for the
m
first step and (12N +M2)/p for the second, respectively.
Let us analyze the efficiency of Algorithm 1,3 without regard to communication interactions. As the
{ }
criterion, we will enter the speedup
S =T /T ,
1 p
whereT isthenumberofoperationsforsolvingoneproblemfromseries(1)byasequentialsweepalgorithm,
1
and T , by Algorithm 1,3 . Assuming p = M, T = 8N, where N is the number of unknowns and
p 1
{ }
T = 12N +2M2 /p, we have 4
p
(cid:0) (cid:1)
8Np
S = (32)
12N +p2
From (32) it follows that the speedup value increases monotonically as the number of PEs grows, and
then starting from some p>p decreases monotonically to zero.
0
Evidently, the minimal time of problem solution is achieved for the number of PEs
8Np
p =max =√6N,
0 p (cid:18)12N +p2(cid:19)
√6N
S = .
max
3
Thus, the efficiency of parallel Algorithm 1,3 for the maximum possible speedup is 30%. The
{ } ≈
remaining 70% computations fall on ”additional” operations for maintaining parallelism. From (30) and
(19) it follows that the volume of these additional computations has order O(p2), where p is the number of
PEs.
For comparison, that difficulty is also characteristic of the algorithm proposed in [10, 20] , where for
computingoffirst Un andlast Un (inourdesignation)elementsitisnecessarytosolvethetridiagonal
{ i} { i}
systemofequationswiththenumberofunknownsequaltothenumberofPEs. Sincetheauthorsproposeto
calculatethesolutionbymeansofasequentialsweepalgorithmversion,thenumberofadditionaloperations
will be of order O(p), but contrary to (19), parallel computing is not allowed.
Thus, in Algorithm 1,3 as well as in the algorithm [10, 20], the time of computing first Un and
{ } { i}
last Un elements depends linearly on the number of PEs.
{ i}
For increasing the efficiency of Algorithm 1,3 ,we will task to reduce the number of arithmetic opera-
{ }
tions in realizing formula.
4Theobtainedestimate(32)isconditionalandrepresentsratherthequalitativebehaviorofthespeedupdependenceonthe
numberofunknowns andPEs.
8
3.4. Parallel dichotomy algorithm for solving tridiagonal SLAEs.
It is required to calculate the components of the vector of solution defined in (18), it is assumed that
p=2p0 1 n, p0 >0.
− ≤
Let us enter into the consideration the sets
i−1
Ωi =n(X)nj(cid:12)(cid:12)(X)nj ∈Ω, j =2⌊log2p⌋+1−ik, k =1,...,2i−1o\j[=1Ωj, (33)
(cid:12)  
where i=1,..., log (p) +1.
⌊ 2 ⌋
It is evident that
⌊log2(p)⌋+1
Ω= Ω , Ω Ω = , i=j
i i j
{∅} 6
i[=1 \
Theorem 3. Let the components of the solution vector from the set Ω , j 1 and the quanities
j
≥
βR,L, ZR,L,G , i=1,...,p, j =1,...,p 1bedetermined. Thenforall(X) Ω , j <m log (p) +1,
ni ni nj − ni ∈ m ≤⌊ 2 ⌋
the following identity holds true
i k2−1
(X) = βR ZR + βL ZL +δ +δ , (34)
ni j=Xk1+1 nj(cid:16) nj(cid:17)ni j=Xi+1 nj(cid:16) nj(cid:17)ni k1 k2
0, k =0
1
δ =
k1
 (X)k1 ZRk1 ni −(Gk1+1)k1(F)k1 ZRk1+1 ni, k1 >0
(cid:0) (cid:1) (cid:0) (cid:1)

0, k =p+1
2
δ =
k2
 (X) ZL , k <p+1
k2 k2 ni 2
(cid:0) (cid:1)
where k and k are defined as follows
1 2
k = min (n n ), k = min (n n )
1 i t 2 t i
t,t<k,(X)nt∈(ΩjS{X0}) − t,t>k,(X)nt∈(ΩjS{Xp+1}) −
Proof. Validity ofthe theoremfollowsfromthe factthatthe knowncomponentsfromthe setΩ partition
j
the initial system according to Lemma 1 into independent subsystems and the solution of each subsystem
canberepresentedasthe sumofgeneralsolutionofahomogeneousequationandapartialnonhomogeneous
equation [21].
Based on Theorem 3, we formulate the efficient parallel algorithm for computing separate components
from the solution vector
Algorithm 4. Dichotomy algorithm. Calculation of M different components of the solution
vector (Xn) , n=1,...,N, m=1,...,M from series of equations (1), requires:
km
1. Preliminary step. (Performed once for the whole series of problems).
1.1 Find Gk , m=1,...,M, from the solution of (31).
m
1.2 Find ZR,L, m=1,...,M from (21).
k
m
2. Step of obtaining solutions. (Performed for each right-hand side F , n=1,...,N.)
n
2.1 Find βR,L, m=1,...,M according to (30) .
km
2.2 Calculate in ascending order of index i = 1,..., log (p) +1 the components of the solution
⌊ 2 ⌋
vector (Xn)km,∈Ωi, using (34).
Remark3. Elements belonging to the same set Ω , can be calculated independently.
i
9
LetusanalyzetheissueofcomputationalstabilityofAlgorithm4. WewillsaythatAlgorithm4isstable
if for (34) for all j
ZR,L 1, (35)
(cid:13) nj (cid:13)C ≤
(cid:13) (cid:13)
where (cid:13) (cid:13)
X =max (X)
k kC i {| i|}
Let us formulate stability criterion of Algorithm 4.
Theorem 4. Let the matrix A have the diagonal dominance [2]
b a + c , i=2,...,N 1, (36)
i i i
| |≥| | | | −
b c , b a , (37)
1 1 N N
| |≥| | | |≥| |
and at least in one of inequalities (36) or (37) , strict inequality holds, then Algorithm 4 is stable.
Proof. IfthematrixAhasthediagonaldominance,thenobviouslythematricesBR,L alsohaveadiagonal
k
dominance. Following the sweep algorithm [7], the solution of system BRZR =eR may be written as
k k
zR = i+1α , i=1,...,k 1,
i i=k i −
Q
c
i
α = − , i=2,...,k,
i
b +a α
i i i−1
α = c /a .
1 1 1
−
From conditions (36),(37)follows the inequality α 1 [2], from where the following estimate takes
i
| | ≤
place
i+1
zR = α 1
(cid:12)(cid:12) i (cid:12)(cid:12) (cid:12)(cid:12)(cid:12)(cid:12)iY=k i(cid:12)(cid:12)(cid:12)(cid:12)≤
Wecansimilarlyshowthat zL 1. Thus,theprese(cid:12)nceof(cid:12)diagonaldominanceentailsstabilityofAlgorithm
i ≤
4. (cid:12) (cid:12)
(cid:12) (cid:12)
The theorem has been proved
Remark 5. Since for calculating all elements from the set Ω to perform O(p) arithmetic operations
i
regardless of the index i, for realizing Algorithm 4, computing p components from the solution vector
requires O(plog p) operations.
2
Comparingthe dependence ofthe speedup onthe computing time (Fig. 2)for the dichotomy Algorithm
1,3 andalgorithm[10,20],weconcludethatthedichotomyalgorithmefficiencyforfewPEsiscomparable
{ }
withthatofAlgorithms 1,3 and[10,20]. ForagreatnumberofPEs,Algorithms 1,3 and[10,20]possess
{ } { }
a nearly zero speedup, whereas the dichotomy algorithm efficiency remains quite high.
10