MATH4063 - SCIENTIFIC COMPUTING AND C++ Coursework 1: GMRES linear algebra solver
AUTUMN SEMESTER 2022-2023
MATH4063 - SCIENTIFIC COMPUTING AND C++
Coursework 1 - Released 30th October 2023, 4pm
Your work should be submitted electronically via the MATH4063 Moodle page by 12noon on Monday 20th November (unless you have arranged an extension). Since this work is assessed, your submission must be entirely your own work (see the University’s policy on Academic Misconduct). Submissions up to five working days late will be marked, but subject to a penalty of 5% of the maximum mark per working day.
The marks for each question are given by means of a figure enclosed by square brackets, eg [20]. There are a total of 100 marks available for the coursework and it contributes 45% to the module. The marking rubric available on Moodle will be applied to each full question to further break down this mark.
You are free to name the functions you write as you wish, but bear in mind these names should be meaningful. Functions should be grouped together in .cpp files and accessed in other files using correspondingly named .hpp files.
All calculations should be done in double precision.
A single zip file containing your full solution should be submitted on Moodle. This zip file should contain three folders called main, source and include, with the following files in them:
main:
• q1d.cpp • q2c.cpp • q3b.cpp • q4b.cpp
source:
• vector.cpp
• dense_matrix.cpp • csr_matrix.cpp
• linear_algebra.cpp • finite_volume.cpp
include:
• vector.hpp
• dense_matrix.hpp
• csr_matrix.hpp
• linear_algebra.hpp
• finite_difference.hpp
Prior to starting the coursework, please download the CW1_templates.zip from Moodle and extract the files. More information about the contents of the files included in the template is given in the questions below.
Hint: When using a C++ struct with header files, the whole struct needs to be defined fully in the header file, and the header file included in the corresponding .cpp file. Include guards should also be used.
MATH4063 Turn Over
2 MATH4063 In this coursework you will build a 2D finite volume solver for the following PDE boundary value problem
−𝛥𝑢+∇⋅(b𝑢)=𝑓 (𝑥,𝑦)∈𝛺, (1) 𝑢=𝑔, (𝑥,𝑦)∈𝜕𝛺, (2)
where 𝑓 ∶ 𝛺 → R, 𝑔 ∶ 𝜕𝛺 → R and b ∶ 𝛺 → R2.
In order to solve this problem, you will first define a sparse matrix structure, then write functions to apply the GMRES linear algebra solver and finally build and solve the linear system arising from the finite volume approximation of (1)-(2).
1. Matrices arising from the discretisation of partial differential equations using, for example, finite volume methods, are generally sparse in the sense that they have many more zero entries than nonzero ones. We would like to avoid storing the zero entries and only store the nonzero ones.
A commonly employed sparse matrix storage format is the Compressed Sparse Row (CSR) format. Here, the nonzero entries of an 𝑛 × 𝑛 matrix are stored in a vector matrix_entries, the vector column_no gives the column position of the corresponding entries in matrix_entries, while the vector row_start of length 𝑛+1 is the list of indices which indicates where each row starts in matrix_entries. For example, consider the following:
⎛8 0 0 2⎞
𝐴=⎜0 3 1 0⎟ ⟶ ⎜0 0 4 0⎟
⎝6 0 0 7⎠
matrix_entries=(8 2 3 1 4 6 7) column_no=(0 3 1 2 2 0 3) row_start=(0 2 4 5 7)
Note, in the above, C++ indexing has been assumed, i.e, indices begin at 0.
-
(a) In csr_matrix.hpp, define a C++ struct called csr_matrix to store a matrix in CSR format. In addition to matrix_entries, column_no and row_start, you should store the number of rows of the matrix explicitly.
-
(b) In csr_matrix.cpp, write a C++ function that will set up the matrix 𝐴 from above in CSR format. Remember, if you are using dynamically allocated memory, then you should also have corresponding functions that will deallocate the memory you have set up.
-
(c) In csr_matrix.cpp, write a C++ function that takes as input a matrix 𝐴 stored in CSR format and a vector x and computes the product 𝐴x. The prototype for your function should be:
void MultiplyMatrixVector(csr_matrix& matrix,double* vector, double* productVector)
Hence, the input vector and the output productVector should be pointers to dynamically allocated arrays. In particular, it should be assumed that productVector has been preallocated to the correct size already.
-
(d) By setting a vector x = (4,−1,3,6)⊤, write a test program in q1d.cpp to compute and print to the screen the product 𝐴x, where 𝐴 is the matrix given above.
[20 marks]
MATH4063
2. Suppose we wish to find x ∈ R𝑛 such that 𝐴x = b,
where𝐴isan𝑛×𝑛matrixandb ∈R𝑛.
(3)
𝑖 𝑖=1
3: Compute q𝑘+1 = 𝐴q𝑘
3 MATH4063
One algorithm for solving this problem is the (restarted) Generalised Minimal RESidual (GMRES) algorithm. The method is too complicated to explain here, but works to quickly find approximations x𝑘 = x0 + y𝑘 where y𝑘 ∈ 𝒦𝑘 ∶= Span{𝐴q0, 𝐴2q0 ... 𝐴𝑘q0} for 𝑘 = 1, 2, .... y𝑘 is chosen to minimise the residual
‖b − 𝐴x𝑘‖2.
Here x0 is some initial guess vector and q0 is the normed initial residual
q0= b−𝐴x0 . ‖b − 𝐴x0‖2
𝒦𝑘 is called a Krylov subspace of 𝐴.
The algorithm stops when ‖b − 𝐴x𝑘‖2 < tol for some termination tolerance tol. As the method becomes
very memory inefficient when 𝑘 is large, the method is restarted every so often and x𝑘 reset to be x0. An incomplete GMRES algorithm function PerformGMRESRestarted() has been written in
linear_algebra.cpp.
A key component of the GMRES algorithm is the Arnoldi iteration that seeks to find an orthonormal basis of 𝒦𝑘. At the 𝑘th step of the iteration, the Arnoldi method constructs the following matrix decomposition of 𝐴:
𝐴𝑄𝑘 = 𝑄𝑘+1𝐻̃𝑘,
where the columns of 𝑄𝑘 (𝑄𝑘+1) contain the orthonormal basis of 𝒦𝑘 (𝒦𝑘+1, resp.) and 𝐻̃𝑘 is a (𝑘+1)×𝑘 upper Hessenberg matrix. That is, a matrix that is nearly upper triangular but has non-zero components on the first subdiagonal.
The 𝑘th step of the Arnoldi algorithm is:
Algorithm 1 One step of the Arnoldi Iteration. Require: 𝑘 > 0, 𝐴, 𝑄𝑘:
-
1: Let q𝑖 be the 𝑖th column of 𝑄𝑘.
-
2: Leth={h}𝑘+1 beavectoroflength𝑘+1.
-
4: for𝑖=1,...,𝑘do
-
5: h𝑖 = q𝑘+1 ⋅ q𝑖.
-
6: q𝑘+1 = q𝑘+1 − h𝑖q𝑖.
-
7: end for
-
8: h𝑘+1 = ‖q𝑘+1‖2.
-
9: q𝑘+1 = q𝑘+1/h𝑘+1.
-
10: 𝑄𝑘+1 = [𝑄𝑘, q𝑘+1].
-
11: return 𝑄𝑘+1 and h.
(a) In linear_algebra.cpp, write a C++ function which implements one step of the Arnoldi iteration
method defined above.
The function should have the following prototype
void PerformArnoldiIteration(csr_matrix& matrix,
dense_matrix& krylov_matrix , int k,
double*
hessenberg)
Here, matrix is 𝐴, k is the step of the iteration to perform, krylov_matrix is the matrix containing
the orthonormal basis, where each row is a basis vector. Upon entry, krylov_matrix should have 𝑘
rows and upon exit it should contain 𝑘 + 1 rows, with the new basis vector in the last row.
Finally, upon exit, hessenberg should contain h, which is the final column of 𝐻̃𝑘. You may assume that
hessenberg has been preallocated to be of length 𝑘 + 1 before the call to PerformArnoldiIteration.
Your function should make use, where possible, of prewritten functions defined in dense_matrix.cpp
and vector.cpp. Your code should also make use of the matrix multiplication function from Q1.
Once you have written PerformArnoldiIteration() the GMRES function should function as intended. Note: Storage of the basis functions in the rows of krylov_matrix, rather than in the columns,
improves efficiency of the code.
-
(b) In csr_matrix.cpp, write a C++ function that will read from a file a matrix already stored in CSR format and a vector. You may assume the file structures are as in matrix1.dat and vector1.dat on Moodle and you may use these data files to test your function.
-
(c) Write a test program in file q2c.cpp that will read in the matrix 𝐴 from matrix2.dat and the vector x from vector2.dat, compute b = 𝐴x, then use PerformGMRESRestarted() with the default input arguments to find an approximation x̂ to x. At the end of the calculation, print to the screen the error ‖x − x̂‖2.
[20 marks]
3. Thefilemesh.hppcontainsastructthatdefinesameshdatastructuremeshforageneralmeshcomprising axis-aligned rectangular cells. In particular, each cell in the mesh has an additional struct called cell_information that contains, among other things, information about the cell neighbours. Familiarise yourself with these data structures by looking in mesh.hpp.
mesh.cpp contains two functions that will generate meshes, they are:
• ConstructRectangularMesh() - this constructs a mesh on the rectangular domain 𝛺𝑅 = [𝑎, 𝑏] ×
[𝑐, 𝑑].
• ConstructLShapedMesh() - this constructs a mesh on the L-shaped domain 𝛺𝐿 = 𝛺𝑅\𝛺𝐶, where
𝛺𝐶 =[(𝑎+𝑏)/2,𝑏]×[(𝑐+𝑑)/2,𝑑].
-
(a) In finite_volume.cpp, write a C++ function that will create the storage for a matrix 𝐴 in CSR format and a RHS vector F required for a cell-centred finite volume method for solving (1)-(2). You should follow the procedure outlined in the Unit 6 lecture notes. As one of the inputs, your function should take in a variable of type mesh.
-
(b) In csr_matrix.cpp, write a C++ function that will output to the screen a matrix stored in CSR format in the same style as in matrix1.dat.
-
(c) In Q3c.cpp, write a program that will ask the user to supply the number of cells in each coordinate direction of a rectangular mesh, sets up the mesh using ConstructRectangularMesh() then calls the function from part (a) to set up the corresponding matrix and finally prints it to the screen using the function from part (b).
[20 marks]
MATH4063
4.
(a) In finite_volume.cpp, write a function that takes in a mesh, uses the function from Q3(a) to construct 𝐴 and F, then populates it with the correct entries to solve problem (1)-(2) using the cell-centred finite volume method, as outlined in the Unit 6 notes. The function should also take as input the functions 𝑓(𝑥, 𝑦), b(𝑥, 𝑦) and the Dirichlet boundary function 𝑔(𝑥, 𝑦).
(b) In Q4b.cpp, write a main program to ask the user to select from the following problems and supply the number of cells in each coordinate direction.
RectangularMesh-𝑎=0,𝑏=1,𝑐=0and𝑑=1;
1. •
• 𝑓(𝑥,𝑦)=1;
-
𝑔(𝑥,𝑦)=0;
-
b = 0.
L-shapedMesh-𝑎=0,𝑏=1,𝑐=0and𝑑=1;
2. •
-
𝑓(𝑥, 𝑦) = 8𝜋2 cos(2𝜋𝑥) cos(2𝜋𝑦);
-
𝑔(𝑥,𝑦)=cos(2𝜋𝑥)cos(2𝜋𝑦);
-
b = 0.
RectangularMesh-𝑎=−1,𝑏=1,𝑐=−1and𝑑=1;
3. •
• 𝑓(𝑥,𝑦)=1;
-
𝑔(𝑥,𝑦)=0;
-
b=(10,10)⊤.
L-ShapedMesh-𝑎=0,𝑏=1,𝑐=0and𝑑=1;
4. •
• 𝑓(𝑥,𝑦)=0;
•
𝑔(𝑥,𝑦)={1, 𝑥=0,0.25<𝑦<0.75, 0, otherwise;
5 MATH4063
50𝑦 −50𝑥 ⊤ • b=(√𝑥2+𝑦2,√𝑥2+𝑦2) .
The code should then set up the linear system arising from the finite volume discretisation and solve the system
𝐴uh = F
using PerformGMRESRestarted().
Finally, print to the screen the maximum value of uh.
Hint: Once you have computed uh you can output it to together with the mesh to a file using OutputSolution() in finite_volume.cpp. plot_solution.py can then be used to plot the solution in Python.
Note, if you are unable to get the iterative solver from Q2 working, then you may create the finite volume matrix 𝐴 as if it were a dense matrix (i.e store all the zero entries) and use the function PerformGaussianElimination() from dense_matrix.cpp to solve the system of equations. This will incur a small penalty. Note, an illustration of the use of PerformGaussianElimination() can be found in the main program inside gaussian_elimination_test.cpp.
