dbltrsylv(3)

📦

version 1.1.2 · 3 versions

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

✨ dbltrsylv ✨

okay bestie, welcome to dbltrsylv(3) — it's a MEPACK page and it is THICC. Like, 17-parts long. This is part 1. Buckle up, hydrate, and let's go 💖

🌸 NAME 🌸

dbltrsylv — Double Precision

Double Precision routines for triangular standard Sylvester.

🦋 SYNOPSIS 🦋

🎀 Functions 🎀

subroutine dla_trsylv2_dag (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info) — the Level-3 Bartels-Stewart Algorithm for the discrete time Sylvester equation, except it's running with DAG based parallelization. parallel queen behavior fr fr.
subroutine dla_trsylv_dag (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info) — Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation, also with DAG parallelization. the generalized one, not the discrete one, don't mix them up 💅
subroutine dla_trsylv2_kernel_44nn (sgn, m, n, a, lda, b, ldb, x, ldx, scale, info) — Solver for a 4x4 discrete time Sylvester equation (TRANSA = N, TRANSB = N).
subroutine dla_trsylv2_kernel_44nt (sgn, m, n, a, lda, b, ldb, x, ldx, scale, info) — Solver for a 4x4 discrete time Sylvester equation (TRANSA = N, TRANSB = T).
subroutine dla_trsylv2_kernel_44tn (sgn, m, n, a, lda, b, ldb, x, ldx, scale, info) — Solver for a 4x4 discrete time Sylvester equation (TRANSA = T, TRANSB = N).
subroutine dla_trsylv2_kernel_44tt (sgn, m, n, a, lda, b, ldb, x, ldx, scale, info) — Solver for a 4x4 discrete time Sylvester equation (TRANSA = T, TRANSB = T).
subroutine dla_trsylv_kernel_44nn (sgn, m, n, a, lda, b, ldb, x, ldx, scale, info) — Solver for a 4x4 Sylvester equation (TRANSA = N, TRANSB = N).
subroutine dla_trsylv_kernel_44nt (sgn, m, n, a, lda, b, ldb, x, ldx, scale, info) — Solver for a 4x4 Sylvester equation (TRANSA = N, TRANSB = T).
subroutine dla_trsylv_kernel_44tn (sgn, m, n, a, lda, b, ldb, x, ldx, scale, info) — Solver for a 4x4 Sylvester equation (TRANSA = T, TRANSB = N).
subroutine dla_trsylv_kernel_44tt (sgn, m, n, a, lda, b, ldb, x, ldx, scale, info) — Solver for a 4x4 Sylvester equation (TRANSA = T, TRANSB = T).
subroutine dla_trsylv2_l2 (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). glow-up included ✨
subroutine dla_trsylv2_l2_local_copy (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).
subroutine dla_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).
subroutine dla_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).
subroutine dla_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).
subroutine dla_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).
subroutine dla_trsylv2_l2_reorder (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).
subroutine dla_trsylv2_l2_unopt (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 (unoptimized). no glow-up, just vibes 😭
subroutine dla_trsylv_l2 (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info) — Level-2 Bartels-Stewart Algorithm for the Sylvester equation (optimized).
subroutine dla_trsylv_l2_local_copy (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info) — Level-2 Bartels-Stewart Algorithm for the Sylvester equation (optimized).
subroutine dla_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).
subroutine dla_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).
subroutine dla_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).
subroutine dla_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).
subroutine dla_trsylv_l2_reorder (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info) — Level-2 Bartels-Stewart Algorithm for the Sylvester equation (optimized).
subroutine dla_trsylv_l2_unopt (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info) — Level-2 Bartels-Stewart Algorithm for the Sylvester equation (unoptimized).
subroutine dla_trsylv2_l3_2s (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info) — Two level blocked Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized).
subroutine dla_trsylv2_l3 (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info) — Level-3 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized).
subroutine dla_trsylv2_l3_unopt (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info) — Level-3 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Not Optimized).
subroutine dla_trsylv_l3_2s (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info) — Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation.
subroutine dla_trsylv_l3 (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info) — Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation.
subroutine dla_trsylv_l3_unopt (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info) — Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation (unoptimized).
recursive subroutine dla_trsylv2_recursive (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info) — Level-3 Recursive Blocked Algorithm for the discrete time Sylvester equation.
recursive subroutine dla_trsylv_recursive (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info) — Level-3 Recursive Blocked Solver for the Sylvester equation.

💅 Detailed Description 💅

Double Precision routines for triangular standard Sylvester.

This section contains the solvers for the standard Sylvester equation with (quasi) triangular coefficient matrices. The coefficient matrices are normally generated with the help of the Schur decomposition from LAPACK. The routines use double precision arithmetic for computation and data.

so real talk: the coefficient matrices don't just appear, bestie — they usually come out of the Schur decomposition over in LAPACK. and everything in here is double precision, both the math and the data. no float shenanigans 💅

👑 Function Documentation 👑

🫶 subroutine dla_trsylv2_dag 🫶

subroutine dla_trsylv2_dag (character, dimension(1) transa, character,
    dimension(1) transb, double precision sgn, integer m, integer n, double
    precision, dimension(lda, m) a, integer lda, double precision,
    dimension(ldb, n) b, integer ldb, double precision, dimension(ldx, n) x,
    integer ldx, double precision scale, double precision, dimension(*)
    work, integer info)

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

🌟 Purpose 🌟

DLA_TRSYLV2_DAG solves a discrete time Sylvester equation of the following forms:

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

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 DGEES from LAPACK.

so like, <pink>A</pink> is M-by-M quasi upper triangular, <pink>B</pink> is N-by-N upper quasi triangular, and <pink>X</pink> + <pink>Y</pink> are the M-by-N girlies doing all the work. and yeah — DGEES from LAPACK is usually the one who made A and B in the first place, iykyk 🦋

🌸 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 the DLA_TRSYLV2_L3 routine.

not me being dramatic but this one MATTERS, bestie 😭 because the algorithm is parallel, the scaling is not applied to the right hand. if your problem is ill-posed and scaling shows up, you gotta run it again with a solver that has complete scaling support, like DLA_TRSYLV2_L3. no cap.

🎀 Parameters 🎀

TRANSA

TRANSA is CHARACTER(1)

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

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

TRANSB

TRANSB is CHARACTER(1)

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

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

SGN

SGN is DOUBLE PRECISION, 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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

...and that's exactly where this part cuts off, bestie — X's whole moment is waiting for you in the next part 🫖✨

X — 💖 X is a DOUBLE PRECISION</> 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) — whichever one TRANSA</>, TRANSB</>, and SGN</> picked, bestie. Right hand side Y and the solution X are M-by-N matrices. 🦋
LDX — LDX</> is INTEGER. The leading dimension of the array X. LDX >= max(1,M). no cap.
SCALE — ✨ SCALE</> is DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0; otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true. If SCALE .NE. 1 the problem is no solved correctly in this case one have to use an other solver. 👑

🚨 fr fr: if SCALE isn't exactly 1, that solve did not go correctly and you gotta go find a different solver. not a vibe, but that's the deal.

WORK — 🛠️ WORK</> is a DOUBLE PRECISION array, dimension LWORK. Workspace for the algorithm. The workspace needs to queried before the running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO. 🫶
INFO — INFO</> is INTEGER.
- On input:
  - == -1 : Perform a workspace query
  - <> -1: normal operation
- On exit, workspace query:
  - < 0 : if INFO = -i, the i-th argument had an illegal value
  - >= 0: The value of INFO is the required number of elements in the workspace.
- On exit, normal operation:
  - == 0: successful exit
  - < 0: if INFO = -i, the i-th argument had an illegal value
  - > 0: The equation is not solved correctly. One of the arising inner system got singular. 💀

Author — Martin Koehler, MPI Magdeburg 🎀

Date — October 2025 🌸

Definition at line 181 of file dla_trsylv2_dag.f90. ⭐

subroutine `dla_trsylv2_kernel_44nn` ✨

subroutine dla_trsylv2_kernel_44nn (double precision sgn, integer m, integer
    n, double precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, *) x,
    integer ldx, double precision scale, integer info)

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

Purpose 🎯

DLA_TRSYLV2_KERNEL_44NN solves a discrete time Sylvester equation of the following form

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 DGEES 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 — SGN</> is DOUBLE PRECISION, allowed values: +/-1. 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 a DOUBLE PRECISION 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 a DOUBLE PRECISION 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 a DOUBLE PRECISION 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. LDB >= max(1,M). 🫶
SCALE — SCALE</> is DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0; 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 at line 135 of file dla_trsylv2_kernel_44_nn.f90. ⭐

subroutine `dla_trsylv2_kernel_44nt` ✨

subroutine dla_trsylv2_kernel_44nt (double precision sgn, integer m, integer
    n, double precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, *) x,
    integer ldx, double precision scale, integer info)

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

Purpose 🎯

DLA_TRSYLV2_KERNEL_44NT solves a discrete time Sylvester equation of the following form

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

Parameters

SGN — SGN</> is DOUBLE PRECISION, allowed values: +/-1. 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 a DOUBLE PRECISION 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 a DOUBLE PRECISION 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 a DOUBLE PRECISION 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. LDB >= max(1,M). 🫶
SCALE — SCALE</> is DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0; 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 at line 134 of file dla_trsylv2_kernel_44_nt.f90. ⭐

subroutine `dla_trsylv2_kernel_44tn` ✨

subroutine dla_trsylv2_kernel_44tn (double precision sgn, integer m, integer
    n, double precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, *) x,
    integer ldx, double precision scale, integer info)

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

Purpose 🎯

DLA_TRSYLV2_KERNEL_44TN solves a discrete time Sylvester equation of the following form

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

Parameters

SGN — SGN</> is DOUBLE PRECISION, allowed values: +/-1. 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 a DOUBLE PRECISION array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular. ⭐
LDA —
LDA — INTEGER, the leading dimension of the array A. Constraint: LDA >= max(1,M). 💅
B — DOUBLE PRECISION array, dimension (LDB,N). Heads up bestie: the matrix B must be (quasi-) upper triangular. Non-negotiable, that's the vibe check. 🌸
LDB — INTEGER, the leading dimension of the array B. Constraint: LDB >= max(1,N). ✨
X — DOUBLE PRECISION 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 glow-up, honestly. 💖
LDX — INTEGER, the leading dimension of the array X. Constraint as written: LDB >= max(1,M). (yeah it says LDB in the source, we are reporting it exactly as it is, no editing the queen's receipts 💀)
SCALE — DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0, otherwise if one of the inner systems could not be solved correctly, 0 < SCAL <= 1 holds true. 🫶
INFO — 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 dla_trsylv2_kernel_44_tn.f90.

🦋 `subroutine dla_trsylv2_kernel_44tt`

subroutine dla_trsylv2_kernel_44tt (double precision sgn, integer m, integer
    n, double precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, *) x,
    integer ldx, double precision scale, integer info)

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

Purpose 👑

!>
!> DLA_TRSYLV2_KERNEL_44TT solves a discrete time Sylvester equation of the following form
!>
!>    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
!> M-by-N matrices. Typically the matrices A and B are created by DGEES 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.
!>
!>

Okay so real talk: this one does not check the input arguments, because it's built for speed. You're on your own for validation, bestie. Speedrun energy. ⚡

Parameters 💅

SGN — DOUBLE PRECISION, allowed values: +/-1. Specifies the sign between the two parts of the Sylvester equation. ✨
M — INTEGER. The order of the matrix A. 4 >= M >= 0.
N — INTEGER. The order of the matrix B. 4 >= N >= 0.
A — DOUBLE PRECISION array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular. No cap. 🌸
LDA — INTEGER, the leading dimension of the array A. Constraint: LDA >= max(1,M).
B — DOUBLE PRECISION array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. 💜
LDB — INTEGER, the leading dimension of the array B. Constraint: LDB >= max(1,N).
X — DOUBLE PRECISION 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 — INTEGER, the leading dimension of the array X. Constraint as written: LDB >= max(1,M). (again, reporting it exactly, not our typo to fix 🙈)
SCALE — DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0, otherwise if one of the inner systems could not be solved correctly, 0 < SCAL <= 1 holds true. 🫶
INFO — 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 dla_trsylv2_kernel_44_tt.f90.

🌈 `subroutine dla_trsylv2_l2`

subroutine dla_trsylv2_l2 (character, dimension(1) transa, character,
    dimension(1) transb, double precision sgn, integer m, integer n, double
    precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, *) x,
    integer ldx, double precision scale, double precision, dimension(*)
    work, integer info)

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

Purpose 👑

!>
!> DLA_TRSYLV2_L2 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 DGEES 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 — DOUBLE PRECISION, 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 — DOUBLE PRECISION array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular.
LDA — INTEGER, the leading dimension of the array A. Constraint: LDA >= max(1,M).
B — DOUBLE PRECISION 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. 👑 Conditional queen behaviour, we love it.
LDB — INTEGER, the leading dimension of the array B. Constraint: LDB >= max(1,N).
X — DOUBLE PRECISION 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. Constraint as written: LDB >= max(1,M). (yes, still LDB, still verbatim, still iconic 💀)
SCALE — DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0, otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true. 🫶
WORK — DOUBLE PRECISION 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 DAXPY operation by Fortran intrinsic
- Replaced all BLAS calls except of DTRMV and DGER by Fortran intrinsic,

Author 🌸

Martin Koehler, MPI Magdeburg

Date 🎀

October 2025

Definition 💖

Definition at line 182 of file dla_trsylv2_l2.f90.

🦋 `subroutine dla_trsylv2_l2_local_copy` 🦋

subroutine dla_trsylv2_l2_local_copy (character, dimension(1) transa,
    character, dimension(1) transb, double precision sgn, integer m, integer
    n, double precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, *) x,
    integer ldx, double precision scale, double precision, dimension(*)
    work, integer info)

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

Purpose 💜

DLA_TRSYLV2_L2_LOCAL_COPY solves a generalized Sylvester equation of the following forms

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

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 DGEES from LAPACK. 🫶

Parameters 🌈

TRANSA — TRANSA is CHARACTER(1) bestie, and 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), and 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 DOUBLE PRECISION, 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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). 🙈 (yup, the docs literally say LDB right there — i'm just copying it exactly, don't come for me)
SCALE — SCALE is DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true. 💯
WORK — WORK is DOUBLE PRECISION 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 DAXPY operation by Fortran intrinsic
Replaced all BLAS calls except of DTRMV and DGER 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 💖

Definition at line 184 of file dla_trsylv2_l2_opt_local_copy.f90.

🧚 `subroutine dla_trsylv2_l2_local_copy_128` 🧚

subroutine dla_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 💜

DLA_TRSYLV2_L2_LOCAL_COPY_128 solves a generalized Sylvester equation of the following forms

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

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

Attention ⚠️

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

Parameters 🌈

TRANSA — TRANSA is CHARACTER(1), and 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), and 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 DOUBLE PRECISION, 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 — the docs said LDB, not LDX, i'm just the messenger)
SCALE — SCALE is DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true. 💯
WORK — WORK is DOUBLE PRECISION 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 DAXPY operation by Fortran intrinsic
Replaced all BLAS calls except of DTRMV 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 dla_trsylv2_l2_opt_local_copy_128.f90.

🦋 `dla_trsylv2_l2_local_copy_32`

subroutine dla_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). yes bestie it's giving optimized, lowercase-optimized, the whole vibe 💅

Purpose:

!>
!> DLA_TRSYLV2_L2_LOCAL_COPY_32 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 DGEES from LAPACK.
!>

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

Parameters

TRANSA — TRANSA is CHARACTER(1). It tells the routine what the system of equations looks like on the A and C side:
- == '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 deal but for the B and D girlies:
- == '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 DOUBLE PRECISION, 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. M >= 0.
N — N is INTEGER. The order of the matrices B and D. N >= 0.
A — A is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 what it says, keep it moving ✨)
SCALE — SCALE is DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — WORK is DOUBLE PRECISION 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 DAXPY operation by Fortran intrinsic
Replaced all BLAS calls except of DTRMV 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 dla_trsylv2_l2_opt_local_copy_32.f90.

🌸 `dla_trsylv2_l2_local_copy_64`

subroutine dla_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). same energy, bigger container, no notes 💖

Purpose:

!>
!> DLA_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 DGEES from LAPACK.
!>

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

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 DOUBLE PRECISION, 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — WORK is DOUBLE PRECISION 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:

ok bestie hold that thought ✋ the rest of this list is living in the next part, don't go anywhere 💜

okay so we were mid-list, continuing the glow-up receipts for the 64-variant, no repeats, let's go 💅

✨ Replaced level-3 by level-2 and level-1 calls,
🔁 Reorder the solution order to column first style,
🔥 Replaced DAXPY operation by Fortran intrinsic
💅 Replaced all BLAS calls except of DTRMV 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 dla_trsylv2_l2_opt_local_copy_64.f90.

🌸 subroutine `dla_trsylv2_l2_local_copy_96`

subroutine dla_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) — yes bestie, the sequel, and she's bigger 👑

🎯 Purpose

!> !> DLA_TRSYLV2_L2_LOCAL_COPY_96 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 DGEES from LAPACK. !>

!>
!> DLA_TRSYLV2_L2_LOCAL_COPY_96 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 DGEES from LAPACK.
!>

⚠️ Attention

The size of the Problem is limited by M,N <= 96 — read that again bestie, she taps out past 96, no cap 😭

🧾 Parameters

TRANSA — TRANSA is CHARACTER(1). Specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)
TRANSB — TRANSB is CHARACTER(1). Specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)
SGN — SGN is DOUBLE PRECISION, 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — WORK is DOUBLE PRECISION 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 DAXPY operation by Fortran intrinsic
💅 Replaced all BLAS calls except of DTRMV 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 dla_trsylv2_l2_opt_local_copy_96.f90.

🦋 subroutine `dla_trsylv2_l2_reorder`

subroutine dla_trsylv2_l2_reorder (character, dimension(1) transa,
    character, dimension(1) transb, double precision sgn, integer m, integer
    n, double precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, *) x,
    integer ldx, double precision scale, double precision, dimension(*)
    work, integer info)

Level-2 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized) — reorder edition, she reshuffled the whole vibe 🫶

🎯 Purpose

!>
!> DLA_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 DGEES from LAPACK.
!>

🧾 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 DOUBLE PRECISION, 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — WORK is DOUBLE PRECISION 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

okay and the receipts for this one roll in riiiight about now, bestie — 💖

Replaced level-3 by level-2 and level-1 calls, bestie. 💅
Reordered the solution order to column first style.
Replaced the DAXPY operation by a Fortran intrinsic. ✨
Replaced all BLAS calls except DTRMV and DGER by Fortran intrinsics.

✨ Author ✨

Martin Koehler, MPI Magdeburg

📅 Date 📅

October 2025

📍 Definition 📍

Definition at line 183 of file dla_trsylv2_l2_opt_reorder.f90.

   subroutine dla_trsylv2_l2_unopt (character, dimension(1) transa, character,
       dimension(1) transb, double precision sgn, integer m, integer n, double
       precision, dimension(lda, *) a, integer lda, double precision,
       dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, *) x,
       integer ldx, double precision scale, double precision, dimension(*)
       work, integer info)

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

🎯 Purpose:

DLA_TRSYLV2_L2_UNOPT solves a generalized Sylvester equation of the following forms:

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

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 DGEES from LAPACK.

🚨 Attention 🚨

The algorithm is implemented using BLAS level 2 operations without further optimizations. For a faster implementation see DLA_TRSYLV2_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 DOUBLE PRECISION, 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — WORK is DOUBLE PRECISION 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.

✨ Optimizations ✨

Nothing. 💀 Like, literally nothing. No cap.

✨ Author ✨

Martin Koehler, MPI Magdeburg

📅 Date 📅

October 2025

📍 Definition 📍

Definition at line 181 of file dla_trsylv2_l2_unopt.f90.

   subroutine dla_trsylv2_l3 (character, dimension(1) transa, character,
       dimension(1) transb, double precision sgn, integer m, integer n, double
       precision, dimension(lda, *) a, integer lda, double precision,
       dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, *) x,
       integer ldx, double precision scale, double precision, dimension(*)
       work, integer info)

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

🎯 Purpose:

okay so bestie, meet DLA_TRSYLV2_L3 💅 she solves a generalized Sylvester equation, and she does it in two possible shapes:

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

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 (yes, the original says "matrices" plural, we love a typo, we move on 🫶). The right hand side Y and the solution X are M-by-N matrices. And like, typically A and B get created by DGEES from LAPACK — that's the whole pipeline, iykyk ✨

🌸 Parameters 🌸

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 — 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 — is DOUBLE PRECISION, allowed values: +/-1. It 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDX,N). On input, the matrix X contains the right hand side Y. On output, the matrix X contains the solution of Equation (1) or (2) as selected by TRANSA, TRANSB, and SGN. Right hand side Y and the solution X are M-by-N matrices. She's doing double duty, main character energy honestly 💖
LDX — is INTEGER. The leading dimension of the array X. LDB >= max(1,M). (that's what the original says, bestie — I'm not gonna gaslight you and "fix" it, fr fr)
SCALE — is DOUBLE PRECISION. It's a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0, otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — is DOUBLE PRECISION 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. Query first, then run — no cap 🔥
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 dla_trsylv2_l3_opt.f90.

🦋 `dla_trsylv2_l3_2s` 🦋

subroutine dla_trsylv2_l3_2s (character, dimension(1) transa, character,
    dimension(1) transb, double precision sgn, integer m, integer n, double
    precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, *) x,
    integer ldx, double precision scale, double precision, dimension(*)
    work, integer info)

Two level blocked Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized) — she's the optimized one, she's booked and busy 💅

✨ Purpose ✨

okay so DLA_TRSYLV2_L3_2S solves a generalized Sylvester equation of the following forms

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

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 ("matrices" again, we're consistent 🌸). The right hand side Y and the solution X are M-by-N matrices. Typically the matrices A and B are created by DGEES from LAPACK. 👑

🌈 Remarks 🌈

The algorithm used level-3 BLAS operations and a DAG scheduled inner solver. We love a parallel pipeline, she's thriving 🧚

💀 Attention 💀

Due to the parallel nature of the inner solvers the scaling is turned off and SCALE is set to ONE. So don't come crying to her about your scaling, bestie — she told you upfront 🙈

🌸 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 DOUBLE PRECISION, 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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). (yep, LDB again, verbatim, don't @ me 🫶)
SCALE — is DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0, otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — is DOUBLE PRECISION array, dimension LWORK. Workspace for the algorithm. The workspace needs to be queried before running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO.
INFO — is INTEGER.
- On input:
  - == -1 : Perform a workspace query
  - <> -1: normal operation
- On exit, workspace query:
  - < 0 : if INFO = -i, the i-th argument had an illegal value
  - >= 0: The value of INFO is the required number of elements in the workspace.
- On exit, normal operation:
  - == 0: successful exit
  - < 0: if INFO = -i, the i-th argument had an illegal value
  - > 0: The equation is not solved correctly. One of the arising inner system got singular.

Author — Martin Koehler, MPI Magdeburg

Date — October 2025

Definition at line 174 of file dla_trsylv2_l3_2stage.f90.

🎀 `dla_trsylv2_l3_unopt` 🎀

subroutine dla_trsylv2_l3_unopt (character, dimension(1) transa, character,
    dimension(1) transb, double precision sgn, integer m, integer n, double
    precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, *) x,
    integer ldx, double precision scale, double precision, dimension(*)
    work, integer info)

Level-3 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Not Optimized) — she's the chill one, no rush, no stress, still gets the job done 🌸

✨ Purpose ✨

`DLA_TRSYLV2_L3_UNOPT` — the "one column at a time" girlie 💅

okay so bestie, DLA_TRSYLV2_L3_UNOPT solves a generalized Sylvester equation of one of these two shapes:

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

   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 are M-by-N matrices. Typically the matrices A and B are created by DGEES from LAPACK. ✨

Attention — this function iterates column first over the result, and from this fact it will be a bit slower than DLA_TRSYLV2_L3. it's giving "doing it the hard way on purpose", 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 DOUBLE PRECISION, 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — is DOUBLE PRECISION 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 175 of file dla_trsylv2_l3_unopt.f90.

`dla_trsylv2_recursive` — the recursive glow-up 🦋

recursive subroutine dla_trsylv2_recursive (character, dimension(1) transa,
    character, dimension(1) transb, double precision sgn, integer m, integer
    n, double precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, *) x,
    integer ldx, double precision scale, double precision, dimension(*)
    work, integer info)

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

Purpose:

okay so DLA_TRSYLV2_RECURSIVE solves a generalized Sylvester equation of the following forms

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

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

Remarks

The algorithm uses recursive blocking instead of the Bartels-Stewart approach. it's giving "we do things differently here" and honestly? iconic 👑

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 DOUBLE PRECISION, 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — is DOUBLE PRECISION 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 dla_trsylv2_recursive.f90.

`dla_trsylv_dag` — DAG parallelization, main character energy 👑

subroutine dla_trsylv_dag (character, dimension(1) transa, character,
    dimension(1) transb, double precision sgn, integer m, integer n, double
    precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, n) x,
    integer ldx, double precision scale, double precision, dimension(*)
    work, integer info)

Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation with DAG parallelization. ate and left no crumbs 💅

Purpose:

okay so DLA_TRSYLV_DAG solves a Sylvester equation, and there are two forms it can handle, bestie:

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

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 DGEES form LAPACK. ✨ slay ✨

💬 Remarks 💬

The algorithm is implemented using BLAS level 3 operations and OpenMP 4.0 DAG Scheduling. it's giving parallel, it's giving modern, periodt 💅

⚠️ 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 DLA_TRSYLV_L3. no cap, this one matters bestie 😭

🎀 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 DOUBLE PRECISION, 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — is DOUBLE PRECISION array, dimension LWORK. Workspace for the algorithm. The workspace needs to queried before the running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO. (yes you HAVE to query it first, don't skip this step bestie 🫶)
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 dla_trsylv_dag.f90.

🦋 `dla_trsylv_kernel_44nn` — Solver for a 4x4 Sylvester equation (TRANSA = N, TRANSB = N) 🦋

subroutine dla_trsylv_kernel_44nn (double precision sgn, integer m, integer n, double precision, dimension(lda, *) a, integer lda, double precision, dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, *) x, integer ldx, double precision scale, integer info)

🎯 Purpose 🎯

okay so DLA_TRSYLV_KERNEL_44NN solves a Sylvester equation of the following form:

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 M-by-N matrices. Typically the matrices A and B are create by DGEES from LAPACK. The algorithm is implemented without 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. so like… be careful, she trusts you completely 🙈

🎀 Parameters 🎀

SGN — is DOUBLE PRECISION, allowed values: +/-1. Specifies the sign between the two parts of the Sylvester equation.
M — is INTEGER. The order of the matrices A and C. 4 >= M >= 0.
N — is INTEGER. The order of the matrices B and D. 4 >= N >= 0.
A — is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0 otherwise if one of the inner systems could not be solved correctly, 0 < SCAL <= 1 holds true.
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 dla_trsylv_kernel_44_nn.f90.

🌈 `dla_trsylv_kernel_44nt` — Solver for a 4x4 Sylvester equation (TRANSA = N, TRANSB = T) 🌈

subroutine dla_trsylv_kernel_44nt (double precision sgn, integer m, integer n, double precision, dimension(lda, *) a, integer lda, double precision, dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, *) x, integer ldx, double precision scale, integer info)

🎯 Purpose 🎯

`dla_trsylv_kernel_44nt` 🦋

okay so bestie, DLA_TRSYLV_KERNEL_44NT solves a Sylvester equation that looks exactly like this:

   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 DGEES from LAPACK. ✨

The algorithm is implemented without BLAS level 2 operations, so the order of M and N is at most 4. And ngl, for fast execution the function does not check the input arguments. Like, at all. 🙈

Parameters

SGN — SGN is DOUBLE PRECISION, allowed values: +/-1. It 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 literally says LDB there, we do not make the rules 💀)
SCALE — SCALE is DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0 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 dla_trsylv_kernel_44_nt.f90.

`dla_trsylv_kernel_44tn` 💜

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

Purpose:

DLA_TRSYLV_KERNEL_44TN solves a Sylvester equation of the following form:

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

The algorithm is implemented without using BLAS level 2 operations, so the order of M and N is at most 4. And for fast execution the function does not check the input arguments. No cap. 🙈

Parameters

SGN — SGN is DOUBLE PRECISION, allowed values: +/-1. It 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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). (same energy as before, iconic 💀)
SCALE — SCALE is DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0 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 dla_trsylv_kernel_44_tn.f90.

`dla_trsylv_kernel_44tt` 🌈

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

Purpose:

DLA_TRSYLV_KERNEL_44TT solves a Sylvester equation of the following form:

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

The algorithm is implemented without using BLAS level 2 operations, so the order of M and N is at most 4. And for fast execution the function does not check the input arguments. Fr fr. 🙈

Parameters

SGN — SGN is DOUBLE PRECISION, allowed values: +/-1. It 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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). (third time's the charm bestie, still says LDB 💀)
SCALE — SCALE is DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0 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 dla_trsylv_kernel_44_tt.f90.

`dla_trsylv_l2` ⭐

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

Here's the signature, bestie, keep it safe:

subroutine dla_trsylv_l2 (character, dimension(1) transa, character,
    dimension(1) transb, double precision sgn, integer m, integer
    n, double precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, *) x,
    integer ldx, double precision scale, double precision, dimension(*)
    work, integer info)

Purpose:

okay so here's the full tea on what DLA_TRSYLV_L2 actually does 💖 bestie it solves a Sylvester equation in one of these two forms:

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

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

💖 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 DOUBLE PRECISION, 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDX,N). On input, the matrix X contains the right hand side Y. On output, the matrix X contains the solution of Equation (1) or (2) as selected by TRANSA, TRANSB, and SGN. Right hand side Y and the solution X are M-by-N matrices. it literally shape-shifts from Y to the answer, main character energy fr 🦋
LDX — LDX is INTEGER. The leading dimension of the array X. LDB >= max(1,M).
SCALE — SCALE is DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true. 🌸
WORK — WORK is DOUBLE PRECISION 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. bestie do NOT skip the query, the workspace will simply not be there for you 🫶
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 DAXPY operation by Fortran intrinsic
Replaced all BLAS calls except of DTRMV and DGER by Fortran intrinsic,

👑 Author

Martin Koehler, MPI Magdeburg

📅 Date

October 2025

Definition at line 180 of file dla_trsylv_l2.f90.

🦋 `dla_trsylv_l2_local_copy` 🦋

subroutine dla_trsylv_l2_local_copy (character, dimension(1) transa,
    character, dimension(1) transb, double precision sgn, integer m, integer
    n, double precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, *) x,
    integer ldx, double precision scale, double precision, dimension(*)
    work, integer info)

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

Purpose:

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

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

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

💜 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 DOUBLE PRECISION, 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true. 🌸
WORK — WORK is DOUBLE PRECISION 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 DAXPY operation by Fortran intrinsic
Replaced all BLAS calls except of DTRMV and DGER 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 dla_trsylv_l2_opt_local_copy.f90.

💀 `dla_trsylv_l2_local_copy_128` 💀

subroutine dla_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:

...and bestie, DLA_TRSYLV_L2_LOCAL_COPY_128 is out here solving Sylvester equations, and it does it in two flavours:

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

or

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

🎀 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 DOUBLE PRECISION, 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — WORK is DOUBLE PRECISION 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 DAXPY operation by Fortran intrinsic
Replaced all BLAS calls except of DTRMV 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

Definition at line 183 of file dla_trsylv_l2_opt_local_copy_128.f90.

🔥 `dla_trsylv_l2_local_copy_32` 🔥

subroutine dla_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:

TRANSA — TRANSA is CHARACTER(1). It picks what we're doing to A in the equation: == 'N' means op1(A) = A, and == 'T' means op1(A) = A**T. that's it, that's the whole vibe check 💅
TRANSB — TRANSB is CHARACTER(1). Same energy, but for B: == 'N' gives you op2(B) = B, and == 'T' gives you op2(B) = B**T. ✨
SGN — SGN is DOUBLE PRECISION, allowed values: +/-1. It's the sign between the two halves of the Sylvester equation, lowkey the whole plot twist: = 1 solves Equation (1), and == -1 solves Equation (2). 🎀
M — M is INTEGER. The order of the matrices A and C. And listen, this one is strict: 64 >= M >= 0. 👑
N — N is INTEGER. The order of the matrices B and D. Same rules bestie: 64 >= N >= 0. 🌸
A — A is a DOUBLE PRECISION array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular, no exceptions, no shortcuts. 🦋
LDA — LDA is INTEGER. The leading dimension of the array A. Requirement: LDA >= max(1,M). ⭐
B — B is a DOUBLE PRECISION array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular, fr fr. 🦋
LDB — LDB is INTEGER. The leading dimension of the array B. Requirement: LDB >= max(1,N). ⭐
X — X is a DOUBLE PRECISION array, dimension (LDX,N). On input, the matrix X contains the right hand side Y. On output, the matrix X contains the solution of Equation (1) or (2) as selected by TRANSA, TRANSB, and SGN. And just so we're clear: right hand side Y and the solution X are M-by-N matrices. it's giving glow-up, honestly 💖
LDX — LDX is INTEGER. The leading dimension of the array X. And the doc says LDB >= max(1,M) — yes bestie, we know, that's literally what's written, we do not edit the sacred texts. 🙈
SCALE — SCALE is DOUBLE PRECISION. It's the scaling factor that stops the result from overflowing, a true protective icon. If INFO == 0 then SCALE is 1.0D0, otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true. 🫶
WORK — WORK is a DOUBLE PRECISION array, dimension LWORK. Workspace for the algorithm. And no, you can't just guess a size: the workspace needs to be queried before running the computation, and 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.

so basically 0 is the slay, negatives are "you handed me nonsense", and positives are "the math said no". 💀

✨ Optimizations ✨

the receipts, bestie, here's what got upgraded:

Replaced level-3 by level-2 and level-1 calls,
Reorder the solution order to column first style,
Replaced DAXPY operation by Fortran intrinsic
Replaced all BLAS calls except of DTRMV 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 dla_trsylv_l2_opt_local_copy_64.f90.

🦋 subroutine dla_trsylv_l2_local_copy_96

subroutine dla_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:

DLA_TRSYLV_L2_LOCAL_COPY_96 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 DGEES 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 DOUBLE PRECISION, 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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. we don't touch it. 🙈)
SCALE — SCALE is DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true. 🫶
WORK — WORK is DOUBLE PRECISION 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 DAXPY operation by Fortran intrinsic
Replaced all BLAS calls except of DTRMV by Fortran intrinsic,
Use local copies of A, B, C, D, and X.
Align local copies and fix problem size to <= 64

okay so yes, the _96 routine says <= 64 in its optimization list. that's literally what's written, bestie. i'm the messenger, not the editor. 💀

👑 Author

Martin Koehler, MPI Magdeburg

📅 Date

October 2025

📍 Definition

At line 181 of file dla_trsylv_l2_opt_local_copy_96.f90.

🌈 subroutine dla_trsylv_l2_reorder

subroutine dla_trsylv_l2_reorder (character, dimension(1) transa, character,
    dimension(1) transb, double precision sgn, integer m, integer n, double
    precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, *) x,
    integer ldx, double precision scale, double precision, dimension(*)
    work, integer info)

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

💫 Purpose:

💅 `DLA_TRSYLV_L2_REORDER` 💅

okay so bestie, <pink>DLA_TRSYLV_L2_REORDER</pink> solves a Sylvester equation of one of these two forms:

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

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 DGEES form LAPACK. <lavender>(iykyk 🦋)

🙈 Attention 🙈

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

🌸 Parameters 🌸

TRANSA — TRANSA is CHARACTER(1). It says what form the system of equations takes with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)
TRANSB — TRANSB is CHARACTER(1). It says what form the system of equations takes with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)
SGN — SGN is DOUBLE PRECISION, 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. M >= 0.
N — N is INTEGER. The order of the matrices B and D. N >= 0.
A — A is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — WORK is DOUBLE PRECISION 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 DAXPY operation by Fortran intrinsic
Replaced all BLAS calls except of DTRMV and DGER by Fortran intrinsic,

🫶 Author 🫶

Martin Koehler, MPI Magdeburg

📅 Date 📅

October 2025

📍 Definition 📍

Definition at line 182 of file dla_trsylv_l2_opt_reorder.f90.

🧚 `dla_trsylv_l2_unopt` 🧚

subroutine dla_trsylv_l2_unopt (character, dimension(1) transa, character,
    dimension(1) transb, double precision sgn, integer m, integer n, double
    precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, *) x,
    integer ldx, double precision scale, double precision, dimension(*)
    work, integer info)

<cyan>Level-2 Bartels-Stewart Algorithm for the Sylvester equation (unoptimized)</cyan> (the raw, un-glow-up version bestie, no filter 😭)

💜 Purpose 💜

<hotpink>DLA_TRSYLV_L2_UNOPT</hotpink> solves a Sylvester equation of one of these two forms:

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

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

🙈 Attention 🙈

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

🌸 Parameters 🌸

TRANSA — TRANSA is CHARACTER(1). It says what form the system of equations takes with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)
TRANSB — TRANSB is CHARACTER(1). It says what form the system of equations takes with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)
SGN — SGN is DOUBLE PRECISION, 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. M >= 0.
N — N is INTEGER. The order of the matrices B and D. N >= 0.
A — A is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — WORK is DOUBLE PRECISION 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 ✨

Nothing.

literally nothing. it's giving "we did not touch a thing and we're proud" 💀

🫶 Author 🫶

Martin Koehler, MPI Magdeburg

📅 Date 📅

October 2025

📍 Definition 📍

Definition at line 180 of file dla_trsylv_l2_unopt.f90.

👑 `dla_trsylv_l3` 👑

subroutine dla_trsylv_l3 (character, dimension(1) transa, character,
    dimension(1) transb, double precision sgn, integer m, integer n, double
    precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, *) x,
    integer ldx, double precision scale, double precision, dimension(*)
    work, integer info)

<rainbow>Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation.</rainbow> ✨ the main character has entered the chat ✨

💜 Purpose 💜

okay so DLA_TRSYLV_L3 is out here solving Sylvester equations in two iconic flavours, and bestie they are literally the same equation with a vibe flip in the middle:

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

   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 DGEES form LAPACK. iykyk 💜

💅 Attention 💅

The algorithm is implemented using BLAS level 3 operations. level 3, we're not playing around ✨

🎀 Parameters 🎀

TRANSA — TRANSA is CHARACTER(1). It says which form of the system of equations we're using with respect to A:
- == 'N': op1(A) = A
- == 'T': op1(A) = A**T
TRANSB — TRANSB is CHARACTER(1). Same energy, but for B:
- == 'N': op2(B) = B,
- == 'T': op2(B) = B**T
SGN — SGN is DOUBLE PRECISION, 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 DOUBLE PRECISION array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular. quasi, bestie, don't forget the quasi 🦋
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — B is DOUBLE PRECISION 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 DOUBLE PRECISION 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 comes in as Y, it leaves as the answer, main character transformation 💅
LDX — LDX is INTEGER. The leading dimension of the array X. LDB >= max(1,M). (yep, it says LDB, we're just the messengers 🙈)
SCALE — SCALE is DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — WORK is DOUBLE PRECISION array, dimension LWORK. Workspace for the algorithm. The workspace needs to queried before the running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO. query first, run later, no cap 🍓
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 172 of file dla_trsylv_l3_opt.f90.

subroutine dla_trsylv_l3_2s

subroutine dla_trsylv_l3_2s (character, dimension(1) transa, character,
    dimension(1) transb, double precision sgn, integer m, integer n, double
    precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, *) x,
    integer ldx, double precision scale, double precision, dimension(*)
    work, integer info)

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

🧚 Purpose 🧚

DLA_TRSYLV_L3_2S solves a Sylvester equation of the following forms:

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

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

🌸 Remarks 🌸

The algorithm is implemented using BLAS level 3 operations and a two level inner solver consisting of a massive parallel DAG scheduled solver from DLA_TRSYLV_DAG and optimized kernel solvers. a whole little ecosystem, obsessed 💜

💀 Attention 💀

Due to the parallel nature of the inner solvers the scaling is turned off and SCALE is set to ONE.

bestie this is important: scaling is OFF here. SCALE is just ONE, periodt. Don't come for the scaling, it's not happening in this version.

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

✨ Author ✨

Martin Koehler, MPI Magdeburg

📅 Date 📅

October 2025

Definition at line 174 of file dla_trsylv_l3_2stage.f90.

subroutine dla_trsylv_l3_unopt

subroutine dla_trsylv_l3_unopt (character, dimension(1) transa, character,
    dimension(1) transb, double precision sgn, integer m, integer n, double
    precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, *) x,
    integer ldx, double precision scale, double precision, dimension(*)
    work, integer info)

Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation (unoptimized) the unoptimized one, she's doing her best 🍓

🦋 Purpose 🦋

okay so bestie, here's where it's at:

okay so DLA_TRSYLV_L3_UNOPT is the unoptimized level-3 solver and it's here to solve a Sylvester equation in two iconic forms, no cap:

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

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 DGEES form LAPACK.

✨ Attention ✨

The algorithm is implemented using BLAS level 3 operations. That's the whole vibe. 💅

🌸 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 DOUBLE PRECISION, 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — WORK is DOUBLE PRECISION 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 📍

Definition at line 172 of file dla_trsylv_l3_unopt.f90.

Next up, bestie, meet dla_trsylv_recursive — the level-3 recursive blocked solver for the Sylvester equation. Its signature is literally this:

recursive subroutine dla_trsylv_recursive (character, dimension(1) transa,
    character, dimension(1) transb, double precision sgn, integer m, integer
    n, double precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, *) x,
    integer ldx, double precision scale, double precision, dimension(*)
    work, integer info)

Level-3 Recursive Blocked Solver for the Sylvester equation.

🦋 Purpose 🦋

DLA_TRSYLV_RECURSIVE solves a Sylvester equation of the following forms:

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

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

💅 Attention 💅

The algorithm is implemented using BLAS level 3 operations and recursive blocking.

🌺 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 DOUBLE PRECISION, 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION. SCALE is a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0 otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — WORK is DOUBLE PRECISION 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 📍

Definition at line 160 of file dla_trsylv_recursive.f90.

✨ Author ✨

Generated automatically by Doxygen for MEPACK from the source code.

This whole thing is MEPACK version 1.1.2, dated Thu Oct 9 2025 00:00:00, and it's serving dbltrsylv(3) energy. That's the tea, bestie. Byeee 💖🦋✨

📦 where this page comes from

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

dbltrsylv(3)

✨ dbltrsylv ✨

🌸 NAME 🌸

🦋 SYNOPSIS 🦋

🎀 Functions 🎀

💅 Detailed Description 💅

👑 Function Documentation 👑

🫶 subroutine dla_trsylv2_dag 🫶

🌟 Purpose 🌟

🌸 Remarks 🌸

⚠️ Attention ⚠️

🎀 Parameters 🎀

subroutine dla_trsylv2_kernel_44nn ✨

subroutine dla_trsylv2_kernel_44nt ✨

subroutine dla_trsylv2_kernel_44tn ✨

Author 🎀

Date ⭐

Definition 📍

🦋 subroutine dla_trsylv2_kernel_44tt

Purpose 👑

Parameters 💅

Author 🎀

Date ⭐

Definition 📍

🌈 subroutine dla_trsylv2_l2

Purpose 👑

Parameters 💅

✨ Optimizations ✨

Author 🌸

Date 🎀

Definition 💖

🦋 subroutine dla_trsylv2_l2_local_copy 🦋

Purpose 💜

Parameters 🌈

Optimizations ⭐

Author 🌸

Date 🎀

Definition 💖

🧚 subroutine dla_trsylv2_l2_local_copy_128 🧚

Purpose 💜

Attention ⚠️

Parameters 🌈

Optimizations ⭐

🦋 dla_trsylv2_l2_local_copy_32

🌸 dla_trsylv2_l2_local_copy_64

👑 Author

🗓️ Date

📍 Definition

🌸 subroutine dla_trsylv2_l2_local_copy_96

🎯 Purpose

⚠️ Attention

🧾 Parameters

✨ Optimizations

👑 Author

🗓️ Date

📍 Definition

🦋 subroutine dla_trsylv2_l2_reorder

🎯 Purpose

🧾 Parameters

✨ Optimizations

✨ Author ✨

📅 Date 📅

📍 Definition 📍

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

🎯 Purpose:

🚨 Attention 🚨

✨ Parameters ✨

✨ Optimizations ✨

✨ Author ✨

📅 Date 📅

📍 Definition 📍

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

🎯 Purpose:

🌸 Parameters 🌸

🦋 dla_trsylv2_l3_2s 🦋

✨ Purpose ✨

🌈 Remarks 🌈

💀 Attention 💀

🌸 Parameters 🌸

🎀 dla_trsylv2_l3_unopt 🎀

subroutine `dla_trsylv2_kernel_44nn` ✨

subroutine `dla_trsylv2_kernel_44nt` ✨

subroutine `dla_trsylv2_kernel_44tn` ✨

🦋 `subroutine dla_trsylv2_kernel_44tt`

🌈 `subroutine dla_trsylv2_l2`

🦋 `subroutine dla_trsylv2_l2_local_copy` 🦋

🧚 `subroutine dla_trsylv2_l2_local_copy_128` 🧚

🦋 `dla_trsylv2_l2_local_copy_32`

🌸 `dla_trsylv2_l2_local_copy_64`

🌸 subroutine `dla_trsylv2_l2_local_copy_96`

🦋 subroutine `dla_trsylv2_l2_reorder`

🦋 `dla_trsylv2_l3_2s` 🦋

🎀 `dla_trsylv2_l3_unopt` 🎀

`DLA_TRSYLV2_L3_UNOPT` — the "one column at a time" girlie 💅

`dla_trsylv2_recursive` — the recursive glow-up 🦋

`dla_trsylv_dag` — DAG parallelization, main character energy 👑

🦋 `dla_trsylv_kernel_44nn` — Solver for a 4x4 Sylvester equation (TRANSA = N, TRANSB = N) 🦋

🌈 `dla_trsylv_kernel_44nt` — Solver for a 4x4 Sylvester equation (TRANSA = N, TRANSB = T) 🌈

`dla_trsylv_kernel_44nt` 🦋

`dla_trsylv_kernel_44tn` 💜

`dla_trsylv_kernel_44tt` 🌈

`dla_trsylv_l2` ⭐

🦋 `dla_trsylv_l2_local_copy` 🦋

💀 `dla_trsylv_l2_local_copy_128` 💀

🔥 `dla_trsylv_l2_local_copy_32` 🔥

💅 `DLA_TRSYLV_L2_REORDER` 💅

🧚 `dla_trsylv_l2_unopt` 🧚