Numerical Linear Algebra
Mathwrist Code Example Series
(January 20, 2023)
Abstract
This document gives some C++ code examples on how to use the numerical linear algebra features in Mathwrist’s C++ Numerical Programming Library (NPL).
1 Matrix, View and Iterator
1
// create a 5 x 5 matrix of zero elements
2
Matrix m(5,5);
4
// Obtain a 3 x 3 leading sub matrix view
5
Matrix::View v = m.sub(0,0,3,3);
7
// Change the value from the view also
8
// changes the content of the original matrix.
9
v(1,1) = 0.5;
11
// Now the (1,1) element is 0.5.
12
assert(std::abs(m(1, 1) - 0.5) < 1.e-11);
14
// Let’s have a constant reference to matrix m.
15
const Matrix& mr = m;
17
// From a const matrix, we can only get a constant view.
18
Matrix::ViewConst vc = mr.row(0);
20
// This triggers compilation error if uncommented.
21
// vc(1,1) = 0.5;
23
// Elements of m are contiguous in memory.
24
// Hence an iterator support is available.
25
Matrix::iterator_support its = m.it_support();
27
// its should be not null.
28
assert(its);
30
// Change the element value of the first row.
31
Matrix::iterator_traits::row_iterator
32
it = its->row_begin(0);
33
Matrix::iterator_traits::row_iterator
34
it_end = its->row_end(0);
36
for (; it != it_end; ++it)
37
(*it) = 0.5;
39
// constant matrix has read-only iterator support.
40
Matrix::const_iterator_support itcs = mr.it_support();
41
assert(itcs);
43
Matrix::const_iterator_traits::row_iterator
44
itc = itcs->col_begin(0);
46
// uncomment this triggers a compilation error
47
// *itc = 1;
2 Arithmetic, Inverse and Factorization
1
Matrix A(2, 2);
2
A(0, 0) = 1;
3
A(0, 1) = 0.5;
4
A(1, 0) = 0.6;
5
A(1, 1) = 1.5;
7
Matrix I;
8
Matrix::identity(2, I);
10
// Matrix muplication expression
11
Matrix B = A * I;
12
assert(A.equals_to(B));
14
// Matrix member function.
15
A.multiply(I, B);
16
assert(A.equals_to(B));
18
A.inv(B, Matrix::GENERAL);
19
assert(I.equals_to(A * B));
21
Matrix U, s, V;
22
A.svd(U, s, V);
24
Matrix S(2, 2);
25
S.diagonalize(s);
27
// B = U^T * S * V^T recovers A.
28
B = U * S * *V;
29
assert(A.equals_to(B));
3 Linear System
1
const double _A[] =
2
{
3
1.2205, 1.1124, -3.8936,
4
1.1124, 3.5822, -3.3333,
5
-3.8936, -3.3333, 19.6172
6
};
8
const double _b[] =
9
{
10
2.7694, -1.3499, 3.0349
11
};
13
// Solution of linear system A x = b
14
const double _x[] =
15
{
16
9.1936, -1.6509, 1.6989
17
};
19
Matrix::ViewConst A(3, 3, _A);
20
Matrix::ViewConst b(3, 1, _b);
21
Matrix::ViewConst x_benchmark(3, 1, _x);
23
Matrix x;
25
// solve the linear system: A x = b
26
LinearSystem::solve(A, b, x);
27
assert(x.equals_to(x_benchmark, 1.e-4));
29
// solve by conjugate gradient step update
30
x.zeros();
32
// p = -(A * x - b) at the first call.
33
// Initially x = 0, so p = b;
34
Matrix p = b;
36
// Initially x = 0, residual r = -b.
37
Matrix r = -b;
38
double r2_prev = 1.e10;
40
// The first step update should reduce residuals.
41
for (int i = 0; i < 3; ++i)
42
{
43
LinearSystem::cg_step_update(A, x, p, r);
44
double r2_next = r.vector_norm_2();
45
assert(r2_next < r2_prev);
46
r2_prev = r2_next;
47
}
49
// 3 step update should solve the equation exactly.
50
assert(r2_prev < 1.e-12);
51
assert(x.equals_to(x_benchmark, 1.e-4));