%
%  Apply GMRES to a example problems
%
%  Matthias Heinkenschloss
%  February 3, 2012 (last modified February 7, 2017)



set(0, 'defaultaxesfontsize',24,'defaultaxeslinewidth',2.0,...
       'defaultlinelinewidth',2.0,'defaultpatchlinewidth',2.0,...
       'defaulttextfontsize',24);

 

%% Random Matrix Example
fprintf('\n\n Random Matrix. Hit return to continue ....\n');
pause
n = 20;
A = rand(n,n);
b = rand(n,1);

[x, res, iter, flag] = mygmres_simple( A, b, [], [] );
[x, res, iter, flag] = mygmres( A, b, [], [] );


figure(1)
fprintf(' GMRES return after %d iterations with flag = %d\n', iter,flag)
semilogy((0:iter),res/res(1))
xlabel('Iteration')
ylabel('GMRES Residual')


%% 1D advection diffusion equation
fprintf('\n\n 1D Advection Diffusion Equation. Hit return to continue ....\n');
pause

% apply GMRES to compute an approximate solution of the FD discretization of the BVP
%
%     - alpha  u''(x) + beta u'(x) + gamma u(x) = f(x)  x in [a,b],
%                          u(a) = r_a,
%                          u(b) = r_b,
% on [a,b] = [0,1]
%


% set problem data
alpha = 1.e-3;
beta  = 1;
gamma = 0;
ra    = 0;
rb    = 0;


f     = inline( '3*(x.^2 - 2*alpha*x)', 'x', 'alpha', 'beta', 'gamma');
uex   = inline( 'x.^3 - (exp((x-1)/alpha)-exp(-1/alpha))./(1 - exp(-1/alpha))', 'x', 'alpha', 'beta', 'gamma');


% set mesh points
n   = 39;
dx  = 1/(n+1);
x   = [0:dx:1]';


% Use upwind finite difference approximation
if( beta < 0 ); error( 'beta < 0 not yet implemented'); end

% set subdiagonal
c    = -(alpha + dx*beta)*ones(n,1);
% set diagonal
d    = (2*alpha + dx*beta + dx*dx*gamma)*ones(n,1);
% set superdiagonal
e    = (-alpha)*ones(n,1);

% generate system matrix
A = spdiags([c d e], -1:1, n, n);

% right hand side
b      = dx*dx*f(x(2:n+1), alpha, beta, gamma);
b(1)   = b(1) + (alpha + 0.5*dx*beta)*dx*ra;
b(n)   = b(n) + (alpha - 0.5*dx*beta)*dx*rb;

% FD  solution
ufd = A\b;


x0 = zeros(n,1);
linestyles = [{'-k'},{':r'},{'--b'},{'-.g'},{':m'}]; 
m          = [2;5;10;20;2*n];
for i = 1:5
    options.restart = m(i);
    options.maxit   = n;
    
    [u, res, iter, flag] = mygmres( A, b, [], [], x0, options );
    
    figure(1)
    fprintf(' GMRES(%d) return after %d iterations with flag = %d\n', m(i), iter,flag)
    semilogy((0:iter),res/res(1),linestyles{i}); hold on
end
xlabel('Iteration')
ylabel('GMRES Residual (Normalized)')
legend('m=2','m=5','m=10','m=20','full','Location','SouthWest')
hold off
%print('-depsc','gmres_bvp.eps')

% plot solution
figure(2)
xx = (0:1/100:1);
plot(x,[ra;ufd;rb],'k-'); hold on
plot(x,[ra;u;rb],'r--'); hold on
plot(xx,uex(xx, alpha, beta, gamma),'b:'); hold off
xlabel('x')
ylabel('u(x)')
legend('FD approximation','GMRES solution', 'exact BVP solution','Location','NorthWest')
axis([-0.1,1.1,-0.1,1.1])
%print('-depsc',filename_sol)


%% Embree Example
fprintf('\n\n Embree Example. Hit return to continue ....\n');
pause
A  = [ 1 1 1; 0 1 3; 0 0 1];
b  = [ 2; -4; 1];
x0 = zeros(3,1);

linestyles = [{'-k'},{':r'},{'--b'}]; 
for m = 1:3
    options.restart = m;
    options.maxit   = 30;
    
    [x, res, iter, flag] = mygmres( A, b, [], [], x0, options );
    
    figure(1)
    fprintf(' GMRES(%d) return after %d iterations with flag = %d\n', m, iter,flag)
    semilogy((0:iter),res/res(1),linestyles{m}); hold on
end
xlabel('Iteration')
ylabel('GMRES Residual (Normalized)')
legend('m=1','m=2','m=3','Location','SouthEast')
hold off

    