dbltglyap(3)

📦

version 1.1.2 · 3 versions

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

💅 dbltglyap — Double Precision 💅

okay bestie, gather round, because this is dbltglyap(3) from MEPACK and it is a lot in the best way. It's the whole double precision family for generalized Lyapunov and Stein equations, but specifically the ones where your coefficient matrices are triangular. 👑

✨ NAME ✨

dbltglyap — Double Precision

Double precision solvers for generalized Lyapunov and Stein equations with triangular coefficient matrices. that's the whole vibe. 💖

🌸 SYNOPSIS 🌸

🎀 Functions

okay so here's the full cast, no one gets left out:

dla_tglyap_dag — Level-3 Bartels-Stewart Algorithm for the generalized Lyapunov equation with DAG scheduling.

subroutine dla_tglyap_dag (trans, m, a, lda, b, ldb, x, ldx, scale,
    work, info)

dla_tgstein_dag — Level-3 Bartels-Stewart Algorithm for the generalized Stein equation with DAG scheduling.

subroutine dla_tgstein_dag (trans, m, a, lda, b, ldb, x, ldx, scale,
    work, info)

dla_tglyap_l2 — Level-2 Bartels-Stewart Algorithm for the generalized Lyapunov equation.

subroutine dla_tglyap_l2 (trans, m, a, lda, b, ldb, x, ldx, scale, work,
    info)

dla_tglyap_l2_unopt — Level-2 Bartels-Stewart Algorithm for the generalized Lyapunov equation (Unoptimized)

subroutine dla_tglyap_l2_unopt (trans, m, a, lda, b, ldb, x, ldx, scale,
    work, info)

dla_tgstein_l2 — Level-2 Bartels-Stewart Algorithm for the generalized Stein equation.

subroutine dla_tgstein_l2 (trans, m, a, lda, b, ldb, x, ldx, scale,
    work, info)

dla_tglyap_l3 — Level-3 Bartels-Stewart Algorithm for the generalized Lyapunov equation.

subroutine dla_tglyap_l3 (trans, m, a, lda, b, ldb, x, ldx, scale, work,
    info)

dla_tglyap_l3_2s — Level-3 Bartels-Stewart Algorithm for the generalized Lyapunov equation two stage.

subroutine dla_tglyap_l3_2s (trans, m, a, lda, b, ldb, x, ldx, scale,
    work, info)

dla_tgstein_l3 — Level-3 Bartels-Stewart Algorithm for the generalized Stein equation.

subroutine dla_tgstein_l3 (trans, m, a, lda, b, ldb, x, ldx, scale, work,
    info)

dla_tgstein_l3_2s — Level-3 Bartels-Stewart Algorithm for the generalized Stein equation (Two Stage)

subroutine dla_tgstein_l3_2s (trans, m, a, lda, b, ldb, x, ldx, scale,
    work, info)

dla_tglyap_recursive — Recursive Blocking Algorithm for the generalized Lyapunov equation.

recursive subroutine dla_tglyap_recursive (trans, m, a, lda, b, ldb, x,
    ldx, scale, work, info)

dla_tgstein_recursive — Recursive Blocking Algorithm for the generalized Stein equation.

recursive subroutine dla_tgstein_recursive (trans, m, a, lda, b, ldb, x,
    ldx, scale, work, info)

dla_ggstein_refine — Iterative Refinement for the Generalized Lyapunov Equations.

subroutine dla_ggstein_refine (trans, guess, m, a, lda, b, ldb, x, ldx,
    y, ldy, as, ldas, bs, ldbs, q, ldq, z, ldz, maxit, tau, convlog,
    work, ldwork, info)

🦋 DETAILED DESCRIPTION 🦋

Double precision solvers for generalized Lyapunov and Stein equations with triangular coefficient matrices. that's the entire thesis statement, and honestly? iconic. ✨

📖 FUNCTION DOCUMENTATION 📖

🌈 `dla_ggstein_refine`

Iterative Refinement for the Generalized Lyapunov Equations.

subroutine dla_ggstein_refine (character, dimension(1) trans, character,
    dimension(1) guess, integer m, double precision, dimension(lda, *) a,
    integer lda, double precision, dimension(ldb, *) b, integer ldb, double
    precision, dimension( ldx, *) x, integer ldx, double precision,
    dimension(ldy,*) y, integer ldy, double precision, dimension(ldas, *)
    as, integer ldas, double precision, dimension(ldbs,*) bs, integer ldbs,
    double precision, dimension(ldq, *) q, integer ldq, double precision,
    dimension(ldz, *) z, integer ldz, integer maxit, double precision tau,
    double precision, dimension(*) convlog, double precision, dimension(*)
    work, integer ldwork, integer info)

💫 Purpose

so DLA_GGSTEIN solves a generalized Stein equation, and it comes in exactly two flavours bestie:

   A * X * A^T - B * X * B^T = SCALE * Y                                              (1)

   A^T * X * A - B^T * X * B =  SCALE * Y                                             (2)

where (A,B) is a M-by-M matrix pencil, and we're solving it with iterative refinement. the right hand side Y and the solution X are both M-by-M matrices. 🫶

and here's the thing — the matrix pencil (A,B) has to hand over both its original data and its generalized Schur decomposition, because the iterative refinement process needs both. no shortcuts, both or nothing. 💅

🔥 Parameters

TRANS — TRANS is CHARACTER. it picks which form of the system we're solving, with respect to A:
- == 'N': Equation (1) is solved
- == 'T': Equation (2) is solved
GUESS — GUESS is CHARACTER. tells us whether X is walking in with an initial guess or not:
- = 'I': X contains an initial guess
- = 'N': No initial guess — X gets set to zero at the beginning of the iteration
M — M is INTEGER. the order of the matrices A and B. and it's M >= 0, periodt.
A — A is DOUBLE PRECISION array, dimension (LDA,M). this array holds the original matrix A that defines the equation. 💖
LDA — LDA is INTEGER. the leading dimension of the array A, and LDA >= max(1,M).
B — B is DOUBLE PRECISION array, dimension (LDB,M). this one holds the original matrix B that defines the equation.
LDB — LDB is INTEGER. leading dimension of the array B, with LDB >= max(1,M).
X — X is DOUBLE PRECISION array, dimension (LDX,M). on input, it holds the initial guess — but only if GUESS = 'I'. on output, it holds the solution X. glow up. ✨
LDX — LDX is INTEGER. leading dimension of the array X, LDX >= max(1,M).
Y — Y is DOUBLE PRECISION array, dimension (LDY,M). on input it holds the right hand side Y, and it stays unchanged for the entire iteration. it's just built different, it doesn't move. 🙈
LDY — LDY is INTEGER. leading dimension of the array Y, LDY >= max(1,M).
AS — AS is DOUBLE PRECISION array, dimension (LDAS,M). holds the generalized Schur decomposition of the A.
LDAS — LDAS is INTEGER. leading dimension of the array AS, LDAS >= max(1,M).
BS — BS is DOUBLE PRECISION array, dimension (LDBS,M). holds the generalized Schur decomposition of the B.
LDBS — LDBS is INTEGER. leading dimension of the array BS, LDBS >= max(1,M).
Q — Q is DOUBLE PRECISION array, dimension (LDQ,M). holds the left generalized Schur vectors for (A,B) as returned by DGGES. 👑
LDQ — LDQ is INTEGER. leading dimension of the array Q, LDQ >= max(1,M).
Z — Z is DOUBLE PRECISION array, dimension (LDZ,M). holds the right generalized Schur vectors for (A,B) as returned by DGGES.
LDZ — LDZ is INTEGER. leading dimension of the array Z, LDZ >= max(1,M).
MAXIT — MAXIT is INTEGER. on input it's the maximum number of iterations that get performed, and MAXIT <= 100. on exit it holds the number of iteration steps the algorithm actually took. she reports back. 💅
TAU — TAU is DOUBLE PRECISION. on input, the additional security factor for the stopping criterion — typical values are 0.1. on exit, it holds the last relative residual from when the stopping criterion kicked in.
CONVLOG — CONVLOG is DOUBLE PRECISION array, dimension (MAXIT). the convergence history of the iterative refinement lives here. CONVLOG(I) holds the maximum relative residual before it is solved for the I-th time. receipts, basically. 📝
WORK — WORK is DOUBLE PRECISION array, dimension (MAX(1,LDWORK)). workspace for the algorithm. the optimal workspace gets returned in LDWORK if LDWORK == -1 on input — and in that case, no computations are performed at all.
LDWORK — LDWORK is INTEGER. if LDWORK == -1, the subroutine returns the required size of the workspace in LDWORK on exit. no computations are performed and none of the arrays are referenced. it's just vibes and math. ✨
INFO — INFO is INTEGER:
- == 0: Success 🎉
- > 0: Iteration failed in step INFO 😭
- < 0: if INFO = -i, the i-th argument had an illegal value
- = -50: Some of the internal settings like NB,... are incorrect. 💀

🧚 See also

DLA_TGSTEIN_L3
DLA_TGSTEIN_L2
DLA_TGSTEIN_L3_2S
DLA_TGSTEIN_DAG
DLA_TGSTEIN_RECURSIVE

🌟 Author

Martin Koehler, MPI Magdeburg

🌙 Date

October 2025

📌 Definition

Definition at line 243 of file dla_ggstein_refine.f90.

💜 `dla_tglyap_dag`

Level-3 Bartels-Stewart Algorithm for the generalized Lyapunov equation with DAG scheduling.

subroutine dla_tglyap_dag (character, dimension(1) trans, integer m, double
    precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldb,*) b, integer ldb, double precision, dimension(ldx, *) x,
    integer ldx, double precision scale, double precision, dimension(*)
    work, integer info)

💫 Purpose

DLA_TGLYAP_DAG solves a generalized Lyapunov equation, and yeah, it also comes in two forms:

   A * X * B^T + A * X * B^T = SCALE * Y                                              (2)

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

where A is a M-by-M quasi upper triangular matrix, B is a M-by-M upper triangular matrix, and X and Y are symmetric M-by-M matrices. 🌸

typically the matrix pencil (A,B) gets created by DGGES from LAPACK. iykyk.

⚠️ Attention

The algorithm is implemented using BLAS level 3 operations and DAG. that's the whole flex.

🔥 Parameters

TRANS — TRANS is CHARACTER. specifies the form of the system of equations with respect to A and B:
- == 'N': Equation (1) is solved.
- == 'T': Equation (2) is solved.
M — M is INTEGER. the order of the matrices A and B. M >= 0.
A
A — DOUBLE PRECISION array, dimension (LDA,M). Bestie, this matrix has to be (quasi-) upper triangular. That's not a suggestion, that's the vibe check it has to pass. 💅
LDA — INTEGER. The leading dimension of the array A. Rule of the land: LDA >= max(1,M). 🎀
B — DOUBLE PRECISION array, dimension (LDB,M). The matrix B must be upper triangular. Note it's upper triangular proper, no quasi-ness here, she's a strict queen. 👑
LDB — INTEGER. The leading dimension of the array B. And yes: LDB >= max(1,M). 🦋
X — DOUBLE PRECISION array, dimension (LDX,N). On input, X holds the right hand side Y. On output, X holds the solution of Equation (1) or (2), whichever one you picked with TRANSA, TRANSB, and SGN. Both the right hand side Y and the solution X are M-by-N matrices. So she's the glow-up container: comes in as a problem, leaves as the answer. ✨
LDX — INTEGER. The leading dimension of the array X. Constraint, say it with me: LDX >= max(1,M). 💖
SCALE — DOUBLE PRECISION. This is the scaling factor whose whole job is stopping overflow from ruining the result. 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. So she's the safety net, no cap. 🌸
WORK — DOUBLE PRECISION array, dimension LWORK. Workspace for the algorithm. And here's the thing, bestie: you need to query the workspace before running the computation. You do that by calling the subroutine with INFO == -1 on input, and the required workspace comes back to you in INFO. Don't skip the query, that's the whole plot. 🙈
INFO — INTEGER. She's the receipts, so let's read her carefully.

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 154 of file dla_tglyap_dag.f90.

subroutine dla_tglyap_l2

subroutine dla_tglyap_l2 (character, dimension(1) trans, integer m, double
    precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, *) x,
    integer ldx, double precision scale, double precision, dimension(*)
    work, integer info)

Level-2 Bartels-Stewart Algorithm for the generalized Lyapunov equation.

Purpose:

fortran !> !> DLA_TGLYAP_L2 solves a generalized Lyapunov equation of the following forms !> !> A * X * B^T + B * X * A^T = SCALE * Y (1) !> !> or !> !> A^T * X * B + B^T * X * A = SCALE * Y (2) !> !> where A is a M-by-M quasi upper triangular matrix, B is a M-by-M upper triangular, !> and X and Y are symmetric M-by-M matrices. !> Typically the matrix pencil (A,B) is created by DGGES from LAPACK. !>

Attention

The algorithm is implemented using BLAS level 2 operations. That's the whole tea. 🍓

Parameters

TRANS — CHARACTER. Specifies the form of the system of equations with respect to A and B:
- == 'N' : Equation (1) is solved.
- == 'T' : Equation (2) is solved.
M — INTEGER. The order of the matrices A and B. Constraint: M >= 0. 🎀
A — DOUBLE PRECISION array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular. No exceptions, bestie. 💅
LDA — INTEGER. The leading dimension of the array A. LDA >= max(1,M). 🦋
B — DOUBLE PRECISION array, dimension (LDB,M). The matrix B must be upper triangular. 💖
LDB — INTEGER. The leading dimension of the array B. LDB >= max(1,M). 🌸
X — DOUBLE PRECISION array, dimension (LDX,N). On input, X contains the right hand side Y. On output, X contains the solution of Equation (1) or (2) as selected by TRANS. Right hand side Y and the solution X are M-by-N matrices. ✨
LDX — INTEGER. The leading dimension of the array X. LDX >= max(1,M). 💜
SCALE — DOUBLE PRECISION. 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. 🌈
WORK — DOUBLE PRECISION array, dimension LWORK. Workspace for the algorithm. The workspace needs to be queried before running the computation — you query by calling the subroutine with INFO == -1 on input, and 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 155 of file dla_tglyap_l2.f90.

subroutine dla_tglyap_l2_unopt

subroutine dla_tglyap_l2_unopt (character, dimension(1) trans, integer m,
    double precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, *) x,
    integer ldx, double precision scale, double precision, dimension(*)
    work, integer info)

Level-2 Bartels-Stewart Algorithm for the generalized Lyapunov equation (Unoptimized)

Purpose:

fortran !> !> DLA_TGLYAP_L2 solves a generalized Lyapunov equation of the following forms !> !> A * X * B^T + B * X * A^T = SCALE * Y (1) !> !> or !> !> A^T * X * B + B^T * X * A = SCALE * Y (2) !> !> where A is a M-by-M quasi upper triangular matrix, B is a M-by-M upper triangular, !> and X and Y are symmetric M-by-M matrices. !> Typically the matrix pencil (A,B) is created by DGGES from LAPACK. !>

Attention

The algorithm is implemented using BLAS level 2 operations but works with inefficient DSYR2 calls in the case TRANS="T". So like, she functions, but she is not thriving in that mode. 😭

Parameters

TRANS — CHARACTER. Specifies the form of the system of equations with respect to A and B:
- == 'N' : Equation (1) is solved.
- == 'T' : Equation (2) is solved.
M — INTEGER. The order of the matrices A and B. Constraint: M >= 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,M). The matrix B must be upper triangular. 💖
LDB — INTEGER. The leading dimension of the array B. LDB >= max(1,M). 🌸
X — DOUBLE PRECISION array, dimension (LDX,N). On input, X contains the right hand side Y. On output, X contains the solution of Equation (1) or (2) as selected by TRANS. Right hand side Y and the solution X are M-by-N matrices. ✨
LDX — INTEGER. The leading dimension of the array X. LDX >= max(1,M). 💜
SCALE — DOUBLE PRECISION. 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. 🌈
WORK — DOUBLE PRECISION array, dimension LWORK. Workspace for the algorithm. The workspace needs to be queried before running the computation — call the subroutine with INFO == -1 on input and the required workspace gets 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 156 of file dla_tglyap_l2_unopt.f90.

💅 `dla_tglyap_l3`

subroutine dla_tglyap_l3 (character, dimension(1) trans, integer m, double
    precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, *) x,
    integer ldx, double precision scale, double precision, dimension(*)
    work, integer info)

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

Purpose:

okay bestie so DLA_TGLYAP_L3 solves a generalized Lyapunov equation, and it comes in two flavors:

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

A^T * X * B + B^T * X * A =  SCALE * Y                                             (2)

where A is a M-by-M quasi upper triangular matrix, B is a M-by-M upper triangular, and X and Y are symmetric M-by-M matrices. Typically the matrix pencil (A,B) is created by DGGES from LAPACK. 🦋

Attention

The algorithm is implemented using BLAS level 3 operations. she's efficient, your honor ✨

Parameters

TRANS

TRANS is CHARACTER. It specifies the form of the system of equations with respect to A and B:

== 'N': Equation (1) is solved.
== 'T': Equation (2) is solved.

M

M is INTEGER. The order of the matrices A and B. M >= 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,M). The matrix B must be upper triangular.

LDB

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

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 154 of file dla_tglyap_l3.f90.

🌸 `dla_tglyap_l3_2s`

subroutine dla_tglyap_l3_2s (character, dimension(1) trans, integer m,
    double precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, *) x,
    integer ldx, double precision scale, double precision, dimension(*)
    work, integer info)

Level-3 Bartels-Stewart Algorithm for the generalized Lyapunov equation two stage.

Purpose:

so DLA_TGLYAP_L3_2S solves a generalized Lyapunov equation, same two forms as always:

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

A^T * X * B + B^T * X * A =  SCALE * Y                                             (2)

Attention

The algorithm is implemented using BLAS level 3 operations with two stage blocking (DAG inner solves). it's giving parallel energy ✨

Parameters

TRANS

TRANS is CHARACTER. It specifies the form of the system of equations with respect to A and B:

== 'N': Equation (1) is solved.
== 'T': Equation (2) is solved.

M

M is INTEGER. The order of the matrices A and B. M >= 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,M). The matrix B must be upper triangular.

LDB

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

X

LDX

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

SCALE

WORK

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 153 of file dla_tglyap_l3_2stage.f90.

🧚 `dla_tglyap_recursive`

recursive subroutine dla_tglyap_recursive (character, dimension(1) trans,
    integer m, double precision, dimension(lda, *) a, integer lda, double
    precision, dimension(ldb,*) b, integer ldb, double precision,
    dimension(ldx, *) x, integer ldx, double precision scale, double
    precision, dimension(*) work, integer info)

Recursive Blocking Algorithm for the generalized Lyapunov equation.

Purpose:

and now for the main character: DLA_TGLYAP_RECURSIVE solves a generalized Lyapunov equation in the usual two forms:

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

A^T * X * B + B^T * X * A =  SCALE * Y                                             (2)

Attention

The algorithm is implemented using recursive blocking. it calls itself, it's self-care fr fr 💅

Parameters

TRANS

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

== 'N': Equation (1) is solved.
== 'T': Equation (2) is solved.

M

M is INTEGER. The order of the matrices A and B. M >= 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,M). The matrix B must be upper triangular.

LDB

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

X

LDX

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

SCALE —

SCALE

okay so SCALE is DOUBLE PRECISION and it's literally the vibe preserver 💅 — it's a scaling factor whose whole job is to stop your result from overflowing. if INFO == 0 then SCALE is 1.0D0, no drama. but if one of the inner systems couldn't be solved correctly, then 0 < SCALE <= 1 holds true. so a SCALE that's not 1.0D0 is the algorithm quietly telling you "bestie i did what i could" 😭

WORK

WORK is a DOUBLE PRECISION array, dimension M*M, and it's just workspace for the algorithm. it's the backstage area, not the main character ✨

INFO

INFO is INTEGER, and here's the tea:

== 0: successful exit, we ate 👑
< 0: if INFO = -i, the i-th argument had an illegal value — you passed something cursed, go fix it 🙈
> 0: the equation is not solved correctly. one of the arising inner system got singular. not a slay. 💀

Author

Martin Koehler, MPI Magdeburg

Date

October 2025

Definition

at line 140 of file dla_tglyap_recursive.f90.

🦋 subroutine dla_tgstein_dag 🦋

okay bestie, next up, and she is giving parallel energy:

subroutine dla_tgstein_dag (character, dimension(1) trans, integer m, double
    precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldb,*) b, integer ldb, double precision, dimension(ldx, *) x,
    integer ldx, double precision scale, double precision, dimension(*)
    work, integer info)

Level-3 Bartels-Stewart Algorithm for the generalized Stein equation with DAG scheduling. the DAG is the whole point, she's out here running tasks in a graph, iykyk 🧚

Purpose

okay so DLA_TGSTEIN_DAG solves a generalized Lyapunov equation, and it comes in two flavors:

A * X * A^T - B * X * B^T = SCALE * Y                                              (1)

A^T * X * A - B^T * X * B =  SCALE * Y                                             (2)

where A is a M-by-M quasi upper triangular matrix, B is a M-by-M upper triangular, and X and Y are symmetric M-by-M matrices. and like, typically the matrix pencil (A,B) is created by DGGES from LAPACK. so if you got your (A,B) from DGGES, you're in the right place 🫶

Attention: the algorithm is implemented using BLAS level 3 operations and DAG . (yes, there's a space before the period, she's quirky like that 🌸)

Parameters

TRANS

TRANS is CHARACTER(1). it specifies the form of the system of equations with respect to A, B, and C: - == 'N': Equation (1) is solved. - == 'T': Equation (2) is solved.

M

M is INTEGER. the order of the matrices A and C. M >= 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,M). the matrix B must be upper triangular.

LDB

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

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. and listen, the workspace needs to be queried before running the computation — the query is performed by calling the subroutine with INFO == -1 on input, and the required workspace is then returned in INFO. so don't just roll up with a random WORK, do the query first bestie 🎀

INFO

INFO is INTEGER. on input: - == -1: Perform a workspace query - <> -1: normal operation on exit, workspace query: - < 0: if INFO = -i, the i-th argument had an illegal value - >= 0: the value of INFO is the required number of elements in the workspace. on exit, normal operation: - == 0: successful exit - < 0: if INFO = -i, the i-th argument had an illegal value - > 0: the equation is not solved correctly. one of the arising inner system got singular.

Author

Martin Koehler, MPI Magdeburg

Date

October 2025

Definition

at line 152 of file dla_tgstein_dag.f90.

🌈 subroutine dla_tgstein_l2 🌈

and now the level-2 sister, same energy, fewer BLAS levels:

subroutine dla_tgstein_l2 (character, dimension(1) trans, integer m, double
    precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, *) x,
    integer ldx, double precision scale, double precision, dimension(*)
    work, integer info)

Level-2 Bartels-Stewart Algorithm for the generalized Stein equation. periodt. 💜

Purpose

okay so DLA_TGSTEIN_L2 solves a generalized Lyapunov equation of the following forms:

A * X * A^T - B * X * B^T = SCALE * Y                                              (1)

A^T * X * A - B^T * X * B =  SCALE * Y                                             (2)

where A is a M-by-M quasi upper triangular matrix, B is a M-by-M upper triangular, and X and Y are symmetric M-by-M matrices. and typically the matrix pencil (A,B) is created by DGGES from LAPACK. iconic behavior honestly ⭐

Attention: the algorithm is implemented using BLAS level 2 operations.

Parameters

TRANS

TRANS is CHARACTER(1). it specifies the form of the system of equations with respect to A, B, and C: - == 'N': Equation (1) is solved. - == 'T': Equation (2) is solved.

M

M is INTEGER. the order of the matrices A and C. M >= 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,M). the matrix B must be upper triangular.

LDB

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

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. and yes, same deal as before: the workspace needs to be queried before running the computation — the query is performed by calling the subroutine with INFO == -1 on input, and the required workspace is then returned in INFO. always query, never guess 💅

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 154 of file dla_tgstein_l2.f90.

💖 subroutine dla_tgstein_l3 💖

and the trilogy concludes with the level-3 baddie:

subroutine dla_tgstein_l3 (character, dimension(1) trans, integer m, double
    precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, *) x,
    integer ldx, double precision scale, double precision, dimension(*)
    work, integer info)

Level-3 Bartels-Stewart Algorithm for the generalized Stein equation. the final boss 🔥

Purpose

okay so DLA_TGSTEIN_L3 solves a generalized Lyapunov equation of the following forms:

A * X * A^T - B * X * B^T = SCALE * Y                                              (1)

A^T * X * A - B^T * X * B =  SCALE * Y                                             (2)

Attention: the algorithm is implemented using BLAS level 3 operations.

Parameters

TRANS

TRANS is CHARACTER(1). it specifies the form of the system of equations with respect to A, B, and C: - == 'N': Equation (1) is solved. - == 'T': Equation (2) is solved.

M

M is INTEGER. the order of the matrices A and C. M >= 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 — B is a DOUBLE PRECISION array, dimension (LDB,M). Bestie, this matrix has to be upper triangular or we're not doing this. 💜
LDB — LDB is INTEGER, the leading dimension of the array B. Requirement: LDB >= max(1,M). No exceptions, she's strict. 🎀
X — X is a DOUBLE PRECISION array, dimension (LDX,N). On input, X is holding your right hand side Y like it's carrying your emotional baggage. On output, X contains the solution of Equation (1) or (2), whichever one you picked with TRANSA, TRANSB, and SGN. Both the right hand side Y and the solution X are M-by-N matrices. ✨
LDX — LDX is INTEGER, the leading dimension of the array X. Requirement: LDX >= max(1,M). Respect the leading dimension, it's the main character. 👑
SCALE — SCALE is DOUBLE PRECISION, and she's here purely to stop 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, then 0 < SCALE <= 1 holds true. Iconic damage-control behaviour honestly. 💅
WORK — WORK is a DOUBLE PRECISION array, dimension LWORK. It's the workspace for the algorithm. Heads up though: the workspace needs to be queried before running the computation. You do that by calling the subroutine with INFO == -1 on input, and then the required workspace comes back to you in INFO. Like a lil fortune teller but for memory. 🔮
INFO — INFO is 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 system got singular. 💀

Author

Martin Koehler, MPI Magdeburg

Date

October 2025

Definition at line 153 of file dla_tgstein_l3.f90.

🌈 `dla_tgstein_l3_2s` — Level-3 Bartels-Stewart Algorithm for the generalized Stein equation (Two Stage)

subroutine dla_tgstein_l3_2s (character, dimension(1) trans, integer m,
    double precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldb, *) b, integer ldb, double precision, dimension(ldx, *) x,
    integer ldx, double precision scale, double precision, dimension(*)
    work, integer info)

✨ Purpose ✨

DLA_TGSTEIN_L3_2S solves a generalized Lyapunov equation of the following forms:

A * X * A^T - B * X * B^T = SCALE * Y                                              (1)

A^T * X * A - B^T * X * B =  SCALE * Y                                             (2)

where A is a M-by-M quasi upper triangular matrix, B is a M-by-M upper triangular, and X and Y are symmetric M-by-M matrices. Typically the matrix pencil (A,B) is created by DGGES from LAPACK. 🦋

👀 Attention

The algorithm is implemented using BLAS level 3 operations with two level blocking.

🎀 Parameters

TRANS — TRANS is CHARACTER(1). Specifies the form of the system of equations with respect to A, B, and C:
- == 'N' : Equation (1) is solved.
- == 'T' : Equation (2) is solved.
M — M is INTEGER. The order of the matrices A and C. M >= 0.
A — A is a DOUBLE PRECISION array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular.
LDA — LDA is INTEGER, the leading dimension of the array A. Requirement: LDA >= max(1,M).
B — B is a DOUBLE PRECISION array, dimension (LDB,M). The matrix B must be upper triangular.
LDB — LDB is INTEGER, the leading dimension of the array B. Requirement: LDB >= max(1,M).
X — X is a DOUBLE PRECISION array, dimension (LDX,N). On input, the matrix X contains the right hand side Y. On output, the matrix X contains the solution of Equation (1) or (2) as selected by TRANSA, TRANSB, and SGN. Right hand side Y and the solution X are M-by-N matrices.
LDX — LDX is INTEGER, the leading dimension of the array X. Requirement: LDX >= max(1,M).
SCALE — SCALE is DOUBLE PRECISION. It's a scaling factor to prevent the overflow in the result. If INFO == 0 then SCALE is 1.0D0, otherwise if one of the inner systems could not be solved correctly, 0 < SCALE <= 1 holds true.
WORK — 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.
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 153 of file dla_tgstein_l3_2stage.f90.

🌸 `dla_tgstein_recursive` — Recursive Blocking Algorithm for the generalized Stein equation

recursive subroutine dla_tgstein_recursive (character, dimension(1) trans,
    integer m, double precision, dimension(lda, *) a, integer lda, double
    precision, dimension(ldb, *) b, integer ldb, double precision,
    dimension(ldx, *) x, integer ldx, double precision scale, double
    precision, dimension(*) work, integer info)

✨ Purpose ✨

DLA_TGSTEIN_RECURSIVE solves a generalized Lyapunov equation of the following forms:

A * X * A^T - B * X * B^T = SCALE * Y                                              (1)

A^T * X * A - B^T * X * B =  SCALE * Y                                             (2)

👀 Attention

The algorithm is implemented using recursive blocking.

🎀 Parameters

TRANS — TRANS is CHARACTER(1). Specifies the form of the system of equations with respect to A, B, and C:
- == 'N' : Equation (1) is solved.
- == 'T' : Equation (2) is solved.
M — M is INTEGER. The order of the matrices A and C. M >= 0.
A — A is a DOUBLE PRECISION array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular.
LDA — LDA is INTEGER, the leading dimension of the array A. Requirement: LDA >= max(1,M).
B — B is a DOUBLE PRECISION array, dimension (LDB,M). The matrix B must be upper triangular.
LDB — LDB is INTEGER, the leading dimension of the array B. Requirement: LDB >= max(1,M).
X — X is a DOUBLE PRECISION array, dimension (LDX,N). On input, the matrix X contains the right hand side Y. On output, the matrix X contains the solution of Equation (1) or (2) as selected by TRANSA, TRANSB, and SGN. Right hand side Y and the solution X are M-by-N matrices.
LDX — LDX is INTEGER, the leading dimension of the array X. Requirement: LDX >= max(1,M).
SCALE — SCALE is DOUBLE PRECISION. 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 a DOUBLE PRECISION array, dimension M*M. Workspace for the algorithm. No query dance needed this time, she just wants M*M. 🙈
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 140 of file dla_tgstein_recursive.f90.

💖 Author

Generated automatically by Doxygen for MEPACK from the source code.

And that's a wrap bestie!! 🫶 This whole lil masterpiece is MEPACK version 1.1.2, stamped Thu Oct 9 2025 00:00:00, and filed under dbltglyap(3). Go forth and solve those generalized Stein equations, you absolute legend. 💅✨👑

📦 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

dbltglyap(3)

💅 dbltglyap — Double Precision 💅

✨ NAME ✨

🌸 SYNOPSIS 🌸

🎀 Functions

🦋 DETAILED DESCRIPTION 🦋

📖 FUNCTION DOCUMENTATION 📖

🌈 dla_ggstein_refine

💫 Purpose

🔥 Parameters

🧚 See also

🌟 Author

🌙 Date

📌 Definition

💜 dla_tglyap_dag

💫 Purpose

⚠️ Attention

🔥 Parameters

Author

Date

Definition at line 154 of file dla_tglyap_dag.f90.

subroutine dla_tglyap_l2

Purpose:

Attention

Parameters

Author

Date

Definition at line 155 of file dla_tglyap_l2.f90.

subroutine dla_tglyap_l2_unopt

Purpose:

Attention

Parameters

Author

Date

Definition at line 156 of file dla_tglyap_l2_unopt.f90.

💅 dla_tglyap_l3

🌸 dla_tglyap_l3_2s

🧚 dla_tglyap_recursive

Author

Date

Definition

🦋 subroutine dla_tgstein_dag 🦋

Purpose

Parameters

Author

Date

Definition

🌈 subroutine dla_tgstein_l2 🌈

Purpose

Parameters

Author

Date

Definition

💖 subroutine dla_tgstein_l3 💖

Purpose

Parameters

Author

Date

🌈 dla_tgstein_l3_2s — Level-3 Bartels-Stewart Algorithm for the generalized Stein equation (Two Stage)

✨ Purpose ✨

👀 Attention

🎀 Parameters

Author

Date

🌸 dla_tgstein_recursive — Recursive Blocking Algorithm for the generalized Stein equation

✨ Purpose ✨

👀 Attention

🎀 Parameters

Author

Date

💖 Author

🌈 `dla_ggstein_refine`

💜 `dla_tglyap_dag`

💅 `dla_tglyap_l3`

🌸 `dla_tglyap_l3_2s`

🧚 `dla_tglyap_recursive`

🌈 `dla_tgstein_l3_2s` — Level-3 Bartels-Stewart Algorithm for the generalized Stein equation (Two Stage)

🌸 `dla_tgstein_recursive` — Recursive Blocking Algorithm for the generalized Stein equation