sgltgsylv(3)

📦

version 1.1.2 · 3 versions

1.1.2 · fedora 43 (latest)
1.1.1 · fedora 43
1.1.2 · fedora 44

💖 sgltgsylv(3) — MEPACK 💖

okay so bestie, gather round, because this is sgltgsylv from the MEPACK manual, section 3, and it is a LOT. It's the whole single-precision corner of the generalized Sylvester equation world. Yes, all of it. Grab a snack. 🍓

✨ NAME ✨

sgltgsylv — Single Precision

Single Precision routines for triangular generalized Sylvester equations.

no cap, that's the entire name section. short but iconic. 💅

🌸 SYNOPSIS 🌸

🎀 Functions 🎀

this is the full roster, in order, every single one, don't skip. 🫶

`sla_tgcsylv_dag`

subroutine sla_tgcsylv_dag (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info)

Level-3 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation, but make it DAG scheduling. parallel bestie behavior. 🔥

`sla_tgcsylv_dual_dag`

subroutine sla_tgcsylv_dual_dag (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info)

OpenMP4 - DAG Algorithm for the dual generalized coupled Sylvester equation. it's giving threads. 🧚

`sla_tgsylv_dag`

subroutine sla_tgsylv_dag (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info)

Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation with DAG based parallelization. 💫

`sla_tgcsylv_dual_kernel_44nn`

subroutine sla_tgcsylv_dual_kernel_44nn (sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, info)

Solver for a 4x4 dual generalized coupled Sylvester equation (TRANSA = N, TRANSB = N). tiny but mighty, the kernel that carries the whole operation. 👑

`sla_tgcsylv_dual_kernel_44nt`

subroutine sla_tgcsylv_dual_kernel_44nt (sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, info)

Solver for a 4x4 dual generalized coupled Sylvester equation (TRANSA = N, TRANSB = T). same energy, different letters. 💅

`sla_tgcsylv_dual_kernel_44tn`

subroutine sla_tgcsylv_dual_kernel_44tn (sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, info)

Solver for a 4x4 dual generalized coupled Sylvester equation (TRANSA = T, TRANSB = N). 🔁

`sla_tgcsylv_dual_kernel_44tt`

subroutine sla_tgcsylv_dual_kernel_44tt (sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, info)

Solver for a 4x4 dual generalized coupled Sylvester equation (TRANSA = T, TRANSB = T). double transposed, double slay. ✨

`sla_tgcsylv_dual_l2`

subroutine sla_tgcsylv_dual_l2 (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info)

Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the dual generalized coupled Sylvester equation. lvl 2 era. 🦋

`sla_tgcsylv_dual_l2_local_copy`

subroutine sla_tgcsylv_dual_l2_local_copy (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info)

Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the dual generalized coupled Sylvester equation (Local Copy). it copies locally, respectfully. 💖

`sla_tgcsylv_dual_l2_local_copy_128`

subroutine sla_tgcsylv_dual_l2_local_copy_128 (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info)

Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the dual generalized coupled Sylvester equation (Local Copy). the _128 one. 🎀

`sla_tgcsylv_dual_l2_local_copy_32`

subroutine sla_tgcsylv_dual_l2_local_copy_32 (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info)

Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the dual generalized coupled Sylvester equation (Local Copy). the _32 one. 🌸

`sla_tgcsylv_dual_l2_local_copy_64`

subroutine sla_tgcsylv_dual_l2_local_copy_64 (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info)

Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the dual generalized coupled Sylvester equation (Local Copy). the _64 one. ⭐

`sla_tgcsylv_dual_l2_local_copy_96`

subroutine sla_tgcsylv_dual_l2_local_copy_96 (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info)

Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the dual generalized coupled Sylvester equation (Local Copy). the _96 one, the last of the local copy quartet. 🙈

`sla_tgcsylv_l2`

subroutine sla_tgcsylv_l2 (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info)

Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation (Optimized). 💚

`sla_tgcsylv_l2_local_copy`

subroutine sla_tgcsylv_l2_local_copy (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info)

Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation (Optimized). yes, same label, different routine, keep up bestie. 🫶

`sla_tgcsylv_l2_local_copy_128`

subroutine sla_tgcsylv_l2_local_copy_128 (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info)

Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation (Optimized). _128. 🎀

`sla_tgcsylv_l2_local_copy_32`

subroutine sla_tgcsylv_l2_local_copy_32 (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info)

Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation (Optimized). _32. 🌸

`sla_tgcsylv_l2_local_copy_64`

subroutine sla_tgcsylv_l2_local_copy_64 (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info)

Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation (Optimized). _64. ⭐

`sla_tgcsylv_l2_local_copy_96`

subroutine sla_tgcsylv_l2_local_copy_96 (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info)

Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation (Optimized). _96. 🙈

`sla_tgcsylv_l2_reorder`

subroutine sla_tgcsylv_l2_reorder (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info)

Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation (Optimized). reordering, as a treat. 🔁

`sla_tgcsylv_l2_unopt`

subroutine sla_tgcsylv_l2_unopt (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info)

Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation (unoptimized). it's the reference-y one, no shame in it. 💛

`sla_tgsylv_gardiner_laub`

subroutine sla_tgsylv_gardiner_laub (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info)

Level-2 Gardiner-Laub algorithm for the generalized Sylvester equation. named algorithm, extra iconic. 👑

`sla_tgsylv_l2`

subroutine sla_tgsylv_l2 (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info)

Level-2 Bartels-Stewart Algorithm for the generalized Sylvester equation (Optimized). 💚

`sla_tgsylv_l2_colwise`

subroutine sla_tgsylv_l2_colwise (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info)

Level-2 Bartels-Stewart Algorithm for the generalized Sylvester equation (unoptimized). going column by column, one at a time, like a lil ladybug. 🐞

`sla_tgsylv_l2_local_copy`

subroutine sla_tgsylv_l2_local_copy (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info)

Level-2 Bartels-Stewart Algorithm for the generalized Sylvester equation (Optimized). 💖

`sla_tgsylv_l2_local_copy_128`

subroutine sla_tgsylv_l2_local_copy_128 (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info)

Level-2 Bartels-Stewart Algorithm for the generalized Sylvester equation (Optimized). _128. 🎀

`sla_tgsylv_l2_local_copy_32`

subroutine sla_tgsylv_l2_local_copy_32 (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info)

Level-2 Bartels-Stewart Algorithm for the generalized Sylvester equation (Optimized). _32. 🌸

`sla_tgsylv_l2_local_copy_64`

subroutine sla_tgsylv_l2_local_copy_64 (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info)

Level-2 Bartels-Stewart Algorithm for the generalized Sylvester equation (Optimized). _64. ⭐

`sla_tgsylv_l2_local_copy_96`

subroutine sla_tgsylv_l2_local_copy_96 (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info)

Level-2 Bartels-Stewart Algorithm for the generalized Sylvester equation (Optimized). _96. 🙈

`sla_tgsylv_l2_reorder`

subroutine sla_tgsylv_l2_reorder (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info)

Level-2 Bartels-Stewart Algorithm for the generalized Sylvester equation (Optimized). 🔁

`sla_tgsylv_l2_rowwise`

subroutine sla_tgsylv_l2_rowwise (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info)

Level-2 Bartels-Stewart Algorithm for the generalized Sylvester equation (unoptimized). row by row this time. we contain multitudes. 🐛

`sla_tgcsylv_dual_l3`

subroutine sla_tgcsylv_dual_l3 (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info)

Level-3 Kagstroem-Westin/Bartels-Stewart Algorithm for the dual generalized coupled Sylvester equation. level 3 era, glow up complete. 🔥

`sla_tgcsylv_dual_l3_2s`

subroutine sla_tgcsylv_dual_l3_2s (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info)

Level-3 Kagstroem-Westin/Bartels-Stewart Algorithm with sub blocking for the dual generalized coupled Sylvester equation. sub blocking, sub slaying. 💅

`sla_tgcsylv_l3_2s`

subroutine sla_tgcsylv_l3_2s (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info)

Level-3 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation (optimized). 🎀

`sla_tgcsylv_l3`

subroutine sla_tgcsylv_l3 (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info)

Level-3 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation (optimized). the main character, lowkey. 👑

`sla_tgcsylv_l3_unopt`

subroutine sla_tgcsylv_l3_unopt (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info)

Level-3 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation (unoptimized). 💛

`sla_tgsylv_l3_2s`

subroutine sla_tgsylv_l3_2s (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info)

Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation. that's the whole description, it's confident. ✨

`sla_tgsylv_l3_colwise`

subroutine sla_tgsylv_l3_colwise (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info)

Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation. column-wise variant, still level-3. 🐞

`sla_tgsylv_l3`

subroutine sla_tgsylv_l3 (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info)

Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation. 🔥

`sla_tgsylv_l3_rowwise`

subroutine sla_tgsylv_l3_rowwise (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info)

Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation. row-wise variant, same deal. 🐛

`sla_tgcsylv_dual_recursive`

recursive subroutine sla_tgcsylv_dual_recursive (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info)

Recursive blocking Algorithm for the dual generalized coupled Sylvester equation. it calls itself, that's the whole vibe, it's a recursive subroutine bestie. 🌀

`sla_tgcsylv_recursive`

recursive subroutine sla_tgcsylv_recursive (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info)

Recursive Algorithm for the generalized coupled Sylvester equation. 🌀

`sla_tgsylv_recursive`

recursive subroutine sla_tgsylv_recursive (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info)

Level-3 Recursive Blocked Solver for the generalized Sylvester equation. 🌀

okay that's the entire SYNOPSIS and I'm not even winded. more sections incoming, don't go anywhere bestie 💖

Detailed Description 💖

okay so bestie, this whole page is about single precision routines for triangular generalized Sylvester equations. i know, i know, the name is giving "professor wrote this at 3am" but stay with me, it's actually so good.

This section contains the solvers for the standard Sylvester equation with (quasi) triangular coefficient matrices. The coefficient matrices are normally generated with the help of the Schur decomposition from LAPACK. And — important — all routines work in single precision arithmetic. no doubles in this house. 💅

✨ Function Documentation ✨

`sla_tgcsylv_dag` 🦋

subroutine sla_tgcsylv_dag (character, dimension(1) transa, character,
    dimension(1) transb, real sgn1, real sgn2, integer m, integer n, real,
    dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer
    ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d,
    integer ldd, real, dimension(lde, *) e, integer lde, real,
    dimension(ldf,*) f, integer ldf, real scale, real, dimension(*) work,
    integer info)

Level-3 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation with DAG scheduling. yes that is a lot of names. they earned them. 👑

Purpose 🎀

!>
!> SLA_TGCSYLV_DAG solves a generalized coupled  Sylvester equation of the following form
!>
!>    op1(A) * R  + SGN1 * L  * op2(B) = SCALE * E                              (1)
!>    op1(C) * R  + SGN2 * L  * op2(D) = SCALE * F
!>
!> where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular
!> matrix and B and D are N-by-N upper triangular matrices. Furthermore either B or
!> D can be quasi upper triangular as well. The right hand side Y and the solution X
!> M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) ( or (D,B) ) are
!> created by SGGES from LAPACK.
!>

Remarks 🌸

The Algorithm uses OpenMP DAG scheduling. that's the whole remark. short and iconic. ✨

Parameters 🫶

TRANSA — okay so TRANSA is CHARACTER(1). It specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)
TRANSB — TRANSB is CHARACTER(1). It specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)
SGN1 — SGN1 is REAL, allowed values: +/-1. It specifies the sign between in the first equation. (yes, "between in". we don't question the phrasing, we just pass the sign 💅)
SGN2 — SGN2 is REAL, allowed values: +/-1. It specifies the sign between in the second equation.
M — M is INTEGER. The order of the matrices A and C. M >= 0.
N — N is INTEGER. The order of the matrices B and D. N >= 0.
A — A is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular.
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — B is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular the matrix B must be upper triangular.
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — C is REAL array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular.
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — E is REAL array, dimension (LDE,N).
- On input, the matrix E contains the right hand side E.
- On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2.
- Right hand side E and the solution R are M-by-N matrices.
LDE — LDE is INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F — F is REAL array, dimension (LDF,N).
- On input, the matrix F contains the right hand side F.
- On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2.
- Right hand side F and the solution L are M-by-N matrices.
LDF — LDF is INTEGER. The leading dimension of the array F. LDE >= max(1,M). (fr fr that's what it says, i'm not touching it, LDE is right there in the source 💀)
SCALE — SCALE is REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — WORK is REAL array, dimension LWORK. Workspace for the algorithm. The workspace needs to queried before the running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO. (query first bestie, do NOT skip the query, the vibes will be off 🙈)
INFO — INFO is INTEGER.
- On input:
  - == -1 : Perform a workspace query
  - <> -1: normal operation
- On exit, workspace query:
  - < 0 : if INFO = -i, the i-th argument had an illegal value
  - >= 0: The value of INFO is the required number of elements in the workspace.
- On exit, normal operation:
  - == 0: successful exit
  - < 0: if INFO = -i, the i-th argument had an illegal value
  - > 0: The equation is not solved correctly. One of the arising inner system got singular. 😭

Author 👑

Martin Koehler, MPI Magdeburg

Date ⭐

October 2025

Definition 📍

Definition at line 217 of file sla_tgcsylv_dag.f90.

`sla_tgcsylv_dual_dag` 🧚

subroutine sla_tgcsylv_dual_dag (character, dimension(1) transa, character,
    dimension(1) transb, real sgn1, real sgn2, integer m, integer n, real,
    dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer
    ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d,
    integer ldd, real, dimension(lde, *) e, integer lde, real,
    dimension(ldf,*) f, integer ldf, real scale, real, dimension(*) work,
    integer info)

OpenMP4 - DAG Algorithm for the dual generalized coupled Sylvester equation. the sequel. the part two. the dual. 💜

Purpose 🦋

!>
!> SLA_TGCSYLV_DUAL_DAG solves a generalized coupled  Sylvester equation of the following form
!>
!>    op1(A)**T * R  + op1(C)**T * L               = SCALE * E                             (1)
!>    SGN1 * R * op2(B)**T + SGN2 * L * op2(D)** T = SCALE * F
!>
!> where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular
!> matrix and B and D are N-by-N upper triangular matrices. Furthermore either B or
!> D can be quasi upper triangular as well. The right hand sides E, F  and the solutions R, L are
!> M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) ( or (D,B) ) are
!> created by SGGES from LAPACK.
!> The equation (1) is the dual to the generalized coupled Sylvester equation
!>
!>    op1(A) * R + SGN1 * L * op2(B)  = SCALE * E                                          (2)
!>    op1(C) * R + SGN2 * L * op2(D)  = SCALE * F
!>
!> The equation (1) is the dual one to equation (2) with respect to the underlying linear system.
!> Let Z be the matrix formed by rewriting (2) into its Kronecker form. This yields
!>
!>         | kron(I, op1(A))   SGN1*kron(op2(B)**T, I) | | Vec R |   | Vec E |
!>   Z X = |                                           |*|       | = |       |
!>         | kron(I, op1(C))   SGN2*kron(op2(D)**T, I) | | Vec L |   | Vec F |
!>
!> Regarding Z**T one obtains
!>
!>            | kron(I, op1(A)**T )    kron(I, op1(C)**T)   | | Vec R |   | Vec E |
!>   Z**T X = |                                             |*|       | = |       |
!>            | SGN1*kron(op2(B), I)   SGN2*kron(op2(D), I) | | Vec L |   | Vec F |
!>
!> which belongs to the Sylvester equation (1). For this reason the parameters TRANSA and TRANSB
!> are expressed in terms of the Sylvester equation (2).
!>
!>

heads up bestie: because of that last line, TRANSA and TRANSB are expressed in terms of the Sylvester equation (2), not (1). that's not a typo, that's the whole point. the transposes are transposed. iykyk. 🔥

Parameters 🫶

TRANSA — TRANSA is CHARACTER(1). It specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)
TRANSB — TRANSB is CHARACTER(1). It specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)
SGN1 — SGN1 is REAL, allowed values: +/-1. It specifies the sign between in the first equation.
SGN2 — SGN2 is REAL, allowed values: +/-1. It specifies the sign between in the second equation.
M — M is INTEGER. The order of the matrices A and C. M >= 0.
N — N is INTEGER. The order of the matrices B and D. N >= 0.
A — A is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular.
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — B is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular the matrix B must be upper triangular.
LDB —
LDB — INTEGER energy, bestie. It's the leading dimension of the array B, and the vibe check says LDB >= max(1,N). 🦋
C — REAL array, dimension (LDC,M). And listen 👑 — the matrix C must be upper triangular. That's not a suggestion, that's the law. 💅
LDC — INTEGER, leading dimension of the array C. Requirement: LDC >= max(1,M). Periodt. ✨
D — REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. And here's the lil plot twist 🫶 — if the matrix B is already quasi-upper triangular, then the matrix D must be upper triangular. It's giving "one of you has to be the normal one." 🙈
LDD — INTEGER, the leading dimension of the array D. LDD >= max(1,N), no cap. 💖
E — REAL array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices. So it walks in as the problem and struts out as the answer — main character behavior fr fr. 🔥
LDE — INTEGER, leading dimension of the array E. LDE >= max(1,M). Slay. 💅
F — REAL array, dimension (LDF,N). On input, the matrix F holds the right hand side F. On output, the matrix F holds the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices. Same glow-up arc as E, honestly. ✨
LDF — INTEGER, leading dimension of the array F. And the doc literally says LDE >= max(1,M) here — not a typo I'm fixing, bestie, that's exactly what's written. 😭
SCALE — REAL. This one is the bodyguard 🫶 — SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0, otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — REAL array, dimension LWORK. Workspace for the algorithm. And you cannot just roll up with an unqueried workspace — the workspace needs to be queried before running the computation. The query is performed by calling the subroutine with INFO == -1 on input, and the required workspace is then returned in INFO. So basically you ask first, then you compute. Communication, queen. 👑
INFO — INTEGER, the lil status report 📋

On input:

== -1 : Perform a workspace query
<> -1: normal operation

On exit, workspace query:

< 0 : if INFO = -i, the i-th argument had an illegal value
>= 0: The value of INFO is the required number of elements in the workspace.

On exit, normal operation:

== 0: successful exit
< 0: if INFO = -i, the i-th argument had an illegal value
> 0: The equation is not solved correctly. One of the arising inner system got singular.

Ngl, that > 0 case is the "it's not me, it's the singular inner system" moment. Not great, not terrible. 💀

Author 🖋️

Martin Koehler, MPI Magdeburg

Date 📅

October 2025

Definition at line 237 of file sla_tgcsylv_dual_dag.f90.

💖 subroutine `sla_tgcsylv_dual_kernel_44nn` 💖

subroutine sla_tgcsylv_dual_kernel_44nn (real sgn1, real sgn2, integer m,
    integer n, real, dimension(lda, *) a, integer lda, real, dimension(ldb,
    *) b, integer ldb, real, dimension(ldc, *) c, integer ldc, real,
    dimension(ldd, *) d, integer ldd, real, dimension(lde, *) e, integer
    lde, real, dimension(ldf,*) f, integer ldf, real scale, integer info)

Solver for a 4x4 dual generalized coupled Sylvester equation (TRANSA = N, TRANSB = N). Yes bestie, a 4x4. Tiny. Adorable. The little baby of the recursion. 🎀

Purpose

!>
!> SLA_TGCSYLV_DUAL_KERNEL_44NN solves a generalized coupled  Sylvester equation of the following form
!>
!>    op1(A)**T * R  + op1(C)**T * L               = SCALE * E                             (1)
!>    SGN1 * R * op2(B)**T + SGN2 * L * op2(D)** T = SCALE * F
!>
!> where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular
!> matrix and B and D are N-by-N upper triangular matrices. Furthermore either B or
!> D can be quasi upper triangular as well. The right hand side Y and the solution X
!> M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) ( or (D,B) ) are
!> created by SGGES from LAPACK.
!> The equation (1) is the dual to the generalized coupled Sylvester equation
!>
!>    op1(A) * R + SGN1 * L * op2(B)  = SCALE * E                                          (2)
!>    op1(C) * R + SGN2 * L * op2(D)  = SCALE * F
!>
!> The equation (1) is the dual one to equation (2) with respect to the underlying linear system.
!> Let Z be the matrix formed by rewriting (2) into its Kronecker form. This yields
!>
!>         | kron(I, op1(A))   SGN1*kron(op2(B)**T, I) | | Vec R |   | Vec E |
!>   Z X = |                                           |*|       | = |       |
!>         | kron(I, op1(C))   SGN2*kron(op2(D)**T, I) | | Vec L |   | Vec F |
!>
!> Regarding Z**T one obtains
!>
!>            | kron(I, op1(A)**T )    kron(I, op1(C)**T)   | | Vec R |   | Vec E |
!>   Z**T X = |                                             |*|       | = |       |
!>            | SGN1*kron(op2(B), I)   SGN2*kron(op2(D), I) | | Vec L |   | Vec F |
!>
!> which belongs to the Sylvester equation (1). For this reason the parameters TRANSA and TRANSB
!> are expressed in terms of the Sylvester equation (2).
!>
!>

Attention ⚠️

This routine only allows to solve equations that are at most of size 4 x 4. The operator arguments TRANSA and TRANSB are fixed to op1(H) = H and op2(H) = H. The primary use of this function is to end the recursion in dla_tgcsylv_dual_recursive.

Okay so this one is literally the "base case" bestie — it exists to stop the recursion from going forever. Iconic, honestly. 🦋

Parameters

SGN1

!>          SGN1 is REAL, allowed values: +/-1
!>          Specifies the sign between in the first equation.
!>

SGN2

!>          SGN2 is REAL, allowed values: +/-1
!>          Specifies the sign between in the second equation.
!>

M

!>          M is INTEGER
!>          The order of the matrices A and C.  128 >= M >= 0.
!>

N

!>          N is INTEGER
!>          The order of the matrices B and D.  128 >= N >= 0.
!>

A

!>          A is REAL array, dimension (LDA,M)
!>          The matrix A must be (quasi-) upper triangular.
!>

LDA

!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,M).
!>

B

!>          B is REAL array, dimension (LDB,N)
!>          The matrix B must be (quasi-) upper triangular. If the matrix D is already
!>          quasi-upper triangular the matrix B must be upper triangular.
!>

LDB

!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>

C

!>          C is REAL array, dimension (LDC,M)
!>          The matrix C must be upper triangular.
!>

LDC

!>          LDC is INTEGER
!>          The leading dimension of the array C.  LDC >= max(1,M).
!>

D

!>          D is REAL array, dimension (LDD,N)
!>          The matrix D must be (quasi-) upper triangular. If the matrix B is already
!>          quasi-upper triangular the matrix D must be upper triangular.
!>

LDD

!>          LDD is INTEGER
!>          The leading dimension of the array D.  LDD >= max(1,N).
!>

E

!>          E is REAL array, dimension (LDE,N)
!>          On input, the matrix E contains the right hand side E.
!>          On output, the matrix E contains the solution R of Equation (1)
!>          as selected by TRANSA, TRANSB, and SGN1/SGN2.
!>          Right hand side E and the solution R are M-by-N matrices.
!>

LDE

!>          LDE is INTEGER
!>          The leading dimension of the array E.  LDE >= max(1,M).
!>

F

!>          F is REAL array, dimension (LDF,N)
!>          On input, the matrix F contains the right hand side F.
!>          On output, the matrix F contains the solution L of Equation (1)
!>          as selected by TRANSA, TRANSB, and SGN1/SGN2.
!>          Right hand side F and the solution L are M-by-N matrices.
!>

LDF

!>          LDF is INTEGER
!>          The leading dimension of the array F.  LDE >= max(1,M).
!>

SCALE

!>          SCALE is REAL
!>          SCALE is a scaling factor to prevent the overflow in the result.
!>          If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems
!>          could not be solved correctly, 0 < SCALE <= 1 holds true.
!>

INFO

!>          INFO is INTEGER
!>
!>          On exit:
!>            == 0:  successful exit
!>            > 0:  The equation is not solved correctly. One of the arising inner
!>                  system got singular.
!>

Author 🖋️

Martin Koehler, MPI Magdeburg

Date 📅

October 2025

Definition at line 206 of file sla_tgcsylv_dual_kernel_44_nn.f90.

🦋 subroutine `sla_tgcsylv_dual_kernel_44nt` 🦋

subroutine sla_tgcsylv_dual_kernel_44nt (real sgn1, real sgn2, integer m,
    integer n, real, dimension(lda, *) a, integer lda, real, dimension(ldb,
    *) b, integer ldb, real, dimension(ldc, *) c, integer ldc, real,
    dimension(ldd, *) d, integer ldd, real, dimension(lde, *) e, integer
    lde, real, dimension(ldf,*) f, integer ldf, real scale, integer info)

Solver for a 4x4 dual generalized coupled Sylvester equation (TRANSA = N, TRANSB = T). Look at her — same 4x4 energy, but this time TRANSB is T. She's the transposed twin and she ate. 💅✨

Purpose

okay bestie, so this whole routine is out here solving a generalized coupled Sylvester equation that looks exactly like this:

   op1(A)**T * R  + op1(C)**T * L               = SCALE * E                             (1)
   SGN1 * R * op2(B)**T + SGN2 * L * op2(D)** T = SCALE * F

where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix, and B and D are N-by-N upper triangular matrices. Also — and this is the little plot twist — either B or D can be quasi upper triangular too. The right hand sides E, F and the solutions R, L are M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) (or (D,B)) get created by SGGES from LAPACK. 💅

Equation (1) is the dual to the generalized coupled Sylvester equation

   op1(A) * R + SGN1 * L * op2(B)  = SCALE * E                                          (2)
   op1(C) * R + SGN2 * L * op2(D)  = SCALE * F

Equation (1) is the dual one to equation (2) with respect to the underlying linear system. Let Z be the matrix formed by rewriting (2) into its Kronecker form. That gives you:

        | kron(I, op1(A))   SGN1*kron(op2(B)**T, I) | | Vec R |   | Vec E |
  Z X = |                                           |*|       | = |       |
        | kron(I, op1(C))   SGN2*kron(op2(D)**T, I) | | Vec L |   | Vec F |

And regarding Z**T you obtain

           | kron(I, op1(A)**T )    kron(I, op1(C)**T)   | | Vec R |   | Vec E |
  Z**T X = |                                             |*|       | = |       |
           | SGN1*kron(op2(B), I)   SGN2*kron(op2(D), I) | | Vec L |   | Vec F |

which belongs to the Sylvester equation (1). For this reason the parameters TRANSA and TRANSB are expressed in terms of the Sylvester equation (2). 🦋

Attention (read this fr) 💀

This routine only allows to solve equations that are at most of size 4 x 4. The operator arguments TRANSA and TRANSB are fixed to op1(H) = H and op2(H) = H**T. The primary use of this function is to end the recursion in dla_tgcsylv_dual_recursive.

Parameters 🧾

SGN1 — SGN1 is REAL, allowed values: +/-1. Specifies the sign between in the first equation. ✨
SGN2 — SGN2 is REAL, allowed values: +/-1. Specifies the sign between in the second equation. ✨
M — M is INTEGER. The order of the matrices A and C. 128 >= M >= 0.
N — N is INTEGER. The order of the matrices B and D. 128 >= N >= 0.
A — A is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular. 👑
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — B is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular the matrix B must be upper triangular. 🫶
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — C is REAL array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular. 🫶
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — E is REAL array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices. 🌸
LDE — LDE is INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F — F is REAL array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices. 🌸
LDF — LDF is INTEGER. The leading dimension of the array F. LDE >= max(1,M). (yes bestie, the source literally says LDE there and not LDF — not me, I'm just the messenger 🙈)
SCALE — SCALE is REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true. 💖
INFO — INFO is INTEGER. On exit:
- == 0: successful exit 🎉
- > 0: The equation is not solved correctly. One of the arising inner system got singular. 😭

Author 👑

Martin Koehler, MPI Magdeburg

Date 📅

October 2025

Definition at line 206 of file sla_tgcsylv_dual_kernel_44_nt.f90.

`sla_tgcsylv_dual_kernel_44tn`

subroutine sla_tgcsylv_dual_kernel_44tn (real sgn1, real sgn2, integer m,
    integer n, real, dimension(lda, *) a, integer lda, real, dimension(ldb,
    *) b, integer ldb, real, dimension(ldc, *) c, integer ldc, real,
    dimension(ldd, *) d, integer ldd, real, dimension(lde, *) e, integer
    lde, real, dimension(ldf,*) f, integer ldf, real scale, integer info)

Solver for a 4x4 dual generalized coupled Sylvester equation (TRANSA = T, TRANSB = N). 🦋

Purpose 🎯

SLA_TGCSYLV_DUAL_KERNEL_44TN solves a generalized coupled Sylvester equation of the following form

   op1(A)**T * R  + op1(C)**T *L               = SCALE * E                             (1)
   SGN1 * R * op2(B)**T + SGN2 * L * op2(D)** T = SCALE * F

where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix and B and D are N-by-N upper triangular matrices. Furthermore either B or D can be quasi upper triangular as well. The right hand sides E, F and the solutions R, L are M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) (or (D,B)) are created by SGGES from LAPACK. The equation (1) is the dual to the generalized coupled Sylvester equation

   op1(A) * R + SGN1 * L * op2(B)  = SCALE * E                                          (2)
   op1(C) * R + SGN2 * L * op2(D)  = SCALE * F

The equation (1) is the dual one to equation (2) with respect to the underlying linear system. Let Z be the matrix formed by rewriting (2) into its Kronecker form. This yields

        | kron(I, op1(A))   SGN1*kron(op2(B)**T, I) | | Vec R |   | Vec E |
  Z X = |                                           |*|       | = |       |
        | kron(I, op1(C))   SGN2*kron(op2(D)**T, I) | | Vec L |   | Vec F |

Regarding Z**T one obtains

           | kron(I, op1(A)**T )    kron(I, op1(C)**T)   | | Vec R |   | Vec E |
  Z**T X = |                                             |*|       | = |       |
           | SGN1*kron(op2(B), I)   SGN2*kron(op2(D), I) | | Vec L |   | Vec F |

which belongs to the Sylvester equation (1). For this reason the parameters TRANSA and TRANSB are expressed in terms of the Sylvester equation (2). ✨

Attention (read this fr) 💀

This routine only allows to solve equations that are at most of size 4 x 4. The operator arguments TRANSA and TRANSB are fixed to op1(H) = H**T and op2(H) = H. The primary use of this function is to end the recursion in dla_tgcsylv_dual_recursive.

Parameters 🧾

SGN1 — SGN1 is REAL, allowed values: +/-1. Specifies the sign between in the first equation. ✨
SGN2 — SGN2 is REAL, allowed values: +/-1. Specifies the sign between in the second equation. ✨
M — M is INTEGER. The order of the matrices A and C. 128 >= M >= 0.
N — N is INTEGER. The order of the matrices B and D. 128 >= N >= 0.
A — A is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular. 👑
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — B is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular the matrix B must be upper triangular. 🫶
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — C is REAL array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular. 🫶
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — E is REAL array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices. 🌸
LDE — LDE is INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F —
F — F is REAL array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices. 💖
LDF — LDF is INTEGER. The leading dimension of the array F. LDE >= max(1,M).

not me noticing it literally says LDE here and not LDF 😭 that's what the page says bestie, we do not touch the sacred text 💅

SCALE — SCALE is REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true. no overflow allowed on this campus 🫶
INFO — INFO is INTEGER.

On exit:

== 0: successful exit ✨
> 0: The equation is not solved correctly. One of the arising inner system got singular. 💀

Author

Martin Koehler, MPI Magdeburg

Date

October 2025

Definition at line 206 of file sla_tgcsylv_dual_kernel_44_tn.f90.

👑 subroutine `sla_tgcsylv_dual_kernel_44tt`

subroutine sla_tgcsylv_dual_kernel_44tt (real sgn1, real sgn2, integer m, integer n, real, dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d, integer ldd, real, dimension(lde, *) e, integer lde, real, dimension(ldf,*) f, integer ldf, real scale, integer info)

Solver for a 4x4 dual generalized coupled Sylvester equation (TRANSA = T, TRANSB = T) 🌈

Purpose

okay so SLA_TGCSYLV_DUAL_KERNEL_44TT solves a generalized coupled Sylvester equation shaped exactly like this 👇

   op1(A)**T * R  + op1(C)**T * L               = SCALE * E                             (1)
   SGN1 * R * op2(B)**T + SGN2 * L * op2(D)** T = SCALE * F

where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix and B and D are N-by-N upper triangular matrices. Furthermore either B or D can be quasi upper triangular as well. The right hand sides E, F and the solutions R, L are M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) ( or (D,B) ) are created by SGGES from LAPACK. ✨

The equation (1) is the dual to the generalized coupled Sylvester equation

   op1(A) * R + SGN1 * L * op2(B)  = SCALE * E                                          (2)
   op1(C) * R + SGN2 * L * op2(D)  = SCALE * F

The equation (1) is the dual one to equation (2) with respect to the underlying linear system. Let Z be the matrix formed by rewriting (2) into its Kronecker form. This yields 🧚

        | kron(I, op1(A))   SGN1*kron(op2(B)**T, I) | | Vec R |   | Vec E |
  Z X = |                                           |*|       | = |       |
        | kron(I, op1(C))   SGN2*kron(op2(D)**T, I) | | Vec L |   | Vec F |

Regarding Z**T one obtains

           | kron(I, op1(A)**T )    kron(I, op1(C)**T)   | | Vec R |   | Vec E |
  Z**T X = |                                             |*|       | = |       |
           | SGN1*kron(op2(B), I)   SGN2*kron(op2(D), I) | | Vec L |   | Vec F |

which belongs to the Sylvester equation (1). For this reason the parameters TRANSA and TRANSB are expressed in terms of the Sylvester equation (2). 🤯

Attention

This routine only allows to solve equations that are at most of size 4 x 4. The operator arguments TRANSA and TRANSB are fixed to op1(H) = H**T and op2(H) = H**T. The primary use of this function is to end the recursion in dla_tgcsylv_dual_recursive. 🚨

Parameters

SGN1 — SGN1 is REAL, allowed values: +/-1. Specifies the sign between in the first equation. ➕➖
SGN2 — SGN2 is REAL, allowed values: +/-1. Specifies the sign between in the second equation. ➕➖
M — M is INTEGER. The order of the matrices A and C. 128 >= M >= 0.
N — N is INTEGER. The order of the matrices B and D. 128 >= N >= 0.
A — A is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular. 🔺
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — B is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular the matrix B must be upper triangular. 🔺
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — C is REAL array, dimension (LDC,M). The matrix C must be upper triangular. 🔺
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular. 🔺
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — E is REAL array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices. 💖
LDE — LDE is INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F — F is REAL array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices. 💖
LDF — LDF is INTEGER. The leading dimension of the array F. LDE >= max(1,M).

same LDE situation as before bestie, it's giving copy-paste energy and honestly? iconic 🙈

SCALE — SCALE is REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true. protecting your floats from the void 🫶
INFO — INFO is INTEGER.

On exit:

== 0: successful exit ✨
> 0: The equation is not solved correctly. One of the arising inner system got singular. 💀

Author

Martin Koehler, MPI Magdeburg

Date

October 2025

Definition at line 206 of file sla_tgcsylv_dual_kernel_44_tt.f90.

🌸 subroutine `sla_tgcsylv_dual_l2`

subroutine sla_tgcsylv_dual_l2 (character, dimension(1) transa, character, dimension(1) transb, real sgn1, real sgn2, integer m, integer n, real, dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d, integer ldd, real, dimension(lde, *) e, integer lde, real, dimension(ldf,*) f, integer ldf, real scale, real, dimension(*) work, integer info)

Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the dual generalized coupled Sylvester equation. 👑

Purpose

SLA_TGCSYLV_DUAL_L2 solves a generalized coupled Sylvester equation that looks like this 👇

   op1(A)**T * R  + op1(C)**T * L               = SCALE * E                             (1)
   SGN1 * R * op2(B)**T + SGN2 * L * op2(D)** T = SCALE * F

where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix and B and D are N-by-N upper triangular matrices. Furthermore either B or D can be quasi upper triangular as well. The right hand sides E, F and the solutions R, L are M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) ( or (D,B) ) are created by SGGES from LAPACK. ✨

The equation (1) is the dual to the generalized coupled Sylvester equation

   op1(A) * R + SGN1 * L * op2(B)  = SCALE * E                                          (2)
   op1(C) * R + SGN2 * L * op2(D)  = SCALE * F

The equation (1) is the dual one to equation (2) with respect to the underlying linear system. Let Z be the matrix formed by rewriting (2) into its Kronecker form. This yields 🧚

        | kron(I, op1(A))   SGN1*kron(op2(B)**T, I) | | Vec R |   | Vec E |
  Z X = |                                           |*|       | = |       |
        | kron(I, op1(C))   SGN2*kron(op2(D)**T, I) | | Vec L |   | Vec F |

Regarding Z**T one obtains

           | kron(I, op1(A)**T )    kron(I, op1(C)**T)   | | Vec R |   | Vec E |
  Z**T X = |                                             |*|       | = |       |
           | SGN1*kron(op2(B), I)   SGN2*kron(op2(D), I) | | Vec L |   | Vec F |

which belongs to the Sylvester equation (1). For this reason the parameters TRANSA and TRANSB are expressed in terms of the Sylvester equation (2). 🤯

Remarks

The algorithm is implemented using BLAS level 2 operations and hand written optimizations. 💅

Parameters

TRANSA — TRANSA is CHARACTER(1). Specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)
TRANSB — TRANSB is CHARACTER(1). Specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)
SGN1 — SGN1 is REAL, allowed values: +/-1. Specifies the sign between in the first equation. ➕➖
SGN2
SGN2 — it's a REAL, bestie, and the only allowed values are +/-1. It specifies the sign between in the second equation. 💅
M — INTEGER. The order of the matrices A and C. M >= 0. 👑
N — INTEGER. The order of the matrices B and D. N >= 0. ✨
A — REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular. 🦋
LDA — INTEGER. The leading dimension of the array A. LDA >= max(1,M). 💖
B — REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular, then the matrix B must be upper triangular. 🌸
LDB — INTEGER. The leading dimension of the array B. LDB >= max(1,N). 💜
C — REAL array, dimension (LDC,M). The matrix C must be upper triangular. 🎀
LDC — INTEGER. The leading dimension of the array C. LDC >= max(1,M). 🫶
D — REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular, then the matrix D must be upper triangular. 💫
LDD — INTEGER. The leading dimension of the array D. LDD >= max(1,N). 🌈
E — REAL array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices. 🔥
LDE — INTEGER. The leading dimension of the array E. LDE >= max(1,M). 😭
F — REAL array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices. 💖
LDF — INTEGER. The leading dimension of the array F. LDE >= max(1,M). (Yes, it really says LDE here — we're keeping it exactly as the original, bestie. 🤷‍♀️)
SCALE — REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true. 🧚
WORK — REAL array, dimension LWORK. Workspace for the algorithm. The workspace needs to queried before the running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO. ⭐
INFO — INTEGER.
- On input:
  - == -1 : Perform a workspace query
  - <> -1: normal operation
- On exit, workspace query:
  - < 0 : if INFO = -i, the i-th argument had an illegal value
  - >= 0: The value of INFO is the required number of elements in the workspace.
- On exit, normal operation:
  - == 0: successful exit
  - < 0: if INFO = -i, the i-th argument had an illegal value
  - > 0: The equation is not solved correctly. One of the arising inner system got singular.

✨ Optimizations ✨

Replaced SAXPY operation by Fortran intrinsic
Replaced all BLAS calls except of SGER by Fortran intrinsic,

Author: Martin Koehler, MPI Magdeburg Date: October 2025 Definition at line 245 of file sla_tgcsylv_dual_l2.f90.

Okay so next up, we have a whole new subroutine: sla_tgcsylv_dual_l2_local_copy! 👑

   subroutine sla_tgcsylv_dual_l2_local_copy (character, dimension(1) transa,
       character, dimension(1) transb, real sgn1, real sgn2, integer m, integer
       n, real, dimension(lda, *) a, integer lda, real, dimension(ldb, *) b,
       integer ldb, real, dimension(ldc, *) c, integer ldc, real,
       dimension(ldd, *) d, integer ldd, real, dimension(lde, *) e, integer
       lde, real, dimension(ldf,*) f, integer ldf, real scale, real,
       dimension(*) work, integer info)

Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the dual generalized coupled Sylvester equation (Local Copy)

Purpose:

...and the whole point of it is to solve a generalized coupled Sylvester equation that looks like this, bestie:

   op1(A)**T * R  + op1(C)**T * L               = SCALE * E                             (1)
   SGN1 * R * op2(B)**T + SGN2 * L * op2(D)** T = SCALE * F

where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix and B and D are N-by-N upper triangular matrices. And plot twist ✨: either B or D can be quasi upper triangular as well. The right hand sides E, F and the solutions R, L are M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) ( or (D,B) ) are created by SGGES from LAPACK. ikyk what that is 💅

The equation (1) is the dual to the generalized coupled Sylvester equation

   op1(A) * R + SGN1 * L * op2(B)  = SCALE * E                                          (2)
   op1(C) * R + SGN2 * L * op2(D)  = SCALE * F

The equation (1) is the dual one to equation (2) with respect to the underlying linear system. Let Z be the matrix formed by rewriting (2) into its Kronecker form. This yields

        | kron(I, op1(A))   SGN1*kron(op2(B)**T, I) | | Vec R |   | Vec E |
  Z X = |                                           |*|       | = |       |
        | kron(I, op1(C))   SGN2*kron(op2(D)**T, I) | | Vec L |   | Vec F |

Regarding Z**T one obtains

           | kron(I, op1(A)**T )    kron(I, op1(C)**T)   | | Vec R |   | Vec E |
  Z**T X = |                                             |*|       | = |       |
           | SGN1*kron(op2(B), I)   SGN2*kron(op2(D), I) | | Vec L |   | Vec F |

which belongs to the Sylvester equation (1). For this reason the parameters TRANSA and TRANSB are expressed in terms of the Sylvester equation (2). no cap, that's the whole reason for the naming scheme 🫶

✨ REMARKS ✨

The algorithm is implemented using BLAS level 2 operations and hand written optimizations and uses local copies of the input data. For this reason M and N are limited to 128. yes literally capped at 128, that's not a vibe, that's the law 📏

✨ PARAMETERS ✨

TRANSA — TRANSA is CHARACTER(1). Specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)
TRANSB — TRANSB is CHARACTER(1). Specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)
SGN1 — SGN1 is REAL, allowed values: +/-1. Specifies the sign between in the first equation. it's just a lil plus-or-minus, she's simple like that 💖
SGN2 — SGN2 is REAL, allowed values: +/-1. Specifies the sign between in the second equation.
M — M is INTEGER. The order of the matrices A and C. 128 >= M >= 0.
N — N is INTEGER. The order of the matrices B and D. 128 >= N >= 0.
A — A is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular. don't try to sneak a random dense matrix in there, she will not be having it 🙈
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — B is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular the matrix B must be upper triangular. read that twice bestie, it's a real constraint ⚠️
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — C is REAL array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular.
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — E is REAL array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices. she goes in as a RHS and comes out as a whole solution, main character energy fr 🎀
LDE — LDE is INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F — F is REAL array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices.
LDF — LDF is INTEGER. The leading dimension of the array F. LDE >= max(1,M).

⚠️ heads up: yes, that last line says LDE and not LDF. That is what the page says, so that is what you get. copied faithfully, no edits, we don't fix other people's typos, we just report them 💅

SCALE — SCALE is REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — WORK is REAL array, dimension LWORK. Workspace for the algorithm. The workspace needs to queried before the running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO. yes you have to ask first, she doesn't just hand out her workspace for free ✨
INFO — INFO is INTEGER.
- On input:
  - == -1 : Perform a workspace query
  - <> -1: normal operation
- On exit, workspace query:
  - < 0 : if INFO = -i, the i-th argument had an illegal value
  - >= 0: The value of INFO is the required number of elements in the workspace.
- On exit, normal operation:
  - == 0: successful exit
  - < 0: if INFO = -i, the i-th argument had an illegal value
  - > 0: The equation is not solved correctly. One of the arising inner system got singular. 💀 yeah that one's the bad one, singular inner system means the tea is spilled and the math said no

✨ OPTIMIZATIONS: ✨

Replaced SAXPY operation by Fortran intrinsic
No BLAS calls.
Input data is copied to local arrays to increase the data locality if used from inside a level-3 solver.

she's got no BLAS calls at all, she's doing it by hand, artisanal, small-batch linear algebra 🧚

✨ AUTHOR ✨

Martin Koehler, MPI Magdeburg

✨ DATE ✨

October 2025

✨ DEFINITION ✨

Definition at line 247 of file sla_tgcsylv_dual_l2_opt_local_copy.f90.

✨ `sla_tgcsylv_dual_l2_local_copy_128` ✨

subroutine sla_tgcsylv_dual_l2_local_copy_128 ( transa,  transb,  sgn1,
    sgn2,  m,  n,  a,  lda,  b,  ldb,  c,  ldc,  d,  ldd,  e,  lde,  f,
    ldf,  scale,  work,  info)

Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the dual generalized coupled Sylvester equation (Local Copy)

okay so meet her sibling 👑 same energy, same local-copy lifestyle, and yes she is also doing the Kagstroem-Westin/Bartels-Stewart thing at level 2.

✨ PURPOSE: ✨

okay so SLA_TGCSYLV_DUAL_L2_LOCAL_COPY_128 solves a generalized coupled Sylvester equation that looks exactly like this, no cap:

op1(A)**T * R  + op1(C)**T * L               = SCALE * E                             (1)
SGN1 * R * op2(B)**T + SGN2 * L * op2(D)** T = SCALE * F

where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix, and B and D are N-by-N upper triangular matrices. And bestie, plot twist: either B or D can be quasi upper triangular as well. The right hand sides E, F and the solutions R, L are all M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) (or (D,B)) are created by SGGES from LAPACK. 💅

Equation (1) is the dual to the generalized coupled Sylvester equation

op1(A) * R + SGN1 * L * op2(B)  = SCALE * E                                          (2)
op1(C) * R + SGN2 * L * op2(D)  = SCALE * F

The equation (1) is the dual one to equation (2) with respect to the underlying linear system. Let Z be the matrix formed by rewriting (2) into its Kronecker form. This yields

        | kron(I, op1(A))   SGN1*kron(op2(B)**T, I) | | Vec R |   | Vec E |
  Z X = |                                           |*|       | = |       |
        | kron(I, op1(C))   SGN2*kron(op2(D)**T, I) | | Vec L |   | Vec F |

Regarding Z**T one obtains

           | kron(I, op1(A)**T )    kron(I, op1(C)**T)   | | Vec R |   | Vec E |
  Z**T X = |                                             |*|       | = |       |
           | SGN1*kron(op2(B), I)   SGN2*kron(op2(D), I) | | Vec L |   | Vec F |

which belongs to the Sylvester equation (1). For this reason the parameters TRANSA and TRANSB are expressed in terms of the Sylvester equation (2). iykyk ✨

🌟 Remarks

The algorithm is implemented using BLAS level 2 operations plus hand written optimizations, and it uses local copies of the input data. Because of that, M and N are limited to 128. Also, the algorithm uses alignment of the local arrays if this is supported by the compiler. At the moment this is used with ifort and xlf but not with gfortran. 🫶

🎀 Parameters

TRANSA — it's a CHARACTER(1). It specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (no transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (transpose A and C)
TRANSB — it's a CHARACTER(1). It specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (no transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (transpose B and D)
SGN1 — it's REAL, and the only allowed values are +/-1. It specifies the sign between in the first equation. 💖
SGN2 — it's REAL, and the only allowed values are +/-1. It specifies the sign between in the second equation.
M — it's INTEGER. The order of the matrices A and C. 128 >= M >= 0.
N — it's INTEGER. The order of the matrices B and D. 128 >= N >= 0.
A — it's a REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular.
LDA — it's INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — it's a REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular, then the matrix B must be upper triangular.
LDB — it's INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — it's a REAL array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — it's INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — it's a REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular, then the matrix D must be upper triangular.
LDD — it's INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — it's a REAL array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. The right hand side E and the solution R are M-by-N matrices. ✨
LDE — it's INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F — it's a REAL array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. The right hand side F and the solution L are M-by-N matrices.
LDF — it's INTEGER. The leading dimension of the array F. LDE >= max(1,M). (yep, that's literally what it says 🫶 we are simply the messenger)
SCALE — it's REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0; otherwise, if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true. 🔥
WORK — it's a REAL array, dimension LWORK. Workspace for the algorithm. The workspace needs to be queried before running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO.
INFO — it's INTEGER. 👑

On input:

== -1 : perform a workspace query
<> -1: normal operation

On exit, workspace query:

< 0 : if INFO = -i, the i-th argument had an illegal value
>= 0: the value of INFO is the required number of elements in the workspace.

On exit, normal operation:

== 0: successful exit
< 0: if INFO = -i, the i-th argument had an illegal value
> 0: the equation is not solved correctly. One of the arising inner system got singular. 💀

🔥 Optimizations:

Replaced SAXPY operation by Fortran intrinsic
No BLAS calls.
Input data is copied to local arrays to increase the data locality if used from inside a level-3 solver.
The local copies are aligned to 64 byte if supported by the compiler.

✍️ Author

Martin Koehler, MPI Magdeburg

📅 Date

October 2025

Definition at line 248 of file sla_tgcsylv_dual_l2_opt_local_copy_128.f90.

subroutine sla_tgcsylv_dual_l2_local_copy_32 ( transa,  transb,  sgn1,
    sgn2,  m,  n,  a,  lda,  b,  ldb,  c,  ldc,  d,  ldd,  e,  lde,  f,
    ldf,  scale,  work,  info)

Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the dual generalized coupled Sylvester equation (Local Copy) 🦋

Purpose:

okay bestie so SLA_TGCSYLV_DUAL_L2_LOCAL_COPY_32 solves a generalized coupled Sylvester equation that looks exactly like this:

   op1(A)**T * R  + op1(C)**T * L               = SCALE * E                             (1)
   SGN1 * R * op2(B)**T + SGN2 * L * op2(D)** T = SCALE * F

Let me set the scene 💅 — A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix, and B and D are N-by-N upper triangular matrices. And plot twist: either B or D can be quasi upper triangular as well. The right hand sides E, F and the solutions R, L are all M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) (or (D,B)) are the ones SGGES from LAPACK hands you. iconic behavior honestly.

So equation (1) is the dual to the generalized coupled Sylvester equation:

   op1(A) * R + SGN1 * L * op2(B)  = SCALE * E                                          (2)
   op1(C) * R + SGN2 * L * op2(D)  = SCALE * F

And it's the dual one to equation (2) with respect to the underlying linear system, fr fr. Let Z be the matrix you get by rewriting (2) into its Kronecker form. That gives you:

         | kron(I, op1(A))   SGN1*kron(op2(B)**T, I) | | Vec R |   | Vec E |
   Z X = |                                           |*|       | = |       |
         | kron(I, op1(C))   SGN2*kron(op2(D)**T, I) | | Vec L |   | Vec F |

And looking at Z**T you get:

            | kron(I, op1(A)**T )    kron(I, op1(C)**T)   | | Vec R |   | Vec E |
   Z**T X = |                                             |*|       | = |       |
            | SGN1*kron(op2(B), I)   SGN2*kron(op2(D), I) | | Vec L |   | Vec F |

which belongs to the Sylvester equation (1). 🎀 That's literally why the parameters TRANSA and TRANSB are expressed in terms of the Sylvester equation (2). iykyk.

✨ Remarks ✨

The algorithm is implemented using BLAS level 2 operations plus hand written optimizations, and it uses local copies of the input data. For this reason M and N are limited to 32. Also, the algorithm uses alignment of the local arrays if this is supported by the compiler. At the moment this is used with ifort and xlf but not with gfortran. 💀

💅 Parameters 💅

TRANSA — is CHARACTER(1). Specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)
TRANSB — is CHARACTER(1). Specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)
SGN1 — is REAL, allowed values: +/-1. Specifies the sign between in the first equation.
SGN2 — is REAL, allowed values: +/-1. Specifies the sign between in the second equation.
M — is INTEGER. The order of the matrices A and C. 32 >= M >= 0.
N — is INTEGER. The order of the matrices B and D. 32 >= N >= 0.
A — is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular.
LDA — is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular the matrix B must be upper triangular.
LDB — is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — is REAL array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — is REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular.
LDD — is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — is REAL array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices.
LDE — is INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F — is REAL array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices.
LDF — is INTEGER. The leading dimension of the array F. LDE >= max(1,M).
SCALE — is REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — is REAL array, dimension LWORK. Workspace for the algorithm. The workspace needs to queried before the running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO.
INFO — is INTEGER.
- On input:
  - == -1 : Perform a workspace query
  - <> -1: normal operation
- On exit, workspace query:
  - < 0 : if INFO = -i, the i-th argument had an illegal value
  - >= 0: The value of INFO is the required number of elements in the workspace.
- On exit, normal operation:
  - == 0: successful exit
  - < 0: if INFO = -i, the i-th argument had an illegal value
  - > 0: The equation is not solved correctly. One of the arising inner system got singular.

🚀 Optimizations: 🚀

Replaced SAXPY operation by Fortran intrinsic
No BLAS calls.
Input data is copied to local arrays to increase the data locality if used from inside a level-3 solver.
The local copies are aligned to 64 byte if supported by the compiler.

👑 Author 👑

Martin Koehler, MPI Magdeburg

📅 Date 📅

October 2025

📍 Definition 📍

At line 247 of file sla_tgcsylv_dual_l2_opt_local_copy_32.f90.

🌈 `sla_tgcsylv_dual_l2_local_copy_64`

   subroutine sla_tgcsylv_dual_l2_local_copy_64 ( transa,  transb,  sgn1,
       sgn2,  m,  n,  a,  lda,  b,  ldb,  c,  ldc,  d,  ldd,  e,  lde,  f,
       ldf,  scale,  work,  info)

Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the dual generalized coupled Sylvester equation (Local Copy) 🦋

Purpose:

Okay so this whole thing is the doc block for SLA_TGCSYLV_DUAL_L2_LOCAL_COPY_64 — lemme finish that equation breakdown for you bestie, we were literally mid-sentence 💅

!> SLA_TGCSYLV_DUAL_L2_LOCAL_COPY_64 solves a generalized coupled  Sylvester equation of the following form
!>
!>    op1(A)**T * R  + op1(C)**T * L               = SCALE * E                             (1)
!>    SGN1 * R * op2(B)**T + SGN2 * L * op2(D)** T = SCALE * F
!>
!> where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular
!> matrix and B and D are N-by-N upper triangular matrices. Furthermore either B or
!> D can be quasi upper triangular as well. The right hand sides E, F and the solutions R, L are
!> M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) ( or (D,B) ) are
!> created by SGGES from LAPACK.
!> The equation (1) is the dual to the generalized coupled Sylvester equation
!>
!>    op1(A) * R + SGN1 * L * op2(B)  = SCALE * E                                          (2)
!>    op1(C) * R + SGN2 * L * op2(D)  = SCALE * F
!>
!> The equation (1) is the dual one to equation (2) with respect to the underlying linear system.
!> Let Z be the matrix formed by rewriting (2) into its Kronecker form. This yields
!>
!>         | kron(I, op1(A))   SGN1*kron(op2(B)**T, I) | | Vec R |   | Vec E |
!>   Z X = |                                           |*|       | = |       |
!>         | kron(I, op1(C))   SGN2*kron(op2(D)**T, I) | | Vec L |   | Vec F |
!>
!> Regarding Z**T one obtains
!>
!>            | kron(I, op1(A)**T )    kron(I, op1(C)**T)   | | Vec R |   | Vec E |
!>   Z**T X = |                                             |*|       | = |       |
!>            | SGN1*kron(op2(B), I)   SGN2*kron(op2(D), I) | | Vec L |   | Vec F |
!>
!> which belongs to the Sylvester equation (1). For this reason the parameters TRANSA and TRANSB
!> are expressed in terms of the Sylvester equation (2).
!>
!>

So basically Z is just equation (2) in its Kron form, and the transpose Z**T is giving... equation (1) ✨ That's the whole "dual" thing, that's why TRANSA and TRANSB get described in terms of equation (2) even though we're out here solving (1). iykyk 🦋

🌸 Remarks 🌸

The algorithm is implemented using BLAS level 2 operations and hand written optimizations, and it uses local copies of the input data. For this reason M and N are limited to 64 😭 (yes, sixty-four, that's the ceiling bestie). Furthermore, the algorithm uses alignment of the local arrays if this is supported by the compiler. At the moment this is used with ifort and xlf but not with gfortran 💀

✨ Parameters ✨

TRANSA — is CHARACTER(1). Specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)
TRANSB — is CHARACTER(1). Specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)
SGN1 — is REAL, allowed values: +/-1. Specifies the sign between in the first equation. 🔥
SGN2 — is REAL, allowed values: +/-1. Specifies the sign between in the second equation. 🔥
M — is INTEGER. The order of the matrices A and C. 64 >= M >= 0.
N — is INTEGER. The order of the matrices B and D. 64 >= N >= 0.
A — is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular.
LDA — is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular the matrix B must be upper triangular. 👑
LDB — is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — is REAL array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — is REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular. 👑
LDD — is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — is REAL array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices. 🫶
LDE — is INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F — is REAL array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices. 🫶
LDF — is INTEGER. The leading dimension of the array F. LDE >= max(1,M). (yes bestie, that says LDE in the original, we don't fix typos in the manpage, we report them 💀)
SCALE — is REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — is REAL array, dimension LWORK. Workspace for the algorithm. The workspace needs to queried before the running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO. 📋
INFO — is INTEGER.
- On input:
  - == -1 : Perform a workspace query
  - <> -1: normal operation
- On exit, workspace query:
  - < 0 : if INFO = -i, the i-th argument had an illegal value
  - >= 0: The value of INFO is the required number of elements in the workspace.
- On exit, normal operation:
  - == 0: successful exit
  - < 0: if INFO = -i, the i-th argument had an illegal value
  - > 0: The equation is not solved correctly. One of the arising inner system got singular. 😭

🎀 Optimizations 🎀

Replaced SAXPY operation by Fortran intrinsic
No BLAS calls.
Input data is copied to local arrays to increase the data locality if used from inside a level-3 solver.
The local copies are aligned to 64 byte if supported by the compiler.

👑 Author 👑

Martin Koehler, MPI Magdeburg

⭐ Date ⭐

October 2025

Definition at line 246 of file sla_tgcsylv_dual_l2_opt_local_copy_64.f90.

okay and now we're meeting the next subroutine, the ninety-six one, bestie, she's coming in hot 🌈

   subroutine sla_tgcsylv_dual_l2_local_copy_96 ( transa,  transb,  sgn1,
       sgn2,  m,  n,  a,  lda,  b,  ldb,  c,  ldc,  d,  ldd,  e,  lde,  f,
       ldf,  scale,  work,  info)

Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the dual generalized coupled Sylvester equation (Local Copy)

Purpose:

💫 SLA_TGCSYLV_DUAL_L2_LOCAL_COPY_96 💫

okay bestie so this one solves a generalized coupled Sylvester equation and it looks like this:

op1(A)**T * R  + op1(C)**T * L               = SCALE * E                             (1)
SGN1 * R * op2(B)**T + SGN2 * L * op2(D)** T = SCALE * F

where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix, and B and D are N-by-N upper triangular matrices. Also — plot twist — either B or D can be quasi upper triangular too 🫶. The right hand sides E, F and the solutions R, L are M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) (or (D,B)) are created by SGGES from LAPACK.

Equation (1) is the dual to the generalized coupled Sylvester equation

op1(A) * R + SGN1 * L * op2(B)  = SCALE * E                                          (2)
op1(C) * R + SGN2 * L * op2(D)  = SCALE * F

and equation (1) is the dual one to equation (2) with respect to the underlying linear system. Let Z be the matrix formed by rewriting (2) into its Kronecker form. That gives you:

        | kron(I, op1(A))   SGN1*kron(op2(B)**T, I) | | Vec R |   | Vec E |
  Z X = |                                           |*|       | = |       |
        | kron(I, op1(C))   SGN2*kron(op2(D)**T, I) | | Vec L |   | Vec F |

and if you look at Z**T you get

           | kron(I, op1(A)**T )    kron(I, op1(C)**T)   | | Vec R |   | Vec E |
  Z**T X = |                                             |*|       | = |       |
           | SGN1*kron(op2(B), I)   SGN2*kron(op2(D), I) | | Vec L |   | Vec F |

which belongs to the Sylvester equation (1). That's literally why the parameters TRANSA and TRANSB are expressed in terms of the Sylvester equation (2) ✨. iykyk.

🌸 Remarks 🌸

The algorithm is implemented using BLAS level 2 operations and hand written optimizations, and it uses local copies of the input data. Because of that, M and N are limited to 96. On top of that, the algorithm uses alignment of the local arrays if this is supported by the compiler. At the moment this is used with ifort and xlf but not with gfortran.

💖 Parameters 💖

TRANSA — TRANSA is CHARACTER(1). Specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)
TRANSB — TRANSB is CHARACTER(1). Specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)
SGN1 — SGN1 is REAL, allowed values: +/-1. Specifies the sign between in the first equation.
SGN2 — SGN2 is REAL, allowed values: +/-1. Specifies the sign between in the second equation.
M — M is INTEGER. The order of the matrices A and C. 96 >= M >= 0.
N — N is INTEGER. The order of the matrices B and D. 96 >= N >= 0.
A — A is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular.
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — B is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular the matrix B must be upper triangular.
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — C is REAL array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular.
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — E is REAL array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices.
LDE — LDE is INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F — F is REAL array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices.
LDF — LDF is INTEGER. The leading dimension of the array F. LDE >= max(1,M).
SCALE — SCALE is REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — WORK is REAL array, dimension LWORK. Workspace for the algorithm. The workspace needs to queried before the running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO.
INFO — INFO is INTEGER.
- On input:
  - == -1 : Perform a workspace query
  - <> -1: normal operation
- On exit, workspace query:
  - < 0 : if INFO = -i, the i-th argument had an illegal value
  - >= 0: The value of INFO is the required number of elements in the workspace.
- On exit, normal operation:
  - == 0: successful exit
  - < 0: if INFO = -i, the i-th argument had an illegal value
  - > 0: The equation is not solved correctly. One of the arising inner system got singular.

🔥 Optimizations: 🔥

Replaced SAXPY operation by Fortran intrinsic
No BLAS calls.
Input data is copied to local arrays to increase the data locality if used from inside a level-3 solver.
The local copies are aligned to 64 byte if supported by the compiler.

🎀 Author 🎀

Martin Koehler, MPI Magdeburg

🌈 Date 🌈

October 2025

Definition at line 246 of file sla_tgcsylv_dual_l2_opt_local_copy_96.f90.

✨ subroutine sla_tgcsylv_dual_l3 ✨

subroutine sla_tgcsylv_dual_l3 (character, dimension(1) transa, character,
    dimension(1) transb, real sgn1, real sgn2, integer m, integer n, real,
    dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer
    ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d,
    integer ldd, real, dimension(lde, *) e, integer lde, real,
    dimension(ldf,*) f, integer ldf, real scale, real, dimension(*) work,
    integer info)

Level-3 Kagstroem-Westin/Bartels-Stewart Algorithm for the dual generalized coupled Sylvester equation.

🎯 Purpose: 🎯

Which is literally the dual of the generalized coupled Sylvester equation

   op1(A) * R + SGN1 * L * op2(B)  = SCALE * E                                          (2)
   op1(C) * R + SGN2 * L * op2(D)  = SCALE * F

and bestie, "dual" here is a technical thing, not a vibe thing 💅. It means equation (1) is the dual of equation (2) with respect to the underlying linear system. Let Z be the matrix you get when you rewrite (2) into its Kronecker form. That gives you

        | kron(I, op1(A))   SGN1*kron(op2(B)**T, I) | | Vec R |   | Vec E |
  Z X = |                                           |*|       | = |       |
        | kron(I, op1(C))   SGN2*kron(op2(D)**T, I) | | Vec L |   | Vec F |

and then if you look at Z**T instead, you get

           | kron(I, op1(A)**T )    kron(I, op1(C)**T)   | | Vec R |   | Vec E |
  Z**T X = |                                             |*|       | = |       |
           | SGN1*kron(op2(B), I)   SGN2*kron(op2(D), I) | | Vec L |   | Vec F |

which — plot twist 🫢 — belongs to the Sylvester equation (1). So that's why the parameters TRANSA and TRANSB are expressed in terms of the Sylvester equation (2). iykyk 🧚

🎀 Parameters 🎀

TRANSA — is CHARACTER(1). It says what form the system of equations takes with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (no transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (transpose A and C)
TRANSB — is CHARACTER(1). It says what form the system of equations takes with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (no transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (transpose B and D)
SGN1 — is REAL, allowed values: +/-1. It specifies the sign between in the first equation. ✨
SGN2 — is REAL, allowed values: +/-1. It specifies the sign between in the second equation. ✨
M — is INTEGER. The order of the matrices A and C. M >= 0.
N — is INTEGER. The order of the matrices B and D. N >= 0.
A — is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular. 👑
LDA — is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular the matrix B must be upper triangular. 🦋
LDB — is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — is REAL array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — is REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular. 🦋
LDD — is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — is REAL array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices. it literally transforms in place, no cap ✨
LDE — is INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F — is REAL array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices. same energy as E, fr fr ✨
LDF — is INTEGER. The leading dimension of the array F. LDE >= max(1,M). 🫶 (yes, bestie, that's exactly what it says — I'm just the messenger)
SCALE — is REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0, otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — is REAL array, dimension LWORK. Workspace for the algorithm. The workspace needs to queried before the running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO. 🔥
INFO — is INTEGER.
- On input:
  - == -1 : Perform a workspace query
  - <> -1: normal operation
- On exit, workspace query:
  - < 0 : if INFO = -i, the i-th argument had an illegal value
  - >= 0: The value of INFO is the required number of elements in the workspace.
- On exit, normal operation:
  - == 0: successful exit 🎉
  - < 0: if INFO = -i, the i-th argument had an illegal value
  - > 0: The equation is not solved correctly. One of the arising inner system got singular. 💀

🌸 Author 🌸

Martin Koehler, MPI Magdeburg

📅 Date 📅

October 2025

Definition at line 236 of file sla_tgcsylv_dual_l3.f90.

🦋 subroutine sla_tgcsylv_dual_l3_2s 🦋

subroutine sla_tgcsylv_dual_l3_2s (character, dimension(1) transa,
     character, dimension(1) transb, real sgn1, real sgn2, integer m, integer
     n, real, dimension(lda, *) a, integer lda, real, dimension(ldb, *) b,
     integer ldb, real, dimension(ldc, *) c, integer ldc, real,
     dimension(ldd, *) d, integer ldd, real, dimension(lde, *) e, integer
     lde, real, dimension(ldf,*) f, integer ldf, real scale, real,
     dimension(*) work, integer info)

Level-3 Kagstroem-Westin/Bartels-Stewart Algorithm with sub blocking for the dual generalized coupled Sylvester equation.

🎯 Purpose 🎯

okay so real quick, bestie, let me tell you what SLA_TGCSYLV_DUAL_L3_2S actually solves 💖 — it takes on a generalized coupled Sylvester equation shaped like this:

op1(A)**T * R  + op1(C)**T * L               = SCALE * E                             (1)
SGN1 * R * op2(B)**T + SGN2 * L * op2(D)** T = SCALE * F

where, and stay with me here:

A is a M-by-M quasi upper triangular matrix 👑
C is a M-by-M upper triangular matrix
B and D are N-by-N upper triangular matrices 🌸
and either B or D can be quasi upper triangular as well, so slay, flexible

The right hand sides E, F and the solutions R, L are all M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) (or (D,B)) are created by SGGES from LAPACK — iykyk 🧚

Equation (1) is the dual to the generalized coupled Sylvester equation:

op1(A) * R + SGN1 * L * op2(B)  = SCALE * E                                          (2)
op1(C) * R + SGN2 * L * op2(D)  = SCALE * F

Like, equation (1) is the dual one to equation (2) with respect to the underlying linear system. Let Z be the matrix formed by rewriting (2) into its Kronecker form. That gives you:

        | kron(I, op1(A))   SGN1*kron(op2(B)**T, I) | | Vec R |   | Vec E |
  Z X = |                                           |*|       | = |       |
        | kron(I, op1(C))   SGN2*kron(op2(D)**T, I) | | Vec L |   | Vec F |

And then, regarding Z**T, you obtain:

           | kron(I, op1(A)**T )    kron(I, op1(C)**T)   | | Vec R |   | Vec E |
  Z**T X = |                                             |*|       | = |       |
           | SGN1*kron(op2(B), I)   SGN2*kron(op2(D), I) | | Vec L |   | Vec F |

which belongs to the Sylvester equation (1). And that's literally why the parameters TRANSA and TRANSB are expressed in terms of the Sylvester equation (2), not (1) — no cap, that's the whole plot twist 🎀

✨ Parameters ✨

TRANSA — TRANSA is CHARACTER(1). It specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (no transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (transpose A and C)
TRANSB — TRANSB is CHARACTER(1). It specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (no transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (transpose B and D)
SGN1 — SGN1 is REAL, allowed values: +/-1. Specifies the sign between in the first equation 💅
SGN2 — SGN2 is REAL, allowed values: +/-1. Specifies the sign between in the second equation 💅
M — M is INTEGER. The order of the matrices A and C. M >= 0.
N — N is INTEGER. The order of the matrices B and D. N >= 0.
A — A is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular. 👑
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — B is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular the matrix B must be upper triangular. (yes, it's a whole trade-off situation 🫶)
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — C is REAL array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular.
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — E is REAL array, dimension (LDE,N).
- On input, the matrix E contains the right hand side E.
- On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2.
- Right hand side E and the solution R are M-by-N matrices. 🦋
LDE — LDE is INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F — F is REAL array, dimension (LDF,N).
- On input, the matrix F contains the right hand side F.
- On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2.
- Right hand side F and the solution L are M-by-N matrices. 🦋
LDF — LDF is INTEGER. The leading dimension of the array F. LDE >= max(1,M).
SCALE — SCALE is REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true. tiny SCALE = the math was fighting for its life 😭
WORK — WORK is REAL array, dimension LWORK. Workspace for the algorithm. The workspace needs to queried before the running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO. (peep the INFO block below, bestie, it's all connected)
INFO — INFO is INTEGER.
- On input:
  - == -1 : Perform a workspace query
  - <> -1: normal operation
- On exit, workspace query:
  - < 0 : if INFO = -i, the i-th argument had an illegal value
  - >= 0: The value of INFO is the required number of elements in the workspace.
- On exit, normal operation:
  - == 0: successful exit ✨
  - < 0: if INFO = -i, the i-th argument had an illegal value
  - > 0: The equation is not solved correctly. One of the arising inner system got singular. 💀

ngl the workspace query thing is iconic — you literally call the routine with INFO == -1, it tells you how much WORK to bring, and then you actually run it. main character energy fr fr 👑

🌈 Author 🌈

Martin Koehler, MPI Magdeburg

🌸 Date 🌸

October 2025

⭐ Definition ⭐

Definition at line 237 of file sla_tgcsylv_dual_l3_2stage.f90.

💜 recursive subroutine sla_tgcsylv_dual_recursive 💜

okay so here's the signature, and yes she's long, she's worth it:

recursive subroutine sla_tgcsylv_dual_recursive (character, dimension(1)
     transa, character, dimension(1) transb, real sgn1, real sgn2, integer m,
     integer n, real, dimension(lda, *) a, integer lda, real, dimension(ldb,
     *) b, integer ldb, real, dimension(ldc, *) c, integer ldc, real,
     dimension(ldd, *) d, integer ldd, real, dimension(lde, *) e, integer
     lde, real, dimension(ldf,*) f, integer ldf, real scale, real,
     dimension(*) work, integer info)

Recursive blocking Algorithm for the dual generalized coupled Sylvester equation. 🌀

Purpose:

okay bestie, we're deep in the SLA_TGCSYLV_DUAL_RECURSIVE zone and the math is giving ✨

SLA_TGCSYLV_DUAL_RECURSIVE solves a generalized coupled Sylvester equation shaped like this 👇

   op1(A)**T * R  + op1(C)**T * L               = SCALE * E                             (1)
   SGN1 * R * op2(B)**T + SGN2 * L * op2(D)** T = SCALE * F

okay so let's break that down, because it's a lot 💅

A is a M-by-M quasi upper triangular matrix
C is a M-by-M upper triangular matrix
B and D are N-by-N upper triangular matrices
and either B or D can be quasi upper triangular as well (not both being weird at the same time, iykyk)
the right hand sides E, F and the solutions R, L are all M-by-N matrices
typically the matrix pencils (A,C) and (B,D) (or (D,B)) are created by SGGES from LAPACK 🫶

Equation (1) is the dual to the generalized coupled Sylvester equation 🌈

   op1(A) * R + SGN1 * L * op2(B)  = SCALE * E                                          (2)
   op1(C) * R + SGN2 * L * op2(D)  = SCALE * F

Equation (1) is the dual one to equation (2) with respect to the underlying linear system. Let Z be the matrix formed by rewriting (2) into its Kronecker form. This yields 👇

        | kron(I, op1(A))   SGN1*kron(op2(B)**T, I) | | Vec R |   | Vec E |
  Z X = |                                           |*|       | = |       |
        | kron(I, op1(C))   SGN2*kron(op2(D)**T, I) | | Vec L |   | Vec F |

Regarding Z**T one obtains 👇

           | kron(I, op1(A)**T )    kron(I, op1(C)**T)   | | Vec R |   | Vec E |
  Z**T X = |                                             |*|       | = |       |
           | SGN1*kron(op2(B), I)   SGN2*kron(op2(D), I) | | Vec L |   | Vec F |

which belongs to the Sylvester equation (1). For this reason the parameters TRANSA and TRANSB are expressed in terms of the Sylvester equation (2). mind = blown 🤯

TRANSA — TRANSA is CHARACTER(1). It specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)
TRANSB — TRANSB is CHARACTER(1). It specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)
SGN1 — SGN1 is REAL, allowed values: +/-1. Specifies the sign between in the first equation. ✨
SGN2 — SGN2 is REAL, allowed values: +/-1. Specifies the sign between in the second equation. ✨
M — M is INTEGER. The order of the matrices A and C. M >= 0. 💖
N — N is INTEGER. The order of the matrices B and D. N >= 0. 💖
A — A is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular.
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — B is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular the matrix B must be upper triangular. 👀
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — C is REAL array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular. 👀
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — E is REAL array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices. 🌸
LDE — LDE is INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F — F is REAL array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices. 🌸
LDF — LDF is INTEGER. The leading dimension of the array F. LDE >= max(1,M).

ngl, that's literally what the page says — it talks about LDE in the LDF entry 🙈 keeping it exactly as written, no cap.

SCALE — SCALE is REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true. 🔥
WORK — WORK is REAL array, dimension LWORK. Workspace for the algorithm. The workspace needs to queried before the running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO. 🛠️
INFO — INFO is INTEGER.
- On input:
  - == -1 : Perform a workspace query
  - <> -1: normal operation
- On exit, workspace query:
  - < 0 : if INFO = -i, the i-th argument had an illegal value
  - >= 0: The value of INFO is the required number of elements in the workspace.
- On exit, normal operation:
  - == 0: successful exit 🎉
  - < 0: if INFO = -i, the i-th argument had an illegal value
  - > 0: The equation is not solved correctly. One of the arising inner system got singular. 💀

👑 Author

Martin Koehler, MPI Magdeburg

🗓️ Date

October 2025

📍 Definition

Definition at line 240 of file sla_tgcsylv_dual_recursive.f90.

and now, plot twist, a whole other subroutine walks in 💅

subroutine sla_tgcsylv_l2 (character, dimension(1) transa, character,
    dimension(1) transb, real sgn1, real sgn2, integer m, integer n, real,
    dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer
    ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d,
    integer ldd, real, dimension(lde, *) e, integer lde, real,
    dimension(ldf,*) f, integer ldf, real scale, real, dimension(*) work,
    integer info)

Level-2 Kagstroem-Westin/Bartels Algorithm for the generalized coupled Sylvester equation (Optimized) ✨

Purpose:

Remarks

The algorithm is implemented using BLAS level 2 operations and hand written optimizations. 💅✨

Parameters

TRANSA — TRANSA is CHARACTER(1). It specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (No transpose for A and C) 🎀
- == 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C) 🔥
TRANSB — TRANSB is CHARACTER(1). It specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (No transpose for B and D) 🌸
- == 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D) 💫
SGN1 — SGN1 is REAL, allowed values: +/-1. Specifies the sign between in the first equation. (Yes, "between in" is literally what it says, bestie. We keep it exactly.) 🫶
SGN2 — SGN2 is REAL, allowed values: +/-1. Specifies the sign between in the second equation. (Same energy.) 💖
M — M is INTEGER. The order of the matrices A and C. M >= 0. 👑
N — N is INTEGER. The order of the matrices B and D. N >= 0. 👑
A — A is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular. (Don't try to sneak in a lower triangular, it won't slay.) 💅
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M). 📏
B — B is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular the matrix B must be upper triangular. (So like, one of them can be quasi, but if D is quasi then B has to be strictly upper triangular. The drama.) 🎭
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N). 📐
C — C is REAL array, dimension (LDC,M). The matrix C must be upper triangular. (No quasi here, bestie. Just upper triangular.) 🦋
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M). 📏
D — D is REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular. (Same energy as B, but for D.) 🎀
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N). 📐
E — E is REAL array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices. (It goes in as E, comes out as R. Glow up.) ✨
LDE — LDE is INTEGER. The leading dimension of the array E. LDE >= max(1,M). 📏
F — F is REAL array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices. 🌈
LDF — LDF is INTEGER. The leading dimension of the array F. LDE >= max(1,M). (Wait, it says LDE >= max(1,M) for LDF? That's a typo in the original, bestie. We keep it exactly.) 💀
SCALE — SCALE is REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true. 🛟
WORK — WORK is REAL array, dimension LWORK. Workspace for the algorithm. The workspace needs to queried before the running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO. 🧰
INFO — INFO is INTEGER. On input:
- == -1 : Perform a workspace query 🔍
- <> -1: normal operation 🛠️ On exit, workspace query:
- < 0 : if INFO = -i, the i-th argument had an illegal value 😱
- >= 0: The value of INFO is the required number of elements in the workspace. 📦 On exit, normal operation:
- == 0: successful exit 🎉
- < 0: if INFO = -i, the i-th argument had an illegal value 😭
- > 0: The equation is not solved correctly. One of the arising inner system got singular. 💔

Optimizations

Replaced level-3 by level-2 and level-1 calls,
Replaced SAXPY operation by Fortran intrinsic
Replaced all BLAS calls except of SGER by Fortran intrinsic,

Author

Martin Koehler, MPI Magdeburg

Date

October 2025

Definition

Definition at line 225 of file sla_tgcsylv_l2.f90.

subroutine sla_tgcsylv_l2_local_copy

subroutine sla_tgcsylv_l2_local_copy (character, dimension(1) transa,
     character, dimension(1) transb, real sgn1, real sgn2, integer m, integer
     n, real, dimension(lda, *) a, integer lda, real, dimension(ldb, *) b,
     integer ldb, real, dimension(ldc, *) c, integer ldc, real,
     dimension(ldd, *) d, integer ldd, real, dimension(lde, *) e, integer
     lde, real, dimension(ldf,*) f, integer ldf, real scale, real,
     dimension(*) work, integer info)

Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation (Optimized) 🚀

Purpose:

okay so SLA_TGCSYLV_L2_LOCAL_COPY solves a generalized coupled Sylvester equation, and it looks exactly like this:

   op1(A) * R  + SGN1 * L  * op2(B) = SCALE * E                              (1)
   op1(C) * R  + SGN2 * L  * op2(D) = SCALE * F

where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix, and B and D are N-by-N upper triangular matrices. And not only that — either B or D can be quasi upper triangular as well. The right hand side Y and the solution X are M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) (or (D,B)) are created by SGGES from LAPACK. bestie the LAPACK pipeline is real 💅

Remarks

The algorithm is implemented using BLAS level 2 operations and hand written optimizations. no cap, hand written, iconic behavior ✨

Parameters

TRANSA — TRANSA is CHARACTER(1). It specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)
TRANSB — TRANSB is CHARACTER(1). It specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)
SGN1 — SGN1 is REAL, allowed values: +/-1. It specifies the sign between in the first equation.
SGN2 — SGN2 is REAL, allowed values: +/-1. It specifies the sign between in the second equation.
M — M is INTEGER. The order of the matrices A and C. M >= 0.
N — N is INTEGER. The order of the matrices B and D. N >= 0.
A — A is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular.
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — B is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular the matrix B must be upper triangular.
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — C is REAL array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular.
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — E is REAL array, dimension (LDE,N).
- On input, the matrix E contains the right hand side E.
- On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2.
- On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. (yes bestie the doc literally says this twice — we are not editing, we are ✨preserving✨ 🫶)
- Right hand side E and the solution R are M-by-N matrices.
LDE — LDE is INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F — F is REAL array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices.
LDF — LDF is INTEGER. The leading dimension of the array F. LDE >= max(1,M). (that is genuinely what it says — LDE, not LDF. we are just the messengers 🙈)
SCALE — SCALE is REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — WORK is REAL array, dimension LWORK. Workspace for the algorithm. The workspace needs to queried before the running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO.
INFO — INFO is INTEGER.
- On input:
  - == -1 : Perform a workspace query
  - <> -1: normal operation
- On exit, workspace query:
  - < 0 : if INFO = -i, the i-th argument had an illegal value
  - >= 0: The value of INFO is the required number of elements in the workspace.
- On exit, normal operation:
  - == 0: successful exit
  - < 0: if INFO = -i, the i-th argument had an illegal value
  - > 0: The equation is not solved correctly. One of the arising inner system got singular.

Optimizations:

Replaced level-3 by level-2 and level-1 calls,
Reorder the solution order to column first style,
Replaced SAXPY operation by Fortran intrinsic
Replaced all BLAS calls except of STRMV and SGER by Fortran intrinsic,
Use local copies of A, B, C, D, E, and F (M, N <=128) .

Author

Martin Koehler, MPI Magdeburg

Date

October 2025

Definition at line 227 of file sla_tgcsylv_l2_opt_local_copy.f90.

🦋 `sla_tgcsylv_l2_local_copy_128` 🦋

subroutine sla_tgcsylv_l2_local_copy_128 ( transa,  transb,  sgn1,  sgn2,
    m,  n,  a,  lda,  b,  ldb,  c,  ldc,  d,  ldd,  e,  lde,  f,  ldf,
    scale,  work,  info)

Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation (Optimized)

Purpose:

okay so this one? this one is the _128 girlie. she does SLA_TGCSYLV_L2_LOCAL_COPY_128 and she solves a generalized coupled Sylvester equation that looks like this: ✨

op1(A) * R  + SGN1 * L  * op2(B) = SCALE * E                              (1)
op1(C) * R  + SGN2 * L  * op2(D) = SCALE * F

where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix, and B and D are N-by-N upper triangular matrices. also, ngl, either B or D can be quasi upper triangular as well. the right hand side is Y and the solution is X, both M-by-N matrices. typically the matrix pencils (A,C) and (B,D) (or (D,B)) are made by SGGES from LAPACK. 💅

🌸 Remarks 🌸

the algorithm is implemented using BLAS level 2 operations plus hand written optimizations and aligned local copies. periodt.

✨ Parameters ✨

here comes the whole cast, bestie. every single one:

TRANSA — TRANSA is CHARACTER(1). it specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (no transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (transpose A and C)
TRANSB — TRANSB is CHARACTER(1). it specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (no transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (transpose B and D)
SGN1 — SGN1 is REAL, allowed values: +/-1. it specifies the sign between in the first equation.
SGN2 — SGN2 is REAL, allowed values: +/-1. it specifies the sign between in the second equation.
M — M is INTEGER. the order of the matrices A and C. 128 => M >= 0.
N — N is INTEGER. the order of the matrices B and D. 128 >= N >= 0.
A — A is REAL array, dimension (LDA,M). the matrix A must be (quasi-) upper triangular.
LDA — LDA is INTEGER. the leading dimension of the array A. LDA >= max(1,M).
B — B is REAL array, dimension (LDB,N). the matrix B must be (quasi-) upper triangular. if the matrix D is already quasi-upper triangular the matrix B must be upper triangular.
LDB — LDB is INTEGER. the leading dimension of the array B. LDB >= max(1,N).
C — C is REAL array, dimension (LDC,M). the matrix C must be upper triangular.
LDC — LDC is INTEGER. the leading dimension of the array C. LDC >= max(1,M).
D — D is REAL array, dimension (LDD,N). the matrix D must be (quasi-) upper triangular. if the matrix B is already quasi-upper triangular the matrix D must be upper triangular.
LDD — LDD is INTEGER. the leading dimension of the array D. LDD >= max(1,N).
E — E is REAL array, dimension (LDE,N). on input, the matrix E contains the right hand side E. on output, the matrix E contains the solution R of Equation (1) — (yeah bestie it says it twice in the source, so it says it twice here too, we don't gatekeep) — on output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. right hand side E and the solution R are M-by-N matrices.
LDE — LDE is INTEGER. the leading dimension of the array E. LDE >= max(1,M).
F — F is REAL array, dimension (LDF,N). on input, the matrix F contains the right hand side F. on output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. right hand side F and the solution L are M-by-N matrices.
LDF — LDF is INTEGER. the leading dimension of the array F. LDE >= max(1,M). (yes. it says LDE here. in the LDF entry. iconic of it honestly.)
SCALE — SCALE is REAL. SCALE is a scaling factor to prevent the overflow in the result. if INFO == 0 then SCALE is 1.0, otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — WORK is REAL array, dimension LWORK. workspace for the algorithm. the workspace needs to be queried before running the computation. the query is performed by calling the subroutine with INFO == -1 on input. the required workspace is then returned in INFO.
INFO — INFO is INTEGER.
- on input:
  - == -1 : perform a workspace query
  - <> -1: normal operation
- on exit, workspace query:
  - < 0 : if INFO = -i, the i-th argument had an illegal value
  - >= 0: the value of INFO is the required number of elements in the workspace.
- on exit, normal operation:
  - == 0: successful exit
  - < 0: if INFO = -i, the i-th argument had an illegal value
  - > 0: the equation is not solved correctly. one of the arising inner system got singular. (that's the "we tried, bestie" exit 💀)

👑 Optimizations: 👑

the glow-up list. she ate:

replaced level-3 by level-2 and level-1 calls,
reorder the solution order to column first style,
replaced SAXPY operation by Fortran intrinsic
replaced all BLAS calls except of STRMV and SGER by Fortran intrinsic,
align local copies and fix problem size to <= 128

🦋 Author 🦋

Martin Koehler, MPI Magdeburg 🫶

🎀 Date 🎀

October 2025 ✨

Definition at line 225 of file sla_tgcsylv_l2_opt_local_copy_128.f90.

and now, entering the chat: the _32 version. same energy, smaller alignment vibe 💅

subroutine sla_tgcsylv_l2_local_copy_32 ( transa,  transb,  sgn1,  sgn2,  m,
     n,  a,  lda,  b,  ldb,  c,  ldc,  d,  ldd,  e,  lde,  f,  ldf,  scale,
     work,  info)

Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation (Optimized)

Purpose:

WORK — WORK is a REAL array, dimension LWORK. It's the workspace for the algorithm, bestie, and you gotta query it before you run the computation. How? Call the subroutine with INFO == -1 on input, and the required workspace gets handed back to you in INFO. It's giving "ask before you order" energy. 💅
INFO — INFO is INTEGER. She's got two modes and you need to know both, no cap:

On input:

== -1 : perform a workspace query 🕵️
<> -1 : normal operation

On exit, workspace query:

< 0 : if INFO = -i, the i-th argument had an illegal value 😭
>= 0 : the value of INFO is the required number of elements in the workspace ✨

On exit, normal operation:

== 0 : successful exit, slay 💖
< 0 : if INFO = -i, the i-th argument had an illegal value
> 0 : the equation is not solved correctly. One of the arising inner system got singular 💀

✨ Optimizations: ✨

Replaced level-3 by level-2 and level-1 calls,
Reorder the solution order to column first style,
Replaced SAXPY operation by Fortran intrinsic
Replaced all BLAS calls except of STRMV and SGER by Fortran intrinsic,
Align local copies and fix problem size to <= 32

👑 Author

Martin Koehler, MPI Magdeburg

🗓️ Date

October 2025

📍 Definition

At line 226 of file sla_tgcsylv_l2_opt_local_copy_32.f90.

subroutine sla_tgcsylv_l2_local_copy_64 ( transa,  transb,  sgn1,  sgn2,  m,
    n,  a,  lda,  b,  ldb,  c,  ldc,  d,  ldd,  e,  lde,  f,  ldf,  scale,
    work,  info)

Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation (Optimized) 🦋

🌸 Purpose:

TRANSA — this one's CHARACTER(1), bestie, and it decides what shape the equations take on the A-and-C side 💅
- 'N' → op1(A) = A, op1(C) = C (no transpose for A and C, we're keeping it natural ✨)
- 'T' → op1(A) = A**T, op1(C) = C **T (transpose A and C, flip it and reverse it 🦋)
TRANSB — also CHARACTER(1), and it does the exact same job but for B and D 🎀
- 'N' → op2(B) = B, op2(D) = D (no transpose for B and D)
- 'T' → op2(B) = B**T, op2(D) = D **T (transpose B and D)
SGN1 — REAL, and the only allowed values are +/-1. It specifies the sign between in the first equation 🔥
SGN2 — REAL, allowed values +/-1 again. It specifies the sign between in the second equation 🔥
M — INTEGER. The order of the matrices A and C. 64 => M >= 0 👑
N — INTEGER. The order of the matrices B and D. 64 >= N >= 0 🌸
A — REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular, no exceptions, that's the whole vibe 💖
LDA — INTEGER. The leading dimension of the array A. LDA >= max(1,M) ⭐
B — REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. And listen: if the matrix D is already quasi-upper triangular, then the matrix B must be upper triangular. That's the deal, okay? 🫶
LDB — INTEGER. The leading dimension of the array B. LDB >= max(1,N) ⭐
C — REAL array, dimension (LDC,M). The matrix C must be upper triangular 💜
LDC — INTEGER. The leading dimension of the array C. LDC >= max(1,M) ⭐
D — REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. And same energy as B: if the matrix B is already quasi-upper triangular, then the matrix D must be upper triangular 🌈
LDD — INTEGER. The leading dimension of the array D. LDD >= max(1,N) ⭐
E — REAL array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1). On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices. It literally swaps roles on you, main character behavior 😭
LDE — INTEGER. The leading dimension of the array E. LDE >= max(1,M) ⭐
F — REAL array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices 💖
LDF — INTEGER. The leading dimension of the array F. LDE >= max(1,M) ⭐
SCALE — REAL. It's a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0; otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true. So if it's not 1.0, something got messy 💀
WORK — REAL array, dimension LWORK. Workspace for the algorithm. And no cap: the workspace needs to be queried before running the computation. You do that query by calling the subroutine with INFO == -1 on input, and the required workspace is then returned in INFO. Do the query, bestie, don't skip it 🙈
INFO — INTEGER

On input:

== -1 → perform a workspace query
<> -1 → normal operation

On exit, workspace query:

< 0 → if INFO = -i, the i-th argument had an illegal value
>= 0 → the value of INFO is the required number of elements in the workspace

On exit, normal operation:

== 0 → successful exit ✨
< 0 → if INFO = -i, the i-th argument had an illegal value
> 0 → the equation is not solved correctly. One of the arising inner system got singular. 💀

✨ Optimizations: ✨

Replaced level-3 by level-2 and level-1 calls 🔥
Reorder the solution order to column first style 🎀
Replaced SAXPY operation by Fortran intrinsic ⚡
Replaced all BLAS calls except of STRMV and SGER by Fortran intrinsic 🌸
Align local copies and fix problem size to <= 64 👑

🌟 Author 🌟

Martin Koehler, MPI Magdeburg 🫶

📅 Date 📅

October 2025

Definition at line 226 of file sla_tgcsylv_l2_opt_local_copy_64.f90.

subroutine sla_tgcsylv_l2_local_copy_96 ( transa,  transb,  sgn1,  sgn2,  m,
    n,  a,  lda,  b,  ldb,  c,  ldc,  d,  ldd,  e,  lde,  f,  ldf,  scale,
    work,  info)

Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation (Optimized) 💅

Purpose:

SLA_TGCSYLV_L2_LOCAL_COPY_96</p> solves a generalized coupled Sylvester equation that looks like this, bestie:

op1(A) * R  + SGN1 * L  * op2(B) = SCALE * E                              (1)
op1(C) * R  + SGN2 * L  * op2(D) = SCALE * F

Where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix, and B and D are N-by-N upper triangular matrices. Also — plot twist — either B or D can be quasi upper triangular as well. The right hand side Y and the solution X are M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) ( or (D,B) ) are created by SGGES from LAPACK. ✨

💅 Remarks 💅

The algorithm is implemented using BLAS level 2 operations and hand written optimizations and aligned local copies. (Yes, she hand-wrote her own optimizations. Iconic behaviour.)

🎀 Parameters 🎀

TRANSA — TRANSA is CHARACTER(1). It specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)
TRANSB — TRANSB is CHARACTER(1). It specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)
SGN1 — SGN1 is REAL, allowed values: +/-1. Specifies the sign between in the first equation. 💫
SGN2 — SGN2 is REAL, allowed values: +/-1. Specifies the sign between in the second equation. 💫
M — M is INTEGER. The order of the matrices A and C. 96 => M >= 0.
N — N is INTEGER. The order of the matrices B and D. 96 >= N >= 0.
A — A is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular.
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — B is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular the matrix B must be upper triangular. (Yeah, only one of them gets to be quasi. The universe said balance. ⚖️)
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — C is REAL array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi upper triangular the matrix D must be upper triangular.
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — E is REAL array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices. 🔁
LDE — LDE is INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F — F is REAL array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices. 🔁
LDF — LDF is INTEGER. The leading dimension of the array F. LDE >= max(1,M).
SCALE — SCALE is REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — WORK is REAL array, dimension LWORK. Workspace for the algorithm. The workspace needs to queried before the running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO. 👑
INFO — INFO is INTEGER.
- On input:
  - == -1 : Perform a workspace query
  - <> -1: normal operation
- On exit, workspace query:
  - < 0 : if INFO = -i, the i-th argument had an illegal value
  - >= 0: The value of INFO is the required number of elements in the workspace.
- On exit, normal operation:
  - == 0: successful exit
  - < 0: if INFO = -i, the i-th argument had an illegal value
  - > 0: The equation is not solved correctly. One of the arising inner system got singular.

bestie if you get INFO > 0 the vibes are rancid — an inner system went singular and the equation did not get solved correctly. Not a slay. 😭

🔥 Optimizations 🔥

Replaced level-3 by level-2 and level-1 calls,
Reorder the solution order to column first style,
Replaced SAXPY operation by Fortran intrinsic
Replaced all BLAS calls except of STRMV and SGER by Fortran intrinsic,
Align local copies and fix problem size to <= 96

Author — Martin Koehler, MPI Magdeburg 🫶

Date — October 2025

Definition at line 225 of file sla_tgcsylv_l2_opt_local_copy_96.f90.

🌈 and now… the next one 🌈

subroutine sla_tgcsylv_l2_reorder (character, dimension(1) transa,
    character, dimension(1) transb, real sgn1, real sgn2, integer m, integer
    n, real, dimension(lda, *) a, integer lda, real, dimension(ldb, *) b,
    integer ldb, real, dimension(ldc, *) c, integer ldc, real,
    dimension(ldd, *) d, integer ldd, real, dimension(lde, *) e, integer
    lde, real, dimension(ldf,*) f, integer ldf, real scale, real,
    dimension(*) work, integer info)

Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation (Optimized)

🦋 Purpose 🦋

okay so bestie, let's talk about what SLA_TGCSYLV_L2_REORDER actually does 💅 — it solves a generalized coupled Sylvester equation that looks like this:

   op1(A) * R  + SGN1 * L  * op2(B) = SCALE * E                              (1)
   op1(C) * R  + SGN2 * L  * op2(D) = SCALE * F

where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix and B and D are N-by-N upper triangular matrices. Furthermore either B or D can be quasi upper triangular as well. The right hand side Y and the solution X M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) ( or (D,B) ) are created by SGGES from LAPACK.

(yeah the original literally says "the right hand side Y and the solution X M-by-N matrices" — we're keeping it exactly like that, no notes 😭)

🎀 Remarks

The algorithm is implemented using BLAS level 2 operations and hand written optimizations. bestie said "I'll do it myself" and honestly? she ate 🔥

✨ Parameters

TRANSA — is CHARACTER(1). Specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)
TRANSB — is CHARACTER(1). Specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)
SGN1 — is REAL, allowed values: +/-1. Specifies the sign between in the first equation.
SGN2 — is REAL, allowed values: +/-1. Specifies the sign between in the second equation.
M — is INTEGER. The order of the matrices A and C. M >= 0.
N — is INTEGER. The order of the matrices B and D. N >= 0.
A — is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular.
LDA — is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular the matrix B must be upper triangular.
LDB — is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — is REAL array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — is REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular.
LDD — is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — is REAL array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices.
LDE — is INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F — is REAL array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices.
LDF — is INTEGER. The leading dimension of the array F. LDE >= max(1,M).

heads up 🫶 — yes, that last one really does say LDE >= max(1,M) in the source, not LDF. we are not fixing it, we are reporting it, periodt.

SCALE — is REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — is REAL array, dimension LWORK. Workspace for the algorithm. The workspace needs to queried before the running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO.
INFO — is INTEGER.
- On input:
  - == -1 : Perform a workspace query
  - <> -1: normal operation
- On exit, workspace query:
  - < 0 : if INFO = -i, the i-th argument had an illegal value
  - >= 0: The value of INFO is the required number of elements in the workspace.
- On exit, normal operation:
  - == 0: successful exit
  - < 0: if INFO = -i, the i-th argument had an illegal value
  - > 0: The equation is not solved correctly. One of the arising inner system got singular.

🔥 Optimizations

Replaced level-3 by level-2 and level-1 calls,
Reorder the solution order to column first style,
Replaced SAXPY operation by Fortran intrinsic
Replaced all BLAS calls except of STRMV and SGER by Fortran intrinsic,
Column first computation of the solution.

look at her go, hand-optimizing like the main character she is 👑

👑 Author

Martin Koehler, MPI Magdeburg

📅 Date

October 2025

Definition at line 226 of file sla_tgcsylv_l2_opt_reorder.f90.

🌸 `sla_tgcsylv_l2_unopt`

and now for the unoptimized sister version, because every queen needs a chill twin 🦋

subroutine sla_tgcsylv_l2_unopt (character, dimension(1) transa, character,
    dimension(1) transb, real sgn1, real sgn2, integer m, integer n, real,
    dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer
    ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d,
    integer ldd, real, dimension(lde, *) e, integer lde, real,
    dimension(ldf,*) f, integer ldf, real scale, real, dimension(*) work,
    integer info)

Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation (unoptimized)

Purpose:

...and the whole thing is basically: it solves the system without any of the fancy optimizations. 💀

🫶 Attention 🫶

The algorithm is implemented using BLAS level 2 operations without further optimizations. For a faster implementation see SLA_TGSYLV_L2. Translation: this one is the slow-and-steady version, bestie. It works! It's just not the glow-up version. ✨

🎀 Parameters 🎀

TRANSA — TRANSA is CHARACTER(1). It says what form the system of equations takes with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (no transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C**T (transpose A and C)
TRANSB — TRANSB is CHARACTER(1). Same vibe check but for B and D:
- == 'N': op2(B) = B, op2(D) = D (no transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D**T (transpose B and D)
SGN1 — SGN1 is REAL, allowed values: +/-1. This is the sign sitting between the two terms in the first equation. 🫡
SGN2 — SGN2 is REAL, allowed values: +/-1. This is the sign sitting between the two terms in the second equation. 🫡
M — M is INTEGER. The order of the matrices A and C. M >= 0. 📏
N — N is INTEGER. The order of the matrices B and D. N >= 0. 📏
A — A is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular. No cap. 🦋
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — B is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular, then the matrix B must be upper triangular. It's giving "one of y'all has to be the chill one."
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — C is REAL array, dimension (LDC,M). The matrix C must be upper triangular. Periodt. 💅
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular, then the matrix D must be upper triangular. Same energy as B, mirrored. Iconic.
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — E is REAL array, dimension (LDE,N).
- On input, the matrix E contains the right hand side E.
- On output, the matrix E contains the solution R of Equation (1), as selected by TRANSA, TRANSB, and SGN1/SGN2.
- Right hand side E and the solution R are M-by-N matrices. Glow-up moment: it goes in as a right hand side and comes out as a whole solution. 🌈
LDE — LDE is INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F — F is REAL array, dimension (LDF,N).
- On input, the matrix F contains the right hand side F.
- On output, the matrix F contains the solution L of Equation (1), as selected by TRANSA, TRANSB, and SGN1/SGN2.
- Right hand side F and the solution L are M-by-N matrices. Same transformation arc, different letter. 🦋
LDF — LDF is INTEGER. The leading dimension of the array F. LDE >= max(1,M). (yes, that's what it says — LDE — we are reporting it exactly, not fixing it, that's the rule bestie 😭)
SCALE — SCALE is REAL. It's a scaling factor to prevent overflow in the result. If INFO == 0 then SCALE is 1.0; otherwise, if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true. So if SCALE isn't 1, something in the inside got a little wobbly. 🙈
WORK — WORK is REAL array, dimension LWORK. Workspace for the algorithm. The workspace needs to be queried before running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO. Ask before you take — always query first, bestie. ✨
INFO — INFO is INTEGER.

On input:

== -1: perform a workspace query
<> -1: normal operation

On exit, workspace query:

< 0: if INFO = -i, the i-th argument had an illegal value
>= 0: the value of INFO is the required number of elements in the workspace

On exit, normal operation:

== 0: successful exit 🎉
< 0: if INFO = -i, the i-th argument had an illegal value
> 0: the equation is not solved correctly. One of the arising inner systems got singular. Not the singular inner system. Not today. 💀

✨ Optimizations: ✨

Nothing. 🫠

That's it. That's the optimizations section. "Nothing." A whole bullet point for it. Slay in the most cursed way possible.

👑 Author 👑

Martin Koehler, MPI Magdeburg

📅 Date 📅

October 2025

📍 Definition 📍

At line 225 of file sla_tgcsylv_l2_unopt.f90.

🦋 And now... presenting... `sla_tgcsylv_l3` 🦋

subroutine sla_tgcsylv_l3 (character, dimension(1) transa, character,
    dimension(1) transb, real sgn1, real sgn2, integer m, integer n, real,
    dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer
    ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d,
    integer ldd, real, dimension(lde, *) e, integer lde, real,
    dimension(ldf,*) f, integer ldf, real scale, real, dimension(*) work,
    integer info)

Level-3 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation (optimized) 🔥

Yes bestie, THIS is the optimized one. The glow-up. The main character. The _l3. 👑

🌟 Purpose: 🌟

👀 Attention

The algorithm is implemented using BLAS level 3 operations. That's the whole vibe. Fast, chunky, blocky, slay. ✨

🎀 Parameters

🌸 TRANSA

TRANSA is CHARACTER(1) and it says what form the system takes with respect to A and C:

== 'N': op1(A) = A, op1(C) = C (no transpose for A and C)
== 'T': op1(A) = A**T, op1(C) = C **T (transpose A and C)

🦋 TRANSB

TRANSB is CHARACTER(1) and it says what form the system takes with respect to B and D:

== 'N': op2(B) = B, op2(D) = D (no transpose for B and D)
== 'T': op2(B) = B**T, op2(D) = D **T (transpose B and D)

✨ SGN1

SGN1 is REAL, allowed values: +/-1. It specifies the sign between in the first equation. Lowkey just a lil plus or minus, bestie. 💅

🌈 SGN2

SGN2 is REAL, allowed values: +/-1. It specifies the sign between in the second equation. Same energy as SGN1, different equation. 🫶

💜 M

M is INTEGER. It's the order of the matrices A and C. M >= 0.

💖 N

N is INTEGER. It's the order of the matrices B and D. N >= 0.

🍓 A

A is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular. No exceptions, fr fr.

⭐ LDA

LDA is INTEGER. It's the leading dimension of the array A. LDA >= max(1,M).

🎀 B

B is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. And here's the plot twist: if the matrix D is already quasi-upper triangular, then the matrix B must be upper triangular. Iconic mutual-exclusion moment. 🙈

🧚 LDB

LDB is INTEGER. It's the leading dimension of the array B. LDB >= max(1,N).

🌸 C

C is REAL array, dimension (LDC,M). The matrix C must be upper triangular.

👑 LDC

LDC is INTEGER. It's the leading dimension of the array C. LDC >= max(1,M).

🦋 D

D is REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. And the mirror of before: if the matrix B is already quasi-upper triangular, then the matrix D must be upper triangular. It's giving "only one of you gets to be quasi at a time." 💀

💅 LDD

LDD is INTEGER. It's the leading dimension of the array D. LDD >= max(1,N).

🌈 E

E is REAL array, dimension (LDE,N).

On input, the matrix E contains the right hand side E.
On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2.

Right hand side E and the solution R are M-by-N matrices. Glow-up arc tbh. ✨

⭐ LDE

LDE is INTEGER. It's the leading dimension of the array E. LDE >= max(1,M).

🍓 F

F is REAL array, dimension (LDF,N).

On input, the matrix F contains the right hand side F.
On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2.

Right hand side F and the solution L are M-by-N matrices. Twin flame behaviour with E. 💖

🧚 LDF

LDF is INTEGER. It's the leading dimension of the array F. LDE >= max(1,M).

Okay bestie, real talk, that last one says LDE >= max(1,M) in the original — not LDF. Not me pointing it out, the page literally says that. Keeping it exactly as written because facts don't change even when the typo has main character energy. 😭

💜 SCALE

SCALE is REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0, otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true. It's the safety net keeping your numbers from absolutely losing it. 🫶

🌸 WORK

WORK is REAL array, dimension LWORK. Workspace for the algorithm.

The workspace needs to be queried before running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO. Yes bestie, you have to ask first. 💅

👑 INFO

INFO is INTEGER.

On input:

== -1: perform a workspace query
<> -1: normal operation

On exit, workspace query:

< 0: if INFO = -i, the i-th argument had an illegal value
>= 0: the value of INFO is the required number of elements in the workspace

On exit, normal operation:

== 0: successful exit 🎉
< 0: if INFO = -i, the i-th argument had an illegal value
> 0: the equation is not solved correctly. One of the arising inner system got singular. 💀

⚡ Optimizations:

Nothing but standard L3 BLAS calls. Clean, no tricks, just blocks.
Changed computation order to column major. ✨

✍️ Author

Martin Koehler, MPI Magdeburg

📅 Date

October 2025

📍 Definition

At line 223 of file sla_tgcsylv_l3_opt.f90.

🫧 subroutine sla_tgcsylv_l3_2s

subroutine sla_tgcsylv_l3_2s (character, dimension(1) transa, character,
    dimension(1) transb, real sgn1, real sgn2, integer m, integer n, real,
    dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer
    ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d,
    integer ldd, real, dimension(lde, *) e, integer lde, real,
    dimension(ldf,*) f, integer ldf, real scale, real, dimension(*) work,
    integer info)

Level-3 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation (optimized)

🌟 Purpose:

okay bestie so SLA_TGCSYLV_L3_2S solves a generalized coupled Sylvester equation that looks exactly like this 💅

op1(A) * R  + SGN1 * L  * op2(B) = SCALE * E                              (1)
op1(C) * R  + SGN2 * L  * op2(D) = SCALE * F

where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix, and B and D are N-by-N upper triangular matrices. And like… either B or D can be quasi upper triangular as well 🫶. The right hand side Y and the solution X are M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) ( or (D,B) ) are created by SGGES from LAPACK. that's just the lore, iykyk ✨

✨ Remarks ✨

This algorithm uses 2 Stage blocking with DAG scheduled solves. that's the whole flex, no cap 💅

✨ Parameters ✨

TRANSA — CHARACTER(1). Tells you the form of the system of equations with respect to A and C: if it's 'N' then op1(A) = A, op1(C) = C (no transpose for A and C), and if it's 'T' then op1(A) = A**T, op1(C) = C **T (transpose A and C). pick your fighter ✨
TRANSB — CHARACTER(1). Same energy but for B and D: 'N' gives you op2(B) = B, op2(D) = D (no transpose for B and D), 'T' gives you op2(B) = B**T, op2(D) = D **T (transpose B and D).
SGN1 — REAL, allowed values: +/-1. It's the sign between in the first equation. literally just a plus or minus, she's simple like that 🎀
SGN2 — REAL, allowed values: +/-1. The sign between in the second equation. same vibe, second equation 🦋
M — INTEGER. The order of the matrices A and C. M >= 0.
N — INTEGER. The order of the matrices B and D. N >= 0.
A — REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular. no exceptions bestie 😭
LDA — INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. But heads up: if the matrix D is already quasi-upper triangular, then the matrix B must be upper triangular 👀
LDB — INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — REAL array, dimension (LDC,M). The matrix C must be upper triangular. periodt.
LDC — INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. And same deal as B: if the matrix B is already quasi-upper triangular, then D must be upper triangular 🫶
LDD — INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — REAL array, dimension (LDE,N). On input, E holds the right hand side E. On output, E holds the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. The right hand side E and the solution R are M-by-N matrices. she does a little glow-up in place 💖
LDE — INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F — REAL array, dimension (LDF,N). On input, F holds the right hand side F. On output, F holds the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices. same transformation arc, slay ✨
LDF — INTEGER. The leading dimension of the array F. LDE >= max(1,M). (yes, it literally says LDE in the source, we do not fix, we do not question 🙈)
SCALE — REAL. It's a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0, otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true. so if you see a tiny scale, something was struggling 😭
WORK — REAL array, dimension LWORK. Workspace for the algorithm. And listen: the workspace needs to queried before the running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO. so you gotta ask before you act, main character energy 💅
INFO — INTEGER.
- On input:
  - == -1 : Perform a workspace query
  - <> -1: normal operation
- On exit, workspace query:
  - < 0 : if INFO = -i, the i-th argument had an illegal value
  - >= 0: The value of INFO is the required number of elements in the workspace.
- On exit, normal operation:
  - == 0: successful exit
  - < 0: if INFO = -i, the i-th argument had an illegal value
  - > 0: The equation is not solved correctly. One of the arising inner system got singular.

bestie if INFO > 0 comes back at you, the equation did not get solved correctly — one of the arising inner systems went singular. do not ignore her, she is telling you something 🔥

✍️ Author

Martin Koehler, MPI Magdeburg

🗓️ Date

October 2025

Definition at line 217 of file sla_tgcsylv_l3_2stage.f90. 💫

🦋 `sla_tgcsylv_l3_unopt`

okay and then there's this one:

subroutine sla_tgcsylv_l3_unopt (character, dimension(1) transa, character,
    dimension(1) transb, real sgn1, real sgn2, integer m, integer n, real,
    dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer
    ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d,
    integer ldd, real, dimension(lde, *) e, integer lde, real,
    dimension(ldf,*) f, integer ldf, real scale, real, dimension(*) work,
    integer info)

Level-3 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation (unoptimized) 🧚

Purpose:

okay bestie</p> so SLA_TGCSYLV_L3_UNOPT 💅 is the one that solves a generalized coupled Sylvester equation, and it looks like this:

op1(A) * R  + SGN1 * L  * op2(B) = SCALE * E                              (1)
op1(C) * R  + SGN2 * L  * op2(D) = SCALE * F

where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix, and B and D are N-by-N upper triangular matrices. And listen — either B or D can be quasi upper triangular too, it's flexible like that ✨. The right hand side Y and the solution X are M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) (or (D,B)) are the ones made by SGGES from LAPACK 🫶.

Attention ⚠️

The algorithm is implemented using BLAS level 3 operations. slay.

Parameters 🎀

TRANSA — TRANSA is CHARACTER(1). It says what form the system of equations takes with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (no transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (transpose A and C)
TRANSB — TRANSB is CHARACTER(1). It says what form the system of equations takes with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (no transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (transpose B and D)
SGN1 — SGN1 is REAL, allowed values: +/-1. Specifies the sign between in the first equation. (yes, that's literally what it says, bestie, it's the sign thing in equation one 😭)
SGN2 — SGN2 is REAL, allowed values: +/-1. Specifies the sign between in the second equation.
M — M is INTEGER. The order of the matrices A and C. M >= 0.
N — N is INTEGER. The order of the matrices B and D. N >= 0.
A — A is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular. no cap.
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — B is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular, the matrix B must be upper triangular. 👑
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — C is REAL array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular, the matrix D must be upper triangular. 💅
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — E is REAL array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices. ✨
LDE — LDE is INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F — F is REAL array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices.
LDF — LDF is INTEGER. The leading dimension of the array F. LDE >= max(1,M). (yes bestie it says LDE there in the original, we do not "fix" it, we simply report it 🙈)
SCALE — SCALE is REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0, otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true. 🫶
WORK — WORK is REAL array, dimension LWORK. Workspace for the algorithm. The workspace needs to be queried before running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO. iykyk 💜
INFO — INFO is INTEGER.
- On input:
  - == -1 : Perform a workspace query
  - <> -1: normal operation
- On exit, workspace query:
  - < 0 : if INFO = -i, the i-th argument had an illegal value
  - >= 0: The value of INFO is the required number of elements in the workspace.
- On exit, normal operation:
  - == 0: successful exit 🎉
  - < 0: if INFO = -i, the i-th argument had an illegal value
  - > 0: The equation is not solved correctly. One of the arising inner system got singular. 💀

Optimizations: 🔥

Nothing but standard L3 BLAS calls.

Author 🌸

Martin Koehler, MPI Magdeburg

Date 🦋

October 2025

Definition 📌

At line 223 of file sla_tgcsylv_l3_unopt.f90.

and now, bestie, meet the next icon — here's the signature, verbatim, no touchy:

recursive subroutine sla_tgcsylv_recursive (character, dimension(1) transa,
    character, dimension(1) transb, real sgn1, real sgn2, integer m, integer
    n, real, dimension(lda, *) a, integer lda, real, dimension(ldb, *) b,
    integer ldb, real, dimension(ldc, *) c, integer ldc, real,
    dimension(ldd, *) d, integer ldd, real, dimension(lde, *) e, integer
    lde, real, dimension(ldf,*) f, integer ldf, real scale, real,
    dimension(*) work, integer info)

Recursive Algorithm for the generalized coupled Sylvester equation. ✨

Purpose:

TRANSA — it's a CHARACTER(1), bestie, and it tells us what form the system is in for A and C:
- 'N': op1(A) = A, op1(C) = C — no transpose, we're keeping it natural 💅
- 'T': op1(A) = A**T, op1(C) = C **T — transposing A and C, flipping the whole vibe
TRANSB — also a CHARACTER(1), but this one handles B and D:
- 'N': op2(B) = B, op2(D) = D — no transpose for B and D, as-is, iconic
- 'T': op2(B) = B**T, op2(D) = D **T — transpose B and D, rotate the fit 🌀
SGN1 — a REAL, allowed values: +/-1. It's the sign between the two terms in the first equation, no biggie
SGN2 — a REAL too, allowed values: +/-1. Same deal but for the second equation 🫶
M — INTEGER, the order of the matrices A and C. So like, M >= 0, okay?
N — INTEGER, the order of the matrices B and D. Also N >= 0 💖
A — a REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular, that's the rule, no exceptions
LDA — INTEGER, the leading dimension of the array A. LDA >= max(1,M), got it? ✨
B — a REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. And listen — if the matrix D is already quasi-upper triangular, then B must be upper triangular. The two of them can't both be doing the quasi thing, that's the arrangement 💀
LDB — INTEGER, the leading dimension of the array B. LDB >= max(1,N)
C — a REAL array, dimension (LDC,M). The matrix C must be upper triangular. Periodt 💅
LDC — INTEGER, the leading dimension of the array C. LDC >= max(1,M)
D — a REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. And mirroring the B situation: if the matrix B is already quasi-upper triangular, then D must be upper triangular. Only one of them gets to be quasi, we've been over this 😭
LDD — INTEGER, the leading dimension of the array D. LDD >= max(1,N)
E — a REAL array, dimension (LDE,N).
- On input, the matrix E contains the right hand side E.
- On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2.
- Right hand side E and the solution R are M-by-N matrices. She works hard, she gets replaced by her own answer 💖
LDE — INTEGER, the leading dimension of the array E. LDE >= max(1,M)
F — a REAL array, dimension (LDF,N).
- On input, the matrix F contains the right hand side F.
- On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2.
- Right hand side F and the solution L are M-by-N matrices ✨
LDF — INTEGER, the leading dimension of the array F. LDE >= max(1,M) — yes, that's what it says, LDE, the original is being a little silly and we're preserving the lore fr fr 🙈
SCALE — a REAL. It's the scaling factor that keeps the result from overflowing, like a little seatbelt for your numbers. If INFO == 0 then SCALE is 1.0, otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true 🫶
WORK — a REAL array, dimension 1. Workspace for the algorithm. Tiny but mighty, one element of pure hustle 👑
INFO — INTEGER, the vibe report:
- == 0: successful exit, we ate 🔥
- < 0: if INFO = -i, the i-th argument had an illegal value
- > 0: The equation is not solved correctly. One of the arising inner system got singular. That's the 💀 outcome, no cap

Author

Martin Koehler, MPI Magdeburg 💜

Date

October 2025

Definition

at line 205 of file sla_tgcsylv_recursive.f90.

🌈 `sla_tgsylv_dag`

Okay bestie, new function, new era. Behold the signature, verbatim, don't touch it:

subroutine sla_tgsylv_dag (character, dimension(1) transa, character,
    dimension(1) transb, real sgn, integer m, integer n, real,
    dimension(lda, m) a, integer lda, real, dimension(ldb, n) b, integer
    ldb, real, dimension(ldc, m) c, integer ldc, real, dimension(ldd, n) d,
    integer ldd, real, dimension(ldx, n) x, integer ldx, real scale, real,
    dimension(*) work, integer info)

Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation with DAG based parallelization. It's giving organized chaos, and we're obsessed 🦋

Purpose

op1(A) * X * op2(B) + op1(C) * X * op2(D) = SCALE * Y                              (1)

op1(A) * X * op2(B) - op1(C) * X * op2(D) = SCALE * Y                              (2)

where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix, and B and D are N-by-N upper triangular matrices. plot twist: either B or D can be quasi upper triangular too 💅. The right hand side Y and the solution X are M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) ( or (D,B) ) are created by SGGES from LAPACK. iykyk 🫶

🎀 Remarks 🎀

The algorithm is implemented using BLAS level 3 operations and OpenMP 4.0 DAG scheduling. she's built different fr fr 🔥

⚠️ Attention ⚠️

Due to the parallel nature of the algorithm the scaling is not applied to the right hand. If the problem is ill-posed and scaling appears you have to solve the equation again with a solver with complete scaling support like SLA_TGSYLV_L3. no cap, don't skip this step bestie 😭

🌸 Parameters 🌸

TRANSA — is CHARACTER(1). Specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)
TRANSB — is CHARACTER(1). Specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)
SGN — is REAL, allowed values: +/-1. Specifies the sign between the two parts of the Sylvester equation.
- = 1 : Solve Equation (1)
- == -1: Solve Equation (2)
M — is INTEGER. The order of the matrices A and C. M >= 0.
N — is INTEGER. The order of the matrices B and D. N >= 0.
A — is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular.
LDA — is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular the matrix B must be upper triangular.
LDB — is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — is REAL array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — is REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular.
LDD — is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
X — is REAL array, dimension (LDX,N). On input, the matrix X contains the right hand side Y. On output, the matrix X contains the solution of Equation (1) or (2) as selected by TRANSA, TRANSB, and SGN. Right hand side Y and the solution X are M-by-N matrices.
LDX — is INTEGER. The leading dimension of the array X. LDX >= max(1,M).
SCALE — is REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true. If SCALE .NE. 1 the problem is no solved correctly in this case one have to use an other solver. red flag, bestie 🚩
WORK — is REAL array, dimension LWORK. Workspace for the algorithm. The workspace needs to queried before the running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO. (yes, you gotta ask first, she's not gonna just tell you 💅)
INFO — is INTEGER
- On input:
  - == -1 : Perform a workspace query
  - <> -1: normal operation
- On exit, workspace query:
  - < 0 : if INFO = -i, the i-th argument had an illegal value
  - >= 0: The value of INFO is the required number of elements in the workspace.
- On exit, normal operation:
  - == 0: successful exit 🎉
  - < 0: if INFO = -i, the i-th argument had an illegal value
  - > 0: The equation is not solved correctly. One of the arising inner system got singular.

👑 Author 👑

Martin Koehler, MPI Magdeburg

📅 Date 📅

October 2025

Definition at line 208 of file sla_tgsylv_dag.f90.

✨ subroutine `sla_tgsylv_gardiner_laub` ✨

subroutine sla_tgsylv_gardiner_laub (character, dimension(1) transa,
    character, dimension(1) transb, real sgn, integer m, integer n, real,
    dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer
    ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d,
    integer ldd, real, dimension(ldx, *) x, integer ldx, real scale, real,
    dimension(2*m,*) work, integer info)

Level-2 Gardiner-Laub algorithm for the generalized Sylvester equation. she's the classic, the blueprint, the moment 🦋

Purpose:

Okay so SLA_TGSYLV_GARDINER_LAUB solves a generalized Sylvester equation of the following forms, bestie, pay attention because there are literally two of them:

   op1(A) * X * op2(B) + op1(C) * X * op2(D) = SCALE * Y                              (1)

   op1(A) * X * op2(B) - op1(C) * X * op2(D) = SCALE * Y                              (2)

where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix and B and D are N-by-N upper triangular matrices. Furthermore either B or D can be quasi upper triangular as well. The right hand side Y and the solution X are M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) ( or (D,B) ) are created by SGGES from LAPACK. ✨

🧠 Remarks

The algorithm implements the ideas by Gardiner et. al. on top of BLAS level 2 operations. It's giving built different, no cap. 💅

🚨 Attention

Scaling is not supported by this algorithm so SCALE == 1. That's it, that's the whole vibe — don't come for me about scaling, it's simply not on the menu. 😭

🎀 Parameters

TRANSA — TRANSA is CHARACTER(1).

Specifies the form of the system of equations with respect to A and C:

== 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
== 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)

TRANSB — TRANSB is CHARACTER(1).

Specifies the form of the system of equations with respect to B and D:

== 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
== 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)

SGN — SGN is REAL, allowed values: +/-1.

Specifies the sign between the two parts of the Sylvester equation.

= 1 : Solve Equation (1)
== -1: Solve Equation (2)

M — M is INTEGER.

The order of the matrices A and C. M >= 0.

N — N is INTEGER.

The order of the matrices B and D. N >= 0.

A — A is REAL array, dimension (LDA,M).

The matrix A must be (quasi-) upper triangular.

LDA — LDA is INTEGER.

The leading dimension of the array A. LDA >= max(1,M).

B — B is REAL array, dimension (LDB,N).

The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular the matrix B must be upper triangular.

LDB — LDB is INTEGER.

The leading dimension of the array B. LDB >= max(1,N).

C — C is REAL array, dimension (LDC,M).

The matrix C must be upper triangular.

LDC — LDC is INTEGER.

The leading dimension of the array C. LDC >= max(1,M).

D — D is REAL array, dimension (LDD,N).

The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular.

LDD — LDD is INTEGER.

The leading dimension of the array D. LDD >= max(1,N).

X — X is REAL array, dimension (LDX,N).

On input, the matrix X contains the right hand side Y. On output, the matrix X contains the solution of Equation (1) or (2) as selected by TRANSA, TRANSB, and SGN. Right hand side Y and the solution X are M-by-N matrices. 💖

LDX — LDX is INTEGER.

The leading dimension of the array X. LDB >= max(1,M).

bestie I'm just the messenger 💀 that's literally what it says, LDB, we do not touch the words, we only read them

SCALE — SCALE is REAL.

SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.

WORK — WORK is REAL array, dimension LWORK.

Workspace for the algorithm. The workspace needs to queried before the running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO. ✨ (yes, you query first, you don't just rawdog the array, iykyk)

INFO — INFO is INTEGER.

On input:

== -1 : Perform a workspace query
<> -1: normal operation

On exit, workspace query:

< 0 : if INFO = -i, the i-th argument had an illegal value
>= 0: The value of INFO is the required number of elements in the workspace.

On exit, normal operation:

== 0: successful exit
< 0: if INFO = -i, the i-th argument had an illegal value
> 0: The equation is not solved correctly. One of the arising inner system got singular.

👑 Author

Martin Koehler, MPI Magdeburg

📅 Date

October 2025

📍 Definition

Definition at line 204 of file sla_tgsylv_gardiner_laub.f90.

subroutine sla_tgsylv_l2 (character, dimension(1) transa, character,
    dimension(1) transb, real sgn, integer m, integer n, real,
    dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer
    ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d,
    integer ldd, real, dimension(ldx, *) x, integer ldx, real scale, real,
    dimension(*) work, integer info)

Level-2 Bartels-Stewart Algorithm for the generalized Sylvester equation (Optimized) 🦋

Purpose:

Okay so bestie, SLA_TGSYLV_L2 is literally out here solving a generalized Sylvester equation, and it comes in two iconic forms:

op1(A) * X * op2(B) + op1(C) * X * op2(D) = SCALE * Y                              (1)

op1(A) * X * op2(B) - op1(C) * X * op2(D) = SCALE * Y                              (2)

where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix and B and D are N-by-N upper triangular matrices. Furthermore either B or D can be quasi upper triangular as well. The right hand side Y and the solution X M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) ( or (D,B) ) are created by SGGES from LAPACK. iykyk 💅

Remarks ✨

The algorithm is implemented using BLAS level 2 operations and hand written optimizations. obsessed with the hand-written glow-up tbh

Parameters 💅 — vibe check every single one:

TRANSA — TRANSA is CHARACTER(1). Specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)
TRANSB — TRANSB is CHARACTER(1). Specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)
SGN — SGN is REAL, allowed values: +/-1. Specifies the sign between the two parts of the Sylvester equation.
- = 1 : Solve Equation (1)
- == -1: Solve Equation (2)
M — M is INTEGER. The order of the matrices A and C. M >= 0.
N — N is INTEGER. The order of the matrices B and D. N >= 0.
A — A is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular.
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — B is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular the matrix B must be upper triangular.
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — C is REAL array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular.
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
X — X is REAL array, dimension (LDX,N). On input, the matrix X contains the right hand side Y. On output, the matrix X contains the solution of Equation (1) or (2) as selected by TRANSA, TRANSB, and SGN. Right hand side Y and the solution X are M-by-N matrices.
LDX — LDX is INTEGER. The leading dimension of the array X. LDB >= max(1,M). not me side-eyeing that LDB in the LDX section 🙈 keeping it exactly as written tho
SCALE — SCALE is REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — WORK is REAL array, dimension LWORK. Workspace for the algorithm. The workspace needs to queried before the running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO.
INFO — INFO is INTEGER. On input:
- == -1 : Perform a workspace query
- <> -1: normal operation

On exit, workspace query:

< 0 : if INFO = -i, the i-th argument had an illegal value
>= 0: The value of INFO is the required number of elements in the workspace.

On exit, normal operation:

== 0: successful exit
< 0: if INFO = -i, the i-th argument had an illegal value
> 0: The equation is not solved correctly. One of the arising inner system got singular.

Optimizations 🚀

Replaced level-3 by level-2 and level-1 calls,

Author: Martin Koehler, MPI Magdeburg

Date: October 2025

Definition at line 207 of file sla_tgsylv_l2.f90.

Next up, the colwise sister function:

   subroutine sla_tgsylv_l2_colwise (character, dimension(1) transa, character,
       dimension(1) transb, real sgn, integer m, integer n, real,
       dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer
       ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d,
       integer ldd, real, dimension(ldx, *) x, integer ldx, real scale, real,
       dimension(*) work, integer info)

Level-2 Bartels-Stewart Algorithm for the generalized Sylvester equation (unoptimized)

Purpose:

bestie okay so real talk — SLA_TGSYLV_L2_UNOPT solves a generalized Sylvester equation, and it comes in exactly two flavors:

op1(A) * X * op2(B) + op1(C) * X * op2(D) = SCALE * Y                              (1)

op1(A) * X * op2(B) - op1(C) * X * op2(D) = SCALE * Y                              (2)

where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix, and B and D are N-by-N upper triangular matrices. And plot twist: either B or D can be quasi upper triangular as well. The right hand side Y and the solution X are M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) (or (D,B)) are created by SGGES from LAPACK. 💅

Attention 🙈

The algorithm is implemented using BLAS level 2 operations without further optimizations. For a faster implementation see SLA_TGSYLV_L2. (yes queen, there's a glow-up version, go get her)

Parameters 🎀

TRANSA — is CHARACTER(1). Specifies the form of the system of equations with respect to A and C:
- 'N': op1(A) = A, op1(C) = C (no transpose for A and C)
- 'T': op1(A) = A**T, op1(C) = C **T (transpose A and C)
TRANSB — is CHARACTER(1). Specifies the form of the system of equations with respect to B and D:
- 'N': op2(B) = B, op2(D) = D (no transpose for B and D)
- 'T': op2(B) = B**T, op2(D) = D **T (transpose B and D)
SGN — is REAL, allowed values: +/-1. Specifies the sign between the two parts of the Sylvester equation.
- = 1 : solve Equation (1)
- == -1: solve Equation (2)
M — is INTEGER. The order of the matrices A and C. M >= 0.
N — is INTEGER. The order of the matrices B and D. N >= 0.
A — is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular.
LDA — is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular, the matrix B must be upper triangular.
LDB — is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — is REAL array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — is REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular, the matrix D must be upper triangular.
LDD — is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
X — is REAL array, dimension (LDX,N).
- On input, the matrix X contains the right hand side Y.
- On output... sorry, on output, the matrix X contains the solution of Equation (1) or (2) as selected by TRANSA, TRANSB, and SGN.
- Right hand side Y and the solution X are M-by-N matrices.
LDX — is INTEGER. The leading dimension of the array X. LDB >= max(1,M). (yes it says LDB in the source, we're just reporting the vibes fr fr 💀)
SCALE — is REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0, otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — is REAL array, dimension LWORK. Workspace for the algorithm. The workspace needs to be queried before running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO.
INFO — is INTEGER.
- On input:
  - == -1 : perform a workspace query
  - <> -1: normal operation
- On exit, workspace query:
  - < 0 : if INFO = -i, the i-th argument had an illegal value
  - >= 0: the value of INFO is the required number of elements in the workspace
- On exit, normal operation:
  - == 0: successful exit 🎉
  - < 0: if INFO = -i, the i-th argument had an illegal value
  - > 0: the equation is not solved correctly. One of the arising inner system got singular. not the slay we wanted, bestie

Optimizations: 🌈

Nothing. (the audacity. the honesty. iconic.)

Author 👑

Martin Koehler, MPI Magdeburg

Date 📅

October 2025

Definition 📌

At line 208 of file sla_tgsylv_l2_colwise.f90.

🦋 `sla_tgsylv_l2_local_copy` 🦋

okay and then we got this other one, same energy, different name:

subroutine sla_tgsylv_l2_local_copy (character, dimension(1) transa,
    character, dimension(1) transb, real sgn, integer m, integer n, real,
    dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer
    ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d,
    integer ldd, real, dimension(ldx, *) x, integer ldx, real scale, real,
    dimension(*) work, integer info)

Level-2 Bartels-Stewart Algorithm for the generalized Sylvester equation (Optimized)

Purpose:

⚠️ Attention ⚠️

The algorithm solves problems of dimension at most 128 ( M,N <=128 ). Like, that's the ceiling bestie, do not test it. 🚫

🌸 Parameters 🌸

TRANSA — TRANSA is CHARACTER(1). It picks the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (no transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C**T (transpose A and C)
TRANSB — TRANSB is CHARACTER(1). Same energy but for B and D:
- == 'N': op2(B) = B, op2(D) = D (no transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D**T (transpose B and D)
SGN — SGN is REAL, allowed values: +/-1. It decides the sign between the two parts of the Sylvester equation, the literal plus-or-minus moment 💅:
- = 1 : Solve Equation (1)
- == -1: Solve Equation (2)
M — M is INTEGER. The order of the matrices A and C. So: 128 >= M >= 0.
N — N is INTEGER. The order of the matrices B and D. And: 128 >= N >= 0.
A — A is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular. Non-negotiable, she's structured like that. 👑
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — B is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. And listen: if the matrix D is already quasi-upper triangular, then the matrix B must be upper triangular. It's giving tradeoff, okay? 🔄
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — C is REAL array, dimension (LDC,M). The matrix C must be upper triangular. No quasi anything, just straight upper triangular, periodt.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. And same deal as before but mirrored: if the matrix B is already quasi-upper triangular, then the matrix D must be upper triangular. 🤝
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
X — X is REAL array, dimension (LDX,N). On input, the matrix X contains the right hand side Y. On output, the matrix X contains the solution of Equation (1) or (2) as selected by TRANSA, TRANSB, and SGN. Right hand side Y and the solution X are M-by-N matrices. So she walks in as the problem and struts out as the answer, we love a glow-up ✨
LDX — LDX is INTEGER. The leading dimension of the array X. LDB >= max(1,M).

Heads up bestie: yes, that LDB is literally what the original says, right there in the LDX entry. I'm not fixing it, I'm not "correcting" it, I'm just reporting it exactly as written. 🙈

SCALE — SCALE is REAL. It's a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0; otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true. She's the safety net and we stan. 🛟
WORK — WORK is REAL array, dimension LWORK. Workspace for the algorithm. And the workspace needs to be queried before running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO. Ask first, then run. Consent culture for memory allocation fr fr. 🙏
INFO — INFO is INTEGER.

On input:

== -1 : Perform a workspace query
<> -1: normal operation

On exit, workspace query:

< 0 : if INFO = -i, the i-th argument had an illegal value
>= 0: The value of INFO is the required number of elements in the workspace.

On exit, normal operation:

== 0: successful exit 🎉
< 0: if INFO = -i, the i-th argument had an illegal value
> 0: The equation is not solved correctly. One of the arising inner system got singular. 💀

🔥 Optimizations: 🔥

Replaced level-3 by level-2 and level-1 calls,
Reorder the solution order to column first style,
Replaced SAXPY operation by Fortran intrinsic
Replaced all BLAS calls except of STRMV and SGER by Fortran intrinsic,
Use local copies of A, B, C, D, and X ( M, N <=128 ) .
STRMV works sequential.

👑 Author 👑

Martin Koehler, MPI Magdeburg

📅 Date 📅

October 2025

📍 Definition 📍

At line 214 of file sla_tgsylv_l2_opt_local_copy.f90.

subroutine sla_tgsylv_l2_local_copy_128 ( transa,  transb,  sgn,  m,  n,  a,
    lda,  b,  ldb,  c,  ldc,  d,  ldd,  x,  ldx,  scale,  work,  info)

Level-2 Bartels-Stewart Algorithm for the generalized Sylvester equation

(Optimized)

Purpose:

💖 `SLA_TGSYLV_L2_LOCAL_COPY_128` 💖

okay so bestie, this one solves a generalized Sylvester equation, and it comes in two flavors:

op1(A) * X * op2(B) + op1(C) * X * op2(D) = SCALE * Y                              (1)

op1(A) * X * op2(B) - op1(C) * X * op2(D) = SCALE * Y                              (2)

where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix, and B and D are N-by-N upper triangular matrices. Also?? either B or D can be quasi upper triangular too. The right hand side Y and the solution X are M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) (or (D,B)) are created by SGGES from LAPACK. iykyk ✨

🌸 Attention 🌸

The size of the Problem is limited by M,N <= 128 — that's the whole vibe check, don't fight it 💀

✨ Parameters ✨

TRANSA — TRANSA is CHARACTER(1). It picks the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)
TRANSB — TRANSB is CHARACTER(1). It picks the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)
SGN — SGN is REAL, allowed values: +/-1. It's the sign between the two parts of the Sylvester equation, literally the hinge of the whole thing.
- = 1 : Solve Equation (1)
- == -1: Solve Equation (2)
M — M is INTEGER. The order of the matrices A and C. 128 >= M >= 0. 👑
N — N is INTEGER. The order of the matrices B and D. 128 >= N >= 0. 👑
A — A is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular. no cap, it has to be.
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — B is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular the matrix B must be upper triangular. it's giving ✨mutual exclusivity✨
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — C is REAL array, dimension (LDC,M). The matrix C must be upper triangular. periodt.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular. see? told you it was mutual 🦋
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
X — X is REAL array, dimension (LDX,N). On input, the matrix X contains the right hand side Y. On output, the matrix X contains the solution of Equation (1) or (2) as selected by TRANSA, TRANSB, and SGN. Right hand side Y and the solution X are M-by-N matrices. she's doing double duty and she's slaying 💅
LDX — LDX is INTEGER. The leading dimension of the array X. LDB >= max(1,M).

bestie it literally says LDB right there in the source, I am just the messenger, I do not make the rules 💀 (you and I both know it means LDX)

SCALE — SCALE is REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true. she's the safety net 🫶
WORK — WORK is REAL array, dimension LWORK. Workspace for the algorithm. The workspace needs to queried before the running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO. so like, ask first, then compute. communication is key 🔥
INFO — INFO is INTEGER

On input:

== -1 : Perform a workspace query
<> -1: normal operation

On exit, workspace query:

< 0 : if INFO = -i, the i-th argument had an illegal value
>= 0: The value of INFO is the required number of elements in the workspace.

On exit, normal operation:

== 0: successful exit 🎉
< 0: if INFO = -i, the i-th argument had an illegal value
> 0: The equation is not solved correctly. One of the arising inner system got singular. 😭

🎀 Optimizations: 🎀

Replaced level-3 by level-2 and level-1 calls,
Reorder the solution order to column first style,
Replaced SAXPY operation by Fortran intrinsic
Replaced all BLAS calls except of STRMV by Fortran intrinsic,
Use local copies of A, B, C, D, and X.
STRMV works purely sequential
Align local copies and fix problem size to <= 128

👑 Author 👑

Martin Koehler, MPI Magdeburg

🦋 Date 🦋

October 2025

⭐ Definition ⭐

Definition at line 215 of file sla_tgsylv_l2_opt_local_copy_128.f90.

🌈 `sla_tgsylv_l2_local_copy_32` 🌈

subroutine sla_tgsylv_l2_local_copy_32 ( transa,  transb,  sgn,  m,  n,  a,
    lda,  b,  ldb,  c,  ldc,  d,  ldd,  x,  ldx,  scale,  work,  info)

Level-2 Bartels-Stewart Algorithm for the generalized Sylvester equation (Optimized) 🧚

💜 Purpose:

INFO — it's an INTEGER, bestie, and it's doing double duty depending on when you look at it 👀

On input:

== -1 → do a workspace query ✨
<> -1 → normal operation, business as usual 💅

On exit (workspace query mode):

< 0 → if INFO = -i, the i-th argument had an illegal value 😭
>= 0 → the value of INFO is the required number of elements in the workspace 🎀

On exit (normal operation mode):

== 0 → successful exit, we ate 🍓
< 0 → if INFO = -i, the i-th argument had an illegal value 🙈
> 0 → the equation is not solved correctly. One of the arising inner systems got singular 💀

🌸 Optimizations: 🌸

The glow-up era, fr fr:

Replaced level-3 by level-2 and level-1 calls
Reorder the solution order to column first style
Replaced SAXPY operation by Fortran intrinsic
Replaced all BLAS calls except of STRMV by Fortran intrinsic
Use local copies of A, B, C, D, and X.
STRMV works purely sequential
Align local copies and fix problem size to <= 32

💜 Author 💜

Martin Koehler, MPI Magdeburg

📅 Date 📅

October 2025

Definition at line 215 of file sla_tgsylv_l2_opt_local_copy_32.f90.

subroutine sla_tgsylv_l2_local_copy_64 ( transa,  transb,  sgn,  m,  n,  a,
    lda,  b,  ldb,  c,  ldc,  d,  ldd,  x,  ldx,  scale,  work,  info)

Level-2 Bartels-Stewart Algorithm for the generalized Sylvester equation

(Optimized) 👑

✨ Purpose: ✨

🚨 Attention 🚨

bestie, read this before you do <em>anything</em> else: the size of the Problem is limited by M,N <= 64. that's the whole ceiling. sixty-four. not sixty-five. do not try it. 🫶

📋 Parameters 📋

TRANSA — okay so TRANSA is CHARACTER(1), and it decides the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)
TRANSB — TRANSB is CHARACTER(1) too, same deal but for B and D:
- == 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)
SGN — SGN is REAL, and the allowed values are +/-1. it's literally the sign between the two parts of the Sylvester equation, the vibe check of the whole thing:
- = 1 : Solve Equation (1)
- == -1: Solve Equation (2)
M — M is INTEGER. it's the order of the matrices A and C. so like, 64 >= M >= 0. ✨
N — N is INTEGER. the order of the matrices B and D. 64 >= N >= 0. yes girl, another ceiling. 💅
A — A is REAL array, dimension (LDA,M). and listen, the matrix A must be (quasi-) upper triangular. non-negotiable.
LDA — LDA is INTEGER. the leading dimension of the array A. LDA >= max(1,M). that's the rule, no notes.
B — B is REAL array, dimension (LDB,N). the matrix B must be (quasi-) upper triangular. but wait — if the matrix D is already quasi-upper triangular, then the matrix B must be upper triangular. so it's giving conditional energy, iykyk.
LDB — LDB is INTEGER. leading dimension of the array B. LDB >= max(1,N). 👑
C — C is REAL array, dimension (LDC,M). the matrix C must be upper triangular. that's it, that's the assignment.
LDC — LDC is INTEGER. leading dimension of the array C. LDC >= max(1,M).
D — D is REAL array, dimension (LDD,N). the matrix D must be (quasi-) upper triangular. and same little plot twist as B: if the matrix B is already quasi-upper triangular, the matrix D must be upper triangular. 🌸
LDD — LDD is INTEGER. leading dimension of the array D. LDD >= max(1,N).
X — X is REAL array, dimension (LDX,N). on input, the matrix X contains the right hand side Y. on output, the matrix X contains the solution of Equation (1) or (2) as selected by TRANSA, TRANSB, and SGN. and for the record: right hand side Y and the solution X are M-by-N matrices. slay. ✨
LDX — LDX is INTEGER. the leading dimension of the array X. and the page says, character for character, LDB >= max(1,M). yes bestie, it literally says LDB there. i am not the one who wrote it, i'm just the one who's telling you. 🙈
SCALE — SCALE is REAL. it's a scaling factor to prevent the overflow in the result. if INFO == 0 then SCALE is 1.0; otherwise, if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true. so basically SCALE is out here doing damage control. 😭
WORK — WORK is REAL array, dimension LWORK. it's the workspace for the algorithm. and the workspace needs to queried before the running the computation. the query is performed by calling the subroutine with INFO == -1 on input, and the required workspace is then returned in INFO. so query first, queen. 👑
INFO — INFO is INTEGER.
- On input:
  - == -1 : Perform a workspace query
  - <> -1: normal operation
- On exit, workspace query:
  - < 0 : if INFO = -i, the i-th argument had an illegal value
  - >= 0: The value of INFO is the required number of elements in the workspace.
- On exit, normal operation:
  - == 0: successful exit
  - < 0: if INFO = -i, the i-th argument had an illegal value
  - > 0: The equation is not solved correctly. One of the arising inner system got singular.

not the singular inner system 😭 if INFO > 0 comes back at you, that's what happened. no cap.

⚡ Optimizations: ⚡

okay so this version is optimized and here's the receipts, bestie:

Replaced level-3 by level-2 and level-1 calls,
Reorder the solution order to column first style,
Replaced SAXPY operation by Fortran intrinsic
Replaced all BLAS calls except of STRMV by Fortran intrinsic,
Use local copies of A, B, C, D, and X.
STRMV works purely sequential
Align local copies and fix problem size to <= 64

literally every one of those bullets is a glow-up. the local copies? iconic. the column-first reorder? main character energy. 🦋

🎀 Author 🎀

Martin Koehler, MPI Magdeburg

📅 Date 📅

October 2025

🔖 Definition 🔖

Definition at line 215 of file sla_tgsylv_l2_opt_local_copy_64.f90.

🌈 subroutine `sla_tgsylv_l2_local_copy_96` 🌈

subroutine sla_tgsylv_l2_local_copy_96 ( transa,  transb,  sgn,  m,  n,  a,
    lda,  b,  ldb,  c,  ldc,  d,  ldd,  x,  ldx,  scale,  work,  info)

Level-2 Bartels-Stewart Algorithm for the generalized Sylvester equation (Optimized)

💫 Purpose: 💫

where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix, and B and D are N-by-N upper triangular matrices. Also — either B or D can be quasi upper triangular as well, bestie, it's flexible like that. The right hand side Y and the solution X are M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) ( or (D,B) ) are created by SGGES from LAPACK. ✨

🚨 Attention 🚨

The size of the Problem is limited by M,N <= 96

No cap, that's the hard ceiling. Ninety-six. Do not come for her with a 97. 🙈

💅 Parameters 💅

TRANSA — TRANSA is CHARACTER(1). It specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)
TRANSB — TRANSB is CHARACTER(1). It specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)
SGN — SGN is REAL, allowed values: +/-1. It specifies the sign between the two parts of the Sylvester equation.
- = 1 : Solve Equation (1)
- == -1: Solve Equation (2)
M — M is INTEGER. The order of the matrices A and C. 96 >= M >= 0.
N — N is INTEGER. The order of the matrices B and D. 96 >= N >= 0.
A — A is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular. 👑
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — B is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular the matrix B must be upper triangular. (yes, it's a whole agreement between the two of them, iykyk 🤝)
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — C is REAL array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular. (same deal as above, just swapped 🫶)
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
X — X is REAL array, dimension (LDX,N). On input, the matrix X contains the right hand side Y. On output, the matrix X contains the solution of Equation (1) or (2) as selected by TRANSA, TRANSB, and SGN. Right hand side Y and the solution X are M-by-N matrices. she goes in as the question and comes out as the answer 💖
LDX — LDX is INTEGER. The leading dimension of the array X. LDB >= max(1,M). 👀 (fr fr that's literally what the source says, we are simply reporting the lore)
SCALE — SCALE is REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — WORK is REAL array, dimension LWORK. Workspace for the algorithm. The workspace needs to queried before the running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO. 📋
INFO — INFO is INTEGER.

On input:

== -1 : Perform a workspace query
<> -1: normal operation

On exit, workspace query:

< 0 : if INFO = -i, the i-th argument had an illegal value
>= 0: The value of INFO is the required number of elements in the workspace.

On exit, normal operation:

== 0: successful exit 🎉
< 0: if INFO = -i, the i-th argument had an illegal value
> 0: The equation is not solved correctly. One of the arising inner system got singular. 💀

⚡ Optimizations: ⚡

Replaced level-3 by level-2 and level-1 calls,
Reorder the solution order to column first style,
Replaced SAXPY operation by Fortran intrinsic
Replaced all BLAS calls except of STRMV by Fortran intrinsic,
Use local copies of A, B, C, D, and X.
STRMV works purely sequential
Align local copies and fix problem size to <= 96

👑 Author 👑

Martin Koehler, MPI Magdeburg

📅 Date 📅

October 2025

📍 Definition 📍

Definition at line 215 of file sla_tgsylv_l2_opt_local_copy_96.f90.

subroutine sla_tgsylv_l2_reorder (character, dimension(1) transa, character,
    dimension(1) transb, real sgn, integer m, integer n, real,
    dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer
    ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d,
    integer ldd, real, dimension(ldx, *) x, integer ldx, real scale, real,
    dimension(*) work, integer info)

Level-2 Bartels-Stewart Algorithm for the generalized Sylvester equation (Optimized) 🦋

🌸 Purpose: 🌸

okay so SLA_TGSYLV_L2_REORDER solves a generalized Sylvester equation, and bestie it comes in two iconic flavors:

op1(A) * X * op2(B) + op1(C) * X * op2(D) = SCALE * Y   (1)

op1(A) * X * op2(B) - op1(C) * X * op2(D) = SCALE * Y   (2)

where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix, and B and D are N-by-N upper triangular matrices. ✨ And like, either B or D can be quasi upper triangular as well — we're not greedy, one quasi is plenty. The right hand side Y and the solution X are M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) ( or (D,B) ) are created by SGGES from LAPACK. 🫶

💅 Attention 💅

The algorithm is implemented using BLAS level 2 operations and hand written optimizations. no cap, it's literally hand-optimized.

🎀 Parameters 🎀

TRANSA — it's a CHARACTER(1) and it decides how we treat A and C: 'N' means op1(A) = A, op1(C) = C (no transpose for A and C, keeping it natural 💅) and 'T' means op1(A) = A**T, op1(C) = C **T (transpose A and C, periodt).
TRANSB — also a CHARACTER(1), same energy but for B and D: 'N' gives you op2(B) = B, op2(D) = D (no transpose for B and D) and 'T' gives you op2(B) = B**T, op2(D) = D **T (transpose B and D).
SGN — this one's REAL, and the allowed values are +/-1 🎯. It sets the sign between the two parts of the Sylvester equation: = 1 → Solve Equation (1), and == -1 → Solve Equation (2).
M — INTEGER. The order of the matrices A and C. M >= 0.
N — INTEGER. The order of the matrices B and D. N >= 0.
A — REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular. no exceptions bestie.
LDA — INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. And here's the lil catch: if the matrix D is already quasi-upper triangular, then the matrix B must be upper triangular. Only one of them gets to be quasi, we share. 🤝
LDB — INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — REAL array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. And the mirror rule: if the matrix B is already quasi-upper triangular, then the matrix D must be upper triangular.
LDD — INTEGER. The leading dimension of the array D. LDD >= max(1,N).
X — REAL array, dimension (LDX,N). On input, the matrix X contains the right hand side Y. On output, the matrix X contains the solution of Equation (1) or (2) as selected by TRANSA, TRANSB, and SGN. Right hand side Y and the solution X are M-by-N matrices. it does a whole glow-up in place ✨
LDX — INTEGER. The leading dimension of the array X. The doc literally says LDB >= max(1,M) right here — yep, LDB, not LDX 🙈 — and we're passing that along exactly as written, character for character, because that is our whole thing.
SCALE — REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0, otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true. it's the seatbelt 🚗
WORK — REAL array, dimension LWORK. Workspace for the algorithm. The workspace needs to queried before the running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO. so like, ask first, allocate second 🎀
INFO — INTEGER.

On input:

== -1 : Perform a workspace query
<> -1: normal operation

On exit, workspace query:

< 0 : if INFO = -i, the i-th argument had an illegal value
>= 0: The value of INFO is the required number of elements in the workspace.

On exit, normal operation:

== 0: successful exit 🎉
< 0: if INFO = -i, the i-th argument had an illegal value
> 0: The equation is not solved correctly. One of the arising inner system got singular. 💀

✨ Optimizations ✨

Replaced level-3 by level-2 and level-1 calls,
Reorder the solution order to column first style,
Replaced SAXPY operation by Fortran intrinsic
Replaced all BLAS calls except of STRMV and SGER by Fortran intrinsic,

👑 Author 👑

Martin Koehler, MPI Magdeburg

📅 Date 📅

October 2025

📍 Definition 📍

Definition at line 211 of file sla_tgsylv_l2_reorder.f90.

subroutine sla_tgsylv_l2_rowwise (character, dimension(1) transa, character,
    dimension(1) transb, real sgn, integer m, integer n, real,
    dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer
    ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d,
    integer ldd, real, dimension(ldx, *) x, integer ldx, real scale, real,
    dimension(*) work, integer info)

Level-2 Bartels-Stewart Algorithm for the generalized Sylvester equation (unoptimized) 🦋

🎯 Purpose 🎯

op1(A) * X * op2(B) + op1(C) * X * op2(D) = SCALE * Y                              (1)

or

op1(A) * X * op2(B) - op1(C) * X * op2(D) = SCALE * Y                              (2)

where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix and B and D are N-by-N upper triangular matrices. And ok bestie, plot twist: either B or D can be quasi upper triangular as well 🦋. The right hand side Y and the solution X are M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) ( or (D,B) ) are created by SGGES from LAPACK. iykyk 💅

✨ Attention ✨

The algorithm is implemented using BLAS level 2 operations without further optimizations. For a faster implementation see SLA_TGSYLV_L2 . 🏃‍♀️💨

✨ Parameters ✨

TRANSA — TRANSA is CHARACTER(1). It says which form the system of equations takes with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)
TRANSB — TRANSB is CHARACTER(1). Same energy, but for B and D:
- == 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)
SGN — SGN is REAL, allowed values: +/-1. It picks the sign between the two parts of the Sylvester equation, no cap:
- = 1 : Solve Equation (1)
- == -1: Solve Equation (2)
M — M is INTEGER. The order of the matrices A and C. M >= 0. 🫶
N — N is INTEGER. The order of the matrices B and D. N >= 0. 🫶
A — A is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular.
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — B is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. And listen — if the matrix D is already quasi-upper triangular, then the matrix B must be upper triangular. Only one of them gets to be the ✨spicy✨ one at a time, bestie.
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — C is REAL array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. Same rule as above but flipped: if the matrix B is already quasi-upper triangular, then the matrix D must be upper triangular 💀
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
X — X is REAL array, dimension (LDX,N). On input, the matrix X contains the right hand side Y. On output, the matrix X contains the solution of Equation (1) or (2) as selected by TRANSA, TRANSB, and SGN. Right hand side Y and the solution X are M-by-N matrices. It literally serves both looks, main character energy 👑
LDX — LDX is INTEGER. The leading dimension of the array X. LDB >= max(1,M).
SCALE — SCALE is REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true. 🫧
WORK — WORK is REAL array, dimension LWORK. Workspace for the algorithm. The workspace needs to queried before the running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO. Yes bestie, you gotta ask first, it's giving ✨consent✨ 🎀
INFO — INFO is INTEGER.
- On input:
  - == -1 : Perform a workspace query
  - <> -1: normal operation
- On exit, workspace query:
  - < 0 : if INFO = -i, the i-th argument had an illegal value
  - >= 0: The value of INFO is the required number of elements in the workspace.
- On exit, normal operation:
  - == 0: successful exit 🎉
  - < 0: if INFO = -i, the i-th argument had an illegal value
  - > 0: The equation is not solved correctly. One of the arising inner system got singular. 😭

✨ Optimizations: ✨

Nothing. 💀

✨ Author ✨

Martin Koehler, MPI Magdeburg 🫶

✨ Date ✨

October 2025 🌸

Definition at line 208 of file sla_tgsylv_l2_rowwise.f90.

✨ subroutine sla_tgsylv_l3 ✨

subroutine sla_tgsylv_l3 (character, dimension(1) transa, character, dimension(1) transb, real sgn, integer m, integer n, real, dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d, integer ldd, real, dimension(ldx, *) x, integer ldx, real scale, real, dimension(*) work, integer info)

Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation. 🔥

✨ Purpose: ✨

okay so SLA_TGSYLV_L3 is out here solving a generalized Sylvester equation and it's giving these two forms:

op1(A) * X * op2(B) + op1(C) * X * op2(D) = SCALE * Y                              (1)

op1(A) * X * op2(B) - op1(C) * X * op2(D) = SCALE * Y                              (2)

where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix and B and D are N-by-N upper triangular matrices. Furthermore either B or D can be quasi upper triangular as well. The right hand side Y and the solution X are M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) ( or (D,B) ) are created by SGGES from LAPACK. 💖

✨ Attention ✨

The algorithm is implemented using BLAS level 3 operations. No biggie, just level 3 realness. 💅

🎀 Remarks 🎀

So the algorithm works without STRMM calls because of performance issues. All STRMM calls are replaced by block-GEMM operations. Therefore it is necessary that A and C contain zeros below the diagonal. 🫶

👑 Parameters 👑

TRANSA — TRANSA is CHARACTER(1). Specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)
TRANSB — TRANSB is CHARACTER(1). Specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)
SGN — SGN is REAL, allowed values: +/-1. Specifies the sign between the two parts of the Sylvester equation.
- = 1 : Solve Equation (1)
- == -1: Solve Equation (2)
M — M is INTEGER. The order of the matrices A and C. M >= 0.
N — N is INTEGER. The order of the matrices B and D. N >= 0.
A — A is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular.
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — B is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular the matrix B must be upper triangular.
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — C is REAL array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular.
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
X — X is REAL array, dimension (LDX,N). On input, the matrix X contains the right hand side Y. On output, the matrix X contains the solution of Equation (1) or (2) as selected by TRANSA, TRANSB, and SGN. Right hand side Y and the solution X are M-by-N matrices.
LDX — LDX is INTEGER. The leading dimension of the array X. LDX >= max(1,M).
SCALE — SCALE is REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — WORK is REAL array, dimension LWORK. Workspace for the algorithm. The workspace needs to queried before the running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO.
INFO — INFO is INTEGER.

On input:

== -1 : Perform a workspace query
<> -1: normal operation

On exit, workspace query:

< 0 : if INFO = -i, the i-th argument had an illegal value
>= 0: The value of INFO is the required number of elements in the workspace.

On exit, normal operation:

== 0: successful exit
< 0: if INFO = -i, the i-th argument had an illegal value
> 0: The equation is not solved correctly. One of the arising inner system got singular.

🫶 Author 🫶

Martin Koehler, MPI Magdeburg

📅 Date 📅

October 2025

Definition at line 205 of file sla_tgsylv_l3_opt.f90.

Next up, bestie, we got another subroutine: sla_tgsylv_l3_2s. Here's the signature:

   subroutine sla_tgsylv_l3_2s (character, dimension(1) transa, character,
       dimension(1) transb, real sgn, integer m, integer n, real,
       dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer
       ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d,
       integer ldd, real, dimension(ldx, *) x, integer ldx, real scale, real,
       dimension(*) work, integer info)

It's the Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation. 💅

Purpose:

💅 SLA_TGSYLV_L3_2S — what it actually solves

okay bestie, so SLA_TGSYLV_L3_2S solves a generalized Sylvester equation, and it comes in two flavours:

op1(A) * X * op2(B) + op1(C) * X * op2(D) = SCALE * Y                              (1)

op1(A) * X * op2(B) - op1(C) * X * op2(D) = SCALE * Y                              (2)

where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix, and B and D are N-by-N upper triangular matrices. and here's the cute lil twist: either B or D can be quasi upper triangular as well. ✨

the right hand side Y and the solution X are M-by-N matrices. and like, typically the matrix pencils (A,C) and (B,D) (or (D,B), same vibe) are created by SGGES from LAPACK. iykyk 🫶

⚠️ Attention

the algorithm is implemented using BLAS level 3 operations. that's it, that's the note. fast girl behaviour 💅

🎀 Parameters

TRANSA — TRANSA is CHARACTER(1). it specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (no transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C**T (transpose A and C)
TRANSB — TRANSB is CHARACTER(1). it specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (no transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D**T (transpose B and D)
SGN — SGN is REAL, allowed values: +/-1. it specifies the sign between the two parts of the Sylvester equation.
- = 1: solve Equation (1)
- == -1: solve Equation (2)
M — M is INTEGER. the order of the matrices A and C. M >= 0. ✨
N — N is INTEGER. the order of the matrices B and D. N >= 0. ✨
A — A is REAL array, dimension (LDA,M). the matrix A must be (quasi-) upper triangular.
LDA — LDA is INTEGER. the leading dimension of the array A. LDA >= max(1,M).
B — B is REAL array, dimension (LDB,N). the matrix B must be (quasi-) upper triangular. but listen — if the matrix D is already quasi-upper triangular, then the matrix B must be upper triangular. no cap.
LDB — LDB is INTEGER. the leading dimension of the array B. LDB >= max(1,N).
C — C is REAL array, dimension (LDC,M). the matrix C must be upper triangular.
LDC — LDC is INTEGER. the leading dimension of the array C. LDC >= max(1,M).
D — D is REAL array, dimension (LDD,N). the matrix D must be (quasi-) upper triangular. and same energy as before: if the matrix B is already quasi-upper triangular, the matrix D must be upper triangular. 💀
LDD — LDD is INTEGER. the leading dimension of the array D. LDD >= max(1,N).
X — X is REAL array, dimension (LDX,N). on input, the matrix X contains the right hand side Y. on output, the matrix X contains the solution of Equation (1) or (2) as selected by TRANSA, TRANSB, and SGN. right hand side Y and the solution X are M-by-N matrices. main character matrix fr fr 👑
LDX — LDX is INTEGER. the leading dimension of the array X. LDX >= max(1,M).
SCALE — SCALE is REAL. SCALE is a scaling factor to prevent the overflow in the result. if INFO == 0 then SCALE is 1.0 — otherwise, if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — WORK is REAL array, dimension LWORK. workspace for the algorithm. 🫶 and bestie, the workspace needs to be queried before the running the computation. the query is performed by calling the subroutine with INFO == -1 on input. the required workspace is then returned in INFO.
INFO — INFO is INTEGER.

on input:

== -1: perform a workspace query
<> -1: normal operation

on exit, workspace query:

< 0: if INFO = -i, the i-th argument had an illegal value
>= 0: the value of INFO is the required number of elements in the workspace

on exit, normal operation:

== 0: successful exit 🎉
< 0: if INFO = -i, the i-th argument had an illegal value
> 0: the equation is not solved correctly. one of the arising inner system got singular. 💀

🌸 Author

Martin Koehler, MPI Magdeburg

📅 Date

October 2025

🔍 Definition

Definition at line 202 of file sla_tgsylv_l3_2stage.f90.

🦋 subroutine sla_tgsylv_l3_colwise

subroutine sla_tgsylv_l3_colwise (character, dimension(1) transa, character,
    dimension(1) transb, real sgn, integer m, integer n, real,
    dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer
    ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d,
    integer ldd, real, dimension(ldx, *) x, integer ldx, real scale, real,
    dimension(*) work, integer info)

Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation.

Purpose:

okay so like, what is this one actually for? pull up a chair bestie, we're about to find out —

okay bestie, so SLA_TGSYLV_L3_COLWISE solves a generalized Sylvester equation, and it comes in exactly two flavours, no more no less:

op1(A) * X * op2(B) + op1(C) * X * op2(D) = SCALE * Y                              (1)

op1(A) * X * op2(B) - op1(C) * X * op2(D) = SCALE * Y                              (2)

where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix, and B and D are N-by-N upper triangular matrices. Also, plot twist: either B or D can be quasi upper triangular as well. The right hand side Y and the solution X are M-by-N matrices. And typically — like, usually — the matrix pencils (A,C) and (B,D) (or (D,B)) are created by SGGES from LAPACK. iykyk 💅

⚠️ Attention ⚠️

The algorithm is implemented using BLAS level 3 operations. periodt.

🧐 Remarks 🧐

ngl, real talk: the algorithm works without STRMM calls because of performance issues. All STRMM calls are replaced by block-GEMM operations. So therefore it is necessary that A and C contain zeros below the diagonal. No zeros below the diagonal, no slay. 🫶

✨ Parameters ✨

🔤 `TRANSA`

TRANSA is CHARACTER(1). It specifies the form of the system of equations with respect to A and C:

== 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
== 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)

🔤 `TRANSB`

TRANSB is CHARACTER(1). It specifies the form of the system of equations with respect to B and D:

== 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
== 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)

➕➖ `SGN`

SGN is REAL, allowed values: +/-1. It specifies the sign between the two parts of the Sylvester equation.

= 1: Solve Equation (1)
== -1: Solve Equation (2)

📏 `M`

M is INTEGER. The order of the matrices A and C. M >= 0.

📏 `N`

N is INTEGER. The order of the matrices B and D. N >= 0.

🅰️ `A`

A is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular.

📐 `LDA`

LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).

🅱️ `B`

B is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular the matrix B must be upper triangular. 👀

📐 `LDB`

LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).

🅾️ `C`

C is REAL array, dimension (LDC,M). The matrix C must be upper triangular.

📐 `LDC`

LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).

🆔 `D`

D is REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular. (see? it's the mirror of the B situation, they're matching 💞)

📐 `LDD`

LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).

❌ `X`

X is REAL array, dimension (LDX,N). On input, the matrix X contains the right hand side Y. On output, the matrix X contains the solution of Equation (1) or (2) as selected by TRANSA, TRANSB, and SGN. Right hand side Y and the solution X are M-by-N matrices. she does both jobs, iconic behavior 👑

📐 `LDX`

LDX is INTEGER. The leading dimension of the array X. LDX >= max(1,M).

⚖️ `SCALE`

SCALE is REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true. she's out here protecting your numbers from overflowing, a true bodyguard 🛡️

🛠️ `WORK`

WORK is REAL array, dimension LWORK. Workspace for the algorithm. You gotta query the workspace before running the computation, no skipping. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO.

ℹ️ `INFO`

INFO is INTEGER.

On input:

== -1 : Perform a workspace query
<> -1: normal operation

On exit, workspace query:

< 0 : if INFO = -i, the i-th argument had an illegal value
>= 0: The value of INFO is the required number of elements in the workspace.

On exit, normal operation:

== 0: successful exit
< 0: if INFO = -i, the i-th argument had an illegal value
> 0: The equation is not solved correctly. One of the arising inner system got singular.

👑 Author 👑

Martin Koehler, MPI Magdeburg

📅 Date 📅

October 2025

📍 Definition 📍

Definition at line 204 of file sla_tgsylv_l3_colwise.f90.

alright bestie, new subroutine, new me ✨ meet the rowwise twin:

subroutine sla_tgsylv_l3_rowwise (character, dimension(1) transa, character,
    dimension(1) transb, real sgn, integer m, integer n, real,
    dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer
    ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d,
    integer ldd, real, dimension(ldx, *) x, integer ldx, real scale, real,
    dimension(*) work, integer info)

Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation. 🔥

🎯 Purpose 🎯

Purpose:

!>
!> SLA_TGSYLV_L3_ROWWISE solves a generalized Sylvester equation of the following forms
!>
!>    op1(A) * X * op2(B) + op1(C) * X * op2(D) = SCALE * Y                              (1)
!>
!> or
!>
!>    op1(A) * X * op2(B) - op1(C) * X * op2(D) = SCALE * Y                              (2)
!>
!> where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular
!> matrix and B and D are N-by-N upper triangular matrices. Furthermore either B or
!> D can be quasi upper triangular as well. The right hand side Y and the solution X
!> M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) ( or (D,B) ) are
!> created by SGGES from LAPACK.
!>

Attention 💅

The algorithm is implemented using BLAS level 3 operations. That's the whole vibe, bestie. No cap.

Parameters 🫶

TRANSA — TRANSA is CHARACTER(1). Specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)
TRANSB — TRANSB is CHARACTER(1). Specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)
SGN — SGN is REAL, allowed values: +/-1. Specifies the sign between the two parts of the Sylvester equation.
- = 1 : Solve Equation (1)
- == -1: Solve Equation (2)
M — M is INTEGER. The order of the matrices A and C. M >= 0.
N — N is INTEGER. The order of the matrices B and D. N >= 0.
A — A is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular.
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — B is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular the matrix B must be upper triangular.
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — C is REAL array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular.
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
X — X is REAL array, dimension (LDX,N). On input, the matrix X contains the right hand side Y. On output, the matrix X contains the solution of Equation (1) or (2) as selected by TRANSA, TRANSB, and SGN. Right hand side Y and the solution X are M-by-N matrices.
LDX — LDX is INTEGER. The leading dimension of the array X. LDX >= max(1,M).
SCALE — SCALE is REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — WORK is REAL array, dimension LWORK. Workspace for the algorithm. The workspace needs to queried before the running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO.
INFO — INFO is INTEGER.
- On input:
  - == -1 : Perform a workspace query
  - <> -1: normal operation
- On exit, workspace query:
  - < 0 : if INFO = -i, the i-th argument had an illegal value
  - >= 0: The value of INFO is the required number of elements in the workspace.
- On exit, normal operation:
  - == 0: successful exit
  - < 0: if INFO = -i, the i-th argument had an illegal value
  - > 0: The equation is not solved correctly. One of the arising inner system got singular.

Author 🎀

Martin Koehler, MPI Magdeburg

Date 🌸

October 2025

Definition at line ... 💖

Definition at line 203 of file sla_tgsylv_l3_rowwise.f90.

🦋 `sla_tgsylv_recursive` — Level-3 Recursive Blocked Solver 🦋

recursive subroutine sla_tgsylv_recursive (character, dimension(1) transa,
    character, dimension(1) transb, real sgn, integer m, integer n, real,
    dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer
    ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d,
    integer ldd, real, dimension(ldx, *) x, integer ldx, real scale, real,
    dimension(*) work, integer info)

Level-3 Recursive Blocked Solver for the generalized Sylvester equation.

Purpose:

!>
!> SLA_TGSYLV_RECURSIVE solves a generalized Sylvester equation of the following forms
!>
!>    op1(A) * X * op2(B) + op1(C) * X * op2(D) = SCALE * Y                              (1)
!>
!> or
!>
!>    op1(A) * X * op2(B) - op1(C) * X * op2(D) = SCALE * Y                              (2)
!>
!> where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular
!> matrix and B and D are N-by-N upper triangular matrices. Furthermore either B or
!> D can be quasi upper triangular as well. The right hand side E and the solution X
!> M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) ( or (D,B) ) are
!> created by SGGES from LAPACK.
!>

Attention 💅

The algorithm is a rewrite of RECSY with some optimizations.

Parameters 🫶

TRANSA — TRANSA is CHARACTER(1). Specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)
TRANSB — TRANSB is CHARACTER(1). Specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)
SGN — SGN is REAL, allowed values: +/-1. Specifies the sign between the two parts of the Sylvester equation.
- = 1 : Solve Equation (1)
- == -1: Solve Equation (2)
M — M is INTEGER. The order of the matrices A and C. M >= 0.
N — N is INTEGER. The order of the matrices B and D. N >= 0.
A — A is REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular.
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — B is REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular the matrix B must be upper triangular.
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — C is REAL array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is REAL array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular.
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
X — X is REAL array, dimension (LDX,N). On input, the matrix X contains the right hand side Y. On output, the matrix X contains the solution X of Equation (1) or (2) as selected by TRANSA, TRANSB, and SGN. Right hand side Y and the solution X are M-by-N matrices.
LDX — LDX is INTEGER. The leading dimension of the array X. LDX >= max(1,M).
SCALE — SCALE is REAL. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — WORK is REAL array, dimension M*N. Workspace for the algorithm
INFO — INFO is INTEGER
- == 0: successful exit
- < 0: if INFO = -i, the i-th argument had an illegal value
- > 0: The equation is not solved correctly. One of the arising inner system got singular.

Author 🎀

Martin Koehler, MPI Magdeburg

Date 🌸

October 2025

Definition at line ... 💖

Definition at line 190 of file sla_tgsylv_recursive.f90.

Author 🎀

Generated automatically by Doxygen for MEPACK from the source code.

And that's a wrap, bestie! 💖 This page is Version 1.1.2, dated Thu Oct 9 2025 00:00:00, and it's the sgltgsylv(3) manual. Slay! 💅✨

📦 where this page comes from

fedora 43: mepack-devel 1.1.2-1.fc43 · GPL-3.0-or-later and LGPL-2.1-or-later and BSD-3-Clause
fedora 43: mepack-devel 1.1.2-1.fc43 · GPL-3.0-or-later and LGPL-2.1-or-later and BSD-3-Clause
💖 girlified by deepseek-ai/DeepSeek-V4.1-Flash on 2026-10-14: a language model's rewrite of the original, under its terms

sgltgsylv(3)

💖 sgltgsylv(3) — MEPACK 💖

✨ NAME ✨

🌸 SYNOPSIS 🌸

🎀 Functions 🎀

sla_tgcsylv_dag

sla_tgcsylv_dual_dag

sla_tgsylv_dag

sla_tgcsylv_dual_kernel_44nn

sla_tgcsylv_dual_kernel_44nt

sla_tgcsylv_dual_kernel_44tn

sla_tgcsylv_dual_kernel_44tt

sla_tgcsylv_dual_l2

sla_tgcsylv_dual_l2_local_copy

sla_tgcsylv_dual_l2_local_copy_128

sla_tgcsylv_dual_l2_local_copy_32

sla_tgcsylv_dual_l2_local_copy_64

sla_tgcsylv_dual_l2_local_copy_96

sla_tgcsylv_l2

sla_tgcsylv_l2_local_copy

sla_tgcsylv_l2_local_copy_128

sla_tgcsylv_l2_local_copy_32

sla_tgcsylv_l2_local_copy_64

sla_tgcsylv_l2_local_copy_96

sla_tgcsylv_l2_reorder

sla_tgcsylv_l2_unopt

sla_tgsylv_gardiner_laub

sla_tgsylv_l2

sla_tgsylv_l2_colwise

sla_tgsylv_l2_local_copy

sla_tgsylv_l2_local_copy_128

sla_tgsylv_l2_local_copy_32

sla_tgsylv_l2_local_copy_64

sla_tgsylv_l2_local_copy_96

sla_tgsylv_l2_reorder

sla_tgsylv_l2_rowwise

sla_tgcsylv_dual_l3

sla_tgcsylv_dual_l3_2s

sla_tgcsylv_l3_2s

sla_tgcsylv_l3

sla_tgcsylv_l3_unopt

sla_tgsylv_l3_2s

sla_tgsylv_l3_colwise

sla_tgsylv_l3

sla_tgsylv_l3_rowwise

sla_tgcsylv_dual_recursive

sla_tgcsylv_recursive

sla_tgsylv_recursive

Detailed Description 💖

✨ Function Documentation ✨

sla_tgcsylv_dag 🦋

Purpose 🎀

Remarks 🌸

Parameters 🫶

Author 👑

Date ⭐

Definition 📍

sla_tgcsylv_dual_dag 🧚

Purpose 🦋

Parameters 🫶

Author 🖋️

Date 📅

💖 subroutine sla_tgcsylv_dual_kernel_44nn 💖

Purpose

Attention ⚠️

Parameters

Author 🖋️

Date 📅

🦋 subroutine sla_tgcsylv_dual_kernel_44nt 🦋

Purpose

Attention (read this fr) 💀

Parameters 🧾

Author 👑

Date 📅

sla_tgcsylv_dual_kernel_44tn

Purpose 🎯

Attention (read this fr) 💀

Parameters 🧾

Author

Date

`sla_tgcsylv_dag`

`sla_tgcsylv_dual_dag`

`sla_tgsylv_dag`

`sla_tgcsylv_dual_kernel_44nn`

`sla_tgcsylv_dual_kernel_44nt`

`sla_tgcsylv_dual_kernel_44tn`

`sla_tgcsylv_dual_kernel_44tt`

`sla_tgcsylv_dual_l2`

`sla_tgcsylv_dual_l2_local_copy`

`sla_tgcsylv_dual_l2_local_copy_128`

`sla_tgcsylv_dual_l2_local_copy_32`

`sla_tgcsylv_dual_l2_local_copy_64`

`sla_tgcsylv_dual_l2_local_copy_96`

`sla_tgcsylv_l2`

`sla_tgcsylv_l2_local_copy`

`sla_tgcsylv_l2_local_copy_128`

`sla_tgcsylv_l2_local_copy_32`

`sla_tgcsylv_l2_local_copy_64`

`sla_tgcsylv_l2_local_copy_96`

`sla_tgcsylv_l2_reorder`

`sla_tgcsylv_l2_unopt`

`sla_tgsylv_gardiner_laub`

`sla_tgsylv_l2`

`sla_tgsylv_l2_colwise`

`sla_tgsylv_l2_local_copy`

`sla_tgsylv_l2_local_copy_128`

`sla_tgsylv_l2_local_copy_32`

`sla_tgsylv_l2_local_copy_64`

`sla_tgsylv_l2_local_copy_96`

`sla_tgsylv_l2_reorder`

`sla_tgsylv_l2_rowwise`

`sla_tgcsylv_dual_l3`

`sla_tgcsylv_dual_l3_2s`

`sla_tgcsylv_l3_2s`

`sla_tgcsylv_l3`

`sla_tgcsylv_l3_unopt`

`sla_tgsylv_l3_2s`

`sla_tgsylv_l3_colwise`

`sla_tgsylv_l3`

`sla_tgsylv_l3_rowwise`

`sla_tgcsylv_dual_recursive`

`sla_tgcsylv_recursive`

`sla_tgsylv_recursive`

`sla_tgcsylv_dag` 🦋

`sla_tgcsylv_dual_dag` 🧚

💖 subroutine `sla_tgcsylv_dual_kernel_44nn` 💖

🦋 subroutine `sla_tgcsylv_dual_kernel_44nt` 🦋

`sla_tgcsylv_dual_kernel_44tn`

👑 subroutine `sla_tgcsylv_dual_kernel_44tt`

🌸 subroutine `sla_tgcsylv_dual_l2`

✨ `sla_tgcsylv_dual_l2_local_copy_128` ✨

🌈 `sla_tgcsylv_dual_l2_local_copy_64`

🦋 `sla_tgcsylv_l2_local_copy_128` 🦋