Developer Reference for Intel® oneAPI Math Kernel Library - C

Developer Reference

  • 2021.1
  • 12/04/2020
  • Public Content
Contents

Direct Method

For solvers that use the direct method, the basic technique employed in finding the solution of the system
Ax
 =
b
 is to first factor
A
 into triangular matrices. That is, find a lower triangular matrix
L
 and an upper triangular matrix
U
, such that
A
=
LU
. Having obtained such a factorization (usually referred to as an
LU
 decomposition or
LU
 factorization), the solution to the original problem can be rewritten as follows.
Ax
 =
b
Equation
LUx
 =
b
Equation
L
(
Ux
) =
b
This leads to the following two-step process for finding the solution to the original system of equations:
  1. Solve the systems of equations
    Ly
     =
    b
    .
  2. Solve the system
    Ux
     =
    y
    .
Solving the systems
Ly
 =
b
 and
Ux
 =
y
 is referred to as a forward solve and a backward solve, respectively.
If a symmetric matrix
A
 is also positive definite, it can be shown that
A
 can be factored as
LL
T
 where
L
 is a lower triangular matrix. Similarly, a Hermitian matrix,
A
, that is positive definite can be factored as
A
 =
LL
H
. For both symmetric and Hermitian matrices, a factorization of this form is called a Cholesky factorization.
In a Cholesky factorization, the matrix
U
 in an
LU
 decomposition is either
L
T
 or
L
H
. Consequently, a solver can increase its efficiency by only storing
L
, and one-half of
A
, and not computing
U
. Therefore, users who can express their application as the solution of a system of positive definite equations will gain a significant performance improvement over using a general representation.
For matrices that are symmetric (or Hermitian) but not positive definite, there are still some significant efficiencies to be had. It can be shown that if
A
 is symmetric but not positive definite, then
A
 can be factored as
A
 =
LDL
T
, where
D
 is a diagonal matrix and
L
 is a lower unit triangular matrix. Similarly, if
A
 is Hermitian, it can be factored as
A
 =
LDL
H
. In either case, we again only need to store
L
,
D
, and half of
A
 and we need not compute
U
. However, the backward solve phases must be amended to solving
L
T
x
 =
D
-1
y
 rather than
L
T
x
 =
y
.

Fill-In and Reordering of Sparse Matrices

Two important concepts associated with the solution of sparse systems of equations are fill-in and reordering. The following example illustrates these concepts.
Consider the system of linear equation
Ax
 =
b
, where
A
 is a symmetric positive definite sparse matrix, and
A
 and
b
 are defined by the following:
Equation
A star (*) is used to represent zeros and to emphasize the sparsity of
A
. The Cholesky factorization of
A
 is:
A
 =
LL
T
, where
L
 is the following:
Equation
Notice that even though the matrix
A
 is relatively sparse, the lower triangular matrix
L
 has no zeros below the diagonal. If we computed
L
 and then used it for the forward and backward solve phase, we would do as much computation as if
A
 had been dense.
The situation of
L
 having non-zeros in places where
A
 has zeros is referred to as fill-in. Computationally, it would be more efficient if a solver could exploit the non-zero structure of
A
 in such a way as to reduce the fill-in when computing
L
. By doing this, the solver would only need to compute the non-zero entries in
L
. Toward this end, consider permuting the rows and columns of
A
. As described in Matrix Fundamentals, the permutations of the rows of
A
 can be represented as a permutation matrix,
P
. The result of permuting the rows is the product of
P
 and
A
. Suppose, in the above example, we swap the first and fifth row of
A
, then swap the first and fifth columns of
A
, and call the resulting matrix
B
. Mathematically, we can express the process of permuting the rows and columns of
A
 to get
B
 as
B
 =
PAP
T
. After permuting the rows and columns of
A
, we see that
B
 is given by the following:
Equation
Since
B
 is obtained from
A
 by simply switching rows and columns, the numbers of non-zero entries in
A
 and
B
 are the same. However, when we find the Cholesky factorization,
B
 =
LL
T
, we see the following:
Equation
The fill-in associated with
B
 is much smaller than the fill-in associated with
A
. Consequently, the storage and computation time needed to factor
B
 is much smaller than to factor
A
. Based on this, we see that an efficient sparse solver needs to find permutation
P
 of the matrix
A
, which minimizes the fill-in for factoring
B
 =
PAP
T
, and then use the factorization of
B
 to solve the original system of equations.
Although the above example is based on a symmetric positive definite matrix and a Cholesky decomposition, the same approach works for a general
LU
 decomposition. Specifically, let
P
 be a permutation matrix,
B
 =
PAP
T
 and suppose that
B
 can be factored as
B
 =
LU
. Then 
Ax = b
Equation    
PA
(
P
-1
P
)
x
 =
Pb
Equation    
PA
(
P
T
P
)
x
 =
Pb
Equation    
(
PAP
T
)(
P
x
) =
Pb
Equation    
B
(
Px
) =
Pb
Equation    
LU
(
Px
) =
Pb
It follows that if we obtain an
LU
 factorization for
B
, we can solve the original system of equations by a three step process:
  1. Solve
    Ly
     =
    Pb
    .
  2. Solve
    Uz
     =
    y
    .
  3. Set
    x
     =
    P
    T
    z
    .
If we apply this three-step process to the current example, we first need to perform the forward solve of the systems of equation
Ly
 =
Pb
:
Equation
This gives:
Equation
The second step is to perform the backward solve,
Uz
 =
y
. Or, in this case, since a Cholesky factorization is used,
L
T
z
 =
y
.
Equation
This gives
Equation
The third and final step is to set
x
 =
P
T
z
. This gives
Equation

Product and Performance Information

1

Performance varies by use, configuration and other factors. Learn more at www.Intel.com/PerformanceIndex.