代写 R C++ algorithm math scala compiler OOP in C++

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 . How you implement the vector is up to you. Define necessary and convenient member functions, and in particular appropriate constructor functions, a destructor and a function that returns the dimension of the vector. Furthermore, make sure that the following operations are properly defined for your class: x = y, z = x + y, z = x – y, x = – y, x = r * yandx += y,x *= r;wherex,y,zareinstancesofmathvector and r is of type T. Finally, also overload the operator* function to return the scalar/dot product between two vectors, i.e. a . b = 􏰁Ni=1 ai bi, where z denotes the complex conjugate of z ∈ C.
2∗. Matrices
Create the abstract (!) base class template mathmatrix that can be used to represent square matrices in KN×N. Create a pure virtual member function to overload the operator* function to return the result of a matrix times mathvector multiplication. Provide also a pure virtual member function y_eq_Ax(mathvector &y, const mathvector &x) that computes and returns Ax in-situ on y, i.e. without creating any temporary objects.
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::power_method(const mathvector &, int), which can be used to find the largest eigenvalue of e.g. the matrices AN ∈ RN×N and B ∈ R4×4, where
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 x_a;
for (int i = 0; i < N; i++) x_a.push_back(i+1); CW5_matrix A(N);
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 vb(b, b+4); mathvector x_b(vb); // x = (1,1,1,1)^T
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 row1(r1, r1+4), row2(r2, r2+4), row3(r3, r3+4), row4(r4,r4+4);
vector > BB;
BB.push_back(row1); BB.push_back(row2); BB.push_back(row3); BB.push_back(row4);
fullmatrix B(BB);
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 D(row4);
lambda = D.power_method(x_b, K);
cout << "largest lambda (in modulus) for D = " << lambda << endl; vector > > CC; // … fill CC appropriately …
fullmatrix > C(CC); // … etc …
}
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::CG_solver() that takes the two mathvectors x0 and b, the parameters double tolerance, int kmax and mathmatrix P. Let the function perform the above algorithm and return the solution in place of x0. The return value of CG_solver() should be the number of iterations k.
[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.