%          Test script for Jacobi and GS
% First define: method = 'Jacobi' or: method = 'GS' .
% Also define: myname = 'xxxxx xxxxx' , when you are 
% ready to print pictures.
% ------------------------------------------------------

if ~exist('method','var'); method = 'Jacobi'; end

fprintf('\n====== CAAM554/454 Jacobi/GS Tests ======\n')
format1 = 'rres = %8.4e  err = %8.4e  cputime = %5.2f  \n'; 
format2 = 'r1 = %8.4e  r2 = %5.2f  \n'; 
fprintf(['\nMethod is ' method '\n\n'])

N = input(' N = ');
A = gallery('poisson',N);
n = size(A,1);
%xs = numgrid('H',N) > 0;
xs = ones(n,1);
b = A*xs(:); 

tol = 1.e-4;
fprintf('\n----- Poisson Matrix, size = %g -----\n',n)  

% ------ run my code ------
   fprintf('Mycode: ')
   t = cputime;
   [myx, myinfo, myrres] = feval(['my' method],A,b,10*n,tol);
   mycput = cputime - t;
   myerr = norm(myx-xs)/norm(xs);
   fprintf(format1,myrres(end),myerr,mycput);

% ------ run instructor's code ------
   fprintf('YZcode: ')
   t = cputime;
   [zx, yzinfo, zrres] = feval(['z' method],A,b,10*n,tol);
   zcput = cputime - t;
   zerr = norm(myx-xs)/norm(xs);
   fprintf(format1,zrres(end),zerr,zcput);

% ------ calculate ratios ------
   fprintf('ratios: ')
   r1 = myerr/zerr;
   r2 = mycput/zcput;
   fprintf(format2,r1,r2);
   disp(['Myinfo: ' myinfo]);
   disp(['YZinfo: ' yzinfo]);
   
   str1 = ['ERR ratio = ' num2str(r1) '.  ']; 
   str2 = ['CPU ratio = ' num2str(r2) '.  '];
   str3 = date;
   if exist('myname','var')
      str3 = [myname '. ' str3];      
   end

% ------ plotting ------
m = length(myrres);
figure(1)
semilogy(1:m, myrres, '.r');
xlabel('Iteration');  
ylabel('Residual')
title([method ' Method Residual Profile'])
a = axis; 
tx = (a(1) + a(2))/2; ty = a(4);
text(tx,ty/2,str1,'fontsize',14); 
text(tx,ty/4,str2,'fontsize',14); 
text(tx,ty/8,str3,'fontsize',14); 
set(gca,'fontsize',16)

figure(2); surfc(reshape(myx,[N N]))
title('My Solution')
set(gca,'fontsize',16)
fprintf('\n')
