dbltrlyap(3)

📦
version 1.1.2 · 3 versions
, in fedora 43

💖 dbltrlyap(3) — MEPACK 💖

okay bestie, sit down, hydrate, because we're about to talk about dbltrlyap from MEPACK and it's honestly giving solver royalty ✨

🎀 NAME 🎀

dbltrlyap — Double Precision routines.

Like, the whole vibe is: double precision routines for standard Lyapunov and Stein equations with triangular coefficient matrices. no cap. that's the whole thing. 🦋

✨ SYNOPSIS ✨

🌸 Functions 🌸

  • subroutine dla_trlyap_dag (trans, m, a, lda, x, ldx, scale, work, info) — DAG Scheduled Bartels-Stewart Algorithm for the standard Lyapunov Equation. 💅
  • subroutine dla_trstein_dag (trans, m, a, lda, x, ldx, scale, work, info) — Level-3 Bartels-Stewart Algorithm for the Stein equation. 👑
  • subroutine dla_trlyap_kernel_44n (m, a, lda, x, ldx, scale, info) — Solver for a 4x4 standard Lyapunov equation (TRANS = N). 🍓
  • subroutine dla_trlyap_kernel_44t (m, a, lda, x, ldx, scale, info) — Solver for a 4x4 standard Lyapunov equation (TRANS = T). 🍓
  • subroutine dla_trlyap_l2 (transa, m, a, lda, x, ldx, scale, work, info) — Level-2 Bartels-Stewart Algorithm for the Lyapunov Equation. 🌈
  • subroutine dla_trlyap_l2_opt (transa, m, a, lda, x, ldx, scale, work, info) — Level-2 Bartels-Stewart Algorithm for the Lyapunov Equation (Optimized). ⚡
  • subroutine dla_trstein_l2 (trans, m, a, lda, x, ldx, scale, work, info) — Level-2 Bartels-Stewart Algorithm for the Stein equation. 💜
  • subroutine dla_trlyap_l3 (trans, m, a, lda, x, ldx, scale, work, info) — Level-3 Bartels-Stewart Algorithm for the standard Lyapunov Equation. 🔥
  • subroutine dla_trlyap_l3_2s (trans, m, a, lda, x, ldx, scale, work, info) — Level-3 Bartels-Stewart Algorithm for the standard Lyapunov Equation with 2 stage blocking. 🧚
  • subroutine dla_trstein_l3 (trans, m, a, lda, x, ldx, scale, work, info) — Level-3 Bartels-Stewart Algorithm for the Stein equation. 🌸
  • subroutine dla_trstein_l3_2s (trans, m, a, lda, x, ldx, scale, work, info) — Level-3 Bartels-Stewart Algorithm for the Stein equation with 2 stage blocking. 🧚
  • recursive subroutine dla_trlyap_recursive (trans, m, a, lda, x, ldx, scale, work, info) — Level-3 Recursive Blocked Algorithm for the standard Lyapunov Equation. 🌀
  • recursive subroutine dla_trstein_recursive (trans, m, a, lda, x, ldx, scale, work, info) — Recursive Blocked Algorithm for the Stein equation. 🌀

💫 Detailed Description 💫

Double precision routines for standard Lyapunov and Stein equations with triangular coefficient matrices. periodt. 💅

🫶 Function Documentation 🫶

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

DAG Scheduled Bartels-Stewart Algorithm for the standard Lyapunov Equation. okay so this one is the DAG scheduled girlie and she's iconic 💖

🌸 Purpose: 🌸

!> DLA_TRLYAP_DAG solves a Lyapunov equation of the following forms
!>
!>    A  * X  +  X * A**T = SCALE * Y                              (1)
!>
!> or
!>
!>    A ** T * X  +  X * A = SCALE * Y                              (2)
!>
!> where A is a M-by-M quasi upper triangular matrix.
!> The right hand side Y and the solution X are M-by-N matrices.  Typically the matrix A
!> is generated by DGEES from LAPACK.
!>

💅 Remarks 💅

The algorithm is implemented using DAG Scheduling. lowkey the whole personality tbh.

🎀 Parameters 🎀

  • TRANS
!>          TRANS is CHARACTER(1)
!>          Specifies the form of the system of equations with respect to A:
!>          == 'N':  op1(A) = A
!>          == 'T':  op1(A) = A**T
!>
  • 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).
!>
  • 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)
!>          Right hand side Y and the solution X are symmetric M-by-M 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.
!>

👑 Author 👑

Martin Koehler, MPI Magdeburg

🗓️ Date 🗓️

October 2025

Definition at line 137 of file dla_trlyap_dag.f90.

⭐ subroutine dla_trlyap_kernel_44n (integer m, double precision, dimension(lda, *) a, integer lda, double precision, dimension(ldx, *) x, integer ldx, double precision scale, integer info) ⭐

Solver for a 4x4 standard Lyapunov equation (TRANS = N). tiny lil 4x4, so cute, so fast 🍓

🌸 Purpose: 🌸

!>
!> DLA_TRLYAP_KERNEL_44N solves a Lyapunov  equation of the following form
!>
!>    A * X  + X * A**T = SCALE * Y                              (1)
!>
!> where A is a M-by-M quasi upper triangular matrix.
!> The right hand side Y and the solution X M-by-M matrices.
!> Typically the matrix  A is created by DGEES from LAPACK.
!> The algorithm is implemented without BLAS level 2
!> operations. Thereby the order of M and N is at most 4. Furthermore, for fast execution
!> the function does not check the input arguments.
!>
!>
⚠️ bestie. it says the function does not check the input arguments. so like, you GOTTA have your inputs correct yourself. no safety net. we die like the fast path we are 💀

🎀 Parameters 🎀

  • M
!>          M is INTEGER
!>          The order of the matrices A and C.  4 >= 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).
!>
  • 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 < SCAL <= 1 holds true.
!>
  • INFO
!>          INFO is INTEGER
!>          On output:
!>          == 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  The equation is not solved correctly. One of the arising inner
!>                system got singular.
!>

👑 Author 👑

Martin Koehler, MPI Magdeburg

🗓️ Date 🗓️

October 2025

Definition at line 109 of file dla_trlyap_kernel_44_n.f90.

⭐ subroutine dla_trlyap_kernel_44t (integer m, double precision, dimension(lda, *) a, integer lda, double precision, dimension(ldx, *) x, integer ldx, double precision scale, integer info) ⭐

Solver for a 4x4 standard Lyapunov equation (TRANS = T). the transposed twin of 44n, serving looks 💅✨

🌸 Purpose: 🌸

!>
!> DLA_TRLYAP_KERNEL_44T solves a Lyapunov  equation of the following form
!>
!>    A **T * X  + X * A = SCALE * Y                              (1)
!>
!> where A is a M-by-M quasi upper triangular matrix.
!> The right hand side Y and the solution X M-by-M matrices.
!> Typically the matrix  A is created by DGEES from LAPACK.
!> The algorithm is implemented without BLAS level 2
!> operations. Thereby the order of M and N is at most 4. Furthermore, for fast execution
!> the function does not check the input arguments.
!>
!>
⚠️ same warning bestie: the function does not check the input arguments. you've been warned fr fr 💀

🎀 Parameters 🎀

  • M
!>          M is INTEGER
!>          The order of the matrices A and C.  4 >= 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).
!>
  • 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 < SCAL <= 1 holds true.
!>
  • INFO

(bestie the page cuts off right here mid-INFO and I am NOT about to invent the rest — part 2 is coming, stay tuned 🦋)

  • INFO — is INTEGER 🫶

On output:

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

✍️ Author

Martin Koehler, MPI Magdeburg

📅 Date

October 2025

📌 Definition

at line 109 of file dla_trlyap_kernel_44_t.f90.

🦋 dla_trlyap_l2 — Level-2 Bartels-Stewart Algorithm for the Lyapunov Equation 🦋

subroutine dla_trlyap_l2 (character, dimension(1) transa, integer m, double
    precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldx, *) x, integer ldx, double precision scale, double
    precision, dimension(*) work, integer info)

okay bestie, this one is the level-2 glow-up of the Bartels-Stewart algorithm for the Lyapunov equation 💅

🎯 Purpose

DLA_TRLYAP_L2 solves a Lyapunov equation of the following forms:

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

or

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

where A is a M-by-M quasi upper triangular matrix. The right hand side Y and the solution X are M-by-N matrices. Typically the matrix A is generated by DGEES from LAPACK. she's giving LAPACK-origin realness, no cap 🔥

📝 Remarks

  • The algorithm is implemented using BLAS level 2 operations. level 2, level iconic
  • The transposed case (2) is optimized w.r.t. to the usage of the DSYR2 operation. ✨

🎀 Parameters

  • TRANSA — is CHARACTER(1)

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

  • == 'N': op1(A) = A
  • == 'T': op1(A) = A**T
  • M — is INTEGER

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

  • A — is DOUBLE PRECISION array, dimension (LDA,M)

The matrix A must be (quasi-) upper triangular. quasi, meaning she's got those lil 2x2 blocks living their truth

  • LDA — is INTEGER

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

  • 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). Right hand side Y and the solution X are symmetric M-by-M matrices. in goes Y, out comes the solution, it's a whole transformation arc 🦋

  • LDX — is INTEGER

The leading dimension of the array X. LDB >= max(1,M). (yes it literally says LDB in there, i'm just the messenger bestie 🙈)

  • 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 scales so your numbers don't explode, very protective, very mindful

  • 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 gotta ASK first, don't just show up 💅

  • 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 138 of file dla_trlyap_l2.f90.

🌈 dla_trlyap_l2_opt — Level-2 Bartels-Stewart Algorithm for the Lyapunov Equation (Optimized) 🌈

subroutine dla_trlyap_l2_opt (character, dimension(1) transa, integer m,
    double precision, dimension(lda, *) a, integer lda, double precision,
    dimension(ldx, *) x, integer ldx, double precision scale, double
    precision, dimension(*) work, integer info)

the optimized sister of the last one — same energy, extra polish 👑

🎯 Purpose

DLA_TRLYAP_L2 solves a Lyapunov equation of the following forms:

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

or

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

where A is a M-by-M quasi upper triangular matrix. The right hand side Y and the solution X are M-by-N matrices. Typically the matrix A is generated by DGEES from LAPACK. same facts, same vibes, just optimized ✨

📝 Remarks

  • The algorithm is implemented using BLAS level 2 operations.
  • The transposed case (2) is optimized w.r.t. to the usage of the DSYR2 operation. 🔥

🎀 Parameters

  • TRANSA — is CHARACTER(1)

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

  • == 'N': op1(A) = A
  • == 'T': op1(A) = A**T
  • M — is INTEGER

The order of the matrices A and C. M >= 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).

  • 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). Right hand side Y and the solution X are symmetric M-by-M matrices.

  • LDX — is INTEGER

The leading dimension of the array X. LDB >= max(1,M). (still says LDB, still not my fault 🙈)

  • 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. anti-overflow queen behavior

  • WORK — is DOUBLE PRECISION array, dimension LWORK

Workspace for the algorithm. The workspace needs to queried before the running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO.

  • INFO — is INTEGER

On input:

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

On exit, workspace query:

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

On exit, normal operation:

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

✍️ Author

Martin Koehler, MPI Magdeburg

📅 Date

October 2025

📌 Definition

at line 138 of file dla_trlyap_l2_opt.f90.

💜 dla_trlyap_l3 — Level-3 Bartels-Stewart Algorithm for the standard Lyapunov Equation 💜

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

now THIS one is level 3, she moved up in the world, she's doing level-3 BLAS things 🎀

🎯 Purpose

DLA_TRLYAP_L3 solves a Lyapunov equation of the following forms:

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

or

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

where A is a M-by-M quasi upper triangular matrix. The right hand side Y and the solution X are M-by-N matrices. Typically the matrix A is generated by DGEES from LAPACK. yes bestie, case (2) is a PLUS here, not a minus — read it again, it's right there 👀

📝 Remarks

  • The algorithm is implemented using BLAS level 3 operations. level 3 baby, we're in the big leagues now 🔥

🎀 Parameters

  • TRANS — is CHARACTER

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

  • == 'N': op1(A) = A
  • == 'T': op1(A) = A**T
  • M — is INTEGER

The order of the matrices A and C. M >= 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).

  • 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). Right hand side Y and the solution X are symmetric M-by-M matrices. same in-and-out deal as her level-2 sisters 🦋

  • LDX — is INTEGER

The leading dimension of the array X. LDB >= max(1,M). (yep. LDB. again. we're all just living with 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. keeping your doubles from going feral

  • 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 later, that's the rule 💅

  • 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. 💀
  • INFO — an INTEGER, and honestly she's the whole narrator of this story 👑

On input:

  • == -1 — bestie, this means do a workspace query ✨
  • <> -1 — normal operation, business as usual

On exit, workspace query flavour:

  • < 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 flavour:

  • == 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 🫠

👑 Author

Martin Koehler, MPI Magdeburg

🗓️ Date

October 2025

📍 Definition

At line 137 of file dla_trlyap_l3.f90.

💖 dla_trlyap_l3_2s

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

Level-3 Bartels-Stewart Algorithm for the standard Lyapunov Equation with 2 stage blocking. Two stages, zero chill. 🦋

🎯 Purpose

DLA_TRLYAP_L3_2S solves a Lyapunov equation of the following forms:

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

or

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

where A is a M-by-M quasi upper triangular matrix.

The right hand side Y and the solution X are M-by-N matrices. Typically the matrix A is generated by DGEES from LAPACK. No cap. 💜

💭 Remarks

The algorithm is implemented using BLAS level 3 operations and DAG schedule block solves. That's the whole vibe — level 3, DAG-scheduled, iconic. ✨

🧩 Parameters

  • TRANS — a CHARACTER, and she picks which form of the system of equations we're solving with respect to A:
    • == 'N' — op1(A) = A
    • == 'T' — op1(A) = A**T
  • M — an INTEGER. The order of the matrices A and C. M >= 0. Periodt. 🎀
  • A — a DOUBLE PRECISION array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular, bestie, don't try anything else 🙈
  • LDA — an INTEGER. The leading dimension of the array A. LDA >= max(1,M).
  • X — 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).
    • Right hand side Y and the solution X are symmetric M-by-M matrices. 💖
  • LDX — an INTEGER. The leading dimension of the array X. LDB >= max(1,M). (yes it says LDB, we are simply passing the message along 💀)
  • 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. She's here to catch you before you overflow 🌸
  • WORK — a DOUBLE PRECISION array, dimension LWORK. Workspace for the algorithm.
    • The workspace needs to queried before the running the computation. (verbatim, she's doing her best 😭)
    • 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

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 137 of file dla_trlyap_l3_2stage.f90.

🌈 dla_trlyap_recursive

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

Level-3 Recursive Blocked Algorithm for the standard Lyapunov Equation. She recurses. She blocks. She slays. 🧚

🎯 Purpose

DLA_TRLYAP_RECURSIVE solves a Lyapunov equation of the following forms:

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

or

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

where A is a M-by-M quasi upper triangular matrix.

The right hand side Y and the solution X are M-by-N matrices. Typically the matrix A is generated by DGEES from LAPACK. 🔥

💭 Remarks

The algorithm is implemented using BLAS level 3 operations and recursive blocking. Lowkey the most elegant thing on this whole page, ngl 💅

🧩 Parameters

  • TRANS — a CHARACTER(1). Specifies the form of the system of equations with respect to A:
    • == 'N' — op1(A) = A
    • == 'T' — op1(A) = A**T
  • M — an INTEGER. The order of the matrices A and C. M >= 0.
  • A — a DOUBLE PRECISION array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular. 🎀
  • LDA — an INTEGER. The leading dimension of the array A. LDA >= max(1,M).
  • X — 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).
    • Right hand side Y and the solution X are symmetric M-by-M matrices. ✨
  • LDX — an INTEGER. The leading dimension of the array X. LDB >= max(1,M). (again, LDB — iconic typo, we report, we do not judge 🙈)
  • 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. 💖
  • WORK — a DOUBLE PRECISION array, dimension 1. Workspace for the algorithm. Just one, bestie. Recursion carries the rest 🦋
  • INFO — an 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.
Heads up: this one has no workspace query mode. The INFO == -1 thing is simply not a thing here, don't come for me 💀

👑 Author

Martin Koehler, MPI Magdeburg

🗓️ Date

October 2025

📍 Definition

At line 125 of file dla_trlyap_recursive.f90.

⭐ dla_trstein_dag

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

Level-3 Bartels-Stewart Algorithm for the Stein equation. New equation, same impeccable energy 👑

🎯 Purpose

DLA_TRSTEIN_DAG solves a Stein equation of the following forms:

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

or

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

where A is a M-by-M quasi upper triangular matrix. And X and Y are symmetric M-by-M matrices. Typically the matrix A is created by DGEES from LAPACK.

Not me pointing out that the labels go (2) then (1). That's just how it is, bestie, I am merely the messenger 🫶

🚨 Attention

The algorithm is implemented using BLAS level 3 operations.

💭 Remarks

The algorithm is implemented using DAG Scheduling. Two whole notes about implementation, she's thorough like that ✨

🧩 Parameters

  • TRANS — a CHARACTER(1). Specifies the form of the system of equations with respect to A and C:
    • == 'N' — Equation (1) is solved.
    • == 'T' — Equation (2) is solved.
  • M — an INTEGER. The order of the matrices A and C. M >= 0. 🎀
  • A — a DOUBLE PRECISION array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular.
  • LDA — an INTEGER. The leading dimension of the array A. LDA >= max(1,M).
  • X — 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. (yes, the page says TRANSA, TRANSB and SGN here — I said what I said 💀)
    • Right hand side Y and the solution X are M-by-N matrices. 💜
  • LDX — an INTEGER. The leading dimension of the array X. LDX >= max(1,M). Look at that, actual LDX this time, growth 💅
  • 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. 🌸
  • WORK — a DOUBLE PRECISION array, dimension LWORK. Workspace for the algorithm.
    • The workspace needs to queried before the running the computation.
    • The query is performed by calling the subroutine with INFO == -1 on input.
    • The required workspace is then returned in INFO. 🦋
  • INFO —
  • INFO — it's an INTEGER, bestie, and it's literally the vibe check of the entire call 💅

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 142 of file dla_trstein_dag.f90.

dla_trstein_l2

Level-2 Bartels-Stewart Algorithm for the Stein equation.

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

Purpose:

DLA_TRSTEIN_L2 solves a generalized Lyapunov equation of the following forms, okurrr 🦋

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

or

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

where A is a M-by-M quasi upper triangular matrix, and X and Y are symmetric M-by-M matrices. Typically the matrix A is created by DGEES from LAPACK. Periodt. ✨

Attention

The algorithm is implemented using BLAS level 2 operations. No cap, that's just how it is 🫶

Parameters

  • TRANS — it's a CHARACTER(1). Specifies the form of the system of equations with respect to A and C:
    • == 'N': Equation (1) is solved.
    • == 'T': Equation (2) is solved.
  • M — it's an INTEGER. The order of the matrices A and C. M >= 0.
  • A — it's a DOUBLE PRECISION array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular.
  • LDA — it's an INTEGER. The leading dimension of the array A. LDA >= max(1,M).
  • X — it's 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 — it's an INTEGER. The leading dimension of the array X. LDX >= max(1,M).
  • SCALE — it's 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.
  • WORK — it's a DOUBLE PRECISION array, dimension LWORK. Workspace for the algorithm. The workspace needs to queried before the running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO.
  • INFO — it's an 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 143 of file dla_trstein_l2.f90.

dla_trstein_l3

Level-3 Bartels-Stewart Algorithm for the Stein equation.

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

Purpose:

DLA_TRSTEIN_L3 solves a standard Stein equation of the following forms, bestie 💖

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

or

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

where A is a M-by-M quasi upper triangular matrix, and X and Y are symmetric M-by-M matrices. Typically the matrix A is created by DGEES from LAPACK. Slay. 🔥

Attention

The algorithm is implemented using BLAS level 3 operations. She leveled up, fr fr 🧚

Parameters

  • TRANS — it's a CHARACTER. Specifies the form of the system of equations with respect to A and C:
    • == 'N': Equation (1) is solved.
    • == 'T': Equation (2) is solved.
  • M — it's an INTEGER. The order of the matrices A and C. M >= 0.
  • A — it's a DOUBLE PRECISION array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular.
  • LDA — it's an INTEGER. The leading dimension of the array A. LDA >= max(1,M).
  • X — it's 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 — it's an INTEGER. The leading dimension of the array X. LDX >= max(1,M).
  • SCALE — it's 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.
  • WORK — it's a DOUBLE PRECISION array, dimension LWORK. Workspace for the algorithm. The workspace needs to queried before the running the computation. The query is performed by calling the subroutine with INFO == -1 on input. The required workspace is then returned in INFO.
  • INFO — it's an 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 141 of file dla_trstein_l3.f90.

dla_trstein_l3_2s

Level-3 Bartels-Stewart Algorithm for the Stein equation with 2 stage blocking.

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

Purpose:

DLA_TRSTEIN_L3_2S solves a generalized Lyapunov equation of the following forms, and she's doing it in two stages because she's that girl 👑

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

or

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

where A is a M-by-M quasi upper triangular matrix, and X and Y are symmetric M-by-M matrices. Typically the matrix A is created by DGEES from LAPACK. Iconic behaviour only. 💅

Attention

The algorithm is implemented using BLAS level 3 operations and a DAG scheduled inner solver. Yes bestie, a whole dependency graph, we're not playing games 😭

Parameters

  • TRANS — it's a CHARACTER(1). Specifies the form of the system of equations with respect to A and C:
    • == 'N': Equation (1) is solved.
    • == 'T': Equation (2) is solved.
  • M — it's an INTEGER. The order of the matrices A and C. M >= 0.
  • A — it's a DOUBLE PRECISION array, dimension (LDA,M). The matrix A must be (quasi-) upper triangular.
  • LDA — it's an INTEGER. The leading dimension of the array A. LDA >= max(1,M).
  • X — it's 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 — it's an INTEGER. The leading dimension of the array X. LDX >= max(1,M).
  • SCALE — it's 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.
  • WORK
  • WORK — DOUBLE PRECISION array, dimension LWORK ✨ it's the workspace for the algorithm, bestie. Heads up: you gotta query that workspace before you run the computation. You do the query by calling the subroutine with INFO == -1 on input, and the required workspace comes right back at you in INFO. 💅
  • INFO — INTEGER 🦋

On input:

  • == -1 — perform a workspace query 🔍
  • <> -1 — normal operation, we're actually doing the thing ✨

On exit, workspace query:

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

On exit, normal operation:

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

Author

Martin Koehler, MPI Magdeburg 👑

Date

October 2025 🌸

Definition at line 141 of file dla_trstein_l3_2stage.f90.

🦋 recursive subroutine dla_trstein_recursive 🦋

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

Recursive Blocked Algorithm for the Stein equation. 🌈

✨ Purpose ✨

DLA_TRSTEIN_RECURSIVE solves a generalized Lyapunov equation of the following forms
A * X * A^T - X  = SCALE * Y                                              (2)

or

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

where A is a M-by-M quasi upper triangular matrix. 🎀 and X and Y are symmetric M-by-M matrices. 💜 Typically the matrix A is created by DGEES from LAPACK. ⭐

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

  • 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 M*N Workspace for the algorithm. ✨

  • INFO

INFO is INTEGER

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

Author

Martin Koehler, MPI Magdeburg 👑

Date

October 2025 🌸

Definition at line 129 of file dla_trstein_recursive.f90.

💫 Author 💫

Generated automatically by Doxygen for MEPACK from the source code. 🤖✨

okay bestie that's the whole entire page, no cap 🫶 this was Version 1.1.2, dated Thu Oct 9 2025 00:00:00, and the page is dbltrlyap(3) — literally iconic, go solve something symmetric, periodt 💅🎀

📦 where this page comes from