8. `LvArray::tensorOps`¶

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.

8.1. Object Representation¶

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 \times 2$ symmetric matrix is represented as vector of length three as $[a_{00}, a_{11}, a_{01}]$ and that a $3 \times 3$ symmetric matrix is represented as vector of length six as $[a_{00}, a_{11}, a_{22}, a_{12}, a_{02}, a_{01}]$ . One consequence of this is that you can use the vector methods to for example add one symmetric matrix to another.

8.2. Common operations¶

Variables
$\alpha$ : A scalar.
$x$ : A vector in $\mathbb{R}^m \times \mathbb{R}^1$ .
$y$ : A vector in $\mathbb{R}^m \times \mathbb{R}^1$ .
$z$ : A vector in $\mathbb{R}^m \times \mathbb{R}^1$ .
$\mathbf{A}$ : A matrix in $\mathbb{R}^m \times \mathbb{R}^n$ .
$\mathbf{B}$ : A matrix in $\mathbb{R}^m \times \mathbb{R}^n$ .
Operation	Function
$x \leftarrow y$	`tensorOps::copy< m >( x, y )`
$\mathbf{A} \leftarrow \mathbf{B}$	`tensorOps::copy< m, n >( A, B )`
$x \leftarrow \alpha x$	`tensorOps::scale< m >( x, alpha )`
$\mathbf{A} \leftarrow \alpha \mathbf{A}$	`tensorOps::scale< m, n >( A, alpha )`
$x \leftarrow \alpha y$	`tensorOps::scaledCopy< m >( x, y, alpha )`
$\mathbf{A} \leftarrow \alpha \mathbf{B}$	`tensorOps::scaledCopy< m, n >( A, B, alpha )`
$x \leftarrow x + y$	`tensorOps::add< m >( x, y )`
$\mathbf{A} \leftarrow \mathbf{A} + \mathbf{B}$	`tensorOps::add< m, n >( A, B, alpha )`
$x \leftarrow x - y$	`tensorOps::subtract< m >( x, y )`
$x \leftarrow x + \alpha y$	`tensorOps::scaledAdd< m >( x, y, alpha )`
$x \leftarrow y \circ 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 ).

8.3. Vector operations¶

Variables
$\alpha$ : A scalar.
$x$ : A vector in $\mathbb{R}^m \times \mathbb{R}^1$ .
$y$ : A vector in $\mathbb{R}^m \times \mathbb{R}^1$ .
$z$ : A vector in $\mathbb{R}^m \times \mathbb{R}^1$ .
Operation	Function
$\|x\|_\infty$	`tensorOps::maxAbsoluteEntry< m >( x )`
$\|x\|_2^2$	`tensorOps::l2NormSquared< m >( x )`
$\|x\|_2$	`tensorOps::l2Norm< m >( x )`
$x \leftarrow \hat{x}$	`tensorOps::normalize< m >( x )`
$x^T y$	`tensorOps::AiBi( x, y )`

If $x$ , $y$ and $z$ are all in $\mathbb{R}^3 \times \mathbb{R}^1$ then you can perform the operation $x \leftarrow y \times z$ as tensorOps::crossProduct( x, y, z ).

8.4. Matrix vector operations¶

Variables
$x$ : A vector in $\mathbb{R}^m \times \mathbb{R}^1$ .
$y$ : A vector in $\mathbb{R}^n \times \mathbb{R}^1$ .
$\mathbf{A}$ : A matrix in $\mathbb{R}^m \times \mathbb{R}^n$ .
Operation	Function
$\mathbf{A} \leftarrow x y^T$	`tensorOps::Rij_eq_AiBj< m, n >( A, x, y )`
$\mathbf{A} \leftarrow \mathbf{A} + x y^T$	`tensorOps::Rij_add_AiBj< m, n >( A, x, y )`
$x \leftarrow \mathbf{A} y$	`tensorOps::Ri_eq_AijBj< m, n >( x, A, y )`
$x \leftarrow x + \mathbf{A} y$	`tensorOps::Ri_add_AijBj< m, n >( x, A, y )`
$y \leftarrow \mathbf{A}^T x$	`tensorOps::Ri_eq_AjiBj< n, m >( y, A, x )`
$y \leftarrow y + \mathbf{A}^T x$	`tensorOps::Ri_add_AjiBj< n, m >( y, A, x )`

8.5. Matrix operations¶

Variables
$\mathbf{A}$ : A matrix in $\mathbb{R}^m \times \mathbb{R}^n$ .
$\mathbf{B}$ : A matrix in $\mathbb{R}^m \times \mathbb{R}^p$ .
$\mathbf{C}$ : A matrix in $\mathbb{R}^p \times \mathbb{R}^n$ .
$\mathbf{D}$ : A matrix in $\mathbb{R}^n \times \mathbb{R}^m$ .
$\mathbf{E}$ : A matrix in $\mathbb{R}^m \times \mathbb{R}^m$ .
Operation Function
$\mathbf{A} \leftarrow \mathbf{D}^T$	`tensorOps::transpose< m, n >( A, D )`
$\mathbf{A} \leftarrow \mathbf{B} \mathbf{C}$	`tensorOps::Rij_eq_AikBkj< m, n, p >( A, B, C )`
$\mathbf{A} \leftarrow \mathbf{A} + \mathbf{B} \mathbf{C}$	`tensorOps::Rij_add_AikBkj< m, n, p >( A, B, C )`
$\mathbf{B} \leftarrow \mathbf{A} \mathbf{C}^T$	`tensorOps::Rij_eq_AikBjk< m, p, n >( B, A, C )`
$\mathbf{B} \leftarrow \mathbf{B} + \mathbf{A} \mathbf{C}^T$	`tensorOps::Rij_add_AikBjk< m, p, n >( B, A, C )`
$\mathbf{E} \leftarrow \mathbf{E} + \mathbf{A} \mathbf{A}^T$	`tensorOps::Rij_add_AikAjk< m, n >( E, A )`
$\mathbf{C} \leftarrow \mathbf{B}^T \mathbf{A}$	`tensorOps::Rij_eq_AkiBkj< p, n, m >( C, B, A )`
$\mathbf{C} \leftarrow \mathbf{C} + \mathbf{B}^T \mathbf{A}$	`tensorOps::Rij_add_AkiBkj< p, n, m >( C, B, A )`

8.6. Square matrix operations¶

Variables
$\alpha$ : A scalar.
$\mathbf{A}$ : A matrix in $\mathbb{R}^m \times \mathbb{R}^m$ .
$\mathbf{B}$ : A matrix in $\mathbb{R}^m \times \mathbb{R}^m$ .
Operation	Function
$\mathbf{A} \leftarrow \mathbf{A}^T$	`tensorOps::transpose< m >( A )`
$\mathbf{A} \leftarrow \mathbf{A} + \alpha \mathbf{I}$	`tensorOps::tensorOps::addIdentity< m >( A, alpha )`
$tr(\mathbf{A})$	`tensorOps::trace< m >( A )`
$\|\mathbf{A}\|$	`tensorOps::determinant< m >( A )`
$\mathbf{A} \leftarrow \mathbf{B}^-1$	`tensorOps::invert< m >( A, B )`
$\mathbf{A} \leftarrow \mathbf{A}^-1$	`tensorOps::invert< m >( A )`

Note

Apart from tensorOps::determinant and tensorOps::invert are only implemented for matrices of size $2 \times 2$ and $3 \times 3$ .

8.7. Symmetric matrix operations¶

Variables
$\alpha$ : A scalar.
$x$ : A vector in $\mathbb{R}^m \times \mathbb{R}^1$ .
$y$ : A vector in $\mathbb{R}^m \times \mathbb{R}^1$ .
$\mathbf{A}$ : A matrix in $\mathbb{R}^m \times \mathbb{R}^m$ .
$\mathbf{B}$ : A matrix in $\mathbb{R}^m \times \mathbb{R}^m$ .
$\mathbf{S}$ : A symmetric matrix in $\mathbb{R}^m \times \mathbb{R}^m$ .
$\mathbf{Q}$ : A symmetric matrix in $\mathbb{R}^m \times \mathbb{R}^m$ .
Operation	Function
$\mathbf{S} \leftarrow \mathbf{S} + \alpha \mathbf{I}$	`tensorOps::symAddIdentity< m >( S, alpha )`
$tr(\mathbf{S})$	`tensorOps::symTrace< m >( S )`
$x \leftarrow \mathbf{S} y$	`tensorOps::Ri_eq_symAijBj< m >( x, S ,y )`
$x \leftarrow x + \mathbf{S} y$	`tensorOps::Ri_add_symAijBj< m >( x, S ,y )`
$\mathbf{A} \leftarrow \mathbf{S} \mathbf{B}^T$	`tensorOps::Rij_eq_symAikBjk< m >( A, S, B )`
$\mathbf{S} \leftarrow \mathbf{A} \mathbf{Q} \mathbf{A}^T$	`tensorOps::Rij_eq_AikSymBklAjl< m >( S, A, Q )`
$\|\mathbf{S}\|$	`tensorOps::symDeterminant< m >( S )`
$\mathbf{S} \leftarrow \mathbf{Q}^-1$	`tensorOps::symInvert< m >( S, Q )`
$\mathbf{S} \leftarrow \mathbf{S}^-1$	`tensorOps::symInvert< m >( S )`
$x \leftarrow \mathbf{S} \mathbf{A}^T = diag(x) \mathbf{A}^T$	`tensorOps::symEigenvalues< M >( x, S )`
$x, \mathbf{A} \leftarrow \mathbf{S} \mathbf{A}^T = diag(x) \mathbf{A}^T$	`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 \times 2$ and $3 \times 3$ .

8.8. Examples¶

TEST( tensorOps, AiBi )
{
  double const x[ 3 ] = { 0, 1, 2 };
  double const y[ 3 ] = { 1, 2, 3 };

  // Works with c-arrays
  EXPECT_EQ( LvArray::tensorOps::AiBi< 3 >( x, y ), 8 );

  LvArray::Array< double,
                  1,
                  camp::idx_seq< 0 >,
                  std::ptrdiff_t,
                  LvArray::MallocBuffer > xArray( 3 );

  xArray( 0 ) = 0;
  xArray( 1 ) = 1;
  xArray( 2 ) = 2;

  // Works with LvArray::Array.
  EXPECT_EQ( LvArray::tensorOps::AiBi< 3 >( xArray, y ), 8 );

  // Works with LvArray::ArrayView.
  EXPECT_EQ( LvArray::tensorOps::AiBi< 3 >( xArray.toView(), y ), 8 );
  EXPECT_EQ( LvArray::tensorOps::AiBi< 3 >( xArray.toViewConst(), y ), 8 );

  // Create a 2D array with fortran layout.
  LvArray::Array< double,
                  2,
                  camp::idx_seq< 1, 0 >,
                  std::ptrdiff_t,
                  LvArray::MallocBuffer > yArray( 1, 3 );

  yArray( 0, 0 ) = 1;
  yArray( 0, 1 ) = 2;
  yArray( 0, 2 ) = 3;

  // Works with non contiguous slices.
  EXPECT_EQ( LvArray::tensorOps::AiBi< 3 >( x, yArray[ 0 ] ), 8 );
  EXPECT_EQ( LvArray::tensorOps::AiBi< 3 >( xArray, yArray[ 0 ] ), 8 );
}

[Source: examples/exampleTensorOps.cpp]

You can mix and match the data types of the objects and also call the tensorOps methods on device.

CUDA_TEST( tensorOps, device )
{
  // Create an array of 5 3x3 symmetric matrices.
  LvArray::Array< int,
                  2,
                  camp::idx_seq< 1, 0 >,
                  std::ptrdiff_t,
                  LvArray::ChaiBuffer > symmetricMatrices( 5, 6 );

  int offset = 0;
  for( int & value : symmetricMatrices )
  {
    value = offset++;
  }

  LvArray::ArrayView< int const,
                      2,
                      0,
                      std::ptrdiff_t,
                      LvArray::ChaiBuffer > symmetricMatricesView = symmetricMatrices.toViewConst();

  // The tensorOps methods work on device.
  RAJA::forall< RAJA::cuda_exec< 32 > >(
    RAJA::TypedRangeSegment< std::ptrdiff_t >( 0, symmetricMatricesView.size( 0 ) ),
    [symmetricMatricesView] __device__ ( std::ptrdiff_t const i )
  {
    double result[ 3 ];
    float x[ 3 ] = { 1.3f, 2.2f, 5.3f };

    // You can mix value types.
    LvArray::tensorOps::Ri_eq_symAijBj< 3 >( result, symmetricMatricesView[ i ], x );

    LVARRAY_ERROR_IF_NE( result[ 0 ], x[ 0 ] * symmetricMatricesView( i, 0 ) +
                         x[ 1 ] * symmetricMatricesView( i, 5 ) +
                         x[ 2 ] * symmetricMatricesView( i, 4 ) );
  }
    );
}

[Source: examples/exampleTensorOps.cpp]

8.9. Bounds checking¶

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.

TEST( tensorOps, boundsCheck )
{
  double x[ 3 ] = { 0, 1, 2 };
  LvArray::tensorOps::normalize< 3 >( x );

  // This would fail to compile since x is not length 4.
  // LvArray::tensorOps::normalize< 4 >( x );

  LvArray::Array< double,
                  1,
                  camp::idx_seq< 0 >,
                  std::ptrdiff_t,
                  LvArray::MallocBuffer > xArray( 8 );

  xArray( 0 ) = 10;

#if defined(LVARRAY_BOUNDS_CHECK)
  // This will fail at runtime.
  EXPECT_DEATH_IF_SUPPORTED( LvArray::tensorOps::normalize< 3 >( xArray ), "" );
#endif

  int matrix[ 2 ][ 2 ] = { { 1, 0 }, { 0, 1 } };
  LvArray::tensorOps::normalize< 2 >( matrix[ 0 ] );

  // This would fail to compile since normalize expects a vector
  // LvArray::tensorOps::normalize< 4 >( matrix );

  LvArray::Array< double,
                  2,
                  camp::idx_seq< 0, 1 >,
                  std::ptrdiff_t,
                  LvArray::MallocBuffer > matrixArray( 2, 2 );

  matrixArray( 0, 0 ) = 1;
  matrixArray( 0, 1 ) = 1;

  LvArray::tensorOps::normalize< 2 >( matrixArray[ 0 ] );

  // This will also fail to compile
  // LvArray::tensorOps::normalize< 2 >( matrixArray );
}

[Source: examples/exampleTensorOps.cpp]

8.10. Doxygen¶

LvArray::tensorOps

8. LvArray::tensorOps¶

8.1. Object Representation¶

8.2. Common operations¶

8.3. Vector operations¶

8.4. Matrix vector operations¶

8.5. Matrix operations¶

8.6. Square matrix operations¶

8.7. Symmetric matrix operations¶

8.8. Examples¶

8.9. Bounds checking¶

8.10. Doxygen¶

8. `LvArray::tensorOps`¶