%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solution:
% Your Homework Assignment
% Iterative Methods for Linear Systems:
% Following the Meandering Way
% 
% This program solves a set of linear systems of various sizes using
% several methods:
% 
%    -- Cholesky factorization
%    -- Gauss-Seidel
%    -- Conjugate gradients
%    -- Preconditioned Conjugate gradients, with preconditioner
%       -- symmetric Gauss-Seidel
%       -- ICCG, with various drop tolerances
% 
% and two orderings:
%
%    -- Original
%    -- Approximate minimum degree (AMD)
%
% The results (time and storage) are compared.
%
% We use two problems, problem1 with 7 refines and problem 2 with 4 refines.
% 
% Dianne P. O'Leary 04/2006
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all
nrefine_array = [7,4];
cut = [.05 .5];
kappa = 0;
maxit = 100;

for myproblem=1:2,

   clear itns timest nnzR n mditns mdtimest
   disp(sprintf('Problem %d',myproblem))
   nrefine = nrefine_array(myproblem);
   mesh = generateproblem(myproblem,nrefine,kappa);
   if myproblem==1
      refines = 2:nrefine;
   else
      refines = 1:nrefine;
   end
   
   for nrefine=refines
      A = mesh(nrefine).A;
      b = mesh(nrefine).b;
      tol = mesh(nrefine).tol;
      [n(nrefine),m]=size(A);
      disp(sprintf('ORIGINAL ORDERING'))
      [itns(nrefine,:),timest(nrefine,:),nnzR] ...
               = runmethods(A,b,tol,cut,maxit,'yes');
   
      p = symamd(A);
      Ap = A(p,p);
      bp = b(p);
      disp(sprintf('AMD ORDERING'))
      [mditns(nrefine,:),mdtimest(nrefine,:),mdnnzR] ...
               = runmethods(Ap,bp,tol,cut,maxit,'yes');
   
   end
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot the results in four figures.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   mysymb='.bogxr+c*mskdbvg^r<c>mpkhb.gorxc+m*ksbdg';
   legenders = 'Chol GS CGGS  cut=.05cut=.5';
   nmethods = length(itns);
   myxlabel = 'Number of unknowns';
   myylabeli = 'Number of iterations';
   myylabelt = 'seconds';
   
   mytitle =  ...
        sprintf('Prob %d: Number of iterations, Original ordering',myproblem);
   makeplot(n,itns,mytitle,myxlabel,myylabeli)
   filename = sprintf('orig%ditns.eps',myproblem);
   feval(@print,'-depsc',filename)

   mytitle = sprintf('Prob %d: Number of iterations, AMD ordering',myproblem);
   makeplot(n,mditns,mytitle,myxlabel,myylabeli)
   filename = sprintf('amd%ditns.eps',myproblem);
   feval(@print,'-depsc',filename)

   mytitle=sprintf('Prob %d: Time, Original ordering',myproblem);
   makeplot(n,timest,mytitle,myxlabel,myylabelt)
   filename = sprintf('orig%dtime.eps',myproblem);
   feval(@print,'-depsc',filename)

   mytitle=sprintf('Prob %d: Time, AMD ordering',myproblem);
   makeplot(n,mdtimest,mytitle,myxlabel,myylabelt)
   filename = sprintf('amd%dtime.eps',myproblem);
   feval(@print,'-depsc',filename)
   
end