function results_table=SOL_Take_Home
% SOL_TAKE_HOME
% Do the take-home part of the exam with Tim's gmres code.
% The gmres code is old and does not offer Gram-Schmidt twice
% as an option. For this exercise it does not matter.
%
% List of demensions, step sizes, and solver parameters.
%
harray=[2^-5, 2^-6, 2^-7];
narray=[31, 63, 127];
tol=1.d-6;
maxit=1000;
params=[tol, maxit, 1];
results=zeros(3,7);
%
% Sweep through the problems.
%
for ip=1:3
    h=harray(ip);
    n=narray(ip);
    u0=zeros(n*n,1);
%
% Make a column vector of grid points.
%
    h=1/(n+1); x=h:h:1-h; x=x';
%
% Create the right side on the n x n geometric grid.
%
    s=sin(pi*x); c=cos(pi*x);
    f=2*pi*pi*s*s' + pi*(c*s' - s*c');
%
% Do the same for the solution. Then map it to a linear array.
%
    ustar2=s*s'; ustar1=reshape(ustar2,n*n,1);
%
% Map the right side to a linear array.
%
    rhs=reshape(f,n*n,1);
%
% I'm using left preconditioning, so I must precondition rhs for 
% the preconditioned solve.
%
    prhs=P2D(rhs);
%
%  Unpreconditioned solve.   
%
    tic;
    [sol, error, total_iters] = gmres(u0, rhs, @pdeop, params);
    uptime=toc;
    delu=norm(sol-ustar1,inf);
%
%  Preconditioned solve.   
%
    tic;
    [psol, perror, ptotal_iters] = gmres(u0, prhs, @ppdeop, params);
    ptime=toc;
    pdelu=norm(psol-ustar1,inf);
%
%   Put the results in a neat and tidy array.
%  
    results_table(ip,:)=[h, total_iters, delu, uptime, ...
                      ptotal_iters, pdelu, ptime];
end
%
% Now format the table in LaTeX for cutting/pasting into the homework report.
%
headers={'h','its','err','time','its','err','time'};
labels=' & \multicolumn{3}{|c|}{Unpreconditioned} & \multicolumn{3}{|c|}{Preconditioned}';
fprintf('\\begin{tabular}{|l|lll|lll|}\n');
fprintf('\\hline \n');
fprintf('%s \\\\ \n',labels);
fprintf('\\hline \n');
for i=1:6
   fprintf('%8s &',headers{i});
end
   fprintf('%8s ',headers{7});
fprintf('\\\\ \n');
fprintf('\\hline \n');
for i=1:3
 fprintf('%8.2e & %6d & %10.2e & %10.2e & %6d & %10.2e & %10.2e \\\\ \n', ...
          results_table(i,:) )
end
fprintf('\\hline \n');
fprintf('\\end{tabular}');


function lu=pdeop(u)
%
% Unpreconditioned operator.
%
lu=lapmf(u);
lu=lu + dxmf(u) - dymf(u);

function plu=ppdeop(u)
%
% Left preconditioned operator.
%
lu=lapmf(u);
lu=lu + dxmf(u) - dymf(u);
plu=P2D(lu);


function dxu = dxmf(u)
%
% matrix-free partial derivative wrt x
% homogeneous Dirichlet BC
%
n2=length(u);
n=sqrt(n2);
h=.5*(n+1);
%
% turn u into a 2D array with the BCs built in
%
uu=zeros(n+2,n+2);
vv=zeros(n,n);
vv(:)=u;
uu(2:n+1,2:n+1)=vv;
%
% compute the partial
%
dxuu=zeros(n,n);
dxuu=uu(3:n+2,2:n+1)-uu(1:n,2:n+1);
%
% fix it up
%
dxu=h*dxuu(:);


function dyu = dymf(u)
%
% matrix-free partial derivative wrt y
% homogeneous Dirichlet BC
%
n2=length(u);
n=sqrt(n2);
h=.5*(n+1);
%
% turn u into a 2D array with the BCs built in
%
uu=zeros(n+2,n+2);
vv=zeros(n,n);
vv(:)=u;
uu(2:n+1,2:n+1)=vv;
%
% compute the partial
%
dyuu=zeros(n,n);
dyuu=uu(2:n+1,3:n+2)-uu(2:n+1,1:n);
%
% fix it up
%
dyu=h*dyuu(:);