Matlab implementation of the algorithms presented in the lecture notes of Numerical methods for large algebraic systems

• Chapter 2 Rounding errors and LU decomposition
• Chapter 3 Special linear systems
• Chapter 4 Least squares problems
• Chapter 5 Basic iterative methods
• Chapter 6 Krylov methods for symmetric matrices
• Chapter 7 Krylov methods for general matrices

2. Rounding errors and LU decomposition

 

poisson
k = 1;
x=a(:,1);
alpha = gauss(x,k);
aprod = gaussxma ( k, alpha, a );

 

1. Poisson Program to construct a test matrix
2. Algoritme 2.2.1 Computes a Gauss vector
3. Algoritme 2.2.2 Computes the product of a Gauss vector with a matrix A
4. Algoritme 2.3.1 Computes an LU decomposition
5. Algoritme 2.4.1 Computes a solution of a lower triangular system
6. Algoritme 2.4.2 Computes a solution of an upper triangular system
7. Condition estimation

 
% This program fills a matrix with the discretized Poisson equation
%
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 28-11-1996

clear a f
disp(' ')
disp(' ')
disp('This program is used to fill the matrix and righthandside')
disp(' ')
disp('This is a discretization of the following equation:')
disp(' ')
disp(' -c  - c   +u1*c  + u2*c  = f  x,y in [0,1]x[0,1],')
disp('   xx   yy      x       y')
disp(' ')
disp(' with Dirichlet boundary conditions everywhere.')
disp(' ')
disp(' ')
n1 = 3;
n2 = 3;
u1 = 0;
u2 = 0;
upwind = 0;
if norm([u1 u2]) > 0
   upwind = input(' Do you want upwind (choose 1) or not (choose 0) :');
  disp(' ')
   disp(' ')
end
m = n1*n2;
h1 = 1/(n1+1);
h2 = 1/(n2+1);
if upwind ==0
for i = 1:m
   a(i,i) = 2/h1^2+2/h2^2;
end
for i = 1:m-n1
   a(i,i+n1) = -1/h2^2+u2/(2*h2);
   a(i+n1,i) = -1/h2^2-u2/(2*h2);
end
for i = 1:m-1
   a(i,i+1) = -1/h1^2+u1/(2*h1);
   a(i+1,i) = -1/h1^2-u1/(2*h1);
end
end
%
% end while loop
%
if upwind ==1
for i = 1:m
   a(i,i) = 2/h1^2+2/h2^2+abs(u1)/h1+abs(u2)/h2;
end
for i = 1:m-n1
 if u2 > 0
   a(i,i+n1) = -1/h2^2;
   a(i+n1,i) = -1/h2^2-u2/(h2);
 else
   a(i,i+n1) = -1/h2^2+u2/(h2);
   a(i+n1,i) = -1/h2^2;
 end
end
for i = 1:m-1
 if u1 > 0
   a(i+1,i) = -1/h1^2-u1/(h1);
   a(i,i+1) = -1/h1^2;
 else
   a(i,i+1) = -1/h1^2+u1/(h1);
   a(i+1,i) = -1/h1^2;
 end
end
end
%
% end while loop
%

for i = 1:n2-1
   a(i*n1,i*n1+1) = 0;
   a(i*n1+1,i*n1) = 0;
end
for i = 1:n1
   for j = 1:n2
      f(i+(j-1)*n1) = 1+sin(pi*i*h1)*sin(pi*j*h2);
   end
end
%a = sparse(a);
 

 
% This function computes a gauss vector when an original vector
% is given and k is given such that the resulting vector should
% be zero from k+1 to the end  Algoritme 2.2.1
%
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 28-11-1996

function alpha = gauss(x,k);

n = max(size(x));

if x(k) ==0
   disp(' The original vector contains a zero element at position k')
else
   for i = k+1: n
      alpha(i,1) = x(i)/x(k);
   end
end
 

 
% This function computes the product of a gauss vector with a matrix A
% Algoritme 2.2.2
%
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 28-11-1996

function aprod = gaussxma ( k, alpha, a );

[n,q] = size(a);

for j = 1: q
   for i = 1 : k
      aprod(i,j) = a(i,j);
   end
   for i = k+1: n
      aprod(i,j) = a(i,j) - alpha(i)*a(k,j);
   end
end
 

 
% This function computes an LU decomposition.
% Algoritme 2.3.1
%
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 28-11-1996

function [l,u] = ludecom(a)

clear l u

n=max(size(a));
u = a;
for k = 1: n-1
   for i = k+1:n
      l(i,k) = u(i,k)/u(k,k);
      for j = k+1:n
         u(i,j) = u(i,j)-l(i,k)*u(k,j);
      end
   end
end
l(n,n) = 0;
l = l+eye(n);
for k = 1: n-1
   for i = k+1:n
      u(i,k) = 0;
   end
end
 

 
% This function computes a solution of a lower triangular system
% of equations.  Algoritme 2.4.1
%
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 25-11-1996

function y = lsol(l,b);

n = max(size(l));
nulpiv = 1;

for i = 1 : n
   if l(i,i) == 0
      disp(' The main diagonal of the lower triangular matrix')
      disp(' contains a zero element')
      nulpiv = 0;
   else
      if nulpiv == 1
         y(i,1) = b(i);
         for j = 1: i-1
            y(i,1) = y(i,1)-l(i,j)*y(j,1);
         end
         y(i,1) = y(i,1)/l(i,i);
      end
   end
end
 

 
% This function computes a solution of an upper triangular system
% of equations.  Algoritme 2.4.2
%
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 25-11-1996

function y = usol(u,b);

n = max(size(u));
nulpiv = 1;

for i = n : -1 : 1
   if u(i,i) == 0
      disp(' The main diagonal of the upper triangular matrix')
      disp(' contains a zero element')
      nulpiv = 0;
   else
      if nulpiv == 1
         y(i,1) = b(i);
         for j = i+1: n
            y(i,1) = y(i,1)-u(i,j)*y(j,1);
         end
         y(i,1) = y(i,1)/u(i,i);
      end
   end
end
 

 
% This function computes an estimate of the condition number of an
% upper triangular matrix
%
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 28-11-1996

function con = ucond(u);

n = max(size(u));
nulpiv = 1;

for i = 1 : n
   p(i,1) = 0;
end

for i = n : -1 : 1
   if u(i,i) == 0
      disp(' The main diagonal of the upper triangular matrix')
      disp(' contains a zero element')
      nulpiv = 0;
   else
      if nulpiv == 1
         if p(i) > 0
            y(i,1) = ( -1-p(i,1))/u(i,i);
         else
            y(i,1) = ( 1-p(i,1))/u(i,i);
         end
         for j = 1: i-1
            p(j,1) = p(j,1)+u(j,i)*y(i,1);
         end
      end
   end
end

con = norm(y,inf)*norm(u,inf);
 

3. Special linear systems

 
% This function computes a Choleski decomposition of a symmetric postive
% definite matrix Algoritme 3.1.1
%
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 28-11-1996

function l = choldecom(a);

n = max(size(a));
nulpiv = 1;

for k = 1 : n
   for p = 1 : k
      l(k,p) = a(k,p);
   end
end

for k = 1 : n
   for p = 1 : k-1
      l(k,k) = l(k,k) - l(k,p)*l(k,p);
   end
   if l(k,k) <= 0
      disp(' The matrix is not postive definite')
      nulpiv = 0;
   else
      if nulpiv == 1
         l(k,k) = sqrt(l(k,k));
         for i = k+1: n
            for p = 1 : k-1
               l(i,k) = l(i,k) - l(i,p)*l(k,p);
            end
            l(i,k) = l(i,k)/l(k,k);
         end
      end
   end
end
 

4. Least squares problems

1. Algoritme 4.1.1 Computes a Householder vector
2. Algoritme 4.1.2 Computes the product of an Householder vector with a matrix A
3. Algoritme 4.1.3 and 4.1.4 Computes an LU decomposition using Householder transformations

 
% This function computes an Householder vector when an original vector
% is given and k, j are given such that the resulting vector is only
% non-zero from k+1 to j.  Algoritme 4.1.1
%
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 02-12-1996
% date      : 16-05-1997 choice of m adapted (thanks to S.M. Roggeband)

function [beta,v] = hous(x,k,j);

m = max(abs(x(k:j)));
v = 0*x;

if m == 0
   disp(' The original vector only contains zero elements')
else
   alpha = 0;
   for i = k: j
      v(i) = x(i)/m;
      alpha = alpha+v(i)*v(i);
   end
   alpha = sqrt(alpha);
   beta  = 1/(alpha*(alpha + abs(v(k))));
   v(k)  = v(k)+sign(v(k))*alpha;
end
 

 
% This function computes the product of an Householder vector with a matrix A
% Algoritme 4.1.2

function aprod = housxma ( k, j, v, beta, a );

[n,q] = size(a);

aprod =a;

for p = 1: q
   s = 0;
   for i = k: j
      s = s + v(i)*a(i,p);
   end
   s = beta*s;
   for i = k: j
      aprod(i,p) = a(i,p) - s*v(i);
   end
end
 

 
% This function computes an LU decomposition using hous and housxma
% It is a combination of Algorithm 4.1.3 and 4.1.4
%
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 02-12-1996

function [q,u] = houslu(a)

clear q u

n=max(size(a));
u = a;
for i = 1: n-1
   [beta(i), v(:,i)] = hous(u(:,i),i,n);
   u = housxma( i, n, v(:,i), beta(i), u );
end
q = eye(n);
for i = n-1: -1 :1
   q = housxma( i, n, v(:,i), beta(i), q );
end
 

5. Basic iterative methods

1. Gauss-Jacobi method
2. Gauss-Seidel method
3. Successive Over-Relaxation method
4. Chebyshev method

 

% This function is an implementation of the Gauss-Jacobi iterative
% method
%
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 10-12-1996

function x=jacobi(a,f,reps);
clear x n1 n2 n m r
format long e
r = f;
disp(' ')
disp(' ')
disp('     iteration                 norm of the residual')
disp('     ---------                 --------------------')
disp(' ')
disp(' ')
bnorm = norm(f);
[n1,n2] = size(a);
for i = 1 : n1
   m(i,i) = a(i,i);
   x(i,1) = 0;
end
n = m-a;
i = 0;

% Below the main iteration takes place

while (norm(r)/bnorm) > reps & i < 1000
  i = i+1;
  x = inv(m)*(n*x+f);
  r = a*x-f;
  res(i) = norm(r);
  disp([i,norm(r)])
end

% end of main iteration

res = log10(res);
plot(res);
title(' The iterations using Gauss-Jacobi')
xlabel(' iterations')
ylabel(' 10log(residu)')
 

 
% This function is an implementation of the Gauss-Seidel iterative
% method
%
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 10-12-1996

function x=seidel(a,f,reps);
clear x n1 n2 n m r
format long e
r = f;
disp(' ')
disp(' ')
disp('     iteration                 norm of the residual')
disp('     ---------                 --------------------')
disp(' ')
disp(' ')
bnorm = norm(f);
[n1,n2] = size(a);
for i = 1 : n1
   m(i,i) = a(i,i);
   x(i,1) = 0;
   for j = 1:(i-1)
      m(i,j) = a(i,j);
   end
end
n = m-a;
i = 0;

% Below the main iteration takes place

while (norm(r)/bnorm) > reps & i < 1000
  i = i+1;
  x = inv(m)*(n*x+f);
  r = a*x-f;
  res(i) = norm(r);
  disp([i,norm(r)])
end

% end of main iteration

res = log10(res);
plot(res);
title(' The iterations using Gauss-Seidel')
xlabel(' iterations')
ylabel(' 10log(residu)')
 

 
% This function is used to implement SOR iterative method
% When the function is called with omega = 0 an optimal omega is
% estimated. The formulas used to estimate omega are only applicable
% to a restricted class of problems (for instance the Poisson equation)
%
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 17-12-1996

function x=sor(a,f,reps,omega);
clear x n1 n2 n m r
format long e
r = f;
disp(' ')
disp(' ')
disp('     iteration                 norm of the residual')
disp('     ---------                 --------------------')
disp(' ')
disp(' ')

bnorm = norm(f);
[n1,n2] = size(a);
omegal = omega;

% estimate of optimal omega

if omega ==0
   for i = 1 : n1
      m(i,i) = a(i,i);
      x(i,1) = 0;
      for j = 1:(i-1)
         m(i,j) = a(i,j);
      end
   end
   n = m-a;
   mu = max(abs(eig(inv(m)*n)));
   omegal = 2/(1+sqrt(1-mu));
end

for i = 1 : n1
   m(i,i) = a(i,i)/omegal;
   x(i,1) = 0;
   for j = 1:(i-1)
      m(i,j) = a(i,j);
   end
end

n = m-a;

i = 0;
while (norm(r)/bnorm) > reps & i < 1000
  i = i+1;
  x = inv(m)*(n*x+f);
  r = a*x-f;
  res(i) = norm(r);
  disp([i,norm(r)])
end
if omega ==0
   disp(' ')
   disp(' ')
   disp(' The value of the optimal omega is:')
   disp(omegal)
   disp(' ')
   disp(' ')
end
res = log10(res);
plot(res);
title(' The iterations using SOR')
xlabel(' iterations')
ylabel(' 10log(residu)')
 

 
% This function is used to implement Chebychev iterative method
% The exact eigenvalues are used. This method can only be used for
% a symmetric and positive definite matrix.
%
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 17-12-1996

function x1=cheby(a,f,reps);
clear x0 x1 n1 n2 n m r u1 un
format long e
r = f;
disp(' ')
disp(' ')
disp('     iteration                 norm of the residual')
disp('     ---------                 --------------------')
disp(' ')
disp(' ')
bnorm = norm(f);
i = 0;
x0 = 0*f;
r = f;
s = f;
u1=min(eig(a));
un=max(eig(a));
i = 1;
x1 = x0 + 2/(u1+un)*r;
c0 = 1;
c1 = -(u1+un)/(un-u1);
r=f-a*x1;
res(i) = norm(r);
disp([i,norm(r)])

while (norm(r)/bnorm) > reps & i < 1000
  i = i+1;
  chulp = c1;
  c1=-2*(u1+un)/(un-u1)*c1-c0;
  c0 = chulp;
  xhulp = x1;
  x1 = x0-2*c0/c1*(u1+un)/(un-u1)*(2/(u1+un)*r+x1-x0);
  x0 = xhulp;
  r = f-a*x1;
  res(i) = norm(r);
  disp([i,norm(r)])
end

for j = 1 : (i-1)
   fac(j) = res(j+1)/res(j);
end
res = log10(res);
subplot(2,1,1);
plot(res);
title(' The iterations using the Chebychev iterative method')
xlabel(' iterations')
ylabel(' 10log(residu)')  
subplot(2,1,2);
plot(fac)
title(' The reduction factor |r_(i+1)|/|r_i)|')
xlabel(' iterations')
ylabel(' factor')
 

6. Krylov methods for symmetric matrices

 
% This function is used to implement the CG iterative method
%
%
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 17-12-1996

function x=cg(a,f,reps);
clear x n1 n2 n m r
format long e
r = f;
disp(' ')
disp(' ')
disp('     iteration                 norm of the residual')
disp('     ---------                 --------------------')
disp(' ')
disp(' ')
bnorm = norm(f);
i = 0;
x = 0*f;
r = f;
s = f;

while (norm(r)/bnorm) > reps & i < 1000
  i = i+1;
  rnormold = r'*r;
  alfa = rnormold/(s'*(a*s));
  x = x+alfa*s;
  r = r-alfa*a*s;
  beta = r'*r/rnormold;
  s = r+beta*s;
  res(i) = norm(r);
  disp([i,norm(r)])
end

for j = 1 : (i-1)
   fac(j) = res(j+1)/res(j);
end
res = log10(res);
subplot(2,1,1);
plot(res);
title(' The iterations using Conjugate Gradients')
xlabel(' iterations')
ylabel(' 10log(residu)')  
subplot(2,1,2);
plot(fac)
title(' The reduction factor |r_(i+1)|/|r_i)|')
xlabel(' iterations')
ylabel(' factor')
 

7. Krylov methods for general matrices

1. Generalized Conjugate Residual (GCR)
2. Generalized Minimal Residual (GMRES)

 
% This function implements the GCR method
%
% call      : x=gcr(a,f,reps)
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 28-06-2004

function x=gcr(a,f,reps);

format long e
r = f;
i = 0;
x = 0*f;
bnorm = norm(f);

disp(' ')
disp(' ')
disp('     iteration                 norm of the residual')
disp('     ---------                 --------------------')
disp(' ')
disp(' ')

while (norm(r)/bnorm) > reps & i < 1000
   i = i+1;
   s(:,i) = r;
   v(:,i) = a*s(:,i);
   for j = 1:i-1
      alpha = v(:,j)'*v(:,i);
      s(:,i) = s(:,i) - alpha*s(:,j);
      v(:,i) = v(:,i) - alpha*v(:,j);
   end
   beta = norm(v(:,i));
   s(:,i) = s(:,i)/beta;
   v(:,i) = v(:,i)/beta;

   gamma = v(:,i)'*r;
   x = x + gamma*s(:,i);  
   r = r - gamma*v(:,i);  
   disp([i,norm(r)])
end

disp(' ')
disp(' ')
disp('  The exact norm of the residual')
disp(norm(f-a*x))
 

 
% This function implements the GMRES method
%
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 20-12-1996

function x=gmres(a,f,reps);
clear v hh cosgm singm rs x 

format long e
r = f;
i = 0;
v(:,1) = r/norm(r);
rs(1) = norm(r);
rnorm = rs(1);
disp(' ')
disp(' ')
disp('     iteration                 norm of the residual')
disp('     ---------                 --------------------')
disp(' ')
disp(' ')
bnorm = norm(f);

while (rnorm/bnorm) > reps & i < 1000
   i = i+1;
   v(:,i+1) = a*v(:,i);
   for j = 1:i
      hh(j,i) = v(:,j)'*v(:,i+1);
      v(:,i+1) = v(:,i+1) - hh(j,i)*v(:,j);
   end
   hh(i+1,i) = norm(v(:,i+1));
   v(:,i+1) = v(:,i+1)/hh(i+1,i);

%  --- Perform plane rotations on new column.

   if  i > 1 
      for k = 2: i
         thulp = hh(k-1,i);
         hh(k-1,i) =  cosgm(k-1)*thulp + singm(k-1)*hh(k,i);
         hh(k,i)   = -singm(k-1)*thulp + cosgm(k-1)*hh(k,i);
      end
   end

%  --- Calculate new plane rotation.

   thulp = sqrt(hh(i,i)^2 + hh(i+1,i)^2);
   cosgm(i)  = hh(i,i)/thulp;
   singm(i)  = hh(i+1,i)/thulp;
   rs(i+1)   = -singm(i)*rs(i);
   rs(i)     = cosgm(i)*rs(i);
   hh(i,i) = cosgm(i)*hh(i,i) + singm(i)*hh(i+1,i);
   res(i) = abs(rs(i+1));
   rnorm = res(i);
   disp([i,res(i)])
end

%  --- Compute new solution. First solve upper triangular system.

rs(i) = rs(i)/hh(i,i);
for k = i-1: -1: 1
   thulp = rs(k);
   for j = k+1: i
      thulp = thulp-hh(k,j)*rs(j);
   end
   rs(k) = thulp/hh(k,k);
end
x = 0*f;
for j = 1:i
   x = x + rs(j)*v(:,j);
end

disp(' ')
disp(' ')
disp('  The exact norm of the residual')
disp(norm(f-a*x))
 

Contact information:

Back to Numerical methods for large algebraic systems or home page or educational page of Kees Vuik

Last modified on 27-01-2003 by Kees Vuik