Line data Source code
1 1 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : /// \file
5 : /// Defines function computing the determinant and inverse of a tensor.
6 :
7 : #pragma once
8 :
9 : #include <type_traits>
10 : #include <utility>
11 :
12 : #include "DataStructures/DataVector.hpp"
13 : #include "DataStructures/Tensor/Tensor.hpp"
14 : #include "DataStructures/Variables.hpp"
15 : #include "Utilities/ErrorHandling/Assert.hpp"
16 : #include "Utilities/Gsl.hpp"
17 : #include "Utilities/MakeWithValue.hpp"
18 : #include "Utilities/Requires.hpp"
19 : #include "Utilities/SetNumberOfGridPoints.hpp"
20 : #include "Utilities/TMPL.hpp"
21 :
22 : namespace determinant_and_inverse_detail {
23 : // Helps to shorten some repeated code:
24 : template <typename Index0, typename Index1>
25 : using inverse_indices =
26 : tmpl::list<change_index_up_lo<Index1>, change_index_up_lo<Index0>>;
27 :
28 : template <typename Symm, typename Index0, typename Index1,
29 : typename = std::nullptr_t>
30 : struct DetAndInverseImpl;
31 :
32 : template <typename Symm, typename Index0, typename Index1>
33 : struct DetAndInverseImpl<Symm, Index0, Index1, Requires<Index0::dim == 1>> {
34 : template <typename T>
35 : static void apply(
36 : const gsl::not_null<Scalar<T>*> det,
37 : const gsl::not_null<Tensor<T, Symm, inverse_indices<Index0, Index1>>*>
38 : inv,
39 : const Tensor<T, Symm, tmpl::list<Index0, Index1>>& tensor) {
40 : const T& t00 = get<0, 0>(tensor);
41 : // inv is non-const so that it can be moved into the std::pair:
42 : get(*det) = t00;
43 : get<0, 0>(*inv) = 1.0 / t00;
44 : }
45 : };
46 :
47 : // The inverse of a 2x2 tensor is computed from Cramer's rule.
48 : template <typename Index0, typename Index1>
49 : struct DetAndInverseImpl<Symmetry<2, 1>, Index0, Index1,
50 : Requires<Index0::dim == 2>> {
51 : template <typename T>
52 : static void apply(
53 : const gsl::not_null<Scalar<T>*> det,
54 : const gsl::not_null<
55 : Tensor<T, Symmetry<2, 1>, inverse_indices<Index0, Index1>>*>
56 : inv,
57 : const Tensor<T, Symmetry<2, 1>, tmpl::list<Index0, Index1>>& tensor) {
58 : const T& t00 = get<0, 0>(tensor);
59 : const T& t01 = get<0, 1>(tensor);
60 : const T& t10 = get<1, 0>(tensor);
61 : const T& t11 = get<1, 1>(tensor);
62 :
63 : get(*det) = t00 * t11 - t01 * t10;
64 : const T one_over_det = make_with_value<T>(*det, 1.0) / get(*det);
65 : get<0, 0>(*inv) = t11 * one_over_det;
66 : get<0, 1>(*inv) = -t01 * one_over_det;
67 : get<1, 0>(*inv) = -t10 * one_over_det;
68 : get<1, 1>(*inv) = t00 * one_over_det;
69 : }
70 : };
71 :
72 : template <typename Index0>
73 : struct DetAndInverseImpl<Symmetry<1, 1>, Index0, Index0,
74 : Requires<Index0::dim == 2>> {
75 : template <typename T>
76 : static void apply(
77 : const gsl::not_null<Scalar<T>*> det,
78 : const gsl::not_null<
79 : Tensor<T, Symmetry<1, 1>, inverse_indices<Index0, Index0>>*>
80 : inv,
81 : const Tensor<T, Symmetry<1, 1>, tmpl::list<Index0, Index0>>& tensor) {
82 : const T& t00 = get<0, 0>(tensor);
83 : const T& t01 = get<0, 1>(tensor);
84 : const T& t11 = get<1, 1>(tensor);
85 :
86 : get(*det) = t00 * t11 - t01 * t01;
87 : const T one_over_det = make_with_value<T>(*det, 1.0) / get(*det);
88 : get<0, 0>(*inv) = t11 * one_over_det;
89 : get<0, 1>(*inv) = -t01 * one_over_det;
90 : get<1, 1>(*inv) = t00 * one_over_det;
91 : }
92 : };
93 :
94 : // The inverse of a 3x3 tensor is computed from Cramer's rule. By reusing some
95 : // terms, the determinant is computed efficiently at the same time.
96 : template <typename Index0, typename Index1>
97 : struct DetAndInverseImpl<Symmetry<2, 1>, Index0, Index1,
98 : Requires<Index0::dim == 3>> {
99 : template <typename T>
100 : static void apply(
101 : const gsl::not_null<Scalar<T>*> det,
102 : const gsl::not_null<
103 : Tensor<T, Symmetry<2, 1>, inverse_indices<Index0, Index1>>*>
104 : inv,
105 : const Tensor<T, Symmetry<2, 1>, tmpl::list<Index0, Index1>>& tensor) {
106 : const T& t00 = get<0, 0>(tensor);
107 : const T& t01 = get<0, 1>(tensor);
108 : const T& t02 = get<0, 2>(tensor);
109 : const T& t10 = get<1, 0>(tensor);
110 : const T& t11 = get<1, 1>(tensor);
111 : const T& t12 = get<1, 2>(tensor);
112 : const T& t20 = get<2, 0>(tensor);
113 : const T& t21 = get<2, 1>(tensor);
114 : const T& t22 = get<2, 2>(tensor);
115 : const T a = t11 * t22 - t12 * t21;
116 : const T b = t12 * t20 - t10 * t22;
117 : const T c = t10 * t21 - t11 * t20;
118 :
119 : get(*det) = t00 * a + t01 * b + t02 * c;
120 : const T one_over_det = make_with_value<T>(*det, 1.0) / get(*det);
121 : get<0, 0>(*inv) = a * one_over_det;
122 : get<0, 1>(*inv) = (t21 * t02 - t22 * t01) * one_over_det;
123 : get<0, 2>(*inv) = (t01 * t12 - t02 * t11) * one_over_det;
124 : get<1, 0>(*inv) = b * one_over_det;
125 : get<1, 1>(*inv) = (t22 * t00 - t20 * t02) * one_over_det;
126 : get<1, 2>(*inv) = (t02 * t10 - t00 * t12) * one_over_det;
127 : get<2, 0>(*inv) = c * one_over_det;
128 : get<2, 1>(*inv) = (t20 * t01 - t21 * t00) * one_over_det;
129 : get<2, 2>(*inv) = (t00 * t11 - t01 * t10) * one_over_det;
130 : }
131 : };
132 :
133 : template <typename Index0>
134 : struct DetAndInverseImpl<Symmetry<1, 1>, Index0, Index0,
135 : Requires<Index0::dim == 3>> {
136 : template <typename T>
137 : static void apply(
138 : const gsl::not_null<Scalar<T>*> det,
139 : const gsl::not_null<
140 : Tensor<T, Symmetry<1, 1>, inverse_indices<Index0, Index0>>*>
141 : inv,
142 : const Tensor<T, Symmetry<1, 1>, tmpl::list<Index0, Index0>>& tensor) {
143 : const T& t00 = get<0, 0>(tensor);
144 : const T& t01 = get<0, 1>(tensor);
145 : const T& t02 = get<0, 2>(tensor);
146 : const T& t11 = get<1, 1>(tensor);
147 : const T& t12 = get<1, 2>(tensor);
148 : const T& t22 = get<2, 2>(tensor);
149 : const T a = t11 * t22 - t12 * t12;
150 : const T b = t12 * t02 - t01 * t22;
151 : const T c = t01 * t12 - t11 * t02;
152 :
153 : get(*det) = t00 * a + t01 * b + t02 * c;
154 : const T one_over_det = make_with_value<T>(*det, 1.0) / get(*det);
155 : get<0, 0>(*inv) = (t11 * t22 - t12 * t12) * one_over_det;
156 : get<0, 1>(*inv) = (t12 * t02 - t22 * t01) * one_over_det;
157 : get<0, 2>(*inv) = (t01 * t12 - t02 * t11) * one_over_det;
158 : get<1, 1>(*inv) = (t22 * t00 - t02 * t02) * one_over_det;
159 : get<1, 2>(*inv) = (t02 * t01 - t00 * t12) * one_over_det;
160 : get<2, 2>(*inv) = (t00 * t11 - t01 * t01) * one_over_det;
161 : }
162 : };
163 :
164 : // The 4x4 tensor inverse is implemented using the formula for inverting a
165 : // partitioned matrix (here, partitioned into 2x2 blocks). This is more
166 : // efficient than Cramer's rule for matrices larger than 3x3.
167 : //
168 : // In this algorithm, the 4x4 input matrix is partitioned into the 2x2 blocks:
169 : // P Q
170 : // R S
171 : // The 4x4 inverse matrix is partitioned into the 2x2 blocks:
172 : // U V
173 : // W X
174 : // Each inverse block is obtained from the blocks {P,Q,R,S} of the input matrix
175 : // as follows:
176 : // X = inv[S - R inv(P) Q] (i.e. inv(X) is the Schur complement of P)
177 : // W = - X R inv(P)
178 : // V = - inv(P) Q X
179 : // U = inv(P) + inv(P) Q X R inv(P) = inv(P) - inv(P) Q W
180 : template <typename Index0, typename Index1>
181 : struct DetAndInverseImpl<Symmetry<2, 1>, Index0, Index1,
182 : Requires<Index0::dim == 4>> {
183 : template <typename T>
184 : static void apply(
185 : const gsl::not_null<Scalar<T>*> det,
186 : const gsl::not_null<
187 : Tensor<T, Symmetry<2, 1>, inverse_indices<Index0, Index1>>*>
188 : inv,
189 : const Tensor<T, Symmetry<2, 1>, tmpl::list<Index0, Index1>>& tensor) {
190 : const T& p00 = get<0, 0>(tensor);
191 : const T& p01 = get<0, 1>(tensor);
192 : const T& p10 = get<1, 0>(tensor);
193 : const T& p11 = get<1, 1>(tensor);
194 : const T& q00 = get<0, 2>(tensor);
195 : const T& q01 = get<0, 3>(tensor);
196 : const T& q10 = get<1, 2>(tensor);
197 : const T& q11 = get<1, 3>(tensor);
198 : const T& r00 = get<2, 0>(tensor);
199 : const T& r01 = get<2, 1>(tensor);
200 : const T& r10 = get<3, 0>(tensor);
201 : const T& r11 = get<3, 1>(tensor);
202 : const T& s00 = get<2, 2>(tensor);
203 : const T& s01 = get<2, 3>(tensor);
204 : const T& s10 = get<3, 2>(tensor);
205 : const T& s11 = get<3, 3>(tensor);
206 :
207 : T& u00 = get<0, 0>(*inv);
208 : T& u01 = get<0, 1>(*inv);
209 : T& u10 = get<1, 0>(*inv);
210 : T& u11 = get<1, 1>(*inv);
211 : T& v00 = get<0, 2>(*inv);
212 : T& v01 = get<0, 3>(*inv);
213 : T& v10 = get<1, 2>(*inv);
214 : T& v11 = get<1, 3>(*inv);
215 : T& w00 = get<2, 0>(*inv);
216 : T& w01 = get<2, 1>(*inv);
217 : T& w10 = get<3, 0>(*inv);
218 : T& w11 = get<3, 1>(*inv);
219 : T& x00 = get<2, 2>(*inv);
220 : T& x01 = get<2, 3>(*inv);
221 : T& x10 = get<3, 2>(*inv);
222 : T& x11 = get<3, 3>(*inv);
223 :
224 : // Temporarily store det(P) in det
225 : get(*det) = p00 * p11 - p01 * p10;
226 : const T& det_p = get(*det);
227 : const T one_over_det_p = make_with_value<T>(det_p, 1.0) / det_p;
228 : const T inv_p00 = p11 * one_over_det_p;
229 : const T inv_p01 = -p01 * one_over_det_p;
230 : const T inv_p10 = -p10 * one_over_det_p;
231 : const T inv_p11 = p00 * one_over_det_p;
232 :
233 : const T r_inv_p00 = r00 * inv_p00 + r01 * inv_p10;
234 : const T r_inv_p01 = r00 * inv_p01 + r01 * inv_p11;
235 : const T r_inv_p10 = r10 * inv_p00 + r11 * inv_p10;
236 : const T r_inv_p11 = r10 * inv_p01 + r11 * inv_p11;
237 :
238 : const T inv_p_q00 = inv_p00 * q00 + inv_p01 * q10;
239 : const T inv_p_q01 = inv_p00 * q01 + inv_p01 * q11;
240 : const T inv_p_q10 = inv_p10 * q00 + inv_p11 * q10;
241 : const T inv_p_q11 = inv_p10 * q01 + inv_p11 * q11;
242 :
243 : const T inv_x00 = s00 - (r_inv_p00 * q00 + r_inv_p01 * q10);
244 : const T inv_x01 = s01 - (r_inv_p00 * q01 + r_inv_p01 * q11);
245 : const T inv_x10 = s10 - (r_inv_p10 * q00 + r_inv_p11 * q10);
246 : const T inv_x11 = s11 - (r_inv_p10 * q01 + r_inv_p11 * q11);
247 :
248 : const T det_inv_x = inv_x00 * inv_x11 - inv_x01 * inv_x10;
249 : const T one_over_det_inv_x = make_with_value<T>(det_inv_x, 1.0) / det_inv_x;
250 : x00 = inv_x11 * one_over_det_inv_x;
251 : x01 = -inv_x01 * one_over_det_inv_x;
252 : x10 = -inv_x10 * one_over_det_inv_x;
253 : x11 = inv_x00 * one_over_det_inv_x;
254 :
255 : w00 = -x00 * r_inv_p00 - x01 * r_inv_p10;
256 : w01 = -x00 * r_inv_p01 - x01 * r_inv_p11;
257 : w10 = -x10 * r_inv_p00 - x11 * r_inv_p10;
258 : w11 = -x10 * r_inv_p01 - x11 * r_inv_p11;
259 :
260 : v00 = -inv_p_q00 * x00 - inv_p_q01 * x10;
261 : v01 = -inv_p_q00 * x01 - inv_p_q01 * x11;
262 : v10 = -inv_p_q10 * x00 - inv_p_q11 * x10;
263 : v11 = -inv_p_q10 * x01 - inv_p_q11 * x11;
264 :
265 : u00 = inv_p00 - (inv_p_q00 * w00 + inv_p_q01 * w10);
266 : u01 = inv_p01 - (inv_p_q00 * w01 + inv_p_q01 * w11);
267 : u10 = inv_p10 - (inv_p_q10 * w00 + inv_p_q11 * w10);
268 : u11 = inv_p11 - (inv_p_q10 * w01 + inv_p_q11 * w11);
269 :
270 : get(*det) *= det_inv_x;
271 : }
272 : };
273 :
274 : template <typename Index0>
275 : struct DetAndInverseImpl<Symmetry<1, 1>, Index0, Index0,
276 : Requires<Index0::dim == 4>> {
277 : template <typename T>
278 : static void apply(
279 : const gsl::not_null<Scalar<T>*> det,
280 : const gsl::not_null<
281 : Tensor<T, Symmetry<1, 1>, inverse_indices<Index0, Index0>>*>
282 : inv,
283 : const Tensor<T, Symmetry<1, 1>, tmpl::list<Index0, Index0>>& tensor) {
284 : const T& p00 = get<0, 0>(tensor);
285 : const T& p01 = get<0, 1>(tensor);
286 : const T& p11 = get<1, 1>(tensor);
287 : const T& q00 = get<0, 2>(tensor);
288 : const T& q01 = get<0, 3>(tensor);
289 : const T& q10 = get<1, 2>(tensor);
290 : const T& q11 = get<1, 3>(tensor);
291 : const T& s00 = get<2, 2>(tensor);
292 : const T& s01 = get<2, 3>(tensor);
293 : const T& s11 = get<3, 3>(tensor);
294 :
295 : T& u00 = get<0, 0>(*inv);
296 : T& u01 = get<0, 1>(*inv);
297 : T& u11 = get<1, 1>(*inv);
298 : T& v00 = get<0, 2>(*inv);
299 : T& v01 = get<0, 3>(*inv);
300 : T& v10 = get<1, 2>(*inv);
301 : T& v11 = get<1, 3>(*inv);
302 : T& x00 = get<2, 2>(*inv);
303 : T& x01 = get<2, 3>(*inv);
304 : T& x11 = get<3, 3>(*inv);
305 :
306 : // Temporarily store det(P) in det
307 : get(*det) = p00 * p11 - p01 * p01;
308 : const T& det_p = get(*det);
309 : const T one_over_det_p = make_with_value<T>(det_p, 1.0) / det_p;
310 : const T inv_p00 = p11 * one_over_det_p;
311 : const T inv_p01 = -p01 * one_over_det_p;
312 : const T inv_p11 = p00 * one_over_det_p;
313 :
314 : const T r_inv_p00 = q00 * inv_p00 + q10 * inv_p01;
315 : const T r_inv_p01 = q00 * inv_p01 + q10 * inv_p11;
316 : const T r_inv_p10 = q01 * inv_p00 + q11 * inv_p01;
317 : const T r_inv_p11 = q01 * inv_p01 + q11 * inv_p11;
318 :
319 : const T inv_x00 = s00 - (r_inv_p00 * q00 + r_inv_p01 * q10);
320 : const T inv_x01 = s01 - (r_inv_p00 * q01 + r_inv_p01 * q11);
321 : const T inv_x11 = s11 - (r_inv_p10 * q01 + r_inv_p11 * q11);
322 :
323 : const T det_inv_x = inv_x00 * inv_x11 - inv_x01 * inv_x01;
324 : const T one_over_det_inv_x = make_with_value<T>(det_inv_x, 1.0) / det_inv_x;
325 : x00 = inv_x11 * one_over_det_inv_x;
326 : x01 = -inv_x01 * one_over_det_inv_x;
327 : x11 = inv_x00 * one_over_det_inv_x;
328 :
329 : v00 = -x00 * r_inv_p00 - x01 * r_inv_p10;
330 : v10 = -x00 * r_inv_p01 - x01 * r_inv_p11;
331 : v01 = -x01 * r_inv_p00 - x11 * r_inv_p10;
332 : v11 = -x01 * r_inv_p01 - x11 * r_inv_p11;
333 :
334 : u00 = inv_p00 - (r_inv_p00 * v00 + r_inv_p10 * v01);
335 : u01 = inv_p01 - (r_inv_p00 * v10 + r_inv_p10 * v11);
336 : u11 = inv_p11 - (r_inv_p01 * v10 + r_inv_p11 * v11);
337 :
338 : get(*det) *= det_inv_x;
339 : }
340 : };
341 : } // namespace determinant_and_inverse_detail
342 :
343 : /// @{
344 : /*!
345 : * \ingroup TensorGroup
346 : * \brief Computes the determinant and inverse of a rank-2 Tensor.
347 : *
348 : * Computes the determinant and inverse together, because this leads to
349 : * fewer operations compared to computing the determinant independently.
350 : *
351 : * \details
352 : * Treats the input rank-2 tensor as a matrix. The first (second) index
353 : * of the tensor corresponds to the rows (columns) of the matrix. The
354 : * determinant is a scalar tensor. The inverse is a rank-2 tensor whose
355 : * indices are reversed and of opposite valence relative to the input
356 : * tensor, i.e. given \f$T_a^b\f$ returns \f$(Tinv)_b^a\f$.
357 : *
358 : * \note
359 : * When inverting a 4x4 spacetime metric, it is typically more efficient
360 : * to use the 3+1 decomposition of the 4-metric in terms of lapse,
361 : * shift, and spatial 3-metric, in which only the spatial 3-metric needs
362 : * to be inverted.
363 : */
364 : template <typename T, typename Symm, typename Index0, typename Index1>
365 1 : void determinant_and_inverse(
366 : const gsl::not_null<Scalar<T>*> det,
367 : const gsl::not_null<Tensor<
368 : T, Symm,
369 : tmpl::list<change_index_up_lo<Index1>, change_index_up_lo<Index0>>>*>
370 : inv,
371 : const Tensor<T, Symm, tmpl::list<Index0, Index1>>& tensor) {
372 : static_assert(Index0::dim == Index1::dim,
373 : "Cannot take the inverse of a Tensor whose Indices are not "
374 : "of the same dimensionality.");
375 : static_assert(Index0::index_type == Index1::index_type,
376 : "Taking the inverse of a mixed Spatial and Spacetime index "
377 : "Tensor is not allowed since it's not clear what that means.");
378 : static_assert(not std::is_integral<T>::value, "Can't invert a Tensor<int>.");
379 :
380 : set_number_of_grid_points(det, tensor);
381 : set_number_of_grid_points(inv, tensor);
382 : determinant_and_inverse_detail::DetAndInverseImpl<Symm, Index0,
383 : Index1>::apply(det, inv,
384 : tensor);
385 : }
386 :
387 : template <typename T, typename Symm, typename Index0, typename Index1>
388 1 : auto determinant_and_inverse(
389 : const Tensor<T, Symm, tmpl::list<Index0, Index1>>& tensor)
390 : -> std::pair<Scalar<T>, Tensor<T, Symm,
391 : tmpl::list<change_index_up_lo<Index1>,
392 : change_index_up_lo<Index0>>>> {
393 : std::pair<Scalar<T>, Tensor<T, Symm,
394 : tmpl::list<change_index_up_lo<Index1>,
395 : change_index_up_lo<Index0>>>>
396 : result{};
397 : determinant_and_inverse_detail::DetAndInverseImpl<
398 : Symm, Index0, Index1>::apply(make_not_null(&result.first),
399 : make_not_null(&result.second), tensor);
400 : return result;
401 : }
402 : /// @}
403 :
404 : /// @{
405 : /*!
406 : * \ingroup TensorGroup
407 : * \brief Computes the determinant and inverse of a rank-2 Tensor.
408 : *
409 : * Computes the determinant and inverse together, because this leads to fewer
410 : * operations compared to computing the determinant independently.
411 : *
412 : * \tparam DetTag the Tag for the determinant of input Tensor.
413 : * \tparam InvTag the Tag for the inverse of input Tensor.
414 : *
415 : * \details
416 : * See determinant_and_inverse().
417 : */
418 : template <typename DetTag, typename InvTag, typename T, typename Symm,
419 : typename Index0, typename Index1>
420 1 : void determinant_and_inverse(
421 : const gsl::not_null<Variables<tmpl::list<DetTag, InvTag>>*> det_and_inv,
422 : const Tensor<T, Symm, tmpl::list<Index0, Index1>>& tensor) {
423 : static_assert(std::is_same_v<typename DetTag::type, Scalar<T>>,
424 : "Type of first return tag must correspond to that of input's "
425 : "determinant.");
426 : static_assert(
427 : std::is_same_v<typename InvTag::type,
428 : Tensor<T, Symm,
429 : tmpl::list<change_index_up_lo<Index1>,
430 : change_index_up_lo<Index0>>>>,
431 : "Type of second return tag must correspond to that of input's inverse.");
432 : const auto number_of_grid_points = get<0, 0>(tensor).size();
433 : if (UNLIKELY(number_of_grid_points != det_and_inv->number_of_grid_points())) {
434 : det_and_inv->initialize(number_of_grid_points);
435 : }
436 : determinant_and_inverse_detail::DetAndInverseImpl<
437 : Symm, Index0, Index1>::apply(make_not_null(&get<DetTag>(*det_and_inv)),
438 : make_not_null(&get<InvTag>(*det_and_inv)),
439 : tensor);
440 : }
441 :
442 : template <typename DetTag, typename InvTag, typename T, typename Symm,
443 : typename Index0, typename Index1>
444 1 : auto determinant_and_inverse(
445 : const Tensor<T, Symm, tmpl::list<Index0, Index1>>& tensor)
446 : -> Variables<tmpl::list<DetTag, InvTag>> {
447 : static_assert(std::is_same_v<typename DetTag::type, Scalar<T>>,
448 : "Type of first return tag must correspond to that of input's "
449 : "determinant.");
450 : static_assert(std::is_same_v<typename InvTag::type,
451 : Tensor<T, Symm,
452 : tmpl::list<change_index_up_lo<Index1>,
453 : change_index_up_lo<Index0>>>>,
454 : "Type of second return tag must correspond to that of input's "
455 : "inverse.");
456 : Variables<tmpl::list<DetTag, InvTag>> result(get<0, 0>(tensor).size());
457 : determinant_and_inverse_detail::DetAndInverseImpl<
458 : Symm, Index0, Index1>::apply(make_not_null(&get<DetTag>(result)),
459 : make_not_null(&get<InvTag>(result)), tensor);
460 : return result;
461 : }
462 : /// @}
|
