From bfa24ea5cba1f73e2a5cc5bb3795b0bb50d9f1ff Mon Sep 17 00:00:00 2001 From: Joseph Lee Date: Fri, 7 Apr 2023 16:57:44 +0100 Subject: [PATCH 1/3] Update eigen lecture & exercises --- AUTHORS | 1 + exercises/eigen/explicit.cpp | 71 +++---- exercises/eigen/implicit.cpp | 56 +++--- exercises/eigen/movie.py | 28 +++ exercises/eigen/sparse.cpp | 58 +++--- lectures/eigen/README.md | 347 ++++++++++++++++++++++++++++++----- lectures/eigen/index.html | 2 +- 7 files changed, 435 insertions(+), 128 deletions(-) create mode 100644 exercises/eigen/movie.py diff --git a/AUTHORS b/AUTHORS index dbb4fd9..2b27819 100644 --- a/AUTHORS +++ b/AUTHORS @@ -1 +1,2 @@ Rupert Nash +Joseph Lee \ No newline at end of file diff --git a/exercises/eigen/explicit.cpp b/exercises/eigen/explicit.cpp index ba27e73..3e0003e 100644 --- a/exercises/eigen/explicit.cpp +++ b/exercises/eigen/explicit.cpp @@ -6,40 +6,43 @@ int main() { - int n = 12; + int n = 20; + int steps = 200; + std::vector Avec(n * n); + + + // Set up matrix A + Eigen::Map A(Avec.data(), n, n); + A = Eigen::MatrixXd::Identity(n, n); + double delta = 0.4; + for (int i = 0; i < n - 1; ++i) + { + A(i + 1, i) += delta; + A(i + 1, i + 1) += -delta; + + A(i, i) += -delta; + A(i, i + 1) += +delta; + } + + std::cout << "A = \n" << A << std::endl + << std::endl; + + + // T_n + Eigen::VectorXd b(n); + b.setZero(); + b.head(n / 2).array() = 1.0; + + std::ofstream f; + f.open("explicit_sim.txt"); + for (int i = 0; i < steps; ++i) + { + f << b.transpose() << std::endl; + // update time-step T_{n+1} = A.T_{n} + b = A * b; + } - std::vector Avec(n * n); - - Eigen::Map A(Avec.data(), n, n); - A = Eigen::MatrixXd::Identity(n, n); - - double delta = 0.4; - - for (int i = 0; i < n - 1; ++i) - { - A(i + 1, i) += delta; - A(i + 1, i + 1) += -delta; - - A(i, i) += -delta; - A(i, i + 1) += +delta; - } - - std::cout << A << std::endl << std::endl; - - Eigen::VectorXd b(n); - b.setZero(); - b.head(n / 2).array() = 1.0; - - std::ofstream f; - f.open("a.txt"); - for (int i = 0; i < 20; ++i) - { std::cout << b.transpose() << std::endl; - f << b << std::endl << std::endl << std::endl; - b = A * b; - } - - std::cout << b.transpose() << std::endl; - return 0; -} + return 0; +} \ No newline at end of file diff --git a/exercises/eigen/implicit.cpp b/exercises/eigen/implicit.cpp index 2ae3a1d..e2b0a39 100644 --- a/exercises/eigen/implicit.cpp +++ b/exercises/eigen/implicit.cpp @@ -5,38 +5,42 @@ int main() { - int n = 12; + int n = 20; + int steps = 100; - Eigen::MatrixXd A(n, n); - A = Eigen::MatrixXd::Identity(n, n); + // Set up matrix A + Eigen::MatrixXd A(n, n); + A = Eigen::MatrixXd::Identity(n, n); - double delta = -0.4; + double delta = -0.4; - for (int i = 0; i < n - 1; ++i) - { - A(i + 1, i) += delta; - A(i + 1, i + 1) += -delta; + for (int i = 0; i < n - 1; ++i) + { + A(i + 1, i) += delta; + A(i + 1, i + 1) += -delta; - A(i, i) += -delta; - A(i, i + 1) += +delta; - } + A(i, i) += -delta; + A(i, i + 1) += +delta; + } - std::cout << A << std::endl << std::endl; + std::cout << "A = \n" + << A << std::endl + << std::endl; - Eigen::VectorXd b(n); - b.setZero(); + Eigen::VectorXd b(n); + b.setZero(); - b.head(n / 2).array() = 1.0; + b.head(n / 2).array() = 1.0; - std::ofstream f; - f.open("a.txt"); - for (int i = 0; i < 20; ++i) - { - std::cout << b.transpose() << std::endl; - f << b << std::endl << std::endl << std::endl; - b = A.colPivHouseholderQr().solve(b); - } + std::ofstream f; + f.open("implicit_sim.txt"); + for (int i = 0; i < 200; ++i) + { + // Solve A.T_{n+1} = T_{n} for T_{n+1} + f << b.transpose() << std::endl; + b = A.colPivHouseholderQr().solve(b); + } - std::cout << b.transpose() << std::endl; - return 0; -} + std::cout << b.transpose() << std::endl; + return 0; +} \ No newline at end of file diff --git a/exercises/eigen/movie.py b/exercises/eigen/movie.py new file mode 100644 index 0000000..01e149d --- /dev/null +++ b/exercises/eigen/movie.py @@ -0,0 +1,28 @@ +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.animation import FuncAnimation + +n = 20 +steps = 200 + +fig, ax = plt.subplots() +xdata = np.arange(0, n) +ydata = [] +ln, = ax.plot([], [], 'ro') + +data = np.loadtxt("implicit_sim.txt") +def init(): + ax.set_xlim(0,n) + ax.set_ylim(0,1.1) + return ln, + +def update(frame): + ydata = data[frame] + ln.set_data(xdata, ydata) + return ln, + +ani = FuncAnimation(fig, update, frames=np.arange(0, 200), + init_func=init, blit=True) + + +plt.show() diff --git a/exercises/eigen/sparse.cpp b/exercises/eigen/sparse.cpp index 66cded5..da6f7d9 100644 --- a/exercises/eigen/sparse.cpp +++ b/exercises/eigen/sparse.cpp @@ -1,42 +1,50 @@ #include #include +#include #include +#include int main() { - int n = 10; + int n = 20; + int steps = 200; - Eigen::SparseMatrix A; - A.resize(n, n); + Eigen::SparseMatrix A; + A.resize(n, n); - double delta = 0.4; + double delta = 0.4; - std::vector> fill; - fill.reserve(n * 4); + std::vector> fill; + fill.reserve(n * 4); - for (int i = 0; i < n - 1; ++i) - { - fill.push_back(Eigen::Triplet(i + 1, i, delta)); - fill.push_back(Eigen::Triplet(i + 1, i + 1, -delta)); - fill.push_back(Eigen::Triplet(i, i, 1.0 - delta)); - fill.push_back(Eigen::Triplet(i, i + 1, delta)); - } - fill.push_back(Eigen::Triplet(n - 1, n - 1, 1.0)); - A.setFromTriplets(fill.begin(), fill.end()); + for (int i = 0; i < n - 1; ++i) + { + fill.push_back(Eigen::Triplet(i + 1, i, delta)); + fill.push_back(Eigen::Triplet(i + 1, i + 1, -delta)); + fill.push_back(Eigen::Triplet(i, i, 1.0 - delta)); + fill.push_back(Eigen::Triplet(i, i + 1, delta)); + } + fill.push_back(Eigen::Triplet(n - 1, n - 1, 1.0)); + A.setFromTriplets(fill.begin(), fill.end()); - std::cout << A << std::endl; + std::cout << A << std::endl; - Eigen::VectorXd b(n); - b.head(n / 2).array() = 1.0; - b.tail(n / 2).array() = 0.0; + Eigen::VectorXd b(n); + b.head(n / 2).array() = 1.0; + b.tail(n / 2).array() = 0.0; - std::cout << b.transpose() << std::endl; + std::cout << b.transpose() << std::endl; - for (int i = 0; i < 100; ++i) - b = A * b; + std::ofstream f; + f.open("sparse_sim.txt"); + for (int i = 0; i < steps; ++i) + { + f << b.transpose() << std::endl; + b = A * b; + } - std::cout << b.transpose() << std::endl; + std::cout << b.transpose() << std::endl; - return 0; -} + return 0; +} \ No newline at end of file diff --git a/lectures/eigen/README.md b/lectures/eigen/README.md index c6fe6dd..ecc1039 100644 --- a/lectures/eigen/README.md +++ b/lectures/eigen/README.md @@ -1,23 +1,111 @@ template: titleslide -# Using Eigen -## Chris Richardson, Rupert Nash -## chris@bpi.cam.ac.uk, r.nash@epcc.ed.ac.uk -## CC-BY +# Linear Algebra for C++ (using Eigen) +## Joseph Lee, EPCC +## j.lee@epcc.ed.ac.uk --- -# What is Eigen3 and why use it? -- C++ library for matrix arithmetic -- “Header only” implementation: no libraries to compile and install (easy) -- Provides some “missing” features needed for scientific computing in C++ -- Contains optimisations to get good performance out of ARM & Intel processors -- Easy to use interface -- Support for dense and sparse matrices, vectors and “arrays” -- Some support for ‘solvers’ (A.x = b) -- Download from https://eigen.tuxfamily.org or e.g. `apt install libeigen3-dev` -- If you know Python, it is a bit like a NumPy for C++ +# Source +Original: +- Chris Richardson (chris@bpi.cam.ac.uk) +- Rupert Nash (r.nash@epcc.ed.ac.uk) --- -# Basics + +# Rundown + +- C++ Linear Algebra libraries: what and why? + +- Eigen3: + + - Getting started + + - Matrices: Dense, Sparse, Geometry + + - Solvers: Dense, Sparse, Iterative + + - General notes + + - Diffusion equation demo + exercise + + +--- + +# Linear Algebra: What do we need? + +__What for?__ + +- Scientific computing, PDE, simulations, finance, graphics, ML and more + +__What we expect?__ + +Numeric types: `complex` `complex` + +Object types: + +- Vectors + +- Matrix: (Dense / Sparse / Diagonal / Symmetric / Column/Row major) + +- Operations: Matrix (Multiplication/transpose/geometric rotation) + +- Solvers (`\(Ax = b\)`, Eigenvalues) : Direct (LU/QR factorization) / Indirect (CG) + +--- + +# Why use a library? + +- Simple & easy to use : clean API + +- Fast Performance + +- Correctness + +- Expressiveness + +- Don’t reinvent the wheel + +--- + +# Linear Algebra on C++ + +- Simple: **Eigen** + +- Big & popular: **PETSc** (C) + +- Standard: **LAPACK/LAPACK++/ScaLAPACK** + +- GPU: Nvidia/AMD **(cu/roc)(BLAS/SPARSE/SOLVER)** + +- CPU: **Intel MKL**, **OpenBLAS** + +- Future: C++26 std library? + +- vs: Python Numpy/Scipy + +- +- +--- + +# Eigen3 + +## Numpy, but C++ + +“Header only” no library to link + - `g++ -I /path/to/eigen/ my_program.cpp -o my_program` + +Contains optimisation for ARM & Intel processors (Vectorization - SIMD instructions) + +Easy to use interface + +Support for dense and sparse matrices, vectors, and “arrays” (coefficient-wise ops) + +Support for some solvers + +Download from or e.g. apt install libeigen3-dev + +--- + +# Matrix: Basic ```C++ #include @@ -51,18 +139,26 @@ std::cout << A.rows() << “ x ” << A.cols() << std::endl; So this is more like a 2D version of `std::vector` --- -# Convenience aliases +# Convenience typedefs ```C++ -using Eigen::Matrix3d = Eigen::Matrix -using Eigen::Matrix3i = Eigen::Matrix -using Eigen::MatrixXd = Eigen::Matrix -using Eigen::VectorXd = Eigen::Matrix -using Eigen::RowVectorXd = Eigen::Matrix - +MatrixNt = Matrix +MatrixXNt = Matrix +MatrixNXt = Matrix +VectorNt = Matrix +RowVectorNt = Matrix ``` +``` +N = {2, 3, 4, X = dynamic} +t = {i = int, f = float, d = double, cf = complex, cd = complex} +``` +e.g.: +- `Matrix3d` = 3x3 double matrix +- `Matrix3i` = 3x3 int matrix +- `MatrixXd` = (Dynamic)x(Dynamic) double matrix +- `VectorXd` = (Dynamic) double vector --- -# You can do Matrix arithmetic... +# Matrix arithmetics ```C++ Eigen::MatrixXd A(5, 10); @@ -73,7 +169,16 @@ Eigen::MatrixXd C = A * B; Eigen::VectorXd w = A * vec; ``` -Also dot and cross products for vectors, transpose, and usual scalar arithmetic `+ - * /` +More: + +Dot product: `v.dot(w)` + +Cross product: `v.cross(w)` + +Transpose: `A.transpose()` + +Set constant: +`A.setZeros(rows, cols)`, `A.setOnes(rows, cols)`, `A.setConstant(rows, cols, value)`, `A.setRandom(rows, cols)` --- # Element-wise ops with `Array`s @@ -107,16 +212,147 @@ a_eigen(10, 0) = 1.0; Eigen::Map a2_eigen(a.data(), 10, 100); ``` +--- +# Sparse matrix: Basic + +When dealing with very large matrices with many zeros (e.g. from Differential Equations), store as sparse matrices + +```C++ +#include + +SparseMatrix > mat(1000,2000); +// declares a 1000x2000 column-major compressed sparse matrix of complex + +SparseMatrix mat(1000,2000); +// declares a 1000x2000 row-major compressed sparse matrix of double + +SparseVector > vec(1000); +// declares a column sparse vector of complex of size 1000 + +SparseVector vec(1000); +// declares a row sparse vector of double of size 1000 +``` + +--- +# Sparse matrix +Simplest way to create a sparse matrix is to build a triplet: +```C++ +typedef Triplet T; +std::vector tripletList; +tripletList.reserve(estimation_of_entries); +for(...) +{ +// ... +tripletList.push_back(T(i,j,v_ij)); +} +SparseMatrixType m(rows,cols); +m.setFromTriplets(tripletList.begin(), tripletList.end()); +// m is ready to go! +``` +Operators: +- Sparse: `SM.transpose()`, `SM.adjoint()`… +- Sparse-Sparse: `SM1 + SM2`, `SM1 *SM2` … +- Sparse-Dense: `DM2 = DM1 + SM1` , `DM2 = SM1* DM2` + +--- +# Matrix: Geometry + +```C++ +#include +``` +- 2D Rotations: +```C++ + Rotation2D rot2(angle_in_radian) + ``` +- 3D Rotaions: +```C++ +AngleAxis AA(angle_in_radian, Vector3f(ax,ay,az)) +``` +- Quarternion, ND-Transformations… + +--- +# Solvers + +Simple example: `\(Ax = b\)` +```C++ +#include +#include + +int main() +{ + Eigen::Matrix3f A; + Eigen::Vector3f b; + A << 1,2,3, 4,5,6, 7,8,10; + b << 3, 3, 4; + std::cout << "Here is the matrix A:\n" << A << std::endl; + std::cout << "Here is the vector b:\n" << b << std::endl; + Eigen::Vector3f x = A.colPivHouseholderQr().solve(b); #HERE + std::cout << "The solution is:\n" << x << std::endl; +} +``` +2 steps: +- Decompose +- Solve + +--- +# Decompositions: + +E.g. +- PartialPivLU + +- FullPivLU + +- HouseholderQR + +- ColPivHousehoulderQR + +Varying matrix requirements, speed, reliability, accuracy, optimization + + + + +--- +# Other solvers +### Singular values +E.g. JacobiSVD, BDCSVD + +### Eigenvalues/vectors +E.g. SelfAdjointEigenSolver, ComplexEigenSolver + +### Sparse Solver +E.g. SparseLU, SparseQR, SimplicialLLT, SimplicialLDLT + +--- + +# Iterative Solvers + +```C++ +#include +``` +Useful for solving `\(Ax = b\)` where `\(A\)` is large and sparse + +E.g. +- ConjugateGradient + +- LeastSquaresConjugateGradient + +- BICGSTAB + +Use with preconditioner --- -# Efficiency: Eigen does lots of checks in debug mode +# General Notes: +Eigen does a lot of checks in debug mode: + +- Turn optimisation on: `-O2` etc -Turn optimisation on: `-O2` etc. +- Turn off debug: `-DNDEBUG` -Turn off debug: `-DNDEBUG` +Can use external BLAS/LAPACK library (e.g. MKL, OpenBLAS) by enabling macro +- e.g. `EIGEN_USE_BLAS` --- -# Walk through example +# Demo: diffusion equation Solve diffusion equation @@ -124,16 +360,16 @@ $$\frac{\partial T}{\partial t} = k \frac{\partial^2 T}{\partial x^2}$$ in 1D using an explicit method -Each timestep can be solved by using `\(T_{n+1} = A T_n\)` +Each timestep can be iterated by using `\(T_{n+1} = A T_n\)` -1. Create an initial vector for `T` -2. Create a dense matrix for `A` -3. Matrix multiply several times -- Convert to an implicit method: `\(A.T_{n+1} = T_n\)` -- Sparse matrices +1. Create an initial vector for `\(T\)` + +2. Create a dense matrix for `\(A\)` + +3. Matrix multiply `\(A\)` several times -Type along with me - see `exercises/eigen` +See `exercises/eigen/explicit.cpp` --- # Diffusion equation (explicit) @@ -150,9 +386,6 @@ Left-hand side is unknown (next time step) Let: `\(\delta = k\Delta t/ \Delta x^2\)` ---- -# Matrix - ``` A = [ 1-ẟ ẟ 0 0 0 … ] [ ẟ 1-2ẟ ẟ 0 0 … ] @@ -161,9 +394,9 @@ A = [ 1-ẟ ẟ 0 0 0 … ] [ 0 0 0 ẟ 1-2ẟ ẟ 0 …] ``` - --- # Diffusion equation (implicit) +More stable: Subscript means time, square brackets position @@ -172,7 +405,37 @@ Subscript means time, square brackets position Left-hand side is unknown (next time step) `$$A.T_{n+1} = T_n$$` -Let: $$\delta = k \Delta t/\Delta x^2$$ - -??? +``` +A = [ 1+ẟ -ẟ 0 0 0 …] + [ -ẟ 1+2ẟ -ẟ 0 0 …] + [ 0 -ẟ 1+2ẟ -ẟ 0 …] + [ 0 0 -ẟ 1+2ẟ -ẟ …] + [ 0 0 0 -ẟ 1+2ẟ …] + +``` The matrix A is very similar - just flip the sign of the delta terms + +--- +# Exercise: Diffusion equation (sparse) + +Hints: + +```C++ +#include +``` +```C++ +std::vector> fill; + fill.reserve(...); +``` +```C++ +for (int i = 0; i < n - 1; ++i) + { + fill.push_back(Eigen::Triplet(...); + ... + } +``` +```C++ +A.setFromTriplets(fill.begin(), fill.end()); +``` + +See `exercises/eigen/sparse.cpp` \ No newline at end of file diff --git a/lectures/eigen/index.html b/lectures/eigen/index.html index 0b675fd..d5bd99d 100644 --- a/lectures/eigen/index.html +++ b/lectures/eigen/index.html @@ -14,7 +14,7 @@