# Source code for pymor.algorithms.sylvester

```# This file is part of the pyMOR project (http://www.pymor.org).

import scipy.linalg as spla

from pymor.algorithms.to_matrix import to_matrix
from pymor.operators.interfaces import OperatorInterface
from pymor.operators.constructions import IdentityOperator

[docs]def solve_sylv_schur(A, Ar, E=None, Er=None, B=None, Br=None, C=None, Cr=None):
r"""Solve Sylvester equation by Schur decomposition.

Solves Sylvester equation

.. math::
A V E_r^T + E V A_r^T + B B_r^T = 0

or

.. math::
A^T W E_r + E^T W A_r + C^T C_r = 0

or both using (generalized) Schur decomposition (Algorithms 3 and 4
in [BKS11]_), if the necessary parameters are given.

Parameters
----------
A
Real |Operator|.
Ar
Real |Operator|.
It is converted into a |NumPy array| using
:func:`~pymor.algorithms.to_matrix.to_matrix`.
E
Real |Operator| or `None` (then assumed to be the identity).
Er
Real |Operator| or `None` (then assumed to be the identity).
It is converted into a |NumPy array| using
:func:`~pymor.algorithms.to_matrix.to_matrix`.
B
Real |Operator| or `None`.
Br
Real |Operator| or `None`.
It is converted into a |VectorArray| using
`Br.as_source_array()`.
C
Real |Operator| or `None`.
Cr
Real |Operator| or `None`.
It is converted into a |VectorArray| using
`Cr.as_range_array()`.

Returns
-------
V
Returned if `B` and `Br` are given, |VectorArray| from
`A.source`.
W
Returned if `C` and `Cr` are given, |VectorArray| from
`A.source`.

Raises
------
ValueError
If `V` and `W` cannot be returned.
"""
# check types
assert isinstance(A, OperatorInterface) and A.linear and A.source == A.range
assert isinstance(Ar, OperatorInterface) and Ar.linear and Ar.source == Ar.range

assert E is None or isinstance(E, OperatorInterface) and E.linear and E.source == E.range == A.source
if E is None:
E = IdentityOperator(A.source)
assert Er is None or isinstance(Er, OperatorInterface) and Er.linear and Er.source == Er.range == Ar.source

compute_V = B is not None and Br is not None
compute_W = C is not None and Cr is not None

if not compute_V and not compute_W:
raise ValueError('Not enough parameters are given to solve a Sylvester equation.')

if compute_V:
assert isinstance(B, OperatorInterface) and B.linear and B.range == A.source
assert isinstance(Br, OperatorInterface) and Br.linear and Br.range == Ar.source
assert B.source == Br.source

if compute_W:
assert isinstance(C, OperatorInterface) and C.linear and C.source == A.source
assert isinstance(Cr, OperatorInterface) and Cr.linear and Cr.source == Ar.source
assert C.range == Cr.range

# convert reduced operators
Ar = to_matrix(Ar, format='dense')
r = Ar.shape
if Er is not None:
Er = to_matrix(Er, format='dense')
if Br is not None:
Br = Br.as_source_array()
if Cr is not None:
Cr = Cr.as_range_array()

# (Generalized) Schur decomposition
if Er is None:
TAr, Z = spla.schur(Ar, output='complex')
Q = Z
else:
TAr, TEr, Q, Z = spla.qz(Ar, Er, output='complex')

# solve for V, from the last column to the first
if compute_V:
V = A.source.empty(reserve=r)

Br2 = Br.lincomb(Q.T)
BBr2 = B.apply(Br2)
for i in range(-1, -r - 1, -1):
rhs = -BBr2[i].copy()
if i < -1:
if Er is not None:
rhs -= A.apply(V.lincomb(TEr[i, :i:-1].conjugate()))
rhs -= E.apply(V.lincomb(TAr[i, :i:-1].conjugate()))
TErii = 1 if Er is None else TEr[i, i]
eAaE = TErii.conjugate() * A + TAr[i, i].conjugate() * E
V.append(eAaE.apply_inverse(rhs))

V = V.lincomb(Z.conjugate()[:, ::-1])
V = V.real

# solve for W, from the first column to the last
if compute_W:
W = A.source.empty(reserve=r)

Cr2 = Cr.lincomb(Z.T)
for i in range(r):
rhs = -CTCr2[i].copy()
if i > 0:
if Er is not None:
TErii = 1 if Er is None else TEr[i, i]
eAaE = TErii.conjugate() * A + TAr[i, i].conjugate() * E