We implement a parallel blocked Cholesky factorization using a static tf::Taskflow task graph. This example shows how a classic dense linear algebra algorithm with a rich, irregular dependency structure maps naturally onto Taskflow's task graph model.
Many scientific and engineering applications need to solve a linear system Ax=b where A is symmetric positive definite (SPD): structural mechanics, Gaussian process regression, Kalman filtering, and interior-point optimization all reduce to this problem at their core. The standard approach is to first factor A into A=L*L^T (Cholesky decomposition), where L is a lower-triangular matrix, and then solve two triangular systems. The factorization dominates the cost, so making it fast is critical.
To understand why task graph parallelism is the right model here, consider a small concrete example. Take a 4x4 matrix tiled into 2x2 blocks of size 2x2 each, giving a 2x2 grid of tiles:
The goal is to compute the Cholesky factor L in-place, one diagonal tile at a time. Here, each elimination step k corresponds to factorizing the k-th diagonal tile and using it to update all tiles below and to the right of it in the lower triangle.
Elimination step k=0: The first step is to factorize the top-left diagonal tile in-place.
Once potrf(A00) completes, tile A10 can be updated independently using the freshly computed L00:
With L10 in hand, the trailing submatrix must be prepared for the next elimination step. The diagonal tile A11 is updated by subtracting the outer product of L10 from itself:
There are no off-diagonal tiles below step k=0 in this 2x2 example, so there is no gemm task here.
Elimination step k=1: Now that syrk has updated A11, the second diagonal tile can be factorized.
There are no further sub-diagonal tiles in a 2x2 grid, so factorization is complete. The lower triangle now holds L:
Notice that syrk(A11, L10) and potrf(A11) are strictly ordered (syrk must finish before potrf can begin), but in a larger tiling with more steps, many trsm, syrk, and gemm tasks across different columns have no dependency on each other at all and can run simultaneously. For a matrix tiled into T blocks, the number of tasks is T potrf + T*(T-1) trsm + T*(T-1) syrk + T*(T-1)*(T-2)/6 gemm tasks, and the graph widens significantly as T grows. As a consequence, a static tf::Taskflow captures this full dependency structure at construction time, giving the scheduler complete visibility to maximize parallelism across all elimination steps. The full DAG is determined by the tile dimensions alone, not by the matrix values, so it can be constructed once before any computation begins and executed by the scheduler without further coordination.
Given an SPD matrix A of size NxN tiled into TxT blocks of size bxb (so N=T*b), the blocked Cholesky algorithm sweeps down the diagonal one elimination step at a time. At step k it performs four kinds of operations on tiles, each of which becomes one task in the graph:
potrf(k, k): Cholesky factorization of the diagonal tile (k, k). This is the only serial bottleneck of each step and must complete before any other work at step k can begin.trsm(i, k): Triangular solve that updates the sub-diagonal tile (i, k) using the just-factored diagonal tile (k, k), for each i>k. All T-k-1 trsm tasks within the same step are independent of each other and can run in parallel.syrk(i, k): Symmetric rank-b update that subtracts the outer product of tile (i, k) from the diagonal tile (i, i), for each i>k. Prepares the diagonal tile for its own future potrf.gemm(i, j, k): General rank-b update that subtracts the product of tiles (i, k) and (j, k) from the off-diagonal tile (i, j), for each i>j>k. Prepares the off-diagonal tile for its own future trsm.The algorithm repeats for k=0, 1,..., T-1.
The dependency graph for a 4-panel blocked Cholesky is shown below. Each node represents one tile-level task. Nodes are coloured by type: potrf (peach), trsm (blue), syrk (green), and gemm (yellow). Edges represent data dependencies: a task may not begin until all its predecessors have written their output tiles.
The dependency graph for a blocked Cholesky with T=4 is shown below. Each node represents one tile-level task. Nodes are coloured by type: potrf (peach), trsm (blue), syrk (green), and gemm (yellow). Edges represent data dependencies: a task may not begin until all its predecessors have written their output tiles. The parallelism grows as factorization progresses down the diagonal and then narrows again as the matrix is exhausted.
We construct the full task graph before execution. Each tile is identified by its (row, column) indices; each task captures references to the relevant tiles and calls the corresponding BLAS routine. For clarity the code below uses placeholder functions in place of actual BLAS calls:
Several design choices here are worth noting, which also apply broadly to other blocked dense linear algebra factorizations:
potrf(0, 0) -> trsm(1, 0) -> syrk(0, 1) -> potrf(1, 1) -> ... -> potrf(T-1, T-1). Every other task is off this path and can be scheduled opportunistically. Taskflow's work-stealing scheduler naturally keeps workers busy with off-path tasks while the critical-path tasks proceed as fast as possible.bxb tile. If b is too small, task scheduling overhead dominates. If b is too large, there is not enough parallelism in the early steps. A tile size of 128 to 256 is typical for modern CPUs; the right value depends on cache size and the number of workers.(i, j) with i>=j are allocated and referenced, which halves memory usage and keeps the tile indexing consistent with the mathematical convention that L is lower triangular.potrf tasks, T*(T-1) trsm tasks, T*(T-1) syrk tasks, and T*(T-1)*(T-2)/6 gemm tasks. The degree of parallelism peaks near the middle steps. Real implementations with T=32 or more generate thousands of tasks and achieve near-linear scaling on multi-core CPUs. 
1.12.0