%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Solution to Problem 4 in the homework assignment,
%  Updating and Downdating Matrix Factorizations:
%     A Change in Plans
%
%  For matrices of various sizes n, find the number of updates
%  that makes the cost of using the Sherman-Morrison-Woodbury
%  formula to solve a modified linear system comparable to 
%  solving the problem with backslash.
%
%  problem4.m       Dianne P. O'Leary  12/2005
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nn = [10:10:50,100:100:1000];

ntries = 5;   % Number of repeats, used to try to make the timings
              % more reproducible.
k = 0;

for nsize=1:length(nn);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Generate a random problem of size n.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   n = nn(nsize);
   A = rand(n,n);
   x = rand(n,1);
   b = A*x;
   [L,U,P]=lu(A);
   P = sparse(P);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Set timefact to be the time for solving using backslash.
%  (Ignore the overhead in the "for" loop.)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   tic
   for j=1:ntries,
      xcomp = A \ b;
   end
   timefact(nsize) = toc/ntries;
   ttime = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Use Sherman-Morrison-Woodbury with a rank k update, and increase
%  k until the time is greater than that for backslash.
%  (Ignore the overhead in the "for" loop.)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   while (ttime < timefact(nsize))&(k < n),
     k = k + 1;
     Z = rand(n,k);
     V = rand(n,k);

     tic
     for j=1:ntries
        xnew = sherman_mw(L,U,b,Z,V,P);
     end
     ttime = toc/ntries;
   end
   cross(nsize) = k;
   k = k - 1;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Plot the results.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

plot(nn,cross,'LineWidth',1.5)
xlabel('Size of matrix')
ylabel('Rank of update')
title('Making Sherman-Morrison-Woodbury time comparable to Backslash')
print -dtiff smw.tif