%------------------------------------------------------------------------------------
% A throw-away script to plot the results of using Bi-CGstab on a second order PDE
% discretized via 5-point centered differences. 
%
% Data was extracted from runs on Laplace operator on uniform  N x N mesh on [0, 1]
% x [0, 1]. I just used grep and vim to extract the data from the output files named
% "results", then inserted it here.
%
% The time per iteration is computed using the matrix-vector product time, not the
% iteration time. Since generally discretization and initialization is trivial, on the
% order of less than 0.2%, that's good enough. I also tried a polynomial fit with a
% quadratic since the time per iteration should be quadratic in N or linear in number
% of unknowns = N^2.
%
% This opens up larger than usual graphics windows, slightly staggered in location.
%
%-----------------
% Randall Bramley
% Department of Computer Science
% Indiana University, Bloomington
%------------------------------------------------------------------------------------

A = [ ...
1000      4.690678E+01
2000      2.406304E+02
3000      9.576622E+02
4000      2.329062E+03
5000      6.308759E+03
6000      6.775227E+03
7000      1.197025E+04
8000      1.578574E+04];  % That is over 4.4 hours!

N = A(:,1);
t = A(:,2);  % total time

its = [ ...
    1512
    1897
    3398
    4625
    8118
    6117
    8058
    8134];

% Time in the matrix-vector products
matvec = [...
   3.321102E+01
   1.672217E+02
   6.636362E+02
   1.618257E+03
   4.374013E+03
   4.673960E+03
   8.197807E+03
   1.080435E+04];

% Time in the iterative phase of Bi-CGstab. Includes the mat-vec product time.
ittime = [ ...
   4.684606E+01
   2.404780E+02
   9.573645E+02
   2.328537E+03
   6.308012E+03
   6.774172E+03
   1.196883E+04
   1.578391E+04];

figure('Units', 'normalized', 'Position', [0.15 0.15 0.5 0.5]);
semilogy(N, t/60, 'r-', N, t/60, 'b*')
xlabel('Mesh size parameter')
ylabel('Total time (minutes)')
title('PDE Timings')
grid minor

figure('Units', 'normalized', 'Position', [0.1 0.1 0.5 0.5]);
plot(N, its, 'r-', N, its, 'b*')
xlabel('Mesh size parameter')
ylabel('Number of iterations')
title('Bicgstab Number of Iterations')

%-------------------------------------------------
% Time per iteration based on matrix-vector time.
%-------------------------------------------------
% Notice the ./ instead of / ; that makes Matlab use componentwise division.
tpi = matvec./its;
% Try a quadratic fit to tpi; it should be a good match:
q = polyfit(N, tpi, 2);
% The polynomial has coefficients in decreasing power order, viz.,
%       s(N) = q(1)*N^2 + q(2)*N + q(3)
s = polyval(q, N);   % Evaluate the polynomial at the sizes in N
% As usual with a polynomial fit, the constant term q(3) may be negative. While 
% that makes no physical sense, it is still the best approximant for the given data.
% If you want to see what the coefficients are, just type "q" without a semicolon
% at the command window prompt.


figure('Units', 'normalized', 'Position', [0.05 0.05 0.5 0.5]);
% The order of the 3 data sets is critical to make the legend correct.
plot(N, tpi, 'b*', N, s, 'go-', N, tpi, 'r-')
xlabel('Mesh size parameter')
ylabel('time (s)')
title('Bicgstab Time per Iteration')
legend('TPI', 'Quad fit', 'Location', 'NorthWest')

%-------------------------------------------------------
% What proportion of the time is spent in the mat-vec?
%-------------------------------------------------------
p = 100*matvec./ittime;
figure('Units', 'normalized', 'Position', [0.01 0.01 0.5 0.5]);
plot(N, p, 'b*-')
xlabel('Mesh size parameter')
ylabel('Percent time')
title('Percent Time in Matrix-Vector products')
% This one looks like the percentage drops rapidly. But notice the scale of the
% y-axis. A more illuminating picture occurs if the axes are changed:
a = axis;
a(3) = 0; a(4) = 100;  % Make the y-axis go from 0% to 100%
axis(a);




