X = Matrices.discreteSylvester(A, B, C);
         X = Matrices.discreteSylvester(A, B, C, AisHess, BTisSchur, sgn, eps);
Function discreteSylvester computes the solution X of the discrete-time Sylvester equation
A*X*B + sgn*X = C.
where sgn = 1 or sgn = -1. The algorithm applies the Hessenberg-Schur method proposed by Golub et al [1]. For sgn = -1, the discrete Sylvester equation is also known as Stein equation:
A*X*B - X + Q = 0.
In a nutshell, the problem is reduced to the corresponding problem
H*Y*S' + sgn*Y = F.
with
H=U'*A*U
is the Hessenberg form of A and
S=V'*B'*V
is the real Schur form of B',
F=U'*C*V
and
Y=U*X*V'
are appropriate transformations of C and
X. This problem is solved sequentially by
exploiting the specific forms of S and
H. Finally the solution of the original problem is
recovered as
X=U'*Y*V.
The Boolean inputs "AisHess" and "BTisSchur" indicate to omit one
or both of the transformation to Hessenberg form or Schur form
respectively in the case that A and/or
B have already Hessenberg form or Schur
respectively.
  [1] Golub, G.H., Nash, S. and Van Loan, C.F.
      A Hessenberg-Schur method for the problem AX + XB = C.
      IEEE Transaction on Automatic Control, AC-24, no. 6, pp. 909-913, 1979.
  A = [1.0,   2.0,   3.0;
       6.0,   7.0,   8.0;
       9.0,   2.0,   3.0];
  B = [7.0,   2.0,   3.0;
       2.0,   1.0,   2.0;
       3.0,   4.0,   1.0];
  C = [271.0,   135.0,   147.0;
       923.0,   494.0,   482.0;
       578.0,   383.0,   287.0];
  X = discreteSylvester(A, B, C);
  results in:
  X = [2.0,   3.0,   6.0;
       4.0,   7.0,   1.0;
       5.0,   3.0,   2.0];
Matrices.continuousSylvester, Matrices.discreteLyapunov
function discreteSylvester extends Modelica.Icons.Function; import Modelica.Math.Matrices; input Real A[:, size(A, 1)] "Square matrix A in A*X*B + sgn*X = C"; input Real B[:, size(B, 1)] "Square matrix B in A*X*B + sgn*X = C"; input Real C[size(A, 2), size(B, 1)] "Rectangular matrix C in A*X*B + sgn*X = C"; input Boolean AisHess = false "True if A has already Hessenberg form"; input Boolean BTisSchur = false "True if B' has already real Schur form"; input Integer sgn = 1 "Specifies the sign in A*X*B + sgn*X = C"; input Real eps = Matrices.norm(A, 1) * 10 * Modelica.Constants.eps "Tolerance"; output Real X[size(A, 2), size(B, 1)] "solution of the discrete Sylvester equation A*X*B + sgn*X = C"; end discreteSylvester;