# sgltgsylv - Man Page

## Name

sgltgsylv — Single Precision

— Single Precision routines for triangular generalized Sylvester equations.

## Synopsis

### Functions

subroutine sla_tgcsylv_dag (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 with DAG scheduling.
subroutine sla_tgcsylv_dual_dag (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info)
OpenMP4 - DAG Algorithm for the dual generalized coupled Sylvester equation.
subroutine sla_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.
subroutine sla_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 (TRANSA = N, TRANSB = N)
subroutine sla_tgcsylv_dual_kernel_44nt (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 (TRANSA = N, TRANSB = T)
subroutine sla_tgcsylv_dual_kernel_44tn (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 (TRANSA = T, TRANSB = N)
subroutine sla_tgcsylv_dual_kernel_44tt (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 (TRANSA = T, TRANSB = T)
subroutine sla_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.
subroutine sla_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 (Local Copy)
subroutine sla_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)
subroutine sla_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)
subroutine sla_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)
subroutine sla_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)
subroutine sla_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 (Optimized)
subroutine sla_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 Algorithm for the generalized coupled Sylvester equation (Optimized)
subroutine sla_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)
subroutine sla_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)
subroutine sla_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)
subroutine sla_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)
subroutine sla_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 Algorithm for the generalized coupled Sylvester equation (Optimized)
subroutine sla_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 Algorithm for the generalized coupled Sylvester equation (unoptimized)
subroutine sla_tgsylv_gardiner_laub (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info)
Level-2 Gardiner-Laub algorithm for the generalized Sylvester equation.
subroutine sla_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 sla_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 Algorithm for the generalized Sylvester equation (unoptimized)
subroutine sla_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 Algorithm for the generalized Sylvester equation (Optimized)
subroutine sla_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)
subroutine sla_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)
subroutine sla_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)
subroutine sla_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)
subroutine sla_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 Algorithm for the generalized Sylvester equation (Optimized)
subroutine sla_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 Algorithm for the generalized Sylvester equation (unoptimized)
subroutine sla_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.
subroutine sla_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.
subroutine sla_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)
subroutine sla_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 sla_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)
subroutine sla_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.
subroutine sla_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.
subroutine sla_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.
subroutine sla_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.
recursive subroutine sla_tgcsylv_dual_recursive (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info)
Recursive blocking Algorithm for the dual generalized coupled Sylvester equation.
recursive subroutine sla_tgcsylv_recursive (transa, transb, sgn1, sgn2, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, work, info)
Recursive Algorithm for the generalized coupled Sylvester equation.
recursive subroutine sla_tgsylv_recursive (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info)
Level-3 Recursive Blocked Solver for the generalized Sylvester equation.

## Detailed Description

Single Precision routines for triangular generalized Sylvester equations.

This section contains the solvers for the standard Sylvester equation with (quasi) triangular coefficient matrices. The coefficient matrices are normally generated with the help of the Schur decomposition from LAPACK. All routines work in single precision arithmetic.

## Function Documentation

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

Level-3 Kagstroem-Westin/Bartels-Stewart Algorithm for the generalized coupled Sylvester equation with DAG scheduling.

Purpose:

``` SLA_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 SGGES from LAPACK.```
Remarks

The Algorithm uses OpenMP DAG scheduling

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 REAL, allowed values: +/-1
Specifies the sign between in the first equation.```

SGN2

```          SGN2 is REAL, allowed values: +/-1
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 REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

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

WORK

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

INFO

```          INFO is INTEGER

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

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

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

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 217 of file sla_tgcsylv_dag.f90.

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

OpenMP4 - DAG Algorithm for the dual generalized coupled Sylvester equation.

Purpose:

``` SLA_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 SGGES 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

```          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 REAL, allowed values: +/-1
Specifies the sign between in the first equation.```

SGN2

```          SGN2 is REAL, allowed values: +/-1
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 REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

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

WORK

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

INFO

```          INFO is INTEGER

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

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

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

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 237 of file sla_tgcsylv_dual_dag.f90.

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

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

Purpose:

``` SLA_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
M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) ( or (D,B) ) are
created by SGGES 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 and op2(H) = H. The primary use of this function is to end the recursion in dla_tgcsylv_dual_recursive.

Parameters

SGN1

```          SGN1 is REAL, allowed values: +/-1
Specifies the sign between in the first equation.```

SGN2

```          SGN2 is REAL, 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 REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

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

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

Januar 2023

Definition at line 206 of file sla_tgcsylv_dual_kernel_44_nn.f90.

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

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

Purpose:

``` SLA_TGCSYLV_DUAL_KERNEL_44NT 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 SGGES 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 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 REAL, allowed values: +/-1
Specifies the sign between in the first equation.```

SGN2

```          SGN2 is REAL, 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 REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

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

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

Januar 2023

Definition at line 206 of file sla_tgcsylv_dual_kernel_44_nt.f90.

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

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

Purpose:

``` SLA_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 SGGES 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 REAL, allowed values: +/-1
Specifies the sign between in the first equation.```

SGN2

```          SGN2 is REAL, 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 REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

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

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

Januar 2023

Definition at line 206 of file sla_tgcsylv_dual_kernel_44_tn.f90.

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

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

Purpose:

``` SLA_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 SGGES 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.

Parameters

SGN1

```          SGN1 is REAL, allowed values: +/-1
Specifies the sign between in the first equation.```

SGN2

```          SGN2 is REAL, 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 REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

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

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

Januar 2023

Definition at line 206 of file sla_tgcsylv_dual_kernel_44_tt.f90.

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

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

Purpose:

``` SLA_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 SGGES 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.

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 REAL, allowed values: +/-1
Specifies the sign between in the first equation.```

SGN2

```          SGN2 is REAL, allowed values: +/-1
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 REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

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

WORK

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

INFO

```          INFO is INTEGER

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

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

On exit, normal operation:
== 0:  successful exit
< 0:  if INFO = -i, the i-th argument had an illegal value
> 0:  The equation is not solved correctly. One of the arising inner
system got singular.```
Optimizations:
• Replaced SAXPY operation by Fortran intrinsic
• Replaced all BLAS calls except of SGER by Fortran intrinsic,
Author

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 245 of file sla_tgcsylv_dual_l2.f90.

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

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

Purpose:

``` SLA_TGCSYLV_DUAL_L2_LOCAL_COPY 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 SGGES 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 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 REAL, allowed values: +/-1
Specifies the sign between in the first equation.```

SGN2

```          SGN2 is REAL, 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 REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

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

WORK

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

INFO

```          INFO is INTEGER

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

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

On exit, normal operation:
== 0:  successful exit
< 0:  if INFO = -i, the i-th argument had an illegal value
> 0:  The equation is not solved correctly. One of the arising inner
system got singular.```
Optimizations:
• Replaced SAXPY 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

Januar 2023

Definition at line 247 of file sla_tgcsylv_dual_l2_opt_local_copy.f90.

### subroutine sla_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:

``` SLA_TGCSYLV_DUAL_L2_LOCAL_COPY_128 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 SGGES 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 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

```          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 REAL, allowed values: +/-1
Specifies the sign between in the first equation.```

SGN2

```          SGN2 is REAL, 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 REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

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

WORK

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

INFO

```          INFO is INTEGER

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

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

On exit, normal operation:
== 0:  successful exit
< 0:  if INFO = -i, the i-th argument had an illegal value
> 0:  The equation is not solved correctly. One of the arising inner
system got singular.```
Optimizations:
• Replaced SAXPY 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

Januar 2023

Definition at line 248 of file sla_tgcsylv_dual_l2_opt_local_copy_128.f90.

### subroutine sla_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:

``` SLA_TGCSYLV_DUAL_L2_LOCAL_COPY_32 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 SGGES 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 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 REAL, allowed values: +/-1
Specifies the sign between in the first equation.```

SGN2

```          SGN2 is REAL, 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 REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

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

WORK

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

INFO

```          INFO is INTEGER

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

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

On exit, normal operation:
== 0:  successful exit
< 0:  if INFO = -i, the i-th argument had an illegal value
> 0:  The equation is not solved correctly. One of the arising inner
system got singular.```
Optimizations:
• Replaced SAXPY 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

Januar 2023

Definition at line 247 of file sla_tgcsylv_dual_l2_opt_local_copy_32.f90.

### subroutine sla_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:

``` SLA_TGCSYLV_DUAL_L2_LOCAL_COPY_64 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 SGGES 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 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.

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 REAL, allowed values: +/-1
Specifies the sign between in the first equation.```

SGN2

```          SGN2 is REAL, 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 REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

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

WORK

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

INFO

```          INFO is INTEGER

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

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

On exit, normal operation:
== 0:  successful exit
< 0:  if INFO = -i, the i-th argument had an illegal value
> 0:  The equation is not solved correctly. One of the arising inner
system got singular.```
Optimizations:
• Replaced SAXPY 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

Januar 2023

Definition at line 246 of file sla_tgcsylv_dual_l2_opt_local_copy_64.f90.

### subroutine sla_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:

``` SLA_TGCSYLV_DUAL_L2_LOCAL_COPY_96 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 SGGES 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 uses local copies of the input data. For this reason M and N are limited to 96. 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 REAL, allowed values: +/-1
Specifies the sign between in the first equation.```

SGN2

```          SGN2 is REAL, 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 REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

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

WORK

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

INFO

```          INFO is INTEGER

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

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

On exit, normal operation:
== 0:  successful exit
< 0:  if INFO = -i, the i-th argument had an illegal value
> 0:  The equation is not solved correctly. One of the arising inner
system got singular.```
Optimizations:
• Replaced SAXPY 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

Januar 2023

Definition at line 246 of file sla_tgcsylv_dual_l2_opt_local_copy_96.f90.

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

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

Purpose:

``` SLA_TGCSYLV_DUAL_L3 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 SGGES 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

```          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 REAL, allowed values: +/-1
Specifies the sign between in the first equation.```

SGN2

```          SGN2 is REAL, allowed values: +/-1
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 REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

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

WORK

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

INFO

```          INFO is INTEGER

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

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

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

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 236 of file sla_tgcsylv_dual_l3.f90.

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

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

Purpose:

``` SLA_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 SGGES 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

```          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 REAL, allowed values: +/-1
Specifies the sign between in the first equation.```

SGN2

```          SGN2 is REAL, allowed values: +/-1
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 REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

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

WORK

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

INFO

```          INFO is INTEGER

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

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

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

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 237 of file sla_tgcsylv_dual_l3_2stage.f90.

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

Recursive blocking Algorithm for the dual generalized coupled Sylvester equation.

Purpose:

``` SLA_TGCSYLV_DUAL_RECURSIVE 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 SGGES 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

```          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 REAL, allowed values: +/-1
Specifies the sign between in the first equation.```

SGN2

```          SGN2 is REAL, allowed values: +/-1
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 REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

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

WORK

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

INFO

```          INFO is INTEGER

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

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

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

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 240 of file sla_tgcsylv_dual_recursive.f90.

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

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

Purpose:

``` SLA_TGCSYLV_L2 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 SGGES from LAPACK.```
Remarks

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

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 REAL, allowed values: +/-1
Specifies the sign between in the first equation.```

SGN2

```          SGN2 is REAL, allowed values: +/-1
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 REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

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

WORK

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

INFO

```          INFO is INTEGER

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

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

On exit, normal operation:
== 0:  successful exit
< 0:  if INFO = -i, the i-th argument had an illegal value
> 0:  The equation is not solved correctly. One of the arising inner
system got singular.```
Optimizations:
• Replaced level-3 by level-2 and level-1 calls,
• Replaced SAXPY operation by Fortran intrinsic
• Replaced all BLAS calls except of SGER by Fortran intrinsic,
Author

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 225 of file sla_tgcsylv_l2.f90.

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

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

Purpose:

``` SLA_TGCSYLV_L2_LOCAL_COPY 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 SGGES from LAPACK.```
Remarks

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

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 REAL, allowed values: +/-1
Specifies the sign between in the first equation.```

SGN2

```          SGN2 is REAL, allowed values: +/-1
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 REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

```          C is REAL 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 REAL 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 REAL 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.```

LDE

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

F

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

WORK

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

INFO

```          INFO is INTEGER

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

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

On exit, normal operation:
== 0:  successful exit
< 0:  if INFO = -i, the i-th argument had an illegal value
> 0:  The equation is not solved correctly. One of the arising inner
system got singular.```
Optimizations:
• Replaced level-3 by level-2 and level-1 calls,
• Reorder the solution order to column first style,
• Replaced SAXPY operation by Fortran intrinsic
• Replaced all BLAS calls except of STRMV and SGER by Fortran intrinsic,
• Use local copies of A, B, C, D, E, and F (M, N <=128) .
Author

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 227 of file sla_tgcsylv_l2_opt_local_copy.f90.

### subroutine sla_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:

``` SLA_TGCSYLV_L2_LOCAL_COPY_128 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 SGGES from LAPACK.```
Remarks

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

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 REAL, allowed values: +/-1
Specifies the sign between in the first equation.```

SGN2

```          SGN2 is REAL, 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 REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

```          C is REAL 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 REAL 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 REAL 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.```

LDE

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

F

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

WORK

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

INFO

```          INFO is INTEGER

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

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

On exit, normal operation:
== 0:  successful exit
< 0:  if INFO = -i, the i-th argument had an illegal value
> 0:  The equation is not solved correctly. One of the arising inner
system got singular.```
Optimizations:
• Replaced level-3 by level-2 and level-1 calls,
• Reorder the solution order to column first style,
• Replaced SAXPY operation by Fortran intrinsic
• Replaced all BLAS calls except of STRMV and SGER by Fortran intrinsic,
• Align local copies and fix problem size to <= 128
Author

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 225 of file sla_tgcsylv_l2_opt_local_copy_128.f90.

### subroutine sla_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:

``` SLA_TGCSYLV_L2_LOCAL_COPY_32 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 SGGES from LAPACK.```
Remarks

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

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 REAL, allowed values: +/-1
Specifies the sign between in the first equation.```

SGN2

```          SGN2 is REAL, 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 REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

```          C is REAL 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 REAL 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 REAL 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.```

LDE

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

F

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

WORK

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

INFO

```          INFO is INTEGER

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

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

On exit, normal operation:
== 0:  successful exit
< 0:  if INFO = -i, the i-th argument had an illegal value
> 0:  The equation is not solved correctly. One of the arising inner
system got singular.```
Optimizations:
• Replaced level-3 by level-2 and level-1 calls,
• Reorder the solution order to column first style,
• Replaced SAXPY operation by Fortran intrinsic
• Replaced all BLAS calls except of STRMV and SGER by Fortran intrinsic,
• Align local copies and fix problem size to <= 32
Author

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 226 of file sla_tgcsylv_l2_opt_local_copy_32.f90.

### subroutine sla_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)

Purpose:

``` SLA_TGCSYLV_L2_LOCAL_COPY_64 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 SGGES from LAPACK.```
Remarks

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

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 REAL, allowed values: +/-1
Specifies the sign between in the first equation.```

SGN2

```          SGN2 is REAL, 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 REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

```          C is REAL 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 REAL 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 REAL 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.```

LDE

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

F

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

WORK

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

INFO

```          INFO is INTEGER

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

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

On exit, normal operation:
== 0:  successful exit
< 0:  if INFO = -i, the i-th argument had an illegal value
> 0:  The equation is not solved correctly. One of the arising inner
system got singular.```
Optimizations:
• Replaced level-3 by level-2 and level-1 calls,
• Reorder the solution order to column first style,
• Replaced SAXPY operation by Fortran intrinsic
• Replaced all BLAS calls except of STRMV and SGER by Fortran intrinsic,
• Align local copies and fix problem size to <= 64
Author

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 226 of file sla_tgcsylv_l2_opt_local_copy_64.f90.

### subroutine sla_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:

``` SLA_TGCSYLV_L2_LOCAL_COPY_96 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 SGGES from LAPACK.```
Remarks

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

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 REAL, allowed values: +/-1
Specifies the sign between in the first equation.```

SGN2

```          SGN2 is REAL, 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 REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

```          C is REAL 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 REAL 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 REAL 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.```

LDE

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

F

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

WORK

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

INFO

```          INFO is INTEGER

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

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

On exit, normal operation:
== 0:  successful exit
< 0:  if INFO = -i, the i-th argument had an illegal value
> 0:  The equation is not solved correctly. One of the arising inner
system got singular.```
Optimizations:
• Replaced level-3 by level-2 and level-1 calls,
• Reorder the solution order to column first style,
• Replaced SAXPY operation by Fortran intrinsic
• Replaced all BLAS calls except of STRMV and SGER by Fortran intrinsic,
• Align local copies and fix problem size to <= 96
Author

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 225 of file sla_tgcsylv_l2_opt_local_copy_96.f90.

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

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

Purpose:

``` SLA_TGCSYLV_L2_REORDER 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 SGGES from LAPACK.```
Remarks

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

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 REAL, allowed values: +/-1
Specifies the sign between in the first equation.```

SGN2

```          SGN2 is REAL, allowed values: +/-1
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 REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

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

WORK

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

INFO

```          INFO is INTEGER

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

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

On exit, normal operation:
== 0:  successful exit
< 0:  if INFO = -i, the i-th argument had an illegal value
> 0:  The equation is not solved correctly. One of the arising inner
system got singular.```
Optimizations:
• Replaced level-3 by level-2 and level-1 calls,
• Reorder the solution order to column first style,
• Replaced SAXPY operation by Fortran intrinsic
• Replaced all BLAS calls except of STRMV and SGER by Fortran intrinsic,
• Column first computation of the solution.
Author

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 226 of file sla_tgcsylv_l2_opt_reorder.f90.

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

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

Purpose:

``` SLA_TGCSYLV_L2_UNOPT 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 SGGES from LAPACK.```
Attention

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

Parameters

TRANSA

```          TRANSA is CHARACTER(1)
Specifies the form of the system of equations with respect to A and C:
== 'N':  op1(A) = A, op1(C) = C (No transpose for A and C)
== 'T':  op1(A) = A**T, op1(C) = C **T (Transpose A and C)```

TRANSB

```          TRANSB is CHARACTER(1)
Specifies the form of the system of equations with respect to B and D:
== 'N':  op2(B) = B, op2(D) = D (No transpose for B and D)
== 'T':  op2(B) = B**T, op2(D) = D **T (Transpose B and D)```

SGN1

```          SGN1 is REAL, allowed values: +/-1
Specifies the sign between in the first equation.```

SGN2

```          SGN2 is REAL, allowed values: +/-1
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 REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

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

WORK

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

INFO

```          INFO is INTEGER

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

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

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

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 225 of file sla_tgcsylv_l2_unopt.f90.

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

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

Purpose:

``` SLA_TGCSYLV_L3 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 SGGES from LAPACK.```
Attention

The algorithm is implemented using BLAS level 3 operations.

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 REAL, allowed values: +/-1
Specifies the sign between in the first equation.```

SGN2

```          SGN2 is REAL, allowed values: +/-1
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 REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

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

WORK

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

INFO

```          INFO is INTEGER

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

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

On exit, normal operation:
== 0:  successful exit
< 0:  if INFO = -i, the i-th argument had an illegal value
> 0:  The equation is not solved correctly. One of the arising inner
system got singular.```
Optimizations:
• Nothing but standard L3 BLAS calls.
• Changed computation order to column major
Author

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 223 of file sla_tgcsylv_l3_opt.f90.

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

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

Purpose:

``` SLA_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 SGGES from LAPACK.```
Remarks

This algorithms uses 2 Stage blocking with DAG scheduled solves.

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 REAL, allowed values: +/-1
Specifies the sign between in the first equation.```

SGN2

```          SGN2 is REAL, allowed values: +/-1
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 REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

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

WORK

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

INFO

```          INFO is INTEGER

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

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

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

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 217 of file sla_tgcsylv_l3_2stage.f90.

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

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

Purpose:

``` SLA_TGCSYLV_L3_UNOPT 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 SGGES from LAPACK.```
Attention

The algorithm is implemented using BLAS level 3 operations.

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 REAL, allowed values: +/-1
Specifies the sign between in the first equation.```

SGN2

```          SGN2 is REAL, allowed values: +/-1
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 REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

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

WORK

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

INFO

```          INFO is INTEGER

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

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

On exit, normal operation:
== 0:  successful exit
< 0:  if INFO = -i, the i-th argument had an illegal value
> 0:  The equation is not solved correctly. One of the arising inner
system got singular.```
Optimizations:
• Nothing but standard L3 BLAS calls.
Author

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 223 of file sla_tgcsylv_l3_unopt.f90.

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

Recursive Algorithm for the generalized coupled Sylvester equation.

Purpose:

``` SLA_TGCSYLV_RECURSIVE 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 SGGES from LAPACK.```
Attention

The algorithm is implemented using BLAS level 3 operations.

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 REAL, allowed values: +/-1
Specifies the sign between in the first equation.```

SGN2

```          SGN2 is REAL, allowed values: +/-1
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 REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

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

WORK

```          WORK is REAL array, dimension 1
Workspace for the algorithm```

INFO

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

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 205 of file sla_tgcsylv_recursive.f90.

### subroutine sla_tgsylv_dag (character, dimension(1) transa, character, dimension(1) transb, real sgn, integer m, integer n, real, dimension(lda, m) a, integer lda, real, dimension(ldb, n) b, integer ldb, real, dimension(ldc, m) c, integer ldc, real, dimension(ldd, n) d, integer ldd, real, dimension(ldx, n) x, integer ldx, real scale, real, dimension(*) work, integer info)

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

Purpose:

``` SLA_TGSYLV_DAG solves a generalized Sylvester equation of the following forms

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)

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 SGGES from LAPACK.```
Remarks

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

Attention

Due to the parallel nature of the algorithm the scaling is not applied to the right hand. If the problem is ill-posed and scaling appears you have to solve the equation again with a solver with complete scaling support like SLA_TGSYLV_L3.

Parameters

TRANSA

```          TRANSA is CHARACTER(1)
Specifies the form of the system of equations with respect to A and C:
== 'N':  op1(A) = A, op1(C) = C (No transpose for A and C)
== 'T':  op1(A) = A**T, op1(C) = C **T (Transpose A and C)```

TRANSB

```          TRANSB is CHARACTER(1)
Specifies the form of the system of equations with respect to B and D:
== 'N':  op2(B) = B, op2(D) = D (No transpose for B and D)
== 'T':  op2(B) = B**T, op2(D) = D **T (Transpose B and D)```

SGN

```          SGN is REAL, allowed values: +/-1
Specifies the sign between the two parts of the Sylvester equation.
= 1 :  Solve Equation (1)
== -1:  Solve Equation (2)```

M

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

N

```          N is INTEGER
The order of the matrices B and D.  N >= 0.```

A

```          A is REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

```          C is REAL 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 REAL 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 REAL array, dimension (LDX,N)
On input, the matrix X contains the right hand side Y.
On output, the matrix X contains the solution of Equation (1) or (2)
as selected by TRANSA, TRANSB, and SGN.
Right hand side Y and the solution X are M-by-N matrices.```

LDX

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

SCALE

```          SCALE is REAL
SCALE is a scaling factor to prevent the overflow in the result.
If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems
could not be solved correctly, 0 < SCALE <= 1 holds true.
If SCALE .NE. 1 the problem is no solved correctly in this case
one have to use an other solver.```

WORK

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

INFO

```          INFO is INTEGER

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

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

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

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 208 of file sla_tgsylv_dag.f90.

### subroutine sla_tgsylv_gardiner_laub (character, dimension(1) transa, character, dimension(1) transb, real sgn, integer m, integer n, real, dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d, integer ldd, real, dimension(ldx, *) x, integer ldx, real scale, real, dimension(2*m,*) work, integer info)

Level-2 Gardiner-Laub algorithm for the generalized Sylvester equation.

Purpose:

``` SLA_TGSYLV_GARDINER_LAUB solves a generalized Sylvester equation of the following forms

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)

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 SGGES from LAPACK.```
Remarks

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

Attention

Scaling is not supported by this algorithm so SCALE == 1.

Parameters

TRANSA

```          TRANSA is CHARACTER(1)
Specifies the form of the system of equations with respect to A and C:
== 'N':  op1(A) = A, op1(C) = C (No transpose for A and C)
== 'T':  op1(A) = A**T, op1(C) = C **T (Transpose A and C)```

TRANSB

```          TRANSB is CHARACTER(1)
Specifies the form of the system of equations with respect to B and D:
== 'N':  op2(B) = B, op2(D) = D (No transpose for B and D)
== 'T':  op2(B) = B**T, op2(D) = D **T (Transpose B and D)```

SGN

```          SGN is REAL, allowed values: +/-1
Specifies the sign between the two parts of the Sylvester equation.
= 1 :  Solve Equation (1)
== -1:  Solve Equation (2)```

M

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

N

```          N is INTEGER
The order of the matrices B and D.  N >= 0.```

A

```          A is REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

```          C is REAL 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 REAL 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 REAL array, dimension (LDX,N)
On input, the matrix X contains the right hand side Y.
On output, the matrix X contains the solution of Equation (1) or (2)
as selected by TRANSA, TRANSB, and SGN.
Right hand side Y and the solution X are M-by-N matrices.```

LDX

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

SCALE

```          SCALE is REAL
SCALE is a scaling factor to prevent the overflow in the result.
If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems
could not be solved correctly, 0 < SCALE <= 1 holds true.```

WORK

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

INFO

```          INFO is INTEGER

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

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

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

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 204 of file sla_tgsylv_gardiner_laub.f90.

### subroutine sla_tgsylv_l2 (character, dimension(1) transa, character, dimension(1) transb, real sgn, integer m, integer n, real, dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d, integer ldd, real, dimension(ldx, *) x, integer ldx, real scale, real, dimension(*) work, integer info)

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

Purpose:

``` SLA_TGSYLV_L2 solves a generalized Sylvester equation of the following forms

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)

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 SGGES from LAPACK.```
Remarks

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

Parameters

TRANSA

```          TRANSA is CHARACTER(1)
Specifies the form of the system of equations with respect to A and C:
== 'N':  op1(A) = A, op1(C) = C (No transpose for A and C)
== 'T':  op1(A) = A**T, op1(C) = C **T (Transpose A and C)```

TRANSB

```          TRANSB is CHARACTER(1)
Specifies the form of the system of equations with respect to B and D:
== 'N':  op2(B) = B, op2(D) = D (No transpose for B and D)
== 'T':  op2(B) = B**T, op2(D) = D **T (Transpose B and D)```

SGN

```          SGN is REAL, allowed values: +/-1
Specifies the sign between the two parts of the Sylvester equation.
= 1 :  Solve Equation (1)
== -1:  Solve Equation (2)```

M

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

N

```          N is INTEGER
The order of the matrices B and D.  N >= 0.```

A

```          A is REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

```          C is REAL 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 REAL 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 REAL array, dimension (LDX,N)
On input, the matrix X contains the right hand side Y.
On output, the matrix X contains the solution of Equation (1) or (2)
as selected by TRANSA, TRANSB, and SGN.
Right hand side Y and the solution X are M-by-N matrices.```

LDX

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

SCALE

```          SCALE is REAL
SCALE is a scaling factor to prevent the overflow in the result.
If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems
could not be solved correctly, 0 < SCALE <= 1 holds true.```

WORK

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

INFO

```          INFO is INTEGER

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

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

On exit, normal operation:
== 0:  successful exit
< 0:  if INFO = -i, the i-th argument had an illegal value
> 0:  The equation is not solved correctly. One of the arising inner
system got singular.```
Optimizations:
• Replaced level-3 by level-2 and level-1 calls,
Author

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 207 of file sla_tgsylv_l2.f90.

### subroutine sla_tgsylv_l2_colwise (character, dimension(1) transa, character, dimension(1) transb, real sgn, integer m, integer n, real, dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d, integer ldd, real, dimension(ldx, *) x, integer ldx, real scale, real, dimension(*) work, integer info)

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

Purpose:

``` SLA_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)

or

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 SGGES from LAPACK.```
Attention

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

Parameters

TRANSA

```          TRANSA is CHARACTER(1)
Specifies the form of the system of equations with respect to A and C:
== 'N':  op1(A) = A, op1(C) = C (No transpose for A and C)
== 'T':  op1(A) = A**T, op1(C) = C **T (Transpose A and C)```

TRANSB

```          TRANSB is CHARACTER(1)
Specifies the form of the system of equations with respect to B and D:
== 'N':  op2(B) = B, op2(D) = D (No transpose for B and D)
== 'T':  op2(B) = B**T, op2(D) = D **T (Transpose B and D)```

SGN

```          SGN is REAL, allowed values: +/-1
Specifies the sign between the two parts of the Sylvester equation.
= 1 :  Solve Equation (1)
== -1:  Solve Equation (2)```

M

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

N

```          N is INTEGER
The order of the matrices B and D.  N >= 0.```

A

```          A is REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. 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 REAL 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 REAL 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 REAL array, dimension (LDX,N)
On input, the matrix X contains the right hand side Y.
On output, the matrix X contains the solution of Equation (1) or (2)
as selected by TRANSA, TRANSB, and SGN.
Right hand side Y and the solution X are M-by-N matrices.```

LDX

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

SCALE

```          SCALE is REAL
SCALE is a scaling factor to prevent the overflow in the result.
If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems
could not be solved correctly, 0 < SCALE <= 1 holds true.```

WORK

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

INFO

```          INFO is INTEGER

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

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

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

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 208 of file sla_tgsylv_l2_colwise.f90.

### subroutine sla_tgsylv_l2_local_copy (character, dimension(1) transa, character, dimension(1) transb, real sgn, integer m, integer n, real, dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d, integer ldd, real, dimension(ldx, *) x, integer ldx, real scale, real, dimension(*) work, integer info)

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

Purpose:

``` SLA_TGSYLV_L2_LOCAL_COPY solves a generalized Sylvester equation of the following forms

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)

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 SGGES from LAPACK.```
Attention

The algorithm solves problems of dimension at most 128 ( M,N <=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)```

SGN

```          SGN is REAL, allowed values: +/-1
Specifies the sign between the two parts of the Sylvester equation.
= 1 :  Solve Equation (1)
== -1:  Solve Equation (2)```

M

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

N

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

A

```          A is REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. 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 REAL 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 REAL 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 REAL array, dimension (LDX,N)
On input, the matrix X contains the right hand side Y.
On output, the matrix X contains the solution of Equation (1) or (2)
as selected by TRANSA, TRANSB, and SGN.
Right hand side Y and the solution X are M-by-N matrices.```

LDX

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

SCALE

```          SCALE is REAL
SCALE is a scaling factor to prevent the overflow in the result.
If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems
could not be solved correctly, 0 < SCALE <= 1 holds true.```

WORK

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

INFO

```          INFO is INTEGER

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

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

On exit, normal operation:
== 0:  successful exit
< 0:  if INFO = -i, the i-th argument had an illegal value
> 0:  The equation is not solved correctly. One of the arising inner
system got singular.```
Optimizations:
• Replaced level-3 by level-2 and level-1 calls,
• Reorder the solution order to column first style,
• Replaced SAXPY operation by Fortran intrinsic
• Replaced all BLAS calls except of STRMV and SGER by Fortran intrinsic,
• Use local copies of A, B, C, D, and X (M, N <=128) .
• STRMV works sequential.
Author

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 214 of file sla_tgsylv_l2_opt_local_copy.f90.

### subroutine sla_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:

``` SLA_TGSYLV_L2_LOCAL_COPY_128 solves a generalized Sylvester equation of the following forms

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)

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 SGGES from LAPACK.```
Attention

The size of the Problem is limited by M,N <= 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)```

SGN

```          SGN is REAL, allowed values: +/-1
Specifies the sign between the two parts of the Sylvester equation.
= 1 :  Solve Equation (1)
== -1:  Solve Equation (2)```

M

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

N

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

A

```          A is REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. 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 REAL 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 REAL 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 REAL array, dimension (LDX,N)
On input, the matrix X contains the right hand side Y.
On output, the matrix X contains the solution of Equation (1) or (2)
as selected by TRANSA, TRANSB, and SGN.
Right hand side Y and the solution X are M-by-N matrices.```

LDX

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

SCALE

```          SCALE is REAL
SCALE is a scaling factor to prevent the overflow in the result.
If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems
could not be solved correctly, 0 < SCALE <= 1 holds true.```

WORK

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

INFO

```          INFO is INTEGER

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

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

On exit, normal operation:
== 0:  successful exit
< 0:  if INFO = -i, the i-th argument had an illegal value
> 0:  The equation is not solved correctly. One of the arising inner
system got singular.```
Optimizations:
• Replaced level-3 by level-2 and level-1 calls,
• Reorder the solution order to column first style,
• Replaced SAXPY operation by Fortran intrinsic
• Replaced all BLAS calls except of STRMV by Fortran intrinsic,
• Use local copies of A, B, C, D, and X.
• STRMV works purely sequential
• Align local copies and fix problem size to <= 128
Author

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 215 of file sla_tgsylv_l2_opt_local_copy_128.f90.

### subroutine sla_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:

``` SLA_TGSYLV_L2_LOCAL_COPY_32 solves a generalized Sylvester equation of the following forms

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)

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 SGGES from LAPACK.```
Attention

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

Parameters

TRANSA

```          TRANSA is CHARACTER(1)
Specifies the form of the system of equations with respect to A and C:
== 'N':  op1(A) = A, op1(C) = C (No transpose for A and C)
== 'T':  op1(A) = A**T, op1(C) = C **T (Transpose A and C)```

TRANSB

```          TRANSB is CHARACTER(1)
Specifies the form of the system of equations with respect to B and D:
== 'N':  op2(B) = B, op2(D) = D (No transpose for B and D)
== 'T':  op2(B) = B**T, op2(D) = D **T (Transpose B and D)```

SGN

```          SGN is REAL, allowed values: +/-1
Specifies the sign between the two parts of the Sylvester equation.
= 1 :  Solve Equation (1)
== -1:  Solve Equation (2)```

M

```          M is INTEGER
The order of the matrices A and C.  32 >= M >= 0.```

N

```          N is INTEGER
The order of the matrices B and D.  32 >= N >= 0.```

A

```          A is REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. 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 REAL 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 REAL 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 REAL array, dimension (LDX,N)
On input, the matrix X contains the right hand side Y.
On output, the matrix X contains the solution of Equation (1) or (2)
as selected by TRANSA, TRANSB, and SGN.
Right hand side Y and the solution X are M-by-N matrices.```

LDX

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

SCALE

```          SCALE is REAL
SCALE is a scaling factor to prevent the overflow in the result.
If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems
could not be solved correctly, 0 < SCALE <= 1 holds true.```

WORK

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

INFO

```          INFO is INTEGER

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

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

On exit, normal operation:
== 0:  successful exit
< 0:  if INFO = -i, the i-th argument had an illegal value
> 0:  The equation is not solved correctly. One of the arising inner
system got singular.```
Optimizations:
• Replaced level-3 by level-2 and level-1 calls,
• Reorder the solution order to column first style,
• Replaced SAXPY operation by Fortran intrinsic
• Replaced all BLAS calls except of STRMV by Fortran intrinsic,
• Use local copies of A, B, C, D, and X.
• STRMV works purely sequential
• Align local copies and fix problem size to <= 32
Author

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 215 of file sla_tgsylv_l2_opt_local_copy_32.f90.

### subroutine sla_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:

``` SLA_TGSYLV_L2_LOCAL_COPY_64 solves a generalized Sylvester equation of the following forms

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)

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 SGGES from LAPACK.```
Attention

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

Parameters

TRANSA

```          TRANSA is CHARACTER(1)
Specifies the form of the system of equations with respect to A and C:
== 'N':  op1(A) = A, op1(C) = C (No transpose for A and C)
== 'T':  op1(A) = A**T, op1(C) = C **T (Transpose A and C)```

TRANSB

```          TRANSB is CHARACTER(1)
Specifies the form of the system of equations with respect to B and D:
== 'N':  op2(B) = B, op2(D) = D (No transpose for B and D)
== 'T':  op2(B) = B**T, op2(D) = D **T (Transpose B and D)```

SGN

```          SGN is REAL, allowed values: +/-1
Specifies the sign between the two parts of the Sylvester equation.
= 1 :  Solve Equation (1)
== -1:  Solve Equation (2)```

M

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

N

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

A

```          A is REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. 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 REAL 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 REAL 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 REAL array, dimension (LDX,N)
On input, the matrix X contains the right hand side Y.
On output, the matrix X contains the solution of Equation (1) or (2)
as selected by TRANSA, TRANSB, and SGN.
Right hand side Y and the solution X are M-by-N matrices.```

LDX

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

SCALE

```          SCALE is REAL
SCALE is a scaling factor to prevent the overflow in the result.
If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems
could not be solved correctly, 0 < SCALE <= 1 holds true.```

WORK

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

INFO

```          INFO is INTEGER

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

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

On exit, normal operation:
== 0:  successful exit
< 0:  if INFO = -i, the i-th argument had an illegal value
> 0:  The equation is not solved correctly. One of the arising inner
system got singular.```
Optimizations:
• Replaced level-3 by level-2 and level-1 calls,
• Reorder the solution order to column first style,
• Replaced SAXPY operation by Fortran intrinsic
• Replaced all BLAS calls except of STRMV by Fortran intrinsic,
• Use local copies of A, B, C, D, and X.
• STRMV works purely sequential
• Align local copies and fix problem size to <= 64
Author

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 215 of file sla_tgsylv_l2_opt_local_copy_64.f90.

### subroutine sla_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:

``` SLA_TGSYLV_L2_LOCAL_COPY_96 solves a generalized Sylvester equation of the following forms

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)

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 SGGES from LAPACK.```
Attention

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

Parameters

TRANSA

```          TRANSA is CHARACTER(1)
Specifies the form of the system of equations with respect to A and C:
== 'N':  op1(A) = A, op1(C) = C (No transpose for A and C)
== 'T':  op1(A) = A**T, op1(C) = C **T (Transpose A and C)```

TRANSB

```          TRANSB is CHARACTER(1)
Specifies the form of the system of equations with respect to B and D:
== 'N':  op2(B) = B, op2(D) = D (No transpose for B and D)
== 'T':  op2(B) = B**T, op2(D) = D **T (Transpose B and D)```

SGN

```          SGN is REAL, allowed values: +/-1
Specifies the sign between the two parts of the Sylvester equation.
= 1 :  Solve Equation (1)
== -1:  Solve Equation (2)```

M

```          M is INTEGER
The order of the matrices A and C.  96 >= M >= 0.```

N

```          N is INTEGER
The order of the matrices B and D.  96 >= N >= 0.```

A

```          A is REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. 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 REAL 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 REAL 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 REAL array, dimension (LDX,N)
On input, the matrix X contains the right hand side Y.
On output, the matrix X contains the solution of Equation (1) or (2)
as selected by TRANSA, TRANSB, and SGN.
Right hand side Y and the solution X are M-by-N matrices.```

LDX

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

SCALE

```          SCALE is REAL
SCALE is a scaling factor to prevent the overflow in the result.
If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems
could not be solved correctly, 0 < SCALE <= 1 holds true.```

WORK

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

INFO

```          INFO is INTEGER

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

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

On exit, normal operation:
== 0:  successful exit
< 0:  if INFO = -i, the i-th argument had an illegal value
> 0:  The equation is not solved correctly. One of the arising inner
system got singular.```
Optimizations:
• Replaced level-3 by level-2 and level-1 calls,
• Reorder the solution order to column first style,
• Replaced SAXPY operation by Fortran intrinsic
• Replaced all BLAS calls except of STRMV by Fortran intrinsic,
• Use local copies of A, B, C, D, and X.
• STRMV works purely sequential
• Align local copies and fix problem size to <= 96
Author

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 215 of file sla_tgsylv_l2_opt_local_copy_96.f90.

### subroutine sla_tgsylv_l2_reorder (character, dimension(1) transa, character, dimension(1) transb, real sgn, integer m, integer n, real, dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d, integer ldd, real, dimension(ldx, *) x, integer ldx, real scale, real, dimension(*) work, integer info)

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

Purpose:

``` SLA_TGSYLV_L2_REORDER solves a generalized Sylvester equation of the following forms

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)

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 SGGES from LAPACK.```
Attention

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

Parameters

TRANSA

```          TRANSA is CHARACTER(1)
Specifies the form of the system of equations with respect to A and C:
== 'N':  op1(A) = A, op1(C) = C (No transpose for A and C)
== 'T':  op1(A) = A**T, op1(C) = C **T (Transpose A and C)```

TRANSB

```          TRANSB is CHARACTER(1)
Specifies the form of the system of equations with respect to B and D:
== 'N':  op2(B) = B, op2(D) = D (No transpose for B and D)
== 'T':  op2(B) = B**T, op2(D) = D **T (Transpose B and D)```

SGN

```          SGN is REAL, allowed values: +/-1
Specifies the sign between the two parts of the Sylvester equation.
= 1 :  Solve Equation (1)
== -1:  Solve Equation (2)```

M

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

N

```          N is INTEGER
The order of the matrices B and D.  N >= 0.```

A

```          A is REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

```          C is REAL 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 REAL 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 REAL array, dimension (LDX,N)
On input, the matrix X contains the right hand side Y.
On output, the matrix X contains the solution of Equation (1) or (2)
as selected by TRANSA, TRANSB, and SGN.
Right hand side Y and the solution X are M-by-N matrices.```

LDX

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

SCALE

```          SCALE is REAL
SCALE is a scaling factor to prevent the overflow in the result.
If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems
could not be solved correctly, 0 < SCALE <= 1 holds true.```

WORK

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

INFO

```          INFO is INTEGER

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

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

On exit, normal operation:
== 0:  successful exit
< 0:  if INFO = -i, the i-th argument had an illegal value
> 0:  The equation is not solved correctly. One of the arising inner
system got singular.```
Optimizations:
• Replaced level-3 by level-2 and level-1 calls,
• Reorder the solution order to column first style,
• Replaced SAXPY operation by Fortran intrinsic
• Replaced all BLAS calls except of STRMV and SGER by Fortran intrinsic,
Author

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 211 of file sla_tgsylv_l2_reorder.f90.

### subroutine sla_tgsylv_l2_rowwise (character, dimension(1) transa, character, dimension(1) transb, real sgn, integer m, integer n, real, dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d, integer ldd, real, dimension(ldx, *) x, integer ldx, real scale, real, dimension(*) work, integer info)

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

Purpose:

``` SLA_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)

or

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 SGGES from LAPACK.```
Attention

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

Parameters

TRANSA

```          TRANSA is CHARACTER(1)
Specifies the form of the system of equations with respect to A and C:
== 'N':  op1(A) = A, op1(C) = C (No transpose for A and C)
== 'T':  op1(A) = A**T, op1(C) = C **T (Transpose A and C)```

TRANSB

```          TRANSB is CHARACTER(1)
Specifies the form of the system of equations with respect to B and D:
== 'N':  op2(B) = B, op2(D) = D (No transpose for B and D)
== 'T':  op2(B) = B**T, op2(D) = D **T (Transpose B and D)```

SGN

```          SGN is REAL, allowed values: +/-1
Specifies the sign between the two parts of the Sylvester equation.
= 1 :  Solve Equation (1)
== -1:  Solve Equation (2)```

M

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

N

```          N is INTEGER
The order of the matrices B and D.  N >= 0.```

A

```          A is REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. 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 REAL 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 REAL 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 REAL array, dimension (LDX,N)
On input, the matrix X contains the right hand side Y.
On output, the matrix X contains the solution of Equation (1) or (2)
as selected by TRANSA, TRANSB, and SGN.
Right hand side Y and the solution X are M-by-N matrices.```

LDX

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

SCALE

```          SCALE is REAL
SCALE is a scaling factor to prevent the overflow in the result.
If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems
could not be solved correctly, 0 < SCALE <= 1 holds true.```

WORK

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

INFO

```          INFO is INTEGER

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

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

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

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 208 of file sla_tgsylv_l2_rowwise.f90.

### subroutine sla_tgsylv_l3 (character, dimension(1) transa, character, dimension(1) transb, real sgn, integer m, integer n, real, dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d, integer ldd, real, dimension(ldx, *) x, integer ldx, real scale, real, dimension(*) work, integer info)

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

Purpose:

``` SLA_TGSYLV_L3 solves a generalized Sylvester equation of the following forms

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)

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 SGGES from LAPACK.```
Attention

The algorithm is implemented using BLAS level 3 operations.

Remarks

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

Parameters

TRANSA

```          TRANSA is CHARACTER(1)
Specifies the form of the system of equations with respect to A and C:
== 'N':  op1(A) = A, op1(C) = C (No transpose for A and C)
== 'T':  op1(A) = A**T, op1(C) = C **T (Transpose A and C)```

TRANSB

```          TRANSB is CHARACTER(1)
Specifies the form of the system of equations with respect to B and D:
== 'N':  op2(B) = B, op2(D) = D (No transpose for B and D)
== 'T':  op2(B) = B**T, op2(D) = D **T (Transpose B and D)```

SGN

```          SGN is REAL, allowed values: +/-1
Specifies the sign between the two parts of the Sylvester equation.
= 1 :  Solve Equation (1)
== -1:  Solve Equation (2)```

M

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

N

```          N is INTEGER
The order of the matrices B and D.  N >= 0.```

A

```          A is REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

```          C is REAL 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 REAL 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 REAL array, dimension (LDX,N)
On input, the matrix X contains the right hand side Y.
On output, the matrix X contains the solution of Equation (1) or (2)
as selected by TRANSA, TRANSB, and SGN.
Right hand side Y and the solution X are M-by-N matrices.```

LDX

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

SCALE

```          SCALE is REAL
SCALE is a scaling factor to prevent the overflow in the result.
If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems
could not be solved correctly, 0 < SCALE <= 1 holds true.```

WORK

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

INFO

```          INFO is INTEGER

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

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

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

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 205 of file sla_tgsylv_l3_opt.f90.

### subroutine sla_tgsylv_l3_2s (character, dimension(1) transa, character, dimension(1) transb, real sgn, integer m, integer n, real, dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d, integer ldd, real, dimension(ldx, *) x, integer ldx, real scale, real, dimension(*) work, integer info)

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

Purpose:

``` SLA_TGSYLV_L3_2S solves a generalized Sylvester equation of the following forms

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)

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 SGGES from LAPACK.```
Attention

The algorithm is implemented using BLAS level 3 operations.

Parameters

TRANSA

```          TRANSA is CHARACTER(1)
Specifies the form of the system of equations with respect to A and C:
== 'N':  op1(A) = A, op1(C) = C (No transpose for A and C)
== 'T':  op1(A) = A**T, op1(C) = C **T (Transpose A and C)```

TRANSB

```          TRANSB is CHARACTER(1)
Specifies the form of the system of equations with respect to B and D:
== 'N':  op2(B) = B, op2(D) = D (No transpose for B and D)
== 'T':  op2(B) = B**T, op2(D) = D **T (Transpose B and D)```

SGN

```          SGN is REAL, allowed values: +/-1
Specifies the sign between the two parts of the Sylvester equation.
= 1 :  Solve Equation (1)
== -1:  Solve Equation (2)```

M

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

N

```          N is INTEGER
The order of the matrices B and D.  N >= 0.```

A

```          A is REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

```          C is REAL 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 REAL 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 REAL array, dimension (LDX,N)
On input, the matrix X contains the right hand side Y.
On output, the matrix X contains the solution of Equation (1) or (2)
as selected by TRANSA, TRANSB, and SGN.
Right hand side Y and the solution X are M-by-N matrices.```

LDX

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

SCALE

```          SCALE is REAL
SCALE is a scaling factor to prevent the overflow in the result.
If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems
could not be solved correctly, 0 < SCALE <= 1 holds true.```

WORK

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

INFO

```          INFO is INTEGER

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

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

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

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 202 of file sla_tgsylv_l3_2stage.f90.

### subroutine sla_tgsylv_l3_colwise (character, dimension(1) transa, character, dimension(1) transb, real sgn, integer m, integer n, real, dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d, integer ldd, real, dimension(ldx, *) x, integer ldx, real scale, real, dimension(*) work, integer info)

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

Purpose:

``` SLA_TGSYLV_L3_COLWISE solves a generalized Sylvester equation of the following forms

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)

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 SGGES from LAPACK.```
Attention

The algorithm is implemented using BLAS level 3 operations.

Remarks

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

Parameters

TRANSA

```          TRANSA is CHARACTER(1)
Specifies the form of the system of equations with respect to A and C:
== 'N':  op1(A) = A, op1(C) = C (No transpose for A and C)
== 'T':  op1(A) = A**T, op1(C) = C **T (Transpose A and C)```

TRANSB

```          TRANSB is CHARACTER(1)
Specifies the form of the system of equations with respect to B and D:
== 'N':  op2(B) = B, op2(D) = D (No transpose for B and D)
== 'T':  op2(B) = B**T, op2(D) = D **T (Transpose B and D)```

SGN

```          SGN is REAL, allowed values: +/-1
Specifies the sign between the two parts of the Sylvester equation.
= 1 :  Solve Equation (1)
== -1:  Solve Equation (2)```

M

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

N

```          N is INTEGER
The order of the matrices B and D.  N >= 0.```

A

```          A is REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

```          C is REAL 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 REAL 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 REAL array, dimension (LDX,N)
On input, the matrix X contains the right hand side Y.
On output, the matrix X contains the solution of Equation (1) or (2)
as selected by TRANSA, TRANSB, and SGN.
Right hand side Y and the solution X are M-by-N matrices.```

LDX

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

SCALE

```          SCALE is REAL
SCALE is a scaling factor to prevent the overflow in the result.
If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems
could not be solved correctly, 0 < SCALE <= 1 holds true.```

WORK

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

INFO

```          INFO is INTEGER

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

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

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

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 204 of file sla_tgsylv_l3_colwise.f90.

### subroutine sla_tgsylv_l3_rowwise (character, dimension(1) transa, character, dimension(1) transb, real sgn, integer m, integer n, real, dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d, integer ldd, real, dimension(ldx, *) x, integer ldx, real scale, real, dimension(*) work, integer info)

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

Purpose:

``` SLA_TGSYLV_L3_ROWWISE solves a generalized Sylvester equation of the following forms

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)

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 SGGES from LAPACK.```
Attention

The algorithm is implemented using BLAS level 3 operations.

Parameters

TRANSA

```          TRANSA is CHARACTER(1)
Specifies the form of the system of equations with respect to A and C:
== 'N':  op1(A) = A, op1(C) = C (No transpose for A and C)
== 'T':  op1(A) = A**T, op1(C) = C **T (Transpose A and C)```

TRANSB

```          TRANSB is CHARACTER(1)
Specifies the form of the system of equations with respect to B and D:
== 'N':  op2(B) = B, op2(D) = D (No transpose for B and D)
== 'T':  op2(B) = B**T, op2(D) = D **T (Transpose B and D)```

SGN

```          SGN is REAL, allowed values: +/-1
Specifies the sign between the two parts of the Sylvester equation.
= 1 :  Solve Equation (1)
== -1:  Solve Equation (2)```

M

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

N

```          N is INTEGER
The order of the matrices B and D.  N >= 0.```

A

```          A is REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

```          C is REAL 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 REAL 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 REAL array, dimension (LDX,N)
On input, the matrix X contains the right hand side Y.
On output, the matrix X contains the solution of Equation (1) or (2)
as selected by TRANSA, TRANSB, and SGN.
Right hand side Y and the solution X are M-by-N matrices.```

LDX

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

SCALE

```          SCALE is REAL
SCALE is a scaling factor to prevent the overflow in the result.
If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems
could not be solved correctly, 0 < SCALE <= 1 holds true.```

WORK

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

INFO

```          INFO is INTEGER

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

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

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

Martin Koehler, MPI Magdeburg

Date

Januar 2023

Definition at line 203 of file sla_tgsylv_l3_rowwise.f90.

### recursive subroutine sla_tgsylv_recursive (character, dimension(1) transa, character, dimension(1) transb, real sgn, integer m, integer n, real, dimension(lda, *) a, integer lda, real, dimension(ldb, *) b, integer ldb, real, dimension(ldc, *) c, integer ldc, real, dimension(ldd, *) d, integer ldd, real, dimension(ldx, *) x, integer ldx, real scale, real, dimension(*) work, integer info)

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

Purpose:

``` SLA_TGSYLV_RECURSIVE solves a generalized Sylvester equation of the following forms

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)

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
M-by-N matrices. Typically the matrix pencils (A,C) and (B,D) ( or (D,B) ) are
created by SGGES from LAPACK.```
Attention

The algorithm is a rewrite of RECSY with some optimizations.

Parameters

TRANSA

```          TRANSA is CHARACTER(1)
Specifies the form of the system of equations with respect to A and C:
== 'N':  op1(A) = A, op1(C) = C (No transpose for A and C)
== 'T':  op1(A) = A**T, op1(C) = C **T (Transpose A and C)```

TRANSB

```          TRANSB is CHARACTER(1)
Specifies the form of the system of equations with respect to B and D:
== 'N':  op2(B) = B, op2(D) = D (No transpose for B and D)
== 'T':  op2(B) = B**T, op2(D) = D **T (Transpose B and D)```

SGN

```          SGN is REAL, allowed values: +/-1
Specifies the sign between the two parts of the Sylvester equation.
= 1 :  Solve Equation (1)
== -1:  Solve Equation (2)```

M

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

N

```          N is INTEGER
The order of the matrices B and D.  N >= 0.```

A

```          A is REAL array, dimension (LDA,M)
The matrix A must be (quasi-) upper triangular.```

LDA

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

B

```          B is REAL array, dimension (LDB,N)
The matrix B must be (quasi-) upper triangular. If the matrix D is already
quasi-upper triangular the matrix B must be upper triangular.```

LDB

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

C

```          C is REAL 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 REAL 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 REAL 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

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

SCALE

```          SCALE is REAL
SCALE is a scaling factor to prevent the overflow in the result.
If INFO == 0 then SCALE is 1.0 otherwise if one of the inner systems
could not be solved correctly, 0 < SCALE <= 1 holds true.```

WORK

```          WORK is REAL array, dimension 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

Januar 2023

Definition at line 190 of file sla_tgsylv_recursive.f90.

## Author

Generated automatically by Doxygen for MEPACK from the source code.

## Info

Tue Mar 7 2023 Version 1.0.3 MEPACK