Description
1. Implementation of the Gauss-Seidel method. (10 marks)
Implement the Gauss-Seidel algorithm of Eq. (4.8) in the course notes in MATLAB using
only primitive commands (i.e. do not appeal to the built-in linear algebra functionality in
MATLAB). You should use the function headers given below and submit the code in separate
m-files to the LEARN drop box on the course webpage. You can implement the algorithm for
a general dense matrix A (not assuming A is sparse), and you can directly use the formulation
as in Eq. (4.8) (there is no need to implement this using the interpretation of a stand-alone
preconditioner, as in Eq. (4.6)).
function u = GaussSeidel(A, u0, f, tol, maxit)
% Usage: u = GSSequence(A, u0, f, tol, maxit)
1
% Perform a sequence of Gauss-Seidel iterations on
% A u = f using the initial estimate u0, until the 2-norm
% of the residual has been reduced by a factor
% of at least tol, or until the number of iterations
% reaches maxit.
You can download VerifyGaussSeidel.m from LEARN to check your implementation. Please
report on the output for input n = 100. (It should be a small number, with magnitude around
1e-14 or so.)
Also, include a textual copy of your code in the assignment package submitted on Crowdmark.
Note: to obtain good performance for large problems (see next question), it may be best to
vectorize the implementation, which allows you to avoid the inner loop of the double loop in
the Gauss-Seidel algorithm.
2. Gauss-Seidel applied to a special dense matrix. (10 marks)
Consider the n × n linear system
n 1 1 · · · 1
1 n 1 · · · 1
1 1 n · · · 1
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1 1 1 · · · n
x1
x2
x3
.
.
.
x4
=
1
2
3
.
.
.
n
. (1)
Write a Matlab script that performs the following computations, and provide the output and
answer the questions in your Crowdmark submission:
(a) Compute the solution to this system for n = 250, 500, 1000, 2000, 4000, 8000, 16000 using
LU decomposition (you can use Matlab’s linsolve). (You can use smaller n if your
computer runs out of memory.) Determine and tabulate the amount of compute time
required for each value n. (Note that Matlab’s linsolve uses a highly optimized implementation in a high-level programming language that is much faster than execution of
code written in the Matlab language itself.)
(b) Compute the solution to this system using your implementation of Gauss-Seidel from
the previous question for n = 250, 500, 1000, 2000, 4000 with x0 = [0 0 0 · · · 0 0]T and
tolrel = 10−8
. Determine the amount of compute time required for this solution. How
many iterations do you need to achieve convergence as a function of n? (Note that what
you find here is very specific for this special type of matrix.)
(c) Compare the results from parts (a) and (b). How do you expect the execution time of
each algorithm to scale asymptotically as a function of n? To which extent does this
agree with what you observe numerically (compute and tabulate ratios of consecutive
execution times as n doubles)? Which algorithm/implementation would be best for
small or medium problem sizes? And for very large problem sizes?
Submit your code to LEARN, and include a textual copy of your code in the assignment
package submitted on Crowdmark.
3. Implementation of the steepest descent and conjugate gradient methods. (15 marks)
You are asked to implement the steepest descent and conjugate gradient methods for A ~x = ~b
in Matlab, and apply the methods to the 2D model problem (you can download
build laplace 2D kron.m from LEARN to build the 2D Laplacian matrix).
2
(a) Implement the following in m-file steepest.m:
function [x res steps]=steepest(A,x0,b,tol,maxit)
% Performs a sequence of steepest descent iterations on
% A x = b using the initial estimate x0, until the 2-norm
% of the residual is smaller than tol (relative
% to the initial residual), or until the number of iterations
% reaches maxit. ‘steps’ is the number of steps
% that were performed, and ‘res’ is a vector with the
% residual size after every interation (the 2-norm of
% the residual vector).
(b) Implement the following in m-file conjGrad.m:
function [x res steps]=conjGrad(A,x0,b,tol,maxit)
% Performs a sequence of conjugate gradient iterations
% on A x = b using the initial estimate u0, until the 2-norm
% of the residual is smaller than tol (relative
% to the initial residual), or until the number of iterations
% reaches maxit. ‘steps’ is the number of steps
% that were performed, and ‘res’ is a vector with the
% residual size after every interation (the 2-norm of
% the residual vector).
Submit your m-files steepest.m and conjGrad.m to the LEARN drop box, and provide
a printout in your assignment package to be submitted in Crowdmark.
(c) Download test steepest cg.m to test the correctness of (a) and (b). Report on the
error and number of steps for each method. (The errors should be smaller than 10−10.)
(d) Make a driver program CG driver.m that applies steepest.m and conjGrad.m to A~x = ~b
with A being the 2D model problem matrix generated by build laplace 2D kron.m,
and ~b a vector with all twos. Use maxit=500 and tol=1e-8. Generate a plot in which you
compare, for N = 32, the residual convergence for the steepest descent and conjugate
gradient methods (starting from a zero initial guess), as a function of the iteration. For
each of the methods, plot the 10-log of the residuals as a function of iteration number.
Which method requires the least iterations to obtain an accurate solution, and is this as
you expected?
(e) Provide a table in which you list, for steepest descent and conjugate gradient, how many
iterations you need to reduce the residual by 1e − 6, for N = 16, 32, 64. (You can
use maxit= 20, 000.) What are the total number of unknowns and the approximate
condition number for each problem size? (You can use κ ≈ 4n/π2
, where n = N2
,
which can be derived from a formula for the eigenvalues of A.) For each method, do
you see the expected behaviour in the number of required iterations as a function of the
total number of unknowns? (Explain.) Using these numerical results, briefly discuss the
computational cost/complexity of each of the methods as a function of the total number
of unknowns (discuss the cost per iteration, the number of iterations required, and the
total cost, as a function of total problem size). Which method is best for large problem
size? (Use what you learned in Example 4.14 of the course notes.)
4. Implementation of the preconditioned CG and GMRES methods. (15 marks)
You are asked to implement the preconditioned CG and GMRES methods for A ~x = ~b in Mat3
lab, and apply the methods to the 2D model problem (download build laplace 2D kron.m
from LEARN).
(a) Implement the following in m-file myGMRES.m:
function [x res steps]=myGMRES(A,x0,b,tol,maxit)
% Performs a sequence of GMRES iterations on
% A x = b using the initial estimate x0, until the 2-norm
% of the residual is smaller than tol (relative
% to the initial residual), or until the number of iterations
% reaches maxit. ‘steps’ is the number of steps
% that were performed, and ‘res’ is a vector with the
% residual size after every interation (the 2-norm of
% the residual vector).
(b) Implement the following in m-file myGMRES SSOR.m:
function [x res steps]=myGMRES_SSOR(A,x0,b,tol,maxit)
% Performs a sequence of right-preconditioned GMRES iterations
% on A x = b with the initial estimate x0, using the SSOR preconditioner,
% until the 2-norm of the residual is smaller than tol (relative
% to the initial residual), or until the number of iterations
% reaches maxit. ‘steps’ is the number of steps
% that were performed, and ‘res’ is a vector with the
% residual size after every interation (the 2-norm of
% the residual vector).
(c) Implement the following in m-file myCG SSOR.m:
function [x res steps]=myCG_SSOR(A,x0,b,tol,maxit)
% Performs a sequence of preconditioned CG iterations
% on A x = b with the initial estimate x0, using the SSOR preconditioner,
% until the 2-norm of the residual is smaller than tol (relative
% to the initial residual), or until the number of iterations
% reaches maxit. ‘steps’ is the number of steps
% that were performed, and ‘res’ is a vector with the
% residual size after every interation (the 2-norm of
% the residual vector).
Implementation notes:
– for GMRES, you can solve the small least-squares problem min kM ~y − ~fk by using
Matlab’s backslash operator as in ~y = M\
~f (which solves over-determined systems
in the least-squares sense using QR decomposition), or you can use the normal
equations
– to build the preconditioner, you can use tril and triu to extract the strictly lower
and upper triangular parts of A; you can use spdiags to extract the diagonal of A,
and again to build a sparse diagonal matrix containing the diagonal of A
– take ω = 1.9 in SSOR
– as was discussed in the course notes, since for SSOR the matrix P contains inverses
of sparse triangular matrices, it is a bad idea to form P explicitly, because inverting
4
the sparse triangular matrices will result in dense matrices, and also, inverting a
matrix is in general substantially more expensive than solving a linear system with
that matrix; so rather than forming the P matrices directly and computing the
preconditioned residual P ~rk as a matrix-vector product, you should compute the
preconditioned residual ~qk in
~qk = P ~rk
by solving the linear system
P
−1
~qk = ~rk
with P
−1 a product of sparse matrices; to solve this system for the SSOR preconditioner, you need to do two consecutive triangular solves, and you can again
use backslash in Matlab for each of these solves (which recognizes sparse and triangular matrices, and solves the system uses forward or backward substitution, as
appropriate), or use forward and backward substitution solvers
Submit your m-files to the LEARN drop box, and provide a printout in your assignment
package to be submitted in Crowdmark.
(d) Download test iterative.m to test the correctness of (a–c). Report on the errors and
number of steps. (The errors should be close to 1e-12.)
(e) Make a driver program driverPCG PGMRES.m that applies myGMRES.m, myGMRES SSOR.m,
and myCG SSOR.m to A~x = ~b with A being the 2D model problem matrix generated by
build laplace 2D kron.m, and ~b a vector with all ones. Use maxit=400 and tol=1e-10.
Generate a plot in which you compare, for N = 32, the residual convergence for the
three methods (starting from a zero initial guess), as a function of the iteration. On this
plot, for each of the methods, plot the 10-log of the residuals as a function of iteration
number. Also compare with CG, using your file conjGrad.m from Question 4 (or, a
simplification of myCG SSOR.m, without preconditioner).
(f) Extending driverPCG PGMRES.m, for problem sizes N equal to [80 90 100 110 120 130
140], make a loglog plot of the number of iterations required for convergence as a
function of total problem size n = N2
, for each of the 4 methods (CG, CG+SSOR,
GMRES, GMRES+SSOR), all on the same figure. Use the polyfit command to fit
a polynomial of degree 1 to the measured (n, iterations) points of problem size versus
iterations, for each of the methods. The fitted slope for CG and GMRES should be close
to 1/2, because the number of iterations is O(
√
n) for these methods. It can be shown
that SSOR is an optimal preconditioner for CG for the 2D model problem in terms of
asymptotic order, with number of iterations O(n
1/4
).
Note: This means that the total work is O(n
5/4
) when using the SSOR preconditioner,
and O(n
3/2
) without it (both for CG and GMRES). GMRES gives very similar results for
the 2D model problem than CG, in terms of the number of required iterations. However,
GMRES is much more versatile since it can also be applied to non-SPD systems.
5. Question on neural network implementation will be added later this week as the topic is being
covered in class. (xx marks)
6. Question on neural network implementation will be added later this week as the topic is being
covered in class.. (xx marks)
5