Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-05-27 07:24:19

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008 
0009 // Test include(s).
0010 #include "detray/test/cpu/algebra_fixture.hpp"
0011 
0012 // Project include(s)
0013 #include "detray/algebra/utils/approximately_equal.hpp"
0014 #include "detray/algebra/utils/casts.hpp"
0015 #include "detray/algebra/utils/print.hpp"
0016 
0017 // GoogleTest include(s).
0018 #include <gtest/gtest.h>
0019 
0020 // System include(s).
0021 #include <array>
0022 #include <cmath>
0023 #include <cstddef>
0024 #include <limits>
0025 
0026 namespace detray::test {
0027 
0028 #if !DETRAY_ALGEBRA_VC_AOS
0029 TEST_F(detray_algebra, matrix_2x3) {
0030   static constexpr dindex_type<algebra_t> ROWS = 2;
0031   static constexpr dindex_type<algebra_t> COLS = 3;
0032 
0033   using matrix_2x3_t = dmatrix<algebra_t, 2, 3>;
0034 
0035   // Test type traits
0036   static_assert(std::is_same_v<detray::traits::index_t<matrix_2x3_t>,
0037                                dindex_type<algebra_t>>);
0038   static_assert(std::is_same_v<detray::traits::value_t<matrix_2x3_t>,
0039                                dscalar<algebra_t>>);
0040   static_assert(std::is_same_v<detray::traits::scalar_t<matrix_2x3_t>,
0041                                dscalar<algebra_t>>);
0042   static_assert(std::is_same_v<detray::traits::vector_t<matrix_2x3_t>,
0043                                detray::get_vector2D_t<algebra_t>>);
0044 
0045   static_assert(detray::traits::rows<matrix_2x3_t> == 2);
0046   static_assert(detray::traits::columns<matrix_2x3_t> == 3);
0047   static_assert(detray::traits::max_rank<matrix_2x3_t> == 2);
0048   static_assert(detray::traits::size<matrix_2x3_t> == 6);
0049   static_assert(!detray::traits::is_square<matrix_2x3_t>);
0050   static_assert(detray::traits::is_square<dmatrix<algebra_t, 2, 2>>);
0051   static_assert(detray::traits::is_square<dmatrix<algebra_t, 3, 3>>);
0052 
0053   // Test concepts
0054   static_assert(detray::concepts::matrix<matrix_2x3_t>);
0055   static_assert(!detray::concepts::scalar<matrix_2x3_t>);
0056   static_assert(!detray::concepts::vector<matrix_2x3_t>);
0057   static_assert(!detray::concepts::square_matrix<matrix_2x3_t>);
0058 
0059   static_assert(detray::concepts::index<detray::traits::index_t<matrix_2x3_t>>);
0060   static_assert(detray::concepts::value<detray::traits::value_t<matrix_2x3_t>>);
0061   static_assert(
0062       detray::concepts::scalar<detray::traits::scalar_t<matrix_2x3_t>>);
0063   static_assert(
0064       detray::concepts::vector<detray::traits::vector_t<matrix_2x3_t>>);
0065 
0066   // Test on matrix - vector operations
0067   detray::get_vector3D_t<algebra_t> vE{1.f, 2.f, 3.f};
0068 
0069   matrix_2x3_t m23;
0070 
0071   detray::getter::element(m23, 0, 0) = 1.f;
0072   detray::getter::element(m23, 0, 1) = 2.f;
0073   detray::getter::element(m23, 0, 2) = 3.f;
0074   detray::getter::element(m23, 1, 0) = 4.f;
0075   detray::getter::element(m23, 1, 1) = 5.f;
0076   detray::getter::element(m23, 1, 2) = 6.f;
0077 
0078   // Cast to (different) precision
0079   const auto m23_cast_f = detray::algebra::cast_to<float>(m23);
0080 
0081   for (std::size_t j = 0; j < 3; ++j) {
0082     for (std::size_t i = 0; i < 2; ++i) {
0083       auto elem_i = detray::getter::element(m23_cast_f, i, j);
0084 
0085       static_assert(std::same_as<decltype(elem_i), float>);
0086       ASSERT_FLOAT_EQ(elem_i,
0087                       static_cast<float>(detray::getter::element(m23, i, j)));
0088     }
0089   }
0090 
0091   const auto m23_cast_d = detray::algebra::cast_to<double>(m23);
0092 
0093   for (std::size_t j = 0; j < 3; ++j) {
0094     for (std::size_t i = 0; i < 2; ++i) {
0095       auto elem_i = detray::getter::element(m23_cast_d, i, j);
0096 
0097       static_assert(std::same_as<decltype(elem_i), double>);
0098       ASSERT_DOUBLE_EQ(elem_i,
0099                        static_cast<double>(detray::getter::element(m23, i, j)));
0100     }
0101   }
0102 
0103   const auto m23_cast_i = detray::algebra::cast_to<int>(m23);
0104 
0105   for (std::size_t j = 0; j < 3; ++j) {
0106     for (std::size_t i = 0; i < 2; ++i) {
0107       auto elem_i = detray::getter::element(m23_cast_i, i, j);
0108 
0109       static_assert(std::same_as<decltype(elem_i), int>);
0110       ASSERT_EQ(elem_i, static_cast<int>(detray::getter::element(m23, i, j)));
0111     }
0112   }
0113 
0114   detray::get_vector2D_t<algebra_t> v2 = m23 * vE;
0115 
0116   ASSERT_NEAR(v2[0], 14, this->epsilon());
0117   ASSERT_NEAR(v2[1], 32, this->epsilon());
0118 
0119   this->template matrix_test_ops_any_matrix<ROWS, COLS>();
0120 }
0121 
0122 TEST_F(detray_algebra, matrix_3x1) {
0123   static constexpr dindex_type<algebra_t> ROWS = 3;
0124   static constexpr dindex_type<algebra_t> COLS = 1;
0125 
0126   // Cross product on vector3 and matrix<3,1>
0127   dmatrix<algebra_t, 3, 1> vF;
0128   detray::getter::element(vF, 0, 0) = 5.f;
0129   detray::getter::element(vF, 1, 0) = 6.f;
0130   detray::getter::element(vF, 2, 0) = 13.f;
0131 
0132   // Test printing
0133   std::cout << vF << std::endl;
0134 
0135   // Cast to (different) precision
0136   const auto vF_cast_f = detray::algebra::cast_to<float>(vF);
0137 
0138   for (std::size_t i = 0; i < 3; ++i) {
0139     auto elem_i = detray::getter::element(vF_cast_f, i, 0);
0140 
0141     static_assert(std::same_as<decltype(elem_i), float>);
0142     ASSERT_FLOAT_EQ(elem_i,
0143                     static_cast<float>(detray::getter::element(vF, i, 0)));
0144   }
0145 
0146   const auto vF_cast_d = detray::algebra::cast_to<double>(vF);
0147 
0148   for (std::size_t i = 0; i < 3; ++i) {
0149     auto elem_i = detray::getter::element(vF_cast_d, i, 0);
0150 
0151     static_assert(std::same_as<decltype(elem_i), double>);
0152     ASSERT_DOUBLE_EQ(elem_i,
0153                      static_cast<double>(detray::getter::element(vF, i, 0)));
0154   }
0155 
0156   const auto vF_cast_i = detray::algebra::cast_to<int>(vF);
0157 
0158   for (std::size_t i = 0; i < 3; ++i) {
0159     auto elem_i = detray::getter::element(vF_cast_i, i, 0);
0160 
0161     static_assert(std::same_as<decltype(elem_i), int>);
0162     ASSERT_EQ(elem_i, static_cast<int>(detray::getter::element(vF, i, 0)));
0163   }
0164 
0165   detray::get_vector3D_t<algebra_t> vD{1.f, 1.f, 1.f};
0166   detray::get_vector3D_t<algebra_t> vG = vector::cross(vD, vF);
0167   ASSERT_NEAR(vG[0], 7.f, this->epsilon());
0168   ASSERT_NEAR(vG[1], -8.f, this->epsilon());
0169   ASSERT_NEAR(vG[2], 1.f, this->epsilon());
0170 
0171   // Dot product on vector3 and matrix<3,1>
0172   auto dot = vector::dot(vG, vF);
0173   ASSERT_NEAR(dot, 0.f, this->epsilon());
0174 
0175   this->template matrix_test_ops_any_matrix<ROWS, COLS>();
0176 }
0177 
0178 TEST_F(detray_algebra, matrix_6x4) {
0179   // Create the matrix.
0180   static constexpr dindex_type<algebra_t> ROWS = 6;
0181   static constexpr dindex_type<algebra_t> COLS = 4;
0182   using matrix_6x4_t = dmatrix<algebra_t, ROWS, COLS>;
0183   matrix_6x4_t m;
0184 
0185   // Test type traits
0186   static_assert(std::is_same_v<detray::traits::index_t<matrix_6x4_t>,
0187                                dindex_type<algebra_t>>);
0188   static_assert(std::is_same_v<detray::traits::value_t<matrix_6x4_t>,
0189                                dscalar<algebra_t>>);
0190   static_assert(std::is_same_v<detray::traits::scalar_t<matrix_6x4_t>,
0191                                dscalar<algebra_t>>);
0192 
0193   static_assert(detray::traits::rows<matrix_6x4_t> == 6);
0194   static_assert(detray::traits::columns<matrix_6x4_t> == 4);
0195   static_assert(detray::traits::max_rank<matrix_6x4_t> == 4);
0196   static_assert(detray::traits::size<matrix_6x4_t> == 24);
0197   static_assert(!detray::traits::is_square<matrix_6x4_t>);
0198   static_assert(detray::traits::is_square<dmatrix<algebra_t, 4, 4>>);
0199   static_assert(detray::traits::is_square<dmatrix<algebra_t, 6, 6>>);
0200 
0201   // Test concepts
0202   static_assert(detray::concepts::matrix<matrix_6x4_t>);
0203   static_assert(!detray::concepts::scalar<matrix_6x4_t>);
0204   static_assert(!detray::concepts::vector<matrix_6x4_t>);
0205   static_assert(!detray::concepts::square_matrix<matrix_6x4_t>);
0206 
0207   // Fill it.
0208   for (std::size_t i = 0; i < ROWS; ++i) {
0209     for (std::size_t j = 0; j < COLS; ++j) {
0210       detray::getter::element(m, i, j) =
0211           0.5f * static_cast<dscalar<algebra_t>>(i + j);
0212     }
0213   }
0214 
0215   // Check its content.
0216   const dmatrix<algebra_t, ROWS, COLS>& m_const_ref = m;
0217   for (std::size_t i = 0; i < ROWS; ++i) {
0218     for (std::size_t j = 0; j < COLS; ++j) {
0219       const dscalar<algebra_t> ref =
0220           0.5f * static_cast<dscalar<algebra_t>>(i + j);
0221       ASSERT_NEAR(detray::getter::element(m, i, j), ref, this->epsilon());
0222       ASSERT_NEAR(detray::getter::element(m_const_ref, i, j), ref,
0223                   this->epsilon());
0224     }
0225   }
0226 
0227   // Test set_zero
0228   detray::matrix::set_zero(m);
0229   for (std::size_t i = 0; i < ROWS; ++i) {
0230     for (std::size_t j = 0; j < COLS; ++j) {
0231       ASSERT_NEAR(detray::getter::element(m, i, j), 0., this->epsilon());
0232     }
0233   }
0234 
0235   // Test set_identity
0236   detray::matrix::set_identity(m);
0237   for (std::size_t i = 0; i < ROWS; ++i) {
0238     for (std::size_t j = 0; j < COLS; ++j) {
0239       if (i == j) {
0240         ASSERT_NEAR(detray::getter::element(m, i, j), 1.f, this->epsilon());
0241       } else {
0242         ASSERT_NEAR(detray::getter::element(m, i, j), 0.f, this->epsilon());
0243       }
0244     }
0245   }
0246 
0247   // Test block operations
0248   auto b13 = detray::getter::block<1, 3>(m, 0, 0);
0249   ASSERT_NEAR(detray::getter::element(b13, 0, 0), 1.f, this->epsilon());
0250   ASSERT_NEAR(detray::getter::element(b13, 0, 1), 0.f, this->epsilon());
0251   ASSERT_NEAR(detray::getter::element(b13, 0, 2), 0.f, this->epsilon());
0252 
0253   auto b13_tp = detray::matrix::transpose(b13);
0254   ASSERT_NEAR(detray::getter::element(b13_tp, 0, 0), 1.f, this->epsilon());
0255   ASSERT_NEAR(detray::getter::element(b13_tp, 1, 0), 0.f, this->epsilon());
0256   ASSERT_NEAR(detray::getter::element(b13_tp, 2, 0), 0.f, this->epsilon());
0257 
0258   auto b32 = detray::getter::block<3, 2>(m, 2, 2);
0259   ASSERT_NEAR(detray::getter::element(b32, 0, 0), 1.f, this->epsilon());
0260   ASSERT_NEAR(detray::getter::element(b32, 0, 1), 0.f, this->epsilon());
0261   ASSERT_NEAR(detray::getter::element(b32, 1, 0), 0.f, this->epsilon());
0262   ASSERT_NEAR(detray::getter::element(b32, 1, 1), 1.f, this->epsilon());
0263   ASSERT_NEAR(detray::getter::element(b32, 2, 0), 0.f, this->epsilon());
0264   ASSERT_NEAR(detray::getter::element(b32, 2, 1), 0.f, this->epsilon());
0265 
0266   detray::getter::element(b32, 0, 0) = 4.f;
0267   detray::getter::element(b32, 0, 1) = 3.f;
0268   detray::getter::element(b32, 1, 0) = 12.f;
0269   detray::getter::element(b32, 1, 1) = 13.f;
0270   detray::getter::element(b32, 2, 0) = 5.f;
0271   detray::getter::element(b32, 2, 1) = 6.f;
0272 
0273   detray::getter::set_block(m, b32, 2, 2);
0274   ASSERT_NEAR(detray::getter::element(m, 2, 2), 4.f, this->epsilon());
0275   ASSERT_NEAR(detray::getter::element(m, 2, 3), 3.f, this->epsilon());
0276   ASSERT_NEAR(detray::getter::element(m, 3, 2), 12.f, this->epsilon());
0277   ASSERT_NEAR(detray::getter::element(m, 3, 3), 13.f, this->epsilon());
0278   ASSERT_NEAR(detray::getter::element(m, 4, 2), 5.f, this->epsilon());
0279   ASSERT_NEAR(detray::getter::element(m, 4, 3), 6.f, this->epsilon());
0280 
0281   detray::get_vector3D_t<algebra_t> v = {10.f, 20.f, 30.f};
0282   detray::getter::set_block(m, v, 0, 2);
0283   ASSERT_NEAR(detray::getter::element(m, 0, 2), 10., this->epsilon());
0284   ASSERT_NEAR(detray::getter::element(m, 1, 2), 20., this->epsilon());
0285   ASSERT_NEAR(detray::getter::element(m, 2, 2), 30., this->epsilon());
0286 
0287   // Test printing
0288   std::cout << m << std::endl;
0289 
0290   this->template matrix_test_ops_any_matrix<ROWS, COLS>();
0291 }
0292 
0293 TEST_F(detray_algebra, matrix_3x3) {
0294   static constexpr dindex_type<algebra_t> N = 3;
0295 
0296   {
0297     detray::get_vector3D_t<algebra_t> v = {10.f, 20.f, 30.f};
0298     dmatrix<algebra_t, 3, 3> m33;
0299     detray::getter::element(m33, 0, 0) = 1;
0300     detray::getter::element(m33, 1, 0) = 2;
0301     detray::getter::element(m33, 2, 0) = 3;
0302     detray::getter::element(m33, 0, 1) = 5;
0303     detray::getter::element(m33, 1, 1) = 6;
0304     detray::getter::element(m33, 2, 1) = 7;
0305     detray::getter::element(m33, 0, 2) = 9;
0306     detray::getter::element(m33, 1, 2) = 10;
0307     detray::getter::element(m33, 2, 2) = 11;
0308 
0309     const detray::get_vector3D_t<algebra_t> v2 = m33 * v;
0310     ASSERT_NEAR(v2[0], 380., this->epsilon());
0311     ASSERT_NEAR(v2[1], 440., this->epsilon());
0312     ASSERT_NEAR(v2[2], 500., this->epsilon());
0313   }
0314 
0315   {
0316     dmatrix<algebra_t, 3, 3> m33;
0317     detray::getter::element(m33, 0, 0) = 1.f;
0318     detray::getter::element(m33, 0, 1) = 5.f;
0319     detray::getter::element(m33, 0, 2) = 7.f;
0320     detray::getter::element(m33, 1, 0) = 3.f;
0321     detray::getter::element(m33, 1, 1) = 5.f;
0322     detray::getter::element(m33, 1, 2) = 6.f;
0323     detray::getter::element(m33, 2, 0) = 2.f;
0324     detray::getter::element(m33, 2, 1) = 8.f;
0325     detray::getter::element(m33, 2, 2) = 9.f;
0326 
0327     // Test 3 X 3 matrix determinant
0328     auto m33_det = detray::matrix::determinant(m33);
0329     ASSERT_NEAR(m33_det, 20.f, this->tolerance());
0330 
0331     // Test 3 X 3 matrix inverse
0332     auto m33_inv = detray::matrix::inverse(m33);
0333     ASSERT_NEAR(detray::getter::element(m33_inv, 0, 0), -3.f / 20.f,
0334                 this->tolerance());
0335     ASSERT_NEAR(detray::getter::element(m33_inv, 0, 1), 11.f / 20.f,
0336                 this->tolerance());
0337     ASSERT_NEAR(detray::getter::element(m33_inv, 0, 2), -5.f / 20.f,
0338                 this->tolerance());
0339     ASSERT_NEAR(detray::getter::element(m33_inv, 1, 0), -15.f / 20.f,
0340                 this->tolerance());
0341     ASSERT_NEAR(detray::getter::element(m33_inv, 1, 1), -5.f / 20.f,
0342                 this->tolerance());
0343     ASSERT_NEAR(detray::getter::element(m33_inv, 1, 2), 15.f / 20.f,
0344                 this->tolerance());
0345     ASSERT_NEAR(detray::getter::element(m33_inv, 2, 0), 14.f / 20.f,
0346                 this->tolerance());
0347     ASSERT_NEAR(detray::getter::element(m33_inv, 2, 1), 2.f / 20.f,
0348                 this->tolerance());
0349     ASSERT_NEAR(detray::getter::element(m33_inv, 2, 2), -10.f / 20.f,
0350                 this->tolerance());
0351   }
0352 
0353   this->template matrix_test_ops_square_matrix<N>();
0354 }
0355 
0356 TEST_F(detray_algebra, matrix_2x2) {
0357   static constexpr dindex_type<algebra_t> N = 2;
0358 
0359   dmatrix<algebra_t, 2, 2> m22;
0360   detray::getter::element(m22, 0, 0) = 4.f;
0361   detray::getter::element(m22, 0, 1) = 3.f;
0362   detray::getter::element(m22, 1, 0) = 12.f;
0363   detray::getter::element(m22, 1, 1) = 13.f;
0364 
0365   // Test 2 X 2 matrix determinant
0366   auto m22_det = detray::matrix::determinant(m22);
0367   ASSERT_NEAR(m22_det, 16.f, this->tolerance());
0368 
0369   // Test 2 X 2 matrix inverse
0370   auto m22_inv = detray::matrix::inverse(m22);
0371   ASSERT_NEAR(detray::getter::element(m22_inv, 0, 0), 13.f / 16.f,
0372               this->tolerance());
0373   ASSERT_NEAR(detray::getter::element(m22_inv, 0, 1), -3.f / 16.f,
0374               this->tolerance());
0375   ASSERT_NEAR(detray::getter::element(m22_inv, 1, 0), -12.f / 16.f,
0376               this->tolerance());
0377   ASSERT_NEAR(detray::getter::element(m22_inv, 1, 1), 4.f / 16.f,
0378               this->tolerance());
0379 
0380   this->template matrix_test_ops_square_matrix<N>();
0381 }
0382 
0383 TEST_F(detray_algebra, matrix_5x5) {
0384   // Test 5 X 5 matrix
0385   dmatrix<algebra_t, 5, 5> m55;
0386   detray::getter::element(m55, 0, 0) = 1.f;
0387   detray::getter::element(m55, 0, 1) = 3.f;
0388   detray::getter::element(m55, 0, 2) = -9.f;
0389   detray::getter::element(m55, 0, 3) = -5.f;
0390   detray::getter::element(m55, 0, 4) = -2.f;
0391 
0392   detray::getter::element(m55, 1, 0) = -6.f;
0393   detray::getter::element(m55, 1, 1) = -3.f;
0394   detray::getter::element(m55, 1, 2) = 1.f;
0395   detray::getter::element(m55, 1, 3) = 0.f;
0396   detray::getter::element(m55, 1, 4) = 2.f;
0397 
0398   detray::getter::element(m55, 2, 0) = 12.f;
0399   detray::getter::element(m55, 2, 1) = 7.f;
0400   detray::getter::element(m55, 2, 2) = -9.f;
0401   detray::getter::element(m55, 2, 3) = 11.f;
0402   detray::getter::element(m55, 2, 4) = 2.f;
0403 
0404   detray::getter::element(m55, 3, 0) = -3.f;
0405   detray::getter::element(m55, 3, 1) = 4.f;
0406   detray::getter::element(m55, 3, 2) = 5.f;
0407   detray::getter::element(m55, 3, 3) = -6.f;
0408   detray::getter::element(m55, 3, 4) = 7.f;
0409 
0410   detray::getter::element(m55, 4, 0) = 9.f;
0411   detray::getter::element(m55, 4, 1) = 6.f;
0412   detray::getter::element(m55, 4, 2) = 3.f;
0413   detray::getter::element(m55, 4, 3) = 0.f;
0414   detray::getter::element(m55, 4, 4) = -3.f;
0415 
0416   auto m55_det = detray::matrix::determinant(m55);
0417   ASSERT_NEAR((m55_det - 17334.f) / 17334.f, 0.f, this->tolerance());
0418 
0419   auto m55_inv = detray::matrix::inverse(m55);
0420 
0421   ASSERT_NEAR(detray::getter::element(m55_inv, 0, 0), -2106.f / 17334.f,
0422               this->tolerance());
0423   ASSERT_NEAR(detray::getter::element(m55_inv, 0, 1), -12312.f / 17334.f,
0424               this->tolerance());
0425   ASSERT_NEAR(detray::getter::element(m55_inv, 0, 2), -486.f / 17334.f,
0426               this->tolerance());
0427   ASSERT_NEAR(detray::getter::element(m55_inv, 0, 3), 864.f / 17334.f,
0428               this->tolerance());
0429   ASSERT_NEAR(detray::getter::element(m55_inv, 0, 4), -5112.f / 17334.f,
0430               this->tolerance());
0431 
0432   ASSERT_NEAR(detray::getter::element(m55_inv, 1, 0), 2754.f / 17334.f,
0433               this->tolerance());
0434   ASSERT_NEAR(detray::getter::element(m55_inv, 1, 1), 13878.f / 17334.f,
0435               this->tolerance());
0436   ASSERT_NEAR(detray::getter::element(m55_inv, 1, 2), 1080.f / 17334.f,
0437               this->tolerance());
0438   ASSERT_NEAR(detray::getter::element(m55_inv, 1, 3), -315.f / 17334.f,
0439               this->tolerance());
0440   ASSERT_NEAR(detray::getter::element(m55_inv, 1, 4), 7401.f / 17334.f,
0441               this->tolerance());
0442 
0443   ASSERT_NEAR(detray::getter::element(m55_inv, 2, 0), -918.f / 17334.f,
0444               this->tolerance());
0445   ASSERT_NEAR(detray::getter::element(m55_inv, 2, 1), 1152.f / 17334.f,
0446               this->tolerance());
0447   ASSERT_NEAR(detray::getter::element(m55_inv, 2, 2), -360.f / 17334.f,
0448               this->tolerance());
0449   ASSERT_NEAR(detray::getter::element(m55_inv, 2, 3), 105.f / 17334.f,
0450               this->tolerance());
0451   ASSERT_NEAR(detray::getter::element(m55_inv, 2, 4), 1385.f / 17334.f,
0452               this->tolerance());
0453 
0454   ASSERT_NEAR(detray::getter::element(m55_inv, 3, 0), 108.f / 17334.f,
0455               this->tolerance());
0456   ASSERT_NEAR(detray::getter::element(m55_inv, 3, 1), 7002.f / 17334.f,
0457               this->tolerance());
0458   ASSERT_NEAR(detray::getter::element(m55_inv, 3, 2), 1062.f / 17334.f,
0459               this->tolerance());
0460   ASSERT_NEAR(detray::getter::element(m55_inv, 3, 3), -1032.f / 17334.f,
0461               this->tolerance());
0462   ASSERT_NEAR(detray::getter::element(m55_inv, 3, 4), 2896.f / 17334.f,
0463               this->tolerance());
0464 
0465   ASSERT_NEAR(detray::getter::element(m55_inv, 4, 0), -1728.f / 17334.f,
0466               this->tolerance());
0467   ASSERT_NEAR(detray::getter::element(m55_inv, 4, 1), -8028.f / 17334.f,
0468               this->tolerance());
0469   ASSERT_NEAR(detray::getter::element(m55_inv, 4, 2), 342.f / 17334.f,
0470               this->tolerance());
0471   ASSERT_NEAR(detray::getter::element(m55_inv, 4, 3), 2067.f / 17334.f,
0472               this->tolerance());
0473   ASSERT_NEAR(detray::getter::element(m55_inv, 4, 4), -4927.f / 17334.f,
0474               this->tolerance());
0475 }
0476 
0477 TEST_F(detray_algebra, matrix_6x6) {
0478   static constexpr dindex_type<algebra_t> N = 6;
0479 
0480   // Test 6 X 6 big matrix determinant
0481   dmatrix<algebra_t, 6, 6> m66_big;
0482   detray::getter::element(m66_big, 0, 0) = 1.f;
0483   detray::getter::element(m66_big, 0, 1) = 0.f;
0484   detray::getter::element(m66_big, 0, 2) = 3.f;
0485   detray::getter::element(m66_big, 0, 3) = 0.f;
0486   detray::getter::element(m66_big, 0, 4) = 0.f;
0487   detray::getter::element(m66_big, 0, 5) = 0.f;
0488 
0489   detray::getter::element(m66_big, 1, 0) = 0.f;
0490   detray::getter::element(m66_big, 1, 1) = -2.f;
0491   detray::getter::element(m66_big, 1, 2) = 4.f;
0492   detray::getter::element(m66_big, 1, 3) = 0.f;
0493   detray::getter::element(m66_big, 1, 4) = 5.f;
0494   detray::getter::element(m66_big, 1, 5) = 0.f;
0495 
0496   detray::getter::element(m66_big, 2, 0) = 0.f;
0497   detray::getter::element(m66_big, 2, 1) = 0.f;
0498   detray::getter::element(m66_big, 2, 2) = 3.f;
0499   detray::getter::element(m66_big, 2, 3) = 0.f;
0500   detray::getter::element(m66_big, 2, 4) = 0.f;
0501   detray::getter::element(m66_big, 2, 5) = 0.f;
0502 
0503   detray::getter::element(m66_big, 3, 0) = 0.f;
0504   detray::getter::element(m66_big, 3, 1) = 0.f;
0505   detray::getter::element(m66_big, 3, 2) = 0.f;
0506   detray::getter::element(m66_big, 3, 3) = 4.f;
0507   detray::getter::element(m66_big, 3, 4) = 0.f;
0508   detray::getter::element(m66_big, 3, 5) = 0.f;
0509 
0510   detray::getter::element(m66_big, 4, 0) = 0.f;
0511   detray::getter::element(m66_big, 4, 1) = 0.f;
0512   detray::getter::element(m66_big, 4, 2) = 0.f;
0513   detray::getter::element(m66_big, 4, 3) = 0.f;
0514   detray::getter::element(m66_big, 4, 4) = 9.f;
0515   detray::getter::element(m66_big, 4, 5) = 0.f;
0516 
0517   detray::getter::element(m66_big, 5, 0) = -1.f;
0518   detray::getter::element(m66_big, 5, 1) = -1.f;
0519   detray::getter::element(m66_big, 5, 2) = -1.f;
0520   detray::getter::element(m66_big, 5, 3) = -1.f;
0521   detray::getter::element(m66_big, 5, 4) = -1.f;
0522   detray::getter::element(m66_big, 5, 5) = -1.f;
0523 
0524   auto m66_big_det = detray::matrix::determinant(m66_big);
0525   ASSERT_NEAR(m66_big_det, 216.f, 2.f * this->tolerance());
0526 
0527   auto m66_big_inv = detray::matrix::inverse(m66_big);
0528 
0529   // Test comparison
0530   constexpr auto epsilon{std::numeric_limits<dscalar<algebra_t>>::epsilon()};
0531 
0532   ASSERT_TRUE(detray::algebra::approx_equal(m66_big, m66_big));
0533   ASSERT_TRUE(detray::algebra::approx_equal(m66_big, m66_big, epsilon));
0534   ASSERT_TRUE(detray::algebra::approx_equal(m66_big_inv, m66_big_inv));
0535   ASSERT_TRUE(detray::algebra::approx_equal(m66_big_inv, m66_big_inv, epsilon));
0536 
0537   dscalar<algebra_t> rel_err{1.f + 10.f * epsilon};
0538   dmatrix<algebra_t, 6, 6> m66_big_err = rel_err * m66_big;
0539   ASSERT_TRUE(
0540       detray::algebra::approx_equal(m66_big, m66_big_err, 11.f * epsilon));
0541   ASSERT_FALSE(
0542       detray::algebra::approx_equal(m66_big, m66_big_err, 9.f * epsilon));
0543 
0544   rel_err = 1.f + 17.f * epsilon;
0545   m66_big_err = rel_err * m66_big;
0546   ASSERT_TRUE(
0547       detray::algebra::approx_equal(m66_big, m66_big_err, 18.f * epsilon));
0548   ASSERT_FALSE(
0549       detray::algebra::approx_equal(m66_big, m66_big_err, 16.f * epsilon));
0550 
0551   dmatrix<algebra_t, 6, 6> prod66 = m66_big_inv * m66_big;
0552   auto I66 = detray::matrix::identity<decltype(m66_big)>();
0553   // Set the max error for Fastor and SMatrix plugins
0554   ASSERT_TRUE(
0555       detray::algebra::approx_equal(prod66, I66, epsilon, 2.f * epsilon));
0556 
0557   ASSERT_NEAR(detray::getter::element(m66_big_inv, 0, 0), 36.f / 36.f,
0558               this->tolerance());
0559   ASSERT_NEAR(detray::getter::element(m66_big_inv, 0, 1), 0.f / 36.f,
0560               this->tolerance());
0561   ASSERT_NEAR(detray::getter::element(m66_big_inv, 0, 2), -36.f / 36.f,
0562               this->tolerance());
0563   ASSERT_NEAR(detray::getter::element(m66_big_inv, 0, 3), 0.f / 36.f,
0564               this->tolerance());
0565   ASSERT_NEAR(detray::getter::element(m66_big_inv, 0, 4), 0.f / 36.f,
0566               this->tolerance());
0567   ASSERT_NEAR(detray::getter::element(m66_big_inv, 0, 5), 0.f / 36.f,
0568               this->tolerance());
0569 
0570   ASSERT_NEAR(detray::getter::element(m66_big_inv, 1, 0), 0.f / 36.f,
0571               this->tolerance());
0572   ASSERT_NEAR(detray::getter::element(m66_big_inv, 1, 1), -18.f / 36.f,
0573               this->tolerance());
0574   ASSERT_NEAR(detray::getter::element(m66_big_inv, 1, 2), 24.f / 36.f,
0575               this->tolerance());
0576   ASSERT_NEAR(detray::getter::element(m66_big_inv, 1, 3), 0.f / 36.f,
0577               this->tolerance());
0578   ASSERT_NEAR(detray::getter::element(m66_big_inv, 1, 4), 10.f / 36.f,
0579               this->tolerance());
0580   ASSERT_NEAR(detray::getter::element(m66_big_inv, 1, 5), 0.f / 36.f,
0581               this->tolerance());
0582 
0583   ASSERT_NEAR(detray::getter::element(m66_big_inv, 2, 0), 0.f / 36.f,
0584               this->tolerance());
0585   ASSERT_NEAR(detray::getter::element(m66_big_inv, 2, 1), 0.f / 36.f,
0586               this->tolerance());
0587   ASSERT_NEAR(detray::getter::element(m66_big_inv, 2, 2), 12.f / 36.f,
0588               this->tolerance());
0589   ASSERT_NEAR(detray::getter::element(m66_big_inv, 2, 3), 0.f / 36.f,
0590               this->tolerance());
0591   ASSERT_NEAR(detray::getter::element(m66_big_inv, 2, 4), 0.f / 36.f,
0592               this->tolerance());
0593   ASSERT_NEAR(detray::getter::element(m66_big_inv, 2, 5), 0.f / 36.f,
0594               this->tolerance());
0595 
0596   ASSERT_NEAR(detray::getter::element(m66_big_inv, 3, 0), 0.f / 36.f,
0597               this->tolerance());
0598   ASSERT_NEAR(detray::getter::element(m66_big_inv, 3, 1), 0.f / 36.f,
0599               this->tolerance());
0600   ASSERT_NEAR(detray::getter::element(m66_big_inv, 3, 2), 0.f / 36.f,
0601               this->tolerance());
0602   ASSERT_NEAR(detray::getter::element(m66_big_inv, 3, 3), 9.f / 36.f,
0603               this->tolerance());
0604   ASSERT_NEAR(detray::getter::element(m66_big_inv, 3, 4), 0.f / 36.f,
0605               this->tolerance());
0606   ASSERT_NEAR(detray::getter::element(m66_big_inv, 3, 5), 0.f / 36.f,
0607               this->tolerance());
0608 
0609   ASSERT_NEAR(detray::getter::element(m66_big_inv, 4, 0), 0.f / 36.f,
0610               this->tolerance());
0611   ASSERT_NEAR(detray::getter::element(m66_big_inv, 4, 1), 0.f / 36.f,
0612               this->tolerance());
0613   ASSERT_NEAR(detray::getter::element(m66_big_inv, 4, 2), 0.f / 36.f,
0614               this->tolerance());
0615   ASSERT_NEAR(detray::getter::element(m66_big_inv, 4, 3), 0.f / 36.f,
0616               this->tolerance());
0617   ASSERT_NEAR(detray::getter::element(m66_big_inv, 4, 4), 4.f / 36.f,
0618               this->tolerance());
0619   ASSERT_NEAR(detray::getter::element(m66_big_inv, 4, 5), 0.f / 36.f,
0620               this->tolerance());
0621 
0622   ASSERT_NEAR(detray::getter::element(m66_big_inv, 5, 0), -36.f / 36.f,
0623               this->tolerance());
0624   ASSERT_NEAR(detray::getter::element(m66_big_inv, 5, 1), 18.f / 36.f,
0625               this->tolerance());
0626   ASSERT_NEAR(detray::getter::element(m66_big_inv, 5, 2), 0.f / 36.f,
0627               this->tolerance());
0628   ASSERT_NEAR(detray::getter::element(m66_big_inv, 5, 3), -9.f / 36.f,
0629               this->tolerance());
0630   ASSERT_NEAR(detray::getter::element(m66_big_inv, 5, 4), -14.f / 36.f,
0631               this->tolerance());
0632   ASSERT_NEAR(detray::getter::element(m66_big_inv, 5, 5), -36.f / 36.f,
0633               this->tolerance());
0634 
0635   // Cast to (different) precision
0636   const auto m66_cast_f = detray::algebra::cast_to<float>(m66_big_inv);
0637 
0638   for (std::size_t j = 0; j < 6; ++j) {
0639     for (std::size_t i = 0; i < 6; ++i) {
0640       auto elem_i = detray::getter::element(m66_cast_f, i, j);
0641 
0642       static_assert(std::same_as<decltype(elem_i), float>);
0643       ASSERT_FLOAT_EQ(elem_i, static_cast<float>(
0644                                   detray::getter::element(m66_big_inv, i, j)));
0645     }
0646   }
0647 
0648   const auto m66_cast_d = detray::algebra::cast_to<double>(m66_big_inv);
0649 
0650   for (std::size_t j = 0; j < 6; ++j) {
0651     for (std::size_t i = 0; i < 6; ++i) {
0652       auto elem_i = detray::getter::element(m66_cast_d, i, j);
0653 
0654       static_assert(std::same_as<decltype(elem_i), double>);
0655       ASSERT_DOUBLE_EQ(elem_i, static_cast<double>(
0656                                    detray::getter::element(m66_big_inv, i, j)));
0657     }
0658   }
0659 
0660   const auto m66_cast_i = detray::algebra::cast_to<int>(m66_big_inv);
0661 
0662   for (std::size_t j = 0; j < 6; ++j) {
0663     for (std::size_t i = 0; i < 6; ++i) {
0664       auto elem_i = detray::getter::element(m66_cast_i, i, j);
0665 
0666       static_assert(std::same_as<decltype(elem_i), int>);
0667       ASSERT_EQ(elem_i,
0668                 static_cast<int>(detray::getter::element(m66_big_inv, i, j)));
0669     }
0670   }
0671 
0672   // Test 6 X 6 small matrix determinant
0673   dmatrix<algebra_t, 6, 6> m66_small;
0674 
0675   detray::getter::element(m66_small, 0, 0) =
0676       static_cast<dscalar<algebra_t>>(10.792386);
0677   detray::getter::element(m66_small, 0, 1) =
0678       static_cast<dscalar<algebra_t>>(0.216181);
0679   detray::getter::element(m66_small, 0, 2) =
0680       static_cast<dscalar<algebra_t>>(0.057650);
0681   detray::getter::element(m66_small, 0, 3) =
0682       static_cast<dscalar<algebra_t>>(-0.002764);
0683   detray::getter::element(m66_small, 0, 4) =
0684       static_cast<dscalar<algebra_t>>(0.000001);
0685   detray::getter::element(m66_small, 0, 5) = static_cast<dscalar<algebra_t>>(0);
0686 
0687   detray::getter::element(m66_small, 1, 0) =
0688       static_cast<dscalar<algebra_t>>(43.909368);
0689   detray::getter::element(m66_small, 1, 1) =
0690       static_cast<dscalar<algebra_t>>(10.372997);
0691   detray::getter::element(m66_small, 1, 2) =
0692       static_cast<dscalar<algebra_t>>(0.231496);
0693   detray::getter::element(m66_small, 1, 3) =
0694       static_cast<dscalar<algebra_t>>(-0.065972);
0695   detray::getter::element(m66_small, 1, 4) =
0696       static_cast<dscalar<algebra_t>>(-0.000002);
0697   detray::getter::element(m66_small, 1, 5) = static_cast<dscalar<algebra_t>>(0);
0698 
0699   detray::getter::element(m66_small, 2, 0) =
0700       static_cast<dscalar<algebra_t>>(0.045474);
0701   detray::getter::element(m66_small, 2, 1) =
0702       static_cast<dscalar<algebra_t>>(-0.001730);
0703   detray::getter::element(m66_small, 2, 2) =
0704       static_cast<dscalar<algebra_t>>(0.000246);
0705   detray::getter::element(m66_small, 2, 3) =
0706       static_cast<dscalar<algebra_t>>(0.000004);
0707   detray::getter::element(m66_small, 2, 4) = static_cast<dscalar<algebra_t>>(0);
0708   detray::getter::element(m66_small, 2, 5) = static_cast<dscalar<algebra_t>>(0);
0709 
0710   detray::getter::element(m66_small, 3, 0) =
0711       static_cast<dscalar<algebra_t>>(-0.255134);
0712   detray::getter::element(m66_small, 3, 1) =
0713       static_cast<dscalar<algebra_t>>(-0.059846);
0714   detray::getter::element(m66_small, 3, 2) =
0715       static_cast<dscalar<algebra_t>>(-0.001345);
0716   detray::getter::element(m66_small, 3, 3) =
0717       static_cast<dscalar<algebra_t>>(0.000383);
0718   detray::getter::element(m66_small, 3, 4) = static_cast<dscalar<algebra_t>>(0);
0719   detray::getter::element(m66_small, 3, 5) = static_cast<dscalar<algebra_t>>(0);
0720 
0721   detray::getter::element(m66_small, 4, 0) =
0722       static_cast<dscalar<algebra_t>>(-0.001490);
0723   detray::getter::element(m66_small, 4, 1) =
0724       static_cast<dscalar<algebra_t>>(-0.000057);
0725   detray::getter::element(m66_small, 4, 2) =
0726       static_cast<dscalar<algebra_t>>(-0.000008);
0727   detray::getter::element(m66_small, 4, 3) =
0728       static_cast<dscalar<algebra_t>>(0.000001);
0729   detray::getter::element(m66_small, 4, 4) =
0730       static_cast<dscalar<algebra_t>>(0.000001);
0731   detray::getter::element(m66_small, 4, 5) = static_cast<dscalar<algebra_t>>(0);
0732 
0733   detray::getter::element(m66_small, 5, 0) = static_cast<dscalar<algebra_t>>(0);
0734   detray::getter::element(m66_small, 5, 1) = static_cast<dscalar<algebra_t>>(0);
0735   detray::getter::element(m66_small, 5, 2) = static_cast<dscalar<algebra_t>>(0);
0736   detray::getter::element(m66_small, 5, 3) = static_cast<dscalar<algebra_t>>(0);
0737   detray::getter::element(m66_small, 5, 4) = static_cast<dscalar<algebra_t>>(0);
0738   detray::getter::element(m66_small, 5, 5) =
0739       static_cast<dscalar<algebra_t>>(89875.517874);
0740 
0741   auto m66_small_det = detray::matrix::determinant(m66_small);
0742   ASSERT_NEAR((m66_small_det - 4.30636e-11f) / 4.30636e-11f, 0.f,
0743               2.f * this->tolerance());
0744 
0745   this->template matrix_test_ops_square_matrix<N>();
0746 }
0747 
0748 TEST_F(detray_algebra, matrix_7x4_4x12) {
0749   this->template matrix_test_ops_inhomogeneous_multipliable_matrices<7, 4,
0750                                                                      12>();
0751 }
0752 
0753 TEST_F(detray_algebra, matrix_17x9_9x4) {
0754   this->template matrix_test_ops_inhomogeneous_multipliable_matrices<17, 9,
0755                                                                      4>();
0756 }
0757 
0758 TEST_F(detray_algebra, matrix_5x2_2x3) {
0759   this->template matrix_test_ops_inhomogeneous_multipliable_matrices<5, 2, 3>();
0760 }
0761 
0762 TEST_F(detray_algebra, matrix_small_mixed) {
0763   dmatrix<algebra_t, 2, 2> m22;
0764   detray::getter::element(m22, 0, 0) = 4.f;
0765   detray::getter::element(m22, 0, 1) = 3.f;
0766   detray::getter::element(m22, 1, 0) = 12.f;
0767   detray::getter::element(m22, 1, 1) = 13.f;
0768 
0769   // Test 2 X 2 matrix inverse
0770   auto m22_inv = detray::matrix::inverse(m22);
0771   ASSERT_NEAR(detray::getter::element(m22_inv, 0, 0), 13.f / 16.f,
0772               this->tolerance());
0773   ASSERT_NEAR(detray::getter::element(m22_inv, 0, 1), -3.f / 16.f,
0774               this->tolerance());
0775   ASSERT_NEAR(detray::getter::element(m22_inv, 1, 0), -12.f / 16.f,
0776               this->tolerance());
0777   ASSERT_NEAR(detray::getter::element(m22_inv, 1, 1), 4.f / 16.f,
0778               this->tolerance());
0779 
0780   dmatrix<algebra_t, 3, 3> m33;
0781   detray::getter::element(m33, 0, 0) = 1.f;
0782   detray::getter::element(m33, 0, 1) = 5.f;
0783   detray::getter::element(m33, 0, 2) = 7.f;
0784   detray::getter::element(m33, 1, 0) = 3.f;
0785   detray::getter::element(m33, 1, 1) = 5.f;
0786   detray::getter::element(m33, 1, 2) = 6.f;
0787   detray::getter::element(m33, 2, 0) = 2.f;
0788   detray::getter::element(m33, 2, 1) = 8.f;
0789   detray::getter::element(m33, 2, 2) = 9.f;
0790 
0791   auto m33_inv = detray::matrix::inverse(m33);
0792 
0793   // Test Zero
0794   auto m23 = detray::matrix::zero<dmatrix<algebra_t, 2, 3>>();
0795   detray::getter::element(m23, 0, 0) += 2.f;
0796   detray::getter::element(m23, 0, 1) += 3.f;
0797   detray::getter::element(m23, 0, 2) += 4.f;
0798   detray::getter::element(m23, 1, 0) += 5.f;
0799   detray::getter::element(m23, 1, 1) += 6.f;
0800   detray::getter::element(m23, 1, 2) += 7.f;
0801 
0802   // Test scalar X Matrix
0803   m23 = 2. * m23;
0804   ASSERT_NEAR(detray::getter::element(m23, 0, 0), 4.f, this->epsilon());
0805   ASSERT_NEAR(detray::getter::element(m23, 0, 1), 6.f, this->epsilon());
0806   ASSERT_NEAR(detray::getter::element(m23, 0, 2), 8.f, this->epsilon());
0807   ASSERT_NEAR(detray::getter::element(m23, 1, 0), 10.f, this->epsilon());
0808   ASSERT_NEAR(detray::getter::element(m23, 1, 1), 12.f, this->epsilon());
0809   ASSERT_NEAR(detray::getter::element(m23, 1, 2), 14.f, this->epsilon());
0810 
0811   // Test Transpose
0812   auto m32 = detray::matrix::transpose(m23);
0813 
0814   // Test Identity and (Matrix + Matrix)
0815   m32 = m32 + detray::matrix::identity<dmatrix<algebra_t, 3, 2>>();
0816 
0817   // Test Matrix X scalar
0818   m32 = m32 * 2.f;
0819 
0820   ASSERT_NEAR(detray::getter::element(m32, 0, 0), 10.f, this->epsilon());
0821   ASSERT_NEAR(detray::getter::element(m32, 0, 1), 20.f, this->epsilon());
0822   ASSERT_NEAR(detray::getter::element(m32, 1, 0), 12.f, this->epsilon());
0823   ASSERT_NEAR(detray::getter::element(m32, 1, 1), 26.f, this->epsilon());
0824   ASSERT_NEAR(detray::getter::element(m32, 2, 0), 16.f, this->epsilon());
0825   ASSERT_NEAR(detray::getter::element(m32, 2, 1), 28.f, this->epsilon());
0826 
0827   // Test Matrix multiplication
0828   m22 = m22_inv * m23 * m33_inv * m32;
0829 
0830   ASSERT_NEAR(detray::getter::element(m22, 0, 0), 6.225f, this->tolerance());
0831   ASSERT_NEAR(detray::getter::element(m22, 0, 1), 14.675f, this->tolerance());
0832   ASSERT_NEAR(detray::getter::element(m22, 1, 0), -3.3f, this->tolerance());
0833   ASSERT_NEAR(detray::getter::element(m22, 1, 1), -7.9f, this->tolerance());
0834 }
0835 #endif
0836 }  // namespace detray::test