Robust Task-Parallel Solution of the Triangular Sylvester Equation

The robust solver drsyl solves the scaled triangular Sylvester equation AX + XB = sC in a tiled and task-parallel fashion. The matrices A and B are assumed to be upper quasi-triangular; the right-hand side matrix C and the solution matrix X are dense. The scalar s is computed alongside with the solution X such that overflow is avoided and the solution matrix X at no point during the computation contains infinities representing overflow. This property qualifies our solver as robust.

Download

The source code has migrated to Github

Prerequisites

1. A compiler that supports OpenMP. The compiler optimization level must be chosen such that associative math is disabled. The Makefile contains a suggested list of optimization flags for the Intel compiler and the GNU compiler.
2. An efficient BLAS implementation.
3. An implementation of the LAPACK routines. Only needed for the validation.

Building and executing

A Makefile for the GNU compiler linked against OpenBLAS and the Intel compiler linked against MKL is included. After unpacking, try make. The executable drsyl takes 8 input parameters. For a quick test, type ./drsyl 5000 5000 400 0.5 0.5 1 0 0. The input parameters are as follows.

• m. The matrix size of A and the row count of C and X.
• n. The matrix size of B and the column count of C and X.
• tlsz. The tile size. A good starting value is 400. For a good performance, the tile size needs to be tuned.
• cmplx-ratio-A. The proportion of 2-by-2 blocks on the diagonal of A.
• cmplx-ratio-B. The proportion of 2-by-2 blocks on the diagonal of B.
• sign. The sign of the matrix B in the Sylvester equation.
• mem-layout. The storage format of the matrices, 0=column major or 1=tile layout.
• seed. The seed for the random number generator to reproduce runs.

The solver uses task parallelism. The number of threads is controlled with OMP_NUM_THREADS.

Remarks

By default, a double-precision number is used for the scaling factor s. Sometimes the systems grow so quickly that overflow protection with a double-precision scaling factor does not suffice. Then setting -DINTSCALING during the build process activates integer scaling factors and allows for solving systems that are not solvable by a double-precision scaling factor. This change requires a complete rebuild (make clean, make).

If you want to use the triangular Sylvester equation solver sylvester.c in another code, make sure to pass appropriately partitioned matrices (see main.c for an example on how to generate such partitionings).

The matrices A and B are partitioned identically, i.e., the partitioning uses the input parameter tlsz for both matrices. The usage of distinct tile sizes should work, but has not been tested.

Other high-performance implementations

• The Fortran library RECSY contains a family of recursion-based solvers for the triangular Sylvester equation.
• The high-performance library libflame comprises two triangular Sylvester equation solvers, FLA_Sylv and FLASH_Sylv.

References and supplementary material

1. A. Schwarz and C. C. Kjelgaard Mikkelsen. Robust Task-Parallel Solution of the Triangular Sylvester Equation. In \emph{Parallel Processing and Applied Mathematics}, Lecture Notes in Computer Science, vol 12043, Springer, 2020.
2. A preprint of this paper is available on arXiv.
3. Presentation slides

Last update: 2020-08-19