format short e; format compact
%Es 1(a)(b)
fprintf('Es 1\n')

for n=100:50:500,
    I=speye(n);
    e=ones(n,1);
   %A=spdiags([-e,-e,-e,4*e,-e,-e,-e], -3:3, n,n);
    A=spdiags([-e,2*e,-e], -1:1, n,n);
    c=ones(n,1); C=c*c';
    tic;
    XL=lyap(A,-C);
    tlyap=toc;
    AA=kron(I,A)+kron(A,I);
    tic
    xD=AA\C(:);
    tdir=toc;
    disp([n,norm(XL(:)-xD)/norm(xD), tlyap,tdir])
end


fprintf('Es 2\n')
% Es.2
% WARNING: this script is only for illustration. In particular:
% 1) the approximation to the solution SHOULD BE determined
% during the construction of the approximation space
% so as to stop building the space as soon as the
% stopping criterion is satisfied
% 2) The approx solution X and the residual matrix SHOULD NOT
% be explicitly computed. See 
% www.dm.unibo.it/~simoncin/software.html for proper implementations
%
tol=1e-8;
m_max=80;
n=50;
h=1/(n-1);
e=ones(n,1);
ee=ones(n*n,1);
T=-1/h^2*spdiags([-e,2*e,-e], -1:1, n,n);
B=1/(2*h)*spdiags([e,0*e,-e], -1:1, n,n);
I=speye(n);
C=ee*ee';
normC=norm(C,'fro');
II=speye(n*n);
for imethod=1:3,
   if imethod==1,    %standard Krylov
       fprintf(' standard Krylov\n')
   elseif imethod==2,  % extended Krylov
       fprintf(' Extended Krylov\n')
   else
       fprintf(' Rational Krylov\n')
   end
   for alpha=[0,50,200],
     A=kron(I,T+alpha*B)+kron(T,I);
     if imethod==1,    %standard Krylov
       [RV]=block_Krylov(A,ee,m_max);
     elseif imethod==2,  % extended Krylov
       [RV]=eksm(A,ee,m_max);
     else
       [RV]=rksm(A,ee,m_max,1e-9,1e-2,10);
     end
     for k=1:m_max,
      H = RV(:,1:k)'*A*RV(:,1:k);
      Cm= RV(:,1:k)'*C*RV(:,1:k);
      Y = lyap(H,-Cm);
      X = RV(:,1:k)*Y*RV(:,1:k)';   % this is resource-expensive, should not be done
      norm_res=norm(C-A*X-X*A','fro')/normC;  % this is resource-expensive, should not be done
      if norm_res<tol,
      disp([alpha,k,norm_res])
      break,
      end
   end
   if k>=m_max, disp([alpha,k,norm_res]),end
end
end


% Es 3
fprintf('Es 3\n')
for n=100:50:500,
    u=randn(n,1);
    v=randn(n,1);
    I=speye(n);
    e=ones(n,1);
    A=spdiags([-e,2*e,-e], -1:1, n,n);
    c=ones(n,1); C=c*c';
    XC=lyap(A,-C);
    XU=lyap(A,-u*u');
    gamma = v'*XC*v/(1+v'*XU*v);
    X=XC-gamma*XU;
    disp([n, norm(A*X+X*A'+u*v'*X*v*u'-C,'fro')/norm(C,'fro')])
end

