%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  timetest: a quicky matlab script to calculate the scaling of matlab
% solves for tridiagonal problems x=A\b and x=inv(A)*b;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
for i=1:8
    N(i)=100*2^i; % number of points
    disp(sprintf('problem size=%d points\n',N(i)));
    A=setA(N(i)); % set the matrix
    r=rand(N(i),1); % generate a random rhs;
    % calculate x by gaussian elimination
    tic; x=A\r;  t_gauss(i)=toc;
    disp(sprintf('\tGaussian elimination=%g s',t_gauss(i)));
    % calculate x by LU decomposition 
    tic; [L,U]=lu(A); x=U\(L\r);  t_lu(i)=toc;
    disp(sprintf('\tLU Factorization=%g s',t_lu(i)));
    if ( i < 6 )
    % calculate x by direct inversion of A
        tic; x=inv(A)*r; t_inv(i)=toc;
        disp(sprintf('\tDirect inversion=%g s',t_inv(i)));
    else
        t_inv(i)=NaN;
    end
    disp(sprintf('\tmax/min ratio=%g\n',nanmax([t_inv(i) t_lu(i) t_gauss(i)])/nanmin([t_inv(i) t_lu(i) t_gauss(i)])));
    clear A x r;
end
figure
h=loglog(N,t_gauss,'r-o',N,t_lu,'g-o',N,t_inv,'b-o');
set(h,'linewidth',2);
xlabel('Vector Length N');
ylabel('Solution time (s)');
legend('GE','LU','INV');