ACTS
Experiment-independent tracking
Loading...
Searching...
No Matches
hard_coded.hpp
1// This file is part of the ACTS project.
2//
3// Copyright (C) 2016 CERN for the benefit of the ACTS project
4//
5// This Source Code Form is subject to the terms of the Mozilla Public
6// License, v. 2.0. If a copy of the MPL was not distributed with this
7// file, You can obtain one at https://mozilla.org/MPL/2.0/.
8
9#pragma once
10
11// Project include(s).
12#include "detray/algebra/concepts.hpp"
13#include "detray/algebra/generic/algorithms/matrix/determinant/hard_coded.hpp"
14#include "detray/algebra/type_traits.hpp"
15#include "detray/definitions/detail/qualifiers.hpp"
16
17// System include(s)
18#include <type_traits>
19
20namespace detray::algebra::generic::matrix::inverse {
21
23template <concepts::square_matrix matrix_t>
24struct hard_coded {
26 using scalar_t = detray::traits::value_t<matrix_t>;
28 using index_t = detray::traits::index_t<matrix_t>;
29
31 using element_getter_t = detray::traits::element_getter_t<matrix_t>;
33 using determinant_getter = determinant::hard_coded<matrix_t>;
34
38 template <typename M = matrix_t>
39 requires(detray::traits::max_rank<M> == 2)
40 DETRAY_HOST_DEVICE constexpr matrix_t operator()(const matrix_t &m) const {
41 constexpr element_getter_t elem{};
42
43 const scalar_t det = determinant_getter()(m);
44
45 matrix_t ret{};
46 elem(ret, 0, 0) = elem(m, 1, 1) / det;
47 elem(ret, 0, 1) = -elem(m, 0, 1) / det;
48 elem(ret, 1, 0) = -elem(m, 1, 0) / det;
49 elem(ret, 1, 1) = elem(m, 0, 0) / det;
50
51 return ret;
52 }
53
57 template <typename M = matrix_t>
58 requires(detray::traits::max_rank<M> == 4)
59 DETRAY_HOST_DEVICE constexpr matrix_t operator()(const matrix_t &m) const {
60 constexpr element_getter_t elem{};
61
62 matrix_t ret{};
63 elem(ret, 0, 0) = elem(m, 1, 2) * elem(m, 2, 3) * elem(m, 3, 1) -
64 elem(m, 1, 3) * elem(m, 2, 2) * elem(m, 3, 1) +
65 elem(m, 1, 3) * elem(m, 2, 1) * elem(m, 3, 2) -
66 elem(m, 1, 1) * elem(m, 2, 3) * elem(m, 3, 2) -
67 elem(m, 1, 2) * elem(m, 2, 1) * elem(m, 3, 3) +
68 elem(m, 1, 1) * elem(m, 2, 2) * elem(m, 3, 3);
69 elem(ret, 0, 1) = elem(m, 0, 3) * elem(m, 2, 2) * elem(m, 3, 1) -
70 elem(m, 0, 2) * elem(m, 2, 3) * elem(m, 3, 1) -
71 elem(m, 0, 3) * elem(m, 2, 1) * elem(m, 3, 2) +
72 elem(m, 0, 1) * elem(m, 2, 3) * elem(m, 3, 2) +
73 elem(m, 0, 2) * elem(m, 2, 1) * elem(m, 3, 3) -
74 elem(m, 0, 1) * elem(m, 2, 2) * elem(m, 3, 3);
75 elem(ret, 0, 2) = elem(m, 0, 2) * elem(m, 1, 3) * elem(m, 3, 1) -
76 elem(m, 0, 3) * elem(m, 1, 2) * elem(m, 3, 1) +
77 elem(m, 0, 3) * elem(m, 1, 1) * elem(m, 3, 2) -
78 elem(m, 0, 1) * elem(m, 1, 3) * elem(m, 3, 2) -
79 elem(m, 0, 2) * elem(m, 1, 1) * elem(m, 3, 3) +
80 elem(m, 0, 1) * elem(m, 1, 2) * elem(m, 3, 3);
81 elem(ret, 0, 3) = elem(m, 0, 3) * elem(m, 1, 2) * elem(m, 2, 1) -
82 elem(m, 0, 2) * elem(m, 1, 3) * elem(m, 2, 1) -
83 elem(m, 0, 3) * elem(m, 1, 1) * elem(m, 2, 2) +
84 elem(m, 0, 1) * elem(m, 1, 3) * elem(m, 2, 2) +
85 elem(m, 0, 2) * elem(m, 1, 1) * elem(m, 2, 3) -
86 elem(m, 0, 1) * elem(m, 1, 2) * elem(m, 2, 3);
87 elem(ret, 1, 0) = elem(m, 1, 3) * elem(m, 2, 2) * elem(m, 3, 0) -
88 elem(m, 1, 2) * elem(m, 2, 3) * elem(m, 3, 0) -
89 elem(m, 1, 3) * elem(m, 2, 0) * elem(m, 3, 2) +
90 elem(m, 1, 0) * elem(m, 2, 3) * elem(m, 3, 2) +
91 elem(m, 1, 2) * elem(m, 2, 0) * elem(m, 3, 3) -
92 elem(m, 1, 0) * elem(m, 2, 2) * elem(m, 3, 3);
93 elem(ret, 1, 1) = elem(m, 0, 2) * elem(m, 2, 3) * elem(m, 3, 0) -
94 elem(m, 0, 3) * elem(m, 2, 2) * elem(m, 3, 0) +
95 elem(m, 0, 3) * elem(m, 2, 0) * elem(m, 3, 2) -
96 elem(m, 0, 0) * elem(m, 2, 3) * elem(m, 3, 2) -
97 elem(m, 0, 2) * elem(m, 2, 0) * elem(m, 3, 3) +
98 elem(m, 0, 0) * elem(m, 2, 2) * elem(m, 3, 3);
99 elem(ret, 1, 2) = elem(m, 0, 3) * elem(m, 1, 2) * elem(m, 3, 0) -
100 elem(m, 0, 2) * elem(m, 1, 3) * elem(m, 3, 0) -
101 elem(m, 0, 3) * elem(m, 1, 0) * elem(m, 3, 2) +
102 elem(m, 0, 0) * elem(m, 1, 3) * elem(m, 3, 2) +
103 elem(m, 0, 2) * elem(m, 1, 0) * elem(m, 3, 3) -
104 elem(m, 0, 0) * elem(m, 1, 2) * elem(m, 3, 3);
105 elem(ret, 1, 3) = elem(m, 0, 2) * elem(m, 1, 3) * elem(m, 2, 0) -
106 elem(m, 0, 3) * elem(m, 1, 2) * elem(m, 2, 0) +
107 elem(m, 0, 3) * elem(m, 1, 0) * elem(m, 2, 2) -
108 elem(m, 0, 0) * elem(m, 1, 3) * elem(m, 2, 2) -
109 elem(m, 0, 2) * elem(m, 1, 0) * elem(m, 2, 3) +
110 elem(m, 0, 0) * elem(m, 1, 2) * elem(m, 2, 3);
111 elem(ret, 2, 0) = elem(m, 1, 1) * elem(m, 2, 3) * elem(m, 3, 0) -
112 elem(m, 1, 3) * elem(m, 2, 1) * elem(m, 3, 0) +
113 elem(m, 1, 3) * elem(m, 2, 0) * elem(m, 3, 1) -
114 elem(m, 1, 0) * elem(m, 2, 3) * elem(m, 3, 1) -
115 elem(m, 1, 1) * elem(m, 2, 0) * elem(m, 3, 3) +
116 elem(m, 1, 0) * elem(m, 2, 1) * elem(m, 3, 3);
117 elem(ret, 2, 1) = elem(m, 0, 3) * elem(m, 2, 1) * elem(m, 3, 0) -
118 elem(m, 0, 1) * elem(m, 2, 3) * elem(m, 3, 0) -
119 elem(m, 0, 3) * elem(m, 2, 0) * elem(m, 3, 1) +
120 elem(m, 0, 0) * elem(m, 2, 3) * elem(m, 3, 1) +
121 elem(m, 0, 1) * elem(m, 2, 0) * elem(m, 3, 3) -
122 elem(m, 0, 0) * elem(m, 2, 1) * elem(m, 3, 3);
123 elem(ret, 2, 2) = elem(m, 0, 1) * elem(m, 1, 3) * elem(m, 3, 0) -
124 elem(m, 0, 3) * elem(m, 1, 1) * elem(m, 3, 0) +
125 elem(m, 0, 3) * elem(m, 1, 0) * elem(m, 3, 1) -
126 elem(m, 0, 0) * elem(m, 1, 3) * elem(m, 3, 1) -
127 elem(m, 0, 1) * elem(m, 1, 0) * elem(m, 3, 3) +
128 elem(m, 0, 0) * elem(m, 1, 1) * elem(m, 3, 3);
129 elem(ret, 2, 3) = elem(m, 0, 3) * elem(m, 1, 1) * elem(m, 2, 0) -
130 elem(m, 0, 1) * elem(m, 1, 3) * elem(m, 2, 0) -
131 elem(m, 0, 3) * elem(m, 1, 0) * elem(m, 2, 1) +
132 elem(m, 0, 0) * elem(m, 1, 3) * elem(m, 2, 1) +
133 elem(m, 0, 1) * elem(m, 1, 0) * elem(m, 2, 3) -
134 elem(m, 0, 0) * elem(m, 1, 1) * elem(m, 2, 3);
135 elem(ret, 3, 0) = elem(m, 1, 2) * elem(m, 2, 1) * elem(m, 3, 0) -
136 elem(m, 1, 1) * elem(m, 2, 2) * elem(m, 3, 0) -
137 elem(m, 1, 2) * elem(m, 2, 0) * elem(m, 3, 1) +
138 elem(m, 1, 0) * elem(m, 2, 2) * elem(m, 3, 1) +
139 elem(m, 1, 1) * elem(m, 2, 0) * elem(m, 3, 2) -
140 elem(m, 1, 0) * elem(m, 2, 1) * elem(m, 3, 2);
141 elem(ret, 3, 1) = elem(m, 0, 1) * elem(m, 2, 2) * elem(m, 3, 0) -
142 elem(m, 0, 2) * elem(m, 2, 1) * elem(m, 3, 0) +
143 elem(m, 0, 2) * elem(m, 2, 0) * elem(m, 3, 1) -
144 elem(m, 0, 0) * elem(m, 2, 2) * elem(m, 3, 1) -
145 elem(m, 0, 1) * elem(m, 2, 0) * elem(m, 3, 2) +
146 elem(m, 0, 0) * elem(m, 2, 1) * elem(m, 3, 2);
147 elem(ret, 3, 2) = elem(m, 0, 2) * elem(m, 1, 1) * elem(m, 3, 0) -
148 elem(m, 0, 1) * elem(m, 1, 2) * elem(m, 3, 0) -
149 elem(m, 0, 2) * elem(m, 1, 0) * elem(m, 3, 1) +
150 elem(m, 0, 0) * elem(m, 1, 2) * elem(m, 3, 1) +
151 elem(m, 0, 1) * elem(m, 1, 0) * elem(m, 3, 2) -
152 elem(m, 0, 0) * elem(m, 1, 1) * elem(m, 3, 2);
153 elem(ret, 3, 3) = elem(m, 0, 1) * elem(m, 1, 2) * elem(m, 2, 0) -
154 elem(m, 0, 2) * elem(m, 1, 1) * elem(m, 2, 0) +
155 elem(m, 0, 2) * elem(m, 1, 0) * elem(m, 2, 1) -
156 elem(m, 0, 0) * elem(m, 1, 2) * elem(m, 2, 1) -
157 elem(m, 0, 1) * elem(m, 1, 0) * elem(m, 2, 2) +
158 elem(m, 0, 0) * elem(m, 1, 1) * elem(m, 2, 2);
159
160 const scalar_t idet = static_cast<scalar_t>(1.) / determinant_getter()(ret);
161 for (unsigned int c = 0; c < 4; ++c) {
162 for (unsigned int r = 0; r < 4; ++r) {
163 elem(ret, c, r) *= idet;
164 }
165 }
166
167 return ret;
168 }
169};
170
171} // namespace detray::algebra::generic::matrix::inverse
detray::algebra::generic::matrix::determinant::hard_coded
"Determinant getter", assuming a N X N matrix
Definition hard_coded.hpp:23
detray::algebra::generic::matrix::inverse::hard_coded
"inverse getter", assuming a N X N matrix
Definition hard_coded.hpp:24
detray::algebra::generic::matrix::inverse::hard_coded::index_t
detray::traits::index_t< matrix_t > index_t
Matrix index type.
Definition hard_coded.hpp:28
detray::algebra::generic::matrix::inverse::hard_coded::scalar_t
detray::traits::value_t< matrix_t > scalar_t
Matrix element type.
Definition hard_coded.hpp:26
detray::algebra::generic::matrix::inverse::hard_coded::element_getter_t
detray::traits::element_getter_t< matrix_t > element_getter_t
Function (object) used for accessing a matrix element.
Definition hard_coded.hpp:31
detray::algebra::generic::matrix::inverse::hard_coded::determinant_getter
determinant::hard_coded< matrix_t > determinant_getter
Direct determinant implementation used to normalize the 4 by 4 inverse.
Definition hard_coded.hpp:33