function [X,Y] = sda(E,F,G,P)
% [X,Y]=SDA(E,F,G,P) applies the Structured Doubling Algorithm starting
% from E, F, G, P
% E, F, G, P: initial matrices
% X: limit of the sequence (P_k)
% Y: limit of the sequence (G_k)
tol = 1e-13;
kmax = 30;
err = 1;
k = 0;
n = size(G,1); m = size(P,1);
while err > tol && k < kmax
IGP = eye(n) - G*P;
IPG = eye(m) - P*G;
E1 = E/IGP;
F1 = F/IPG;
G = G + E1*G*F;
P = P + F1*P*E;
E = E1*E;
F = F1*F;
err = min(norm(E,1),norm(F,1));
k = k + 1;
end
X = P; Y = G;
if k == kmax
disp('Warning: reached the maximum number of iterations')
end