dbltgsylv(3)

📦

version 1.1.2 · 3 versions

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

💖 `dbltgsylv` 💖 — double precision, triangular generalized Sylvester era ✨

okay so bestie, you just opened dbltgsylv(3) from the MEPACK girlies and I am so glad you're here 🫶 this page is the entire double precision lineup — every solver, kernel, and recursion situation for triangular generalized Sylvester equations. it's a lot. it's giving ✨algorithm catalog✨. we're gonna get through it together, no cap.

✨ NAME ✨

dbltgsylv — Double Precision 💖 — the Double Precision routines for triangular generalized Sylvester equations. that's the whole vibe, periodt.

✨ SYNOPSIS ✨

🌸 Functions 🌸

okay so every single one of these is a subroutine, bestie — except the last three, which are recursive subroutines and yes that distinction matters 💅 go off.

subroutine dla_tgcsylv_dag (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info) — the Level-3 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation, and it does the whole thing with DAG scheduling 🔥 she's organized, she's parallel, she ate.
subroutine dla_tgcsylv_dual_dag (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info) — the OpenMP4 - DAG Algorithm for the dual generalized coupled Sylvester equation ✨ threads said let's go.
subroutine dla_tgsylv_dag (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info) — Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation, with DAG based parallelization. iconic behaviour honestly.
subroutine dla_tgcsylv_dual_kernel_44nn (sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, info) — solver for a 4x4 dual generalized coupled Sylvester equation when TRANSA = N, TRANSB = N 💅
subroutine dla_tgcsylv_dual_kernel_44nt (sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, info) — same 4x4 dual generalized coupled Sylvester equation solver energy, but TRANSA = N, TRANSB = T 🦋
subroutine dla_tgcsylv_dual_kernel_44tn (sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, info) — 4x4 dual generalized coupled Sylvester equation solver, TRANSA = T, TRANSB = N 🌸
subroutine dla_tgcsylv_dual_kernel_44tt (sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, info) — and the last of the quartet: 4x4 dual generalized coupled Sylvester equation, TRANSA = T, TRANSB = T ⭐
subroutine dla_tgcsylv_dual_l2 (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info) — Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the dual generalized coupled Sylvester equation. level 2 but make it dual ✨
subroutine dla_tgcsylv_dual_l2_local_copy (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info) — Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the dual generalized coupled Sylvester equation, this one tagged (Local Copy) 💖
subroutine dla_tgcsylv_dual_l2_local_copy_128 (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info) — Level-2 Kagstroem-Westin/Bartels-Stewart for the dual generalized coupled Sylvester equation (Local Copy), the _128 flavour 🎀
subroutine dla_tgcsylv_dual_l2_local_copy_32 (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info) — Level-2 Kagstroem-Westin/Bartels-Stewart for the dual generalized coupled Sylvester equation (Local Copy), _32 edition 🍓
subroutine dla_tgcsylv_dual_l2_local_copy_64 (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info) — Level-2 Kagstroem-Westin/Bartels-Stewart for the dual generalized coupled Sylvester equation (Local Copy), _64 girlie 💜
subroutine dla_tgcsylv_dual_l2_local_copy_96 (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info) — Level-2 Kagstroem-Westin/Bartels-Stewart for the dual generalized coupled Sylvester equation (Local Copy), _96 and that's the whole local-copy set 🔥
subroutine dla_tgcsylv_l2 (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info) — Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation, and it's marked (Optimized) 💅
subroutine dla_tgcsylv_l2_local_copy (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info) — Level-2 Kagstroem-Westin/Bartels-Stewart for the generalized coupled Sylvester equation (Optimized) ✨
subroutine dla_tgcsylv_l2_local_copy_128 (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info) — Level-2 Kagstroem-Westin/Bartels-Stewart for the generalized coupled Sylvester equation (Optimized), _128 variant 🦋
subroutine dla_tgcsylv_l2_local_copy_32 (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info) — Level-2 Kagstroem-Westin/Bartels-Stewart for the generalized coupled Sylvester equation (Optimized), _32 variant 🌸
subroutine dla_tgcsylv_l2_local_copy_64 (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info) — Level-2 Kagstroem-Westin/Bartels-Stewart for the generalized coupled Sylvester equation (Optimized), _64 variant ⭐
subroutine dla_tgcsylv_l2_local_copy_96 (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info) — Level-2 Kagstroem-Westin/Bartels-Stewart for the generalized coupled Sylvester equation (Optimized), _96 variant 💜
subroutine dla_tgcsylv_l2_reorder (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info) — Level-2 Kagstroem-Westin/Bartels-Stewart for the generalized coupled Sylvester equation (Optimized), and this one reorders, hence the name, iykyk 💅
subroutine dla_tgcsylv_l2_unopt (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info) — Level-2 Kagstroem-Westin/Bartels-Stewart for the generalized coupled Sylvester equation, (unoptimized) — yes bestie, the unoptimized one exists too and we respect her 🫶
subroutine dla_tgsylv_gardiner_laub (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info) — the Level-2 Gardiner-Laub algorithm for the generalized Sylvester equation 🌈 a whole different algorithm family, we love range.
subroutine dla_tgsylv_l2 (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info) — Level-2 Bartels-Stewart Algorithm for the generalized Sylvester equation, (Optimized) ✨
subroutine dla_tgsylv_l2_colwise (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info) — Level-2 Bartels-Stewart for the generalized Sylvester equation, (unoptimized), going column-wise 💅
subroutine dla_tgsylv_l2_local_copy (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info) — Level-2 Bartels-Stewart for the generalized Sylvester equation, (Optimized) 💖
subroutine dla_tgsylv_l2_local_copy_128 (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info) — Level-2 Bartels-Stewart for the generalized Sylvester equation, (Optimized), _128 🎀
subroutine dla_tgsylv_l2_local_copy_32 (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info) — Level-2 Bartels-Stewart for the generalized Sylvester equation, (Optimized), _32 🍓
subroutine dla_tgsylv_l2_local_copy_64 (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info) — Level-2 Bartels-Stewart for the generalized Sylvester equation, (Optimized), _64 💜
subroutine dla_tgsylv_l2_local_copy_96 (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info) — Level-2 Bartels-Stewart for the generalized Sylvester equation, (Optimized), _96 🔥
subroutine dla_tgsylv_l2_reorder (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info) — Level-2 Bartels-Stewart for the generalized Sylvester equation, (Optimized), reordering edition ✨
subroutine dla_tgsylv_l2_rowwise (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info) — Level-2 Bartels-Stewart for the generalized Sylvester equation, (unoptimized), row-wise this time 🦋
subroutine dla_tgcsylv_dual_l3 (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info) — Level-3 Kagstroem-Westin/Bartels-Stewart Algorithm for the dual generalized coupled Sylvester equation 💅 level 3, dual, she's back and she's bigger.
subroutine dla_tgcsylv_dual_l3_2s (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info) — Level-3 Kagstroem-Westin/Bartels-Stewart Algorithm with sub blocking for the dual generalized coupled Sylvester equation 🌸 the 2s means sub blocking, and it's giving structure.
subroutine dla_tgcsylv_l3_2s (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info) — Level-3 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation, (optimized), sub-blocking style 🎀
subroutine dla_tgcsylv_l3 (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info) — Level-3 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation, (optimized) ✨
subroutine dla_tgcsylv_l3_unopt (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info) — Level-3 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation, (unoptimized) 💜 still valid, still serving.
subroutine dla_tgsylv_l3_2s (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info) — Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation, 2s flavour 🔥
subroutine dla_tgsylv_l3_colwise (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info) — Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation, going column-wise 🦋
subroutine dla_tgsylv_l3 (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info) — Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation ⭐ the classic.
subroutine dla_tgsylv_l3_rowwise (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info) — Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation, row-wise 🌈
recursive subroutine dla_tgcsylv_dual_recursive (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info) — the Recursive blocking Algorithm for the dual generalized coupled Sylvester equation 💅 note the recursive keyword bestie, she calls herself.
recursive subroutine dla_tgcsylv_recursive (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info) — the Recursive Algorithm for the generalized coupled Sylvester equation 🌸 recursion girlies rise.
recursive subroutine dla_tgsylv_recursive (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info) — the Level-3 Recursive Blocked Solver for the generalized Sylvester equation 👑 and that's the whole recursive trio.

and that's the SYNOPSIS roster for this part bestie 💖 there's so much more coming — don't go anywhere, we're not even close to done 👀

💖 Detailed Description 💖

okay so bestie, this whole section is about double precision routines for triangular generalized Sylvester equations. no cap. We're solving the standard Sylvester equation where the coefficient matrices are (quasi) triangular. Those matrices are normally generated with the help of the Schur decomposition from LAPACK. And every single routine here? Works in double precision arithmetic. periodt. 💅

✨ Function Documentation ✨

🦋 `dla_tgcsylv_dag` 🦋

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

This is the Level-3 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation, and it's got DAG scheduling on top. iconic behavior. 🔥

✨ Purpose ✨

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

So basically you've got these two equations, and A and C are M-by-M, B and D are N-by-N. A is quasi upper triangular, C is upper triangular, and B and D are upper triangular — but either B or D can be quasi upper triangular too. The right hand side Y and the solution X are M-by-N matrices (yeah, the original says Y and X, even though the equations use E, F, R, L — we're keeping it exactly as written, bestie 💀). The matrix pencils (A,C) and (B,D) (or (D,B)) usually come from DGGES in LAPACK. 💖

💬 Remarks 💬

The Algorithm uses OpenMP DAG scheduling. that's it, that's the remark. periodt. ✨

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

not me catching that the original literally says LDE >= max(1,M) for LDF 😭 we keep it exactly as written, bestie. no cap.

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.

Author: Martin Koehler, MPI Magdeburg

Date: October 2025

Definition at line 217 of file dla_tgcsylv_dag.f90.

🌈 `dla_tgcsylv_dual_dag` 🌈

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

OpenMP4 - DAG Algorithm for the dual generalized coupled Sylvester equation. it's giving duality, bestie. 🔥

✨ Purpose ✨

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

So yeah, this is the dual to the generalized coupled Sylvester equation (2). The whole Kronecker form thing shows why TRANSA and TRANSB are expressed in terms of equation (2), not (1). iykyk. 💅

🌸 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). ✨
SGN1 — DOUBLE PRECISION, allowed values: +/-1. Specifies the sign between in the first equation. 💅
SGN2 — DOUBLE PRECISION, allowed values: +/-1. Specifies the sign between in the second equation. 💅
M — INTEGER. The order of the matrices A and C. M >= 0. 👑
N — INTEGER. The order of the matrices B and D. N >= 0. 👑
A — DOUBLE PRECISION array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular. 🦋
LDA
LDA — INTEGER ✨ it's the leading dimension of the array A, so the rule is LDA >= max(1,M). No cap.
B — DOUBLE PRECISION array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. And heads up bestie: if the matrix D is already quasi-upper triangular, then B has to be plain upper triangular. it's giving teamwork 🤝
LDB — INTEGER 💅 leading dimension of the array B. Constraint: LDB >= max(1,N).
C — DOUBLE PRECISION array, dimension (LDC,M). The matrix C must be upper triangular, periodt.
LDC — INTEGER 🌸 leading dimension of the array C. Constraint: LDC >= max(1,M).
D — DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. And mirror-mirror: if the matrix B is already quasi-upper triangular, then D must be upper triangular. 🔁
LDD — INTEGER ⭐ leading dimension of the array D. Constraint: LDD >= max(1,N).
E — DOUBLE PRECISION array, dimension (LDE,N). On input, the matrix E holds the right hand side E. On output, the matrix E holds the solution R of Equation (1), as selected by TRANSA, TRANSB, and SGN1/SGN2. Both the right hand side E and the solution R are M-by-N matrices. 🦋 it literally transforms in place.
LDE — INTEGER 👑 leading dimension of the array E. Constraint: LDE >= max(1,M).
F — DOUBLE PRECISION array, dimension (LDF,N). On input, the matrix F holds the right hand side F. On output, the matrix F holds the solution L of Equation (1), as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and solution L are M-by-N matrices. ✨
LDF — INTEGER 🎀 leading dimension of the array F. Constraint, and I'm quoting this exactly: LDE >= max(1,M). yes it says LDE 🙈 we report, we don't rewrite.
SCALE — DOUBLE PRECISION. SCALE is a scaling factor that exists to prevent overflow in the result. If INFO == 0 then SCALE is 1.0D0; otherwise, if one of the inner systems could not be solved correctly, then 0 < SCALE <= 1 holds true. 🔥 protecting your numbers fr fr.
WORK — DOUBLE PRECISION array, dimension LWORK. Workspace for the algorithm. You have to query the workspace before running the computation — you do that by calling the subroutine with INFO == -1 on input, and then the required workspace comes back in INFO. 💅 don't skip the query bestie.
INFO — INTEGER 👑

On input:

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

On exit, workspace query:

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

On exit, normal operation:

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

Author 🌸

Martin Koehler, MPI Magdeburg

Date 📅

October 2025

Definition 🎀

Definition at line 237 of file dla_tgcsylv_dual_dag.f90.

🦋 subroutine `dla_tgcsylv_dual_kernel_44nn` 🦋

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

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

Purpose 💖

DLA_TGCSYLV_DUAL_KERNEL_44NN solves a generalized coupled Sylvester equation of the following form 👇

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

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

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

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

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

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

Regarding Z**T one obtains 🔁

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

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

Attention ⚠️

This routine only allows to solve equations that are at most of size 4 x 4. The operator arguments TRANSA and TRANSB are fixed to op1(H) = H and op2(H) = H. The primary use of this function is to end the recursion in dla_tgcsylv_dual_recursive. 🫶 so like, don't come here with your giant matrix, it is literally not that girl.

Parameters 🎀

SGN1 — DOUBLE PRECISION, allowed values: +/-1. Specifies the sign between in the first equation. 💖
SGN2 — DOUBLE PRECISION, allowed values: +/-1. Specifies the sign between in the second equation. 💖
M — INTEGER. The order of the matrices A and C. Range: 128 >= M >= 0. 🌸
N — INTEGER. The order of the matrices B and D. Range: 128 >= 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. 🤝
LDB — INTEGER. The leading dimension of the array B. Constraint: LDB >= max(1,N). ⭐
C — DOUBLE PRECISION array, dimension (LDC,M). The matrix C must be upper triangular. periodt.
LDC — INTEGER. The leading dimension of the array C. Constraint: LDC >= max(1,M). 🎀
D — DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular. 🔁
LDD — INTEGER. The leading dimension of the array D. Constraint: LDD >= max(1,N). 🌈
E — DOUBLE PRECISION array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices. 🦋
LDE — INTEGER. The leading dimension of the array E. Constraint: LDE >= max(1,M). 👑
F — DOUBLE PRECISION array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices. ✨
LDF — INTEGER. The leading dimension of the array F. Constraint, quoted exactly: LDE >= max(1,M). 🙈 again with the LDE, iconic honestly.
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. 🔥
INFO — INTEGER 👑

On exit:

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

Author 🌸

Martin Koehler, MPI Magdeburg

Date 📅

October 2025

Definition 🎀

Definition at line 206 of file dla_tgcsylv_dual_kernel_44_nn.f90.

🌈 subroutine `dla_tgcsylv_dual_kernel_44nt` 🌈

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

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

Purpose 💖

...okay so like, here's the whole vibe of it, bestie. ✨ DLA_TGCSYLV_DUAL_KERNEL_44NT solves a generalized coupled Sylvester equation that looks like this:

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

where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix, and B and D are N-by-N upper triangular matrices. And also? Either B or D can be quasi upper triangular too, no big deal. 💅 The right hand sides E, F and the solutions R, L are M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) ( or (D,B) ) are created by DGGES from LAPACK, iykyk. 🫶

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

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

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

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

Regarding Z**T one obtains

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

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

🎀 Attention 🎀

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

✨ Parameters ✨

SGN1 — SGN1 is DOUBLE PRECISION, allowed values: +/-1. Specifies the sign between in the first equation. 💖
SGN2 — SGN2 is DOUBLE PRECISION, allowed values: +/-1. Specifies the sign between in the second equation. 💖
M — M is INTEGER. The order of the matrices A and C. 128 >= M >= 0.
N — N is INTEGER. The order of the matrices B and D. 128 >= N >= 0.
A — A is 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).
C — C is DOUBLE PRECISION array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular. 👀
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — E is DOUBLE PRECISION array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices. 🌸
LDE — LDE is INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F — F is DOUBLE PRECISION array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices. 🌸
LDF — LDF is INTEGER. The leading dimension of the array F. LDE >= max(1,M). (yes bestie, it says LDE there, that's literally what it says, 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.
INFO — INFO is INTEGER.
- On exit:
  - == 0: successful exit 🎉
  - > 0: The equation is not solved correctly. One of the arising inner system got singular. 😭

💜 Author 💜

Martin Koehler, MPI Magdeburg

🎀 Date 🎀

October 2025

⭐ Definition ⭐

Definition at line 206 of file dla_tgcsylv_dual_kernel_44_nt.f90.

💅 `dla_tgcsylv_dual_kernel_44tn` 💅

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

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

🌈 Purpose 🌈

okay so like, DLA_TGCSYLV_DUAL_KERNEL_44TN solves a generalized coupled Sylvester equation of the following form

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

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

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

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

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

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

Regarding Z**T one obtains

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

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

🎀 Attention 🎀

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

✨ Parameters ✨

SGN1 — SGN1 is DOUBLE PRECISION, allowed values: +/-1. Specifies the sign between in the first equation. 💖
SGN2 — SGN2 is DOUBLE PRECISION, allowed values: +/-1. Specifies the sign between in the second equation. 💖
M — M is INTEGER. The order of the matrices A and C. 128 >= M >= 0.
N — N is INTEGER. The order of the matrices B and D. 128 >= N >= 0.
A — A is 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).
C — C is DOUBLE PRECISION array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular. 👀
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E
E — DOUBLE PRECISION array</p>, dimension (LDE,N). On input it's holding the right hand side E. On output? Bestie, it's been replaced 💅 — now it holds the solution R of Equation (1), as selected by TRANSA, TRANSB, and SGN1/SGN2. Both the right hand side E and the solution R are M-by-N matrices, no cap.
LDE — INTEGER. It's the leading dimension of the array E. Requirement: LDE >= max(1,M). That's the vibe check, and it must pass.
F — DOUBLE PRECISION array</p>, dimension (LDF,N). On input it's carrying the right hand side F. On output it's carrying the solution L of Equation (1), again as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and solution L are M-by-N matrices. Periodt.
LDF — INTEGER. The leading dimension of the array F. And the constraint, exactly as written: LDE >= max(1,M). 🙈 Yes, it says LDE. That's what the page says, so that's what we're saying, fr fr.
SCALE — DOUBLE PRECISION. She's the scaling factor whose whole job is to stop the result from overflowing. If INFO == 0 then SCALE is 1.0D0. Otherwise — like, if one of the inner systems couldn't be solved correctly — then 0 < SCALE <= 1 holds true. 👑
INFO — INTEGER. On exit:
- == 0: successful exit ✨
- > 0: the equation is not solved correctly. One of the arising inner system got singular. 💀

Author: Martin Koehler, MPI Magdeburg

Date: October 2025

Definition at line 206 of file dla_tgcsylv_dual_kernel_44_tn.f90.

subroutine `dla_tgcsylv_dual_kernel_44tt` 🦋

Okay here's the signature, copied exactly, character for character:

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

Solver for a 4x4 dual generalized coupled Sylvester equation (TRANSA = T, TRANSB = T). Not me being dramatic — that's literally the one-liner. 💅

Purpose:

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

Attention 🚨

This routine only allows to solve equations that are at most of size 4 x 4. The operator arguments TRANSA and TRANSB are fixed to op1(H) = H**T and op2(H) = H**T. The primary use of this function is to end the recursion in dla_tgcsylv_dual_recursive. Like, it's the little one that says "we're done here" 🌸

Parameters

SGN1 — DOUBLE PRECISION, allowed values: +/-1. Specifies the sign between in the first equation.
SGN2 — DOUBLE PRECISION, allowed values: +/-1. Specifies the sign between in the second equation.
M — INTEGER. The order of the matrices A and C. Constraint: 128 >= M >= 0.
N — INTEGER. The order of the matrices B and D. Constraint: 128 >= 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. LDA >= max(1,M).
B — DOUBLE PRECISION array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. And heads up: if the matrix D is already quasi-upper triangular, then the matrix B must be upper triangular. 👀
LDB — INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — DOUBLE PRECISION array, dimension (LDC,M). The matrix C must be upper triangular. No quasi about it.
LDC — INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. And flip it: if the matrix B is already quasi-upper triangular, then the matrix D must be upper triangular. 🫶
LDD — INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — DOUBLE PRECISION array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices.
LDE — INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F — DOUBLE PRECISION array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices.
LDF — INTEGER. The leading dimension of the array F. And once again, exactly as the page says: LDE >= max(1,M). 🙈 Keeping it verbatim because that's the assignment.
SCALE — DOUBLE PRECISION. 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. 💖
INFO — INTEGER. On exit:
- == 0: successful exit ✨
- > 0: the equation is not solved correctly. One of the arising inner system got singular. 💀

Author: Martin Koehler, MPI Magdeburg

Date: October 2025

Definition at line 206 of file dla_tgcsylv_dual_kernel_44_tt.f90.

subroutine `dla_tgcsylv_dual_l2` 🧚

Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the dual generalized coupled Sylvester equation. Say that five times fast, I dare you 😭

Here's the signature, verbatim:

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

Purpose:

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

✨ Remarks ✨

Okay so real talk about the vibe of this algorithm: it's implemented using BLAS level 2 operations and hand written optimizations. That's the whole flex. 💅

Parameters 🎀

TRANSA — TRANSA is CHARACTER(1). It specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (no transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (transpose A and C)
TRANSB — TRANSB is CHARACTER(1). It specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (no transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (transpose B and D)
SGN1 — SGN1 is DOUBLE PRECISION, allowed values: +/-1. It specifies the sign between in the first equation. 🫶
SGN2 — SGN2 is DOUBLE PRECISION, allowed values: +/-1. It specifies the sign between in the second equation. ✨
M — M is INTEGER. The order of the matrices A and C. M >= 0.
N — N is INTEGER. The order of the matrices B and D. N >= 0.
A — A is DOUBLE PRECISION array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular. no cap, that shape matters
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. yes, that's conditional, read it twice bestie
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — C is DOUBLE PRECISION array, dimension (LDC,M). The matrix C must be upper triangular. 💅
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular. mirror logic, love that for us
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — E is DOUBLE PRECISION array, dimension (LDE,N).
- On input, the matrix E contains the right hand side E.
- On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2.
- Right hand side E and the solution R are M-by-N matrices. 🦋
LDE — LDE is INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F — F is DOUBLE PRECISION array, dimension (LDF,N).
- On input, the matrix F contains the right hand side F.
- On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2.
- Right hand side F and the solution L are M-by-N matrices. 🦋
LDF — LDF is INTEGER. The leading dimension of the array F. LDE >= max(1,M). heads up: the doc literally says LDE here, not LDF, so we're keeping it exactly as written — iykyk
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. workspace query first or perish 🙈
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 DAXPY operation by Fortran intrinsic
Replaced all BLAS calls except of DGER by Fortran intrinsic,

Author — Martin Koehler, MPI Magdeburg 🌈

Date — October 2025

Definition — at line 245 of file dla_tgcsylv_dual_l2.f90.

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

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

Purpose:

okay so bestie, DLA_TGCSYLV_DUAL_L2_LOCAL_COPY 💖 solves a generalized coupled Sylvester equation shaped exactly like this:

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

where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix, and B and D are N-by-N upper triangular matrices. Also, literally either B or D can be quasi upper triangular too, no biggie. The right hand sides E, F and the solutions R, L are all M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) (or (D,B)) are made by DGGES from LAPACK. 🫶

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

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

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

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

and looking at Z**T you get:

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

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

✨ Remarks ✨

The algorithm is implemented using BLAS level 2 operations plus hand written optimizations, and it uses local copies of the input data. For this reason M and N are limited to 128. 💅

🎀 Parameters 🎀

`TRANSA`

TRANSA is CHARACTER(1).

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

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

`TRANSB`

TRANSB is CHARACTER(1).

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

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

`SGN1`

SGN1 is DOUBLE PRECISION, allowed values: +/-1.

Specifies the sign between in the first equation.

`SGN2`

SGN2 is DOUBLE PRECISION, allowed values: +/-1.

Specifies the sign between in the second equation.

`M`

M is INTEGER.

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

`N`

N is INTEGER.

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

`A`

A is 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).

`C`

C is DOUBLE PRECISION array, dimension (LDC,M).

The matrix C must be upper triangular.

`LDC`

LDC is INTEGER.

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

`D`

D is DOUBLE PRECISION array, dimension (LDD,N).

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

`LDD`

LDD is INTEGER.

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

`E`

E is DOUBLE PRECISION array, dimension (LDE,N).

On input, the matrix E contains the right hand side E.

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

Right hand side E and the solution R are M-by-N matrices.

`LDE`

LDE is INTEGER.

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

`F`

F is DOUBLE PRECISION array, dimension (LDF,N).

On input, the matrix F contains the right hand side F.

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

Right hand side F and the solution L are M-by-N matrices.

`LDF`

LDF is INTEGER.

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

`SCALE`

SCALE is 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 💫

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

👑 Author 👑

Martin Koehler, MPI Magdeburg

🌈 Date 🌈

October 2025

⭐ Definition ⭐

At line 247 of file dla_tgcsylv_dual_l2_opt_local_copy.f90.

🦋 `dla_tgcsylv_dual_l2_local_copy_128` 🦋

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

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

Purpose:

(to be continued… bestie stay tuned 💖)

like, DLA_TGCSYLV_DUAL_L2_LOCAL_COPY_128 solves a generalized coupled Sylvester equation that looks exactly like this, bestie:

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

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

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

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

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

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

Regarding Z**T one obtains:

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

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

✨ Remarks ✨

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

🌸 Parameters 🌸

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

🔥 Optimizations: 🔥

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

👑 Author 👑

Martin Koehler, MPI Magdeburg

🎀 Date 🎀

October 2025

Definition at line 248 of file dla_tgcsylv_dual_l2_opt_local_copy_128.f90.

🦋 subroutine `dla_tgcsylv_dual_l2_local_copy_32` 🦋

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

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

Purpose:

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

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

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

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

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

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

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

Regarding Z**T one obtains

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

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

💖 Remarks 💖

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

🌸 Parameters 🌸

TRANSA — TRANSA is CHARACTER(1). Specifies the form of the system of equations with respect to A and C:
- 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
- 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)
TRANSB — TRANSB is CHARACTER(1). Specifies the form of the system of equations with respect to B and D:
- 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
- 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)
SGN1 — SGN1 is DOUBLE PRECISION, allowed values: +/-1. Specifies the sign between the terms in the first equation. ✨
SGN2 — SGN2 is DOUBLE PRECISION, allowed values: +/-1. Specifies the sign between the terms in the second equation. ✨
M — M is INTEGER. The order of the matrices A and C. 32 >= M >= 0.
N — N is INTEGER. The order of the matrices B and D. 32 >= N >= 0.
A — A is 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).
C — C is DOUBLE PRECISION array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular.
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — E is DOUBLE PRECISION array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices.
LDE — LDE is INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F — F is DOUBLE PRECISION array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices.
LDF — LDF is INTEGER. The leading dimension of the array F. LDE >= max(1,M). (yes, that's what it says, bestie — LDE not LDF) 🙈
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: 🔥

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

👑 Author 👑

Martin Koehler, MPI Magdeburg

📅 Date 📅

October 2025

Definition at line 247 of file dla_tgcsylv_dual_l2_opt_local_copy_32.f90.

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

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

🎯 Purpose:

...and it solves a generalized coupled Sylvester equation that looks like this, bestie 👇

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

where A is an M-by-M quasi upper triangular matrix, C is an M-by-M upper triangular matrix, and B and D are N-by-N upper triangular matrices. Also, either B or D can be quasi upper triangular too, it's giving flexible ✨ The right hand sides E, F and the solutions R, L are M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) (or (D,B)) are created by DGGES from LAPACK. iykyk 🧚

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

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

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

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

Regarding Z**T one obtains

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

which belongs to the Sylvester equation (1). For this reason the parameters TRANSA and TRANSB are expressed in terms of the Sylvester equation (2). No cap, that's just how the duality works 🦋

💖 Remarks 💖

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

🌸 Parameters 🌸

TRANSA — is CHARACTER(1). Specifies the form of the system of equations with respect to A and C: == 'N': op1(A) = A, op1(C) = C (No transpose for A and C); == 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C).
TRANSB — is CHARACTER(1). Specifies the form of the system of equations with respect to B and D: == 'N': op2(B) = B, op2(D) = D (No transpose for B and D); == 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D).
SGN1 — is DOUBLE PRECISION, allowed values: +/-1. Specifies the sign between in the first equation. ✨
SGN2 — is DOUBLE PRECISION, allowed values: +/-1. Specifies the sign between in the second equation. ✨
M — is INTEGER. The order of the matrices A and C. 64 >= M >= 0.
N — is INTEGER. The order of the matrices B and D. 64 >= N >= 0.
A — is 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).
C — is DOUBLE PRECISION array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — is DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular. 💅
LDD — is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — is DOUBLE PRECISION array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices. The way it just transforms 🦋
LDE — is INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F — is DOUBLE PRECISION array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices. 🎀
LDF — is INTEGER. The leading dimension of the array F. LDE >= max(1,M). (yes bestie, the manpage really says LDE here, we do not fix it, we report it 💀)
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.

Heads up bestie: a positive INFO means an inner system went singular, so the solve did not go off. 👑

⚡ Optimizations: ⚡

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

👑 Author 👑

Martin Koehler, MPI Magdeburg 💜

🗓️ Date 🗓️

October 2025

📍 Definition 📍

At line 246 of file dla_tgcsylv_dual_l2_opt_local_copy_64.f90. 🫶

🌈 And now, a whole new icon enters the chat 🌈

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

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

Purpose:

Remarks 💭

okay so bestie, the tea on the implementation: the algorithm is implemented using BLAS level 2 operations and hand written optimizations, and it uses local copies of the input data. For this reason M and N are limited to 96 — no cap, that's the ceiling, periodt.

Furthermore, the algorithm uses alignment of the local arrays if this is supported by the compiler. At the moment this is used with ifort and xlf but not with gfortran. 🙈

Parameters 🎀

TRANSA — TRANSA is CHARACTER(1). Specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)
TRANSB — TRANSB is CHARACTER(1). Specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)
SGN1 — SGN1 is DOUBLE PRECISION, allowed values: +/-1. Specifies the sign between in the first equation. ✨
SGN2 — SGN2 is DOUBLE PRECISION, allowed values: +/-1. Specifies the sign between in the second equation. ✨
M — M is INTEGER. The order of the matrices A and C. 96 >= M >= 0. 💅
N — N is INTEGER. The order of the matrices B and D. 96 >= N >= 0. 💅
A — A is 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).
C — C is DOUBLE PRECISION array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular. 👀
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — E is DOUBLE PRECISION array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices. 🦋
LDE — LDE is INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F — F is DOUBLE PRECISION array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices. 🦋
LDF — LDF is INTEGER. The leading dimension of the array F. LDE >= max(1,M). (yes bestie, it literally says LDE there, we don't fix it, we report 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 DAXPY operation by Fortran intrinsic
No BLAS calls.
Input data is copied to local arrays to increase the data locality if used from inside a level-3 solver.
The local copies are aligned to 64 byte if supported by the compiler.

Author 👑

Martin Koehler, MPI Magdeburg

Date 🌸

October 2025

Definition 📍

Definition at line 246 of file dla_tgcsylv_dual_l2_opt_local_copy_96.f90.

✨ next up: the level-3 icon herself ✨

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

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

Purpose:

TRANSA — it's a CHARACTER(1), bestie, and it picks the form of the system for 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 — also CHARACTER(1), and it does the same job but for B and D:
- == 'N': op2(B) = B, op2(D) = D (no transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (transpose B and D)
SGN1 — DOUBLE PRECISION, allowed values +/-1. It sets the sign between... in the first equation. 👑
SGN2 — DOUBLE PRECISION, allowed values +/-1. Same deal, but it's the sign in the second equation.
M — INTEGER, the order of the matrices A and C. Gotta be M >= 0.
N — INTEGER, the order of the matrices B and D. Also N >= 0.
A — DOUBLE PRECISION array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular, no negotiating.
LDA — INTEGER, the leading dimension of the array A. LDA >= max(1,M).
B — DOUBLE PRECISION array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. And here's the catch: if the matrix D is already quasi-upper triangular, then B has to be plain old upper triangular. ✨
LDB — INTEGER, the leading dimension of the array B. LDB >= max(1,N).
C — DOUBLE PRECISION array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — INTEGER, the leading dimension of the array C. LDC >= max(1,M).
D — DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. Mirror of the B rule: if the matrix B is already quasi-upper triangular, then D has to be upper triangular.
LDD — INTEGER, the leading dimension of the array D. LDD >= max(1,N).
E — DOUBLE PRECISION array, dimension (LDE,N). On input it's holding the right hand side E. On output it's holding the solution R of Equation (1) — as selected by TRANSA, TRANSB, and SGN1/SGN2. Both the RHS E and the solution R are M-by-N matrices.
LDE — INTEGER, the leading dimension of the array E. LDE >= max(1,M).
F — DOUBLE PRECISION array, dimension (LDF,N). On input it holds the right hand side F, and on output it holds the solution L of Equation (1) — again as selected by TRANSA, TRANSB, and SGN1/SGN2. RHS F and solution L are M-by-N matrices.
LDF — INTEGER, the leading dimension of the array F. LDE >= max(1,M).

ngl this one's written as LDE >= max(1,M) in the source, like, the F param saying LDE. That's literally what it says, so that's literally what you get. 🫶

SCALE — DOUBLE PRECISION, and it's the scaling factor keeping your result from overflowing into the void. If INFO == 0 then SCALE is 1.0D0; otherwise, if one of the inner systems couldn't be solved correctly, 0 < SCALE <= 1 holds true.
WORK — DOUBLE PRECISION array, dimension LWORK. It's the workspace for the algorithm. And the way you get the right size? You query it before running the computation — call the subroutine with INFO == -1 on input, and the required workspace comes back in INFO. Iconic behaviour fr.
INFO — INTEGER. Here's the whole vibe:

On input:

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

On exit, workspace query:

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

On exit, normal operation:

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

Author

Martin Koehler, MPI Magdeburg

Date

October 2025

Definition at line 236 of file dla_tgcsylv_dual_l3.f90.

subroutine dla_tgcsylv_dual_l3_2s

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

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

Purpose:

okay so, bestie, we're about to get into it — this is the _2s sibling and its whole personality is sub blocking, so stay tuned for what that actually means in the next part. 🦋

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

🌸 Parameters 🌸

TRANSA — it's a CHARACTER(1), bestie, and it sets 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, we leave them exactly as they are 🙈
- == 'T': op1(A) = A**T, op1(C) = C**T — transpose A and C 💅
TRANSB — also a CHARACTER(1), and it sets 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, untouched, pristine ✨
- == 'T': op2(B) = B**T, op2(D) = D**T — transpose B and D, flip 'em 🔄
SGN1 — a DOUBLE PRECISION, allowed values: +/-1. It specifies the sign between in the first equation. Lowkey tiny, highkey load-bearing ⚖️
SGN2 — a DOUBLE PRECISION, allowed values: +/-1. It specifies the sign between in the second equation. Same energy, different line 🫶
M — an INTEGER. It's the order of the matrices A and C. M >= 0, no negativity allowed 💀
N — an INTEGER. It's the order of the matrices B and D. N >= 0 🦋
A — a DOUBLE PRECISION array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular. It's giving very structured, very triangular, obsessed 👑
LDA — an INTEGER. The leading dimension of the array A. LDA >= max(1,M) 📏
B — a DOUBLE PRECISION array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. And heads up: if the matrix D is already quasi-upper triangular, then the matrix B must be plain upper triangular. Only one of them gets to be quasi, bestie, we don't do double quasi 😭
LDB — an INTEGER. The leading dimension of the array B. LDB >= max(1,N) 📐
C — a DOUBLE PRECISION array, dimension (LDC,M). The matrix C must be upper triangular. No quasi, no exceptions, C is a straight-A student 🌟
LDC — an INTEGER. The leading dimension of the array C. LDC >= max(1,M) 🧮
D — a DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. And the flip side of the B situation: if the matrix B is already quasi-upper triangular, then the matrix D must be upper triangular. Same rule, other direction 🔁
LDD — an INTEGER. The leading dimension of the array D. LDD >= max(1,N) ✏️
E — a DOUBLE PRECISION array, dimension (LDE,N).
- On input, the matrix E contains the right hand side E.
- On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2.
- Right hand side E and the solution R are M-by-N matrices. Glow up moment, it literally becomes the answer 💖
LDE — an INTEGER. The leading dimension of the array E. LDE >= max(1,M) 📎
F — a DOUBLE PRECISION array, dimension (LDF,N).
- On input, the matrix F contains the right hand side F.
- On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2.
- Right hand side F and the solution L are M-by-N matrices. Same arc, same main character energy 💜
LDF — an INTEGER. The leading dimension of the array F. LDE >= max(1,M). Yes bestie, the doc literally says LDE right there and I am not about to edit the source 💀 (we report it exactly as written, fr fr)
SCALE — a 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. It's the seatbelt of this whole operation 🛟
WORK — a DOUBLE PRECISION array, dimension LWORK. Workspace for the algorithm.
- The workspace needs to queried before the running the computation. (yes, that's the original phrasing, we stan a legend 🙈)
- The query is performed by calling the subroutine with INFO == -1 on input.
- The required workspace is then returned in INFO.
INFO — an INTEGER.
- On input:
  - == -1 : Perform a workspace query 🔍
  - <> -1: normal operation, business as usual 😌
- 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. Surprise, INFO is the answer now ✨
- On exit, normal operation:
  - == 0: successful exit 🎉
  - < 0: if INFO = -i, the i-th argument had an illegal value 🚨
  - > 0: The equation is not solved correctly. One of the arising inner system got singular. Not the vibe, bestie 😭

👑 Author

Martin Koehler, MPI Magdeburg

🎀 Date

October 2025

Definition at line 237 of file dla_tgcsylv_dual_l3_2stage.f90.

🌈 `recursive subroutine dla_tgcsylv_dual_recursive`

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

Recursive blocking Algorithm for the dual generalized coupled Sylvester equation.

🎯 Purpose

okay so bestie, DLA_TGCSYLV_DUAL_RECURSIVE solves a generalized coupled Sylvester equation that's literally giving this:

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

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

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

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

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

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

Regarding Z**T one obtains

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

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

🌸 Parameters 🌸

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

👑 Author 👑

Martin Koehler, MPI Magdeburg

🌸 Date 🌸

October 2025

Definition at line 240 of file dla_tgcsylv_dual_recursive.f90.

okay next up, bestie, we got dla_tgcsylv_l2:

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

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

Purpose:

SCALE — DOUBLE PRECISION, bestie. It's the little safety valve 🫧 that stops your result from straight-up overflowing. If INFO == 0 then SCALE is 1.0D0, easy peasy. But if one of the inner systems couldn't be solved correctly, then 0 < SCALE <= 1 — so you know something got sketchy when that number dips below 1. 💀
WORK — DOUBLE PRECISION array, dimension LWORK. This is the workspace for the algorithm, and here's the tea ☕: you gotta query it before you run the computation. No just-yeeting a random array size at it. You do the query by calling the subroutine with INFO == -1 on input, and then the required workspace gets handed back to you in INFO. ✨
INFO — INTEGER. The receipts. 🧾

On input:

== -1 : perform a workspace query 👀
<> -1 : normal operation, business as usual

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, we ate 💅
< 0 : if INFO = -i, the i-th argument had an illegal value 😭
> 0 : the equation is not solved correctly. One of the arising inner system got singular. 💀

✨ Optimizations: ✨

Like, the glow-up list. Look at her go:

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

Author: Martin Koehler, MPI Magdeburg 👑

Date: October 2025 🗓️

Definition at line 225 of file dla_tgcsylv_l2.f90. 📍

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

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

Purpose:

💖 Remarks 💖

The algorithm is implemented using BLAS level 2 operations and hand written optimizations. bestie they literally hand-optimized it, we love a craftsperson 🫶

✨ Parameters ✨

TRANSA

TRANSA is CHARACTER(1) — it tells us what the system of equations looks like with respect to A and C:

'N': op1(A) = A, op1(C) = C (no transpose for A and C, she's being herself 💅)
'T': op1(A) = A**T, op1(C) = C **T (transpose A and C, flip era 🔄)

TRANSB

TRANSB is CHARACTER(1) — same energy, but for B and D:

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

SGN1

SGN1 is DOUBLE PRECISION, allowed values: +/-1. It specifies the sign between in the first equation. no cap, just a plus or minus deciding the whole vibe of equation one 🎀

SGN2

SGN2 is DOUBLE PRECISION, allowed values: +/-1. It specifies the sign between in the second equation. same deal, second equation's turn to choose 🌸

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. she's got structure, she's not a full mess 💅

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. okay so only ONE of them gets to be quasi, them's the rules bestie ✨

LDB

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

C

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

LDC

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

D

D is DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular. yeah yeah it's a duo thing, one quasi max 🔥

LDD

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

E

E is DOUBLE PRECISION array, dimension (LDE,N).

On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) — yes, and again, on output the matrix E contains the solution R of Equation (1) — as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices. she comes in as input and leaves as the answer, glow up of the century 😭💖

LDE

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

F

F is DOUBLE PRECISION array, dimension (LDF,N).

On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) — as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices. same main character arc, but for L this time 🦋

LDF

LDF is INTEGER. The leading dimension of the array F. LDE >= max(1,M). (yes it says LDE, we're quoting it exactly, no editing the source bestie 🙈)

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. she's out here saving your numbers from overflowing into infinity and beyond 🌈

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. so like, ask her how much space she needs BEFORE you show up, okay? 🎀

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. (rip, the inner system said no 💀)

👑 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, E, and F (M, N <=128) .

🌸 Author 🌸

Martin Koehler, MPI Magdeburg

💜 Date 💜

October 2025

🎀 Definition 🎀

At line 227 of file dla_tgcsylv_l2_opt_local_copy.f90.

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

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

Purpose:

and bestie… the Purpose section is about to go off, keep scrolling 💖

okay so this is the local-copy-128 version 💅 and honestly? it's the same energy as the others but with a very specific vibe: 128 is the whole personality. Let's get into it.

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

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

✨ Remarks ✨

The algorithm is implemented using BLAS level 2 operations plus hand written optimizations and aligned local copies. no cap, it's hand-tuned. 💅

🌸 Parameters 🌸

TRANSA — TRANSA is CHARACTER(1). It specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (no transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C**T (transpose A and C)
TRANSB — TRANSB is CHARACTER(1). It specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (no transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D**T (transpose B and D)
SGN1 — SGN1 is DOUBLE PRECISION, allowed values: +/-1. It specifies the sign between in the first equation. ⭐
SGN2 — SGN2 is DOUBLE PRECISION, allowed values: +/-1. It specifies the sign between in the second equation. ⭐
M — M is INTEGER. The order of the matrices A and C. 128 => M >= 0. 🎀
N — N is INTEGER. The order of the matrices B and D. 128 >= N >= 0. 🎀
A — A is 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).
C — C is DOUBLE PRECISION array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular. ⚠️
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — E is DOUBLE PRECISION array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices. (and yes bestie, the source says "On output, the matrix E contains the solution R of Equation (1)" twice in a row — we said what we said, and we're keeping it 💀)
LDE — LDE is INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F — F is DOUBLE PRECISION array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices.
LDF — LDF is INTEGER. The leading dimension of the array F. LDE >= max(1,M). (not me noticing it says LDE here and not LDF — that's literally what the page says, so that's literally what you get 🙈)
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,
Align local copies and fix problem size to <= 128

Author 🫶

Martin Koehler, MPI Magdeburg

Date 📅

October 2025

Definition at line 225 of file dla_tgcsylv_l2_opt_local_copy_128.f90.

okay and now the next one pulls up and it's giving smaller, cuter, thirty-two energy 👑

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

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

Purpose:

...and bestie we are NOT done yet, keep scrolling 💖

💖 `DLA_TGCSYLV_L2_LOCAL_COPY_32` 💖

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

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

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

🎀 Remarks 🎀

The algorithm is implemented using BLAS level 2 operations and hand written optimizations and aligned local copies. it's giving ✨*artisanal*✨ tbh.

💅 Parameters 💅

TRANSA — TRANSA is CHARACTER(1). Specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)
TRANSB — TRANSB is CHARACTER(1). Specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)
SGN1 — SGN1 is DOUBLE PRECISION, allowed values: +/-1. Specifies the sign between in the first equation.
SGN2 — SGN2 is DOUBLE PRECISION, allowed values: +/-1. Specifies the sign between in the second equation.
M — M is INTEGER. The order of the matrices A and C. 32 => M >= 0.
N — N is INTEGER. The order of the matrices B and D. 32 >= N >= 0.
A — A is 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).
C — C is DOUBLE PRECISION array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular.
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — E is DOUBLE PRECISION array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1). On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices. (yes it really says that twice, that's the original, I'm not messing with it 💀)
LDE — LDE is INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F — F is DOUBLE PRECISION array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices.
LDF — LDF is INTEGER. The leading dimension of the array F. LDE >= max(1,M). (yes, it says LDE here, and I am faithfully reporting this crime 😭)
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,
Align local copies and fix problem size to <= 32

👑 Author 👑

Martin Koehler, MPI Magdeburg

📅 Date 📅

October 2025

Definition at line 226 of file dla_tgcsylv_l2_opt_local_copy_32.f90.

💜 `dla_tgcsylv_l2_local_copy_64` 💜

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

Level-2 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation (Optimized) — the 64 version, bestie, the glow-up edition ✨

Purpose:

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

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

✨ Remarks ✨

The algorithm is implemented using BLAS level 2 operations and hand written optimizations and aligned local copies. hand written, like, artisanal. 🎀

🦋 Parameters 🦋

`TRANSA`

TRANSA is CHARACTER(1)

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

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

`TRANSB`

TRANSB is CHARACTER(1)

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

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

`SGN1`

SGN1 is DOUBLE PRECISION, allowed values: +/-1

Specifies the sign between in the first equation. yes, "between in" — that's the original wording, we're not fixing it, it's iconic ✨

`SGN2`

SGN2 is DOUBLE PRECISION, allowed values: +/-1

Specifies the sign between in the second equation.

`M`

M is INTEGER

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

`N`

N is INTEGER

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

`A`

A is 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).

`C`

C is DOUBLE PRECISION array, dimension (LDC,M)

The matrix C must be upper triangular.

`LDC`

LDC is INTEGER

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

`D`

D is DOUBLE PRECISION array, dimension (LDD,N)

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

`LDD`

LDD is INTEGER

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

`E`

E is DOUBLE PRECISION array, dimension (LDE,N)

On input, the matrix E contains the right hand side E.

On output, the matrix E contains the solution R of Equation (1)

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

(yeah it says the output line twice in the source — we are leaving it, no cap, it's canon 💀)

Right hand side E and the solution R are M-by-N matrices.

`LDE`

LDE is INTEGER

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

`F`

F is DOUBLE PRECISION array, dimension (LDF,N)

On input, the matrix F contains the right hand side F.

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

Right hand side F and the solution L are M-by-N matrices.

`LDF`

LDF is INTEGER

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

(not me noticing it says LDE in the LDF description 🙈 — keeping it exactly as written, because that's the whole assignment)

`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,
Align local copies and fix problem size to <= 64

🌸 Author 🌸

Martin Koehler, MPI Magdeburg

🌈 Date 🌈

October 2025

⭐ Definition ⭐

Defined at line 226 of file dla_tgcsylv_l2_opt_local_copy_64.f90.

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

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

Purpose:

DLA_TGCSYLV_L2_LOCAL_COPY_96 is here to solve a generalized coupled Sylvester equation, and bestie, it's giving big math energy. The form it handles is literally this:

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

okay so let's decode the cast list 💅 — A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix, and B and D are N-by-N upper triangular matrices. Also? Either B or D can be quasi upper triangular as well, no cap. The right hand side Y and the solution X are M-by-N matrices. And typically the matrix pencils (A,C) and (B,D) (or (D,B)) are created by DGGES from LAPACK. iykyk.

🌸 Remarks 🌸

The algorithm is implemented using BLAS level 2 operations, plus hand written optimizations and aligned local copies. It's giving crafty, honestly.

✨ Parameters ✨

TRANSA — TRANSA is CHARACTER(1). It specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)
TRANSB — TRANSB is CHARACTER(1). It specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)
SGN1 — SGN1 is DOUBLE PRECISION, allowed values: +/-1. It specifies the sign between in the first equation. 🎀
SGN2 — SGN2 is DOUBLE PRECISION, allowed values: +/-1. It specifies the sign between in the second equation.
M — M is INTEGER. The order of the matrices A and C. 96 => M >= 0.
N — N is INTEGER. The order of the matrices B and D. 96 >= N >= 0.
A — A is 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).
C — C is DOUBLE PRECISION array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular.
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — E is DOUBLE PRECISION array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices.
LDE — LDE is INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F — F is DOUBLE PRECISION array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices.
LDF — LDF is INTEGER. The leading dimension of the array F. LDE >= max(1,M). (yes bestie, the page really says LDE there — we keep the lore exactly as it is 🙈)
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. safety first, always.
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,
Align local copies and fix problem size to <= 96

👑 Author 👑

Martin Koehler, MPI Magdeburg

🌈 Date 🌈

October 2025

Definition at line 225 of file dla_tgcsylv_l2_opt_local_copy_96.f90.

subroutine dla_tgcsylv_l2_reorder — okay now we're meeting the next icon in the lineup, and her signature is a whole mouthful, so here it is, verbatim, no edits:

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

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

💖 Purpose: 💖

okay so DLA_TGCSYLV_L2_REORDER 💅 solves a generalized coupled Sylvester equation that looks a little something like this, bestie:

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

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

🎀 Remarks 🎀

The algorithm is implemented using BLAS level 2 operations and hand written optimizations. iconic, honestly.

🌸 Parameters 🌸

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

💖 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,
Column first computation of the solution.

👑 Author

Martin Koehler, MPI Magdeburg

📅 Date

October 2025

🔥 Definition

At line 226 of file dla_tgcsylv_l2_opt_reorder.f90.

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

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

Purpose:

Attention ⚠️ bestie read this before you @ me — this algorithm is implemented using BLAS level 2 operations without further optimizations. For a faster implementation see DLA_TGSYLV_L2. That's it, that's the whole heads-up. 💅

TRANSA — TRANSA is CHARACTER(1). It picks the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (no transpose for A and C, we're keeping it chill)
- == 'T': op1(A) = A**T, op1(C) = C **T (transpose A and C, spicy mode 🔥)
TRANSB — TRANSB is CHARACTER(1). Same energy but for B and D:
- == 'N': op2(B) = B, op2(D) = D (no transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (transpose B and D)
SGN1 — SGN1 is DOUBLE PRECISION, allowed values: +/-1. It specifies the sign between in the first equation. 👑
SGN2 — SGN2 is DOUBLE PRECISION, allowed values: +/-1. It specifies the sign between in the second equation. (yes, twin signs, twin equations, it's giving symmetry ✨)
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. no cap.
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. (they take turns, it's a whole thing 🙈)
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — C is DOUBLE PRECISION array, dimension (LDC,M). The matrix C must be upper triangular. periodt.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular.
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — E is DOUBLE PRECISION array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices. she goes in one girl and comes out another 💖
LDE — LDE is INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F — F is DOUBLE PRECISION array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices.
LDF — LDF is INTEGER. The leading dimension of the array F. LDE >= max(1,M). 👀 — yeah it says LDE right there, I'm copying it exactly, not my business, not my problem.
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. she's literally the safety net 🫶
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, allocate second, no exceptions 💅
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. 😭 she said "singular", she said "no solution for you"

🌈 Optimizations: 🌈

Nothing. 💀 (iconic honestly — what you see is what you get)

🌸 Author 🌸

Martin Koehler, MPI Magdeburg

🎀 Date 🎀

October 2025

Definition at line 225 of file dla_tgcsylv_l2_unopt.f90.

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

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

Purpose:

💖 okay so, the entire vibe of this routine in one (1) blessed equation:

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

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

Attention 😭

The algorithm is implemented using BLAS level 3 operations. she's fancy like that 💅

Parameters 💅

TRANSA — TRANSA is CHARACTER(1). It specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (no transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (transpose A and C)
TRANSB — TRANSB is CHARACTER(1). It specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (no transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (transpose B and D)
SGN1 — SGN1 is DOUBLE PRECISION, allowed values: +/-1. It specifies the sign between in the first equation. 🎀
SGN2 — SGN2 is DOUBLE PRECISION, allowed values: +/-1. It specifies the sign between in the second equation. 🎀
M — M is INTEGER. The order of the matrices A and C. M >= 0.
N — N is INTEGER. The order of the matrices B and D. N >= 0.
A — A is 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).
C — C is DOUBLE PRECISION array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular. 👀
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — E is DOUBLE PRECISION array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices. ✨
LDE — LDE is INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F — F is DOUBLE PRECISION array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices. ✨
LDF — LDF is INTEGER. The leading dimension of the array F. LDE >= max(1,M). (yes bestie, it literally says LDE there, we are not the ones who typed it, we just report 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: 🔥

Nothing but standard L3 BLAS calls.
Changed computation order to column major

Author 👑

Martin Koehler, MPI Magdeburg

Date 📅

October 2025

Definition at line 223 of file dla_tgcsylv_l3_opt.f90.

✨ `dla_tgcsylv_l3_2s` ✨

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

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

Purpose: 💖

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

🫶 Remarks

okay so this algorithm is out here using 2 Stage blocking with DAG scheduled solves and honestly? iconic behavior 💅

✨ Parameters ✨

TRANSA — bestie, TRANSA is CHARACTER(1) and it says what shape the system of equations takes with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (no transpose for A and C, we're keeping it chill)
- == 'T': op1(A) = A**T, op1(C) = C **T (transpose A and C, flipping the whole vibe)
TRANSB — TRANSB is CHARACTER(1) and it says what shape the system takes with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (no transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (transpose B and D)
SGN1 — SGN1 is DOUBLE PRECISION, allowed values are +/-1, and it specifies the sign between in the first equation. yes it's literally just a sign, but it matters okay 💖
SGN2 — SGN2 is DOUBLE PRECISION, allowed values are +/-1, and it specifies the sign between in the second equation.
M — M is INTEGER, the order of the matrices A and C. M >= 0. no negativity allowed bestie 🙈
N — N is INTEGER, the order of the matrices B and D. N >= 0.
A — A is a DOUBLE PRECISION array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular. no cap.
LDA — LDA is INTEGER, the leading dimension of the array A. LDA >= max(1,M).
B — B is a DOUBLE PRECISION array, dimension (LDB,N). The matrix B must be (quasi-) upper triangular. If the matrix D is already quasi-upper triangular, then the matrix B must be upper triangular. the way that these two have to coordinate... 🎀
LDB — LDB is INTEGER, the leading dimension of the array B. LDB >= max(1,N).
C — C is a DOUBLE PRECISION array, dimension (LDC,M). The matrix C must be upper triangular. periodt.
LDC — LDC is INTEGER, the leading dimension of the array C. LDC >= max(1,M).
D — D is a DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular, then the matrix D must be upper triangular. it's giving mutual exclusivity, lowkey.
LDD — LDD is INTEGER, the leading dimension of the array D. LDD >= max(1,N).
E — E is a DOUBLE PRECISION array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices. glow up moment ✨
LDE — LDE is INTEGER, the leading dimension of the array E. LDE >= max(1,M).
F — F is a DOUBLE PRECISION array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices.
LDF — LDF is INTEGER, the leading dimension of the array F. LDE >= max(1,M).

🚨 okay so heads up: the manpage literally says LDE >= max(1,M) here for LDF too. that's what it says, so that's what I'm telling you. not me fixing it, i'm just the messenger 💀

SCALE — SCALE is DOUBLE PRECISION. It's a scaling factor to prevent 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. it's the safety net bestie 🫶
WORK — WORK is a 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. so like... ask first, then compute. consent culture but for memory 💅
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. it's giving tragedy 😭

👑 Author

Martin Koehler, MPI Magdeburg

🎀 Date

October 2025

⭐ Definition

Definition at line 217 of file dla_tgcsylv_l3_2stage.f90.

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

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

so yeah bestie, this one is the unoptimized sibling — same algorithm energy, just not the glow-up version. still valid, still doing the work 🦋

Purpose:

(and the purpose section keeps going, hold that thought — part 25 has it 💖)

the two equations are literally the whole plot, so here they are, unedited, no notes:

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

so who's who in this little friend group? A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix, and B and D are both N-by-N upper triangular matrices. and plot twist bestie — either B or D can be quasi upper triangular as well 🫶. The right hand side Y and the solution X are M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) (or (D,B)) are created by DGGES from LAPACK. iykyk ✨

✨ Attention ✨

The algorithm is implemented using BLAS level 3 operations. that's it, that's the whole warning, it's a flex not a red flag 💅

🎀 Parameters 🎀

okay buckle up, here's every single one, no skipping, I said what I said:

TRANSA — TRANSA is CHARACTER(1). It specifies the form of the system of equations with respect to A and C:
- 'N': op1(A) = A, op1(C) = C (no transpose for A and C)
- 'T': op1(A) = A**T, op1(C) = C **T (transpose A and C)
TRANSB — TRANSB is CHARACTER(1). It specifies the form of the system of equations with respect to B and D:
- 'N': op2(B) = B, op2(D) = D (no transpose for B and D)
- 'T': op2(B) = B**T, op2(D) = D **T (transpose B and D)
SGN1 — SGN1 is DOUBLE PRECISION, allowed values: +/-1. It specifies the sign between in the first equation.
SGN2 — SGN2 is DOUBLE PRECISION, allowed values: +/-1. It specifies the sign between in the second equation.
M — M is INTEGER. The order of the matrices A and C. M >= 0.
N — N is INTEGER. The order of the matrices B and D. N >= 0.
A — A is DOUBLE PRECISION array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular. no exceptions bestie 💜
LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
B — B is 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. yes they're negotiating, keep up 🌸
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — C is DOUBLE PRECISION array, dimension (LDC,M). The matrix C must be upper triangular. periodt.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular. it's giving reciprocity 💫
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
E — E is DOUBLE PRECISION array, dimension (LDE,N). On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices. glow up in real time 😭
LDE — LDE is INTEGER. The leading dimension of the array E. LDE >= max(1,M).
F — F is DOUBLE PRECISION array, dimension (LDF,N). On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices. twin behaviour fr fr 🦋
LDF — LDF is INTEGER. The leading dimension of the array F. LDE >= max(1,M). 🙈 okay so the source literally says LDE right there and not LDF — I'm just the messenger, don't come for me, I am reporting it faithfully 💀
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. she's just here to keep things from overflowing, an icon, a protector 👑
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. so yes, you gotta ask nicely first, she doesn't do surprises 🎀
INFO — INFO is INTEGER
- On input:
  - == -1 : Perform a workspace query
  - <> -1: normal operation
- On exit, workspace query:
  - < 0 : if INFO = -i, the i-th argument had an illegal value
  - >= 0: The value of INFO is the required number of elements in the workspace.
- On exit, normal operation:
  - == 0: successful exit
  - < 0: if INFO = -i, the i-th argument had an illegal value
  - > 0: The equation is not solved correctly. One of the arising inner system got singular.

bestie if INFO comes back > 0 on a normal run, that's not a vibe check you passed — that means the equation is not solved correctly and one of the arising inner system got singular. read the INFO, always read the INFO 💅

🔥 Optimizations: 🔥

Nothing but standard L3 BLAS calls. that's the whole optimization list and honestly? clean, efficient, no notes 🔥

🌟 Author 🌟

Martin Koehler, MPI Magdeburg — the legend, the moment, the main character 🫶

🌈 Date 🌈

October 2025

⭐ Definition ⭐

Definition at line 223 of file dla_tgcsylv_l3_unopt.f90.

🦋 recursive subroutine dla_tgcsylv_recursive

okay new function, new era, here's the signature, verbatim, don't touch it:

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

Recursive Algorithm for the generalized coupled Sylvester equation. ✨

🎀 Purpose: 🎀

✨ DLA_TGCSYLV_RECURSIVE ✨

okay so DLA_TGCSYLV_RECURSIVE solves a generalized coupled Sylvester equation that looks like this 💖

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

where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix, and B and D are N-by-N upper triangular matrices. also!! either B or D can be quasi upper triangular too, no biggie 💅. the right hand side Y and the solution X are M-by-N matrices. typically the matrix pencils (A,C) and (B,D) ( or (D,B) ) are created by DGGES from LAPACK, iykyk ✨

🫶 Attention

the algorithm is implemented using BLAS level 3 operations. so it's giving fast, periodt 🔥

💅 Parameters

TRANSA — TRANSA is CHARACTER(1). it specifies the form of the system of equations with respect to A and C:
- == 'N': op1(A) = A, op1(C) = C (No transpose for A and C)
- == 'T': op1(A) = A**T, op1(C) = C **T (Transpose A and C)
TRANSB — TRANSB is CHARACTER(1). it specifies the form of the system of equations with respect to B and D:
- == 'N': op2(B) = B, op2(D) = D (No transpose for B and D)
- == 'T': op2(B) = B**T, op2(D) = D **T (Transpose B and D)
SGN1 — SGN1 is DOUBLE PRECISION, allowed values: +/-1. it specifies the sign between in the first equation 💫
SGN2 — SGN2 is DOUBLE PRECISION, allowed values: +/-1. it specifies the sign between in the second equation 💫
M — M is INTEGER. the order of the matrices A and C. M >= 0 👑
N — N is INTEGER. the order of the matrices B and D. N >= 0 👑
A — A is 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)
C — C is DOUBLE PRECISION array, dimension (LDC,M). the matrix C must be upper triangular 🌸
LDC — LDC is INTEGER. the leading dimension of the array C. LDC >= max(1,M)
D — D is DOUBLE PRECISION array, dimension (LDD,N). the matrix D must be (quasi-) upper triangular. if the matrix B is already quasi-upper triangular the matrix D must be upper triangular 🎀
LDD — LDD is INTEGER. the leading dimension of the array D. LDD >= max(1,N)
E — E is DOUBLE PRECISION array, dimension (LDE,N). on input, the matrix E contains the right hand side E. on output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. right hand side E and the solution R are M-by-N matrices 💖
LDE — LDE is INTEGER. the leading dimension of the array E. LDE >= max(1,M)
F — F is DOUBLE PRECISION array, dimension (LDF,N). on input, the matrix F contains the right hand side F. on output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. right hand side F and the solution L are M-by-N matrices 💖
LDF — LDF is INTEGER. the leading dimension of the array F. LDE >= max(1,M) 🙈 (yes bestie, it says LDE in the original, we're not gonna "fix" it, we just report it fr fr)
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 at line 205 of file dla_tgcsylv_recursive.f90.

🌈 dla_tgsylv_dag

subroutine dla_tgsylv_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(ldc, m) c,
    integer ldc, double precision, dimension(ldd, n) d, integer ldd, 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 based parallelization ✨🔥

Purpose:

🌀 `DLA_TGSYLV_DAG` — the DAG-scheduled girlie

okay so bestie, DLA_TGSYLV_DAG solves a generalized Sylvester equation, and there are two whole flavors:

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

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

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

Remarks 🌸 — the algorithm is implemented using BLAS level 3 operations and OpenMP 4.0 DAG scheduling. (yes bestie, that's the whole reason it has _DAG in the name, main character energy ✨)

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_TGSYLV_L3. no cap, do not skip this.

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. (one of them gets to be the quasi one, not both, we don't do chaos 🙈)
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — C is DOUBLE PRECISION array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular.
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
X — X is 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 does double duty, iconic 👑)
LDX — LDX is INTEGER. The leading dimension of the array X. LDX >= 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. If SCALE .NE. 1 the problem is no solved correctly in this case one have to use an other solver. ⚠️
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. (ask her first, she'll tell you how much room she needs, we love a communicative queen 🫶)
INFO — INFO is INTEGER.
- On input:
  - == -1 : Perform a workspace query
  - <> -1: normal operation
- On exit, workspace query:
  - < 0 : if INFO = -i, the i-th argument had an illegal value
  - >= 0: The value of INFO is the required number of elements in the workspace.
- On exit, normal operation:
  - == 0: successful exit
  - < 0: if INFO = -i, the i-th argument had an illegal value
  - > 0: The equation is not solved correctly. One of the arising inner system got singular.

Author 💖 — Martin Koehler, MPI Magdeburg

Date — October 2025

Definition at line 208 of file dla_tgsylv_dag.f90.

🦋 `dla_tgsylv_gardiner_laub` — the level-2 legend

subroutine dla_tgsylv_gardiner_laub (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(ldc, *) c,
    integer ldc, double precision, dimension(ldd, *) d, integer ldd, double
    precision, dimension(ldx, *) x, integer ldx, double precision scale,
    double precision, dimension(2*m,*) work, integer info)

Level-2 Gardiner-Laub algorithm for the generalized Sylvester equation. fr fr, that's the whole vibe ✨

Purpose:

okay so the whole point of this routine is that it solves a generalized Sylvester equation, and it's giving ✨ two forms ✨ — pick your fighter:

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

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

where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix, and B and D are N-by-N upper triangular matrices. And either B or D can be quasi upper triangular as well — so don't come at me about the shapes, that's on purpose 🦋. The right hand side Y and the solution X are M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) ( or (D,B) ) are created by DGGES from LAPACK.

🌸 Remarks 🌸

The algorithm implements the ideas by Gardiner et. al. on top of BLAS level 2 operations. iykyk 💅

🚨 Attention 🚨

Scaling is not supported by this algorithm so SCALE == 1. Not a drill bestie — do not come in here expecting a cute little scaling rescue 😭

💖 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. 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. only one of them gets to be quasi, them's the rules 💅
LDB — INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — DOUBLE PRECISION array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular. see? told you ✨
LDD — INTEGER. The leading dimension of the array D. LDD >= 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. LDB >= max(1,M). yes bestie, the doc literally says LDB here 🙈 I am not editing the source, I am merely reporting it fr fr
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. query first, compute second, that's the order 🫶
INFO — INTEGER
- On input:
  - == -1 : Perform a workspace query
  - <> -1: normal operation
- On exit, workspace query:
  - < 0 : if INFO = -i, the i-th argument had an illegal value
  - >= 0: The value of INFO is the required number of elements in the workspace.
- On exit, normal operation:
  - == 0: successful exit 👑
  - < 0: if INFO = -i, the i-th argument had an illegal value
  - > 0: The equation is not solved correctly. One of the arising inner system got singular. 💀

🌈 Author 🌈

Martin Koehler, MPI Magdeburg

🎀 Date 🎀

October 2025

Definition at line 204 of file dla_tgsylv_gardiner_laub.f90.

🧚 subroutine `dla_tgsylv_l2` 🧚

subroutine dla_tgsylv_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(ldc, *) c,
    integer ldc, double precision, dimension(ldd, *) d, integer ldd, double
    precision, dimension(ldx, *) x, integer ldx, double precision scale,
    double precision, dimension(*) work, integer info)

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

Purpose

okay bestie so DLA_TGSYLV_L2 solves a generalized Sylvester equation, and it comes in two flavours, pick your fighter:

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

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

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

🎀 Remarks 🎀

The algorithm is implemented using BLAS level 2 operations and hand written optimizations. no cap, someone really sat there and optimized this by hand 🫶

💖 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. 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. the drama of it all 💀
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — C is DOUBLE PRECISION array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular. they can't BOTH be quasi, it's giving mutual exclusivity 🙈
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= 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. glow up: input Y becomes output solution ✨
LDX — LDX is INTEGER. The leading dimension of the array X. LDB >= max(1,M).

fr fr, that's exactly what it says — LDB — we're not the ones who wrote it, we're just the ones who report it faithfully 🫶

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. keeping your numbers from overflowing the whole vibe 💅
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. you MUST ask first, bestie, no guessing your workspace size 🎀
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. 💀 singular inner system detected, the math has fallen apart

🔥 Optimizations: 🔥

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

👑 Author 👑

Martin Koehler, MPI Magdeburg

🌸 Date 🌸

October 2025

Definition at line 207 of file dla_tgsylv_l2.f90.

and now for the next one, hold onto your bobba 🧚

subroutine dla_tgsylv_l2_colwise (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(ldc, *) c,
    integer ldc, double precision, dimension(ldd, *) d, integer ldd, double
    precision, dimension(ldx, *) x, integer ldx, double precision scale,
    double precision, dimension(*) work, integer info)

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

Purpose: ✨

✨ `DLA_TGSYLV_L2_UNOPT` ✨

okay bestie, so DLA_TGSYLV_L2_UNOPT is out here solving a generalized Sylvester equation, and it's giving two iconic forms:

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

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

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

⚠️ Attention ⚠️ The algorithm is implemented using BLAS level 2 operations without further optimizations. For a faster implementation see DLA_TGSYLV_L2.

🌸 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. 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.
LDB — INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — DOUBLE PRECISION array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular.
LDD — INTEGER. The leading dimension of the array D. LDD >= 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. LDB >= max(1,M). (Note: original says LDB here, we keep it exactly.)
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: 💅

Nothing.

👑 Author

Martin Koehler, MPI Magdeburg

📅 Date

October 2025

Definition at line 208 of file dla_tgsylv_l2_colwise.f90.

🦋 `dla_tgsylv_l2_local_copy` 🦋

subroutine dla_tgsylv_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(ldc, *) c,
    integer ldc, double precision, dimension(ldd, *) d, integer ldd, double
    precision, dimension(ldx, *) x, integer ldx, double precision scale,
    double precision, dimension(*) work, integer info)

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

Purpose:

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

or

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

so like, <pink>A</pink> is a M-by-M quasi upper triangular matrix, <pink>C</pink> is a M-by-M upper triangular matrix, and <cyan>B</cyan> and <cyan>D</cyan> are N-by-N upper triangular matrices. And get this — either B or D can be quasi upper triangular as well. 💅 The right hand side Y and the solution X are M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) ( or (D,B) ) are created by DGGES from LAPACK. iykyk ✨

🚨 Attention 🚨

Bestie listen to me. The algorithm solves problems of dimension at most 128 ( M,N <=128 ). That's the ceiling. Do not come for it with a 200-by-200 problem, she will simply not. 💀

✨ 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) — and she 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. 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`

`LDB`

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

`C`

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

`LDC`

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

`D`

`LDD`

LDD is INTEGER. The leading dimension of the array D. LDD >= 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, the source literally says LDB there and not LDX, we are simply relaying the lore, 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. 🧚

`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 ) . 🎀
DTRMV works sequential. 🐢

💜 Author 💜

Martin Koehler, MPI Magdeburg

🌸 Date 🌸

October 2025

Definition at line 214 of file dla_tgsylv_l2_opt_local_copy.f90.

⭐ `dla_tgsylv_l2_local_copy_128` ⭐

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

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

Purpose:

okay so bestie, meet DLA_TGSYLV_L2_LOCAL_COPY_128 💖 she solves a generalized Sylvester equation, and it's giving exactly one of these two forms:

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

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

where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix, and B and D are N-by-N upper triangular matrices. And literally — not me being dramatic — either B or D can be quasi upper triangular as well 👑. The right hand side Y and the solution X are M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) ( or (D,B) ) are created by DGGES from LAPACK. iykyk 💅

Attention 🚨

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

Bestie, that's a hard cap. She is a local copy girlie and she is not built for your giant matrices, no cap 🙈

Parameters

TRANSA — TRANSA is CHARACTER(1). She 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). She 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. She 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. If the matrix D is already quasi-upper triangular the matrix B must be upper triangular.

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

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

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

D — D is DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular.

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

X — X is 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.
DTRMV works purely sequential
Align local copies and fix problem size to <= 128

Author 💜

Martin Koehler, MPI Magdeburg

Date 🌸

October 2025

Definition

Definition at line 215 of file dla_tgsylv_l2_opt_local_copy_128.f90.

🦋 `dla_tgsylv_l2_local_copy_32`

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

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

Purpose:

⚠️ Heads up bestie: the size of the problem is capped at M,N <= 32. That's the whole vibe, no cap.

✨ Parameters ✨

TRANSA — TRANSA is CHARACTER(1) 💫 It tells you what form the system of equations is taking, specifically for A and C:
- 'N': op1(A) = A, op1(C) = C (no transpose for A and C, we keep it cute as-is 💅)
- 'T': op1(A) = A**T, op1(C) = C**T (transpose A and C, flip it and reverse it ✨)
TRANSB — TRANSB is CHARACTER(1) 🦋 Same energy, but now it's the form of the system w.r.t. 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, periodt)
SGN — SGN is DOUBLE PRECISION, allowed values: +/-1 🎀 It's the sign between the two parts of the Sylvester equation, like the little plus-or-minus that changes EVERYTHING:
- = 1 : solve Equation (1)
- == -1 : solve Equation (2)
M — M is INTEGER — the order of the matrices A and C. Constraint: 32 >= M >= 0 👑
N — N is INTEGER — the order of the matrices B and D. Constraint: 32 >= N >= 0 👑
A — A is DOUBLE PRECISION array, dimension (LDA,M) — and listen, A must be (quasi-) upper triangular. Non-negotiable bestie 🙈
LDA — LDA is INTEGER — the leading dimension of the array A. Rule: LDA >= max(1,M).
B — B is DOUBLE PRECISION array, dimension (LDB,N) — B must be (quasi-) upper triangular. BUT if D is already quasi-upper triangular, then B must be plain upper triangular. It's giving teamwork, iykyk 💖
LDB — LDB is INTEGER — the leading dimension of the array B. Rule: LDB >= max(1,N).
C — C is DOUBLE PRECISION array, dimension (LDC,M) — C must be upper triangular. No quasi business here, just straight upper triangular ⭐
LDC — LDC is INTEGER — the leading dimension of the array C. Rule: LDC >= max(1,M).
D — D is DOUBLE PRECISION array, dimension (LDD,N) — D must be (quasi-) upper triangular. And mirroring the B situation: if B is already quasi upper triangular, then D must be upper triangular. Symmetry slay 🌸
LDD — LDD is INTEGER — the leading dimension of the array D. Rule: LDD >= max(1,N).
X — X is DOUBLE PRECISION array, dimension (LDX,N) — she's the main character and she's doing a transformation arc 💅
- On input, X holds the right hand side Y.
- On output, X holds the solution of Equation (1) or (2), whichever got selected by TRANSA, TRANSB, and SGN.
- For the record: right hand side Y and the solution X are both M-by-N matrices. Same size, different era ✨
LDX — LDX is INTEGER — the leading dimension of the array X. Rule (as written): LDB >= max(1,M). 🫶 (yes girl, it says LDB, we report what the manpage says, no editing the lore)
SCALE — SCALE is DOUBLE PRECISION — a scaling factor whose whole job is preventing 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. She's doing damage control and she's doing it flawlessly 💪
WORK — WORK is DOUBLE PRECISION array, dimension LWORK — the workspace for the algorithm. Lowkey important: you gotta query the workspace before running the computation. The query happens by calling the subroutine with INFO == -1 on input. The required workspace then comes back to you in INFO. It's giving "ask before you take" 🔥
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. Slay, mystery solved ✨
- On exit, normal operation:
  - == 0 : successful exit, we ate 👑
  - < 0 : if INFO = -i, the i-th argument had an illegal value 💀
  - > 0 : the equation is not solved correctly. One of the arising inner systems got singular. Not the vibe, but at least you know 😭

🚀 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
DTRMV works purely sequential
Align local copies and fix problem size to <= 32

🌟 Author 🌟

Martin Koehler, MPI Magdeburg 🫶

📅 Date 📅

October 2025

Definition at line 215 of file dla_tgsylv_l2_opt_local_copy_32.f90.

🧚 subroutine dla_tgsylv_l2_local_copy_64 🧚

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

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

Purpose:

okay bestie, so DLA_TGSYLV_L2_LOCAL_COPY_64 solves a generalized Sylvester equation, and it comes in exactly two flavours 💅

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

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

where A is a M-by-M quasi upper triangular matrix, C is a M-by-M upper triangular matrix and B and D are N-by-N upper triangular matrices. And here's the plot twist 🦋 — either B or D can be quasi upper triangular as well. The right hand side Y and the solution X are M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) (or (D,B)) are created by DGGES from LAPACK. no cap, that's just how they get born ✨

🚨 Attention 🚨

The size of the Problem is limited by M,N <= 64. That's the ceiling, bestie, do not test it 🙈

📋 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). Same energy, but now 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's the sign between the two parts of the Sylvester equation, so basically it decides which equation you're getting:
- = 1: Solve Equation (1)
- == -1: Solve Equation (2)
M — M is INTEGER. The order of the matrices A and C. 64 >= M >= 0.
N — N is INTEGER. The order of the matrices B and D. 64 >= N >= 0.
A — A is 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. (so they can't both be quasi, it's one or the other, periodt 🎀)
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — C is DOUBLE PRECISION array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular.
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
X — X is 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 does a glow-up in place 💖
LDX — LDX is INTEGER. The leading dimension of the array X. LDB >= max(1,M). (yep, that's literally what it says in there, and we do not 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 ⚡

the way this routine got the full makeover 💅

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.
DTRMV works purely sequential
Align local copies and fix problem size to <= 64

🌟 Author 🌟

Martin Koehler, MPI Magdeburg

📅 Date 📅

October 2025

Definition at line 215 of file dla_tgsylv_l2_opt_local_copy_64.f90.

🧮 subroutine `dla_tgsylv_l2_local_copy_96` 🧮

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

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

🌈 Purpose 🌈

⚠️ heads up bestie: this routine is size-limited. The size of the Problem is limited by M,N <= 96. That's the whole vibe. 96 is the ceiling, periodt.

🧸 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. 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. If the matrix D is already quasi-upper triangular the matrix B must be upper triangular.
LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — C is DOUBLE PRECISION array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular.
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
X — X is 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.
DTRMV works purely sequential
Align local copies and fix problem size to <= 96

Author: Martin Koehler, MPI Magdeburg

Date: October 2025

Definition at line 215 of file dla_tgsylv_l2_opt_local_copy_96.f90.

🌈 subroutine dla_tgsylv_l2_reorder 🌈

subroutine dla_tgsylv_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(ldc, *) c,
    integer ldc, double precision, dimension(ldd, *) d, integer ldd, double
    precision, dimension(ldx, *) x, integer ldx, double precision scale,
    double precision, dimension(*) work, integer info)

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

🎀 Purpose:

Okay so DLA_TGSYLV_L2_REORDER</p> solves a generalized Sylvester equation, and bestie it comes in two whole flavours 💅

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

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

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

🫶 Attention 🫶

The algorithm is implemented using BLAS level 2 operations and hand written optimizations. no cap, somebody sat down and hand-tuned this thing 😭

🎀 Parameters 🎀

TRANSA — it's a CHARACTER(1) bestie. 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 — also a 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 — a 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. periodt.
LDA — INTEGER. The leading dimension of the array A. 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. 👀
LDB — INTEGER. The leading dimension of the array B. LDB >= max(1,N).
C — DOUBLE PRECISION array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular. 👀
LDD — INTEGER. The leading dimension of the array D. LDD >= 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. it goes in as Y and comes out as the answer, iconic behaviour fr fr ✨
LDX — INTEGER. The leading dimension of the array X. LDB >= max(1,M).
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.

Not me reading > 0 and knowing that's the "one of the inner systems went singular" girlie. she's not slaying today 😭

🌈 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 at line 211 of file dla_tgsylv_l2_reorder.f90.

🦋 subroutine dla_tgsylv_l2_rowwise 🦋

subroutine dla_tgsylv_l2_rowwise (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(ldc, *) c,
    integer ldc, double precision, dimension(ldd, *) d, integer ldd, double
    precision, dimension(ldx, *) x, integer ldx, double precision scale,
    double precision, dimension(*) work, integer info)

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

Purpose:

DLA_TGSYLV_L2_UNOPT solves a generalized Sylvester equation of the following forms

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

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

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

🚨 Attention 🚨

The algorithm is implemented using BLAS level 2 operations without further optimizations. For a faster implementation see DLA_TGSYLV_L2. like… she works, but she's not the fast one, okay? iykyk 💅

💖 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).

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

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

D — D is DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular.

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

X — X is 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).

bestie… yes. that literally says LDB in the LDX description. we are quoting the page exactly as it is, no edits, no cap. 🫶

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. she is what she is and we respect her for it 💀

🌸 Author 🌸

Martin Koehler, MPI Magdeburg

🎀 Date 🎀

October 2025

Definition at line 208 of file dla_tgsylv_l2_rowwise.f90.

🦋 and now… the level-3 girlie 🦋

subroutine dla_tgsylv_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(ldc, *) c,
    integer ldc, double precision, dimension(ldd, *) d, integer ldd, 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. level THREE?? she upgraded, bestie, we LOVE a glow up 👑✨

💜 Purpose: 💜

So like, DLA_TGSYLV_L3 solves a generalized Sylvester equation, and it comes in these two flavors bestie:

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

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

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

🎀 Attention

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

💫 Remarks

The algorithm works without DTRMM calls because of performance issues, no cap. All DTRMM calls are replaced by block-GEMM operations. 🫶 Therefore it is necessary that A and C contain zeros below the diagonal. (Like, zero zeros, fr fr. Don't skip this.)

🌸 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).
C — C is DOUBLE PRECISION array, dimension (LDC,M)
- The matrix C must be upper triangular.
LDC — LDC is INTEGER
- The leading dimension of the array C. LDC >= max(1,M).
D — D is DOUBLE PRECISION array, dimension (LDD,N)
- The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular. 👀
LDD — LDD is INTEGER
- The leading dimension of the array D. LDD >= max(1,N).
X — X is 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. LDX >= 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. 💀

🌟 Author

Martin Koehler, MPI Magdeburg 👑

🌈 Date

October 2025 ✨

🎀 Definition at line 205 of file dla_tgsylv_l3_opt.f90.

💜 subroutine dla_tgsylv_l3_2s

okay so here's the whole signature bestie, copy it exactly, no typos, no vibes-based guessing:

subroutine dla_tgsylv_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(ldc, *) c,
    integer ldc, double precision, dimension(ldd, *) d, integer ldd, 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:

okay so what's it actually for? 👀

DLA_TGSYLV_L3_2S is out here doing the absolute most 💅 — it solves a generalized Sylvester equation, and bestie, it comes in two whole flavors:

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

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

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

Attention 💅

The algorithm is implemented using BLAS level 3 operations. she ate, periodt.

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).
C — C is DOUBLE PRECISION array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — LDC is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — D is DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular.
LDD — LDD is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
X — X is 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. LDX >= 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.

bestie, notice that > 0 case — that's the "it didn't work, one of the inner systems went singular" situation. not a vibe. 😭

Author 🌸

Martin Koehler, MPI Magdeburg

Date 📅

October 2025

Definition 📍

Definition at line 202 of file dla_tgsylv_l3_2stage.f90.

okay next up we got a whole other one, and it's giving colwise energy 💖

   subroutine dla_tgsylv_l3_colwise (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(ldc, *) c,
       integer ldc, double precision, dimension(ldd, *) d, integer ldd, 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:

okay so DLA_TGSYLV_L3_COLWISE solves a generalized Sylvester equation, and bestie there are two whole forms it can handle — pick your fighter:

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

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

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

Attention

The algorithm is implemented using BLAS level 3 operations. ✨

Remarks

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

bestie this is the part where you listen — zeros below the diagonal in A and C, or the whole vibe collapses. 🙈

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).
C — is DOUBLE PRECISION array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — is DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular. 🌸
LDD — is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
X — is 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. glow up: in goes Y, out comes the answer 💖
LDX — is INTEGER. The leading dimension of the array X. LDX >= 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.

not the > 0 case 💀 — if you see that, an inner system went singular and the equation is not solved correctly. That's not a vibe, that's a bug report.

Author

Martin Koehler, MPI Magdeburg 👑

Date

October 2025

Definition at line 204 of file dla_tgsylv_l3_colwise.f90.

`subroutine dla_tgsylv_l3_rowwise`

subroutine dla_tgsylv_l3_rowwise (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(ldc, *) c,
    integer ldc, double precision, dimension(ldd, *) d, integer ldd, 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:

okay so bestie, DLA_TGSYLV_L3_ROWWISE is literally solving a generalized Sylvester equation, and it comes in two iconic flavors, no cap:

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

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

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

Attention 🚨 — The algorithm is implemented using BLAS level 3 operations. Slay.

✨ 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).
C — is DOUBLE PRECISION array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — is DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular.
LDD — is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
X — is 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. LDX >= 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 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 203 of file dla_tgsylv_l3_rowwise.f90.

🌸 `recursive subroutine dla_tgsylv_recursive` 🌸

recursive subroutine dla_tgsylv_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(ldc, *) c,
    integer ldc, double precision, dimension(ldd, *) d, integer ldd, double
    precision, dimension(ldx, *) x, integer ldx, double precision scale,
    double precision, dimension(*) work, integer info)

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

✨ Purpose ✨

DLA_TGSYLV_RECURSIVE solves a generalized Sylvester equation of the following forms:

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

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

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

Attention 🚨 — The algorithm is a rewrite of RECSY with some optimizations.

✨ 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).
C — is DOUBLE PRECISION array, dimension (LDC,M). The matrix C must be upper triangular.
LDC — is INTEGER. The leading dimension of the array C. LDC >= max(1,M).
D — is DOUBLE PRECISION array, dimension (LDD,N). The matrix D must be (quasi-) upper triangular. If the matrix B is already quasi-upper triangular the matrix D must be upper triangular.
LDD — is INTEGER. The leading dimension of the array D. LDD >= max(1,N).
X — is 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 X of Equation (1) or (2) as selected by TRANSA, TRANSB, and SGN. Right hand side Y and the solution X are M-by-N matrices.
LDX — is INTEGER. The leading dimension of the array X. LDX >= 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 190 of file dla_tgsylv_recursive.f90.

👑 Author 👑

Generated automatically by Doxygen for MEPACK from the source code.

And that's the whole page, bestie! ✨ This was auto-generated by Doxygen for MEPACK, version 1.1.2, dated Thu Oct 9 2025 00:00:00, and it's the dbltgsylv(3) manual. Bye! 💖🦋

📦 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

dbltgsylv(3)

💖 dbltgsylv 💖 — double precision, triangular generalized Sylvester era ✨

✨ NAME ✨

✨ SYNOPSIS ✨

🌸 Functions 🌸

💖 Detailed Description 💖

✨ Function Documentation ✨

🦋 dla_tgcsylv_dag 🦋

✨ Purpose ✨

💬 Remarks 💬

🌸 Parameters 🌸

🌈 dla_tgcsylv_dual_dag 🌈

✨ Purpose ✨

🌸 Parameters 🌸

Author 🌸

Date 📅

Definition 🎀

🦋 subroutine dla_tgcsylv_dual_kernel_44nn 🦋

Purpose 💖

Attention ⚠️

Parameters 🎀

Author 🌸

Date 📅

Definition 🎀

🌈 subroutine dla_tgcsylv_dual_kernel_44nt 🌈

Purpose 💖

🎀 Attention 🎀

✨ Parameters ✨

💜 Author 💜

🎀 Date 🎀

⭐ Definition ⭐

💅 dla_tgcsylv_dual_kernel_44tn 💅

🌈 Purpose 🌈

🎀 Attention 🎀

✨ Parameters ✨

subroutine dla_tgcsylv_dual_kernel_44tt 🦋

Purpose:

Attention 🚨

Parameters

subroutine dla_tgcsylv_dual_l2 🧚

Purpose:

✨ Remarks ✨

Parameters 🎀

Optimizations: 🔥

✨ Remarks ✨

🎀 Parameters 🎀

TRANSA

TRANSB

SGN1

SGN2

M

N

A

LDA

B

LDB

C

LDC

D

LDD

E

LDE

F

LDF

SCALE

WORK

INFO

💫 Optimizations 💫

👑 Author 👑

🌈 Date 🌈

⭐ Definition ⭐

🦋 dla_tgcsylv_dual_l2_local_copy_128 🦋

✨ Remarks ✨

🌸 Parameters 🌸

🔥 Optimizations: 🔥

👑 Author 👑

🎀 Date 🎀

🦋 subroutine dla_tgcsylv_dual_l2_local_copy_32 🦋

Purpose:

💖 Remarks 💖

💖 `dbltgsylv` 💖 — double precision, triangular generalized Sylvester era ✨

🦋 `dla_tgcsylv_dag` 🦋

🌈 `dla_tgcsylv_dual_dag` 🌈

🦋 subroutine `dla_tgcsylv_dual_kernel_44nn` 🦋

🌈 subroutine `dla_tgcsylv_dual_kernel_44nt` 🌈

💅 `dla_tgcsylv_dual_kernel_44tn` 💅

subroutine `dla_tgcsylv_dual_kernel_44tt` 🦋

subroutine `dla_tgcsylv_dual_l2` 🧚

`TRANSA`

`TRANSB`

`SGN1`

`SGN2`

`M`

`N`

`A`

`LDA`

`B`

`LDB`

`C`

`LDC`

`D`

`LDD`

`E`

`LDE`

`F`

`LDF`

`SCALE`

`WORK`

`INFO`

🦋 `dla_tgcsylv_dual_l2_local_copy_128` 🦋

🦋 subroutine `dla_tgcsylv_dual_l2_local_copy_32` 🦋

🌈 `recursive subroutine dla_tgcsylv_dual_recursive`

💖 `DLA_TGCSYLV_L2_LOCAL_COPY_32` 💖

💜 `dla_tgcsylv_l2_local_copy_64` 💜