%% Iterative algorithm for computing a generalizad inverse
%
function [G1,G2,G3,t1,t2,t20,kit,t3,t30,lit] = GenInv(A)

[m,n]=size(A);
%% Via SVD
tic
[U,S,V]=svd(A,0);
G1=V*inv(S)*U';
t1=toc;
%A*G1*A-A

%% via an iterative procedure
% Find the method parameter alpha.
% To this end we need an estimate (lower bound) of the largest eigenvalue of A'*A

tic
AA = A'*A;
EA = sort(eig(AA));
alpha=0.95*2/(max(EA));
t20=toc;

I2=2*eye(m);
tic
G2=alpha*A';
G0=ones(size(A'));
kit = 0;
while norm(G0-G2)>1e-8
    G0 = G2;
    G2 = G0*(I2-A*G0);
    kit  = kit + 1;
end
t2=toc;

tic
% Alternative way to estimate max(EA) via Gershgorin's Theorem
% Compute the row-sum of abs(A)
cs = sum(abs(A),2);
for k=1:n
    w(k)=abs(A(:,k))'*cs;
end
wm = max(abs(w));
beta=2/wm;
t30=toc;

tic
G3=beta*A';
G0=ones(size(A'));
lit = 0;
while norm(G0-G3)>1e-8
    G0 = G3;
    G3 = G0*(I2-A*G0);
    lit  = lit + 1;
end
t3=toc;

disp('')
disp(['Time(SVD)    ' num2str(t1)])
disp(['Time(Iter-1) ' num2str(t2), ' Time prep. ' num2str(t20), ' Iter: ' int2str(kit)])
disp(['Time(Iter-2) ' num2str(t3), ' Time prep. ' num2str(t30), ' Iter: ' int2str(lit)])

return
%