This is a light weight Linear Algebra C++17 library I am currently developing for my own stuff.
Early development stage -> Do not use it right now!
The goals are:
- a good compromise between performance and library complexity
- shorter compilation time than libs like Eigen or Blaze.
- a concise wrapping of libraries like Blas, Blis, Lapacke, etc.
As this contructions is only available for “storable” vectors and not
views they are defined using functions (and not usual constructors).
#include "LinearAlgebra/vector.hpp"
#include <array>
#include <vector>
using namespace LinearAlgebra;
int
main()
{
// Construction of a dense vector using move(std::vector)
//
std::vector<int> data_1(4, 1);
auto v_1 = create_vector(std::move(data_1));
assert(data_1.size() == 0);
static_assert(std::is_same_v<Vector<int>, decltype(v_1)>);
std::cout << v_1 << std::endl;
//================
// Construction of a Tiny_Vector from a std::array copy
//
std::array<int, 4> data_2;
for (auto& data_2_i : data_2) data_2_i = 2;
auto v_2 = create_vector(data_2);
static_assert(std::is_same_v<Tiny_Vector<int, 4>, decltype(v_2)>);
std::cout << v_2 << std::endl;
//================
// Construction of a Tiny_Vector from a C array
//
auto v_3 = create_vector({3, 3, 3, 3});
static_assert(std::is_same_v<Tiny_Vector<int, 4>, decltype(v_3)>);
std::cout << v_3 << std::endl;
}prints:
1 1 1 1 2 2 2 2 3 3 3 3[2020-05-08 Fri 12:03] <- lower/upper bound
#include "LinearAlgebra/vector.hpp"
#include "LinearAlgebra/utils/lower_upper_bound.hpp"
using namespace LinearAlgebra;
int
main()
{
double data[13] = {1, 1, 2, 3, 3, 3, 3, 4, 4, 4, 5, 5, 6};
auto X = create_vector_view(data, 13);
std::cout << X << std::endl;
auto X_3 = create_vector_view(X, lower_bound(X, 3), upper_bound(X, 3));
X_3=10;
std::cout << X << std::endl;
}prints:
1 1 2 3 3 3 3 4 4 4 5 5 6 1 1 2 10 10 10 10 4 4 4 5 5 6
#include "LinearAlgebra/matrix.hpp"
using namespace LinearAlgebra;
int
main()
{
Tiny_Matrix<int, 5, 6> mat;
mat = 1;
auto view_mat = create_matrix_view(mat, 1, 3, 2, 5);
view_mat = -1;
std::cout << "full matrix:" << mat << std::endl;
//================
Tiny_Symmetric_Matrix<int, 6> mat_S;
create_matrix_view_full(mat_S) = 0; // fill the full matrix
mat_S = 1; // here as mat_S is symmetric
// only fill lower triangular part is filled
// Take a subview of a _symmetric matrix_
auto view_mat_S = create_matrix_view(mat_S, 1, 3, 1, 3);
view_mat_S = -1; // here the subview mat_S is symmetric too and
// only the lower triangular part will be filled
std::cout << "full matrix:" << create_matrix_view_full(mat_S) << std::endl;
std::cout << "submatrix:" << mat_S << std::endl;
}prints:
full matrix:
1 1 1 1 1 1
1 1 -1 -1 -1 1
1 1 -1 -1 -1 1
1 1 1 1 1 1
1 1 1 1 1 1
full matrix:
1 0 0 0 0 0
1 -1 0 0 0 0
1 -1 -1 0 0 0
1 1 1 1 0 0
1 1 1 1 1 0
1 1 1 1 1 1
submatrix:
1 X X X X X
1 -1 X X X X
1 -1 -1 X X X
1 1 1 1 X X
1 1 1 1 1 X
1 1 1 1 1 1
You can define your own matrix type specialization:
#include "LinearAlgebra/matrix.hpp"
using namespace LinearAlgebra;
template <typename T>
using Matrix_nx2 =
Default_Matrix<T,
Matrix_Special_Structure_Enum::None, // Dense matrix
Matrix_Storage_Mask_Enum::None, //
std::size_t, // Dynamic number of rows
std::integral_constant<std::size_t, 2>, // Static number of columns, here 2
std::integral_constant<std::size_t, 10>>; // Static leading dimension, here 10
// note: having a static leading dimension=10 allows to statically
// allocate matrix, but this also limits max row size to 10.
static_assert(not Has_Static_I_Size_v<Matrix_nx2<int>>);
static_assert(Has_Static_J_Size_v<Matrix_nx2<int>>);
static_assert(Has_Static_Capacity_v<Matrix_nx2<int>>);
static_assert(Has_Static_Capacity_v<Tiny_Matrix<int, 2, 3>>);
static_assert(Has_Static_Capacity_v<Tiny_Vector<int, 2>>);
int
main()
{
// Note:
//
// 1. M(5,3) would lead to run-time error in debug mode
// Reason: request 3 columns, 2 expected
//
// 2. M(15, 2) would lead to run-time error in debug mode
// Reason: number of rows > statically defined leading dimension
//
Matrix_nx2<double> M(5, 2);
M = 0;
// Note: one can also write
// "create_vector_view_matrix_row(M,2)=2;"
auto row_2 = create_vector_view_matrix_row(M, 2);
row_2 = 2;
auto col_1 = create_vector_view_matrix_column(M, 1);
col_1 = 1;
static_assert(Has_Static_Dimension_v<decltype(row_2)>);
static_assert(not Has_Static_Dimension_v<decltype(col_1)>);
std::cout << M;
return 0;
}Easier interfacing with non generic code (.cpp files etc…).
- [ ] add some examples
This is only for demo purpose as a real implementation makes more sense with sparse matrices + a preconditioner.
#include "LinearAlgebra/matrix.hpp"
#include "LinearAlgebra/vector.hpp"
using namespace LinearAlgebra;
// Basic CG implementation
// https://en.wikipedia.org/wiki/Conjugate_gradient_method
//
template <typename A_IMPL, typename X0_IMPL, typename B_IMPL>
bool
cg(const Matrix_Crtp<A_IMPL>& A, Dense_Vector_Crtp<X0_IMPL>& X0, const Dense_Vector_Crtp<B_IMPL>& b)
{
// Sanity checks
//
assert(all_sizes_are_equal_p(A.I_size(), A.J_size(), X0.size(), b.size()));
static_assert(Is_Symmetric_Matrix_v<A_IMPL> or Is_Hermitian_Matrix_v<A_IMPL>);
// Parameters
//
const double eps = 1e-6;
const double squared_eps = eps * eps;
const size_t max_iter = 100;
// Working vector type
//
using element_type = Common_Element_Type_t<A_IMPL, X0_IMPL, B_IMPL>;
auto r = similar(Type_v<element_type>, X0);
auto p = similar(r);
auto Ap = similar(r);
// Initialization
//
r = b - A * X0;
auto squared_norm_r_old = dot(r, r);
if (squared_norm_r_old < squared_eps)
{
return true;
}
p = r;
// Main loop
//
for (size_t i = 0; i < max_iter; i++)
{
Ap = A * p;
auto alpha = squared_norm_r_old / dot(p, Ap);
X0 = X0 + alpha * p;
r = r - alpha * Ap;
auto squared_norm_r_new = dot(r, r);
std::cout << "iter " << i << " residue " << squared_norm_r_new << std::endl;
if (squared_norm_r_new < squared_eps)
{
return true;
}
p = r + squared_norm_r_new / squared_norm_r_old * p;
squared_norm_r_old = squared_norm_r_new;
}
return false;
}
int
main()
{
Symmetric_Matrix<double> M(10, 10);
Vector<double> X0(10);
Vector<double> b(10);
M = 1;
M(6, 5) = 5; // Note: in debug mode M(5, 6) = 2 would lead to an
// assert failure as by default symmetric matrices are
// stored into their lower part.
create_vector_view_matrix_diagonal(M) = 10;
b = 1;
X0 = 0;
bool status = cg(M, X0, b);
std::cout << X0 << std::endl;
std::cout << std::boolalpha << status << std::endl;
}prints:
iter 0 residue 0.0652995 iter 1 residue 1.03037e-32 0.0543933 0.0543933 0.0543933 0.0543933 0.0543933 0.0376569 0.0376569 0.0543933 0.0543933 0.0543933 true
#include "LinearAlgebra/lapack/lapack.hpp"
#include "LinearAlgebra/matrix.hpp"
#include "LinearAlgebra/vector.hpp"
using namespace LinearAlgebra;
int
main()
{
const size_t n = 5;
Vector<double> v(n), w(n);
v = 1;
Symmetric_Matrix<double> M(n, n);
M = 0;
create_vector_view_matrix_diagonal(M) = 10;
M(4, 0) = 5;
w = M * v;
std::cout << "symmetric M :" << M << std::endl << std::endl;
std::cout << "v :" << v << std::endl << std::endl;
std::cout << "w = M.v :" << w << std::endl << std::endl;
// Low level call of lapack: L*L^t decomposition
//
int info = Lapack::potrf(M);
assert(info == 0);
// Create a constant view defining L
// (symmetric matrix uses lower part)
static_assert(Is_Lower_Matrix_Storage_v<decltype(M)>);
//
auto L = create_matrix_view_lower_triangular(M.as_const());
std::cout << "L :" << L << std::endl << std::endl;
// inplace solve of w = M.v = L.L^t.v ...
// ... at the end w "contains" v
//
w = inverse(L) * w; // L^(-1).w = L^t.v
w = inverse(transpose(L)) * w; // L^(-t).L^(-1).w = v
std::cout << "v such that w=M.v :" << w << std::endl;
}prints:
symmetric M :
10 X X X X
0 10 X X X
0 0 10 X X
0 0 0 10 X
5 0 0 0 10
v :
1
1
1
1
1
w = M.v :
15
10
10
10
15
L :
3.16228 X X X X
0 3.16228 X X X
0 0 3.16228 X X
0 0 0 3.16228 X
1.58114 0 0 0 2.73861
v such that w=M.v :
1
1
1
1
1
For the moment I only have defined dense matrices (BLAS compatible with column major order):
#include "LinearAlgebra/dense/matrix.hpp"
#include "LinearAlgebra/dense/vector.hpp"
#include <iostream>
using namespace LinearAlgebra;
int
main()
{
Matrix<double> M_1(4, 5);
Symmetric_Matrix<int> M_2(4, 4);
Tiny_Strict_Lower_Triangular_Matrix<float, 4, 7> M_3;
std::cout << M_1 << std::endl;
std::cout << M_2 << std::endl;
std::cout << M_3 << std::endl;
}prints
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 X X X 0 0 X X 0 0 0 X 0 0 0 0 X X X X X X X 0 X X X X X X 0 0 X X X X X 0 0 0 X X X X
The generic definition for these matrix types is:
template <T, // is the component type
SPECIAL_STRUCTURE, // is in {None, Symmetric, Hermitian, Triangular,
// Unit_Triangular, Triangular_Strict}
MASK, // is in {None, Upper, Upper_Strict, Lower, Lower_Strict }
N_TYPE, // std::size_t or a std::integral_constant<std::size_t,N>
M_TYPE, // std::size_t or a std::integral_constant<std::size_t,M>
LEADING_DIMENSION // std::size_t or a std::integral_constant<std::size_t,LD>
> //
class Default_Matrix; // or {Default_Matrix_View, Default_Matrix_Const_View}There are some alias covering the usual cases:
| Dynamic | Static |
|---|---|
| Matrix<T> M(I_size, J_size); | Tiny_Matrix<T, I_SIZE, J_SIZE> M; |
| Symmetric_Matrix<T> M(I_size, J_size); | Tiny_Symmetric_Matrix<T, SIZE> M; |
| Hermitian_Matrix<T> M(I_size, J_size); | Tiny_Hermitian_Matrix<T, SIZE> M; |
| Lower_Triangular_Matrix<T> M(I_size, J_size); | Tiny_Lower_Triangular_Matrix<T, I_SIZE, J_SIZE> M; |
| Upper_Triangular_Matrix<T> M(I_size, J_size); | Tiny_Upper_Triangular_Matrix<T, I_SIZE, J_SIZE> M; |
| Lower_Triangular_Strict_Matrix<T> M(I_size, J_size); | Tiny_Lower_Triangular_Strict_Matrix<T, I_SIZE, J_SIZE> M; |
| Upper_Triangular_Strict_Matrix<T> M(I_size, J_size); | Tiny_Upper_Triangular_Strict_Matrix<T, I_SIZE, J_SIZE> M; |
| Lower_Unit_Triangular_Matrix<T> M(I_size, J_size); | Tiny_Lower_Unit_Triangular_Matrix<T, I_SIZE, J_SIZE> M; |
| Upper_Unit_Triangular_Matrix<T> MI_size, J_size); | Tiny_Upper_Unit_Triangular_Matrix<T, I_SIZE, J_SIZE> M; |
Please note that by default Symmetric/Hermitian matrices are stored in their Lower part.
For each case you can also use views, there are two types of view: mutable one and constant one. For instance:
Matrix<double> M(10, 5);
auto view = view_as_lower_triangular_strict(M.as_const());will return a constant view (a lightweight matrix where only pointers are stored and not owned).
File: LinearAlgebra/utils/size_utils.hpp
Has_Static_Capacity_v:check if capacity is static or not (a static capacity means that the object can be created without dynamic memory allocation)static_assert(not Has_Static_Capacity_v<Vector<double>>); static_assert(Has_Static_Capacity_v<Tiny_Matrix<int, 3, 4>>);
Has_Static_Size_v:check if a vector has a static sizestatic_assert(Has_Static_Size_v<Tiny_Vector<double, 3>>); static_assert(not Has_Static_Size_v<Vector<double>>); static_assert(not Has_Static_Size_v<Matrix<double>>); // safely // returns // false, even // if not a // vector type
Any_Has_Static_Size_v:check if any vector has a static sizestatic_assert(Any_Has_Static_Size_v<Vector<double>, Tiny_Vector<double, 3>>);
Has_Static_I_Size_v:check if a matrix has a staticI_size(number of rows):static_assert(Has_Static_I_Size_v<Tiny_Matrix<int, 3, 4>>); static_assert(not Has_Static_I_Size_v<Matrix<int>>);
Has_Static_I_Size_v:check if a matrix has a staticI_size(number of rows):static_assert(Has_Static_I_Size_v<Tiny_Matrix<int, 3, 4>>); static_assert(not Has_Static_I_Size_v<Matrix<int>>);
Any_Static_I_Size_v:check if any matrix has a staticI_size(number of rows):static_assert(Any_Has_Static_I_Size_v<Vector<int>, Tiny_Matrix<int, 3, 4>>);
Has_Static_J_Size_v:same thanI_sizebut for columnsAny_Static_J_Size_v:same thanI_sizebut for columns
TODO: move to function with links:
Two functions, declared as constexpr:
get_static_size_if_any(...)that returns a static size if any.CAVEAT: it does not check if all sizes are equal!
all_sizes_are_equal_p(...)check that all sizes are equal
File: LinearAlgebra/utils/sfinae_vmt_helpers.hpp
Note: these predicates will also work for sparse matrices (they are not restricted to dense one).
Checks if the matrix is defined by its lower or upper part:
template <typename MATRIX>
constexpr bool Is_Upper_Matrix_Storage_v = Is_Upper_Matrix_Storage<MATRIX>::value;
template <typename MATRIX>
constexpr bool Is_Upper_Strict_Matrix_Storage_v = Is_Upper_Strict_Matrix_Storage<MATRIX>::value;
template <typename MATRIX>
constexpr bool Is_Lower_Matrix_Storage_v = Is_Lower_Matrix_Storage<MATRIX>::value;
template <typename MATRIX>
constexpr bool Is_Lower_Strict_Matrix_Storage_v = Is_Lower_Strict_Matrix_Storage<MATRIX>::value;Checks for special structure:
template <typename MATRIX>
constexpr bool Is_Full_Matrix_v = Is_Full_Matrix<MATRIX>::value;
template <typename MATRIX>
constexpr bool Is_Symmetric_Matrix_v = Is_Symmetric_Matrix<MATRIX>::value;
template <typename MATRIX>
constexpr bool Is_Hermitian_Matrix_v = Is_Hermitian_Matrix<MATRIX>::value;
template <typename MATRIX>
constexpr bool Is_Triangular_Matrix_v = Is_Triangular_Matrix<MATRIX>::value;
template <typename MATRIX>
constexpr bool Is_Unit_Triangular_Matrix_v = Is_Unit_Triangular_Matrix<MATRIX>::value;Usage example:
static_assert(Is_Symmetric_Matrix_v<Symmetric_Matrix<double>>);
static_assert(Is_Lower_Matrix_Storage_v<Symmetric_Matrix<double>>);Test if two objects are aliased or not.
bool status = are_not_aliased_p(vector0, matrix1)Note: it is easy/fast to check if two objects are not aliased (memory blocks don’t overlap). It can be more tricky (with all the possible memory increments/strides) to check if two objects are aliased. That’s the reason why we also have this function:
bool status = are_maybe_aliased_p(vector0, matrix1)Also see: Same mathematical object
A predicate that checks if two objects represent exactly the same “mathematical” object. By exactly we mean:
- same memory
- same dimension
This is possible, think to a matrix and its view.
bool status = same_mathematical_object_p(matrix_1,matrix_2);Note: according to the definition, two identical mathematical objects are trivially Aliased.
[2020-04-24 Fri 12:50] <- CopyCreate an uninitialized object of the same type and same dimension:
auto u = similar(v);a variant is to change type:
Vector<double> v;
auto int_u = similar(Type_v<int>,v);Like Similar but also performs a copy of the element
auto u = copy(v);
auto int_u = copy(Type_v<int>,v);These functions mimic std::lower_bound, std::upper_bound but work with indices.
// Returns the first index *idx* such "value <= x[idx]"
// or X.size() if such element index does not exist
//
template <typename IMPL>
std::size_t
lower_bound(const Dense_Vector_Crtp<IMPL>& X, const Element_Type_t<IMPL>& value);
// Returns the first index *idx* such "value < x[idx]"
// or X.size() if such element index does not exist
//
template <typename IMPL>
std::size_t
upper_bound(const Dense_Vector_Crtp<IMPL>& X, const Element_Type_t<IMPL>& value);They can be used to create view, see View from lower/upper bounds
In general vectors or matrices cannot be resized.
This avoids introducing an asymmetry in the code between dynamic & static size objects. This asymmetry would have come with some extra complications both for the developer and the user who want to implement some generic routines.
Only a reduced number of expressions are supported (TODO: list them!).
By example you can write
V=2*transpose(M)*U+2*Vas this expression can be directly mapped to a Blas subroutine.
However, you cannot write, in full generality, things like:
V=2*transpose(M)*M**M*U+2*VPlease note that, all in all, this constraint has some positive side effects as it reduces the “chance” of introducing hidden temporary creations.
Also note that beside Expression Template you can call available expressions using reverse polish notation, by example
V = 2 * transpose(M) * U + 3 * Vcan be computed by calling:
assign(V, _plus_, _product_, _product_, 2, _transpose_, M, U, _product_, 3, _lhs_);For the moment I have not introduced all the machinery that “manually”
generates simd code. This would have required a lot to be introduced,
like cpu dependent simd definitions, aligned or packed template
parameters etc… Really not my priority for the moment… Some other
libs do a great job here.