cuLME - Linear Matrix Equation Solvers using CUDA

Copyright: Maximilian Behr

License: The software is licensed under under MIT. See LICENSE for details.

cuLME (Linear Matrix Equations) is a CUDA library implementing GPU solvers for the following Sylvester-like linear matrix equations:

Equation	Form
SYLV	$op_1(A) X + \mathrm{sgn}, X, op_2(B) = Y$
SYLV2	$op_1(A) X op_2(B) + \mathrm{sgn}, X = Y$
LYAP	$op_1(A) X + X, op_1(A)^H = Y$
STEIN	$op_1(A) X op_1(A)^H - X = Y$
GSYLV	$op_1(A) X op_2(B) + \mathrm{sgn}, op_1(C) X op_2(D) = Y$
GLYAP	$op_1(A) X op_1(B)^H + op_1(B) X op_1(A)^H = Y$
GSTEIN	$op_1(A) X op_1(A)^H - op_1(B) X op_1(B)^H = Y$
CSYLV	$op_1(A) R + \mathrm{sgn}_1 L, op_2(B) = E$, $\quad op_1(C) R + \mathrm{sgn}_2 L, op_2(D) = F$
CSYLV_DUAL	$op_1(A)^H R + op_1(C)^H L = E$, $\quad \mathrm{sgn}_1 R, op_2(B)^H + \mathrm{sgn}_2 L, op_2(D)^H = F$

Each equation is available in two flavors:

Triangular solvers (Tr*, Tg*): the coefficient matrices are expected in (generalized) Schur form as computed by LAPACK GEES/GGES. Real matrices may be quasi upper triangular with 1x1/2x2 diagonal blocks, complex matrices are strictly upper triangular.
General solvers (Ge*, Gg*): the coefficient matrices are general dense matrices. Since cuSOLVER provides no non-symmetric (generalized) Schur decomposition, these solvers use GPU-native iterative algorithms (matrix sign function Newton iteration and squared Smith iteration) built entirely from gemm, geam and LU factorizations.

cuLME supports real and complex, single and double precision matrices (float, double, cuComplex, cuDoubleComplex).

cuLME is purely GPU based: all matrices are device pointers in column-major storage, the workspace is queried with the *BufferSize functions and provided by the caller, no hidden allocations are performed by the solvers.

Available Functions

The functions follow the naming scheme culme_<p><Solver> with the precision prefix p in s (float), d (double), c (cuComplex) and z (cuDoubleComplex). The right-hand side d_X (d_E/d_F for the coupled equations) is overwritten by the solution. transA/transB are 'N' and 'T' for real types and 'N' and 'C' for complex types ('T' is accepted as a synonym for 'C'), sgn, sgn1, sgn2 are +1 or -1.

Double Precision Functions

/* triangular solvers, coefficients in (generalized) Schur form */
int culme_dTrSylvBufferSize(const int m, const int n, size_t *d_bufferSize, size_t *h_bufferSize);
int culme_dTrSylv(const char transA, const char transB, const double sgn, const int m, const int n, const double *d_A, const int ldA, const double *d_B, const int ldB, double *d_X, const int ldX, void *d_buffer, void *h_buffer);

int culme_dTrSylv2BufferSize(const int m, const int n, size_t *d_bufferSize, size_t *h_bufferSize);
int culme_dTrSylv2(const char transA, const char transB, const double sgn, const int m, const int n, const double *d_A, const int ldA, const double *d_B, const int ldB, double *d_X, const int ldX, void *d_buffer, void *h_buffer);

int culme_dTrLyapBufferSize(const int m, size_t *d_bufferSize, size_t *h_bufferSize);
int culme_dTrLyap(const char transA, const int m, const double *d_A, const int ldA, double *d_X, const int ldX, void *d_buffer, void *h_buffer);

int culme_dTrSteinBufferSize(const int m, size_t *d_bufferSize, size_t *h_bufferSize);
int culme_dTrStein(const char transA, const int m, const double *d_A, const int ldA, double *d_X, const int ldX, void *d_buffer, void *h_buffer);

int culme_dTgSylvBufferSize(const int m, const int n, size_t *d_bufferSize, size_t *h_bufferSize);
int culme_dTgSylv(const char transA, const char transB, const double sgn, const int m, const int n, const double *d_A, const int ldA, const double *d_B, const int ldB, const double *d_C, const int ldC, const double *d_D, const int ldD, double *d_X, const int ldX, void *d_buffer, void *h_buffer);

int culme_dTgLyapBufferSize(const int m, size_t *d_bufferSize, size_t *h_bufferSize);
int culme_dTgLyap(const char transA, const int m, const double *d_A, const int ldA, const double *d_B, const int ldB, double *d_X, const int ldX, void *d_buffer, void *h_buffer);

int culme_dTgSteinBufferSize(const int m, size_t *d_bufferSize, size_t *h_bufferSize);
int culme_dTgStein(const char transA, const int m, const double *d_A, const int ldA, const double *d_B, const int ldB, double *d_X, const int ldX, void *d_buffer, void *h_buffer);

int culme_dTgCSylvBufferSize(const int m, const int n, size_t *d_bufferSize, size_t *h_bufferSize);
int culme_dTgCSylv(const char transA, const char transB, const double sgn1, const double sgn2, const int m, const int n, const double *d_A, const int ldA, const double *d_B, const int ldB, const double *d_C, const int ldC, const double *d_D, const int ldD, double *d_E, const int ldE, double *d_F, const int ldF, void *d_buffer, void *h_buffer);

int culme_dTgCSylvDualBufferSize(const int m, const int n, size_t *d_bufferSize, size_t *h_bufferSize);
int culme_dTgCSylvDual(const char transA, const char transB, const double sgn1, const double sgn2, const int m, const int n, const double *d_A, const int ldA, const double *d_B, const int ldB, const double *d_C, const int ldC, const double *d_D, const int ldD, double *d_E, const int ldE, double *d_F, const int ldF, void *d_buffer, void *h_buffer);

/* general solvers, dense coefficient matrices */
int culme_dGeLyapBufferSize(const int m, size_t *d_bufferSize, size_t *h_bufferSize);
int culme_dGeLyap(const char transA, const int m, const double *d_A, const int ldA, double *d_X, const int ldX, void *d_buffer, void *h_buffer);

int culme_dGeLyapBatchedBufferSize(const int m, const int batchCount, size_t *d_bufferSize, size_t *h_bufferSize);
int culme_dGeLyapBatched(const char transA, const int m, const double *d_A, const int ldA, const long long strideA, double *d_X, const int ldX, const long long strideX, const int batchCount, void *d_buffer, void *h_buffer);

int culme_dGeSylvBatchedBufferSize(const int m, const int n, const int batchCount, size_t *d_bufferSize, size_t *h_bufferSize);
int culme_dGeSylvBatched(const char transA, const char transB, const double sgn, const int m, const int n, const double *d_A, const int ldA, const long long strideA, const double *d_B, const int ldB, const long long strideB, double *d_X, const int ldX, const long long strideX, const int batchCount, void *d_buffer, void *h_buffer);

int culme_dGeSteinBatchedBufferSize(const int m, const int batchCount, size_t *d_bufferSize, size_t *h_bufferSize);
int culme_dGeSteinBatched(const char transA, const int m, const double *d_A, const int ldA, const long long strideA, double *d_X, const int ldX, const long long strideX, const int batchCount, void *d_buffer, void *h_buffer);

int culme_dGeSylv2BatchedBufferSize(const int m, const int n, const int batchCount, size_t *d_bufferSize, size_t *h_bufferSize);
int culme_dGeSylv2Batched(const char transA, const char transB, const double sgn, const int m, const int n, const double *d_A, const int ldA, const long long strideA, const double *d_B, const int ldB, const long long strideB, double *d_X, const int ldX, const long long strideX, const int batchCount, void *d_buffer, void *h_buffer);

int culme_dGeSteinBufferSize(const int m, size_t *d_bufferSize, size_t *h_bufferSize);
int culme_dGeStein(const char transA, const int m, const double *d_A, const int ldA, double *d_X, const int ldX, void *d_buffer, void *h_buffer);

int culme_dGeSylvBufferSize(const int m, const int n, size_t *d_bufferSize, size_t *h_bufferSize);
int culme_dGeSylv(const char transA, const char transB, const double sgn, const int m, const int n, const double *d_A, const int ldA, const double *d_B, const int ldB, double *d_X, const int ldX, void *d_buffer, void *h_buffer);

int culme_dGeSylv2BufferSize(const int m, const int n, size_t *d_bufferSize, size_t *h_bufferSize);
int culme_dGeSylv2(const char transA, const char transB, const double sgn, const int m, const int n, const double *d_A, const int ldA, const double *d_B, const int ldB, double *d_X, const int ldX, void *d_buffer, void *h_buffer);

int culme_dGgLyapBufferSize(const int m, size_t *d_bufferSize, size_t *h_bufferSize);
int culme_dGgLyap(const char transA, const int m, const double *d_A, const int ldA, const double *d_B, const int ldB, double *d_X, const int ldX, void *d_buffer, void *h_buffer);

int culme_dGgSteinBufferSize(const int m, size_t *d_bufferSize, size_t *h_bufferSize);
int culme_dGgStein(const char transA, const int m, const double *d_A, const int ldA, const double *d_B, const int ldB, double *d_X, const int ldX, void *d_buffer, void *h_buffer);

int culme_dGgSylvBufferSize(const int m, const int n, size_t *d_bufferSize, size_t *h_bufferSize);
int culme_dGgSylv(const char transA, const char transB, const double sgn, const int m, const int n, const double *d_A, const int ldA, const double *d_B, const int ldB, const double *d_C, const int ldC, const double *d_D, const int ldD, double *d_X, const int ldX, void *d_buffer, void *h_buffer);

int culme_dGgCSylvBufferSize(const int m, const int n, size_t *d_bufferSize, size_t *h_bufferSize);
int culme_dGgCSylv(const char transA, const char transB, const double sgn1, const double sgn2, const int m, const int n, const double *d_A, const int ldA, const double *d_B, const int ldB, const double *d_C, const int ldC, const double *d_D, const int ldD, double *d_E, const int ldE, double *d_F, const int ldF, void *d_buffer, void *h_buffer);

int culme_dGgCSylvDualBufferSize(const int m, const int n, size_t *d_bufferSize, size_t *h_bufferSize);
int culme_dGgCSylvDual(const char transA, const char transB, const double sgn1, const double sgn2, const int m, const int n, const double *d_A, const int ldA, const double *d_B, const int ldB, const double *d_C, const int ldC, const double *d_D, const int ldD, double *d_E, const int ldE, double *d_F, const int ldF, void *d_buffer, void *h_buffer);

Single, Complex Single and Complex Double Precision Functions

All functions are also available with the prefixes s (float, sgn as float), c (cuComplex, sgn as float) and z (cuDoubleComplex, sgn as double), e.g.

int culme_sTrSylv(const char transA, const char transB, const float sgn, const int m, const int n, const float *d_A, const int ldA, const float *d_B, const int ldB, float *d_X, const int ldX, void *d_buffer, void *h_buffer);
int culme_cTrSylv(const char transA, const char transB, const float sgn, const int m, const int n, const cuComplex *d_A, const int ldA, const cuComplex *d_B, const int ldB, cuComplex *d_X, const int ldX, void *d_buffer, void *h_buffer);
int culme_zTrSylv(const char transA, const char transB, const double sgn, const int m, const int n, const cuDoubleComplex *d_A, const int ldA, const cuDoubleComplex *d_B, const int ldB, cuDoubleComplex *d_X, const int ldX, void *d_buffer, void *h_buffer);

See culme.h for the complete list.

Algorithms

Triangular Solvers

The triangular solvers implement the substitution phase of the Bartels-Stewart algorithm as a blocked level-3 scheme:

Bartels, R. H., & Stewart, G. W. (1972). Solution of the matrix equation AX + XB = C. Communications of the ACM, 15(9), 820-826.

The matrices are partitioned into diagonal blocks (block size 16-32 depending on the data type, 2x2 bumps of the real Schur form are never split). Each diagonal block equation is solved by a single fused kernel which loads the coefficient tiles into shared memory, decomposes the block into 1x1/2x2 micro blocks, solves the tiny Kronecker systems (size at most 4, or 8 for the coupled equations) by Gaussian elimination with partial pivoting in registers and applies the in-tile updates with one thread per tile entry. The remaining right-hand side is updated by cublas<t>gemm, so the bulk of the arithmetic runs at level-3 speed. A single generic two-sided core $op(L_1),X,op(R_1) + \mathrm{sgn},op(L_2),X,op(R_2) = Y$ instantiates all non-coupled equations; a second core handles the coupled CSYLV/CSYLV_DUAL forms.

The solvers detect singular (or numerically singular) micro systems and return a positive error code in that case.

General Solvers

Since cuSOLVER provides no non-symmetric (generalized) Schur decomposition, the general solvers replace the reduction phase of the Bartels-Stewart approach with GPU-friendly iterations built from gemm, geam and LU factorizations (cusolverDnXgetrf/getrs):

Solver	Algorithm	Requirement
`GeLyap`	sign function Newton iteration	$\Lambda(A)$ in the open left or right half plane
`GeSylv`	sign function Newton iteration	$\Lambda(op_1(A))$ and $\Lambda(\mathrm{sgn},op_2(B))$ in the same open half plane
`GeStein`	squared Smith iteration	$\rho(A) < 1$
`GeSylv2`	squared Smith iteration	$\rho(A),\rho(B) < 1$
`GgLyap`	reduction by $B^{-1}$ + sign	$B$ invertible, $\Lambda(A, B)$ in the open left or right half plane
`GgStein`	reduction by $B^{-1}$ + Smith	$B$ invertible, $\Lambda(A, B)$ inside the unit circle
`GgSylv`	reduction by $C^{-1}$, $B^{-1}$ + sign	$B, C$ invertible, $\Lambda(A, C)$ and $\mathrm{sgn},\Lambda(D, B)$ in the same open half plane
`GgCSylv`	elimination of $L$ + Smith	$B, C$ invertible, $\rho(C^{-1}A),\rho(D B^{-1}) < 1$
`GgCSylvDual`	elimination of $R$ + Smith	$A, D$ invertible, $\rho(A^{-H}C^{H}),\rho(B^{H}D^{-H}) < 1$

The sign function Newton iteration with determinantal scaling

$$A_{k+1} = \frac{1}{2}\left(\frac{A_k}{c_k} + c_k A_k^{-1}\right), \qquad c_k = |\det A_k|^{1/m}$$

converges globally quadratically to $\mathrm{sign}(A)$; the right-hand side is iterated alongside and the solution is recovered from the limit, see

Roberts, J. D. (1980). Linear model reduction and solution of the algebraic Riccati equation by use of the sign function. International Journal of Control, 32(4), 677-687.

Benner, P., & Quintana-Ortí, E. S. (1999). Solving stable generalized Lyapunov equations with the matrix sign function. Numerical Algorithms, 20(1), 75-100.

The explicit inverses inside the sign iterations avoid the partial pivoting bottleneck of getrf on GPUs: for medium problem sizes the iterate is transformed with a depth-2 random butterfly transform (RBT), factorized without pivoting, and the accuracy is restored by a single Newton refinement step of the inverse plus a cheap residual probe with a partial-pivoted fallback. The final iterations near convergence exploit that the iterate approximates its own inverse ($A_k \approx \mathrm{sign}(A) = A_k^{-1}$): a few factorization-free Schulz steps $X \leftarrow X (2I - A_k X)$ starting from $X = A_k$ square the residual and are accepted through the same probe, with the pivoted reference path as fallback. The one-shot reductions of the generalized (Gg*) solvers reuse the same probed inversion helper, so their solves collapse into gemms. See

Parker, D. S. (1995). Random butterfly transformations with applications in computational linear algebra. UCLA Technical Report CSD-950023.

Baboulin, M., Dongarra, J., Hermann, J., Tomov, S. (2013). Accelerating linear system solutions using randomization techniques. ACM Transactions on Mathematical Software, 39(2).

The squared Smith iteration doubles the partial sums of $X = \mathrm{sgn} \sum_{j \geq 0} (-\mathrm{sgn})^{j} A^{j} Y B^{j}$ by repeated squaring of the coefficients and converges quadratically for $\rho(A)\rho(B) < 1$, see

Smith, R. A. (1968). Matrix equation XA + BX = C. SIAM Journal on Applied Mathematics, 16(1), 198-201.

The determinantal scaling factors and the convergence control require one host synchronization per iteration; everything else runs asynchronously on the configured stream. The iterations return a positive error code if an iterate is singular (spectral assumptions violated, e.g. eigenvalues on the imaginary axis) or if the iteration does not converge.

Batched Solvers

For many small independent equations, individual Ge* calls cannot saturate the device. The batched solvers drive the whole batch with one iteration in strided batched storage (layout as in cublas<t>gemmStridedBatched: problem i lives at d_A + i * strideA, d_B + i * strideB and d_X + i * strideX):

Solver	Equation	Iteration
`GeLyapBatched`	$op_1(A_i)X_i + X_i op_1(A_i)^H = Y_i$	sign function
`GeSylvBatched`	$op_1(A_i)X_i + \text{sgn } X_i op_2(B_i) = Y_i$	sign function
`GeSteinBatched`	$op_1(A_i)X_i op_1(A_i)^H - X_i = Y_i$	squared Smith
`GeSylv2Batched`	$op_1(A_i)X_i op_2(B_i) + \text{sgn } X_i = Y_i$	squared Smith

int culme_dGeLyapBatchedBufferSize(const int m, const int batchCount, size_t *d_bufferSize, size_t *h_bufferSize);
int culme_dGeLyapBatched(const char transA, const int m, const double *d_A, const int ldA, const long long strideA, double *d_X, const int ldX, const long long strideX, const int batchCount, void *d_buffer, void *h_buffer);

For the sign function solvers the inversions run through cublas<t>getrf/getriBatched, the per-problem determinantal scaling is computed from the LU diagonals on the device and the rational updates collapse into elementwise kernels plus two strided batched gemms, so every Newton step costs a fixed number of kernel launches for the whole batch. Convergence is tracked per problem: an entry that reaches its tolerance runs one final accurate pass and is then frozen (its iterate and right-hand side stop changing) while the remaining problems keep iterating, so slowly converging entries cannot degrade the already converged ones. The squared Smith solvers consist of strided batched gemms only, converged entries decay on their own (their update contributions shrink quadratically below round-off). The spectral requirements of the corresponding single-problem solvers apply per batch entry. Intended for many small problems (dimensions up to a few hundred): on an H100, 512 problems with m = 32 in double precision solve in 0.7 ms (GeLyapBatched) and 0.4 ms (GeSteinBatched) versus 400 ms for the equivalent GeLyap loop.

Streams

All solvers run on the CUDA stream set by

int culme_setStream(cudaStream_t stream);
cudaStream_t culme_stream(void);

The stream is stored per thread, so different host threads can solve independent equations concurrently on different streams.

CUDA Graphs (experimental)

Repeated solves of the same problem shape can opt into CUDA graph capture:

int culme_setGraphMode(int enable);
int culme_graphMode(void);

When enabled, the substitution phase of the triangular solvers (Tr*, Tg*) is captured into a CUDA graph on first use and replayed on subsequent calls with the same configuration (equation, dimensions, block structure and device pointers), which removes the per-call kernel launch overhead. A non-default stream set via culme_setStream is required, otherwise the solvers fall back to the regular execution. The graph cache is stored per thread and owns a small amount of additional device memory per cached configuration; disabling the graph mode clears the cache.

Reproducibility

The randomized components (the random butterfly transforms of the sign iterations and the randomized range sampling of the Schur divide and conquer) draw their random numbers from a thread-local generator:

int culme_setSeed(unsigned long long seed);
unsigned long long culme_seed(void);

Calling culme_setSeed before a solve makes the run bitwise reproducible on the same hardware and toolkit version.

Error Codes

All functions return 0 on success. Shared return codes: -1 invalid argument, -2/-3/-4 CUDA / cuBLAS / cuSOLVER errors (details on stderr), 1 singular input or iterate (e.g. violated spectral assumptions), 2 no convergence, 3 the Schur prototype could not split a clustered spectrum. A human readable description is returned by

const char *culme_errorString(int error);

Algebraic Riccati Equations

The continuous-time and discrete-time algebraic Riccati equations with general dense coefficients are solved for all four data types (G and Q Hermitian positive semidefinite, the standard stabilizability/detectability assumptions apply). The right-hand side Q is passed in d_X and overwritten by the stabilizing solution:

/* CARE: A^H X + X A - X G X + Q = 0          (trans = 'N')
 *       A X + X A^H - X G X + Q = 0          (trans = 'T'/'C') */
int culme_dGeCARE(const char trans, const int m, const double *d_A, const int ldA,
                  const double *d_G, const int ldG, double *d_X, const int ldX,
                  void *d_buffer, void *h_buffer);

/* DARE: X = A^H X (I + G X)^-1 A + Q         (trans = 'N')
 *       X = A X (I + G X)^-1 A^H + Q         (trans = 'T'/'C') */
int culme_dGeDARE(const char trans, const int m, const double *d_A, const int ldA,
                  const double *d_G, const int ldG, double *d_X, const int ldX,
                  void *d_buffer, void *h_buffer);

(plus the s, c, z variants and the *BufferSize queries). The CARE solver computes the matrix sign function of the $2m \times 2m$ Hamiltonian $\begin{pmatrix} A & -G \ -Q & -A^H \end{pmatrix}$, extracts the stabilizing solution from the stable invariant subspace and polishes it with Newton steps, each solved by the sign function Lyapunov solver (the closed loop $A - GX$ is stable by construction). The DARE solver runs the quadratically convergent structure preserving doubling algorithm, which consists almost entirely of GEMMs and probed RBT inversions. For the control form $G = B R^{-1} B^H$ build G with one factorization of the small $p \times p$ matrix R beforehand.

Separation Estimation

The sensitivity of a Sylvester or Lyapunov solution is governed by the separation $\mathrm{sep}(T) = \sigma_{\min}(T)$ of the underlying linear operator, e.g. $T: X \mapsto op_1(A)X + \mathrm{sgn},X,op_2(B)$: a backward error $\eta$ in the coefficients perturbs the solution by up to $\eta,(\lVert A\rVert + \lVert B\rVert)/\mathrm{sep}$. cuLME estimates the separation of the triangular operators with a few power iteration steps on the inverse operator (two triangular solves per step):

int culme_dTrSylvSep(const char transA, const char transB, const double sgn, const int m, const int n,
                     const double *d_A, const int ldA, const double *d_B, const int ldB,
                     double *sep, void *d_buffer, void *h_buffer);
int culme_dTrLyapSep(const char transA, const int m, const double *d_A, const int ldA,
                     double *sep, void *d_buffer, void *h_buffer);

(plus the s, c, z variants and the *BufferSize queries). The estimate approaches the true separation from above and reports *sep = 0 for singular operators. For general coefficients reduce to triangular form first (culme_dSchur), the separation is invariant under the orthogonal reduction.

Tuning

The block sizes and algorithmic windows are resolved once per device from the compute capability: data center GPUs (full rate FP64) use the default windows, on consumer and visualization parts (throttled FP64) the RBT inversion window is narrowed because its O(m³) FP64 refinement outgrows the saved pivoting earlier. The defaults can be overridden at run time without a rebuild:

Environment variable	Effect
`CULME_RBT_MIN`	lower bound of the RBT inversion window (default 384)
`CULME_RBT_MAX`	upper bound for real types, complex types use half of it
`CULME_SCHUR_NB`	Schur base case size, clamped to the kernel capacity (8..80, 64 for `cuDoubleComplex`)
`CULME_TRRECTH`	recursive blocking threshold of the triangular solvers (default 4096)
`CULME_NO_RBT`	disable the RBT fast path (always use pivoted LU)
`CULME_NO_NS`	disable the Newton-Schulz endgame of the sign iteration

Schur decomposition (experimental prototype)

A prototype of the Schur decomposition $A = Q T Q^\mathsf{H}$ on the GPU is available for all four data types:

int culme_dSchurBufferSize(const int m, size_t *d_bufferSize, size_t *h_bufferSize);
int culme_dSchur(const int m, double *d_A, const int ldA, double *d_Q, const int ldQ, void *d_buffer, void *h_buffer);

(plus the s, c, z variants). d_A is overwritten by the Schur form T (real types: upper quasi-triangular with 1x1/2x2 blocks, complex types: upper triangular, directly usable by the Tr*/Tg* solvers), d_Q receives the unitary Schur vectors. The implementation uses spectral divide and conquer over the matrix sign function (sign Newton iteration with determinantal scaling, randomized range sampling and Householder QR for the spectral projector bases, and a block-parallel base case kernel for diagonal blocks up to size 80 — Hessenberg/Francis double shift QR for real types, Hessenberg/single shift QR with Givens rotations for complex types). The split lines are placed into the largest gaps of the diagonal value distribution. Real types split along vertical lines first and separate spectra clustered on a single vertical line in a second stage through the squared shifted block ((M - aI)^2), whose sign splits sort the eigenvalues by their imaginary parts. Complex types split along vertical and horizontal lines (via the rotated iterate (-i(M - i\tau I))) before falling back to the squared shift. A single (m)-fold eigenvalue beyond the base case size (e.g. one large Jordan block) admits no spectral split at all and the routine returns 3. See example_culme_schur for the accuracy evaluation.

On top of the prototype, Bartels-Stewart style frontends for the non-generalized equations are available for all four data types:

int culme_dGeLyapSchur(const char fact, const char transA, const int m, double *d_A, const int ldA, double *d_Q, const int ldQ, double *d_X, const int ldX, void *d_buffer, void *h_buffer);
int culme_dGeSteinSchur(const char fact, const char transA, const int m, double *d_A, const int ldA, double *d_Q, const int ldQ, double *d_X, const int ldX, void *d_buffer, void *h_buffer);
int culme_dGeSylvSchur(const char factA, const char factB, const char transA, const char transB, const double sgn, const int m, const int n, double *d_A, const int ldA, double *d_QA, const int ldQA, double *d_B, const int ldB, double *d_QB, const int ldQB, double *d_X, const int ldX, void *d_buffer, void *h_buffer);
int culme_dGeSylv2Schur(const char factA, const char factB, const char transA, const char transB, const double sgn, const int m, const int n, double *d_A, const int ldA, double *d_QA, const int ldQA, double *d_B, const int ldB, double *d_QB, const int ldQB, double *d_X, const int ldX, void *d_buffer, void *h_buffer);

(plus the s, c, z variants).

They reduce the coefficients to Schur form, transform the right-hand side, run the corresponding triangular solver and transform back. Unlike the iterative Ge* solvers they carry no spectral requirements beyond unique solvability (the sign function iteration needs all eigenvalues in a common open half plane, the Smith iteration needs spectral radii below one). With fact = 'F' the Schur factors of a previous call ('N' overwrites d_A with T and fills d_Q) are reused, so repeated solves with further right-hand sides skip the expensive reduction and only cost the triangular solve plus four gemms.

Project Layout

include/     public API header culme.h
src/         library implementation (.cu) and internal headers
examples/    test suite (example_culme) and Schur evaluation (example_culme_schur)
tests/       gtest unit tests of the internal numerics and the coverage report script
benchmarks/  performance benchmark (benchmark_culme)
thirdParty/  optional vendored dependencies (googletest), not part of the repository
cmake/       CMake package configuration template

Code Tour

Shared infrastructure (internal headers):

src/checkcuda.h       error propagation macros CHECK_* (early-return chain, maps to the public error codes)
src/culme_math.h      host/device scalar arithmetic for the four element types (cm_add, cm_mul, cm_conj, ...)
src/culme_traits.h    culme_traits<T>: real scalar type, eps, block sizes, cuBLAS/cuSolver type ids and dispatch
src/culme_handles.h   thread-local state: stream, RNG seed, pinned staging scratch, per-device tuning, cached handles
src/culme_graphs.h    CUDA graph cache of the triangular substitution phase (opt-in graph mode)
src/culme_trcore.h    level-3 wavefront substitution engine of all (quasi-)triangular equations:
                      diagonal block micro kernels, grouped batched gemm updates, graph capture
src/culme_gecore.h    workspace arena, sign function Newton and squared Smith iteration cores,
                      RBT/non-pivoted fast inversion with Newton refinement and residual probe
src/culme_rescore.h   shared residual evaluation core

Solver translation units (each ends with the macro-stamped s/d/c/z C API wrappers):

src/culme_trlyap.cu    TrLyap, TrStein            triangular standard equations
src/culme_trsylv.cu    TrSylv, TrSylv2            triangular Sylvester equations
src/culme_tgsylv.cu    TgSylv, TgLyap, TgStein    generalized triangular equations
src/culme_tgcsylv.cu   TgCSylv, TgCSylvDual       coupled generalized triangular equations
src/culme_gelyap.cu    GeLyap, GeStein            general dense, sign/Smith iteration
src/culme_gesylv.cu    GeSylv, GeSylv2            general dense, sign/Smith iteration
src/culme_gglyap.cu    GgLyap, GgStein            generalized dense, reduction with B^-1
src/culme_ggsylv.cu    GgSylv                     generalized dense, reduction with C^-1 and B^-1
src/culme_ggcsylv.cu   GgCSylv, GgCSylvDual       coupled generalized dense, elimination + Smith
src/culme_gebatched.cu Ge*Batched                 batched sign/Smith iterations on strided stacks
src/culme_riccati.cu   GeCARE, GeDARE             Riccati: Hamiltonian sign + Newton, SDA doubling
src/culme_schur.cu     Schur                      spectral divide and conquer prototype
src/culme_geschur.cu   Ge*Schur                   Bartels-Stewart frontends over Schur + Tr*
src/culme_sep.cu       TrSylvSep, TrLyapSep       separation estimation (inverse power iteration)
src/culme_residual.cu  *Residual                  relative residuals of all equations
src/culme_stream.cu    stream/graph/seed state and culme_errorString

Every file starts with a block comment describing the algorithm, the mathematical background with references, and the workspace layout; the internal headers additionally document the maintainer invariants (e.g. the *BufferSize/solver arena contract in culme_gecore.h).

Installation

Prerequisites:

CMake >= 3.23
CUDA >= 11.4.2

  mkdir build && cd build
  cmake ..
  make
  make install

The build can be configured with the following options:

Option	Default	Description
`CULME_BUILD_EXAMPLE`	`ON`	build the example `example_culme`
`CULME_BUILD_BENCHMARK`	`OFF`	build the benchmark `benchmark_culme`
`CULME_BUILD_TESTS`	`ON`	build the gtest unit tests `test_culme`
`CULME_COVERAGE`	`OFF`	instrument the host code paths for gcov

If CMAKE_CUDA_ARCHITECTURES is not specified, the library is built for all supported architectures.

cuLME installs a CMake package configuration, so downstream projects can simply use

find_package(culme REQUIRED)
target_link_libraries(myproject PRIVATE culme::culme)

Usage and Examples

The example example_culme.cu covers all 18 solvers and all supported matrix formats (real and complex, single and double precision). For each equation a known solution is generated, the right-hand side is assembled on the host and the relative forward error of the GPU solution is checked:

  ./example_culme            # full suite including the size sweep (~1 min)
  ./example_culme --quick    # default sizes only

A minimal usage example for the double precision Lyapunov solver with a general stable matrix A:

size_t d_bufferSize, h_bufferSize;
culme_dGeLyapBufferSize(m, &d_bufferSize, &h_bufferSize);

void *d_buffer, *h_buffer;
cudaMalloc(&d_buffer, d_bufferSize);
h_buffer = malloc(h_bufferSize);

/* d_X contains the right-hand side Y on entry and the solution X on exit */
int ret = culme_dGeLyap('N', m, d_A, m, d_X, m, d_buffer, h_buffer);

Testing

The test suite is registered with CTest and consists of five runs:

Test	Description
`culme.full`	all solvers, forward error and GPU residual checks, plus a size sweep and edge case tests
`culme.schur`	accuracy evaluation of the Schur prototype (backward error, orthogonality, known spectra)
`culme.no_rbt`	quick run with `CULME_NO_RBT=1`, forces the pivoted-LU path instead of the RBT fast path
`culme.no_ns`	quick run with `CULME_NO_NS=1`, disables the Newton-Schulz endgame of the sign iteration
`culme.unit`	gtest unit tests of the internal numerical building blocks

  cd build && ctest --output-on-failure

The full suite runs every solver (including the Schur frontends and the four batched solvers in all four data types) over 30 (m, n) pairs covering the tiny sizes 1 to 5, odd and non-power-of-two sizes, sizes straddling the internal block boundaries (triangular kernels NB = 16, Schur base case NB = 80) and one pair inside the RBT window (m >= 384), in both the default mode and the CUDA graph mode. Edge case tests check empty problems (m = 0), argument validation and defect inputs (singular coefficients) for clean error codes — about 11,000 checks in total. The reference right-hand sides and residuals of the example are computed with cuBLAS/cuSolver, so the whole suite runs on the GPU. A GPU is required to run the tests.

Unit Tests

The unit tests in tests/test_internals.cu check mathematical invariants of the internal primitives directly instead of only end-to-end solver errors, for all four data types:

LogDetMatchesHostReference — log|det| of the pivoted LU against a host reference (catches modulus errors like |Re| + |Im| instead of hypot underneath the logarithm)
RbtInvertResidual — ||A * Inv - I||_F of all inversion paths (pivoted, RBT window, Schulz endgame with probe fallback) plus the returned log|det|
SignIterationInvariants — the sign function Newton iteration on a spectrum with known half-plane split: iteration budget (quadratic convergence), ||S^2 - I||_F and the exact trace
PrecisionScaling — the same Lyapunov problem solved in single and double precision must show errors scaling with the machine epsilon, a precision independent error floor signals an algorithmic bug

GoogleTest is acquired via FetchContent: a vendored clone under thirdParty/googletest is preferred (offline friendly), otherwise an installed package is searched and finally the sources are downloaded:

  git clone --branch v1.14.0 https://github.com/google/googletest thirdParty/googletest

Coverage

A coverage build instruments the host compiler passes of the library (device kernel bodies are not covered by gcov):

  cmake -S . -B build-cov -DCULME_COVERAGE=ON
  cmake --build build-cov -j
  (cd build-cov && ctest)
  tests/coverage_report.sh build-cov

The report aggregates the line coverage per source file over all translation units; the annotated sources in build-cov/coverage/ show the uncovered branches.

Benchmark

The benchmark benchmark_culme.cu (built with -DCULME_BUILD_BENCHMARK=ON) times the triangular Sylvester and Lyapunov solvers and the general Lyapunov solver over a sweep of problem sizes:

  ./benchmark_culme             # sweep m = 512, 1024, 2048 in single precision (s, c)
  ./benchmark_culme 4096        # single problem size
  ./benchmark_culme 1024 sdcz   # single problem size, selected precisions

Double precision is opt-in because many consumer and visualization GPUs only provide severely throttled FP64 throughput.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Repository files navigation

cuLME - Linear Matrix Equation Solvers using CUDA

Available Functions

Double Precision Functions

Single, Complex Single and Complex Double Precision Functions

Algorithms

Triangular Solvers

General Solvers

Batched Solvers

Streams

CUDA Graphs (experimental)

Reproducibility

Error Codes

Algebraic Riccati Equations

Separation Estimation

Tuning

Schur decomposition (experimental prototype)

Project Layout

Code Tour

Installation

Usage and Examples

Testing

Unit Tests

Coverage

Benchmark

About

Releases

Packages

Contributors

Languages

Name		Name	Last commit message	Last commit date
Latest commit History 1 Commit
benchmarks		benchmarks
cmake		cmake
examples		examples
include		include
src		src
tests		tests
.clang-format		.clang-format
.gitignore		.gitignore
CMakeLists.txt		CMakeLists.txt
LICENSE		LICENSE
README.md		README.md

Folders and files

Latest commit

History

Repository files navigation

cuLME - Linear Matrix Equation Solvers using CUDA

Available Functions

Double Precision Functions

Single, Complex Single and Complex Double Precision Functions

Algorithms

Triangular Solvers

General Solvers

Batched Solvers

Streams

CUDA Graphs (experimental)

Reproducibility

Error Codes

Algebraic Riccati Equations

Separation Estimation

Tuning

Schur decomposition (experimental prototype)

Project Layout

Code Tour

Installation

Usage and Examples

Testing

Unit Tests

Coverage

Benchmark

About

Resources

Stars

Watchers

Forks

Releases

Packages

Contributors

Languages