%% Numerik D-MAVT, RW (Nipp), SS06
%
% This is the M-file solution to the fourth exercise sheet.
% To run, save this file to your account and ensure that the "Current directory"
% in MATLAB is the same as that where the file is saved.
%
% This file requires that model.m, rhs.m and sparsemodel.m are placed 
% in the same directory.
%
% Then to run the file type "s4" in the command window
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
function s4
% Serie 4, A3
% b) test with model problem: check if for smaller tolerances 
%    errors get smaller
disp(['Aufgabe 3b'])
clear;

disp(['k=40']);
k=40; A=sparsemodel(k); b=rhs(k); exact=A\b; x0=ones(size(b));
for tol=[1e-7]
    tic;
  [x,k,t,r]=cg(A,-b,x0,tol,1000); % cg solves A*x=-b
  time=toc;
  disp(['CPU time = '  , num2str(time)]);
  disp(['tol=',num2str(tol),' final residual=',num2str(r),...
        ' error=',num2str(norm(exact-x))]);
end
% output:
%Aufgabe 3b
%k=40
%Number of iteration for cg   72
%CPU time = 0.028176
%tol=1e-07 final residual=6.9176e-08 error=5.3456e-08

disp(['k=100']);
k=100; A=sparsemodel(k); b=rhs(k); exact=A\b; x0=ones(size(b));
for tol=[1e-7]
    tic;
  [x,k,t,r]=cg(A,-b,x0,tol,1000); % cg solves A*x=-b
  time = toc;
  disp(['CPU time = ' , num2str(time)]);
  disp(['tol=',num2str(tol),' final residual=',num2str(r),...
        ' error=',num2str(norm(exact-x))]);
end
%output
%k=100
%Number of iteration for cg   194
%CPU time = 0.49504
%tol=1e-07 final residual=8.799e-08 error=1.7019e-07

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% c) use now MATLAB routine pcg to solve A*x=b (note: other sign...)
disp(['Aufgabe 3c'])
clear;
k=40;A=sparsemodel(k); b=rhs(k);	
tic;
[x,flag,relres,iter]=pcg(A,b,1e-7,1000);
time=toc;
kappa=condest(A);		% estimate for condition number
disp(['k=',int2str(k),': ',num2str(time),' seconds, ',...
      num2str(iter),' iterations, kappa=',num2str(condest(A))]);
k=100;A=sparsemodel(k); b=rhs(k);
tic;
[x,flag,relres,iter,resvec]=pcg(A,b,1e-7,1000);
time=toc;
disp(relres);
kappa=condest(A);	
disp(['k=',int2str(k),': ',num2str(time),' seconds, ',...
      num2str(iter),' iterations, kappa=',num2str(condest(A))]);
% output:
%Aufgabe 3c
%k=40: 0.10101 seconds, 66 iterations, kappa=894.9698
%k=100: 1.7655 seconds, 167 iterations, kappa=5774.9592
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% d) same as c, but now with preconditioning (incomplete Cholesky)
disp(['Aufgabe 3d'])
clear;
k=40; A=sparsemodel(k); b=rhs(k);  
tic;
R=cholinc(A,1e-3);H=R';
[x,flag,relres,iter]=pcg(A,b,1e-7,1000,H,H');
time=toc;
kappa=condest(A);	
disp(['k=',int2str(k),': ',num2str(time),' seconds, ',...
      num2str(iter),' iterations, kappa=',num2str(condest(A))]);
k=100;A=sparsemodel(k); b=rhs(k);     
tic;
R=cholinc(A,1e-3);H=R';
[x,flag,relres,iter]=pcg(A,b,1e-7,1000,H,H');
time=toc;
kappa=condest(A);		% estimate for condition number
disp(['k=',int2str(k),': ',num2str(time),' seconds, ',...
      num2str(iter),' iterations, kappa=',num2str(condest(A))]);
% output:
%Aufgabe 3d
%k=40: 0.073602 seconds, 7 iterations, kappa=894.9698
%k=100: 1.0632 seconds, 14 iterations, kappa=5774.9592

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% e) memory test
disp(['Aufgabe 3e'])

% the sie of matrix which you are able to consider will depend on the
% RAM on your computer....
clear; format compact; k=68; A=model(k); whos; 
% output:
%   Name      Size         Bytes  Class
%
%   A      4356x4356    151797888  double array
%   k         1x1              8  double array 
%
% Grand total is 18974737 elements using 151797896 bytes
 
clear; format compact; k=69; A=model(k); whos; 
% output:
% ??? Error using ==> full
% Out of memory. Type HELP MEMORY for your options.
% 
% Error in ==> /hg/u/fab/assi/numerik_D-MAVT/solutions/s4/model.m
% On line 6  ==> A = full(delsq(numgrid('S',n)));
clear; format compact; k=299; A=sparsemodel(k); whos; 
% output:
%   Name      Size         Bytes  Class
% 
%   A     88209x88209    5631124  sparse array
%   k         1x1              8  double array
% 
% Grand total is 439858 elements using 5631132 bytes
% 
clear; format compact; k=300; A=sparsemodel(k); whos; 
% output:
% 
% ??? Error using ==> sparse
% Out of memory. Type HELP MEMORY for your options.
% 
% Error in ==> /hg/s/solaris/2.5.1/app/matlab/5.3/toolbox/matlab/demos/delsq.m
% On line 34  ==> D = sparse(i,j,s);
% 
% Error in ==> /hg/u/fab/assi/numerik_D-MAVT/solutions/s4/sparsemodel.m
% On line 7  ==> A = delsq(numgrid('S',n));
return;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Serie 4, A3
% a) implement cg-method
function [u,k,time,resid]=cg(A,b,u0,tol,maxIter)

% conjugate gradients to solve A*u+b=0 (or A*u=-b, mind the sign ;-))
% in:	A 	matrix of linear system
%	b	(negative) right hand side of system
% 	u0	initial guess
% 	tol 	tolerance: if norm(residual)<tol ==> stop
%	maxIter	maximal number of iterations
% out:	u	approximate solution
%	k	number of iterations needed
%	time	cpu time needed
%	resid	2-norm of residual for final approximativ solution
% auth: fab@math.ethz.ch

tic;		% start measuring time
% initialisation
number_iteration=0;
u=u0;
rm1=A*u+b; p=-rm1; n2m1=rm1'*rm1; % rm1=r^{k-1}, n2m1=<rm1,rm1>
for k=1:maxIter
  if k>=2
    e = n2m1/n2m2;
    p = -rm1+e*p;
  end
  z    = A*p;
  rho  = n2m1/(p'*z);
  u    = u + rho*p;
  rm2  = rm1;
  rm1  = rm1+rho*z;
  n2m2 = n2m1;
  n2m1 = rm1'*rm1;
  number_iteration = number_iteration+1;
  if (n2m1<tol*tol) break; end;	% converged!
end
time=toc;	% stop measuring time
if (k==maxIter) 
  error(['not converged after ',int2str(maxIter),' iterations']); 
end
disp(['Number of iteration for cg   ', int2str(  number_iteration)]);
resid=norm(rm1);
return;



