<divclass="textblock"><p>We implement a parallel blocked Cholesky factorization using a static <aclass="el" href="classtf_1_1Taskflow.html" title="class to create a taskflow object">tf::Taskflow</a> 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.</p>
<h1><aclass="anchor" id="CholeskyMotivation"></a>
Motivation</h1>
<p>Many scientific and engineering applications need to solve a linear system <code>Ax=b</code> where <code>A</code> 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 <code>A</code> into <code>A=L*L^T</code> (Cholesky decomposition), where <code>L</code> is a lower-triangular matrix, and then solve two triangular systems. The factorization dominates the cost, so making it fast is critical.</p>
<p>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:</p>
<divclass="fragment"><divclass="line">Tile layout (lower triangle only, since A is symmetric):</div>
<divclass="line"></div>
<divclass="line"> col 0 col 1</div>
<divclass="line">row 0 [ A00 ]</div>
<divclass="line">row 1 [ A10 ] [ A11 ]</div>
</div><!-- fragment --><p>The goal is to compute the Cholesky factor <code>L</code> in-place, one diagonal tile at a time. Here, each elimination step <code>k</code> corresponds to factorizing the <code>k-th</code> diagonal tile and using it to update all tiles below and to the right of it in the lower triangle.</p>
<p><b>Elimination step <code>k=0</code>:</b> The first step is to factorize the top-left diagonal tile in-place.</p>
</div><!-- fragment --><p>Once <code>potrf(A00)</code> completes, tile <code>A10</code> can be updated independently using the freshly computed <code>L00:</code></p>
</div><!-- fragment --><p>With <code>L10</code> in hand, the trailing submatrix must be prepared for the next elimination step. The diagonal tile <code>A11</code> is updated by subtracting the outer product of <code>L10</code> from itself:</p>
</div><!-- fragment --><p>There are no off-diagonal tiles below step <code>k=0</code> in this 2x2 example, so there is no <code>gemm</code> task here.</p>
<p><b>Elimination step <code>k=1</code>:</b> Now that <code>syrk</code> has updated <code>A11</code>, the second diagonal tile can be factorized.</p>
</div><!-- fragment --><p>There are no further sub-diagonal tiles in a 2x2 grid, so factorization is complete. The lower triangle now holds <code>L:</code></p>
</div><!-- fragment --><dlclass="section user"><dt>Why this creates a DAG, not a serial chain?</dt><dd></dd></dl>
<p>Notice that <code>syrk(A11, L10)</code> and <code>potrf(A11)</code> are strictly ordered (<code>syrk</code> must finish before <code>potrf</code> can begin), but in a larger tiling with more steps, many <code>trsm</code>, <code>syrk</code>, and <code>gemm</code> 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 <code>potrf</code> + T*(T-1) <code>trsm</code> + T*(T-1) <code>syrk</code> + T*(T-1)*(T-2)/6 <code>gemm</code> tasks, and the graph widens significantly as T grows. As a consequence, a static <aclass="el" href="classtf_1_1Taskflow.html" title="class to create a taskflow object">tf::Taskflow</a> 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.</p>
<h1><aclass="anchor" id="CholeskyAlgorithm"></a>
The Blocked Algorithm</h1>
<p>Given an SPD matrix <code>A</code> of size <code>NxN</code> tiled into <code>TxT</code> blocks of size <code>bxb</code> (so <code>N=T*b</code>), the blocked Cholesky algorithm sweeps down the diagonal one elimination step at a time. At step <code>k</code> it performs four kinds of operations on tiles, each of which becomes one task in the graph:</p>
<ul>
<li><code>potrf(k, k)</code>: Cholesky factorization of the diagonal tile <code>(k, k)</code>. This is the only serial bottleneck of each step and must complete before any other work at step <code>k</code> can begin.</li>
<li><code>trsm(i, k)</code>: Triangular solve that updates the sub-diagonal tile <code>(i, k)</code> using the just-factored diagonal tile <code>(k, k)</code>, for each <code>i>k</code>. All T-k-1 <code>trsm</code> tasks within the same step are independent of each other and can run in parallel.</li>
<li><code>syrk(i, k)</code>: Symmetric rank-b update that subtracts the outer product of tile <code>(i, k)</code> from the diagonal tile <code>(i, i)</code>, for each <code>i>k</code>. Prepares the diagonal tile for its own future <code>potrf</code>.</li>
<li><code>gemm(i, j, k)</code>: General rank-b update that subtracts the product of tiles <code>(i, k)</code> and <code>(j, k)</code> from the off-diagonal tile <code>(i, j)</code>, for each <code>i>j>k</code>. Prepares the off-diagonal tile for its own future <code>trsm</code>.</li>
</ul>
<p>The algorithm repeats for <code>k=0, 1,..., T-1</code>.</p>
<h1><aclass="anchor" id="CholeskyTaskGraph"></a>
Task Graph Structure</h1>
<p>The dependency graph for a 4-panel blocked Cholesky is shown below. Each node represents one tile-level task. Nodes are coloured by type: <code>potrf</code> (peach), <code>trsm</code> (blue), <code>syrk</code> (green), and <code>gemm</code> (yellow). Edges represent data dependencies: a task may not begin until all its predecessors have written their output tiles.</p>
<divclass="dotgraph">
<iframescrolling="no" frameborder="0" src="dot_cholesky.svg" width="675" height="520"><p><b>This browser is not able to show SVG: try Firefox, Chrome, Safari, or Opera instead.</b></p></iframe></div>
<p>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: <code>potrf</code> (peach), <code>trsm</code> (blue), <code>syrk</code> (green), and <code>gemm</code> (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.</p>
<p>We construct the full task graph before execution. Each tile is identified by its <code>(row, column)</code> 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:</p>
<divclass="line"><spanclass="keywordtype">void</span> syrk (Tile& A, Tile& C); <spanclass="comment">// C <- C - A * A^T</span></div>
<divclass="line"><spanclass="keywordtype">void</span> gemm (Tile& A, Tile& B, Tile& C); <spanclass="comment">// C <- C - A * B^T</span></div>
<divclass="line"><spanclass="keyword">const</span><spanclass="keywordtype">int</span> T = 4; <spanclass="comment">// number of tile blocks along each dimension</span></div>
<divclass="line"><spanclass="keyword">const</span><spanclass="keywordtype">int</span> b = 64; <spanclass="comment">// tile size (b x b elements)</span></div>
<divclass="ttc" id="aclasstf_1_1Executor_html"><divclass="ttname"><ahref="classtf_1_1Executor.html">tf::Executor</a></div><divclass="ttdoc">class to create an executor</div><divclass="ttdef"><b>Definition</b> executor.hpp:62</div></div>
<divclass="ttc" id="aclasstf_1_1Executor_html_a519777f5783981d534e9e53b99712069"><divclass="ttname"><ahref="classtf_1_1Executor.html#a519777f5783981d534e9e53b99712069">tf::Executor::run</a></div><divclass="ttdeci">tf::Future< void > run(Taskflow &taskflow)</div><divclass="ttdoc">runs a taskflow once</div></div>
<divclass="ttc" id="aclasstf_1_1FlowBuilder_html_a4d52a7fe2814b264846a2085e931652c"><divclass="ttname"><ahref="classtf_1_1FlowBuilder.html#a4d52a7fe2814b264846a2085e931652c">tf::FlowBuilder::emplace</a></div><divclass="ttdeci">Task emplace(C &&callable)</div><divclass="ttdoc">creates a static task</div><divclass="ttdef"><b>Definition</b> flow_builder.hpp:1781</div></div>
<divclass="ttc" id="aclasstf_1_1Task_html_a8c78c453295a553c1c016e4062da8588"><divclass="ttname"><ahref="classtf_1_1Task.html#a8c78c453295a553c1c016e4062da8588">tf::Task::precede</a></div><divclass="ttdeci">Task & precede(Ts &&... tasks)</div><divclass="ttdoc">adds precedence links from this to other tasks</div><divclass="ttdef"><b>Definition</b> task.hpp:1305</div></div>
<divclass="ttc" id="aclasstf_1_1Taskflow_html"><divclass="ttname"><ahref="classtf_1_1Taskflow.html">tf::Taskflow</a></div><divclass="ttdoc">class to create a taskflow object</div><divclass="ttdef"><b>Definition</b> taskflow.hpp:64</div></div>
</div><!-- fragment --><h1><aclass="anchor" id="CholeskyDesignPoints"></a>
Design Points</h1>
<p>Several design choices here are worth noting, which also apply broadly to other blocked dense linear algebra factorizations:</p>
<ul>
<li>Static task graph construction: The entire task graph is wired before <aclass="el" href="classtf_1_1Executor.html#a519777f5783981d534e9e53b99712069" title="runs a taskflow once">tf::Executor::run</a> is called. Because the dependency structure depends only on T (the tile count), not on the matrix values, all tasks and edges are known at construction time. This makes static <aclass="el" href="classtf_1_1Taskflow.html" title="class to create a taskflow object">tf::Taskflow</a> the right tool: the scheduler sees the full graph from the start and can make globally optimal scheduling decisions, such as prioritizing the critical path down the diagonal.</li>
<li>Critical path awareness: The critical path runs through the sequence <code>potrf(0, 0)</code> -> <code>trsm(1, 0)</code> -> <code>syrk(0, 1)</code> -> <code>potrf(1, 1)</code> -> ... -> <code>potrf(T-1, T-1)</code>. 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.</li>
<li>Task granularity: Each task operates on one <code>bxb</code> tile. If <code>b</code> is too small, task scheduling overhead dominates. If <code>b</code> 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.</li>
<li>Lower-triangular tile storage: Only tiles <code>(i, j)</code> with <code>i>=j</code> are allocated and referenced, which halves memory usage and keeps the tile indexing consistent with the mathematical convention that <code>L</code> is lower triangular.</li>
</ul>
<dlclass="section note"><dt>Note</dt><dd>The task graph shown in the figure is for T=4 and contains 20 tasks and 30 dependency edges. For T tile blocks in general, the total task count is T <code>potrf</code> tasks, T*(T-1) <code>trsm</code> tasks, T*(T-1) <code>syrk</code> tasks, and T*(T-1)*(T-2)/6 <code>gemm</code> 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. </dd></dl>
</div></div><!-- contents -->
</div><!-- PageDoc -->
</div><!-- doc-content -->
<!-- HTML footer for doxygen 1.13.1-->
<!-- start footer part -->
<divid="nav-path" class="navpath"><!-- id is needed for treeview function! -->
<ul>
<liclass="navelem"><aclass="el" href="Examples.html">Learning from Examples</a></li>
<liclass="footer">
Maintained by <ahref="https://tsung-wei-huang.github.io/">Dr. Tsung-Wei Huang</a>
—
Generated by <ahref="https://www.doxygen.org/index.html"><imgclass="footer" src="doxygen.svg" width="104" height="31" alt="doxygen"/></a> 1.13.1