Robust Task-Parallel Solution of the Triangular Sylvester Equation

The Bartels-Stewart algorithm is a standard approach to solving the dense Sylvester equation. It reduces the problem to the solution of the triangular Sylvester equation. The triangular Sylvester equation is solved with a variant of backward substitution. Backward substitution is prone to overflow. Overflow can be avoided by dynamic scaling of the solution matrix. An algorithm which prevents overflow is said to be robust. The standard library LAPACK contains the robust scalar sequential solver dtrsyl. This paper derives a robust, level-3 BLAS-based task-parallel solver. By adding overflow protection, our robust solver closes the gap between problems solvable by LAPACK and problems solvable by existing non-robust task-parallel solvers. We demonstrate that our robust solver achieves a similar performance as non-robust solvers.

💡 Research Summary

The paper addresses a long‑standing gap between the robustness of LAPACK’s scalar Sylvester solver (dtrsyl) and the performance of modern task‑parallel implementations for the triangular Sylvester equation that arises in the Bartels‑Stewart reduction of the dense Sylvester problem AX + XB = C. While dtrsyl prevents overflow by dynamically scaling the right‑hand side with a global factor α∈(0,1], it relies on level‑2 BLAS operations and therefore cannot exploit the full parallel potential of multicore architectures. Non‑robust task‑parallel solvers (e.g., FLASH‑Sylv) achieve high performance but may overflow on ill‑conditioned or large‑scale problems because they lack any scaling mechanism.

The authors propose a two‑stage solution. The first stage, called RobustSyl, augments the classic backward‑substitution algorithm with two protective primitives introduced in earlier work: ProtectUpdate and ProtectDivision. ProtectUpdate computes a scaling factor ζ based on the ∞‑norms of the matrices involved in an update (‖C‖∞, ‖A‖∞, ‖Y‖∞) so that the updated quantity never exceeds the overflow threshold Ω. ProtectDivision guards the Gaussian‑elimination steps required to solve the small (1×1 or 2×2) Sylvester sub‑problems that appear on the diagonal. Each diagonal block is solved as a tiny linear system with complete pivoting, and a local scaling β is applied to the block’s right‑hand side. The global scaling α is updated multiplicatively as α←α·β·ζ₁·ζ₂ after each block, guaranteeing that every intermediate result stays within the representable range.

The second stage transforms RobustSyl into a tiled, level‑3 BLAS‑centric algorithm, referred to as RobustTile (implemented as drsylv). The matrices ˜A and ˜B are partitioned into M×M and N×N tiles, respectively, such that 2×2 diagonal blocks are not split. Each tile of the solution Y carries its own local scaling factor αₖℓ, forming an “augmented tile” ⟨αₖℓ, Yₖℓ⟩. Tile updates are performed by the RobustUpdate routine, which receives three scaling factors (α, β, γ) from the operands and computes a common factor η = min{α,β,γ}. It then calls ProtectUpdate to obtain an additional factor ζ, yielding a final scaling δ = η·ζ. The update itself is expressed as a matrix‑matrix multiplication:  D = (δ/γ) C −