sgltrsylv(3)

✨ sgltrsylv ✨ — straight out of the MEPACK manual, section 3 💅

okay so bestie, pull up a chair 🪑 because this is the sgltrsylv(3) page and it is long. This is part 1 of 17, so we're gonna take it slow and cute. Let's go 🦋

---

🌸 NAME 🌸

sgltrsylv — Single Precision

• Single Precision routines for triangular standard Sylvester equations.

---

🎀 SYNOPSIS 🎀

✨ Functions ✨

• sla_trsylv2_dag — Level-3 Bartels-Stewart Algorithm for the discrete time Sylvester equation with DAG based parallelization.

subroutine sla_trsylv2_dag (transa, transb, sgn, m, n, a, lda, b, ldb, x,
    ldx, scale, work, info)

• sla_trsylv_dag — Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation with DAG parallelization.

subroutine sla_trsylv_dag (transa, transb, sgn, m, n, a, lda, b, ldb, x,
    ldx, scale, work, info)

• sla_trsylv2_kernel_44nn — Solver for a 4x4 discrete time Sylvester equation (TRANSA = N, TRANSB = N)

subroutine sla_trsylv2_kernel_44nn (sgn, m, n, a, lda, b, ldb, x, ldx,
    scale, info)

• sla_trsylv2_kernel_44nt — Solver for a 4x4 discrete time Sylvester equation (TRANSA = N, TRANSB = T)

subroutine sla_trsylv2_kernel_44nt (sgn, m, n, a, lda, b, ldb, x, ldx,
    scale, info)

• sla_trsylv2_kernel_44tn — Solver for a 4x4 discrete time Sylvester equation (TRANSA = T, TRANSB = N)

subroutine sla_trsylv2_kernel_44tn (sgn, m, n, a, lda, b, ldb, x, ldx,
    scale, info)

• sla_trsylv2_kernel_44tt — Solver for a 4x4 discrete time Sylvester equation (TRANSA = T, TRANSB = T)

subroutine sla_trsylv2_kernel_44tt (sgn, m, n, a, lda, b, ldb, x, ldx,
    scale, info)

• sla_trsylv_kernel_44nn — Solver for a 4x4 Sylvester equation (TRANSA = N, TRANSB = N)

subroutine sla_trsylv_kernel_44nn (sgn, m, n, a, lda, b, ldb, x, ldx,
    scale, info)

• sla_trsylv_kernel_44nt — Solver for a 4x4 Sylvester equation (TRANSA = N, TRANSB = T)

subroutine sla_trsylv_kernel_44nt (sgn, m, n, a, lda, b, ldb, x, ldx,
    scale, info)

• sla_trsylv_kernel_44tn — Solver for a 4x4 Sylvester equation (TRANSA = T, TRANSB = N)

subroutine sla_trsylv_kernel_44tn (sgn, m, n, a, lda, b, ldb, x, ldx,
    scale, info)

• sla_trsylv_kernel_44tt — Solver for a 4x4 Sylvester equation (TRANSA = T, TRANSB = T)

subroutine sla_trsylv_kernel_44tt (sgn, m, n, a, lda, b, ldb, x, ldx,
    scale, info)

• sla_trsylv2_l2 — Level-2 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized)

subroutine sla_trsylv2_l2 (transa, transb, sgn, m, n, a, lda, b, ldb, x,
    ldx, scale, work, info)

• sla_trsylv2_l2_local_copy — Level-2 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized)

subroutine sla_trsylv2_l2_local_copy (transa, transb, sgn, m, n, a, lda,
    b, ldb, x, ldx, scale, work, info)

• sla_trsylv2_l2_local_copy_128 — Level-2 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized)

subroutine sla_trsylv2_l2_local_copy_128 (transa, transb, sgn, m, n, a,
    lda, b, ldb, x, ldx, scale, work, info)

• sla_trsylv2_l2_local_copy_32 — Level-2 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized)

subroutine sla_trsylv2_l2_local_copy_32 (transa, transb, sgn, m, n, a,
    lda, b, ldb, x, ldx, scale, work, info)

• sla_trsylv2_l2_local_copy_64 — Level-2 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized)

subroutine sla_trsylv2_l2_local_copy_64 (transa, transb, sgn, m, n, a,
    lda, b, ldb, x, ldx, scale, work, info)

• sla_trsylv2_l2_local_copy_96 — Level-2 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized)

subroutine sla_trsylv2_l2_local_copy_96 (transa, transb, sgn, m, n, a,
    lda, b, ldb, x, ldx, scale, work, info)

• sla_trsylv2_l2_reorder — Level-2 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized)

subroutine sla_trsylv2_l2_reorder (transa, transb, sgn, m, n, a, lda, b,
    ldb, x, ldx, scale, work, info)

• sla_trsylv2_l2_unopt — Level-2 Bartels-Stewart Algorithm for the discrete time Sylvester equation (unoptimized)

subroutine sla_trsylv2_l2_unopt (transa, transb, sgn, m, n, a, lda, b,
    ldb, x, ldx, scale, work, info)

• sla_trsylv_l2 — Level-2 Bartels-Stewart Algorithm for the Sylvester equation (optimized)

subroutine sla_trsylv_l2 (transa, transb, sgn, m, n, a, lda, b, ldb, x,
    ldx, scale, work, info)

• sla_trsylv_l2_local_copy — Level-2 Bartels-Stewart Algorithm for the Sylvester equation (optimized)

subroutine sla_trsylv_l2_local_copy (transa, transb, sgn, m, n, a, lda,
    b, ldb, x, ldx, scale, work, info)

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

subroutine sla_trsylv_l2_local_copy_128 (transa, transb, sgn, m, n, a,
    lda, b, ldb, x, ldx, scale, work, info)

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

subroutine sla_trsylv_l2_local_copy_32 (transa, transb, sgn, m, n, a,
    lda, b, ldb, x, ldx, scale, work, info)

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

subroutine sla_trsylv_l2_local_copy_64 (transa, transb, sgn, m, n, a,
    lda, b, ldb, x, ldx, scale, work, info)

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

subroutine sla_trsylv_l2_local_copy_96 (transa, transb, sgn, m, n, a,
    lda, b, ldb, x, ldx, scale, work, info)

• sla_trsylv_l2_reorder — Level-2 Bartels-Stewart Algorithm for the Sylvester equation (optimized)

subroutine sla_trsylv_l2_reorder (transa, transb, sgn, m, n, a, lda, b,
    ldb, x, ldx, scale, work, info)

• sla_trsylv_l2_unopt — Level-2 Bartels-Stewart Algorithm for the Sylvester equation (unoptimized)

subroutine sla_trsylv_l2_unopt (transa, transb, sgn, m, n, a, lda, b,
    ldb, x, ldx, scale, work, info)

• sla_trsylv2_l3_2s — Two level blocked Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized)

subroutine sla_trsylv2_l3_2s (transa, transb, sgn, m, n, a, lda, b, ldb,
    x, ldx, scale, work, info)

• sla_trsylv2_l3 — Level-3 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized)

subroutine sla_trsylv2_l3 (transa, transb, sgn, m, n, a, lda, b, ldb, x,
    ldx, scale, work, info)

• sla_trsylv2_l3_unopt — Level-3 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Not Optimized)

subroutine sla_trsylv2_l3_unopt (transa, transb, sgn, m, n, a, lda, b,
    ldb, x, ldx, scale, work, info)

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

subroutine sla_trsylv_l3_2s (transa, transb, sgn, m, n, a, lda, b, ldb,
    x, ldx, scale, work, info)

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

subroutine sla_trsylv_l3 (transa, transb, sgn, m, n, a, lda, b, ldb, x,
    ldx, scale, work, info)

• sla_trsylv_l3_unopt — Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation (unoptimized)

subroutine sla_trsylv_l3_unopt (transa, transb, sgn, m, n, a, lda, b,
    ldb, x, ldx, scale, work, info)

• sla_trsylv2_recursive — Level-3 Recursive Blocked Algorithm for the discrete time Sylvester equation.

recursive subroutine sla_trsylv2_recursive (transa, transb, sgn, m, n,
    a, lda, b, ldb, x, ldx, scale, work, info)

• sla_trsylv_recursive — Level-3 Recursive Blocked Solver for the Sylvester equation.

recursive subroutine sla_trsylv_recursive (transa, transb, sgn, m, n, a,
    lda, b, ldb, x, ldx, scale, work, info)

---

🌈 Detailed Description 🌈

Single Precision routines for triangular standard Sylvester equations.

okay so this whole section? it's the solvers for the standard Sylvester equation with (quasi) triangular coefficient matrices. 💫

and like, where do those coefficient matrices come from? bestie, they're normally generated with the help of the Schur decomposition from LAPACK. that's the origin story 🫶

and every single code in here uses single precision arithmetic. no double, no mixing, just single. periodt. 💅

---

💖 Function Documentation 💖

🦋 subroutine sla_trsylv2_dag 🦋

subroutine sla_trsylv2_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(ldx, n) x, integer ldx, real scale, real,
    dimension(*) work, integer info)

Level-3 Bartels-Stewart Algorithm for the discrete time Sylvester equation with DAG based parallelization.

✨ Purpose ✨

!>
!> SLA_TRSYLV2_DAG solves a discrete time Sylvester equation of the following forms
!>
!>    op1(A) * X * op2(B) +  X  = SCALE * Y                              (1)
!>
!> or
!>
!>    op1(A) * X * op2(B) -  X  = SCALE * Y                              (2)
!>
!> where A is a M-by-M quasi upper triangular matrix and B is a N-by-N upper
!> quasi triangular matrices. The right hand side Y and the solution X
!> M-by-N matrices. Typically the matrices A and B are created by SGEES from LAPACK.
!>

so basically: A is a M-by-M quasi upper triangular matrix, B is a N-by-N upper quasi triangular matrix (yes the original says "matrices", we don't question it 😌), and the right hand side Y and the solution X are M-by-N matrices. and the glow-up moment? A and B are typically created by SGEES from LAPACK. iconic. 👑

🌸 Remarks 🌸

The algorithm is implemented using BLAS level 3 operations and OpenMP 4.0 DAG scheduling. yes, DAG scheduling, she's doing the most and we love her for it 💅

💀 Attention 💀

bestie, listen. due to the parallel nature of the algorithm, the scaling is not applied to the right hand. so if the problem is ill-posed and scaling shows up, you have to solve the equation again with a solver that has complete scaling support — like the SLA_TRSYLV2_L3 routine. no shortcuts, no cap. 🚨

🎀 Parameters 🎀

• TRANSA — 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 — 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)
• SGN — REAL, and the allowed values are +/-1. it specifies the sign between the two parts of the Sylvester equation:
  • = 1: solve Equation (1)
  • == -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.
• 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. heads up: if the matrix D is already quasi-upper triangular, the matrix B must be upper triangular.
• LDB — INTEGER. the leading dimension of the array B. LDB >= 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. 🦋
• LDX — it's an INTEGER, bestie, and it's the leading dimension of the array X. Rule: LDX >= max(1,M). Don't test it. It WILL test you back. 💅
• SCALE — a REAL, and it's the little hero that stops your result from overflowing. If INFO == 0 then SCALE is 1.0. Otherwise, if one of the inner systems couldn't be solved correctly, then 0 < SCALE <= 1 holds true. And if SCALE .NE. 1? Then the problem is not solved correctly and you have to go use another solver. Not me finding this out at runtime 😭
• WORK — a REAL array, dimension LWORK. It's the workspace for the algorithm. You must query this workspace before running the computation, no skipping. You do the query by calling the subroutine with INFO == -1 on input. The required workspace then comes back to you in INFO. It's giving "ask before you take" energy 🫶
• INFO — an INTEGER, and it's the main character of diagnostics honestly 👑

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. 💀

✨ Author ✨

Martin Koehler, MPI Magdeburg

🎀 Date 🎀

October 2025

🦋 Definition 🦋

Definition at line 181 of file sla_trsylv2_dag.f90.

---

👑 subroutine sla_trsylv2_kernel_44nn

subroutine sla_trsylv2_kernel_44nn (real sgn, integer m, integer n, real,
    dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer
    ldb, real, dimension(ldx, *) x, integer ldx, real scale, integer info)

Solver for a 4x4 discrete time Sylvester equation (TRANSA = N, TRANSB = N). Okay so this one is specifically the 4x4 girlies, no cap.

💖 Purpose 💖

SLA_TRSYLV2_KERNEL_44NN solves a discrete time Sylvester equation that looks like this:

A * X * B + SGN * X  = SCALE * Y                              (1)

where A is a M-by-M quasi upper triangular matrix and B is a N-by-N quasi upper triangular matrix. The right hand side Y and the solution X are M-by-N matrices. Typically the matrices A and B are created by SGEES from LAPACK. The algorithm is implemented using BLAS level 2 operations. Thereby the order of M and N is at most 4. Furthermore, for fast execution the function does not check the input arguments. 🫣 (iconic of it, honestly, but you HAVE been warned.)

🌸 Parameters 🌸

• SGN — a REAL, allowed values: +/-1. It specifies the sign between the two parts of the Sylvester equation.
• M — an INTEGER, the order of the matrix A. Constraint: 4 >= M >= 0.
• N — an INTEGER, the order of the matrix B. Constraint: 4 >= N >= 0.
• A — a REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular. No exceptions bestie.
• LDA — an INTEGER, the leading dimension of the array A. Rule: LDA >= max(1,M).
• B — a REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular too. Matching set, we love to see it 💅
• LDB — an INTEGER, the leading dimension of the array B. Rule: LDB >= max(1,N).
• X — a 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) as selected by SGN. Right hand side Y and the solution X are M-by-N matrices. It's giving "same array, new identity" 🌈
• LDX — an INTEGER, the leading dimension of the array X. Rule: LDB >= max(1,M). (Yes, that says LDB. We're reporting it exactly as written, bestie, not our typo 😭)
• SCALE — a REAL, the 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 < SCAL <= 1 holds true. (And yes, SCAL — again, exactly as written 💀)
• INFO — an INTEGER.

On input:

• = 1 — skip the input data checks
• <> 1 — check input data like normal LAPACK like routines

On output:

• == 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 135 of file sla_trsylv2_kernel_44_nn.f90.

---

💜 subroutine sla_trsylv2_kernel_44nt

subroutine sla_trsylv2_kernel_44nt (real sgn, integer m, integer n, real,
    dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer
    ldb, real, dimension(ldx, *) x, integer ldx, real scale, integer info)

Solver for a 4x4 discrete time Sylvester equation (TRANSA = N, TRANSB = T). Same family, different transpose, we move 🌸

💖 Purpose 💖

SLA_TRSYLV2_KERNEL_44NT solves a discrete time Sylvester equation shaped like this:

A * X * B**T + SGN * X  = SCALE * Y                              (1)

where A is a M-by-M quasi upper triangular matrix and B is a N-by-N quasi upper triangular matrix. The right hand side Y and the solution X are M-by-N matrices. Typically the matrices A and B are created by SGEES from LAPACK. The algorithm is implemented using BLAS level 2 operations. Thereby the order of M and N is at most 4. Furthermore, for fast execution the function does not check the input arguments. 👀

🌸 Parameters 🌸

• SGN — a REAL, allowed values: +/-1. It specifies the sign between the two parts of the Sylvester equation.
• M — an INTEGER, the order of the matrix A. Constraint: 4 >= M >= 0.
• N — an INTEGER, the order of the matrix B. Constraint: 4 >= N >= 0.
• A — a REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular.
• LDA — an INTEGER, the leading dimension of the array A. Rule: LDA >= max(1,M).
• B — a REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular.
• LDB — an INTEGER, the leading dimension of the array B. Rule: LDB >= max(1,N).
• X — a 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) as selected by SGN. Right hand side Y and the solution X are M-by-N matrices.
• LDX — an INTEGER, the leading dimension of the array X. Rule: LDB >= max(1,M). (Still LDB. Still reporting faithfully. 🫶)
• SCALE — a REAL, the 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 < SCAL <= 1 holds true.
• INFO — an INTEGER.

On input:

• = 1 — skip the input data checks
• <> 1 — check input data like normal LAPACK like routines

On output:

• == 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 134 of file sla_trsylv2_kernel_44_nt.f90.

---

🍓 subroutine sla_trsylv2_kernel_44tn

subroutine sla_trsylv2_kernel_44tn (real sgn, integer m, integer n, real,
    dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer
    ldb, real, dimension(ldx, *) x, integer ldx, real scale, integer info)

Solver for a 4x4 discrete time Sylvester equation (TRANSA = T, TRANSB = N). Transpose moved to the other side, she's flexible like that 💅

💖 Purpose 💖

SLA_TRSYLV2_KERNEL_44TN solves a discrete time Sylvester equation of this exact form:

A**T * X * B + SGN * X  = SCALE * Y                              (1)

where A is a M-by-M quasi upper triangular matrix and B is a N-by-N quasi upper triangular matrix. The right hand side Y and the solution X are M-by-N matrices. Typically the matrices A and B are created by SGEES from LAPACK. The algorithm is implemented using BLAS level 2 operations. Thereby the order of M and N is at most 4. Furthermore, for fast execution the function does not check the input arguments. (Same disclaimer, same energy, still true 😤)

🌸 Parameters 🌸

• SGN — a REAL, allowed values: +/-1. It specifies the sign between the two parts of the Sylvester equation.
• M — an INTEGER, the order of the matrix A. Constraint: 4 >= M >= 0.
• N — an INTEGER, the order of the matrix B. Constraint: 4 >= N >= 0.
• A — a REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular.
• LDA — an INTEGER, the leading dimension of the array A. Rule: LDA >= max(1,M).
• B — a REAL array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular.
• LDB — an INTEGER, the leading dimension of the array B. Rule: LDB >= max(1,N).
• X — a 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) as selected by SGN. Right hand side Y and the solution X are M-by-N matrices.
• LDX — an INTEGER, the leading dimension of the array X. Rule: LDB >= max(1,M). (You already know the drill 💀)
• LDX — LDX is INTEGER. It's the leading dimension of the array X. And bestie, the page says LDB >= max(1,M) here — yes, it says LDB, not LDX, we're keeping it exactly as written, no cap. 👑
• SCALE — SCALE is REAL. It's a scaling factor whose whole job is to stop the result from overflowing. If INFO == 0 then SCALE is 1.0, otherwise — if one of the inner systems couldn't be solved correctly — 0 < SCAL <= 1 holds true. 💅
• INFO — INFO is INTEGER. On input:
  • = 1 : skip the input data checks ✨
  • <> 1: check input data like normal LAPACK like routines do.

On output:

• == 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 135 of file sla_trsylv2_kernel_44_tn.f90.

---

🌸 subroutine sla_trsylv2_kernel_44tt

subroutine sla_trsylv2_kernel_44tt (real sgn, integer m, integer n, real,
    dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer
    ldb, real, dimension(ldx, *) x, integer ldx, real scale, integer info)

Solver for a 4x4 discrete time Sylvester equation (TRANSA = T, TRANSB = T) 💖

Purpose ✨

Okay so SLA_TRSYLV2_KERNEL_44TT solves a discrete time Sylvester equation that looks like this:

A**T * X * B**T + SGN * X  = SCALE * Y                              (1)

Where A is a M-by-M quasi upper triangular matrix and B is a N-by-N quasi upper triangular matrix. The right hand side Y and the solution X are M-by-N matrices. Typically the matrices A and B are created by SGEES from LAPACK, iykyk. 👑

The algorithm is implemented using BLAS level 2 operations. Thereby the order of M and N is at most 4. And — this is the part you gotta know, bestie — for fast execution, the function does not check the input arguments. It just goes. 🏃‍♀️💨

Parameters 🌈

• SGN — SGN is REAL, allowed values: +/-1. It specifies the sign between the two parts of the Sylvester equation. 💫
• M — M is INTEGER. The order of the matrix A. 4 >= M >= 0.
• N — N is INTEGER. The order of the matrix B. 4 >= 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.
• LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= 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) as selected by 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. And again, the page literally says LDB >= max(1,M) — we are NOT correcting it, we keep it fr fr. 🫶
• 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 < SCAL <= 1 holds true.
• INFO — INFO is INTEGER. On input:
  • = 1 : skip the input data checks ✨
  • <> 1: check input data like normal LAPACK like routines.

On output:

• == 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 135 of file sla_trsylv2_kernel_44_tt.f90.

---

💜 subroutine sla_trsylv2_l2

subroutine sla_trsylv2_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(ldx, *) x, integer ldx, real scale, real,
    dimension(*) work, integer info)

Level-2 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized) 🔥

Purpose ✨

So SLA_TRSYLV2_L2 solves a generalized Sylvester equation, and it comes in two flavors:

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

or

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

Where A is a M-by-M quasi upper triangular matrix and B is a N-by-N upper quasi triangular matrix. The right hand side Y and the solution X are M-by-N matrices. And yep, typically the matrices A and B are created by SGEES from LAPACK. 👑

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. 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).
• 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. And once again the page says LDB >= max(1,M), and once again we are keeping it exactly. Not me fixing a manpage, never. 🙈
• 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 STRMV and SGER by Fortran intrinsic,

Author 🎀

Martin Koehler, MPI Magdeburg

Date 🦋

October 2025

Definition 📌

Definition at line 182 of file sla_trsylv2_l2.f90.

---

🦋 subroutine sla_trsylv2_l2_local_copy

subroutine sla_trsylv2_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(ldx, *) x, integer ldx, real scale, real,
    dimension(*) work, integer info)

Level-2 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized) 💖

Purpose ✨

✨ sla_trsylv2_l2_local_copy ✨

okay bestie, so SLA_TRSYLV2_L2_LOCAL_COPY is out here solving a generalized Sylvester equation and honestly it's giving main character energy. 💅 It handles two forms:

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

or

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

where A is a M-by-M quasi upper triangular matrix and B is a N-by-N upper quasi triangular matrix. The right hand side Y and the solution X are M-by-N matrices. Typically the matrices A and B are created by SGEES from LAPACK. ✨

Parameters 🎀

• TRANSA — 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 — 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 — REAL, allowed values: +/-1. Specifies the sign between the two parts of the Sylvester equation.
  • = 1: Solve Equation (1)
  • == -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.
• 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 the matrix B must be upper triangular.
• LDB — INTEGER. The leading dimension of the array B. LDB >= 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.
• LDX — INTEGER. The leading dimension of the array X. LDB >= max(1,M). (yes, it says LDB, not LDX, we're just the messenger 💀)
• 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 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 — 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) .

Author — Martin Koehler, MPI Magdeburg Date — October 2025 Definition at line 184 of file sla_trsylv2_l2_opt_local_copy.f90.

---

✨ sla_trsylv2_l2_local_copy_128 ✨

subroutine sla_trsylv2_l2_local_copy_128 ( transa,  transb,  sgn,  m,  n,
    a,  lda,  b,  ldb,  x,  ldx,  scale,  work,  info)

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

Purpose:

SLA_TRSYLV2_L2_LOCAL_COPY_128 solves a generalized Sylvester equation of the following forms

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

or

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

where A is a M-by-M quasi upper triangular matrix and B is a N-by-N upper quasi triangular matrices. The right hand side Y and the solution X M-by-N matrices. Typically the matrices A and B are created by SGEES from LAPACK.

Attention ⚠️

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

Parameters (same as above, bestie, but we're listing them again because rules are rules 💖)

• TRANSA — 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 — 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 — REAL, allowed values: +/-1. Specifies the sign between the two parts of the Sylvester equation.
  • = 1: Solve Equation (1)
  • == -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.
• 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 the matrix B must be upper triangular.
• LDB — INTEGER. The leading dimension of the array B. LDB >= 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.
• LDX — INTEGER. The leading dimension of the array X. LDB >= max(1,M). (again with the LDB, iconic typo energy 💀)
• 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 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 — 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.
• Align local copies and fix problem size to <= 128

Author — Martin Koehler, MPI Magdeburg Date — October 2025 Definition at line 184 of file sla_trsylv2_l2_opt_local_copy_128.f90.

---

✨ sla_trsylv2_l2_local_copy_32 ✨

subroutine sla_trsylv2_l2_local_copy_32 ( transa,  transb,  sgn,  m,  n,  a,
    lda,  b,  ldb,  x,  ldx,  scale,  work,  info)

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

Purpose:

okay so it solves a generalized Sylvester equation in one of these two flavors:

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

or

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

where A is a M-by-M quasi upper triangular matrix and B is a N-by-N upper quasi triangular matrix. The right hand side Y and the solution X are M-by-N matrices. Typically the matrices A and B are created by SGEES from LAPACK. iykyk 🦋

Attention 🌸 — The size of the Problem is limited by M,N <= 32

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. (not me clocking that conditional, but okay bestie, noted) 🎀
• LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= 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). (yes, that's literally what it says — copying it exactly, 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. she's just here to keep the numbers from going feral 💅
• 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. (workspace query first, bestie — always)
• 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.
• Align local copies and fix problem size to <= 32

Author — Martin Koehler, MPI Magdeburg

Date — October 2025

Definition at line 184 of file sla_trsylv2_l2_opt_local_copy_32.f90.

---

sla_trsylv2_l2_local_copy_64

subroutine sla_trsylv2_l2_local_copy_64 ( transa,  transb,  sgn,  m,  n,  a,
    lda,  b,  ldb,  x,  ldx,  scale,  work,  info)

Level-2 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized) ✨ the sixty-four era ✨

Purpose:

!>
!> SLA_TRSYLV2_L2_LOCAL_COPY_64 solves a generalized Sylvester equation of the following forms
!>
!>    op1(A) * X * op2(B) +  X  = SCALE * Y                              (1)
!>
!> or
!>
!>    op1(A) * X * op2(B) -  X  = SCALE * Y                              (2)
!>
!> where A is a M-by-M quasi upper triangular matrix and B is a N-by-N upper
!> quasi triangular matrices. The right hand side Y and the solution X
!> M-by-N matrices. Typically the matrices A and B are created by SGEES from LAPACK.
!>

okay so same two equations as before, same A, same B, same SGEES origin story:

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

or

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

where A is a M-by-M quasi upper triangular matrix and B is a N-by-N upper quasi triangular matrix. The right hand side Y and the solution X are M-by-N matrices. Typically the matrices A and B are created by SGEES from LAPACK. 💜

Attention 🦋 — The size of the Problem is limited by M,N <= 64

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. (same conditional, still clocking it 🎀)
• LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= 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). (again, verbatim, don't come for me 💀)
• 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. overflow protection, she's so considerate 💅
• 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.
• Align local copies and fix problem size to <= 64

Author — Martin Koehler, MPI Magdeburg

Date — October 2025

Definition at line 184 of file sla_trsylv2_l2_opt_local_copy_64.f90.

---

sla_trsylv2_l2_local_copy_96

subroutine sla_trsylv2_l2_local_copy_96 ( transa,  transb,  sgn,  m,  n,  a,
    lda,  b,  ldb,  x,  ldx,  scale,  work,  info)

Level-2 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized) ✨ and now we're in the ninety-six timeline ✨

Purpose:

so bestie, SLA_TRSYLV2_L2_LOCAL_COPY_96 solves a generalized Sylvester equation, and it shows up in two lil flavors:

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

or

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

where A is a M-by-M quasi upper triangular matrix and B is a N-by-N upper quasi triangular matrices. The right hand side Y and the solution X M-by-N matrices. Typically the matrices A and B are created by SGEES from LAPACK.

---

🚨 Attention 🚨

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

like, that's the whole vibe of the "96" in the name, ok? don't come for her 💅

---

🌸 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.
  • = 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. no exceptions bestie ✨
• 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).
• 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 really said "glow up in place" 💖
• LDX — LDX is INTEGER. The leading dimension of the array X. LDB >= max(1,M).
• 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.
• 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 gotta ask first, it's called ✨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:

• 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.
• Align local copies and fix problem size to <= 96

---

👑 Author

Martin Koehler, MPI Magdeburg

📅 Date

October 2025

📍 Definition

at line 184 of file sla_trsylv2_l2_opt_local_copy_96.f90.

---

subroutine sla_trsylv2_l2_reorder

subroutine sla_trsylv2_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(ldx, *) x, integer ldx, real scale, real,
    dimension(*) work, integer info)

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

🎯 Purpose:

so this one, SLA_TRSYLV2_L2_UNOPT, solves a generalized Sylvester equation of the following forms:

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

or

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

where A is a M-by-M quasi upper triangular matrix and B is a N-by-N upper quasi triangular matrices. The right hand side Y and the solution X M-by-N matrices. Typically the matrices A and B are created by SGEES from LAPACK.

ik the name says _UNOPT and the title says "Optimized" and that's a whole identity crisis but bestie that's literally what it says, I don't make the rules 💀

---

🌸 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).
• 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).
• 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,

---

👑 Author

Martin Koehler, MPI Magdeburg

📅 Date

October 2025

📍 Definition

at line 183 of file sla_trsylv2_l2_opt_reorder.f90.

---

subroutine sla_trsylv2_l2_unopt

subroutine sla_trsylv2_l2_unopt (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(ldx, *) x, integer ldx, real scale, real,
    dimension(*) work, integer info)

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

🎯 Purpose:

SLA_TRSYLV2_L2_UNOPT solves a generalized Sylvester equation of the following forms ✨

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

or

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

where A is a M-by-M quasi upper triangular matrix and B is a N-by-N upper quasi triangular matrices. The right hand side Y and the solution X M-by-N matrices. Typically the matrices A and B are created by SGEES from LAPACK. No cap 💅

---

Attention 😭

The algorithm is implemented using BLAS level 2 operations without further optimizations. For a faster implementation see SLA_TRSYLV2_L2. Bestie, this one is the no-frills version — it works, it's just not the glow-up 💀

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).
• 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).
• 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, and 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: 🌈

• Nothing.

Author 👑

Martin Koehler, MPI Magdeburg

Date 🗓️

October 2025

Definition at line 181 of file sla_trsylv2_l2_unopt.f90.

---

sla_trsylv2_l3 🦋

subroutine sla_trsylv2_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(ldx, *) x, integer ldx, real scale, real,
    dimension(*) work, integer info)

Level-3 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized) ✨

Purpose: 💜

SLA_TRSYLV2_L3 solves a generalized Sylvester equation of the following forms

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

or

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

where A is a M-by-M quasi upper triangular matrix and B is a N-by-N upper quasi triangular matrices. The right hand side Y and the solution X M-by-N matrices. Typically the matrices A and B are created by SGEES from LAPACK. Same energy as before, just the optimized flavor 💅

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).
• 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).
• 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 — you do that by calling the subroutine with INFO == -1 on input, and the required workspace comes back 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 173 of file sla_trsylv2_l3_opt.f90.

---

sla_trsylv2_l3_2s 🌸

subroutine sla_trsylv2_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(ldx, *) x, integer ldx, real scale, real,
    dimension(*) work, integer info)

Two level blocked Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized) 🔥

Purpose: 💜

...it solves a generalized Sylvester equation of these two iconic shapes:

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

or

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

where A is a M-by-M quasi upper triangular matrix and B is a N-by-N upper quasi triangular matrix. The right hand side Y and the solution X are M-by-N matrices. Typically the matrices A and B are created by SGEES from LAPACK. iykyk ✨

🎀 Remarks

The algorithm used level-3 BLAS operations and a DAG scheduled inner solver. lowkey iconic behavior.

💀 Attention

Because of the parallel nature of the inner solvers, the scaling is turned off and SCALE is set to ONE. No scaling era for this one, bestie — that's just how it is. 🫶

🌸 Parameters

• TRANSA — TRANSA is CHARACTER(1). It tells you 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 tells you 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, and the allowed values are +/-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. no exceptions, bestie.
• 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).
• 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. glow up fr fr 💅
• 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 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 act, periodt.
• 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 174 of file sla_trsylv2_l3_2stage.f90.

---

🌈 sla_trsylv2_l3_unopt

subroutine sla_trsylv2_l3_unopt (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(ldx, *) x, integer ldx, real scale, real,
    dimension(*) work, integer info)

Level-3 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Not Optimized). yes bestie, the not optimized one — she's doing her best and honestly that's valid. 🫶

Purpose

SLA_TRSYLV2_L3_UNOPT solves a generalized Sylvester equation of the following forms:

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

or

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

where A is a M-by-M quasi upper triangular matrix and B is a N-by-N upper quasi triangular matrix. The right hand side Y and the solution X are M-by-N matrices. Typically the matrices A and B are created by SGEES from LAPACK. ✨

💀 Attention

This function iterates column first over the result, and from this fact it will be a bit slower than SLA_TRSYLV2_L3. not a flaw, just a lifestyle. 😭

🌸 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).
• 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).
• 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. 💀

👑 Author

Martin Koehler, MPI Magdeburg

🦋 Date

October 2025

Definition at line 175 of file sla_trsylv2_l3_unopt.f90.

---

🧚 sla_trsylv2_recursive

recursive subroutine sla_trsylv2_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(ldx, *) x, integer ldx, real scale, real,
    dimension(*) work, integer info)

Level-3 Recursive Blocked Algorithm for the discrete time Sylvester equation.

Purpose

okay so SLA_TRSYLV2_RECURSIVE solves a generalized Sylvester equation that looks like this, bestie:

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

or

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

where A is a M-by-M quasi upper triangular matrix and B is a N-by-N upper quasi triangular matrix. The right hand side Y and the solution X are M-by-N matrices. Typically the matrices A and B are created by SGEES from LAPACK. iykyk 💅

🌸 Remarks 🌸

The algorithm uses recursive blocking instead of the Bartels-Stewart approach. it's giving main character energy, no cap 🫶

✨ 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).
• 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. LDB >= 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 M*N. Workspace for the algorithm.
• 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 163 of file sla_trsylv2_recursive.f90.

---

🦋 sla_trsylv_dag 🦋

subroutine sla_trsylv_dag (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(ldx, n) x, integer ldx, real scale, real,
    dimension(*) work, integer info)

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

🌈 Purpose 🌈

SLA_TRSYLV_DAG solves a Sylvester equation of the following forms:

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

or

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

where A is a M-by-M quasi upper triangular matrix, B is a N-by-N quasi upper triangular matrix. The right hand side Y and the solution X are M-by-N matrices. Typically the matrices A and B are generated via SGEES form LAPACK. 😌

🌸 Remarks 🌸

The algorithm is implemented using BLAS level 3 operations and OpenMP 4.0 DAG Scheduling.

💀 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_TRSYLV_L3. fr fr, don't skip this one. 😭

✨ Parameters ✨

• TRANSA — is CHARACTER(1). Specifies the form of the system of equations with respect to A:

== 'N':  op1(A) = A
== 'T':  op1(A) = A**T

• TRANSB — is CHARACTER(1). Specifies the form of the system of equations with respect to B:

== 'N':  op2(B) = B,
== 'T':  op2(B) = B**T

• 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.
• LDB — is INTEGER. The leading dimension of the array B. LDB >= 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. LDB >= 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.

👑 Author 👑

Martin Koehler, MPI Magdeburg

🎀 Date 🎀

October 2025

Definition at line 177 of file sla_trsylv_dag.f90.

---

💜 sla_trsylv_kernel_44nn 💜

subroutine sla_trsylv_kernel_44nn (real sgn, integer m, integer n, real,
    dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer
    ldb, real, dimension(ldx, *) x, integer ldx, real scale, integer info)

Solver for a 4x4 Sylvester equation (TRANSA = N, TRANSB = N) 🍓

🌈 Purpose 🌈

okay so SLA_TRSYLV_KERNEL_44NN solves a Sylvester equation that looks exactly like this:

A * X  + SGN * X * B = SCALE * Y                              (1)

where A is a M-by-M quasi upper triangular matrix and B is N-by-N quasi upper triangular matrices. The right hand side Y and the solution X are M-by-N matrices. Typically the matrices A and B are create by SGEES from LAPACK. The algorithm is implemented without using BLAS level 2 operations. Because of that, the order of M and N is at most 4. And, for fast execution, the function does not check the input arguments. 💀 (yes bestie, it just trusts you — don't make it regret that)

Parameters

• SGN — SGN is REAL, allowed values: +/-1. Specifies the sign between the two parts of the Sylvester equation.
• M — M is INTEGER. The order of the matrices A and C. 4 >= M >= 0.
• N — N is INTEGER. The order of the matrices B and D. 4 >= 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).
• 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).
• 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 < SCAL <= 1 holds true.
• INFO — INFO is INTEGER. On output:
  • == 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 134 of file sla_trsylv_kernel_44_nn.f90.

---

sla_trsylv_kernel_44nt 🦋

subroutine sla_trsylv_kernel_44nt (real sgn, integer m, integer n, real,
    dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer
    ldb, real, dimension(ldx, *) x, integer ldx, real scale, integer info)

Solver for a 4x4 Sylvester equation (TRANSA = N, TRANSB = T)

Purpose:

SLA_TRSYLV_KERNEL_44NT solves a Sylvester equation of the following form:

A * X  + SGN * X * B**T = SCALE * Y                              (1)

where A is a M-by-M quasi upper triangular matrix and B is N-by-N quasi upper triangular matrices. The right hand side Y and the solution X are M-by-N matrices. Typically the matrices A and B are create by SGEES from LAPACK. The algorithm is implemented using without BLAS level 2 operations. Because of that, the order of M and N is at most 4. And for fast execution, the function does not check the input arguments. ✨

Parameters

• SGN — SGN is REAL, allowed values: +/-1. Specifies the sign between the two parts of the Sylvester equation.
• M — M is INTEGER. The order of the matrices A and C. 4 >= M >= 0.
• N — N is INTEGER. The order of the matrices B and D. 4 >= 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).
• 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).
• 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 < SCAL <= 1 holds true.
• INFO — INFO is INTEGER. On output:
  • == 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 133 of file sla_trsylv_kernel_44_nt.f90.

---

sla_trsylv_kernel_44tn 🌈

subroutine sla_trsylv_kernel_44tn (real sgn, integer m, integer n, real,
    dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer
    ldb, real, dimension(ldx, *) x, integer ldx, real scale, integer info)

Solver for a 4x4 Sylvester equation (TRANSA = T, TRANSB = N)

Purpose:

SLA_TRSYLV_KERNEL_44TN solves a Sylvester equation of the following form:

A**T * X  + SGN * X * B = SCALE * Y                              (1)

where A is a M-by-M quasi upper triangular matrix and B is N-by-N quasi upper triangular matrices. The right hand side Y and the solution X are M-by-N matrices. Typically the matrices A and B are create by SGEES from LAPACK. The algorithm is implemented without using BLAS level 2 operations. Because of that, the order of M and N is at most 4. And for fast execution, the function does not check the input arguments. 💅

Parameters

• SGN — SGN is REAL, allowed values: +/-1. Specifies the sign between the two parts of the Sylvester equation.
• M — M is INTEGER. The order of the matrices A and C. 4 >= M >= 0.
• N — N is INTEGER. The order of the matrices B and D. 4 >= 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).
• 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).
• 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 < SCAL <= 1 holds true.
• INFO — INFO is INTEGER. On output:
  • == 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 133 of file sla_trsylv_kernel_44_tn.f90.

---

sla_trsylv_kernel_44tt 👑

subroutine sla_trsylv_kernel_44tt (real sgn, integer m, integer n, real,
    dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer
    ldb, real, dimension(ldx, *) x, integer ldx, real scale, integer info)

Solver for a 4x4 Sylvester equation (TRANSA = T, TRANSB = T)

Purpose:

• SCALE — REAL, the lil' safety-scaling factor that keeps your result from overflowing into the void ✨ if INFO == 0 then SCALE is 1.0; otherwise, when one of the inner systems couldn't be solved correctly, 0 < SCAL <= 1 holds true (yes bestie, that's exactly how it's written in the docs — SCAL, no E — we report, we don't edit 💅)
• INFO — INTEGER, the receipts. On output:
  • == 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 133 of file sla_trsylv_kernel_44_tt.f90.

---

👑 sla_trsylv_l2

subroutine sla_trsylv_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(ldx, *) x, integer ldx, real scale, real,
    dimension(*) work, integer info)

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

✨ Purpose:

!> !> SLA_TRSYLV_L2 solves a Sylvester equation of the following forms !> !> op1(A) * X + X * op2(B) = SCALE * Y (1) !> !> or !> !> op1(A) * X - X * op2(B) = SCALE * Y (2) !> !> where A is a M-by-M quasi upper triangular matrix, B is a N-by-N quasi upper triangular !> matrix. The right hand side Y and the solution X are M-by-N matrices. Typically the matrices !> A and B are generated via SGEES form LAPACK. !>

✨ Parameters

• TRANSA — CHARACTER(1). Tells us what form the system of equations takes with respect to A:
  • == 'N': op1(A) = A
  • == 'T': op1(A) = A**T
• TRANSB — CHARACTER(1). Same energy, but for B:
  • == 'N': op2(B) = B,
  • == 'T': op2(B) = B**T
• SGN — REAL, allowed values: +/-1. Specifies the sign between the two parts of the Sylvester equation.
  • = 1 : Solve Equation (1)
  • == -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. ⭐
• 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. ⭐
• LDB — INTEGER. The leading dimension of the array B. LDB >= 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. 💖
• LDX — INTEGER. The leading dimension of the array X. LDB >= max(1,M) (yup, that's what it says — LDB, on the LDX entry. Copying it down faithfully, no cap 👀)
• SCALE — REAL. 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 level-3 by level-2 and level-1 calls,
• 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 180 of file sla_trsylv_l2.f90.

---

👑 sla_trsylv_l2_local_copy

subroutine sla_trsylv_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(ldx, *) x, integer ldx, real scale, real,
    dimension(*) work, integer info)

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

✨ Purpose:

okay this is where the next part swoops in, bestie — the purpose section for this one is right there and we're picking it up next 💜

okay so SLA_TRSYLV_L2_LOCAL_COPY solves a Sylvester equation that looks like this:

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

or

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

where A is a M-by-M quasi upper triangular matrix and B is a N-by-N quasi upper triangular matrix. The right hand side Y and the solution X are M-by-N matrices. Typically the matrices A and B are generated via SGEES form LAPACK. no cap, that's the whole vibe 💖

🌸 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.
• LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= 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).
• 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 literally have to ask first, it's giving "let me check my calendar" energy 🙈)
• 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 (we're going smaller, and honestly? it's serving)
• 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)

👑 Author

Martin Koehler, MPI Magdeburg

📅 Date

October 2025

Definition at line 181 of file sla_trsylv_l2_opt_local_copy.f90.

---

🦋 subroutine sla_trsylv_l2_local_copy_128

subroutine sla_trsylv_l2_local_copy_128 ( transa,  transb,  sgn,  m,  n,  a,
    lda,  b,  ldb,  x,  ldx,  scale,  work,  info)

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

💫 Purpose 💫

SLA_TRSYLV_L2_LOCAL_COPY_128 solves a Sylvester equation of the following forms

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

or

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

where A is a M-by-M quasi upper triangular matrix and B is a N-by-N quasi upper triangular matrix. The right hand side Y and the solution X are M-by-N matrices. Typically the matrices A and B are generated via SGEES form LAPACK. 💖

🌸 Parameters 🌸

• TRANSA — TRANSA is CHARACTER(1). Specifies the form of the system of equations with respect to A:
  • == 'N': op1(A) = A
  • == 'T': op1(A) = A**T
• TRANSB — TRANSB is CHARACTER(1). Specifies the form of the system of equations with respect to B:
  • == 'N': op2(B) = B,
  • == 'T': op2(B) = B**T
• 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. 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.
• LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= 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).
• 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.
• Align local copies and fix problem size to <= 64

👑 Author

Martin Koehler, MPI Magdeburg

📅 Date

October 2025

Definition at line 183 of file sla_trsylv_l2_opt_local_copy_128.f90.

---

🧚 subroutine sla_trsylv_l2_local_copy_32

subroutine sla_trsylv_l2_local_copy_32 ( transa,  transb,  sgn,  m,  n,  a,
    lda,  b,  ldb,  x,  ldx,  scale,  work,  info)

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

💫 Purpose 💫

okay so SLA_TRSYLV_L2_LOCAL_COPY_32 solves a Sylvester equation, and it comes in these two flavors:

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

or

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

where A is a M-by-M quasi upper triangular matrix and B is a N-by-N quasi upper triangular matrix. The right hand side Y and the solution X are M-by-N matrices. Typically bestie, the matrices A and B get generated via SGEES form LAPACK. 🦋

---

🎀 Parameters 🎀

• TRANSA — TRANSA is CHARACTER(1). It says which form the system of equations takes with respect to A: if it's 'N' then op1(A) = A, and if it's 'T' then op1(A) = A**T. 💅
• TRANSB — TRANSB is CHARACTER(1). Same energy but for B: 'N' means op2(B) = B, and 'T' means op2(B) = B**T. ✨
• SGN — SGN is REAL, allowed values: +/-1. It picks 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. 32 >= M >= 0.
• N — N is INTEGER. The order of the matrices B and D. 32 >= 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.
• LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= 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 😭 the doc literally says LDB in that LDX description. i am NOT fixing it, that's what it says, periodt.

• 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. 💫
• WORK — WORK is REAL array, dimension LWORK. Workspace for the algorithm. Lowkey you have to query the workspace before running the computation — you do the query by calling the subroutine with INFO == -1 on input, and then the required workspace gets 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.
• Align local copies and fix problem size to <= 64

---

👑 Author 👑

Martin Koehler, MPI Magdeburg

🌸 Date 🌸

October 2025

Definition at line 183 of file sla_trsylv_l2_opt_local_copy_32.f90.

---

✨ subroutine sla_trsylv_l2_local_copy_64 ( transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info) ✨

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

🎀 Purpose: 🎀

SLA_TRSYLV_L2_LOCAL_COPY_64 solves a Sylvester equation, same two forms as before bestie:

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

or

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

where A is a M-by-M quasi upper triangular matrix and B is a N-by-N quasi upper triangular matrix. The right hand side Y and the solution X are M-by-N matrices. Typically the matrices A and B get generated via SGEES form LAPACK. 🦋

💜 Parameters 💜

• TRANSA — TRANSA is CHARACTER(1). Which form the system of equations takes with respect to A: 'N' means op1(A) = A, 'T' means op1(A) = A**T. ✨
• TRANSB — TRANSB is CHARACTER(1). Which form with respect to B: 'N' means op2(B) = B, and 'T' means op2(B) = B**T. ✨
• SGN — SGN is REAL, allowed values: +/-1. Picks 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. 64 >= M >= 0.
• N — N is INTEGER. The order of the matrices B and D. 64 >= 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.
• LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= 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).

yeah bestie it says LDB again. copy-paste is a lifestyle apparently 💀

• SCALE — SCALE is REAL. 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. 💫
• WORK — WORK is REAL array, dimension LWORK. Workspace for the algorithm. You gotta query the workspace before running the computation — call the subroutine with INFO == -1 on input and the required workspace comes back 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.
• Align local copies and fix problem size to <= 64

---

👑 Author 👑

Martin Koehler, MPI Magdeburg

🌸 Date 🌸

October 2025

Definition at line 182 of file sla_trsylv_l2_opt_local_copy_64.f90.

---

🦋 subroutine sla_trsylv_l2_local_copy_96 ( transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info) 🦋

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

🎀 Purpose: 🎀

so basically SLA_TRSYLV_L2_LOCAL_COPY_96 solves a Sylvester equation that looks like one of these two glow-ups:

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

or

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

where A is a M-by-M quasi upper triangular matrix and B is a N-by-N quasi upper triangular matrix 🦋. The right hand side Y and the solution X are M-by-N matrices. Typically the matrices A and B are generated via SGEES form LAPACK. iykyk ✨

---

🎀 Parameters 🎀

• TRANSA — TRANSA is CHARACTER(1). Specifies the form of the system of equations with respect to A:
  • == 'N': op1(A) = A
  • == 'T': op1(A) = A**T
• TRANSB — TRANSB is CHARACTER(1). Specifies the form of the system of equations with respect to B:
  • == 'N': op2(B) = B,
  • == 'T': op2(B) = B**T
• 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. 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.
• LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= 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).
• 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.
• Align local copies and fix problem size to <= 64

🌸 Author 🌸

Martin Koehler, MPI Magdeburg

⭐ Date ⭐

October 2025

👑 Definition 👑

At line 181 of file sla_trsylv_l2_opt_local_copy_96.f90.

---

💖 sla_trsylv_l2_reorder 💖

   subroutine sla_trsylv_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(ldx, *) x, integer ldx, real scale, real,
       dimension(*) work, integer info)

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

🦋 Purpose:

okay so SLA_TRSYLV_L2_REORDER solves a Sylvester equation of the following forms:

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

or

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

where A is a M-by-M quasi upper triangular matrix, B is a N-by-N quasi upper triangular matrix. The right hand side Y and the solution X are M-by-N matrices. Typically the matrices A and B are generated via SGEES form LAPACK. fr fr 💜

Attention 🙈 The algorithm is implemented using BLAS level 2 operations without further optimizations. For a faster implementation see SLA_TRSYLV_L2.

🎀 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.
• LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= 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).
• 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,

🌸 Author 🌸

Martin Koehler, MPI Magdeburg

⭐ Date ⭐

October 2025

👑 Definition 👑

At line 182 of file sla_trsylv_l2_opt_reorder.f90.

---

✨ sla_trsylv_l2_unopt ✨

   subroutine sla_trsylv_l2_unopt (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(ldx, *) x, integer ldx, real scale, real,
       dimension(*) work, integer info)

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

🦋 Purpose:

Parameters (continued from the l2_unopt version, bestie, we're still in it)

• TRANSA — it's a CHARACTER(1). Tells us the form of the system of equations with respect to A and C, and the valid vibes are:
  • 'N': op1(A) = A, op1(C) = C (no transpose for A and C, they're serving as-is 💅)
  • 'T': op1(A) = A**T, op1(C) = C **T (transpose A and C, flip 'em)
• TRANSB — also a 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 — a REAL, and the only allowed values are +/-1, no in-between, pick a side bestie. It specifies the sign between the two parts of the Sylvester equation:
  • = 1 : solve Equation (1)
  • == -1: solve Equation (2)
• M — INTEGER. The order of the matrices A and C. Constraint: M >= 0.
• N — INTEGER. The order of the matrices B and D. Constraint: N >= 0.
• A — REAL array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular. That's not a suggestion, that's the law 🫶
• 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.
• LDB — INTEGER. The leading dimension of the array B. LDB >= 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. So yeah she's doing double duty, input then output, main character behaviour honestly 👑
• LDX — INTEGER. The leading dimension of the array X. LDB >= max(1,M). (yes it literally says LDB in there, we are reporting what the page says, no notes 📝)
• 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. She's holding your hand so nothing explodes 💖
• WORK — REAL array, dimension LWORK. Workspace for the algorithm. And listen: 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. Do the query first, I'm begging 😭
• INFO — INTEGER. Buckle up:
  • 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.

Like. Nothing. Zero. Not a single one. It's giving "we're just happy to be here" 🎀

---

Author

Martin Koehler, MPI Magdeburg

Date

October 2025

Definition at line 180 of file sla_trsylv_l2_unopt.f90.

---

subroutine sla_trsylv_l3

subroutine sla_trsylv_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(ldx, *) x, integer ldx, real scale, real,
    dimension(*) work, integer info)

Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation. Okay so this is the same energy as the L2 one but with BLAS level 3, and we love an upgrade 🦋

Purpose:

!> SLA_TRSYLV_L3 solves a Sylvester equation of the following forms
!>
!>    op1(A) * X  +  X * op2(B) = SCALE * Y                              (1)
!>
!> or
!>
!>    op1(A) * X  -  X * op2(B) = SCALE * Y                              (2)
!>
!> where A is a M-by-M quasi upper triangular matrix, B is a N-by-N quasi upper triangular
!> matrix. The right hand side Y and the solution X are M-by-N matrices.  Typically the matrices
!> A and B are generated via SGEES form LAPACK.
!>

Attention

The algorithm is implemented using BLAS level 3 operations. That's the whole glow-up, no further notes 💅

Parameters

• TRANSA — CHARACTER(1). Specifies the form of the system of equations with respect to A:
  • 'N': op1(A) = A
  • 'T': op1(A) = A**T
• TRANSB — CHARACTER(1). Specifies the form of the system of equations with respect to B:
  • 'N': op2(B) = B,
  • 'T': op2(B) = B**T
• SGN — REAL, allowed values: +/-1. Specifies the sign between the two parts of the Sylvester equation:
  • = 1 : solve Equation (1)
  • == -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.
• 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.
• LDB — INTEGER. The leading dimension of the array B. LDB >= 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.
• LDX — INTEGER. The leading dimension of the array X. LDB >= max(1,M). (still says LDB, still quoting it exactly, iykyk 🙈)
• 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.

---

Author

Martin Koehler, MPI Magdeburg

Date

October 2025

Definition at line 172 of file sla_trsylv_l3_opt.f90.

---

subroutine sla_trsylv_l3_2s

subroutine sla_trsylv_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(ldx, *) x, integer ldx, real scale, real,
    dimension(*) work, integer info)

Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation. Yeah, same description as the one above — the "2s" is the plot twist and we're about to find out what it means 👀

Purpose:

SLA_TRSYLV_L3_2S is out here solving a Sylvester equation, and the forms are literally these two, bestie:

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

or

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

where A is a M-by-M quasi upper triangular matrix, B is a N-by-N quasi upper triangular matrix, and the right hand side Y plus the solution X are both M-by-N matrices. Typically A and B are generated via SGEES form LAPACK. 🦋

✨ Remarks ✨

The algorithm is implemented using BLAS level 3 operations and a two level inner solver — that's a massive parallel DAG scheduled solver from SLA_TRSYLV_DAG plus optimized kernel solvers. Two layers of solving. Main character energy fr fr. 👑

💀 Attention 💀

Due to the parallel nature of the inner solvers the scaling is turned off and SCALE is set to ONE. So do not come looking for a scale factor here, she doesn't do that. 🫶

🌸 Parameters 🌸

• TRANSA — TRANSA is CHARACTER(1). It specifies the form of the system of equations with respect to A:
  • == 'N': op1(A) = A
  • == 'T': op1(A) = A**T
• TRANSB — TRANSB is CHARACTER(1). It specifies the form of the system of equations with respect to B:
  • == 'N': op2(B) = B,
  • == 'T': op2(B) = B**T
• 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.
• LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= 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).
• 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 174 of file sla_trsylv_l3_2stage.f90.

---

okay NEXT bestie, new function, same energy 💖

subroutine sla_trsylv_l3_unopt (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(ldx, *) x, integer ldx, real scale, real,
    dimension(*) work, integer info)

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

Purpose:

SLA_TRSYLV_L3_UNOPT solves a Sylvester equation of the following forms

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

or

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

where A is a M-by-M quasi upper triangular matrix, B is a N-by-N quasi upper triangular matrix. The right hand side Y and the solution X are M-by-N matrices. Typically the matrices A and B are generated via SGEES form LAPACK. ✨

💅 Attention 💅

The algorithm is implemented using BLAS level 3 operations.

🌈 Parameters 🌈

• TRANSA — TRANSA is CHARACTER(1). It specifies the form of the system of equations with respect to A:
  • == 'N': op1(A) = A
  • == 'T': op1(A) = A**T
• TRANSB — TRANSB is CHARACTER(1). It specifies the form of the system of equations with respect to B:
  • == 'N': op2(B) = B,
  • == 'T': op2(B) = B**T
• 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.
• LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= 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). (yes bestie, it says LDB here, that's literally what the page says, 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.
• 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 172 of file sla_trsylv_l3_unopt.f90.

---

and now for something recursive, periodt 💜

recursive subroutine sla_trsylv_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(ldx, *) x, integer ldx, real scale, real,
    dimension(*) work, integer info)

💖 Level-3 Recursive Blocked Solver for the Sylvester equation 💖

Purpose:

she recurses, she blocks, she solves — details incoming, stay tuned bestie 🌸

okay bestie, so SLA_TRSYLV_RECURSIVE is literally here to solve a Sylvester equation, and it comes in two iconic forms: 💅

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

or

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

where A is a M-by-M quasi upper triangular matrix, B is a N-by-N quasi upper triangular matrix. The right hand side Y and the solution X are M-by-N matrices. Typically the matrices A and B are generated via SGEES form LAPACK. 💅✨

✨ ATTENTION ✨

The algorithm is implemented using BLAS level 3 operations and recursive blocking. So it's fast, bestie. No cap. 🔥👑

✨ PARAMETERS ✨

• TRANSA — TRANSA is CHARACTER(1). Specifies the form of the system of equations with respect to A:
  • == 'N': op1(A) = A
  • == 'T': op1(A) = A**T
• TRANSB — TRANSB is CHARACTER(1). Specifies the form of the system of equations with respect to B:
  • == 'N': op2(B) = B,
  • == 'T': op2(B) = B**T
• 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.
• LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= 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). (yes bestie, it says LDB 💀)
• 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 1. 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 160 of file sla_trsylv_recursive.f90.

✨ AUTHOR ✨

Generated automatically by Doxygen for MEPACK from the source code.

---

And that's the whole vibe, bestie! This page is from Version 1.1.2, dated Thu Oct 9 2025 00:00:00, for sgltrsylv(3). Stay iconic! 💖✨🦋🌈