Implement GPU-compatible transformations with ORANGE

src/corecel/cont/Span.hh

@@ -177,7 +177,7 @@ template<class T, std::size_t N>
 constexpr std::size_t Span<T, N>::extent;
 
 //---------------------------------------------------------------------------//
-// HELPER FUNCTIONS
+// FREE FUNCTIONS
 //---------------------------------------------------------------------------//
 //! Get a mutable fixed-size view to an array
 template<class T, std::size_t N>
@@ -218,5 +218,19 @@ CELER_FUNCTION Span<const typename T::value_type> make_span(T const& cont)
     return {cont.data(), cont.size()};
 }
 
+//---------------------------------------------------------------------------//
+//! Construct an array from a fixed-size span
+template<class T, std::size_t N>
+CELER_CONSTEXPR_FUNCTION Array<std::remove_cv_t<T>, N>
+make_array(Span<T, N> const& s)
+{
+    Array<std::remove_cv_t<T>, N> result;
+    for (std::size_t i = 0; i < N; ++i)
+    {
+        result[i] = s[i];
+    }
+    return result;
+}
+
 //---------------------------------------------------------------------------//
 }  // namespace celeritas

src/orange/CMakeLists.txt

@@ -22,6 +22,7 @@ list(APPEND SOURCES
   detail/RectArrayInserter.cc
   detail/UnitInserter.cc
   surf/SurfaceIO.cc
+  transform/Transformation.cc
 )
 
 if(CELERITAS_USE_JSON)

src/orange/MatrixUtils.cc

@@ -7,7 +7,10 @@
 //---------------------------------------------------------------------------//
 #include "MatrixUtils.hh"
 
+#include <cmath>
+
 #include "corecel/math/ArrayUtils.hh"
+#include "corecel/math/SoftEqual.hh"
 
 using Mat3 = celeritas::SquareMatrixReal3;
 
@@ -66,11 +69,12 @@ gemm(SquareMatrix<T, N> const& a, SquareMatrix<T, N> const& b)
 /*!
  * Normalize and orthogonalize a small, dense matrix.
  *
- * \return normalization (zero if already orthonormal)
- *
  * 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.
  */
 template<class T, size_type N>
 void orthonormalize(SquareMatrix<T, N>* mat)
@@ -95,19 +99,8 @@ void orthonormalize(SquareMatrix<T, N>* mat)
         }
     }
 
-    // Check result for orthonormality: if the matrix is singular this will
-    // fail
-    if (CELERITAS_DEBUG)
-    {
-        for (size_type i = 0; i != N; ++i)
-        {
-            for (size_type ip = 0; ip != i; ++ip)
-            {
-                CELER_ENSURE(soft_equal(dot_product((*mat)[i], (*mat)[ip]),
-                                        static_cast<T>(i == ip ? 1 : 0)));
-            }
-        }
-    }
+    // Check result for orthonormality
+    CELER_ENSURE(soft_equal(std::fabs(determinant(*mat)), T{1}));
 }
 
 //---------------------------------------------------------------------------//

src/orange/MatrixUtils.hh

@@ -15,6 +15,18 @@
 
 namespace celeritas
 {
+namespace matrix
+{
+//---------------------------------------------------------------------------//
+//!@{
+//! Policy tags
+struct TransposePolicy
+{
+};
+inline constexpr TransposePolicy transpose{};
+//!@}
+}  // namespace matrix
+
 //---------------------------------------------------------------------------//
 // Apply a matrix to an array
 template<class T, size_type N>
@@ -25,14 +37,32 @@ inline CELER_FUNCTION Array<T, N> gemv(T alpha,
                                        Array<T, N> const& y);
 
 //---------------------------------------------------------------------------//
-//! Apply a matrix to an array without scaling or addition
+// Apply the transpose of a matrix to an array
+template<class T, size_type N>
+inline CELER_FUNCTION Array<T, N> gemv(matrix::TransposePolicy,
+                                       T alpha,
+                                       SquareMatrix<T, N> const& a,
+                                       Array<T, N> const& x,
+                                       T beta,
+                                       Array<T, N> const& y);
+
+//---------------------------------------------------------------------------//
+//!@{
+//! Apply a matrix or its transpose to an array, without scaling or addition
 template<class T, size_type N>
 inline CELER_FUNCTION Array<T, N>
 gemv(SquareMatrix<T, N> const& a, Array<T, N> const& x)
 {
     return gemv(T{1}, a, x, T{0}, x);
 }
 
+template<class T, size_type N>
+inline CELER_FUNCTION Array<T, N>
+gemv(matrix::TransposePolicy, SquareMatrix<T, N> const& a, Array<T, N> const& x)
+{
+    return gemv(matrix::transpose, T{1}, a, x, T{0}, x);
+}
+//!@}
 //---------------------------------------------------------------------------//
 // Host-only declarations
 // (double and float (and some int) for N=3 are instantiated in MatrixUtils.cc)
@@ -67,7 +97,7 @@ SquareMatrixReal3 make_rotation(Axis ax, Turn rev, SquareMatrixReal3 const&);
  * z \gets \alpha A x + \beta y
  * \f]
  *
- * This should be equivalent to BLAS' GEMV without the option to transpose. All
+ * 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.
@@ -91,5 +121,41 @@ CELER_FUNCTION Array<T, N> gemv(T alpha,
     return result;
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * Naive transposed generalized matrix-vector multiply.
+ *
+ * \f[
+ * z \gets \alpha A^T x + \beta y
+ * \f]
+ *
+ * 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.
+ */
+template<class T, size_type N>
+CELER_FUNCTION Array<T, N> gemv(matrix::TransposePolicy,
+                                T alpha,
+                                SquareMatrix<T, N> const& a,
+                                Array<T, N> const& x,
+                                T beta,
+                                Array<T, N> const& y)
+{
+    Array<T, N> result;
+    for (size_type i = 0; i != N; ++i)
+    {
+        result[i] = beta * y[i];
+    }
+    for (size_type j = 0; j != N; ++j)
+    {
+        for (size_type i = 0; i != N; ++i)
+        {
+            result[i] += alpha * (a[j][i] * x[j]);
+        }
+    }
+    return result;
+}
+
 //---------------------------------------------------------------------------//
 }  // namespace celeritas

src/orange/OrangeData.hh

@@ -283,7 +283,7 @@ struct OrangeParamsData
     Items<detail::BIHLeafNode> bih_leaf_nodes;
 
     Items<Daughter> daughters;
-    Items<Translation> translations;
+    Items<Real3> translations;
 
     UniverseIndexerData<W, M> universe_indexer_data;
 

src/orange/OrangeTrackView.hh

@@ -14,9 +14,9 @@
 
 #include "OrangeData.hh"
 #include "OrangeTypes.hh"
-#include "Translator.hh"
 #include "detail/LevelStateAccessor.hh"
 #include "detail/UniverseIndexer.hh"
+#include "transform/Translation.hh"
 #include "univ/SimpleUnitTracker.hh"
 #include "univ/UniverseTypeTraits.hh"
 #include "univ/detail/Types.hh"
@@ -263,9 +263,8 @@ OrangeTrackView::operator=(Initializer_t const& init)
         if (daughter_id)
         {
             auto const& daughter = params_.daughters[daughter_id];
-            auto const& trans = params_.translations[daughter.translation_id];
-            TranslatorDown td(trans);
-            local.pos = td(local.pos);
+            Translation trans{params_.translations[daughter.translation_id]};
+            local.pos = trans.transform_down(local.pos);
 
             uid = daughter.universe_id;
             ++level;
@@ -544,10 +543,8 @@ CELER_FUNCTION void OrangeTrackView::move_internal(Real3 const& pos)
             auto daughter_id = tracker.daughter(lsa.vol());
             CELER_ASSERT(daughter_id);
             auto const& daughter = params_.daughters[daughter_id];
-            auto const& trans = params_.translations[daughter.translation_id];
-
-            TranslatorDown td(trans);
-            local_pos = td(pos);
+            Translation trans{params_.translations[daughter.translation_id]};
+            local_pos = trans.transform_down(pos);
         }
     }
 
@@ -625,16 +622,15 @@ CELER_FUNCTION void OrangeTrackView::cross_boundary()
         auto daughter = params_.daughters[daughter_id];
         // Get the translator at the parent level, in order to translate into
         // daughter
-        TranslatorDown translator(
-            params_.translations[daughter.translation_id]);
+        Translation trans{params_.translations[daughter.translation_id]};
 
         // Make the current level the daughter level
         ++level;
         universe_id = daughter.universe_id;
         auto tracker = this->make_tracker(universe_id);
 
         // Create local state on the daughter level
-        local.pos = translator(local.pos);
+        local.pos = trans.transform_down(local.pos);
         local.volume = {};
         local.surface = {};
         local.temp_sense = this->make_temp_sense();

src/orange/OrangeTypes.hh

@@ -64,11 +64,8 @@ using SimpleUnitId = OpaqueId<struct SimpleUnitRecord>;
 //! Opaque index for rectilinear array data
 using RectArrayId = OpaqueId<struct RectArrayRecord>;
 
-//! Translation of a single embedded universe
-using Translation = Real3;
-
 //! Identifier for a translation of a single embedded universe
-using TranslationId = OpaqueId<Translation>;
+using TranslationId = OpaqueId<Real3>;
 
 //! Identifier for a relocatable set of volumes
 using UniverseId = OpaqueId<struct Universe>;
@@ -124,6 +121,17 @@ enum class SurfaceType : unsigned char
     size_  //!< Sentinel value for number of surface types
 };
 
+//---------------------------------------------------------------------------//
+/*!
+ * Enumeration for mapping transform implementations to integers.
+ */
+enum class TransformType : unsigned char
+{
+    translation,  //!< Translation only
+    transformation,  //!< Translation plus rotation
+    size_
+};
+
 //---------------------------------------------------------------------------//
 /*!
  * Enumeration for type-deleted universe storage.