include/orange/MatrixUtils.hh

0001 //------------------------------- -*- C++ -*- -------------------------------//
0002 // Copyright Celeritas contributors: see top-level COPYRIGHT file for details
0003 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0004 //---------------------------------------------------------------------------//
0005 //! \file orange/MatrixUtils.hh
0006 //! \todo Split into BLAS and host-only utils
0007 //! \todo Investigate Eigen for setup-time computations if necessary
0008 //---------------------------------------------------------------------------//
0009 #pragma once
0010
0011 #include <cmath>
0012
0013 #include "corecel/Macros.hh"
0014 #include "corecel/Types.hh"
0015 #include "corecel/cont/Array.hh"
0016 #include "corecel/math/Algorithms.hh"
0017 #include "corecel/math/Turn.hh"
0018 #include "geocel/Types.hh"
0019
0020 namespace celeritas
0021 {
0022 //! Policy tags for matrix operations
0023 namespace matrix
0024 {
0025 //---------------------------------------------------------------------------//
0026 struct TransposePolicy
0027 {
0028 };
0029 //! Indicate that the input matrix is transposed
0030 inline constexpr TransposePolicy transpose{};
0031 }  // namespace matrix
0032
0033 //---------------------------------------------------------------------------//
0034 // Apply a matrix to an array
0035 template<class T, size_type N>
0036 inline CELER_FUNCTION Array<T, N> gemv(T alpha,
0037                                        SquareMatrix<T, N> const& a,
0038                                        Array<T, N> const& x,
0039                                        T beta,
0040                                        Array<T, N> const& y);
0041
0042 //---------------------------------------------------------------------------//
0043 // Apply the transpose of a matrix to an array
0044 template<class T, size_type N>
0045 inline CELER_FUNCTION Array<T, N> gemv(matrix::TransposePolicy,
0046                                        T alpha,
0047                                        SquareMatrix<T, N> const& a,
0048                                        Array<T, N> const& x,
0049                                        T beta,
0050                                        Array<T, N> const& y);
0051
0052 //---------------------------------------------------------------------------//
0053 //!@{
0054 //! Apply a matrix or its transpose to an array, without scaling or addition
0055 template<class T, size_type N>
0056 inline CELER_FUNCTION Array<T, N>
0057 gemv(SquareMatrix<T, N> const& a, Array<T, N> const& x)
0058 {
0059     return gemv(T{1}, a, x, T{0}, x);
0060 }
0061
0062 template<class T, size_type N>
0063 inline CELER_FUNCTION Array<T, N>
0064 gemv(matrix::TransposePolicy, SquareMatrix<T, N> const& a, Array<T, N> const& x)
0065 {
0066     return gemv(matrix::transpose, T{1}, a, x, T{0}, x);
0067 }
0068 //!@}
0069 //---------------------------------------------------------------------------//
0070 // Host-only declarations
0071 // (double and float (and some int) for N=3 are instantiated in MatrixUtils.cc)
0072 //---------------------------------------------------------------------------//
0073
0074 // Calculate the determinant of a matrix
0075 template<class T>
0076 T determinant(SquareMatrix<T, 3> const& mat);
0077
0078 // Calculate the trace of a matrix
0079 template<class T>
0080 T trace(SquareMatrix<T, 3> const& mat);
0081
0082 // Perform a matrix-matrix multiply
0083 template<class T, size_type N>
0084 SquareMatrix<T, N>
0085 gemm(SquareMatrix<T, N> const& a, SquareMatrix<T, N> const& b);
0086
0087 // Perform a matrix-matrix multiply with A transposed
0088 template<class T, size_type N>
0089 SquareMatrix<T, N> gemm(matrix::TransposePolicy,
0090                         SquareMatrix<T, N> const& a,
0091                         SquareMatrix<T, N> const& b);
0092
0093 // Normalize and orthogonalize a small, dense matrix
0094 template<class T, size_type N>
0095 void orthonormalize(SquareMatrix<T, N>* mat);
0096
0097 // Create an identity 3x3 matrix
0098 SquareMatrixReal3 make_identity();
0099
0100 // Create a C-ordered rotation matrix about an arbitrary axis
0101 SquareMatrixReal3 make_rotation(Real3 const& ax, Turn rev);
0102
0103 // Create a C-ordered rotation matrix about a Cartesian axis
0104 SquareMatrixReal3 make_rotation(Axis ax, Turn rev);
0105
0106 // Apply a rotation to an existing C-ordered rotation matrix
0107 SquareMatrixReal3 make_rotation(Axis ax, Turn rev, SquareMatrixReal3 const&);
0108
0109 // Scale uniformly
0110 SquareMatrixReal3 make_scaling(real_type scale);
0111
0112 // Scale along an axis
0113 SquareMatrixReal3 make_scaling(Axis ax, real_type scale);
0114
0115 // Scale along all three Cartesian axes
0116 SquareMatrixReal3 make_scaling(Real3 const& scale);
0117
0118 // Reflect across a plane perpendicular to the axis
0119 SquareMatrixReal3 make_reflection(Axis ax);
0120
0121 // Construct a transposed matrix
0122 SquareMatrixReal3 make_transpose(SquareMatrixReal3 const&);
0123
0124 //---------------------------------------------------------------------------//
0125 // INLINE DEFINITIONS
0126 //---------------------------------------------------------------------------//
0127 /*!
0128  * Naive generalized matrix-vector multiply.
0129  *
0130  * \f[
0131  * z \gets \alpha A x + \beta y
0132  * \f]
0133  *
0134  * This should be equivalent to BLAS' GEMV without transposition. All
0135  * matrix orderings are C-style: mat[i][j] is for row i, column j .
0136  *
0137  * \warning This implementation is limited and slow.
0138  */
0139 template<class T, size_type N>
0140 CELER_FUNCTION Array<T, N> gemv(T alpha,
0141                                 SquareMatrix<T, N> const& a,
0142                                 Array<T, N> const& x,
0143                                 T beta,
0144                                 Array<T, N> const& y)
0145 {
0146     Array<T, N> result;
0147     for (size_type i = 0; i != N; ++i)
0148     {
0149         result[i] = beta * y[i];
0150         for (size_type j = 0; j != N; ++j)
0151         {
0152             result[i] = fma(alpha, a[i][j] * x[j], result[i]);
0153         }
0154     }
0155     return result;
0156 }
0157
0158 //---------------------------------------------------------------------------//
0159 /*!
0160  * Naive transposed generalized matrix-vector multiply.
0161  *
0162  * \f[
0163  * z \gets \alpha A^T x + \beta y
0164  * \f]
0165  *
0166  * This should be equivalent to BLAS' GEMV with the 't' option. All
0167  * matrix orderings are C-style: mat[i][j] is for row i, column j .
0168  *
0169  * \warning This implementation is limited and slow.
0170  */
0171 template<class T, size_type N>
0172 CELER_FUNCTION Array<T, N> gemv(matrix::TransposePolicy,
0173                                 T alpha,
0174                                 SquareMatrix<T, N> const& a,
0175                                 Array<T, N> const& x,
0176                                 T beta,
0177                                 Array<T, N> const& y)
0178 {
0179     Array<T, N> result;
0180     for (size_type i = 0; i != N; ++i)
0181     {
0182         result[i] = beta * y[i];
0183     }
0184     for (size_type j = 0; j != N; ++j)
0185     {
0186         for (size_type i = 0; i != N; ++i)
0187         {
0188             result[i] = fma(alpha, a[j][i] * x[j], result[i]);
0189         }
0190     }
0191     return result;
0192 }
0193
0194 //---------------------------------------------------------------------------//
0195 }  // namespace celeritas