function [L, kstep] = smith(A, B, mu, tol)
%
% This is a simple function using the modified Smith method to
% approximately compute a low-rank approximation to the solution of the
% Lyapunov equation
%
%                  ALL' + LL'A' = BB',
%
% where A is a real n-by-n matrix
%       B is a real n-by-p matrix
%       mu is a negative real shift
%       tol is a stopping tolerance
%
% R. Nong, 2009-04-01
%

stol = sqrt(eps); % for SVD truncation
n = size(A, 1);

rho = sqrt(2 * abs(mu));
Bk = (A + mu * eye(n)) \ B;
Bk = rho * Bk;

[V, S, W] = svd(Bk, 0);
s = diag(S);
m = length(s);
k = 1;
while (k < m & s(k + 1) > stol * s(1))
    k = k + 1;
end
V = V(:, 1:k);
S = diag(s(1:k));

% there are different ways to choose the number of steps
specRad = max(abs(eig((A + mu * eye(n)) \ (A - mu * eye(n)))));
kstep = ceil(log10(tol) / log(specRad));

for i = 1:kstep
    Bk = (A + mu * eye(n)) \ ((A - mu * eye(n)) * Bk);
    G = V' * Bk;
    Vhat = Bk - V * G;
        % DGKS step
        C = V' * Vhat;
        G = G + C;
        Vhat = Vhat - V * C;
    [Vhat, R] = qr(Vhat, 0);
    m2 = size(S, 2);
    m1 = size(R, 1);
    [T, S, Y] = svd([S G; zeros(m1, m2) R], 0);
    s = diag(S);
    m = length(s);
    k = 1;
    while(k < m & s(k + 1) > stol * s(1))
        k = k + 1;
    end
    V = [V Vhat] * T(:, 1:k);
    S = diag(s(1:k));
end

L = V * S;