Skip to content

Commit 2b1137f

Browse files
rrsettgastCusiniMFrancois HamonMatteo Cusini
authored
Add symmetric operations to tensors (#216)
* Change Rij_eq_AiBj to have symmetric Voigt result * add unscaled symmetric outer product * scaledAdd for matrices. * update unitTests * Adding Rij_eq_AkiAkj * added unit test for Rij_eq_AkiAkj Co-authored-by: Matteo Cusini <[email protected]> Co-authored-by: Francois Hamon <[email protected]> Co-authored-by: Matteo Cusini <[email protected]>
1 parent aa7765b commit 2b1137f

File tree

7 files changed

+455
-78
lines changed

7 files changed

+455
-78
lines changed

‎src/fixedSizeSquareMatrixOps.hpp‎

Lines changed: 43 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -176,6 +176,49 @@ void Rij_eq_AikSymBklAjl( DST_SYM_MATRIX && LVARRAY_RESTRICT_REF dstSymMatrix,
176176
symMatrixB );
177177
}
178178

179+
180+
/**
181+
* @brief Perform the outer product of @p vectorA with itself writing the result to @p dstMatrix.
182+
* @tparam M The size of both dimensions of @p dstMatrix and the length of @p vectorA.
183+
* @tparam DST_MATRIX The type of @p dstMatrix.
184+
* @tparam VECTOR_A The type of @p vectorA.
185+
* @param dstMatrix The matrix the result is written to, of size M x M.
186+
* @param vectorA The first vector in the outer product, of length M.
187+
* @details Performs the operations @code dstMatrix[ i ][ j ] = vectorA[ i ] * vectorA[ j ] @endcode
188+
*/
189+
template< std::ptrdiff_t M, typename DST_SYM_MATRIX, typename VECTOR_A >
190+
LVARRAY_HOST_DEVICE CONSTEXPR_WITHOUT_BOUNDS_CHECK inline
191+
void symRij_eq_AiAj( DST_SYM_MATRIX && LVARRAY_RESTRICT_REF dstMatrix,
192+
VECTOR_A const & LVARRAY_RESTRICT_REF vectorA )
193+
{
194+
internal::SquareMatrixOps< M >::symRij_eq_AiAj( std::forward< DST_SYM_MATRIX >( dstMatrix ),
195+
vectorA );
196+
}
197+
198+
/**
199+
* @brief Perform the unscaled symmetric outer product of @p vectorA and
200+
* @p vectorB with itself writing the result to @p dstMatrix.
201+
* @tparam M The size of both dimensions of @p dstMatrix and the length of @p vectorA.
202+
* @tparam DST_MATRIX The type of @p dstMatrix.
203+
* @tparam VECTOR_A The type of @p vectorA.
204+
* @tparam VECTOR_B The type of @p vectorB.
205+
* @param dstMatrix The matrix the result is written to, of size M x N.
206+
* @param vectorA The first vector in the outer product, of length M.
207+
* @param vectorB The second vector in the outer product, of length M.
208+
* @details Performs the operations @code dstMatrix[ i ][ j ] = vectorA[ i ] * vectorB[ j ] + vectorA[ j ] * vectorB[ i ] @endcode
209+
*/
210+
template< std::ptrdiff_t M, typename DST_SYM_MATRIX, typename VECTOR_A, typename VECTOR_B >
211+
LVARRAY_HOST_DEVICE CONSTEXPR_WITHOUT_BOUNDS_CHECK inline
212+
void symRij_eq_AiBj_plus_AjBi( DST_SYM_MATRIX && LVARRAY_RESTRICT_REF dstMatrix,
213+
VECTOR_A const & LVARRAY_RESTRICT_REF vectorA,
214+
VECTOR_B const & LVARRAY_RESTRICT_REF vectorB )
215+
{
216+
internal::SquareMatrixOps< M >::symRij_eq_AiBj_plus_AjBi( std::forward< DST_SYM_MATRIX >( dstMatrix ),
217+
vectorA,
218+
vectorB );
219+
}
220+
221+
179222
/**
180223
* @return Return the determinant of the symmetric matrix @p symMatrix.
181224
* @tparam SYM_MATRIX The type of @p symMatrix.

‎src/fixedSizeSquareMatrixOpsImpl.hpp‎

Lines changed: 103 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -331,6 +331,55 @@ struct SquareMatrixOps< 2 >
331331
matrixA[ 1 ][ 1 ] * symMatrixB[ 1 ] * matrixA[ 0 ][ 1 ];
332332
}
333333

334+
335+
/**
336+
* @brief Perform the outer product of @p vectorA with itself writing the result to @p dstMatrix.
337+
* @tparam M The size of both dimensions of @p dstMatrix and the length of @p vectorA.
338+
* @tparam DST_MATRIX The type of @p dstMatrix.
339+
* @tparam VECTOR_A The type of @p vectorA.
340+
* @param dstMatrix The matrix the result is written to, of size M x N.
341+
* @param vectorA The first vector in the outer product, of length M.
342+
* @details Performs the operations @code dstMatrix[ i ][ j ] = vectorA[ i ] * vectorA[ j ] @endcode
343+
*/
344+
template< typename DST_MATRIX, typename VECTOR_A >
345+
LVARRAY_HOST_DEVICE CONSTEXPR_WITHOUT_BOUNDS_CHECK inline
346+
static void symRij_eq_AiAj( DST_MATRIX && LVARRAY_RESTRICT_REF dstMatrix,
347+
VECTOR_A const & LVARRAY_RESTRICT_REF vectorA )
348+
{
349+
internal::checkSizes< 3 >( dstMatrix );
350+
internal::checkSizes< 2 >( vectorA );
351+
352+
dstMatrix[ 0 ] = vectorA[ 0 ] * vectorA[ 0 ];
353+
dstMatrix[ 1 ] = vectorA[ 1 ] * vectorA[ 1 ];
354+
dstMatrix[ 2 ] = vectorA[ 0 ] * vectorA[ 1 ];
355+
}
356+
357+
/**
358+
* @brief Perform the unscaled symmetric outer product of @p vectorA and
359+
* @p vectorB writing the result to @p dstMatrix.
360+
* @tparam DST_MATRIX The type of @p dstMatrix.
361+
* @tparam VECTOR_A The type of @p vectorA.
362+
* @tparam VECTOR_B The type of @p vectorB.
363+
* @param dstMatrix The matrix the result is written to, of size M x N.
364+
* @param vectorA The first vector in the outer product, of length M.
365+
* @param vectorB The second vector in the outer product, of length M.
366+
* @details Performs the operations @code dstMatrix[ i ][ j ] = vectorA[ i ] * vectorB[ j ] + vectorA[ j ] * vectorB[ i ] @endcode
367+
*/
368+
template< typename DST_SYM_MATRIX, typename VECTOR_A, typename VECTOR_B >
369+
LVARRAY_HOST_DEVICE CONSTEXPR_WITHOUT_BOUNDS_CHECK inline
370+
static void symRij_eq_AiBj_plus_AjBi( DST_SYM_MATRIX && LVARRAY_RESTRICT_REF dstMatrix,
371+
VECTOR_A const & LVARRAY_RESTRICT_REF vectorA,
372+
VECTOR_B const & LVARRAY_RESTRICT_REF vectorB )
373+
{
374+
internal::checkSizes< 3 >( dstMatrix );
375+
internal::checkSizes< 2 >( vectorA );
376+
internal::checkSizes< 2 >( vectorB );
377+
378+
dstMatrix[ 0 ] = 2 * vectorA[ 0 ] * vectorB[ 0 ];
379+
dstMatrix[ 1 ] = 2 * vectorA[ 1 ] * vectorB[ 1 ];
380+
dstMatrix[ 2 ] = vectorA[ 0 ] * vectorB[ 1 ] + vectorA[ 1 ] * vectorB[ 0 ];
381+
}
382+
334383
/**
335384
* @brief Compute the eigenvalues of the symmetric matrix @p symMatrix.
336385
* @tparam DST_VECTOR The type of @p eigenvalues.
@@ -890,6 +939,60 @@ struct SquareMatrixOps< 3 >
890939
matrixA[ 0 ][ 2 ] * symMatrixB[ 2 ] * matrixA[ 1 ][ 2 ];
891940
}
892941

942+
/**
943+
* @brief Perform the outer product of @p vectorA with itself writing the result to @p dstMatrix.
944+
* @tparam M The size of both dimensions of @p dstMatrix and the length of @p vectorA.
945+
* @tparam DST_MATRIX The type of @p dstMatrix.
946+
* @tparam VECTOR_A The type of @p vectorA.
947+
* @param dstMatrix The matrix the result is written to, of size M x N.
948+
* @param vectorA The first vector in the outer product, of length M.
949+
* @details Performs the operations @code dstMatrix[ i ][ j ] = vectorA[ i ] * vectorA[ j ] @endcode
950+
*/
951+
template< typename DST_MATRIX, typename VECTOR_A >
952+
LVARRAY_HOST_DEVICE CONSTEXPR_WITHOUT_BOUNDS_CHECK inline
953+
static void symRij_eq_AiAj( DST_MATRIX && LVARRAY_RESTRICT_REF dstMatrix,
954+
VECTOR_A const & LVARRAY_RESTRICT_REF vectorA )
955+
{
956+
internal::checkSizes< 6 >( dstMatrix );
957+
internal::checkSizes< 3 >( vectorA );
958+
959+
dstMatrix[ 0 ] = vectorA[ 0 ] * vectorA[ 0 ];
960+
dstMatrix[ 1 ] = vectorA[ 1 ] * vectorA[ 1 ];
961+
dstMatrix[ 2 ] = vectorA[ 2 ] * vectorA[ 2 ];
962+
dstMatrix[ 3 ] = vectorA[ 1 ] * vectorA[ 2 ];
963+
dstMatrix[ 4 ] = vectorA[ 0 ] * vectorA[ 2 ];
964+
dstMatrix[ 5 ] = vectorA[ 0 ] * vectorA[ 1 ];
965+
}
966+
967+
/**
968+
* @brief Perform the unscaled symmetric outer product of @p vectorA and
969+
* @p vectorB writing the result to @p dstMatrix.
970+
* @tparam DST_MATRIX The type of @p dstMatrix.
971+
* @tparam VECTOR_A The type of @p vectorA.
972+
* @tparam VECTOR_B The type of @p vectorB.
973+
* @param dstMatrix The matrix the result is written to, of size M x N.
974+
* @param vectorA The first vector in the outer product, of length M.
975+
* @param vectorB The second vector in the outer product, of length M.
976+
* @details Performs the operations @code dstMatrix[ i ][ j ] = vectorA[ i ] * vectorB[ j ] + vectorA[ j ] * vectorB[ i ] @endcode
977+
*/
978+
template< typename DST_SYM_MATRIX, typename VECTOR_A, typename VECTOR_B >
979+
LVARRAY_HOST_DEVICE CONSTEXPR_WITHOUT_BOUNDS_CHECK inline
980+
static void symRij_eq_AiBj_plus_AjBi( DST_SYM_MATRIX && LVARRAY_RESTRICT_REF dstMatrix,
981+
VECTOR_A const & LVARRAY_RESTRICT_REF vectorA,
982+
VECTOR_B const & LVARRAY_RESTRICT_REF vectorB )
983+
{
984+
internal::checkSizes< 6 >( dstMatrix );
985+
internal::checkSizes< 3 >( vectorA );
986+
internal::checkSizes< 3 >( vectorB );
987+
988+
dstMatrix[ 0 ] = 2 * vectorA[ 0 ] * vectorB[ 0 ];
989+
dstMatrix[ 1 ] = 2 * vectorA[ 1 ] * vectorB[ 1 ];
990+
dstMatrix[ 2 ] = 2 * vectorA[ 2 ] * vectorB[ 2 ];
991+
dstMatrix[ 3 ] = vectorA[ 1 ] * vectorB[ 2 ] + vectorA[ 2 ] * vectorB[ 1 ];
992+
dstMatrix[ 4 ] = vectorA[ 0 ] * vectorB[ 2 ] + vectorA[ 2 ] * vectorB[ 0 ];
993+
dstMatrix[ 5 ] = vectorA[ 0 ] * vectorB[ 1 ] + vectorA[ 1 ] * vectorB[ 0 ];
994+
}
995+
893996
/**
894997
* @brief Compute the eigenvalues of the symmetric matrix @p symMatrix.
895998
* @tparam DST_VECTOR The type of @p eigenvalues.

‎src/genericTensorOps.hpp‎

Lines changed: 69 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -613,6 +613,38 @@ void scaledAdd( DST_VECTOR && LVARRAY_RESTRICT_REF dstVector,
613613
}
614614
}
615615

616+
617+
/**
618+
* @brief Add @p srcMatrix scaled by @p scale to @p dstMatrix.
619+
* @tparam ISIZE The size of the first dimension of @p dstMatrix and @p srcMatrix.
620+
* @tparam JSIZE The size of the second dimension of @p dstMatrix and @p srcMatrix.
621+
* @tparam DST_MATRIX The type of @p dstMatrix.
622+
* @tparam SRC_MATRIX The type of @p srcMatrix.
623+
* @param dstMatrix The destination matrix, of size ISIZE x N.
624+
* @param srcMatrix The source matrix, of size ISIZE x N.
625+
* @param scale The value to scale the entries of @p srcMatrix by.
626+
* @details Performs the operation @code dstMatrix[ i ][ j ] += srcMatrix[ i ][ j ] @endcode
627+
*/
628+
template< std::ptrdiff_t ISIZE, std::ptrdiff_t JSIZE, typename DST_MATRIX, typename SRC_MATRIX >
629+
LVARRAY_HOST_DEVICE CONSTEXPR_WITHOUT_BOUNDS_CHECK inline
630+
void scaledAdd( DST_MATRIX && LVARRAY_RESTRICT_REF dstMatrix,
631+
SRC_MATRIX const & LVARRAY_RESTRICT_REF srcMatrix,
632+
std::remove_reference_t< decltype( srcMatrix[ 0 ][ 0 ] ) > const scale )
633+
{
634+
static_assert( ISIZE > 0, "ISIZE must be greater than zero." );
635+
static_assert( JSIZE > 0, "JSIZE must be greater than zero." );
636+
internal::checkSizes< ISIZE, JSIZE >( dstMatrix );
637+
internal::checkSizes< ISIZE, JSIZE >( srcMatrix );
638+
639+
for( std::ptrdiff_t i = 0; i < ISIZE; ++i )
640+
{
641+
for( std::ptrdiff_t j = 0; j < JSIZE; ++j )
642+
{
643+
dstMatrix[ i ][ j ] += scale * srcMatrix[ i ][ j ];
644+
}
645+
}
646+
}
647+
616648
/**
617649
* @brief Multiply the elements of @p vectorA and @p vectorB putting the result into @p dstVector.
618650
* @tparam ISIZE The length of @p dstVector and @p srcVector.
@@ -796,33 +828,6 @@ void Rij_eq_AiBj( DST_MATRIX && LVARRAY_RESTRICT_REF dstMatrix,
796828
}
797829
}
798830

799-
/**
800-
* @brief Perform the outer product of @p vectorA with itself writing the result to @p dstMatrix.
801-
* @tparam M The size of both dimensions of @p dstMatrix and the length of @p vectorA.
802-
* @tparam DST_MATRIX The type of @p dstMatrix.
803-
* @tparam VECTOR_A The type of @p vectorA.
804-
* @param dstMatrix The matrix the result is written to, of size M x N.
805-
* @param vectorA The first vector in the outer product, of length M.
806-
* @details Performs the operations @code dstMatrix[ i ][ j ] = vectorA[ i ] * vectorA[ j ] @endcode
807-
*/
808-
template< std::ptrdiff_t M, typename DST_MATRIX, typename VECTOR_A >
809-
LVARRAY_HOST_DEVICE CONSTEXPR_WITHOUT_BOUNDS_CHECK inline
810-
void Rij_eq_AiAj( DST_MATRIX && LVARRAY_RESTRICT_REF dstMatrix,
811-
VECTOR_A const & LVARRAY_RESTRICT_REF vectorA )
812-
{
813-
static_assert( M > 0, "M must be greater than zero." );
814-
internal::checkSizes< M, M >( dstMatrix );
815-
internal::checkSizes< M >( vectorA );
816-
817-
for( std::ptrdiff_t i = 0; i < M; ++i )
818-
{
819-
for( std::ptrdiff_t j = 0; j < M; ++j )
820-
{
821-
dstMatrix[ i ][ j ] = vectorA[ i ] * vectorA[ j ];
822-
}
823-
}
824-
}
825-
826831
/**
827832
* @brief Perform the outer product of @p vectorA and @p vectorB adding the result to @p dstMatrix.
828833
* @tparam ISIZE The size of the first dimension of @p dstMatrix and the length of @p vectorA.
@@ -1290,6 +1295,43 @@ void Rij_eq_AkiBkj( DST_MATRIX && LVARRAY_RESTRICT_REF dstMatrix,
12901295
}
12911296
}
12921297

1298+
/**
1299+
* @brief Multiply the transpose of @p matrixA with @p matrixA and put the result into @p dstMatrix.
1300+
* @tparam ISIZE The size of both dimensions of @p dstMatrix and the second dimension of @p matrixA.
1301+
* @tparam JSIZE The size of the first dimension of @p matrixA.
1302+
* @tparam DST_MATRIX The type of @p dstMatrix.
1303+
* @tparam MATRIX_A The type of @p matrixA.
1304+
* @param dstMatrix The matrix the result is written to, of size ISIZE x ISIZE.
1305+
* @param matrixA The left matrix in the multiplication, of size JSIZE x ISIZE.
1306+
* @details Performs the operation
1307+
* @code dstMatrix[ i ][ j ] = matrixA[ k ][ i ] * matrixA[ k ][ j ] @endcode
1308+
*/
1309+
template< std::ptrdiff_t ISIZE,
1310+
std::ptrdiff_t JSIZE,
1311+
typename DST_MATRIX,
1312+
typename MATRIX_A >
1313+
LVARRAY_HOST_DEVICE CONSTEXPR_WITHOUT_BOUNDS_CHECK inline
1314+
void Rij_eq_AkiAkj( DST_MATRIX && LVARRAY_RESTRICT_REF dstMatrix,
1315+
MATRIX_A const & LVARRAY_RESTRICT_REF matrixA )
1316+
{
1317+
static_assert( ISIZE > 0, "ISIZE must be greater than zero." );
1318+
static_assert( JSIZE > 0, "JSIZE must be greater than zero." );
1319+
internal::checkSizes< ISIZE, ISIZE >( dstMatrix );
1320+
internal::checkSizes< JSIZE, ISIZE >( matrixA );
1321+
1322+
for( std::ptrdiff_t i = 0; i < ISIZE; ++i )
1323+
{
1324+
for( std::ptrdiff_t j = 0; j < ISIZE; ++j )
1325+
{
1326+
dstMatrix[ i ][ j ] = matrixA[ 0 ][ i ] * matrixA[ 0 ][ j ];
1327+
for( std::ptrdiff_t k = 1; k < JSIZE; ++k )
1328+
{
1329+
dstMatrix[ i ][ j ] = dstMatrix[ i ][ j ] + matrixA[ k ][ i ] * matrixA[ k ][ j ];
1330+
}
1331+
}
1332+
}
1333+
}
1334+
12931335
/**
12941336
* @brief Multiply the transpose of @p matrixA with @p matrixB and add the result into @p dstMatrix.
12951337
* @tparam ISIZE The size of the first dimension of @p dstMatrix and the second dimension of @p matrixA.

0 commit comments

Comments
 (0)