SPARSE DIRECT METHODS Direct Sparse Matrix Methods

• Building blocks for sparse direct solvers

Problem addressed: Linear systems
• SPD case. Sparse Column Cholesky/
Ax = b
• Elimination Trees - Symbolic factorization

ä We will consider mostly Cholesky –

ä We will consider some implementation details and tricks used to
develop efficient solvers
Basic principles:

• Separate computation of structure from rest [symbolic factoriza-

tion]
• Do as much work as possible statically
• Take advantage of clique formation (supernodes, mass-elimination).
8-2 Davis: Chap. 4 – Direct

Sparse Column Cholesky The four essential stages of a solve

1. Reordering: A −→ A := P AP T
For j = 1, . . . , n Do:
l(j : n, j) = a(j : n, j) T
ä Preprocessing: uses graph [Min. deg, AMD, Nested Dissection]
L1
For k = 1, . . . , j − 1 Do:
L1 2. Symbolic Factorization: Build static data structure.
// cmod(k,j):
lj:n,j := lj:n,j − lj,k ∗ lj:n,k
ä Exploits ’elimination tree’, uses graph only.
EndDo
// cdivp (j) [Scale] ä Also: ’supernodes’
A
lj,j = lj,j
lj+1:n,j := lj+1:n,j /ljj
L 2 3. Numerical Factorization: Actual factorization A = LLT
EndDo ä Pattern of L is known. Uses static data structure. Exploits
supernodes (blas3)
4. Triangular solves: Solve Ly = b then LT x = y

8-3 Davis: Chap. 4 – Direct 8-4 Davis: Chap. 4 – Direct

 The notion of elimination tree

ä Elimination trees are useful in many different ways [theory, sym-

bolic factorization, etc..]
ä For a matrix whose graph is a tree, parent of column j < n is
defined by
ELIMINATION TREES
P arent(j) = i, where aij 6= 0 and i >j

ä For a general matrix matrix, consider A = LLT , and GF =

‘filled’ graph = graph of L + LT . Then
P arent(j) = min(i) s.t. aij 6= 0 and i>j

ä Defines a tree rooted at column n (Elimintion tree).

8-6 Davis: Chap. 4 – Direct

Example: Original matrix and Graph Filled matrix+graph

   
1 ? ? ? 1 5 1 ? ? ? 1 5
? 2 ?  ? 2 
   ? 
 ? 


3 ? 
 2 3  3 ? ? 3
 4 ? ?    2
   4 ? ? 
? ? 5 ? 8  
  ? ? 5 ? 8
 ? ? 6  6  
   ? ? 6  6
 ? 7 ?  
 ? 7 ?
? ? ? ? 8 4 7 ? ? ? ? 8
4 7

8-7 Davis: Chap. 4 – Direct 8-8 Davis: Chap. 4 – Direct

 Corresponding Elimination Tree Where does the elimination tree come from?

1 5 8 ä Answer in the form of an excercise.

1
7 1 5
3 Consider the elimination steps for the 2 3
2 1
previous example. A directed edge 5
−→ 6 1 3
8 means a row (column) modification. It 2 5
1
6 5 4 shows the task dependencies. There are 2 3
8
2 unnecessary dependencies. For example: 6
3 6
1 → 5 can be removed because it is 6
7
4 7 1 4
subsumed by the path 1 → 2 → 5.
ä Parent(i) = ’first nonzero entry in L(i+1:n,i)’ 4 7
4
ä Parent(i) = min {j > i | j ∈ AdjGF (i)} To do: Remove all the redundant dependencies.. What is the
result?

8-9 Davis: Chap. 4 – Direct 8-10 Davis: Chap. 4 – Direct

Facts about elimination trees Elim. tree depends on ordering (Not just the graph)

ä Elimination Tree defines dependencies between columns. Example: 3 × 3 grid for 5-point stencil [natural ordering]

ä The root of a subtree cannot be used as pivot before any of its 9


descendents is processed. 
8

ä Elimination tree depends on ordering; ++ + 
7 8 9 +++ + 7


ä Can be used to define ‘parallel’ tasks. ++ + 

4 5 6 + ++ + 6

+ +++ + 
ä For parallelism: flat and wide trees → good; thin and tall (e.g. + ++ 5

of tridiagonal systems) → Bad. 1 2 3 + ++ 
+ +++ 4

ä For parallel executions, Nested Dissection gives better trees than + ++ 

Minimun Degree ordering. 3



2


1


8-11 Davis: Chap. 4 – Direct 8-12 Davis: Chap. 4 – Direct

ä Same example with nested dissection ordering Properties

ä The elimination tree is a spanning tree of the filled graph [a tree

containing all vertices] - obtained by removing edges.

+ + + 9
 k
+ ++ 
1 7 2 + + + ä If lik 6= 0 then i is an
8

parent(k)
+ + + ancestor of k in the tree

2
8 + ++ + parent (k)
5 6 ,
7 - 1 In the previous exam- 3
+ + + + 
,
,
, @
@
@ ple: follow the creation of parent (k)
9 ++ ++ 
5 
6
3 4 +++++  A
S the fill-in (6,8).
   
 A
S
 A
S i
++ ++ 1
  
3 2 4

In particular: if aik 6= 0, k < i then i k

ä Consequence: no fill-in between branches of the same subtree

8-13 Davis: Chap. 4 – Direct 8-14 Davis: Chap. 4 – Direct

Elimination trees and the pattern of L 8

ä It is easy to determine the sparsity pattern of L because the In theory: To construct the 7
pattern of a given column is “inherited” by the ancestors in the tree. pattern of L, go up the
nz(L:5) := nz(L:5)
tree and accumulate the
k patterns of the columns. U nz(L:2) 6

Initially L has the same U nz(L:3)

parent(k) 5 4
Theorem: For i > j, lij 6= pattern as T RIL(A).
2
0 iff j is an ancestor of some parent (k) 2
j 3 3
k ∈ AdjA(i) in the elimina- parent (k)
1
tion tree.
i ä However: Let us assume tree is not available ahead of time
j
ä Solution: Parents can be obtained dynamically as the pattern is
In other words: being built.
∃k ∈ AdjA(i)s.t. ä This is the basis of symbolic factorization.
lij 6= 0, i > j iff
j k
8-16 Davis: Chap. 4 – Direct
Notation : ALGORITHM : 1 Symbolic factorization
ä nz(X) is the pattern of X (matrix or column, or row). A set 1. Set: nz(L) = tril(nz(A)),
of pairs (i, j) 2. Set: list(j) = ∅, j = 1, · · · , n
ä tril(X) = Lower triangular part of pattern [matlab notation] 3. For j = 1 : n
{(i, j) ∈ X |i > j} 4. for k ∈ list(j) do
5. nz(L:,j ) := nz(L:,j ) ∪ nz(L:,k )
ä Idea: dynamically create the list of nodes needed to update L:,j . 6. end
7. p = min{i > j | Li,j 6= 0}
8. list(p) := list(p) ∪ {j}
9. End

8-17 Davis: Chap. 4 – Direct 8-18 Davis: Chap. 4 – Direct

5 List= empty List= empty

1 L={2,6,8} , p = 2
L = {2,5,8}, p = 2 List={2,3}
1 5
3 1 5 List={2,3}
2
Example: Consider List={1} 3
8 2
the earlier example: 6 List = {1} 3 L={5,6,8} List=empty
2 p=5 8 L={5,8},p=5
L={5,6,8} List=empty 6 List={4}
4 7 p=5 8 L={5,8}
6 p={5} 4 7
List=empty List= empty
List=empty
L={2,5,8},p=2 L: = {2,5,8} p= 2 L={6,7}
1 5 1 5 List={2} 4 7 p=6

3 List={1} 2 3
List={1} 2 L={5,6,8}
p=5
8 8
6 6

4 7 4 7
 List= empty List= empty
L={2,5,8} , p = 2 List={2,3},L={6,8}, p=6 L={2,5,8} , p = 2
List={2,3},L={6,8},p=6 Multifrontal methods
1 5 5
1
List=empty List={1} List=empty ä Start with the frontal method.
List={1} 3 L:={5,8}
L={5,6,8}
2
p=5
L={5,6,8}
2 3 L:={5,8}
p=5 8
p=5 p=5 ä Recall: Finite element matrix:
6 8 P [e]
List={4,5} List={4,5}
6 A= A
L={7,8}
4 7 p=7 List={6}
List=empty
List=empty 4 7
L={6,7}
p=6 L={6,7} p = 6 A[e] = element matrix associated with element e.
ä An old idea: Execute Gaussian elimination as the elements are
being assembled
ä Dependency: variabes ↔ elements, creates an assembly tree.
ä Method is called the frontal method
ä Very popular among finite element users: saves storage

8-21 Davis: Chap. 4 – Direct 8-22 – Direct2

1 7 4
Multifrontal methods: extension to general matrices 9

8
ä Elimination tree replaces assembly tree 8
6
3 7
ä Proceed in post-order traversal of elimination tree in order not 3 6
2 5
to violate task dependencies. 1 2 4 5
9
ä When a node is eliminated an update matrix is created.
ä This matrix is passed to the parent which adds it to its frontal 1 3 7 2 3 9
matrix. 1 U1 = U2 =

2
3 7
ä Requires a stack of pending update matrices 3 9
3

3
3 7

3
ä Update matrices popped out as they are needed
7

9
Frontal Update Frontal Update
ä Often implemented with nested dissection-type ordering Matrix Matrix Matrix Matrix
ä More complex than a left-looking algorithm

8-23 – Direct2
 9 7 4
Eliminating nodes 1 and 2: What happens on matrix
8
7 8
6  
3
3 6 1 ? ?

4 5 5  2 ? ?
 
9  ? ? 3 ?  ← U1(3, :) ← U2(3, :)
 
 4 ? ? 
 
A3 + U 1 + U2  5 ? ?
 
3 7 8 9  ? ? 6 ? 
 
U1 U3  ? ? 7 ?  ← U (7, :)
3
  1
7 8 9  
7

 ? ? ? 8 ?

7 8 9
? ? ? 9 ← U2(9, :)
8
9

U2
Frontal Update
Matrix Matrix

8-25 – Direct2 8-26 – Direct2

Supernodes

ä Contiguous columns tend to inherit the pattern of the columns

from they are updated → Many columns with same sparsity pattern.
Supernode = a set of contiguous columns in the Cholesky factor L
that have the same sparsity pattern.

ä The set {j, j + 1, ..., j + s} is a supernode if

S
N Z(L∗,k ) = N Z(L∗,k+1) {k + 1} j ≤ k<j + s

where N Z(L∗,k ) is nonzero set of column k of L.

ä Other terms used: Mass elimination, indistinguishible nodes,
active variables in front, subscript compression,...
ä Gain in performance due to savings in Gather-Scatter operations.
8-27 – Direct2