--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mftoolbox/sylvsol.m	Wed May 06 14:56:53 2015 +0200
@@ -0,0 +1,39 @@
+function X = sylvsol(A,B,C)
+%SYLVSOL  Solve Sylvester equation.
+%   X = SYLSOL(A, B, C) is the solution to the Sylvester equation
+%   AX + XB = C, where A is m-by-m, B is n-by-n and C is m-by-n.
+%   Schur decompositions are used to convert to a (quasi)-triangular
+%   system.
+
+%   Reference:
+%   R. H. Bartels and G. W. Stewart.
+%   Algorithm 432: Solution of the matrix equation AX+XB=C.
+%   Comm. ACM, 15(9):820-826, 1972.
+
+[m,m] = size(A);
+[n,n] = size(B);
+
+realdata = (isreal(A) && isreal(B) && isreal(C));
+if ~isequal(A,triu(A)) || ~isequal(B,triu(B))
+
+   [Q, T] = schur(A,'complex');
+   [P, U] = schur(B,'complex');
+   C = Q'*C*P;
+   schur_red = 1;
+
+else
+
+   schur_red = 0;
+   T = A; U = B;
+
+end
+
+X = zeros(m,n);
+
+% Forward substitution.
+for i = 1:n
+    X(:,i) = (T + U(i,i)*eye(m)) \ (C(:,i) - X(:,1:i-1)*U(1:i-1,i));
+end
+
+if schur_red, X = Q*X*P'; end
+if realdata, X = real(X); end