<< . .

. 16
( : 95)



. . >>




FIGURE 3.5. Reordered graphs using the direct (left) and reverse (right)
Cuthill-McKee algorithm. The label of each vertex, before reordering is per-
formed, is reported in braces



xxx x x xxx
xxxxxx xxxx
xxx x xxxx x
x xx
x xx xx
xxx
xx
xx xx x x
x xx x x
xx x
x xx xx xxxx xx
xx x xx xx x
xx
xx x xx x
xx
xxx x xx xx
xxxx x xx xx
x x xx x xx
x

FIGURE 3.6. Factors L and U after the direct (left) and reverse (right)
Cuthill-McKee reordering. In the second case, ¬ll-in is absent



3.9.2 Decomposition into Substructures
These methods have been developed in the framework of numerical ap-
proximation of partial di¬erential equations. Their basic strategy consists
of splitting the solution of the original linear system into subsystems of
smaller size which are almost independent from each other and can be
easily interpreted as a reordering technique.
We describe the methods on a special example, referring for a more
comprehensive presentation to [BSG96]. Consider the linear system Ax=b,
where A is a symmetric positive de¬nite matrix whose pattern is shown in
Figure 3.3. To help develop an intuitive understanding of the method, we
draw the graph of A in the form as in Figure 3.7.
3.9 Sparse Matrices 101

We then partition the graph of A into the two subgraphs (or substruc-
tures) identi¬ed in the ¬gure and denote by xk , k = 1, 2, the vectors of the
unknowns relative to the nodes that belong to the interior of the k-th sub-
structure. We also denote by x3 the vector of the unknowns that lie along
the interface between the two substructures. Referring to the decomposi-
T T
tion in Figure 3.7, we have x1 = (2, 3, 4, 6) , x2 = (8, 9, 10, 11, 12)
T
and x3 = (1, 5, 7) .
As a result of the decomposition of the unknowns, matrix A will be
partitioned in blocks, so that the linear system can be written in the form

11 10
12


5 1
7 8 9
substructure II
3
4

substructure I
2
6
FIGURE 3.7. Decomposition into two substructures


® ® ® 
A11 0 A13 x1 b1
°0 A23 » ° x2 » = ° b2 » ,
A22
AT AT A33 x3 b3
13 23

having reordered the unknowns and partitioned accordingly the right hand
side of the system. Suppose that A33 is decomposed into two parts, A33
and A33 , which represent the contributions to A33 of each substructure.
Similarly, let the right hand side b3 be decomposed as b3 +b3 . The original
linear system is now equivalent to the following pair

A11 A13 x1 b1
= ,
b3 + γ 3
AT A33 x3
13

A22 A23 x2 b2
=
b3 ’ γ 3
AT A33 x3
23

having denoted by γ 3 a vector that takes into account the coupling between
the substructures. A typical way of proceeding in decomposition techniques
consists of eliminating γ 3 to end up with independent systems, one for each
102 3. Direct Methods for the Solution of Linear Systems

substructure. Let us apply this strategy to the example at hand. The linear
system for the ¬rst substructure is

A11 A13 x1 b1
= . (3.62)
b3 + γ 3
AT A33 x3
13

Let us now factorize A11 as HT H11 and proceed with the reduction method
11
already described in Section 3.8.3 for block tridiagonal matrices. We obtain
the system

H11 H21 x1 c1
=
A33 ’ H21 HT b3 + γ 3 ’ H21 c1
0 x3
21

where H21 = H’T A13 and c1 = H’T b1 . The second equation of this system
11 11
yields γ 3 explicitly as

γ 3 = A33 ’ HT H21 x3 ’ b3 + HT c1 .
21 21

Substituting this equation into the system for the second substructure, one
ends up with a system only in the unknowns x2 and x3

A22 A23 x2 b2
= , (3.63)
AT A33 x3 b3
23

where A33 = A33 ’ HT H21 and b3 = b3 ’ HT c1 . Once (3.63) has been
21 21
solved, it will be possible, by backsubstitution into (3.62), to compute also
x1 .
The technique described above can be easily extended to the case of
several substructures and its e¬ciency will increase the more the substruc-
tures are mutually independent. It reproduces in nuce the so-called frontal
method (introduced by Irons [Iro70]), which is quite popular in the solution
of ¬nite element systems (for an implementation, we refer to the UMF-
PACK library [DD95]).

Remark 3.6 (The Schur complement) An approach that is dual to
the above method consists of reducing the starting system to a system
acting only on the interface unknowns x3 , passing through the assembling
of the Schur complement of matrix A, de¬ned in the 3—3 case at hand as

S = A33 ’ AT A’1 A13 ’ AT A’1 A23 .
13 11 23 22

The original problem is thus equivalent to the system

Sx3 = b3 ’ AT A’1 b1 ’ AT A’1 b2 .
13 11 23 22

This system is full (even if the matrices Aij were sparse) and can be solved
using either a direct or an iterative method, provided that a suitable pre-
conditioner is available. Once x3 has been computed, one can get x1 and
3.10 Accuracy of the Solution Achieved Using GEM 103

x2 by solving two systems of reduced size, whose matrices are A11 and A22 ,
respectively.
We also notice that if the block matrix A is symmetric and positive
de¬nite, then the linear system on the Schur complement S is no more
ill-conditioned than the original system on A, since

K2 (S) ¤ K2 (A)

(for a proof, see Lemma 3.12, [Axe94]. See also [CM94] and [QV99]).


3.9.3 Nested Dissection
This is a renumbering technique quite similar to substructuring. In practice,
it consists of repeating the decomposition process several times at each
substructure level, until the size of each single block is made su¬ciently
small. In Figure 3.8 a possible nested dissection is shown in the case of the
matrix considered in the previous section. Once the subdivision procedure
has been completed, the vertices are renumbered starting with the nodes
belonging to the latest substructuring level and moving progressively up to
the ¬rst level. In the example at hand, the new node ordering is 11, 9, 7,
6, 12, 8, 4, 2, 1, 5, 3.
This procedure is particularly e¬ective if the problem has a large size and
the substructures have few connections between them or exhibit a repetitive
pattern [Geo73].



3.10 Accuracy of the Solution Achieved Using
GEM
Let us analyze the e¬ects of rounding errors on the accuracy of the solution
yielded by GEM. Suppose that A and b are a matrix and a vector of
¬‚oating-point numbers. Denoting by L and U, respectively, the matrices
of the LU factorization induced by GEM and computed in ¬‚oating-point
arithmetic, the solution x yielded by GEM can be regarded as being the
solution (in exact arithmetic) of the perturbed system (A + δA)x = b,
where δA is a perturbation matrix such that

|δA| ¤ nu 3|A| + 5|L||U| + O(u2 ), (3.64)

where u is the roundo¬ unit and the matrix absolute value notation has
been used (see [GL89], Section 3.4.6). As a consequence, the entries of δA
will be small in size if the entries of L and U are small. Using partial
pivoting allows for bounding below 1 the module of the entries of L in such
a way that, passing to the in¬nity norm and noting that L ∞ ¤ n, the
104 3. Direct Methods for the Solution of Linear Systems


2 A
1
0000000
1111111
0000000000000000
1111111111111111
1111111111111111
0000000000000000
1111111
0000000
0000000000000000
1111111111111111
1111111111111111
0000000000000000
1 0000000
1111111
0000000000000000
1111111111111111
1111111111111111
0000000000000000
1111111
0000000
1111111111111111
0000000000000000
0000000000000000
1111111111111111
1111111
0000000
0000000000000000
1111111111111111
2 0000000000000000
1111111111111111
0000000
1111111
0000000000000000
1111111111111111
1111111111111111
0000000000000000
1111111
0000000
00000
11111
1111111111111111
0000000000000000
1111111111111111
0000000000000000
A
00000
11111
0000000000000000
1111111111111111
1111111111111111
0000000000000000
00000
11111
0000000000000000
1111111111111111
0000000000000000
1111111111111111
2 00000
11111
0000000000000000
1111111111111111
1111111111111111
0000000000000000
00000
11111
0000000000000000
1111111111111111
0000000000000000
1111111111111111
1111111111111111
0000000000000000
0000000000000000
1111111111111111
1111111111111111
0000000000000000
0000000000000000
1111111111111111
1 1111111111111111
0000000000000000
A 0000000000000000
1111111111111111
1111111111111111
0000000000000000
0000000000000000
1111111111111111
1111111111111111
0000000000000000


3 4 C 5 6B A
3 000
111 0000000
1111111
000000000000000
111111111111111
11
00
C
1111111
000000000000000
111111111111111
0000000
111
000
00
11 111111111111111
1111111
000000000000000
0000000
11
00 000000000000000
0000000
1111111
111111111111111
111
000 111111111111111
1111111
0000000
000000000000000
000
111 111111111111111
000000000000000
1111111
0000000
0000000
111111111111111
000000000000000
1111111
111
000
4 000000000000000
111111111111111
1111111
0000000
4
3 1111111
0000000
000000000000000
111111111111111
000
111 0000000
111111111111111
000000000000000
1111111
1111111
0000000
111111111111111
000000000000000
C 0000000
111111111111111
000000000000000
1111111
000
111
111111111111111
0000000
1111111
000000000000000
A 11111
00000
000000000000000
111111111111111
00000
11111
000
111
000000000000000
111111111111111
11111
00000
5 000000000000000
111111111111111
00000
11111
000
111
111111111111111
000000000000000
00000
11111
11
00
000000000000000
111111111111111
00000
11111
111111111111111
000000000000000
00000
11111
00
11

<< . .

. 16
( : 95)



. . >>