Coding a Cache Optimized Parallel LU Factorization

For a general \(m \times n\) matrix \(A\) the factorization computes \(P\), \(L\) and \(U\) such that

where

• \(P\) is a \(m \times m\) permutation matrix. Our implementation will represent it a permutation vector \(p\). This vector will only contain integers and encodes the permutation as follows: interchange row \(i\) with \(p_i\) for \(1\leq i \leq m\).

  The permutation matrix is determined by a pivoting strategy. Such a strategy is necessary to prevent divisions by zero or division by very small numbers. So it enhances the numerical stability of the algorithm. Pivoting strategies differ in the amount of effort you want to invest for stability. We will use partial pivoting in our factorization.

• \(L\) is a lower triangular \(m \times m\) matrix with ones on the diagonal.

  The strict lower triangular part of \(A\) gets overwritten with \(L\). So no additional matrix is needed to store \(L\).

• \(U\) is an upper triangular \(m \times n\) matrix.

  The upper triangular part of \(A\) gets overwritten with \(U\). So no additional matrix is needed to store \(U\).

Every student who took an undergraduate class on numerical linear algebra knows how to implement this element wise using 3 loops We will derive two variants that use matrix/vector operations instead of element wise operations. Hence, we can use BLAS for the heavy computations:

• Unblocked Variant lu_unblk

  This variant uses only BLAS Level 1 and Level 2 functions. So here we can not expect high performance. But we will take advantage of it. Low performance means that the CPU is waiting for work most of the time. So we can do some extra computations if it improves numerical stability.

• Blocked Variant lu_blk

  This variant uses the unblocked variant lu_unblk only on a rather small block. Using BLAS Level 3 functions the remaining matrix gets updated with the results. We we achieve high performance. Furthermore, it dominates the overall performance and leads to an immense speedup.

  Comparing the efficiency (measured in MFLOPS) you see that the blocked variant is superior for all problem sizes. For large systems of linear equations it almost reaches 8 GFLOPS. The unblock variant drops from 1 GFLOP to 500 MFLOPS when the problem size exceeds the cache size. So efficiency is increased by a factor of 14. Note that this is achieved without multi threading just by using an algorithm that allows to exploit cache hierarchies:

  

  Of course in practise not efficiency but effectiveness matters. Effectiveness means how long does it take to compute the factorisation. So less is better:

  

  Note that both variants still have a complexity of \(\mathcal{O}(N^3)\). Only the constant for the \(N^3\) differs significantly. This can be visualized by a log-scaled plot:

  

• Parallel Tile Variant lu_tile

  Once we can exploit a single CPU core we can consider parallelization. The matrix is converted to a block format. For this blocked or tiled matrix tasks for the LU factorization are done in parallel if possible.

  The following benchmarks were done on a Quad-Core Intel Xeon Woodcrest with 2.66 GHz:

  

  Using the timing measurements of lu_blk we can see that we have an effective per-core speedup:

  

Unblocked Algorithm

The idea to develop an algorithm in terms of matrix/vector operations is based on considering partitions of matrices.

Partitioning \(A\), \(L\) and \(U\)

We partition \(A \in \mathbb{R}^{m \times n}\) along the first row and first column as follows:

So \(a_{11}\) is a scalar, \(\vec{a}_{21}\) a column vector, \(\vec{a}_{12}^T\) a row vector and \(A_{22}\) a matrix block.

We first assume that the permutation matrix \(P\) is just the identity matrix, i.e. that no pivoting is needed. In this case the factorization reads \(A = LU\). We now partition \(L\) and \(U\) accordingly.

As \(L \in \mathbb{R}^{m \times m}\) is lower triangular with unit diagonal the partition reads:

Analogously we get for \(U\):

Recursive Algorithm for a Factorization \(A = LU\)

In the matrix product equation \(A=LU\) we now simply plug in the partitioned matrices and get

Comparing partition blocks leads to the equations

So except for \(L_{22}\) and \(U_{22}\) all entries of the partitioned matrices \(L\) and \(U\) can be computed directly. The blocks \(L_{22}\) and \(U_{22}\) can be computed by computing the \(LU\) factorization of the \((m-1) \times (n-1)\) matrix \(A_{22} - \vec{l}_{21} \vec{u}_{12}^T\).

Overwriting the partitions of \(A\) with partitions of \(L\) and \(U\) reads:

So obviously this requires that \(a_{11} \neq 0\) If \(a_{11} = 0\) we have to use a pivot strategy that interchanges two rows of \(A\) first. We will add this next. For the moment we derived the well-known Gauss algorithm (without row interchanges):

Recursive Algorithm with Pivoting: \(A = PLU\)

In the above algorithm we have to avoid a division by zero. For numerical stability it is favorable that in a floating point division x/y the absolute value of y is larger than the absolute value of x.

So the idea to achieve both is simple and known as partial pivoting. Before each step we first determine the component in \(\vec{a}_{21}\) with maximal absolute value. We interchange then complete rows of the original matrix \(A\) such that \(a_{11}\) always contains the maximal absolute value. Note that is important the rows of the original matrix get interchanged.

In order to remember what rows of the original matrix where interchanged we use a pivot vector \(\vec{p}\) that initially contains values \(1, \dots, m\). If row \(i\) gets interchanged with row \(j\) we set \(p_i = j\).

Iterative Algorithm with Pivoting: \(A = PLU\)

Of course we prefer an iterative algorithm for an implementation. This is simply achieved by considering the partitioning within the original matrix:

So in the \(i\)-th step we are working on the sub-matrix \(\tilde{A}\). In Matlab notation this sub-matrix is related to \(A\) by \(\tilde{A} = A(i:m, i:n)\). Parts of \(\tilde{A}\) are therefore:

This gives the iterative algorithm for computing the \(LU\) factorization with partial pivoting:

The algorithm aborts if in the \(i\)-th step no non-zero element can be found. So yes, we continue even if \(A_{i,i}\) is almost zero. Numerical stability can be optimized if scaling is done by division if \(A_{i,i}\) is almost zero.

Improving Numerical Stability

Floating point divisions are more expensive (i.e. need more CPU cycles) than floating point multiplications. That makes it favorable to compute first the reciprocal \(\alpha = \frac{1}{A_{i,i}}\) and then the scaling \(A_{i+1:m, i} \leftarrow \alpha A_{i+1:m, i}\).

If the absolute value of \(A_{i,i}\) is small the division \(1/A_{i,i}\) can cause an overflow. Even without an overflow scaling by multiplying with the reciprocal can cause a considerable loss of precision. Because of the partial pivoting the element wise divisions \(A_{j,i}/A_{i,i}\) with \(i < j \leq m\) still give results within machine accuracy. Using the FPU for the division makes the division considerably more expensive than multiplications using SSE or AVX. However, using the FPU the division gets internally done with 80 bits floating point numbers instead of 64 or 32 bits.

Note that the dominant operations in the above algorithm are Vector Scaling (BLAS Level 1) and Rank 1 Update (BLAS Level 2). We don't have any BLAS Level 3 function. So this algorithm can not exploit caches. This means the CPU will wait most of the time for data from the memory. So once the CPU finally gets something to do we can do some extra computations if it improves roundoff errors. So in each step we will check whether the absolute value of \(A_{i,i}\) is smaller than a threshold. Based on that we decide to scale by \(1/A_{i,i}\) (possibly using internally SSE or AVX) or doing element wise divisions using the FPU.

Numerical stability can be further improved by using full pivoting for the \(LU\) factorization. If this is not sufficient more stable (and more expensive) factorizations like \(QR\) factorization or a singular value decomposition can be considered.

Implementation of lu_unblk using Operators

Implementation of lu_unblk using flens::blas functions

In case you find explicit function calls more expressive you can do so. For example, it makes it more obvious that only BLAS Level 1 and BLAS Level 2 functions get used. Actually we usually start an implementation with explicit function calls and only replace them with operator notations if it make code more expressive.

The previous previous code is equivalent to:

Auxiliary Function for Applying the Permutation

Applying the permutation \(P\) to a matrix \(X\) means computing \(PX\). The permutation matrix \(P\) is represented by a pivot vector \(p\) that stores row interchanges that happened during the \(LU\) factorization. In this case applying \(P\) means reversing the interchanges such that \(PLU\) gives the original matrix \(A\). This is done by the auxiliary function apply_perm:

Auxiliary Wall Time Function

For measuring the elapsed time for doing a certain task we use a simple wall-time clock. The term 'wall-time' means that for a task we want to measure the real time that elapses from start to end, i.e. including time used for system calls etc. It can be used following the pattern:

The code for ATL_walltime was taken from ATLAS:

Code for Testing

Note:

Uncomment for playing with small dimensions:

and

Compile and Run

$shell> cd flens/examples                                                                          
$shell> g++ -Wall -std=c++11 -DWITH_OPERATORS -I../.. -o lu_unblk_operators_test lu_unblk_test.cc  
$shell> g++ -Wall -std=c++11 -I../.. -o lu_unblk_flensblas_test lu_unblk_test.cc                   
$shell> ./lu_unblk_operators_test                                                                  
Time elapsed: 13.96
Checking results:
|| A-P*L*U ||_1 = 1.66341e-08
$shell> ./lu_unblk_flensblas_test                                                                  
Time elapsed: 13.85
Checking results:
|| A-P*L*U ||_1 = 1.66341e-08

Blocked Algorithm

We now derive an algorithm that takes advantage of BLAS Level 3 functions and applies the unblocked \(LU\) factorization only on a block that fits into the L2 cache.

Block-Partitioning A, L and U

We now consider a block-wise partition of \(A \in \mathbb{R}^{m \times n}\)

where block \(A_{11}\) is square, i.e. \(A_{11} \in \mathbb{R}^{\text{BS} \,\times\, \text{BS}}\), \(A_{21} \in \mathbb{R}^{(m-\text{BS}) \,\times\, \text{BS}}\), \(A_{12} \in \mathbb{R}^{\text{BS} \,\times\, (n-\text{BS})}\) and \(A_{22} \in \mathbb{R}^{(m-\text{BS}) \,\times\, (n-\text{BS})}\) for some block size \(\text{BS} \in \mathbb{N}\).

Analogously we consider for \(L\) the block-partition

where \(L_{11}\) adn \(L_{22}\) are lower triangular matrix with unit diagonal. For \(U\) the block-partition

blocks \(U_{11}\) and \(U_{22}\) are obviously upper triangular matrices.

So in case no pivoting gets applied the \(LU\) factorization is determined by

Comparing partition blocks leads to the equations

Note that because \(L_{11}\) and \(U_{11}^T\) are both lower triangular computing \(L_{11}^{-1} A_{12}\) and \(U_{11}^{-T} A_{21}^T\) can be done using a triangular solver (forward substitution).

Recursive Block Algorithm for a Factorization \(A = LU\)

So as we again want to overwrite \(A\) with \(L\) and \(U\) we end up with the algorithm

The triangular solvers (blas::sm) as well as the matrix-matrix product (blas::mm) are BLAS Level 3 functions.

Iterative Block Algorithm with Pivoting: \(A = PLU\)

In order to derive an iterative algorithm we again consider the partitioning within the original matrix:

In the \(i\)-th iteration the dimension of \(\tilde{A}\) is \(\bigl(m - (i-1)\cdot\text{BS}\bigr) \times \bigl(n - (i-1)\cdot\text{BS}\bigr)\). Therefore blocks of \(\tilde{A}\) have dimensions

This gives the iterative algorithm for computing block wise the \(A=PLU\) factorization with partial pivoting:

Note

What does updating \(P\) with \(\tilde{P}\) mean? The permutation matrix \(P\) is represented by a pivot vector \(p = (p_1, \dots, p_m)\) and \(\tilde{P}\) by a pivot vector \(\tilde{p} = (\tilde{p}_1, \dots, \tilde{p}_{\text{bs}})\). So an update means to set \(p_{i+k-1} = i - 1 + \tilde{p}_{k}\) for \(1 \leq k \leq \text{bs}\).

Auxiliary Function for Applying the Inverse Permutation

Applying the inverse permutation \(P^{-1}=P^T\) to a matrix \(X\) means computing \(P^T X\). Recall, the permutation matrix \(P\) is represented by a pivot vector \(p\) that stores row interchanges that happened during the \(LU\) factorization. So applying \(P\) means reversing the interchanges such that \(PLU\) gives the original matrix \(A\). Hence, applying \(P^{-1}\) to \(X\) means to apply the same row interchanges that were done during the \(LU\) factorization. This is done by the auxiliary function apply_perm_inv:

Implementation of lu_blk using flens::blas functions

Implementation of lu_blk using Operators

For applying triangular solvers we do not provide overloaded operators. So you might be disappointed if you a big fan of it. In that case we leave it as a simple exercise to the dear user to support for a triangular matrix L something like X = inverse(L)*B that gets evaluated as blas::sm(Left, NoTrans, 1.0, L, X) if X and B are identical and as X = B; blas::sm(Left, NoTrans, 1.0, L, X) otherwise.

Code for Testing

Compile and Run

$shell> cd flens/examples                                                                          
$shell> g++ -Wall -std=c++11 -DWITH_OPERATORS -I../.. -o lu_blk_operators_test lu_blk_test.cc      
$shell> g++ -Wall -std=c++11 -I../.. -o lu_blk_flensblas_test lu_blk_test.cc                       
$shell> ./lu_blk_operators_test                                                                    
Time elapsed: 0.95
Checking results:
|| A-P*L*U ||_1 = 3.56198e-07
$shell> ./lu_blk_flensblas_test                                                                    
Time elapsed: 0.95
Checking results:
|| A-P*L*U ||_1 = 3.56198e-07

Having an algorithm that is capable to exploit caches and BLAS it now makes also sense to compile without debug code and optimizations:

$shell> cd flens/examples                                                                          
$shell> g++ -Wall -DNDEBUG -O3 -std=c++11 -DWITH_OPERATORS -I../.. -o lu_blk_operators_test lu_blk_test.cc                                                                         
$shell> g++ -Wall -DNDEBUG -O3 -std=c++11 -I../.. -o lu_blk_flensblas_test lu_blk_test.cc          
$shell> ./lu_blk_operators_test                                                                    
Time elapsed: 0.76
Checking results:
|| A-P*L*U ||_1 = 3.56198e-07
$shell> ./lu_blk_flensblas_test                                                                    
Time elapsed: 0.76
Checking results:
|| A-P*L*U ||_1 = 3.56198e-07

Block Algorithm with Parallel BLAS

The easiest way to parallelize the \(LU\) factorization is using a multi-threaded BLAS implementation for the blocked \(LU\) factorization. In this case 'easy' means that you don't have to change function lu_blk you just link against a parallel BLAS library. So the multi-threading part is hidden in BLAS.

In the near future we will extend the ulmBLAS tutorial for parallelization...

Parallel Tile Algorithm

Having a parallel computer architecture we hope that using 2 CPU cores for the \(LU\) factorization cuts the wall-time in half, i.e. gives a speedup factor of 2. And of course, with 3 cores we hope for a speedup factor of 3, with 4 cores a speedup factor or 4 etc. This means we hope for a perfect scalability of the problem. However, this usually requires that for \(k\) CPU cores we can work on \(k\) disjoint sub-problems with perfect load balancing. The later means each sub-problem has the same complexity and the same memory bandwidth can be used for solving each of these sub-problems.

For BLAS Level 3 functions this can be achieved in certain cases. For example, in the matrix product \(C = A\cdot B\) each thread computes a different block of \(C\). However, a good speedup factor requires that the computation done by a thread compensates the induced overhead, e.g. scheduling threads or copying buffered results.

The PLASMA library tries to achieve a much higher granularity of sub-problems by considering parallelization already in the algorithms. These algorithms basically push forward ideas from blocked algorithms and are known as tiled algorithms.

Example: Tasks for Factorizing a 2x2 Tile/Block Matrix

We again consider the \(LU\) factorization of a \(2 \times 2\) blocked matrix:

In the blocked algorithm we identified a sequence of necessary tasks that were done in the certain fixed order. In this case the total sequence of tasks done were:

Obviously these tasks can not be done in an arbitrary order. And therefore some tasks can not be done in parallel. However, analyzing the dependencies will show what tasks could be done in parallel.

Task Dependencies

We illustrated the dependencies through a directed graph. Hereby \(T_1 \rightarrow T_2\) denotes that task \(T_1\) must be completed before task \(T_2\) can be started.

So once task \(A_11 = P_1 L_{11} U_{11}\) is done we could

• Compute \(A_{12} \leftarrow P_1^T A_{12}\) by one thread and

• \(A_{12} \leftarrow A_{12} U_{11}^{-1}\) by another thread.

Considering more tiles shows that more and more independent task can be identified.

Example: Tasks for Factorizing a 3x3 Tile/Block Matrix

Listing all operations needed for a \(3 \times 3\) tile matrix is already lengthy.

Gladly we will develop in this section a tool that automatically generates the task dependencies:

Example: Tasks for Factorizing a 10x10 Tile/Block Matrix

Just for fun (and because we have the tool) we also have a graph generated for a 10x10 tile matrix.

Converting Matrices to Tile Matrices

Implementation of lu_tiled

For the implementation of the tile algorithm we need a scheduler. This scheduler can be fed with tasks (i.e. functions) that have dependency constraints, e.g. task A must be completed before task B can be started. Tasks can be done in parallel if constraints permit and if enough threads are available.

A task is just a function pointer. Using C++11 lambdas makes is a convenient way of creating such functions. Furthermore it is required to define unique keys for each task.

We will first look at the implementation of the tile algorithm and discuss the design of the scheduler and key functions later:

Task Scheduler and Thread Pool

Class Scheduler provides an abstract interface a thread pool and scheduler:

Internally task are identified by an integer.

And a task is just a function pointer.

The scheduler and thread pool get initialized when the application gets started, e.g. in main(). Threads will be spawned and initialized. They will immediately wait on a condition variable until a task is ready.

For illustration purposes it is nice to identify task by a string. In our example these strings happen to be Latex code. So we are able to dump the scheduled directed graph and generate nice pictures.

Adding a task identified by keyStr can depend on a list of other tasks (pre) that must have been completed before it can be started.

In some cases it is convenient to add only some dependencies, i.e. add an arc to the directed graph.

Note that if you want to add 'pre must be done before succ gets started' it is necessary that succ has not been scheduled already. Otherwise task succ could already be running.

Again tasks pre and succ are identified by string.

Once all tasks are scheduled we wait until all tasks are completed.

If a task calls abort then no new tasks will be stared. join will wait until all running tasks are completed and returns errorCode.

The implementation uses typical C++11 features for dealing with threads, locks and condition variables. The directed graph is realized using containers from the C++ standard library:

Auxiliary Functions for Task Keys

For creating task identifiers we simply create strings that contain some Latex code describing what the task is supposed to do. For convenience we define a function such that we can use the beloved format strings from C:

Now we simply can define a bunch of functions for generating different keys:

Code for Testing

The test code first initializes the scheduler. Copies and converts the matrix into tile format and calls the implementation of the tile algorithm. Once the factorization is completed the tiled matrix is copied and converted back.

Compile and Run

$shell> cd flens/examples                                               
$shell> g++ -O3 -Wall -std=c++11 -I../.. -c lu/scheduler.cc             
$shell> g++ -O3 -Wall -std=c++11 -I../.. scheduler.o lu_tiled_test.cc   
$shell> ./a.out                                                         
total time elapsed:          0.0999999
time elapsed for 'lu_tiled': 0.0699999
time elapsed for conversion: 0.03
Checking results:
|| A-P*L*U ||_1 = 4.7791e-09

View document source

© 2011-2015 Michael Lehn