cpu/algebra/columnwise_operations.cpp

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 // Detray test include(s)
0010 #include "detray/test/cpu/algebra_fixture.hpp"
0011
0012 // GTest include(s)
0013 #include <gtest/gtest.h>
0014
0015 namespace detray::test {
0016
0017 TEST_F(detray_algebra, column_wise_cross) {
0018   using matrix3_type = matrix<3, 3>;
0019
0020   auto P = detray::matrix::zero<matrix3_type>();
0021
0022   getter::element(P, 0u, 0u) = 0.f;
0023   getter::element(P, 0u, 1u) = 1.f;
0024   getter::element(P, 0u, 2u) = 2.f;
0025   getter::element(P, 1u, 0u) = 3.f;
0026   getter::element(P, 1u, 1u) = 4.f;
0027   getter::element(P, 1u, 2u) = 5.f;
0028   getter::element(P, 2u, 0u) = 6.f;
0029   getter::element(P, 2u, 1u) = 7.f;
0030   getter::element(P, 2u, 2u) = 8.f;
0031
0032   const vector3 u{1.f, 2.f, 3.f};
0033
0034   const matrix3_type Q = detray::matrix::column_wise_cross(P, u);
0035
0036   EXPECT_NEAR(getter::element(Q, 0u, 0u), -3.f, this->tolerance());
0037   EXPECT_NEAR(getter::element(Q, 1u, 0u), 6.f, this->tolerance());
0038   EXPECT_NEAR(getter::element(Q, 2u, 0u), -3.f, this->tolerance());
0039   EXPECT_NEAR(getter::element(Q, 0u, 1u), -2.f, this->tolerance());
0040   EXPECT_NEAR(getter::element(Q, 1u, 1u), 4.f, this->tolerance());
0041   EXPECT_NEAR(getter::element(Q, 2u, 1u), -2.f, this->tolerance());
0042   EXPECT_NEAR(getter::element(Q, 0u, 2u), -1.f, this->tolerance());
0043   EXPECT_NEAR(getter::element(Q, 1u, 2u), 2.f, this->tolerance());
0044   EXPECT_NEAR(getter::element(Q, 2u, 2u), -1.f, this->tolerance());
0045 }
0046
0047 TEST_F(detray_algebra, column_wise_multiply) {
0048   using matrix3_type = matrix<3, 3>;
0049
0050   auto P = detray::matrix::zero<matrix3_type>();
0051
0052   getter::element(P, 0u, 0u) = 0.f;
0053   getter::element(P, 0u, 1u) = 1.f;
0054   getter::element(P, 0u, 2u) = 2.f;
0055   getter::element(P, 1u, 0u) = 3.f;
0056   getter::element(P, 1u, 1u) = 4.f;
0057   getter::element(P, 1u, 2u) = 5.f;
0058   getter::element(P, 2u, 0u) = 6.f;
0059   getter::element(P, 2u, 1u) = 7.f;
0060   getter::element(P, 2u, 2u) = 8.f;
0061
0062   const vector3 u{1.f, 2.f, 3.f};
0063
0064   const matrix3_type Q = detray::matrix::column_wise_multiply(P, u);
0065
0066   EXPECT_NEAR(getter::element(Q, 0u, 0u), 0.f, this->tolerance());
0067   EXPECT_NEAR(getter::element(Q, 0u, 1u), 1.f, this->tolerance());
0068   EXPECT_NEAR(getter::element(Q, 0u, 2u), 2.f, this->tolerance());
0069   EXPECT_NEAR(getter::element(Q, 1u, 0u), 6.f, this->tolerance());
0070   EXPECT_NEAR(getter::element(Q, 1u, 1u), 8.f, this->tolerance());
0071   EXPECT_NEAR(getter::element(Q, 1u, 2u), 10.f, this->tolerance());
0072   EXPECT_NEAR(getter::element(Q, 2u, 0u), 18.f, this->tolerance());
0073   EXPECT_NEAR(getter::element(Q, 2u, 1u), 21.f, this->tolerance());
0074   EXPECT_NEAR(getter::element(Q, 2u, 2u), 24.f, this->tolerance());
0075 }
0076
0077 TEST_F(detray_algebra, cross_matrix) {
0078   using matrix3_type = matrix<3, 3>;
0079
0080   const vector3 u{1.f, 2.f, 3.f};
0081   const vector3 v{3.f, 4.f, 5.f};
0082
0083   const matrix3_type u_cross = detray::matrix::cross_matrix(u);
0084   const matrix3_type v_cross = detray::matrix::cross_matrix(v);
0085
0086   EXPECT_NEAR(getter::element(u_cross, 0u, 0u), 0.f, this->tolerance());
0087   EXPECT_NEAR(getter::element(u_cross, 0u, 1u), -3.f, this->tolerance());
0088   EXPECT_NEAR(getter::element(u_cross, 0u, 2u), 2.f, this->tolerance());
0089   EXPECT_NEAR(getter::element(u_cross, 1u, 0u), 3.f, this->tolerance());
0090   EXPECT_NEAR(getter::element(u_cross, 1u, 1u), 0.f, this->tolerance());
0091   EXPECT_NEAR(getter::element(u_cross, 1u, 2u), -1.f, this->tolerance());
0092   EXPECT_NEAR(getter::element(u_cross, 2u, 0u), -2.f, this->tolerance());
0093   EXPECT_NEAR(getter::element(u_cross, 2u, 1u), 1.f, this->tolerance());
0094   EXPECT_NEAR(getter::element(u_cross, 2u, 2u), 0.f, this->tolerance());
0095
0096   // [u]_cross * v = [v]_cross^T * u
0097   const vector3 u_cross_v = u_cross * v;
0098   const vector3 v_cross_u = detray::matrix::transpose(v_cross) * u;
0099
0100   EXPECT_NEAR(getter::element(u_cross_v, 0u), -2.f, this->tolerance());
0101   EXPECT_NEAR(getter::element(u_cross_v, 1u), 4.f, this->tolerance());
0102   EXPECT_NEAR(getter::element(u_cross_v, 2u), -2.f, this->tolerance());
0103   EXPECT_NEAR(getter::element(u_cross_v, 0u), getter::element(v_cross_u, 0u),
0104               this->tolerance());
0105   EXPECT_NEAR(getter::element(u_cross_v, 1u), getter::element(v_cross_u, 1u),
0106               this->tolerance());
0107   EXPECT_NEAR(getter::element(u_cross_v, 2u), getter::element(v_cross_u, 2u),
0108               this->tolerance());
0109 }
0110
0111 TEST_F(detray_algebra, outer_product) {
0112   using matrix3_type = matrix<3, 3>;
0113
0114   const vector3 u{1.f, 2.f, 3.f};
0115   const vector3 v{3.f, 4.f, 5.f};
0116
0117   const matrix3_type m33 = detray::matrix::outer_product(u, v);
0118
0119   EXPECT_NEAR(getter::element(m33, 0u, 0u), 3.f, this->tolerance());
0120   EXPECT_NEAR(getter::element(m33, 0u, 1u), 4.f, this->tolerance());
0121   EXPECT_NEAR(getter::element(m33, 0u, 2u), 5.f, this->tolerance());
0122   EXPECT_NEAR(getter::element(m33, 1u, 0u), 6.f, this->tolerance());
0123   EXPECT_NEAR(getter::element(m33, 1u, 1u), 8.f, this->tolerance());
0124   EXPECT_NEAR(getter::element(m33, 1u, 2u), 10.f, this->tolerance());
0125   EXPECT_NEAR(getter::element(m33, 2u, 0u), 9.f, this->tolerance());
0126   EXPECT_NEAR(getter::element(m33, 2u, 1u), 12.f, this->tolerance());
0127   EXPECT_NEAR(getter::element(m33, 2u, 2u), 15.f, this->tolerance());
0128 }
0129
0130 TEST_F(detray_algebra, cholesky_decomposition) {
0131   using matrix3_type = matrix<3, 3>;
0132
0133   // Define A
0134   auto A = detray::matrix::zero<matrix3_type>();
0135   getter::element(A, 0u, 0u) = 4.f;
0136   getter::element(A, 0u, 1u) = 12.f;
0137   getter::element(A, 0u, 2u) = -16.f;
0138   getter::element(A, 1u, 0u) = 12.f;
0139   getter::element(A, 1u, 1u) = 37.f;
0140   getter::element(A, 1u, 2u) = -43.f;
0141   getter::element(A, 2u, 0u) = -16.f;
0142   getter::element(A, 2u, 1u) = -43.f;
0143   getter::element(A, 2u, 2u) = 98.f;
0144
0145   // Get L that satisfies A = L * L^T and check if it is the expected value
0146   const matrix3_type L = detray::matrix::cholesky_decomposition(A);
0147
0148   EXPECT_FLOAT_EQ(static_cast<float>(getter::element(L, 0u, 0u)), 2.f);
0149   EXPECT_FLOAT_EQ(static_cast<float>(getter::element(L, 0u, 1u)), 0.f);
0150   EXPECT_FLOAT_EQ(static_cast<float>(getter::element(L, 0u, 2u)), 0.f);
0151   EXPECT_FLOAT_EQ(static_cast<float>(getter::element(L, 1u, 0u)), 6.f);
0152   EXPECT_FLOAT_EQ(static_cast<float>(getter::element(L, 1u, 1u)), 1.f);
0153   EXPECT_FLOAT_EQ(static_cast<float>(getter::element(L, 1u, 2u)), 0.f);
0154   EXPECT_FLOAT_EQ(static_cast<float>(getter::element(L, 2u, 0u)), -8.f);
0155   EXPECT_FLOAT_EQ(static_cast<float>(getter::element(L, 2u, 1u)), 5.f);
0156   EXPECT_FLOAT_EQ(static_cast<float>(getter::element(L, 2u, 2u)), 3.f);
0157
0158   // Compare A and L * L^T
0159   const matrix3_type B = L * detray::matrix::transpose(L);
0160
0161   for (unsigned int i = 0u; i < 3u; i++) {
0162     for (unsigned int j = 0u; j < 3u; j++) {
0163       EXPECT_FLOAT_EQ(static_cast<float>(getter::element(A, i, j)),
0164                       static_cast<float>(getter::element(B, i, j)));
0165     }
0166   }
0167 }
0168
0169 }  // namespace detray::test