function X = liap(A, B, C)
%   X = LIAP(A,B,C) solves the general form of the Liapunov matrix
%   equation (also called Sylvester equation):
%   (A and B must be square and C must have the rows of A and columns of B)
%
%           A*X + X*B + C = 0
%
%   Q = LIAP(A,B) solves the Right Liapunov matrix equation:
%       (A and B have the same number of rows.)
%
%           A*Q + Q*A' + B*B' = 0
%
%   P = LIAP(A',C) solves the Left Liapunov matrix equation:
%       (A and C have the same number of columns.)
%
%           A'*P + P*A + C'*C = 0

ni = nargin;

if (ni ~= 3 && ni ~= 2)
    print_usage ();
end

[ma,na] = size(A);
if ( ma ~= na )
    error ('liap: A is not square');
else
     n = ma;
end

if (ni == 2)  %% Transform Liapunov equation to Sylvester equation form.

    [mb,nb] = size(B);
    if ( mb ~= nb )             % b is not square
      if (mb == n)              % solve a x + x a' + b b' = 0   (right)
        C = -B * B';
        %disp('liap: right LE')
      elseif (nb == n)          % solve a'x + x a + b' b = 0     (left)
        C = -B' * B;
        %disp('liap: left LE')
      else
        error ('liap: A, B not conformably dimensioned');
      end
    else                        % b is square
      C = -B;
    end
    B =  A';
end

% Solve the Sylvester equation.  A X + X B = C
X = sylvester(A,B,C);
% find(abs(A*X + X*B - C) > 1e-14)   %% test

% Ignore complex part if real inputs (imag parts better be small)
if isreal(A) && isreal(B) && isreal(C),
   X = real(X);
end
