OOP in C++
Dr Robert Nu ̈rnberg
Exercise 5
Tasks marked with a ∗ are assessed coursework. Hand in your solutions to these via email to rn@ic.ac.uk. (Resit students do not need to submit coursework.) Use the subject line “ C++ CW: surname firstname CW5”, where surname firstname CW5.cpp is the attached file that contains your solution. The course will be assessed based on 5 pieces of coursework (25%) and an end of term driving test (75%). Your submission must be your own work (submissions will be checked for plagiarism), and it should compile (and run) with the GNU C++ compiler g++. The deadline for submitting the coursework is 10pm on 24/03/2019. fan.qu15 fan.qu15
1. Mathematical vectors
Create a template mathvector<> that can be used to represent vectors in KN, where K = R or K = C. For the latter you may #include
2∗. Matrices
Create the abstract (!) base class template mathmatrix
From this base class derive the class fullmatrix, which stores all elements of a matrix explicitly, and the class diagmatrix, which represents a diagonal matrix and hence stores only N elements. Implement the matrix times mathvector function for both of these classes.
Note that you can now also derive any implicitly declared matrix from the base class, i.e. a matrix for which you do not store all N2 elements explicitly. All you have to do, is to declare the action of that matrix on a vector.
3. Power method
The largest eigenvalue (in magnitude) of a matrix A ∈ KN×N can be found using the power method:
1. Choose an initial guess x0 ∈ KN 2.Letq:= x0 ,where|x0|2=x0.x0
|x0 | 3. For k = 1 : K
4. z := A q
5. q:= z
|z| 6. end
7. Returnλ=q.(Aq)
4∗. Power method example
Write a member function mathmatrix
2−1 0… 0−1
3 1 2 5 1137 4×4
−1
AN=
−1 0 … 0 −1 2
Here AN should be defined implicitly. To this end, define a separate class for it, called CW5_matrix. In particular, your code should be able to execute all of the following statements correctly.
[Hint: To test your code for acceptable speed, compile it with the GNU compiler with the command: g++ -Wall -O3 surname_firstname_CW5.cpp.]
2 −1 0 … 0
. . N×N .. .. ∈R
and B= ∈R . 0 0
0… 0−1 2−1 2324 5742
OOP in C++
Dr Robert Nu ̈rnberg
int main() {
int N = 10000, K = 4000;
mathvector
for (int i = 0; i < N; i++) x_a.push_back(i+1);
CW5_matrix
double lambda = A.power_method(x_a, K);
cout << "largest lambda (in modulus) for A = " << lambda << endl;
double b[] = {1.0, 1.0, 1.0, 1.0};
vector
cout << " x . x = " << x_b * x_b << endl;
double r1[] = {3.0, 1.0, 2.0, 5.0}, r2[] = {1.0, 1.0, 3.0, 7.0},
r3[] = {2.0, 3.0, 2.0, 4.0}, r4[] = {5.0, 7.0, 4.0, 2.0};
vector
vector
BB.push_back(row1); BB.push_back(row2); BB.push_back(row3); BB.push_back(row4);
fullmatrix
B.y_eq_Ax(x_a, x_b); cout << x_a << endl;
lambda = B.power_method(x_b, K);
B -= lambda;
double l2 = B.power_method(x_b, K);
B += lambda;
cout << "The spectrum of B lies between " << lambda << " and " << l2 + lambda << endl;
diagmatrix
lambda = D.power_method(x_b, K);
cout << "largest lambda (in modulus) for D = " << lambda << endl;
vector
fullmatrix
}
5. Conjugate gradient solver
For a symmetric and positive definite matrix A ∈ RN×N the system Ax = b can be solved using the conjugate gradient algorithm:
1. Setk:=0andchooseaninitialguessx0 ∈RN Letr0 :=b−Ax0,ρ0 :=|r0|2 =r0.r0
2. While k < kmax and √ρk > tolerance ∗ |b| do
3. z:=Prk, τk:=z.rk
4. If k = 0 then set β := 0 and v := 0, else set β := τk/τk−1.
// x = (1,2,…,N)^T
// matrix A
// matrix B
// y = B * x
// B = B – lambda * Id
// B = B + lambda * Id
5. v := z + β v
6. w := A v
7. γ := τk /v . w
8. xk+1 :=xk +γv,
9. ρk+1 := |rk+1|2
10. k := k + 1
11. Return xk
rk+1 :=rk −γw
Observe that the above algorithm uses only one matrix times vector multiplication per iteration (6.). The classical conjugate gradient method uses P = I for the preconditioner P in 3., where I is the identity matrix. Otherwise P is usually a very simple matrix approximating A−1.
Write a member function mathmatrix
[Hint: Although efficiency is not the priority in this exercise, try to increase the efficiency of your code, by using the method y_eq_Ax() in 6. and making use of the operators += and *= throughout.]
6. Conjugate gradient example
Use your solver from 5. to compute the solution of some linear systems. Include at least the example matrices A10000 and B := 10 I + B from 4. In each case, solve the system A x = b for an arbitrary right hand side b ̸= 0 with i bi = 0. Once the solution x⋆ is obtained, compute its residual r = |A x⋆ − b|. Also compare the number of iterations when using the preconditioner P = [diag(A)]−1 compared to P = I.