LvArray::tensorOps is a collection of free template functions which perform common linear algebra operations on compile time sized matrices and vectors. All of the functions work on device.
The LvArray::tensorOps methods currently operate on four different types of object; scalars, vectors, matrices, and symmetric matrices.
Vectors are represented by either a one dimensional c-array or a one dimensional LvArray::Array object. Examples of acceptable vector types are
double[ 4 ]int[ 55 ]LvArray::Array< float, 1, camp::idx_seq< 0 >, std::ptrdiff_t, LvArray::MallocBuffer >LvArray::ArrayView< long, 1, 0, std::ptrdiff_t, LvArray::MallocBuffer >LvArray::ArraySlice< double, 1, 0, std::ptrdiff_t >LvArray::ArraySlice< double, 1, -1, std::ptrdiff_t >
Matrices are represented by either a two dimensional c-array or a two dimensional LvArray::Array object. Examples of acceptable matrix types
double[ 3 ][ 3 ]int[ 5 ][ 2 ]LvArray::Array< float, 2, camp::idx_seq< 0, 1 >, std::ptrdiff_t, LvArray::MallocBuffer >LvArray::ArrayView< long, 2, 0, std::ptrdiff_t, LvArray::MallocBuffer >:LvArray::ArraySlice< double, 2, 1, std::ptrdiff_t >LvArray::ArraySlice< double, 2, 0, std::ptrdiff_t >LvArray::ArraySlice< double, 2, -1, std::ptrdiff_t >
Symmetric matrices are represented in Voigt notation as a vector. This means that a 2 × 2 symmetric matrix is represented as vector of length three as [a00, a11, a01] and that a 3 × 3 symmetric matrix is represented as vector of length six as [a00, a11, a22, a12, a02, a01]. One consequence of this is that you can use the vector methods to for example add one symmetric matrix to another.
| Variables |
|---|
| α: A scalar. |
| --------------------------------------------------------------------------------------------------------------- |
| x: A vector in ℝm × ℝ1. |
| --------------------------------------------------------------------------------------------------------------- |
| y: A vector in ℝm × ℝ1. |
| --------------------------------------------------------------------------------------------------------------- |
| z: A vector in ℝm × ℝ1. |
| --------------------------------------------------------------------------------------------------------------- |
| A: A matrix in ℝm × ℝn. |
| --------------------------------------------------------------------------------------------------------------- |
| B: A matrix in ℝm × ℝn. |
| --------------------------------------------------------------------------------------------------------------- |
| Operation Function |
| ====================================================== ======================================================== |
x ← y tensorOps::copy< m >( x, y ) |
A ← B tensorOps::copy< m, n >( A, B ) |
x ← αx tensorOps::scale< m >( x, alpha ) |
A ← αA tensorOps::scale< m, n >( A, alpha ) |
x ← αy tensorOps::scaledCopy< m >( x, y, alpha ) |
A ← αB tensorOps::scaledCopy< m, n >( A, B, alpha ) |
x ← x + y tensorOps::add< m >( x, y ) |
A ← A + B tensorOps::add< m, n >( A, B, alpha ) |
x ← x − y tensorOps::subtract< m >( x, y ) |
x ← x + αy tensorOps::scaledAdd< m >( x, y, alpha ) |
x ← y ∘ z tensorOps::hadamardProduct< m >( x, y, z ) |
There is also a function fill that will set all the entries of the object to a specific value. For a vector it is called as tensorOps::fill< n >( x ) and for a matrix as tensorOps::fill< m, n >( A ).
| Variables |
|---|
| α: A scalar. |
| ------------------------------------------------------------------------------ |
| x: A vector in ℝm × ℝ1. |
| ------------------------------------------------------------------------------ |
| y: A vector in ℝm × ℝ1. |
| ------------------------------------------------------------------------------ |
| z: A vector in ℝm × ℝ1. |
| ------------------------------------------------------------------------------ |
| Operation Function |
| ==================================== ========================================= |
|x|∞ tensorOps::maxAbsoluteEntry< m >( x ) |
|x|22 tensorOps::l2NormSquared< m >( x ) |
|x|2 tensorOps::l2Norm< m >( x ) |
x ← x̂ tensorOps::normalize< m >( x ) |
xTy tensorOps::AiBi( x, y ) |
If x, y and z are all in ℝ3 × ℝ1 then you can perform the operation x ← y × z as tensorOps::crossProduct( x, y, z ).
| Variables |
|---|
| x: A vector in ℝm × ℝ1. |
| -------------------------------------------------------------------------------------------------- |
| y: A vector in ℝn × ℝ1. |
| -------------------------------------------------------------------------------------------------- |
| A: A matrix in ℝm × ℝn. |
| -------------------------------------------------------------------------------------------------- |
| Operation Function |
| ====================================================== =========================================== |
A ← xyT tensorOps::Rij_eq_AiBj< m, n >( A, x, y ) |
A ← A + xyT tensorOps::Rij_add_AiBj< m, n >( A, x, y ) |
x ← Ay tensorOps::Ri_eq_AijBj< m, n >( x, A, y ) |
x ← x + Ay tensorOps::Ri_add_AijBj< m, n >( x, A, y ) |
y ← ATx tensorOps::Ri_eq_AjiBj< n, m >( y, A, x ) |
y ← y + ATx tensorOps::Ri_add_AjiBj< n, m >( y, A, x ) |
| Variables |
|---|
| A: A matrix in ℝm × ℝn. |
| ------------------------------------------------------------------------------------------------------------------- |
| B: A matrix in ℝm × ℝp. |
| ------------------------------------------------------------------------------------------------------------------- |
| C: A matrix in ℝp × ℝn. |
| ------------------------------------------------------------------------------------------------------------------- |
| D: A matrix in ℝn × ℝm. |
| ------------------------------------------------------------------------------------------------------------------- |
| E: A matrix in ℝm × ℝm. |
| ------------------------------------------------------------------------------------------------------------------- |
| Operation Function |
| =================================================================== =============================================== |
A ← DT tensorOps::transpose< m, n >( A, D ) |
A ← BC tensorOps::Rij_eq_AikBkj< m, n, p >( A, B, C ) |
A ← A + BC tensorOps::Rij_add_AikBkj< m, n, p >( A, B, C ) |
B ← ACT tensorOps::Rij_eq_AikBjk< m, p, n >( B, A, C ) |
B ← B + ACT tensorOps::Rij_add_AikBjk< m, p, n >( B, A, C ) |
E ← E + AAT tensorOps::Rij_add_AikAjk< m, n >( E, A ) |
C ← BTA tensorOps::Rij_eq_AkiBkj< p, n, m >( C, B, A ) |
C ← C + BTA tensorOps::Rij_add_AkiBkj< p, n, m >( C, B, A ) |
| Variables |
|---|
| α: A scalar. |
| ----------------------------------------------------------------------------------------------------------------------- |
| A: A matrix in ℝm × ℝm. |
| ----------------------------------------------------------------------------------------------------------------------- |
| B: A matrix in ℝm × ℝm. |
| ----------------------------------------------------------------------------------------------------------------------- |
| Operation Function |
| =============================================================== ======================================================= |
A ← AT tensorOps::transpose< m >( A ) |
A ← A + αI tensorOps::tensorOps::addIdentity< m >( A, alpha ) |
tr(A) tensorOps::trace< m >( A ) |
|A| tensorOps::determinant< m >( A ) |
A ← B−1 tensorOps::invert< m >( A, B ) |
A ← A−1 tensorOps::invert< m >( A ) |
Note
Apart from tensorOps::determinant and tensorOps::invert are only implemented for matrices of size 2 × 2 and 3 × 3.
| Variables |
|---|
| α: A scalar. |
| ------------------------------------------------------------------------------------------------------------------------------- |
| x: A vector in ℝm × ℝ1. |
| ------------------------------------------------------------------------------------------------------------------------------- |
| y: A vector in ℝm × ℝ1. |
| ------------------------------------------------------------------------------------------------------------------------------- |
| A: A matrix in ℝm × ℝm. |
| ------------------------------------------------------------------------------------------------------------------------------- |
| B: A matrix in ℝm × ℝm. |
| ------------------------------------------------------------------------------------------------------------------------------- |
| S: A symmetric matrix in ℝm × ℝm. |
| ------------------------------------------------------------------------------------------------------------------------------- |
| Q: A symmetric matrix in ℝm × ℝm. |
| ------------------------------------------------------------------------------------------------------------------------------- |
| Operation Function |
| ================================================================================ ============================================== |
S ← S + αI tensorOps::symAddIdentity< m >( S, alpha ) |
tr(S) tensorOps::symTrace< m >( S ) |
x ← Sy tensorOps::Ri_eq_symAijBj< m >( x, S ,y ) |
x ← x + Sy tensorOps::Ri_add_symAijBj< m >( x, S ,y ) |
A ← SBT tensorOps::Rij_eq_symAikBjk< m >( A, S, B ) |
S ← AQAT tensorOps::Rij_eq_AikSymBklAjl< m >( S, A, Q ) |
|S| tensorOps::symDeterminant< m >( S ) |
S ← Q−1 tensorOps::symInvert< m >( S, Q ) |
S ← S−1 tensorOps::symInvert< m >( S ) |
x ← SAT = diag(x)AT tensorOps::symEigenvalues< M >( x, S ) |
x, A ← SAT = diag(x)AT tensorOps::symEigenvectors< M >( x, S ) |
There are also two function tensorOps::denseToSymmetric and tensorOps::symmetricToDense which convert between dense and symmetric matrix representation.
Note
Apart from tensorOps::symAddIdentity and tensorOps::symTrace the symmetric matrix operations are only implemented for matrices of size 2 × 2 and 3 × 3.
../../examples/exampleTensorOps.cpp
[Source: examples/exampleTensorOps.cpp]
You can mix and match the data types of the objects and also call the tensorOps methods on device.
../../examples/exampleTensorOps.cpp
[Source: examples/exampleTensorOps.cpp]
Whatever the argument type the number of dimensions is checked at compile time. For example if you pass a double[ 3 ][ 3 ] or a three dimensional LvArray::ArraySlice to LvArray::tensorOps::crossProduct you will get a compilation error since that function is only implemented for vectors. When passing a c-array as an argument the size of the array is checked at compile time. For example if you pass int[ 2 ][ 3 ] to LvArray::tensorOps::addIdentity you will get a compilation error because that function only operates on square matrices. However when passing an LvArray::Array* object the size is only checked at runtime if LVARRAY_BOUNDS_CHECK is defined.
../../examples/exampleTensorOps.cpp
[Source: examples/exampleTensorOps.cpp]