стр. 16 |

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 |