dbltrsylv - Man Page

subroutine dla_trsylv2_dag (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info)
Level-3 Bartels-Stewart Algorithm for the discrete time Sylvester equation with DAG based parallelization.
subroutine dla_trsylv_dag (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info)
Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation with DAG parallelization.
subroutine dla_trsylv2_kernel_44nn (sgn, m, n, a, lda, b, ldb, x, ldx, scale, info)
Solver for a 4x4 discrete time Sylvester equation (TRANSA = N, TRANSB = N)
subroutine dla_trsylv2_kernel_44nt (sgn, m, n, a, lda, b, ldb, x, ldx, scale, info)
Solver for a 4x4 discrete time Sylvester equation (TRANSA = N, TRANSB = T)
subroutine dla_trsylv2_kernel_44tn (sgn, m, n, a, lda, b, ldb, x, ldx, scale, info)
Solver for a 4x4 discrete time Sylvester equation (TRANSA = T, TRANSB = N)
subroutine dla_trsylv2_kernel_44tt (sgn, m, n, a, lda, b, ldb, x, ldx, scale, info)
Solver for a 4x4 discrete time Sylvester equation (TRANSA = T, TRANSB = T)
subroutine dla_trsylv_kernel_44nn (sgn, m, n, a, lda, b, ldb, x, ldx, scale, info)
Solver for a 4x4 Sylvester equation (TRANSA = N, TRANSB = N)
subroutine dla_trsylv_kernel_44nt (sgn, m, n, a, lda, b, ldb, x, ldx, scale, info)
Solver for a 4x4 Sylvester equation (TRANSA = N, TRANSB = T)
subroutine dla_trsylv_kernel_44tn (sgn, m, n, a, lda, b, ldb, x, ldx, scale, info)
Solver for a 4x4 Sylvester equation (TRANSA = T, TRANSB = N)
subroutine dla_trsylv_kernel_44tt (sgn, m, n, a, lda, b, ldb, x, ldx, scale, info)
Solver for a 4x4 Sylvester equation (TRANSA = T, TRANSB = T)
subroutine dla_trsylv2_l2 (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info)
Level-2 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized)
subroutine dla_trsylv2_l2_local_copy (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info)
Level-2 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized)
subroutine dla_trsylv2_l2_local_copy_128 (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info)
Level-2 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized)
subroutine dla_trsylv2_l2_local_copy_32 (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info)
Level-2 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized)
subroutine dla_trsylv2_l2_local_copy_64 (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info)
Level-2 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized)
subroutine dla_trsylv2_l2_local_copy_96 (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info)
Level-2 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized)
subroutine dla_trsylv2_l2_reorder (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info)
Level-2 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized)
subroutine dla_trsylv2_l2_unopt (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info)
Level-2 Bartels-Stewart Algorithm for the discrete time Sylvester equation (unoptimized)
subroutine dla_trsylv_l2 (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info)
Level-2 Bartels-Stewart Algorithm for the Sylvester equation (optimized)
subroutine dla_trsylv_l2_local_copy (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info)
Level-2 Bartels-Stewart Algorithm for the Sylvester equation (optimized)
subroutine dla_trsylv_l2_local_copy_128 (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info)
Level-2 Bartels-Stewart Algorithm for the generalized Sylvester equation (Optimized)
subroutine dla_trsylv_l2_local_copy_32 (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info)
Level-2 Bartels-Stewart Algorithm for the generalized Sylvester equation (Optimized)
subroutine dla_trsylv_l2_local_copy_64 (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info)
Level-2 Bartels-Stewart Algorithm for the generalized Sylvester equation (Optimized)
subroutine dla_trsylv_l2_local_copy_96 (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info)
Level-2 Bartels-Stewart Algorithm for the generalized Sylvester equation (Optimized)
subroutine dla_trsylv_l2_reorder (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info)
Level-2 Bartels-Stewart Algorithm for the Sylvester equation (optimized)
subroutine dla_trsylv_l2_unopt (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info)
Level-2 Bartels-Stewart Algorithm for the Sylvester equation (unoptimized)
subroutine dla_trsylv2_l3_2s (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info)
Two level blocked Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized)
subroutine dla_trsylv2_l3 (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info)
Level-3 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized)
subroutine dla_trsylv2_l3_unopt (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info)
Level-3 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Not Optimized)
subroutine dla_trsylv_l3_2s (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info)
Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation.
subroutine dla_trsylv_l3 (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info)
Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation.
subroutine dla_trsylv_l3_unopt (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info)
Level-3 Bartels-Stewart Algorithm for the generalized Sylvester equation (unoptimized)
recursive subroutine dla_trsylv2_recursive (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info)
Level-3 Recursive Blocked Algorithm for the discrete time Sylvester equation.
recursive subroutine dla_trsylv_recursive (transa, transb, sgn, m, n, a, lda, b, ldb, x, ldx, scale, work, info)
Level-3 Recursive Blocked Solver for the Sylvester equation.

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

Purpose:

!>
!> DLA_TRSYLV2_DAG solves a discrete time Sylvester equation of the following forms
!>
!>    op1(A) * X * op2(B) +  X  = SCALE * Y                              (1)
!>
!> or
!>
!>    op1(A) * X * op2(B) -  X  = SCALE * Y                              (2)
!>
!> where A is a M-by-M quasi upper triangular matrix and B is a N-by-N upper
!> quasi triangular matrices. The right hand side Y and the solution X
!> M-by-N matrices. Typically the matrices A and B are created by DGEES 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 the DLA_TRSYLV2_L3 routine.

Parameters

TRANSA

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

TRANSB

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

SGN

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

M

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

N

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

A

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

LDA

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

B

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

LDB

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

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

!>          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 181 of file dla_trsylv2_dag.f90.

Solver for a 4x4 discrete time Sylvester equation (TRANSA = N, TRANSB = N)

Purpose:

!>
!> DLA_TRSYLV2_KERNEL_44NN solves a discrete time Sylvester equation of the following form
!>
!>    A * X * B + SGN * X  = SCALE * Y                              (1)
!>
!> where A is a M-by-M quasi upper triangular matrix and B is a N-by-N quasi upper
!> triangular matrix.  The right hand side Y and the solution X
!> M-by-N matrices. Typically the matrices A and B are created by DGEES from LAPACK.
!> The algorithm is implemented using BLAS level 2
!> operations. Thereby the order of M and N is at most 4. Furthermore, for fast execution
!> the function does not check the input arguments.
!>
!>

Parameters

SGN

!>          SGN is DOUBLE PRECISION, allowed values: +/-1
!>          Specifies the sign between the two parts of the Sylvester equation.
!>

M

!>          M is INTEGER
!>          The order of the matrix A.  4 >= M >= 0.
!>

N

!>          N is INTEGER
!>          The order of the matrix B.  4 >= 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.
!>

LDB

!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= 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)
!>          as selected by SGN.
!>          Right hand side Y and the solution X are M-by-N matrices.
!>

LDX

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

SCALE

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

INFO

!>          INFO is INTEGER
!>          On input:
!>          = 1 :  Skip the input data checks
!>          <> 1:  Check input data like normal LAPACK like routines.
!>          On output:
!>          == 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  The equation is not solved correctly. One of the arising inner
!>                system got singular.
!>

Author

Martin Koehler, MPI Magdeburg

Date

January 2024

Definition at line 135 of file dla_trsylv2_kernel_44_nn.f90.

Solver for a 4x4 discrete time Sylvester equation (TRANSA = N, TRANSB = T)

Purpose:

!>
!> DLA_TRSYLV2_KERNEL_44NT solves a discrete time Sylvester equation of the following form
!>
!>    A * X * B**T + SGN * X  = SCALE * Y                              (1)
!>
!> where A is a M-by-M quasi upper triangular matrix and B is a N-by-N quasi upper
!> triangular matrix.  The right hand side Y and the solution X
!> M-by-N matrices. Typically the matrices A and B are created by DGEES from LAPACK.
!> The algorithm is implemented using BLAS level 2
!> operations. Thereby the order of M and N is at most 4. Furthermore, for fast execution
!> the function does not check the input arguments.
!>
!>

Parameters

SGN

!>          SGN is DOUBLE PRECISION, allowed values: +/-1
!>          Specifies the sign between the two parts of the Sylvester equation.
!>

M

!>          M is INTEGER
!>          The order of the matrix A.  4 >= M >= 0.
!>

N

!>          N is INTEGER
!>          The order of the matrix B.  4 >= 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.
!>

LDB

!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= 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)
!>          as selected by SGN.
!>          Right hand side Y and the solution X are M-by-N matrices.
!>

LDX

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

SCALE

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

INFO

!>          INFO is INTEGER
!>          On input:
!>          = 1 :  Skip the input data checks
!>          <> 1:  Check input data like normal LAPACK like routines.
!>          On output:
!>          == 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  The equation is not solved correctly. One of the arising inner
!>                system got singular.
!>

Author

Martin Koehler, MPI Magdeburg

Date

January 2024

Definition at line 134 of file dla_trsylv2_kernel_44_nt.f90.

Solver for a 4x4 discrete time Sylvester equation (TRANSA = T, TRANSB = N)

Purpose:

!>
!> DLA_TRSYLV2_KERNEL_44TN solves a discrete time Sylvester equation of the following form
!>
!>    A**T * X * B + SGN * X  = SCALE * Y                              (1)
!>
!> where A is a M-by-M quasi upper triangular matrix and B is a N-by-N quasi upper
!> triangular matrix.  The right hand side Y and the solution X
!> M-by-N matrices. Typically the matrices A and B are created by DGEES from LAPACK.
!> The algorithm is implemented using BLAS level 2
!> operations. Thereby the order of M and N is at most 4. Furthermore, for fast execution
!> the function does not check the input arguments.
!>
!>

Parameters

SGN

!>          SGN is DOUBLE PRECISION, allowed values: +/-1
!>          Specifies the sign between the two parts of the Sylvester equation.
!>

M

!>          M is INTEGER
!>          The order of the matrix A.  4 >= M >= 0.
!>

N

!>          N is INTEGER
!>          The order of the matrix B.  4 >= 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.
!>

LDB

!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= 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)
!>          as selected by SGN.
!>          Right hand side Y and the solution X are M-by-N matrices.
!>

LDX

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

SCALE

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

INFO

!>          INFO is INTEGER
!>          On input:
!>          = 1 :  Skip the input data checks
!>          <> 1:  Check input data like normal LAPACK like routines.
!>          On output:
!>          == 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  The equation is not solved correctly. One of the arising inner
!>                system got singular.
!>

Author

Martin Koehler, MPI Magdeburg

Date

January 2024

Definition at line 135 of file dla_trsylv2_kernel_44_tn.f90.

Solver for a 4x4 discrete time Sylvester equation (TRANSA = T, TRANSB = T)

Purpose:

!>
!> DLA_TRSYLV2_KERNEL_44TT solves a discrete time Sylvester equation of the following form
!>
!>    A**T * X * B**T + SGN * X  = SCALE * Y                              (1)
!>
!> where A is a M-by-M quasi upper triangular matrix and B is a N-by-N quasi upper
!> triangular matrix.  The right hand side Y and the solution X
!> M-by-N matrices. Typically the matrices A and B are created by DGEES from LAPACK.
!> The algorithm is implemented using BLAS level 2
!> operations. Thereby the order of M and N is at most 4. Furthermore, for fast execution
!> the function does not check the input arguments.
!>
!>

Parameters

SGN

!>          SGN is DOUBLE PRECISION, allowed values: +/-1
!>          Specifies the sign between the two parts of the Sylvester equation.
!>

M

!>          M is INTEGER
!>          The order of the matrix A.  4 >= M >= 0.
!>

N

!>          N is INTEGER
!>          The order of the matrix B.  4 >= 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.
!>

LDB

!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= 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)
!>          as selected by SGN.
!>          Right hand side Y and the solution X are M-by-N matrices.
!>

LDX

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

SCALE

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

INFO

!>          INFO is INTEGER
!>          On input:
!>          = 1 :  Skip the input data checks
!>          <> 1:  Check input data like normal LAPACK like routines.
!>          On output:
!>          == 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  The equation is not solved correctly. One of the arising inner
!>                system got singular.
!>

Author

Martin Koehler, MPI Magdeburg

Date

January 2024

Definition at line 135 of file dla_trsylv2_kernel_44_tt.f90.

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

Purpose:

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

Parameters

TRANSA

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

TRANSB

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

SGN

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

M

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

N

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

A

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

LDA

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

B

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

LDB

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

X

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

LDX

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

SCALE

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

WORK

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

INFO

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

Optimizations:

• Replaced level-3 by level-2 and level-1 calls,
• 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 182 of file dla_trsylv2_l2.f90.

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

Purpose:

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

Parameters

TRANSA

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

TRANSB

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

SGN

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

M

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

N

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

A

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

LDA

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

B

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

LDB

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

X

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

LDX

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

SCALE

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

WORK

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

INFO

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

Optimizations:

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

Author

Martin Koehler, MPI Magdeburg

Date

January 2024

Definition at line 184 of file dla_trsylv2_l2_opt_local_copy.f90.

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

Purpose:

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

Attention

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

Parameters

TRANSA

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

TRANSB

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

SGN

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

M

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

N

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

A

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

LDA

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

B

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

LDB

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

X

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

LDX

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

SCALE

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

WORK

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

INFO

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

Optimizations:

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

Author

Martin Koehler, MPI Magdeburg

Date

January 2024

Definition at line 184 of file dla_trsylv2_l2_opt_local_copy_128.f90.

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

Purpose:

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

Attention

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

Parameters

TRANSA

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

TRANSB

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

SGN

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

M

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

N

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

A

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

LDA

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

B

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

LDB

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

X

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

LDX

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

SCALE

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

WORK

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

INFO

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

Optimizations:

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

Author

Martin Koehler, MPI Magdeburg

Date

January 2024

Definition at line 184 of file dla_trsylv2_l2_opt_local_copy_32.f90.

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

Purpose:

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

Attention

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

Parameters

TRANSA

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

TRANSB

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

SGN

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

M

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

N

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

A

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

LDA

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

B

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

LDB

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

X

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

LDX

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

SCALE

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

WORK

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

INFO

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

Optimizations:

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

Author

Martin Koehler, MPI Magdeburg

Date

January 2024

Definition at line 184 of file dla_trsylv2_l2_opt_local_copy_64.f90.

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

Purpose:

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

Attention

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

Parameters

TRANSA

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

TRANSB

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

SGN

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

M

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

N

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

A

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

LDA

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

B

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

LDB

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

X

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

LDX

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

SCALE

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

WORK

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

INFO

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

Optimizations:

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

Author

Martin Koehler, MPI Magdeburg

Date

January 2024

Definition at line 184 of file dla_trsylv2_l2_opt_local_copy_96.f90.

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

Purpose:

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

Parameters

TRANSA

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

TRANSB

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

SGN

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

M

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

N

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

A

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

LDA

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

B

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

LDB

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

X

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

LDX

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

SCALE

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

WORK

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

INFO

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

Optimizations:

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

Author

Martin Koehler, MPI Magdeburg

Date

January 2024

Definition at line 183 of file dla_trsylv2_l2_opt_reorder.f90.

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

Purpose:

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

Attention

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

Parameters

TRANSA

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

TRANSB

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

SGN

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

M

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

N

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

A

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

LDA

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

B

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

LDB

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

X

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

LDX

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

SCALE

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

WORK

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

INFO

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

Optimizations:

• Nothing.

Author

Martin Koehler, MPI Magdeburg

Date

January 2024

Definition at line 181 of file dla_trsylv2_l2_unopt.f90.

Level-3 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized)

Purpose:

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

Parameters

TRANSA

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

TRANSB

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

SGN

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

M

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

N

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

A

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

LDA

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

B

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

LDB

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

X

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

LDX

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

SCALE

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

WORK

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

INFO

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

Author

Martin Koehler, MPI Magdeburg

Date

January 2024

Definition at line 173 of file dla_trsylv2_l3_opt.f90.

Two level blocked Bartels-Stewart Algorithm for the discrete time Sylvester equation (Optimized)

Purpose:

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

Remarks

The algorithm used level-3 BLAS operations and a DAG scheduled inner solver.

Attention

Due to the parallel nature of the inner solvers the scaling is turned off and SCALE is set to ONE.

Parameters

TRANSA

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

TRANSB

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

SGN

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

M

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

N

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

A

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

LDA

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

B

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

LDB

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

X

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

LDX

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

SCALE

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

WORK

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

INFO

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

Author

Martin Koehler, MPI Magdeburg

Date

January 2024

Definition at line 174 of file dla_trsylv2_l3_2stage.f90.

Level-3 Bartels-Stewart Algorithm for the discrete time Sylvester equation (Not Optimized)

Purpose:

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

Attention

This function iterates column first over the result and from this fact it will be a bit slower than DLA_TRSYLV2_L3.

Parameters

TRANSA

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

TRANSB

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

SGN

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

M

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

N

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

A

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

LDA

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

B

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

LDB

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

X

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

LDX

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

SCALE

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

WORK

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

INFO

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

Author

Martin Koehler, MPI Magdeburg

Date

January 2024

Definition at line 175 of file dla_trsylv2_l3_unopt.f90.

Level-3 Recursive Blocked Algorithm for the discrete time Sylvester equation.

Purpose:

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

Remarks

The algorithm uses recursive blocking instead of the Bartels-Stewart approach.

Parameters

TRANSA

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

TRANSB

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

SGN

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

M

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

N

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

A

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

LDA

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

B

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

LDB

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

X

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

LDX

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

SCALE

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

WORK

!>          WORK is DOUBLE PRECISION array, dimension 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 163 of file dla_trsylv2_recursive.f90.

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

Purpose:

!> DLA_TRSYLV_DAG solves a Sylvester equation of the following forms
!>
!>    op1(A) * X  +  X * op2(B) = SCALE * Y                              (1)
!>
!> or
!>
!>    op1(A) * X  -  X * op2(B) = SCALE * Y                              (2)
!>
!> where A is a M-by-M quasi upper triangular matrix, B is a N-by-N quasi upper triangular
!> matrix. The right hand side Y and the solution X are M-by-N matrices.  Typically the matrices
!> A and B are generated via DGEES form 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_TRSYLV_L3.

Parameters

TRANSA

!>          TRANSA is CHARACTER(1)
!>          Specifies the form of the system of equations with respect to A:
!>          == 'N':  op1(A) = A
!>          == 'T':  op1(A) = A**T
!>

TRANSB

!>          TRANSB is CHARACTER(1)
!>          Specifies the form of the system of equations with respect to B:
!>          == 'N':  op2(B) = B,
!>          == 'T':  op2(B) = B**T
!>

SGN

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

M

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

N

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

A

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

LDA

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

B

!>          B is DOUBLE PRECISION array, dimension (LDB,N)
!>          The matrix B must be (quasi-) upper triangular.
!>

LDB

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

X

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

LDX

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

SCALE

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

WORK

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

INFO

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

Author

Martin Koehler, MPI Magdeburg

Date

January 2024

Definition at line 177 of file dla_trsylv_dag.f90.

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

Purpose:

!>
!> DLA_TRSYLV_KERNEL_44NN solves a Sylvester equation of the following form
!>
!>    A * X  + SGN * X * B = SCALE * Y                              (1)
!>
!> where A is a M-by-M quasi upper triangular matrix and B is N-by-N quasi
!> upper triangular matrices. The right hand side Y and the solution X
!> M-by-N matrices. Typically the matrices A and B are create by DGEES from LAPACK.
!> The algorithm is implemented without using BLAS level 2
!> operations. Thereby the order of M and N is at most 4. Furthermore, for fast execution
!> the function does not check the input arguments.
!>
!>

Parameters

SGN

!>          SGN is DOUBLE PRECISION, allowed values: +/-1
!>          Specifies the sign between the two parts of the Sylvester equation.
!>

M

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

N

!>          N is INTEGER
!>          The order of the matrices B and D.  4 >= 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).
!>

X

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

LDX

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

SCALE

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

INFO

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

Author

Martin Koehler, MPI Magdeburg

Date

January 2024

Definition at line 134 of file dla_trsylv_kernel_44_nn.f90.

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

Purpose:

!>
!> DLA_TRSYLV_KERNEL_44NT solves a Sylvester equation of the following form
!>
!>    A * X  + SGN * X * B**T = SCALE * Y                              (1)
!>
!> where A is a M-by-M quasi upper triangular matrix and B is N-by-N quasi
!> upper triangular matrices. The right hand side Y and the solution X
!> M-by-N matrices. Typically the matrices A and B are create by DGEES from LAPACK.
!> The algorithm is implemented using without BLAS level 2
!> operations. Thereby the order of M and N is at most 4. Furthermore, for fast execution
!> the function does not check the input arguments.
!>
!>

Parameters

SGN

!>          SGN is DOUBLE PRECISION, allowed values: +/-1
!>          Specifies the sign between the two parts of the Sylvester equation.
!>

M

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

N

!>          N is INTEGER
!>          The order of the matrices B and D.  4 >= 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).
!>

X

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

LDX

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

SCALE

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

INFO

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

Author

Martin Koehler, MPI Magdeburg

Date

January 2024

Definition at line 133 of file dla_trsylv_kernel_44_nt.f90.

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

Purpose:

!>
!> DLA_TRSYLV_KERNEL_44TN solves a Sylvester equation of the following form
!>
!>    A**T * X  + SGN * X * B = SCALE * Y                              (1)
!>
!> where A is a M-by-M quasi upper triangular matrix and B is N-by-N quasi
!> upper triangular matrices. The right hand side Y and the solution X
!> M-by-N matrices. Typically the matrices A and B are create by DGEES from LAPACK.
!> The algorithm is implemented without using BLAS level 2
!> operations. Thereby the order of M and N is at most 4. Furthermore, for fast execution
!> the function does not check the input arguments.
!>
!>

Parameters

SGN

!>          SGN is DOUBLE PRECISION, allowed values: +/-1
!>          Specifies the sign between the two parts of the Sylvester equation.
!>

M

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

N

!>          N is INTEGER
!>          The order of the matrices B and D.  4 >= 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).
!>

X

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

LDX

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

SCALE

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

INFO

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

Author

Martin Koehler, MPI Magdeburg

Date

January 2024

Definition at line 133 of file dla_trsylv_kernel_44_tn.f90.

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

Purpose:

!>
!> DLA_TRSYLV_KERNEL_44TT solves a Sylvester equation of the following form
!>
!>    A**T * X  + SGN * X * B**T = SCALE * Y                              (1)
!>
!> where A is a M-by-M quasi upper triangular matrix and B is N-by-N quasi
!> upper triangular matrices. The right hand side Y and the solution X
!> M-by-N matrices. Typically the matrices A and B are create by DGEES from LAPACK.
!> The algorithm is implemented without using BLAS level 2
!> operations. Thereby the order of M and N is at most 4. Furthermore, for fast execution
!> the function does not check the input arguments.
!>
!>

Parameters

SGN

!>          SGN is DOUBLE PRECISION, allowed values: +/-1
!>          Specifies the sign between the two parts of the Sylvester equation.
!>

M

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

N

!>          N is INTEGER
!>          The order of the matrices B and D.  4 >= 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).
!>

X

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

LDX

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

SCALE

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

INFO

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

Author

Martin Koehler, MPI Magdeburg

Date

January 2024

Definition at line 133 of file dla_trsylv_kernel_44_tt.f90.

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

Purpose:

!>
!> DLA_TRSYLV_L2 solves a Sylvester equation of the following forms
!>
!>    op1(A) * X  +  X * op2(B) = SCALE * Y                              (1)
!>
!> or
!>
!>    op1(A) * X  -  X * op2(B) = SCALE * Y                              (2)
!>
!> where A is a M-by-M quasi upper triangular matrix, B is a N-by-N quasi upper triangular
!> matrix. The right hand side Y and the solution X are M-by-N matrices.  Typically the matrices
!> A and B are generated via DGEES form LAPACK.
!>

Parameters

TRANSA

!>          TRANSA is CHARACTER(1)
!>          Specifies the form of the system of equations with respect to A:
!>          == 'N':  op1(A) = A
!>          == 'T':  op1(A) = A**T
!>

TRANSB

!>          TRANSB is CHARACTER(1)
!>          Specifies the form of the system of equations with respect to B:
!>          == 'N':  op2(B) = B,
!>          == 'T':  op2(B) = B**T
!>

SGN

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

M

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

N

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

A

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

LDA

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

B

!>          B is DOUBLE PRECISION array, dimension (LDB,N)
!>          The matrix B must be (quasi-) upper triangular.
!>

LDB

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

X

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

LDX

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

SCALE

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

WORK

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

INFO

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

Optimizations:

• Replaced level-3 by level-2 and level-1 calls,
• 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 180 of file dla_trsylv_l2.f90.

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

Purpose:

!>
!> DLA_TRSYLV_L2_LOCAL_COPY solves a Sylvester equation of the following forms
!>
!>    op1(A) * X  +  X * op2(B) = SCALE * Y                              (1)
!>
!> or
!>
!>    op1(A) * X  -  X * op2(B) = SCALE * Y                              (2)
!>
!> where A is a M-by-M quasi upper triangular matrix, B is a N-by-N quasi upper triangular
!> matrix. The right hand side Y and the solution X are M-by-N matrices.  Typically the matrices
!> A and B are generated via DGEES form LAPACK.
!>

Parameters

TRANSA

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

TRANSB

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

SGN

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

M

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

N

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

A

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

LDA

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

B

!>          B is DOUBLE PRECISION array, dimension (LDB,N)
!>          The matrix B must be (quasi-) upper triangular.
!>

LDB

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

X

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

LDX

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

SCALE

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

WORK

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

INFO

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

Optimizations:

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

Author

Martin Koehler, MPI Magdeburg

Date

January 2024

Definition at line 181 of file dla_trsylv_l2_opt_local_copy.f90.

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

Purpose:

!>
!> DLA_TRSYLV_L2_LOCAL_COPY_128 solves a Sylvester equation of the following forms
!>
!>    op1(A) * X  +  X * op2(B) = SCALE * Y                              (1)
!>
!> or
!>
!>    op1(A) * X  -  X * op2(B) = SCALE * Y                              (2)
!>
!> where A is a M-by-M quasi upper triangular matrix, B is a N-by-N quasi upper triangular
!> matrix. The right hand side Y and the solution X are M-by-N matrices.  Typically the matrices
!> A and B are generated via DGEES form LAPACK.
!>

Parameters

TRANSA

!>          TRANSA is CHARACTER(1)
!>          Specifies the form of the system of equations with respect to A:
!>          == 'N':  op1(A) = A
!>          == 'T':  op1(A) = A**T
!>

TRANSB

!>          TRANSB is CHARACTER(1)
!>          Specifies the form of the system of equations with respect to B:
!>          == 'N':  op2(B) = B,
!>          == 'T':  op2(B) = B**T
!>

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.  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.
!>

LDB

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

X

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

LDX

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

SCALE

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

WORK

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

INFO

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

Optimizations:

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

Author

Martin Koehler, MPI Magdeburg

Date

January 2024

Definition at line 183 of file dla_trsylv_l2_opt_local_copy_128.f90.

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

Purpose:

!>
!> DLA_TRSYLV_L2_LOCAL_COPY_32 solves a Sylvester equation of the following forms
!>
!>    op1(A) * X  +  X * op2(B) = SCALE * Y                              (1)
!>
!> or
!>
!>    op1(A) * X  -  X * op2(B) = SCALE * Y                              (2)
!>
!> where A is a M-by-M quasi upper triangular matrix, B is a N-by-N quasi upper triangular
!> matrix. The right hand side Y and the solution X are M-by-N matrices.  Typically the matrices
!> A and B are generated via DGEES form LAPACK.
!>

Parameters

TRANSA

!>          TRANSA is CHARACTER(1)
!>          Specifies the form of the system of equations with respect to A:
!>          == 'N':  op1(A) = A
!>          == 'T':  op1(A) = A**T
!>

TRANSB

!>          TRANSB is CHARACTER(1)
!>          Specifies the form of the system of equations with respect to B:
!>          == 'N':  op2(B) = B,
!>          == 'T':  op2(B) = B**T
!>

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.  32 >= M >= 0.
!>

N

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

A

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

LDA

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

B

!>          B is DOUBLE PRECISION array, dimension (LDB,N)
!>          The matrix B must be (quasi-) upper triangular.
!>

LDB

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

X

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

LDX

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

SCALE

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

WORK

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

INFO

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

Optimizations:

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

Author

Martin Koehler, MPI Magdeburg

Date

January 2024

Definition at line 183 of file dla_trsylv_l2_opt_local_copy_32.f90.

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

Purpose:

!>
!> DLA_TRSYLV_L2_LOCAL_COPY_64 solves a Sylvester equation of the following forms
!>
!>    op1(A) * X  +  X * op2(B) = SCALE * Y                              (1)
!>
!> or
!>
!>    op1(A) * X  -  X * op2(B) = SCALE * Y                              (2)
!>
!> where A is a M-by-M quasi upper triangular matrix, B is a N-by-N quasi upper triangular
!> matrix. The right hand side Y and the solution X are M-by-N matrices.  Typically the matrices
!> A and B are generated via DGEES form LAPACK.
!>

Parameters

TRANSA

!>          TRANSA is CHARACTER(1)
!>          Specifies the form of the system of equations with respect to A:
!>          == 'N':  op1(A) = A
!>          == 'T':  op1(A) = A**T
!>

TRANSB

!>          TRANSB is CHARACTER(1)
!>          Specifies the form of the system of equations with respect to B:
!>          == 'N':  op2(B) = B,
!>          == 'T':  op2(B) = B**T
!>

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.  64 >= M >= 0.
!>

N

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

A

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

LDA

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

B

!>          B is DOUBLE PRECISION array, dimension (LDB,N)
!>          The matrix B must be (quasi-) upper triangular.
!>

LDB

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

X

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

LDX

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

SCALE

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

WORK

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

INFO

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

Optimizations:

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

Author

Martin Koehler, MPI Magdeburg

Date

January 2024

Definition at line 182 of file dla_trsylv_l2_opt_local_copy_64.f90.

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

Purpose:

!>
!> DLA_TRSYLV_L2_LOCAL_COPY_96 solves a Sylvester equation of the following forms
!>
!>    op1(A) * X  +  X * op2(B) = SCALE * Y                              (1)
!>
!> or
!>
!>    op1(A) * X  -  X * op2(B) = SCALE * Y                              (2)
!>
!> where A is a M-by-M quasi upper triangular matrix, B is a N-by-N quasi upper triangular
!> matrix. The right hand side Y and the solution X are M-by-N matrices.  Typically the matrices
!> A and B are generated via DGEES form LAPACK.
!>

Parameters

TRANSA

!>          TRANSA is CHARACTER(1)
!>          Specifies the form of the system of equations with respect to A:
!>          == 'N':  op1(A) = A
!>          == 'T':  op1(A) = A**T
!>

TRANSB

!>          TRANSB is CHARACTER(1)
!>          Specifies the form of the system of equations with respect to B:
!>          == 'N':  op2(B) = B,
!>          == 'T':  op2(B) = B**T
!>

SGN

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

M

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

N

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

A

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

LDA

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

B

!>          B is DOUBLE PRECISION array, dimension (LDB,N)
!>          The matrix B must be (quasi-) upper triangular.
!>

LDB

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

X

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

LDX

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

SCALE

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

WORK

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

INFO

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

Optimizations:

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

Author

Martin Koehler, MPI Magdeburg

Date

January 2024

Definition at line 181 of file dla_trsylv_l2_opt_local_copy_96.f90.

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

Purpose:

!>
!> DLA_TRSYLV_L2_REORDER solves a Sylvester equation of the following forms
!>
!>    op1(A) * X  +  X * op2(B) = SCALE * Y                              (1)
!>
!> or
!>
!>    op1(A) * X  -  X * op2(B) = SCALE * Y                              (2)
!>
!> where A is a M-by-M quasi upper triangular matrix, B is a N-by-N quasi upper triangular
!> matrix. The right hand side Y and the solution X are M-by-N matrices.  Typically the matrices
!> A and B are generated via DGEES form LAPACK.
!>

Attention

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

Parameters

TRANSA

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

TRANSB

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

SGN

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

M

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

N

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

A

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

LDA

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

B

!>          B is DOUBLE PRECISION array, dimension (LDB,N)
!>          The matrix B must be (quasi-) upper triangular.
!>

LDB

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

X

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

LDX

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

SCALE

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

WORK

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

INFO

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

Optimizations:

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

Author

Martin Koehler, MPI Magdeburg

Date

January 2024

Definition at line 182 of file dla_trsylv_l2_opt_reorder.f90.

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

Purpose:

!>
!> DLA_TRSYLV_L2_UNOPT solves a Sylvester equation of the following forms
!>
!>    op1(A) * X  +  X * op2(B) = SCALE * Y                              (1)
!>
!> or
!>
!>    op1(A) * X  -  X * op2(B) = SCALE * Y                              (2)
!>
!> where A is a M-by-M quasi upper triangular matrix, B is a N-by-N quasi upper triangular
!> matrix. The right hand side Y and the solution X are M-by-N matrices.  Typically the matrices
!> A and B are generated via DGEES form LAPACK.
!>

Attention

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

Parameters

TRANSA

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

TRANSB

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

SGN

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

M

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

N

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

A

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

LDA

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

B

!>          B is DOUBLE PRECISION array, dimension (LDB,N)
!>          The matrix B must be (quasi-) upper triangular.
!>

LDB

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

X

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

LDX

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

SCALE

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

WORK

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

INFO

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

Optimizations:

• Nothing.

Author

Martin Koehler, MPI Magdeburg

Date

January 2024

Definition at line 180 of file dla_trsylv_l2_unopt.f90.

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

Purpose:

!> DLA_TRSYLV_L3 solves a Sylvester equation of the following forms
!>
!>    op1(A) * X  +  X * op2(B) = SCALE * Y                              (1)
!>
!> or
!>
!>    op1(A) * X  -  X * op2(B) = SCALE * Y                              (2)
!>
!> where A is a M-by-M quasi upper triangular matrix, B is a N-by-N quasi upper triangular
!> matrix. The right hand side Y and the solution X are M-by-N matrices.  Typically the matrices
!> A and B are generated via DGEES form LAPACK.
!>

Attention

The algorithm is implemented using BLAS level 3 operations.

Parameters

TRANSA

!>          TRANSA is CHARACTER(1)
!>          Specifies the form of the system of equations with respect to A:
!>          == 'N':  op1(A) = A
!>          == 'T':  op1(A) = A**T
!>

TRANSB

!>          TRANSB is CHARACTER(1)
!>          Specifies the form of the system of equations with respect to B:
!>          == 'N':  op2(B) = B,
!>          == 'T':  op2(B) = B**T
!>

SGN

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

M

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

N

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

A

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

LDA

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

B

!>          B is DOUBLE PRECISION array, dimension (LDB,N)
!>          The matrix B must be (quasi-) upper triangular.
!>

LDB

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

X

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

LDX

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

SCALE

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

WORK

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

INFO

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

Author

Martin Koehler, MPI Magdeburg

Date

January 2024

Definition at line 172 of file dla_trsylv_l3_opt.f90.

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

Purpose:

!> DLA_TRSYLV_L3_2S solves a Sylvester equation of the following forms
!>
!>    op1(A) * X  +  X * op2(B) = SCALE * Y                              (1)
!>
!> or
!>
!>    op1(A) * X  -  X * op2(B) = SCALE * Y                              (2)
!>
!> where A is a M-by-M quasi upper triangular matrix, B is a N-by-N quasi upper triangular
!> matrix. The right hand side Y and the solution X are M-by-N matrices.  Typically the matrices
!> A and B are generated via DGEES form LAPACK.
!>

Remarks

The algorithm is implemented using BLAS level 3 operations and a two level inner solver consisting of a massive parallel DAG scheduled solver from DLA_TRSYLV_DAG and optimized kernel solvers.

Attention

Due to the parallel nature of the inner solvers the scaling is turned off and SCALE is set to ONE.

Parameters

TRANSA

!>          TRANSA is CHARACTER(1)
!>          Specifies the form of the system of equations with respect to A:
!>          == 'N':  op1(A) = A
!>          == 'T':  op1(A) = A**T
!>

TRANSB

!>          TRANSB is CHARACTER(1)
!>          Specifies the form of the system of equations with respect to B:
!>          == 'N':  op2(B) = B,
!>          == 'T':  op2(B) = B**T
!>

SGN

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

M

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

N

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

A

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

LDA

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

B

!>          B is DOUBLE PRECISION array, dimension (LDB,N)
!>          The matrix B must be (quasi-) upper triangular.
!>

LDB

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

X

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

LDX

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

SCALE

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

WORK

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

INFO

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

Author

Martin Koehler, MPI Magdeburg

Date

January 2024

Definition at line 174 of file dla_trsylv_l3_2stage.f90.

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

Purpose:

!> DLA_TRSYLV_L3_UNOPT solves a Sylvester equation of the following forms
!>
!>    op1(A) * X  +  X * op2(B) = SCALE * Y                              (1)
!>
!> or
!>
!>    op1(A) * X  -  X * op2(B) = SCALE * Y                              (2)
!>
!> where A is a M-by-M quasi upper triangular matrix, B is a N-by-N quasi upper triangular
!> matrix. The right hand side Y and the solution X are M-by-N matrices.  Typically the matrices
!> A and B are generated via DGEES form LAPACK.
!>

Attention

The algorithm is implemented using BLAS level 3 operations.

Parameters

TRANSA

!>          TRANSA is CHARACTER(1)
!>          Specifies the form of the system of equations with respect to A:
!>          == 'N':  op1(A) = A
!>          == 'T':  op1(A) = A**T
!>

TRANSB

!>          TRANSB is CHARACTER(1)
!>          Specifies the form of the system of equations with respect to B:
!>          == 'N':  op2(B) = B,
!>          == 'T':  op2(B) = B**T
!>

SGN

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

M

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

N

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

A

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

LDA

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

B

!>          B is DOUBLE PRECISION array, dimension (LDB,N)
!>          The matrix B must be (quasi-) upper triangular.
!>

LDB

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

X

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

LDX

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

SCALE

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

WORK

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

INFO

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

Author

Martin Koehler, MPI Magdeburg

Date

January 2024

Definition at line 172 of file dla_trsylv_l3_unopt.f90.

Level-3 Recursive Blocked Solver for the Sylvester equation.

Purpose:

!> DLA_TRSYLV_RECURSIVE solves a Sylvester equation of the following forms
!>
!>    op1(A) * X  +  X * op2(B) = SCALE * Y                              (1)
!>
!> or
!>
!>    op1(A) * X  -  X * op2(B) = SCALE * Y                              (2)
!>
!> where A is a M-by-M quasi upper triangular matrix, B is a N-by-N quasi upper triangular
!> matrix. The right hand side Y and the solution X are M-by-N matrices.  Typically the matrices
!> A and B are generated via DGEES form LAPACK.
!>

Attention

The algorithm is implemented using BLAS level 3 operations and recursive blocking.

Parameters

TRANSA

!>          TRANSA is CHARACTER(1)
!>          Specifies the form of the system of equations with respect to A:
!>          == 'N':  op1(A) = A
!>          == 'T':  op1(A) = A**T
!>

TRANSB

!>          TRANSB is CHARACTER(1)
!>          Specifies the form of the system of equations with respect to B:
!>          == 'N':  op2(B) = B,
!>          == 'T':  op2(B) = B**T
!>

SGN

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

M

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

N

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

A

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

LDA

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

B

!>          B is DOUBLE PRECISION array, dimension (LDB,N)
!>          The matrix B must be (quasi-) upper triangular.
!>

LDB

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

X

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

LDX

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

SCALE

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

WORK

!>          WORK is DOUBLE PRECISION array, dimension 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 160 of file dla_trsylv_recursive.f90.