# dbltgsylv - Man Page

## Name

dbltgsylv — Double Precision

— Double Precision routines for triangular generalized Sylvester equations.

## Synopsis

### Functions

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

subroutine **dla_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 **dla_tgsylv_l2** (transa, transb, sgn, m, n, a, lda, b, ldb, c, ldc, d, ldd, x, ldx, scale, work, info)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Recursive blocking Algorithm for the dual generalized coupled Sylvester equation.

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

Recursive Algorithm for the generalized coupled Sylvester equation.

recursive subroutine **dla_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

Double Precision routines for triangular generalized Sylvester equations.

This section contains the solver 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 double precision arithmetic.

## Function Documentation

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

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

**Purpose:**

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

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

*SGN2*SGN2 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the second equation.

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

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

*E*E is DOUBLE PRECISION array, dimension (LDE,N) On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices.

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

*F*F is DOUBLE PRECISION array, dimension (LDF,N) On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices.

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

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

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

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

**Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **217** of file **dla_tgcsylv_dag.f90**.

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

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

**Purpose:**

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

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

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

*SGN1*SGN1 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the first equation.

*SGN2*SGN2 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the second equation.

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

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

*E*E is DOUBLE PRECISION array, dimension (LDE,N) On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices.

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

*F*F is DOUBLE PRECISION array, dimension (LDF,N) On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices.

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

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

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

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

**Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **237** of file **dla_tgcsylv_dual_dag.f90**.

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

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

**Purpose:**

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

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

**dla_tgcsylv_dual_recursive**.**Parameters***SGN1*SGN1 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the first equation.

*SGN2*SGN2 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the second equation.

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

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

*E*E is DOUBLE PRECISION array, dimension (LDE,N) On input, the matrix E contains the right hand side E. On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices.

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

*F*F is DOUBLE PRECISION array, dimension (LDF,N) On input, the matrix F contains the right hand side F. On output, the matrix F contains the solution L of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side F and the solution L are M-by-N matrices.

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

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

*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**January 2024

Definition at line **206** of file **dla_tgcsylv_dual_kernel_44_nn.f90**.

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

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

**Purpose:**

DLA_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 DGGES from LAPACK. The equation (1) is the dual to the generalized coupled Sylvester equation op1(A) * R + SGN1 * L * op2(B) = SCALE * E (2) op1(C) * R + SGN2 * L * op2(D) = SCALE * F The equation (1) is the dual one to equation (2) with respect to the underlying linear system. Let Z be the matrix formed by rewriting (2) into its Kronecker form. This yields | kron(I, op1(A)) SGN1*kron(op2(B)**T, I) | | Vec R | | Vec E | Z X = | |*| | = | | | kron(I, op1(C)) SGN2*kron(op2(D)**T, I) | | Vec L | | Vec F | Regarding Z**T one obtains | kron(I, op1(A)**T ) kron(I, op1(C)**T) | | Vec R | | Vec E | Z**T X = | |*| | = | | | SGN1*kron(op2(B), I) SGN2*kron(op2(D), I) | | Vec L | | Vec F | which belongs to the Sylvester equation (1). For this reason the parameters TRANSA and TRANSB are expressed in terms of the Sylvester equation (2).

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

**dla_tgcsylv_dual_recursive**.**Parameters***SGN1*SGN1 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the first equation.

*SGN2*SGN2 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the second equation.

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

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

*SCALE**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**January 2024

Definition at line **206** of file **dla_tgcsylv_dual_kernel_44_nt.f90**.

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

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

**Purpose:**

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

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

**dla_tgcsylv_dual_recursive**.**Parameters***SGN1*SGN1 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the first equation.

*SGN2*SGN2 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the second equation.

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

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

*SCALE**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**January 2024

Definition at line **206** of file **dla_tgcsylv_dual_kernel_44_tn.f90**.

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

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

**Purpose:**

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

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

**dla_tgcsylv_dual_recursive**.**Parameters***SGN1*SGN1 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the first equation.

*SGN2*SGN2 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the second equation.

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

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

*SCALE**INFO***Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **206** of file **dla_tgcsylv_dual_kernel_44_tt.f90**.

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

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

**Purpose:**

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

**Remarks**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 DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the first equation.

*SGN2*SGN2 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the second equation.

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

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

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

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

**Optimizations:**- Replaced DAXPY operation by Fortran intrinsic
- Replaced all BLAS calls except of DGER by Fortran intrinsic,

**Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **245** of file **dla_tgcsylv_dual_l2.f90**.

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

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

**Purpose:**

DLA_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 DGGES from LAPACK. The equation (1) is the dual to the generalized coupled Sylvester equation op1(A) * R + SGN1 * L * op2(B) = SCALE * E (2) op1(C) * R + SGN2 * L * op2(D) = SCALE * F The equation (1) is the dual one to equation (2) with respect to the underlying linear system. Let Z be the matrix formed by rewriting (2) into its Kronecker form. This yields | kron(I, op1(A)) SGN1*kron(op2(B)**T, I) | | Vec R | | Vec E | Z X = | |*| | = | | | kron(I, op1(C)) SGN2*kron(op2(D)**T, I) | | Vec L | | Vec F | Regarding Z**T one obtains | kron(I, op1(A)**T ) kron(I, op1(C)**T) | | Vec R | | Vec E | Z**T X = | |*| | = | | | SGN1*kron(op2(B), I) SGN2*kron(op2(D), I) | | Vec L | | Vec F | which belongs to the Sylvester equation (1). For this reason the parameters TRANSA and TRANSB are expressed in terms of the Sylvester equation (2).

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

**Parameters***TRANSA**TRANSB**SGN1*SGN1 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the first equation.

*SGN2*SGN2 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the second equation.

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

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

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

**Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **247** of file **dla_tgcsylv_dual_l2_opt_local_copy.f90**.

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

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

**Purpose:**

DLA_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 DGGES from LAPACK. The equation (1) is the dual to the generalized coupled Sylvester equation op1(A) * R + SGN1 * L * op2(B) = SCALE * E (2) op1(C) * R + SGN2 * L * op2(D) = SCALE * F The equation (1) is the dual one to equation (2) with respect to the underlying linear system. Let Z be the matrix formed by rewriting (2) into its Kronecker form. This yields | kron(I, op1(A)) SGN1*kron(op2(B)**T, I) | | Vec R | | Vec E | Z X = | |*| | = | | | kron(I, op1(C)) SGN2*kron(op2(D)**T, I) | | Vec L | | Vec F | Regarding Z**T one obtains | kron(I, op1(A)**T ) kron(I, op1(C)**T) | | Vec R | | Vec E | Z**T X = | |*| | = | | | SGN1*kron(op2(B), I) SGN2*kron(op2(D), I) | | Vec L | | Vec F | which belongs to the Sylvester equation (1). For this reason the parameters TRANSA and TRANSB are expressed in terms of the Sylvester equation (2).

**Remarks**The algorithm is implemented using BLAS level 2 operations and hand written optimizations and 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**TRANSB**SGN1*SGN1 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the first equation.

*SGN2*SGN2 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the second equation.

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

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

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

**Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **248** of file **dla_tgcsylv_dual_l2_opt_local_copy_128.f90**.

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

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

**Purpose:**

DLA_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 DGGES from LAPACK. The equation (1) is the dual to the generalized coupled Sylvester equation op1(A) * R + SGN1 * L * op2(B) = SCALE * E (2) op1(C) * R + SGN2 * L * op2(D) = SCALE * F The equation (1) is the dual one to equation (2) with respect to the underlying linear system. Let Z be the matrix formed by rewriting (2) into its Kronecker form. This yields | kron(I, op1(A)) SGN1*kron(op2(B)**T, I) | | Vec R | | Vec E | Z X = | |*| | = | | | kron(I, op1(C)) SGN2*kron(op2(D)**T, I) | | Vec L | | Vec F | Regarding Z**T one obtains | kron(I, op1(A)**T ) kron(I, op1(C)**T) | | Vec R | | Vec E | Z**T X = | |*| | = | | | SGN1*kron(op2(B), I) SGN2*kron(op2(D), I) | | Vec L | | Vec F | which belongs to the Sylvester equation (1). For this reason the parameters TRANSA and TRANSB are expressed in terms of the Sylvester equation (2).

**Remarks**The algorithm is implemented using BLAS level 2 operations and hand written optimizations and 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**TRANSB**SGN1*SGN1 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the first equation.

*SGN2*SGN2 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the second equation.

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

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

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

**Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **247** of file **dla_tgcsylv_dual_l2_opt_local_copy_32.f90**.

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

**Purpose:**

DLA_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 DGGES from LAPACK. The equation (1) is the dual to the generalized coupled Sylvester equation op1(A) * R + SGN1 * L * op2(B) = SCALE * E (2) op1(C) * R + SGN2 * L * op2(D) = SCALE * F The equation (1) is the dual one to equation (2) with respect to the underlying linear system. Let Z be the matrix formed by rewriting (2) into its Kronecker form. This yields | kron(I, op1(A)) SGN1*kron(op2(B)**T, I) | | Vec R | | Vec E | Z X = | |*| | = | | | kron(I, op1(C)) SGN2*kron(op2(D)**T, I) | | Vec L | | Vec F | Regarding Z**T one obtains | kron(I, op1(A)**T ) kron(I, op1(C)**T) | | Vec R | | Vec E | Z**T X = | |*| | = | | | SGN1*kron(op2(B), I) SGN2*kron(op2(D), I) | | Vec L | | Vec F | which belongs to the Sylvester equation (1). For this reason the parameters TRANSA and TRANSB are expressed in terms of the Sylvester equation (2).

**Remarks**The algorithm is implemented using BLAS level 2 operations and hand written optimizations and 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**TRANSB**SGN1*SGN1 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the first equation.

*SGN2*SGN2 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the second equation.

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

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

*SCALE**WORK**INFO***Optimizations:**- Replaced DAXPY operation by Fortran intrinsic
- No BLAS calls.
- The local copies are aligned to 64 byte if supported by the compiler.

**Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **246** of file **dla_tgcsylv_dual_l2_opt_local_copy_64.f90**.

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

**Purpose:**

DLA_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 DGGES from LAPACK. The equation (1) is the dual to the generalized coupled Sylvester equation op1(A) * R + SGN1 * L * op2(B) = SCALE * E (2) op1(C) * R + SGN2 * L * op2(D) = SCALE * F The equation (1) is the dual one to equation (2) with respect to the underlying linear system. Let Z be the matrix formed by rewriting (2) into its Kronecker form. This yields | kron(I, op1(A)) SGN1*kron(op2(B)**T, I) | | Vec R | | Vec E | Z X = | |*| | = | | | kron(I, op1(C)) SGN2*kron(op2(D)**T, I) | | Vec L | | Vec F | Regarding Z**T one obtains | kron(I, op1(A)**T ) kron(I, op1(C)**T) | | Vec R | | Vec E | Z**T X = | |*| | = | | | SGN1*kron(op2(B), I) SGN2*kron(op2(D), I) | | Vec L | | Vec F | which belongs to the Sylvester equation (1). For this reason the parameters TRANSA and TRANSB are expressed in terms of the Sylvester equation (2).

**Remarks**The algorithm is implemented using BLAS level 2 operations and hand written optimizations and 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**TRANSB**SGN1*SGN1 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the first equation.

*SGN2*SGN2 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the second equation.

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

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

*SCALE**WORK**INFO***Optimizations:**- Replaced DAXPY operation by Fortran intrinsic
- No BLAS calls.
- The local copies are aligned to 64 byte if supported by the compiler.

**Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **246** of file **dla_tgcsylv_dual_l2_opt_local_copy_96.f90**.

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

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

**Purpose:**

DLA_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 DGGES from LAPACK. The equation (1) is the dual to the generalized coupled Sylvester equation op1(A) * R + SGN1 * L * op2(B) = SCALE * E (2) op1(C) * R + SGN2 * L * op2(D) = SCALE * F The equation (1) is the dual one to equation (2) with respect to the underlying linear system. Let Z be the matrix formed by rewriting (2) into its Kronecker form. This yields | kron(I, op1(A)) SGN1*kron(op2(B)**T, I) | | Vec R | | Vec E | Z X = | |*| | = | | | kron(I, op1(C)) SGN2*kron(op2(D)**T, I) | | Vec L | | Vec F | Regarding Z**T one obtains | kron(I, op1(A)**T ) kron(I, op1(C)**T) | | Vec R | | Vec E | Z**T X = | |*| | = | | | SGN1*kron(op2(B), I) SGN2*kron(op2(D), I) | | Vec L | | Vec F | which belongs to the Sylvester equation (1). For this reason the parameters TRANSA and TRANSB are expressed in terms of the Sylvester equation (2).

**Parameters***TRANSA**TRANSB**SGN1*SGN1 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the first equation.

*SGN2*SGN2 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the second equation.

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

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

*SCALE**WORK**INFO***Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **236** of file **dla_tgcsylv_dual_l3.f90**.

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

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

**Purpose:**

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

**Parameters***TRANSA**TRANSB**SGN1*SGN1 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the first equation.

*SGN2*SGN2 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the second equation.

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

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

*SCALE**WORK**INFO***Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **237** of file **dla_tgcsylv_dual_l3_2stage.f90**.

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

Recursive blocking Algorithm for the dual generalized coupled Sylvester equation.

**Purpose:**

DLA_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 DGGES from LAPACK. The equation (1) is the dual to the generalized coupled Sylvester equation op1(A) * R + SGN1 * L * op2(B) = SCALE * E (2) op1(C) * R + SGN2 * L * op2(D) = SCALE * F The equation (1) is the dual one to equation (2) with respect to the underlying linear system. Let Z be the matrix formed by rewriting (2) into its Kronecker form. This yields | kron(I, op1(A)) SGN1*kron(op2(B)**T, I) | | Vec R | | Vec E | Z X = | |*| | = | | | kron(I, op1(C)) SGN2*kron(op2(D)**T, I) | | Vec L | | Vec F | Regarding Z**T one obtains | kron(I, op1(A)**T ) kron(I, op1(C)**T) | | Vec R | | Vec E | Z**T X = | |*| | = | | | SGN1*kron(op2(B), I) SGN2*kron(op2(D), I) | | Vec L | | Vec F | which belongs to the Sylvester equation (1). For this reason the parameters TRANSA and TRANSB are expressed in terms of the Sylvester equation (2).

**Parameters***TRANSA**TRANSB**SGN1*SGN1 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the first equation.

*SGN2*SGN2 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the second equation.

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

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

*SCALE**WORK**INFO***Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **240** of file **dla_tgcsylv_dual_recursive.f90**.

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

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

**Purpose:**

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

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

**Parameters***TRANSA**TRANSB**SGN1*SGN1 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the first equation.

*SGN2*SGN2 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the second equation.

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

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

*SCALE**WORK**INFO***Optimizations:**- Replaced level-3 by level-2 and level-1 calls,
- Replaced DAXPY operation by Fortran intrinsic
- Replaced all BLAS calls except of DGER by Fortran intrinsic,

**Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **225** of file **dla_tgcsylv_l2.f90**.

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

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

**Purpose:**

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

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

**Parameters***TRANSA**TRANSB**SGN1*SGN1 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the first equation.

*SGN2*SGN2 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the second equation.

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

*E*E is DOUBLE PRECISION array, dimension (LDE,N) On input, the matrix E contains the right hand side E. !> On output, the matrix E contains the solution R of Equation (1) On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices.

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

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

*SCALE**WORK**INFO***Optimizations:**- Replaced level-3 by level-2 and level-1 calls,
- Reorder the solution order to column first style,
- Replaced DAXPY operation by Fortran intrinsic
- Replaced all BLAS calls except of DTRMV and DGER by Fortran intrinsic,
- Use local copies of A, B, C, D, E, and F (M, N <=128) .

**Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **227** of file **dla_tgcsylv_l2_opt_local_copy.f90**.

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

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

**Purpose:**

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

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

**Parameters***TRANSA**TRANSB**SGN1*SGN1 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the first equation.

*SGN2*SGN2 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the second equation.

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

*E*E is DOUBLE PRECISION array, dimension (LDE,N) On input, the matrix E contains the right hand side E. !> On output, the matrix E contains the solution R of Equation (1) On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices.

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

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

*SCALE**WORK**INFO***Optimizations:**- Replaced level-3 by level-2 and level-1 calls,
- Reorder the solution order to column first style,
- Replaced DAXPY operation by Fortran intrinsic
- Replaced all BLAS calls except of DTRMV and DGER by Fortran intrinsic,
- Align local copies and fix problem size to <= 128

**Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **225** of file **dla_tgcsylv_l2_opt_local_copy_128.f90**.

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

**Purpose:**

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

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

**Parameters***TRANSA**TRANSB**SGN1*SGN1 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the first equation.

*SGN2*SGN2 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the second equation.

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

*E*E is DOUBLE PRECISION array, dimension (LDE,N) On input, the matrix E contains the right hand side E. !> On output, the matrix E contains the solution R of Equation (1) On output, the matrix E contains the solution R of Equation (1) as selected by TRANSA, TRANSB, and SGN1/SGN2. Right hand side E and the solution R are M-by-N matrices.

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

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

*SCALE**WORK**INFO***Optimizations:**- Replaced level-3 by level-2 and level-1 calls,
- Reorder the solution order to column first style,
- Replaced DAXPY operation by Fortran intrinsic
- Replaced all BLAS calls except of DTRMV and DGER by Fortran intrinsic,
- Align local copies and fix problem size to <= 32

**Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **226** of file **dla_tgcsylv_l2_opt_local_copy_32.f90**.

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

**Purpose:**

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

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

**Parameters***TRANSA**TRANSB**SGN1*SGN1 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the first equation.

*SGN2*SGN2 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the second equation.

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

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

*SCALE**WORK**INFO***Optimizations:**- Replaced level-3 by level-2 and level-1 calls,
- Reorder the solution order to column first style,
- Replaced DAXPY operation by Fortran intrinsic
- Replaced all BLAS calls except of DTRMV and DGER by Fortran intrinsic,
- Align local copies and fix problem size to <= 64

**Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **226** of file **dla_tgcsylv_l2_opt_local_copy_64.f90**.

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

**Purpose:**

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

**Remarks****Parameters***TRANSA**TRANSB**SGN1*SGN1 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the first equation.

*SGN2*SGN2 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the second equation.

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

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

*SCALE**WORK**INFO***Optimizations:**- Replaced level-3 by level-2 and level-1 calls,
- Reorder the solution order to column first style,
- Replaced DAXPY operation by Fortran intrinsic
- Replaced all BLAS calls except of DTRMV and DGER by Fortran intrinsic,
- Align local copies and fix problem size to <= 96

**Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **225** of file **dla_tgcsylv_l2_opt_local_copy_96.f90**.

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

**Purpose:**

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

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

**Parameters***TRANSA**TRANSB**SGN1*SGN1 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the first equation.

*SGN2*SGN2 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the second equation.

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

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

*SCALE**WORK**INFO***Optimizations:**- Replaced level-3 by level-2 and level-1 calls,
- Reorder the solution order to column first style,
- Replaced DAXPY operation by Fortran intrinsic
- Replaced all BLAS calls except of DTRMV and DGER by Fortran intrinsic,
- Column first computation of the solution.

**Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **226** of file **dla_tgcsylv_l2_opt_reorder.f90**.

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

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

**Purpose:**

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

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

**Parameters***TRANSA**TRANSB**SGN1*SGN1 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the first equation.

*SGN2*SGN2 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the second equation.

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

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

*SCALE**WORK**INFO***Optimizations:**- Nothing.

**Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **225** of file **dla_tgcsylv_l2_unopt.f90**.

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

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

**Purpose:**

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

**Attention**The algorithm is implemented using BLAS level 3 operations.

**Parameters***TRANSA**TRANSB**SGN1*SGN1 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the first equation.

*SGN2*SGN2 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the second equation.

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

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

*SCALE**WORK**INFO***Optimizations:**- Nothing but standard L3 BLAS calls.
- Changed computation order to column major

**Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **223** of file **dla_tgcsylv_l3_opt.f90**.

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

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

**Purpose:**

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

**Remarks**This algorithms uses 2 Stage blocking with DAG scheduled solves.

**Parameters***TRANSA**TRANSB**SGN1*SGN1 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the first equation.

*SGN2*SGN2 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the second equation.

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

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

*SCALE**WORK**INFO***Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **217** of file **dla_tgcsylv_l3_2stage.f90**.

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

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

**Purpose:**

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

**Attention**The algorithm is implemented using BLAS level 3 operations.

**Parameters***TRANSA**TRANSB**SGN1*SGN1 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the first equation.

*SGN2*SGN2 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the second equation.

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

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

*SCALE**WORK**INFO***Optimizations:**- Nothing but standard L3 BLAS calls.

**Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **223** of file **dla_tgcsylv_l3_unopt.f90**.

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

Recursive Algorithm for the generalized coupled Sylvester equation.

**Purpose:**

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

**Attention**The algorithm is implemented using BLAS level 3 operations.

**Parameters***TRANSA**TRANSB**SGN1*SGN1 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the first equation.

*SGN2*SGN2 is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between in the second equation.

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

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

*SCALE**WORK*WORK is DOUBLE PRECISION array, dimension 1 Workspace for the algorithm

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

**Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **205** of file **dla_tgcsylv_recursive.f90**.

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

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

**Purpose:**

DLA_TGSYLV_DAG 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 DGGES 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 DLA_TGSYLV_L3.

**Parameters***TRANSA**TRANSB**SGN*SGN is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between the two parts of the Sylvester equation. = 1 : Solve Equation (1) == -1: Solve Equation (2)

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

*X*X is DOUBLE PRECISION array, dimension (LDX,N) On input, the matrix X contains the right hand side Y. On output, the matrix X contains the solution of Equation (1) or (2) as selected by TRANSA, TRANSB, and SGN. Right hand side Y and the solution X are M-by-N matrices.

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

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

*WORK**INFO***Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **208** of file **dla_tgsylv_dag.f90**.

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

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

**Purpose:**

DLA_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 DGGES 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**TRANSB**SGN*SGN is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between the two parts of the Sylvester equation. = 1 : Solve Equation (1) == -1: Solve Equation (2)

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

*X*X is DOUBLE PRECISION array, dimension (LDX,N) On input, the matrix X contains the right hand side Y. On output, the matrix X contains the solution of Equation (1) or (2) as selected by TRANSA, TRANSB, and SGN. Right hand side Y and the solution X are M-by-N matrices.

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

*SCALE**WORK**INFO***Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **204** of file **dla_tgsylv_gardiner_laub.f90**.

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

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

**Purpose:**

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

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

**Parameters***TRANSA**TRANSB**SGN*SGN is DOUBLE PRECISION, allowed values: +/-1 Specifies the sign between the two parts of the Sylvester equation. = 1 : Solve Equation (1) == -1: Solve Equation (2)

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

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

*X*X is DOUBLE PRECISION array, dimension (LDX,N) On input, the matrix X contains the right hand side Y. On output, the matrix X contains the solution of Equation (1) or (2) as selected by TRANSA, TRANSB, and SGN. Right hand side Y and the solution X are M-by-N matrices.

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

*SCALE**WORK**INFO***Optimizations:**- Replaced level-3 by level-2 and level-1 calls,

**Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **207** of file **dla_tgsylv_l2.f90**.

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

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

**Purpose:**

DLA_TGSYLV_L2_UNOPT 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 DGGES from LAPACK.

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

**Parameters***TRANSA**TRANSB**SGN**M*M is INTEGER The order of the matrices A and C. M >= 0.

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

*SCALE**WORK**INFO***Optimizations:**- Nothing.

**Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **208** of file **dla_tgsylv_l2_colwise.f90**.

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

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

**Purpose:**

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

**Attention**The algorithm solves problems of dimension at most 128 ( M,N <=128)

**Parameters***TRANSA**TRANSB**SGN**M*M is INTEGER The order of the matrices A and C. 128 >= M >= 0.

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

*SCALE**WORK**INFO***Optimizations:**- Replaced level-3 by level-2 and level-1 calls,
- Reorder the solution order to column first style,
- Replaced DAXPY operation by Fortran intrinsic
- Replaced all BLAS calls except of DTRMV and DGER by Fortran intrinsic,
- Use local copies of A, B, C, D, and X (M, N <=128) .
- DTRMV works sequential.

**Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **214** of file **dla_tgsylv_l2_opt_local_copy.f90**.

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

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

**Purpose:**

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

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

**Parameters***TRANSA**TRANSB**SGN**M*M is INTEGER The order of the matrices A and C. 128 >= M >= 0.

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

*SCALE**WORK**INFO***Optimizations:**- Replaced level-3 by level-2 and level-1 calls,
- Reorder the solution order to column first style,
- Replaced DAXPY operation by Fortran intrinsic
- Replaced all BLAS calls except of DTRMV by Fortran intrinsic,
- Use local copies of A, B, C, D, and X.
- DTRMV works purely sequential
- Align local copies and fix problem size to <= 128

**Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **215** of file **dla_tgsylv_l2_opt_local_copy_128.f90**.

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

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

**Purpose:**

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

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

**Parameters***TRANSA**TRANSB**SGN**M*M is INTEGER The order of the matrices A and C. 32 >= M >= 0.

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

*SCALE**WORK**INFO***Optimizations:**- Replaced level-3 by level-2 and level-1 calls,
- Reorder the solution order to column first style,
- Replaced DAXPY operation by Fortran intrinsic
- Replaced all BLAS calls except of DTRMV by Fortran intrinsic,
- Use local copies of A, B, C, D, and X.
- DTRMV works purely sequential
- Align local copies and fix problem size to <= 32

**Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **215** of file **dla_tgsylv_l2_opt_local_copy_32.f90**.

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

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

**Purpose:**

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

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

**Parameters***TRANSA**TRANSB**SGN**M*M is INTEGER The order of the matrices A and C. 64 >= M >= 0.

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

*SCALE**WORK**INFO***Optimizations:**- Replaced level-3 by level-2 and level-1 calls,
- Reorder the solution order to column first style,
- Replaced DAXPY operation by Fortran intrinsic
- Replaced all BLAS calls except of DTRMV by Fortran intrinsic,
- Use local copies of A, B, C, D, and X.
- DTRMV works purely sequential
- Align local copies and fix problem size to <= 64

**Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **215** of file **dla_tgsylv_l2_opt_local_copy_64.f90**.

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

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

**Purpose:**

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

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

**Parameters***TRANSA**TRANSB**SGN**M*M is INTEGER The order of the matrices A and C. 96 >= M >= 0.

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

*SCALE**WORK**INFO***Optimizations:**- Replaced level-3 by level-2 and level-1 calls,
- Reorder the solution order to column first style,
- Replaced DAXPY operation by Fortran intrinsic
- Replaced all BLAS calls except of DTRMV by Fortran intrinsic,
- Use local copies of A, B, C, D, and X.
- DTRMV works purely sequential
- Align local copies and fix problem size to <= 96

**Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **215** of file **dla_tgsylv_l2_opt_local_copy_96.f90**.

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

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

**Purpose:**

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

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

**Parameters***TRANSA**TRANSB**SGN**M*M is INTEGER The order of the matrices A and C. M >= 0.

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

*SCALE**WORK**INFO***Optimizations:**- Replaced level-3 by level-2 and level-1 calls,
- Reorder the solution order to column first style,
- Replaced DAXPY operation by Fortran intrinsic
- Replaced all BLAS calls except of DTRMV and DGER by Fortran intrinsic,

**Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **211** of file **dla_tgsylv_l2_reorder.f90**.

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

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

**Purpose:**

DLA_TGSYLV_L2_UNOPT solves a generalized Sylvester equation of the following forms op1(A) * X * op2(B) + op1(C) * X * op2(D) = SCALE * Y (1) 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 DGGES from LAPACK.

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

**Parameters***TRANSA**TRANSB**SGN**M*M is INTEGER The order of the matrices A and C. M >= 0.

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

*SCALE**WORK**INFO***Optimizations:**- Nothing.

**Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **208** of file **dla_tgsylv_l2_rowwise.f90**.

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

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

**Purpose:**

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

**Attention**The algorithm is implemented using BLAS level 3 operations.

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

**Parameters***TRANSA**TRANSB**SGN**M*M is INTEGER The order of the matrices A and C. M >= 0.

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

*SCALE**WORK**INFO***Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **205** of file **dla_tgsylv_l3_opt.f90**.

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

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

**Purpose:**

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

**Attention**The algorithm is implemented using BLAS level 3 operations.

**Parameters***TRANSA**TRANSB**SGN**M*M is INTEGER The order of the matrices A and C. M >= 0.

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

*SCALE**WORK**INFO***Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **202** of file **dla_tgsylv_l3_2stage.f90**.

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

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

**Purpose:**

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

**Attention**The algorithm is implemented using BLAS level 3 operations.

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

**Parameters***TRANSA**TRANSB**SGN**M*M is INTEGER The order of the matrices A and C. M >= 0.

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

*SCALE**WORK**INFO***Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **204** of file **dla_tgsylv_l3_colwise.f90**.

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

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

**Purpose:**

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

**Attention**The algorithm is implemented using BLAS level 3 operations.

**Parameters***TRANSA**TRANSB**SGN**M*M is INTEGER The order of the matrices A and C. M >= 0.

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

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

*SCALE**WORK**INFO***Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **203** of file **dla_tgsylv_l3_rowwise.f90**.

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

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

**Purpose:**

DLA_TGSYLV_RECURSIVE solves a generalized Sylvester equation of the following forms op1(A) * X * op2(B) + op1(C) * X * op2(D) = SCALE * Y (1) 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 DGGES from LAPACK.

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

**Parameters***TRANSA**TRANSB**SGN**M*M is INTEGER The order of the matrices A and C. M >= 0.

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

*A*A is DOUBLE PRECISION array, dimension (LDA,M) The matrix A must be (quasi-) upper triangular.

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

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

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

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

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

*X*X is DOUBLE PRECISION array, dimension (LDX,N) On input, the matrix X contains the right hand side Y. On output, the matrix X contains the solution 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**WORK*WORK is DOUBLE PRECISION array, dimension M*N Workspace for the algorithm

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

**Author**Martin Koehler, MPI Magdeburg

**Date**January 2024

Definition at line **190** of file **dla_tgsylv_recursive.f90**.

## Author

Generated automatically by Doxygen for MEPACK from the source code.