MatrixUtils.hh File Reference

#include <cmath>
#include "corecel/Macros.hh"
#include "corecel/Types.hh"
#include "corecel/cont/Array.hh"
#include "corecel/math/Algorithms.hh"
#include "corecel/math/Turn.hh"
#include "geocel/Types.hh"

This graph shows which files directly or indirectly include this file:

Classes
struct	celeritas::matrix::TransposePolicy

Namespaces
namespace	celeritas::matrix
	Policy tags for matrix operations.

Functions
template<class T , size_type N>
CELER_FUNCTION Array< T, N >	celeritas::gemv (T alpha, SquareMatrix< T, N > const &a, Array< T, N > const &x, T beta, Array< T, N > const &y)
	Naive generalized matrix-vector multiply.
template<class T , size_type N>
CELER_FUNCTION Array< T, N >	celeritas::gemv (matrix::TransposePolicy, T alpha, SquareMatrix< T, N > const &a, Array< T, N > const &x, T beta, Array< T, N > const &y)
	Naive transposed generalized matrix-vector multiply.
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.
SquareMatrixReal3	celeritas::make_identity ()
	Create an identity 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_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<class T , size_type N>
CELER_FUNCTION Array< T, N >	celeritas::gemv (SquareMatrix< T, N > const &a, Array< T, N > const &x)
template<class T , size_type N>
CELER_FUNCTION Array< T, N >	celeritas::gemv (matrix::TransposePolicy, SquareMatrix< T, N > const &a, Array< T, N > const &x)

Variables
constexpr TransposePolicy	celeritas::matrix::transpose {}
	Indicate that the input matrix is transposed.

Detailed Description

Todo:

Split into BLAS and host-only utils

Investigate Eigen for setup-time computations if necessary

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.

gemv() [1/4]

template<class T , size_type N>

CELER_FUNCTION Array< T, N > celeritas::gemv ( matrix::TransposePolicy  , SquareMatrix< T, N > const &  a, Array< T, N > const &  x  )

inline

Apply a matrix or its transpose to an array, without scaling or addition

gemv() [2/4]

template<class T , size_type N>

CELER_FUNCTION Array< T, N > celeritas::gemv ( matrix::TransposePolicy  , T  alpha, SquareMatrix< T, N > const &  a, Array< T, N > const &  x, T  beta, Array< T, N > const &  y  )

inline

Naive transposed generalized matrix-vector multiply.

\[ z \gets \alpha A^T x + \beta y \]

This should be equivalent to BLAS' GEMV with the 't' option. All matrix orderings are C-style: mat[i][j] is for row i, column j .

Warning

This implementation is limited and slow.

gemv() [3/4]

template<class T , size_type N>

CELER_FUNCTION Array< T, N > celeritas::gemv ( SquareMatrix< T, N > const &  a, Array< T, N > const &  x  )

inline

Apply a matrix or its transpose to an array, without scaling or addition

gemv() [4/4]

template<class T , size_type N>

CELER_FUNCTION Array< T, N > celeritas::gemv ( T  alpha, SquareMatrix< T, N > const &  a, Array< T, N > const &  x, T  beta, Array< T, N > const &  y  )

inline

Naive generalized matrix-vector multiply.

\[ z \gets \alpha A x + \beta y \]

This should be equivalent to BLAS' GEMV without transposition. All matrix orderings are C-style: mat[i][j] is for row i, column j .

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.