MatrixUtils.cc File Reference

#include "MatrixUtils.hh"
#include <cmath>
#include "corecel/math/Algorithms.hh"
#include "corecel/math/ArrayUtils.hh"
#include "corecel/math/SoftEqual.hh"

Typedefs
using	Mat3 = celeritas::SquareMatrixReal3

Functions
template<class T >
T	celeritas::determinant (SquareMatrix< T, 3 > const &mat)
	Calculate the determiniant of a 3x3 matrix.
template<class T >
T	celeritas::trace (SquareMatrix< T, 3 > const &mat)
	Calculate the trace of a 3x3 matrix. More...
template<class T , size_type N>
SquareMatrix< T, N >	celeritas::gemm (SquareMatrix< T, N > const &a, SquareMatrix< T, N > const &b)
	Naive square matrix-matrix multiply. More...
template<class T , size_type N>
SquareMatrix< T, N >	celeritas::gemm (matrix::TransposePolicy, SquareMatrix< T, N > const &a, SquareMatrix< T, N > const &b)
	Naive square matrix-matrix multiply with the first matrix transposed. More...
template<class T , size_type N>
void	celeritas::orthonormalize (SquareMatrix< T, N > *mat)
	Normalize and orthogonalize a small, dense matrix. More...
Mat3	celeritas::make_rotation (Real3 const &ax, Turn theta)
	Create a C-ordered rotation matrix from an arbitrary rotation. More...
Mat3	celeritas::make_rotation (Axis ax, Turn theta)
	Create a C-ordered rotation matrix.
Mat3	celeritas::make_rotation (Axis ax, Turn theta, Mat3 const &other)
	Rotate a C-ordered rotation matrix. More...
SquareMatrixReal3	celeritas::make_transpose (SquareMatrixReal3 const &mat)
	Construct a transposed matrix. More...
template int	celeritas::determinant (SquareMatrix< int, 3 > const &)
template float	celeritas::determinant (SquareMatrix< float, 3 > const &)
template double	celeritas::determinant (SquareMatrix< double, 3 > const &)
template int	celeritas::trace (SquareMatrix< int, 3 > const &)
template float	celeritas::trace (SquareMatrix< float, 3 > const &)
template double	celeritas::trace (SquareMatrix< double, 3 > const &)
template void	celeritas::orthonormalize (SquareMatrix< float, 3 > *)
template void	celeritas::orthonormalize (SquareMatrix< double, 3 > *)
template SquareMatrix< int, 3 >	celeritas::gemm (SquareMatrix< int, 3 > const &, SquareMatrix< int, 3 > const &)
template SquareMatrix< float, 3 >	celeritas::gemm (SquareMatrix< float, 3 > const &, SquareMatrix< float, 3 > const &)
template SquareMatrix< double, 3 >	celeritas::gemm (SquareMatrix< double, 3 > const &, SquareMatrix< double, 3 > const &)
template SquareMatrix< int, 3 >	celeritas::gemm (matrix::TransposePolicy, SquareMatrix< int, 3 > const &, SquareMatrix< int, 3 > const &)
template SquareMatrix< float, 3 >	celeritas::gemm (matrix::TransposePolicy, SquareMatrix< float, 3 > const &, SquareMatrix< float, 3 > const &)
template SquareMatrix< double, 3 >	celeritas::gemm (matrix::TransposePolicy, SquareMatrix< double, 3 > const &, SquareMatrix< double, 3 > const &)
template SquareMatrix< real_type, 4 >	celeritas::gemm (SquareMatrix< real_type, 4 > const &, SquareMatrix< real_type, 4 > const &)
template SquareMatrix< real_type, 4 >	celeritas::gemm (matrix::TransposePolicy, SquareMatrix< real_type, 4 > const &, SquareMatrix< real_type, 4 > const &)

Function Documentation

gemm() [1/2]

template<class T , size_type N>

SquareMatrix< T, N > celeritas::gemm	(	matrix::TransposePolicy	,
		SquareMatrix< T, N > const &	a,
		SquareMatrix< T, N > const &	b
	)

Naive square matrix-matrix multiply with the first matrix transposed.

\[ C \gets A^T * B \]

Note

The first argument is a "tag" that alters the behavior of this function versus the non-transposed one.

gemm() [2/2]

template<class T , size_type N>

SquareMatrix< T, N > celeritas::gemm	(	SquareMatrix< T, N > const &	a,
		SquareMatrix< T, N > const &	b
	)

Naive square matrix-matrix multiply.

\[ C \gets A * B \]

This should be equivalent to BLAS' GEMM without the option to transpose, use strides, or multiply by constants. All matrix orderings are C-style: mat[i][j] is for row i, column j .

Note that this uses celeritas::fma which supports types other than floating point.

Warning

This implementation is limited and slow.

make_rotation() [1/2]

SquareMatrixReal3 celeritas::make_rotation	(	Axis	ax,
		Turn	theta,
		Mat3 const &	other
	)

Rotate a C-ordered rotation matrix.

This applies the new axis + turn as a rotation operator to the left of the matrix.

For example, to rotate first by 135 degrees about the Z axis, then 90 degrees about the X axis:

auto r = make_rotation(Axes::x, Turn{0.25}, make_rotation(Axes::z, 0.375));

make_rotation() [2/2]

SquareMatrixReal3 celeritas::make_rotation	(	Real3 const &	ax,
		Turn	theta
	)

Create a C-ordered rotation matrix from an arbitrary rotation.

This is equation (38) in "Rotation Matrices in Two, Three, and Many Dimensions", Physics 116A, UC Santa Cruz, http://scipp.ucsc.edu/~haber/ph116A/.

Parameters

ax	Axis of rotation (unit vector)
theta	Rotation

make_transpose()

SquareMatrixReal3 celeritas::make_transpose

(

SquareMatrixReal3 const &

mat

)

Construct a transposed matrix.

This should only be used for preprocessing. Prefer methods that transpose on the fly.

orthonormalize()

template<class T , size_type N>

void celeritas::orthonormalize

(

SquareMatrix< T, N > *

mat

)

Normalize and orthogonalize a small, dense matrix.

This is used for constructing rotation matrices from user-given matrices that may only have a few digits of precision (e.g. were read from an XML file). It uses the modified Gram-Schmidt orthogonalization algorithm.

If debug assertions are enabled, the normality of the resulting matrix will be checked. A singular matrix will fail.

trace()

template<class T >

T celeritas::trace

(

SquareMatrix< T, 3 > const &

mat

)

Calculate the trace of a 3x3 matrix.

The trace is just the sum of the diagonal elements.