sol8.tar

ex8_1.m0000644000175100017510000000226013203315137011504 0ustar  arbenzarbenzfunction ex8_1()
%FEM17 - Exercise 8.1: Steepest descent for solving Ax=b with SPD A
%
% Nov 2017

g = 'squareg';
bond = 'squareb1';
c = 1; a = 0; f = 1;

[p,e,t] = initmesh(g,'hmax',1);  % generate a coarse mesh
[p,e,t] = refinemesh(g,p,e,t);  % refine the mesh
[p,e,t] = refinemesh(g,p,e,t);  % refine the mesh
[p,e,t] = refinemesh(g,p,e,t);  % refine the mesh
[p,e,t] = refinemesh(g,p,e,t);  % refine the mesh

% check the number of nodes by size(p,2), it should be between 1000 to
% 3000.

[A,b] = assempde(bond,p,e,t,c,a,f);

D = diag(diag(A));  % diagonal of A
L = tril(A,-1);     % strictly lower triangle of A

x0 = zeros(size(A,1),1); tol = 1e-6; maxit = 20000;
[x,nit] = sd(A,[],[],b,x0,tol,maxit);     % unpreconditioned steepest descent
fprintf('\nUnpreconditioned steepest descent: num. of its = %d\n',nit)

[x,nit] = sd(A,D,[],b,x0,tol,maxit); % Jacobi preconditioned steepest descent
fprintf('\nJacobi preconditioned steepest descent: num. of its = %d\n\n',nit)

for omega=[1:.1:1.9]
   M=sqrt(omega/(2-omega))*(D/omega+L)/sqrt(D);
   [x,nit] = sd(A,M,M',b,x0,tol,maxit);
   fprintf(['SSOR(%4.2f) preconditioned steepest descent: ' ...
            'num. of its = %d\n'],omega,nit)
end
ex8_2.m0000644000175100017510000000220713203310346011503 0ustar  arbenzarbenzfunction ex8_2()
% Exercise 8.2: Conjugate gradient method for solving Ax=b with SPD A
%
% Nov 2017

g = 'squareg';
bond = 'squareb1';
c = 1; a = 0; f = 1;

[p,e,t] = initmesh(g,'hmax',1);  % generate a coarse mesh
[p,e,t] = refinemesh(g,p,e,t);  % refine the mesh
[p,e,t] = refinemesh(g,p,e,t);  % refine the mesh
[p,e,t] = refinemesh(g,p,e,t);  % refine the mesh
[p,e,t] = refinemesh(g,p,e,t);  % refine the mesh

% check the number of nodes by size(p,2), it should be between 1000 to
% 3000.

[A,b] = assempde(bond,p,e,t,c,a,f);

D = diag(diag(A));  % diagonal of A
L = tril(A,-1);     % strictly lower triangle of A

x0 = zeros(size(A,1),1); tol = 1e-6;
[x,nit] = mypcg(A,[],[],b,x0,tol);   % unpreconditioned CG algorithm
fprintf('\nUnpreconditioned CG algorithm: num. of its = %d\n',nit)

[x,nit] = mypcg(A,D,[],b,x0,tol);    % Jacobi preconditioned CG algorithm
fprintf('\nJacobi preconditioned CG algorithm: num. of its = %d\n\n',nit)

for omega=[1:.1:1.9]
   M=sqrt(omega/(2-omega))*(D/omega+L)/sqrt(D);
   [x,nit] = mypcg(A,M,M',b,x0,tol);
   fprintf(['SSOR(%4.2f) preconditioned CG algorithm: ' ...
            'num. of its = %d\n'],omega,nit)
end
ex8_3.m0000644000175100017510000000614313203531755011520 0ustar  arbenzarbenz% FEM17 - Exercise 8-3
% 
% Study the convergence rates of the steepest descent method and CG method
%
% Remarks:
%
% 1. PCG usually converges faster than the SD method.
%
% 2. For the ill-conditioned matrix, i.e. the corresponding condition number
% is large, then the convergence is slow for both SD and PCG method.  the
% condition number "kappa=lambda_max/lambda_min".  In this exercise,
% kappa=10^p.
%  
% 3. Usually, both SD and PCG solver converge faster for the matrix with
% clustered eigenvalue distribution than with equally eigenvalue
% distribution. Here is a reason:
%
% According to the theorem in slide 21 of lecture 8, at step m the residual
% r^(m) is bounded by the value of |p_m(lambda)|, where p_m is a polynomial
% function with degree of m.  Usually m is much smaller than N (A is a
% N-by-N matrix).
% 
% In the case of equally distribution of eigenvalues, it's more difficult to
% find a p_m to be small evaluated on all eigenvalues. That is because the
% eigenvalues are apart from each other. According to the theorem, the
% residual r^(m) may also become large, meaning that solver may converge
% slower.
%
% In the case of clustered distribution of eigenvalues, it's much easier to
% make p_m to be small evaluated on all eigenvalues. In fact in the
% clustered region, if p_m is zero evaluated for one eigenvalue, then for
% all the other eigenvalues nearby, p_m is close to zero.  According to the
% theorem, the residual r^(m) may be bounded and small and the solver
% converges faster.
%
% Hua Guo, 19th Nov. 2010.


N=100;
rand('seed',100);
b=rand(N,1);
tol=1e-6;
i=[0:N-1]';
for p=1:3
   fprintf('----------------------------------------------------\n')
   fprintf('with parameter p=%d\n',p);
   fprintf('----------------------------------------------------\n')
 
   % (1) equal distribution
   lambda1=10^(-p)+i*(1-10^(-p))/(N-1);
   A1=diag(lambda1);

   % (2) clustered distribution
   lambda2=10^(-p)+0.1*(1-10^(-p))*cos(i*3.141592654/(2*(N-1)));
   A2=diag(lambda2);

   % 1
   fprintf('matrix with equal eigenvalue distribution.\n\n');

   [x1,it,res]=sd(A1, [], [], b, zeros(N,1), tol, 1000*p);
   if(it==1000*p)
      fprintf('sd method failed to converge!\n');
      relres=norm(b-A1*x1)/norm(b);
   else
      relres=norm(b-A1*x1)/norm(b);
      fprintf('sd converged to the desired tolerance in %d steps!\n',it);
   end

   % 2
   [x2,flag1,rr1,iter1,rv1] = pcg(A1, b, tol, 1000);
   if (flag1==0)
      fprintf('pcg converged to the desired tolerance in %d steps!\n\n',iter1);
   else
      fprintf('pcg method failed to converge!\n\n');
   end

   % 3
   fprintf('matrix with clustered eigenvalue distribution.\n\n');

   [x3,it]=sd(A2, [], [], b, zeros(N,1), tol, 1000*p);
   if(it==1000*p)
      fprintf('sd method failed to converge!\n');
   else
      fprintf('sd converged to the desired tolerance in %d steps!\n',it);
   end

   %4
   [x4,flag1,rr1,iter1,rv1] = pcg(A2, b, tol, 1000);
   if (flag1==0)
      fprintf('pcg converged to the desired tolerance in %d steps!\n',iter1);
   else
      fprintf('pcg method failed to converge!\n');
   end
end
fprintf('----------------------------------------------------\n')
sd.m0000644000175100017510000000211313203303057011161 0ustar  arbenzarbenzfunction [x,k,reshist] = sd(A,M,M2,b,x,tol,maxit)
%
%STATIT Steepest descent for solving A x = b
%
%       [x,k] = sd(A,M1,M2,b,x,tol)
%	Input:	A  system matrix
%		M1,M2  M = M1*M2 preconditioner
%                      (M1/M2 = [] indicates M1/M2=identity)
%		b  right hand side
%		x  initial vector x^{0} (default x = 0)
%		tol (default tol = eps)
%	Output: x approximate solution
%	        k number of iteration until convergence
%	convergence criterion
%	norm(b - A*x) <= tol*norm(b - A*x0) 

%       12.11.2010/14.11.2017 P. Arbenz

if (nargin < 6), tol = eps; end
if (nargin < 5), x = zeros(size(A,1),1); end

r = b - A*x;
rnrm0 = norm(r);  rnrm = rnrm0;  reshist = rnrm;
if isempty(M),
   z = r;
else
   if isempty(M2),
      z = M\r;
   else
      z = M2\(M\r);
   end
end
q = A*z;
for k=1:maxit
   alpha = (z'*r)/(z'*q);
   x = x + alpha*z;
   r = r - alpha*q;

   rnrm = norm(r);  reshist = [reshist;rnrm];
   if rnrm < tol*rnrm0, return, end

   if isempty(M),
      z = r;
   else
      if isempty(M2),
         z = M\r;
      else
         z = M2\(M\r);
      end
   end
   q = A*z;
end
mypcg.m0000644000175100017510000000223113202552114011672 0ustar  arbenzarbenzfunction [x,k,reshist] = mypcg(A,M,M2,b,x,tol)
%
%MYPCG Conjugate Gradient Method for solving A x = b
%
%       [x,k] = mypcg(A,M1,M2,b,x,tol)
%	Input:	A  system matrix
%		M1,M2  M = M1*M2 preconditioner
%                      (M1/M2 = [] indicates M1/M2=identity)
%		b  right hand side
%		x  initial vector x^{0} (default x = 0)
%		tol (default tol = eps)
%	Output: x approximate solution
%	        k number of iteration until convergence
%	convergence criterion
%	norm(b - A*x) <= tol*norm(b - A*x0) 

%       14.11.2017 P. Arbenz

if (nargin < 6), tol = eps; end
if (nargin < 5), x = zeros(size(A,1),1); end

r = b - A*x;
rnrm0 = norm(r);  rnrm = rnrm0;  reshist = rnrm;
if isempty(M),
   z = r;
else
   if isempty(M2),
      z = M\r;
   else
      z = M2\(M\r);
   end
end
rho = z'*r;
p = z;
for k=1:1000
   q = A*p;
   alpha = (z'*r)/(z'*q);
   x = x + alpha*z;
   r = r - alpha*q;

   rnrm = norm(r);  reshist = [reshist;rnrm];
   if rnrm < tol*rnrm0, return, end
   
   if isempty(M),
      z = r;
   else
      if isempty(M2),
         z = M\r;
      else
         z = M2\(M\r);
      end
   end
   rho0 = rho;
   rho = z'*r;
   beta = rho/rho0;
   p = z + beta*p;
end