%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solution to Problem 5 of
% For Your Homework
% Fast Solvers and Sylvester Equations: Both Sides Now
% Dianne P. O'Leary
%
% In this program, we compare three methods for solving Poisson's equation
% -u_{xx} - u_{yy} = f(x,y) for (x,y) in the unit square
% with u=0 on the boundary.
%
% METHOD 1:  Form the finite difference matrix with mesh-spacing h=1/n
%            and solve using sparse Cholesky factors.
% METHOD 2.  Use the Schur decomposition to solve the Sylvester
%            equation B_y U + U B_x = F. The complexity is O(n^3).
% METHOD 3:  Use the fact that the 2nd derivative operators have
%            an eigenvector matrix equal to the matrix in the discrete
%            sine transform to solve the Sylvester equation.  
%            The complexity is O(n^2 log_2 n) when n is a power of 2.
%
% Compare the times and the accuracies for n=2^p, p=2:9.
% problem5.m  Dianne O'Leary 08/2005
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

p = 2:9;
nrange = 2.^p;
k=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the solutions for each value of n.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for n = nrange

   k = k + 1;
   [K,rhs,h] = laplace2d(n);

   u = rand(n^2,1);
   U = reshape(u,n,n);
   b = K*u;
   B = reshape(b,n,n);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% METHOD 1:
% Compute the solution using sparse Cholesky, implemented in
% Matlab's backslash command, saving 
%          the elapsed time (t1),
%          the norm of the error divided by n (e1) 
%      and the residual norm divided by n (r1).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   tic
   uc = K \ b;
   t1(k) = toc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% METHOD 2:
% Compute the solution using the Schur decomposition to
%    solve the Sylvester equation, saving
%          the elapsed time (t2),
%          the norm of the error divided by n (e2)
%      and the residual norm divided by n (r2).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   tic
   e = ones(n,1)*(n+1)^2;
   T = full(spdiags([-e 2*e -e], -1:1, n, n));
   Us = sylvesterschur(T,T,B);
   t2(k) = toc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% METHOD 3:
%  Compute the solution using the fast sin transform, implemented in
%  helmsolve, saving
%          the elapsed time (t3),
%          the norm of the error divided by n (e3) 
%      and the residual norm divided by n (r3).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   tic
   Uc = helmsolve(0,B);
   t3(k) = toc;
   
   e1(k) = norm(u-uc)/n;
   e2(k) = norm(U-Us,'fro')/n;
   e3(k) = norm(U-Uc,'fro')/n;

   r1(k) = norm(b - K*uc)/n;
   r2(k) = norm(b - K*reshape(Us,n^2,1))/n;
   r3(k) = norm(b - K*reshape(Uc,n^2,1))/n;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot the results.  Figure 1 illustrates the timings
% and Figure 2 illustrates the accuracy.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1)
subplot(2,1,1)
plot(nrange,t1,nrange,t2,nrange,t3)
legend('Matlab backslash','Schur algorithm','Fast sin transform')
xlabel('n')
ylabel('time(sec)')
subplot(2,1,2)
semilogy(nrange,t1,nrange,t2,nrange,t3)
legend('Matlab backslash','Schur algorithm','Fast sin transform')
xlabel('n')
ylabel('time(sec)')
print -depsc time.eps

figure(2)
subplot(2,1,1)
plot(nrange,e1,nrange,e2,nrange,e3)
legend('Matlab backslash','Schur algorithm','Fast sin transform')
xlabel('n')
ylabel('(error norm) / n')
subplot(2,1,2)
plot(nrange,r1,nrange,r2,nrange,r3)
legend('Matlab backslash','Schur algorithm','Fast sin transform')
xlabel('n')
ylabel('(residual norm) / n')
print -depsc accy.eps