corecel/math/SpanUtils.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 corecel/math/SpanUtils.hh
0006 //! \brief Math functions using celeritas::Span
0007 //! \sa SpanUtils.test.cc
0008 //---------------------------------------------------------------------------//
0009 #pragma once
0010
0011 #include <cstddef>
0012
0013 #include "corecel/cont/Span.hh"
0014
0015 #include "ArrayUtils.hh"
0016
0017 #include "detail/SpanUtilsImpl.hh"
0018
0019 namespace celeritas
0020 {
0021 //! \cond (CELERITAS_DOC_DEV)
0022 //---------------------------------------------------------------------------//
0023 /*!
0024  * Copy all elements from src to dst.
0025  *
0026  * \pre src and dst must not overlap.
0027  *
0028  * \internal The current implementation does not take advantage of the lack of
0029  * overlap.
0030  */
0031 template<class T, std::size_t N>
0032 inline CELER_FUNCTION void copy(Span<T const, N> src, Span<T, N> dst)
0033 {
0034     for (std::size_t i = 0; i != N; ++i)
0035         dst[i] = src[i];
0036 }
0037
0038 //---------------------------------------------------------------------------//
0039 /*!
0040  * Fill all elements of dst with a constant value.
0041  */
0042 template<class T, std::size_t N>
0043 inline CELER_FUNCTION void fill(T value, Span<T, N> dst)
0044 {
0045     for (std::size_t i = 0; i != N; ++i)
0046         dst[i] = value;
0047 }
0048
0049 //---------------------------------------------------------------------------//
0050 /*!
0051  * Increment a vector span by another vector span multiplied by a scalar.
0052  *
0053  * \f[
0054  * y \gets \alpha x + y
0055  * \f]
0056  *
0057  * \pre x and y must not overlap.
0058  */
0059 template<class T, std::size_t N>
0060 inline CELER_FUNCTION void axpy(T a, Span<T const, N> x, Span<T, N> y)
0061 {
0062     Array<T, N> xa = detail::load_array(x);
0063     Array<T, N> ya = detail::load_array(Span<T const, N>{y.data(), N});
0064     axpy(a, xa, &ya);
0065     detail::store_array(ya, y);
0066 }
0067
0068 //---------------------------------------------------------------------------//
0069 /*!
0070  * Dot product of two vector spans.
0071  */
0072 template<class T, std::size_t N>
0073 [[nodiscard]] inline CELER_FUNCTION T dot_product(Span<T const, N> x,
0074                                                   Span<T const, N> y)
0075 {
0076     return dot_product(detail::load_array(x), detail::load_array(y));
0077 }
0078
0079 //---------------------------------------------------------------------------//
0080 /*!
0081  * Cross product of two space vector spans, written into dst.
0082  *
0083  * \pre x, y, and dst must not overlap.
0084  */
0085 template<class T>
0086 inline CELER_FUNCTION void cross_product(
0087     Span<T const, 3> x, Span<T const, 3> y, Span<T, 3> dst)
0088 {
0089     detail::store_array(
0090         cross_product(detail::load_array(x), detail::load_array(y)), dst);
0091 }
0092
0093 //---------------------------------------------------------------------------//
0094 /*!
0095  * Calculate the Euclidean (2) norm of a vector span.
0096  */
0097 template<class T, std::size_t N>
0098 [[nodiscard]] inline CELER_FUNCTION T norm(Span<T const, N> v)
0099 {
0100     return norm(detail::load_array(v));
0101 }
0102
0103 //---------------------------------------------------------------------------//
0104 /*!
0105  * Write a unit-magnitude version of v into dst.
0106  *
0107  * \pre v and dst must not overlap.
0108  */
0109 template<class T, std::size_t N>
0110 inline CELER_FUNCTION void make_unit_vector(Span<T const, N> v, Span<T, N> dst)
0111 {
0112     detail::store_array(make_unit_vector(detail::load_array(v)), dst);
0113 }
0114
0115 //---------------------------------------------------------------------------//
0116 /*!
0117  * Write the component of x orthogonal to the unit vector y into dst.
0118  *
0119  * In this implementation, y must be normalized, and the result is not
0120  * normalized.
0121  *
0122  * \f[
0123    \mathbf{x}' \gets \mathbf{x} - (\mathbf{x} \cdot \mathbf{y}) \mathbf{y}
0124    \, , \quad \|\mathbf{y}\| = 1
0125  * \f]
0126  *
0127  * \pre x, y, and dst must not overlap.
0128  */
0129 template<class T, std::size_t N>
0130 inline CELER_FUNCTION void make_orthogonal(
0131     Span<T const, N> x, Span<T const, N> y, Span<T, N> dst)
0132 {
0133     detail::store_array(
0134         make_orthogonal(detail::load_array(x), detail::load_array(y)), dst);
0135 }
0136
0137 //---------------------------------------------------------------------------//
0138 /*!
0139  * Check whether two unit vector spans are approximately orthogonal.
0140  */
0141 template<class T, std::size_t N>
0142 inline CELER_FUNCTION bool is_soft_orthogonal(Span<T const, N> x,
0143                                               Span<T const, N> y)
0144 {
0145     return is_soft_orthogonal(detail::load_array(x), detail::load_array(y));
0146 }
0147
0148 //---------------------------------------------------------------------------//
0149 /*!
0150  * Check whether two unit vector spans are approximately collinear.
0151  *
0152  * \pre Spans must represent normalized vectors.
0153  */
0154 template<class T, std::size_t N>
0155 inline CELER_FUNCTION bool is_soft_collinear(Span<T const, N> x,
0156                                              Span<T const, N> y)
0157 {
0158     return is_soft_collinear(detail::load_array(x), detail::load_array(y));
0159 }
0160
0161 //---------------------------------------------------------------------------//
0162 /*!
0163  * Calculate the Euclidean (2) distance between two point spans.
0164  */
0165 template<class T, std::size_t N>
0166 [[nodiscard]] inline CELER_FUNCTION T distance(Span<T const, N> x,
0167                                                Span<T const, N> y)
0168 {
0169     return distance(detail::load_array(x), detail::load_array(y));
0170 }
0171
0172 //---------------------------------------------------------------------------//
0173 /*!
0174  * Write a Cartesian unit vector from spherical coordinates into dst.
0175  *
0176  * Theta is the angle between the z axis and the outgoing vector, and phi is
0177  * the angle between the x axis and the projection onto the x-y plane.
0178  */
0179 template<class T>
0180 inline CELER_FUNCTION void from_spherical(T costheta, T phi, Span<T, 3> dst)
0181 {
0182     detail::store_array(from_spherical(costheta, phi), dst);
0183 }
0184
0185 //---------------------------------------------------------------------------//
0186 /*!
0187  * Rotate a direction span by the rotation defined by rot, writing to dst.
0188  *
0189  * \pre dir, rot, and dst must not overlap.
0190  * \sa celeritas::rotate in ArrayUtils.hh for the full description.
0191  */
0192 template<class T>
0193 inline CELER_FUNCTION void rotate(
0194     Span<T const, 3> dir, Span<T const, 3> rot, Span<T, 3> dst)
0195 {
0196     detail::store_array(
0197         rotate(detail::load_array(dir), detail::load_array(rot)), dst);
0198 }
0199
0200 //---------------------------------------------------------------------------//
0201 //! \endcond
0202
0203 }  // namespace celeritas