MatrixUtils.cc File Reference

#include "MatrixUtils.hh"
#include <algorithm>
#include <cmath>
#include "corecel/cont/Range.hh"
#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.
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.
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.
template<class T , size_type N>
void	celeritas::orthonormalize (SquareMatrix< T, N > *mat)
	Normalize and orthogonalize a small, dense matrix.
Mat3	celeritas::make_rotation (Real3 const &ax, Turn theta)
	Create a C-ordered rotation matrix from an arbitrary rotation.
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.
SquareMatrixReal3	celeritas::make_identity ()
	Create an identity matrix.
SquareMatrixReal3	celeritas::make_scaling (real_type scale)
	Create a uniform scaling matrix.
SquareMatrixReal3	celeritas::make_scaling (Axis ax, real_type scale)
	Create a scaling matrix along a given axis.
SquareMatrixReal3	celeritas::make_scaling (Real3 const &scale)
	Create a scaling matrix along all three Cartesian axes.
SquareMatrixReal3	celeritas::make_reflection (Axis ax)
	Create a reflection matrix perpendicular to a given axis.
SquareMatrixReal3	celeritas::make_transpose (SquareMatrixReal3 const &mat)
	Construct a transposed matrix.
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

determinant()

template<class T >

T celeritas::determinant

(

SquareMatrix< T, 3 > const &

mat

)

Calculate the determiniant of a 3x3 matrix.

This loop typically gets unrolled and uses literals instead of a modulo.

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_reflection()

SquareMatrixReal3 celeritas::make_reflection

(

Axis

ax

)

Create a reflection matrix perpendicular to a given axis.

This creates a matrix that reflects across a plane on the origin, normal to the specified axis. The sign of the coordinate on that axis is reversed.

Parameters

ax	Axis to reflect

Returns

Matrix that represents reflection across the given axis

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_scaling() [1/3]

SquareMatrixReal3 celeritas::make_scaling	(	Axis	ax,
		real_type	scale
	)

Create a scaling matrix along a given axis.

This creates a matrix that scales by the given factor along the specified axis while preserving the other dimensions.

Parameters

ax	Axis along which to scale
scale	Scale factor to apply

make_scaling() [2/3]

SquareMatrixReal3 celeritas::make_scaling

(

Real3 const &

scale

)

Create a scaling matrix along all three Cartesian axes.

Parameters

scale

Scale factor for each axis

make_scaling() [3/3]

SquareMatrixReal3 celeritas::make_scaling

(

real_type

scale

)

Create a uniform scaling matrix.

This creates a matrix that scales by the given factor along all axes.

Parameters

scale

Scale factor to apply

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.