MATH4063
MATH4063
MATH4063
2 MATH4063
In this coursework you will build a 2D finite volume solver for the following PDE boundary value problem
????? + ? ? (b??) = ?? (?? , ??) ∈ ?? , (1)
?? = ?? , (?? , ??) ∈ ???? , (2)
where ?? ∶ ?? → ?, ?? ∶ ???? → ? and b ∶ ?? → ?2
.
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
3 MATH4063
2. Suppose we wish to find x ∈ ???
such that
?? x = b, (3)
where ?? is an ?? × ?? matrix and b ∈ ???
.
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 fin