Should Aither Use Eigen?
Aither, like many CFD codes requires matrix and vector operations when calculating the flow solution. When using block implicit methods, Aither uses 5x5 matrices for the flow equations and 2x2 matrices if a turbulence model is selected. These matrices are mutiplied with other matrices, multiplied with vectors, scaled by scalar values, added, subtracted, and inverted. If these operations could be performed more efficiently, it would result in a great performance improvement. There are many third party linear algebra libraries available that could be used by Aither such as Eigen, Armadillo, and PETSc. PETSc is widely used and has the ability to run in parallel. It contains many robust matrix solvers as well. However it is written in C, must be linked to, and is probably overkill for small matrix/vector operations as described above. Eigen is also widely used and writen in C++. It has the added advantage of being entirely header-only, so there is no need to link to anything. It could therefore be distributed with the Aither source code eliminating the issue of findind a dependency on various computer systems. Armadillo is another linear algebra library written in C++. It links with libraries such as LAPACK, OpenBLAS, MKL, or ATLAS. Since Eigen claims comprable or better performance than many linear algebra libraries, and it has the best ease of use, it was choosen for evalution. Armadillo may be evaluated at a later point.
The following five test cases were used to evaluate Eigen versus the linear algebra code already in Aither.
- Matrix-matrix multiplication
- Matrix-vector multiplication
- Matrix multiplication with a scalar and addition
- Vector multiplication with a scalar and addition
- Matrix inverse
The above tests were repeated 10 million times each using 5x5 and 2x2 matrices. Eigen has predefined classes for small (< 4)
matrices and vectors which statically allocate their memory on the stack. These can be quite a bit faster than the more general n-dimensional
matrices and vectors which dynamically allocate their memory on the heap. Aither requires general sized matrices because the matrix size is
determined at run time. For scalar implicit methods like LU-SGS and DPLUR the matrix size is 1. For block implicit methods like
BLU-SGS and BDPLUR the matrix size is 5. For the tests using the 2x2 matrices both the static
Eigen::Matrix2d and dynamic
versions were used, but the comparison to Aither was made using the heap version. Aither uses dynamic allocation for all matrices
squareMatrix(2), but uses static allocation for vectors
genArray. In Aither, all vectors are of size 7 which is the
maximum number of equations solved.
The timing results for each of the tests are shown below in seconds. There is no data for the vector multiplication with a scalar and addtion test for Aither for the size 2 vector because in Aither all vectors are the same length. Therefore the test would be the same as for the larger vector.
|Matrix Type||Size||Matrix-Matrix Multiplication||Matrix-Vector Multiplication||Matrix Scale & Addition||Vector Scale & Addition||Matrix Inverse|
Surprisingly Aither outperforms Eigen in all tests with the exception of the matrix multiplication with a scalar and addition test. The statically allocated Eigen matrix and vector classes far outperform both the dynamically allocated Eigen classes and the Aither classes. This is expected as memory access is faster to the stack than it is to the heap. Looking at the Eigen benchmarks, the best performance is expected for larger matrices than the ones tested here. The results of these tests indicate that it would not be worthwhile to integrate Eigen into Aither. The source code for these tests can be found here.