Mercurial > matrix-functions
view mftoolbox/sylvsol.m @ 8:a587712dcf5f draft default tip
funm_atom.m: rename fun_atom to funm_atom
* funm_atom.m: rename fun_atom to funm_atom.
author | Antonio Pino Robles <data.script93@gmail.com> |
---|---|
date | Fri, 29 May 2015 09:48:36 +0200 |
parents | 8f23314345f4 |
children |
line wrap: on
line source
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