DeterminantAndInverse.hpp
Go to the documentation of this file.
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"
14 #include "ErrorHandling/Assert.hpp"
16 #include "Utilities/Requires.hpp"
17 
19 // Helps to shorten some repeated code:
20 template <typename Index0, typename Index1>
21 using inverse_indices =
22  tmpl::list<change_index_up_lo<Index1>, change_index_up_lo<Index0>>;
23 template <typename T, typename Symm, typename Index0, typename Index1>
24 using determinant_inverse_pair =
25  std::pair<Scalar<T>, Tensor<T, Symm, inverse_indices<Index0, Index1>>>;
26 
27 template <typename Symm, typename Index0, typename Index1,
28  typename = std::nullptr_t>
29 struct DetAndInverseImpl;
30 
31 template <typename Symm, typename Index0, typename Index1>
32 struct DetAndInverseImpl<Symm, Index0, Index1, Requires<Index0::dim == 1>> {
33  template <typename T>
34  static determinant_inverse_pair<T, Symm, Index0, Index1> apply(
35  const Tensor<T, Symm, tmpl::list<Index0, Index1>>& tensor) noexcept {
36  const T& t00 = get<0, 0>(tensor);
37  // inv is non-const so that it can be moved into the std::pair:
38  Tensor<T, Symm, inverse_indices<Index0, Index1>> inv{
39  make_with_value<T>(t00, 1.0) / t00};
40  return std::make_pair(Scalar<T>{t00}, std::move(inv));
41  }
42 };
43 
44 // The inverse of a 2x2 tensor is computed from Cramer's rule.
45 template <typename Index0, typename Index1>
46 struct DetAndInverseImpl<Symmetry<2, 1>, Index0, Index1,
47  Requires<Index0::dim == 2>> {
48  template <typename T>
49  static determinant_inverse_pair<T, Symmetry<2, 1>, Index0, Index1> apply(
50  const Tensor<T, Symmetry<2, 1>, tmpl::list<Index0, Index1>>&
51  tensor) noexcept {
52  const T& t00 = get<0, 0>(tensor);
53  const T& t01 = get<0, 1>(tensor);
54  const T& t10 = get<1, 0>(tensor);
55  const T& t11 = get<1, 1>(tensor);
56  // det is non-const so that it can be moved into the std::pair:
57  Scalar<T> det{t00 * t11 - t01 * t10};
58  const T one_over_det = make_with_value<T>(det, 1.0) / det.get();
59  Tensor<T, Symmetry<2, 1>, inverse_indices<Index0, Index1>> inv{};
60  get<0, 0>(inv) = t11 * one_over_det;
61  get<0, 1>(inv) = -t01 * one_over_det;
62  get<1, 0>(inv) = -t10 * one_over_det;
63  get<1, 1>(inv) = t00 * one_over_det;
64  return std::make_pair(std::move(det), std::move(inv));
65  }
66 };
67 
68 template <typename Index0>
69 struct DetAndInverseImpl<Symmetry<1, 1>, Index0, Index0,
70  Requires<Index0::dim == 2>> {
71  template <typename T>
72  static determinant_inverse_pair<T, Symmetry<1, 1>, Index0, Index0> apply(
73  const Tensor<T, Symmetry<1, 1>, tmpl::list<Index0, Index0>>&
74  tensor) noexcept {
75  const T& t00 = get<0, 0>(tensor);
76  const T& t01 = get<0, 1>(tensor);
77  const T& t11 = get<1, 1>(tensor);
78  // det is non-const so that it can be moved into the std::pair:
79  Scalar<T> det{t00 * t11 - t01 * t01};
80  const T one_over_det = make_with_value<T>(det, 1.0) / det.get();
81  Tensor<T, Symmetry<1, 1>, inverse_indices<Index0, Index0>> inv{};
82  get<0, 0>(inv) = t11 * one_over_det;
83  get<0, 1>(inv) = -t01 * one_over_det;
84  get<1, 1>(inv) = t00 * one_over_det;
85  return std::make_pair(std::move(det), std::move(inv));
86  }
87 };
88 
89 // The inverse of a 3x3 tensor is computed from Cramer's rule. By reusing some
90 // terms, the determinant is computed efficiently at the same time.
91 template <typename Index0, typename Index1>
92 struct DetAndInverseImpl<Symmetry<2, 1>, Index0, Index1,
93  Requires<Index0::dim == 3>> {
94  template <typename T>
95  static determinant_inverse_pair<T, Symmetry<2, 1>, Index0, Index1> apply(
96  const Tensor<T, Symmetry<2, 1>, tmpl::list<Index0, Index1>>&
97  tensor) noexcept {
98  const T& t00 = get<0, 0>(tensor);
99  const T& t01 = get<0, 1>(tensor);
100  const T& t02 = get<0, 2>(tensor);
101  const T& t10 = get<1, 0>(tensor);
102  const T& t11 = get<1, 1>(tensor);
103  const T& t12 = get<1, 2>(tensor);
104  const T& t20 = get<2, 0>(tensor);
105  const T& t21 = get<2, 1>(tensor);
106  const T& t22 = get<2, 2>(tensor);
107  const T a = t11 * t22 - t12 * t21;
108  const T b = t12 * t20 - t10 * t22;
109  const T c = t10 * t21 - t11 * t20;
110  // det is non-const so that it can be moved into the std::pair:
111  Scalar<T> det{t00 * a + t01 * b + t02 * c};
112  const T one_over_det = make_with_value<T>(det, 1.0) / det.get();
113  Tensor<T, Symmetry<2, 1>, inverse_indices<Index0, Index1>> inv{};
114  get<0, 0>(inv) = a * one_over_det;
115  get<0, 1>(inv) = (t21 * t02 - t22 * t01) * one_over_det;
116  get<0, 2>(inv) = (t01 * t12 - t02 * t11) * one_over_det;
117  get<1, 0>(inv) = b * one_over_det;
118  get<1, 1>(inv) = (t22 * t00 - t20 * t02) * one_over_det;
119  get<1, 2>(inv) = (t02 * t10 - t00 * t12) * one_over_det;
120  get<2, 0>(inv) = c * one_over_det;
121  get<2, 1>(inv) = (t20 * t01 - t21 * t00) * one_over_det;
122  get<2, 2>(inv) = (t00 * t11 - t01 * t10) * one_over_det;
123  return std::make_pair(std::move(det), std::move(inv));
124  }
125 };
126 
127 template <typename Index0>
128 struct DetAndInverseImpl<Symmetry<1, 1>, Index0, Index0,
129  Requires<Index0::dim == 3>> {
130  template <typename T>
131  static determinant_inverse_pair<T, Symmetry<1, 1>, Index0, Index0> apply(
132  const Tensor<T, Symmetry<1, 1>, tmpl::list<Index0, Index0>>&
133  tensor) noexcept {
134  const T& t00 = get<0, 0>(tensor);
135  const T& t01 = get<0, 1>(tensor);
136  const T& t02 = get<0, 2>(tensor);
137  const T& t11 = get<1, 1>(tensor);
138  const T& t12 = get<1, 2>(tensor);
139  const T& t22 = get<2, 2>(tensor);
140  const T a = t11 * t22 - t12 * t12;
141  const T b = t12 * t02 - t01 * t22;
142  const T c = t01 * t12 - t11 * t02;
143  // det is non-const so that it can be moved into the std::pair:
144  Scalar<T> det{t00 * a + t01 * b + t02 * c};
145  const T one_over_det = make_with_value<T>(det, 1.0) / det.get();
146  Tensor<T, Symmetry<1, 1>, inverse_indices<Index0, Index0>> inv{};
147  get<0, 0>(inv) = (t11 * t22 - t12 * t12) * one_over_det;
148  get<0, 1>(inv) = (t12 * t02 - t22 * t01) * one_over_det;
149  get<0, 2>(inv) = (t01 * t12 - t02 * t11) * one_over_det;
150  get<1, 1>(inv) = (t22 * t00 - t02 * t02) * one_over_det;
151  get<1, 2>(inv) = (t02 * t01 - t00 * t12) * one_over_det;
152  get<2, 2>(inv) = (t00 * t11 - t01 * t01) * one_over_det;
153  return std::make_pair(std::move(det), std::move(inv));
154  }
155 };
156 
157 // The 4x4 tensor inverse is implemented using the formula for inverting a
158 // partitioned matrix (here, partitioned into 2x2 blocks). This is more
159 // efficient than Cramer's rule for matrices larger than 3x3.
160 //
161 // In this algorithm, the 4x4 input matrix is partitioned into the 2x2 blocks:
162 // P Q
163 // R S
164 // The 4x4 inverse matrix is partitioned into the 2x2 blocks:
165 // U V
166 // W X
167 // Each inverse block is obtained from the blocks {P,Q,R,S} of the input matrix
168 // as follows:
169 // X = inv[S - R inv(P) Q] (i.e. inv(X) is the Schur complement of P)
170 // W = - X R inv(P)
171 // V = - inv(P) Q X
172 // U = inv(P) + inv(P) Q X R inv(P) = inv(P) - inv(P) Q W
173 template <typename Index0, typename Index1>
174 struct DetAndInverseImpl<Symmetry<2, 1>, Index0, Index1,
175  Requires<Index0::dim == 4>> {
176  template <typename T>
177  static determinant_inverse_pair<T, Symmetry<2, 1>, Index0, Index1> apply(
178  const Tensor<T, Symmetry<2, 1>, tmpl::list<Index0, Index1>>&
179  tensor) noexcept {
180  const T& p00 = get<0, 0>(tensor);
181  const T& p01 = get<0, 1>(tensor);
182  const T& p10 = get<1, 0>(tensor);
183  const T& p11 = get<1, 1>(tensor);
184  const T& q00 = get<0, 2>(tensor);
185  const T& q01 = get<0, 3>(tensor);
186  const T& q10 = get<1, 2>(tensor);
187  const T& q11 = get<1, 3>(tensor);
188  const T& r00 = get<2, 0>(tensor);
189  const T& r01 = get<2, 1>(tensor);
190  const T& r10 = get<3, 0>(tensor);
191  const T& r11 = get<3, 1>(tensor);
192  const T& s00 = get<2, 2>(tensor);
193  const T& s01 = get<2, 3>(tensor);
194  const T& s10 = get<3, 2>(tensor);
195  const T& s11 = get<3, 3>(tensor);
196 
197  Tensor<T, Symmetry<2, 1>, inverse_indices<Index0, Index1>> inv{};
198  T& u00 = get<0, 0>(inv);
199  T& u01 = get<0, 1>(inv);
200  T& u10 = get<1, 0>(inv);
201  T& u11 = get<1, 1>(inv);
202  T& v00 = get<0, 2>(inv);
203  T& v01 = get<0, 3>(inv);
204  T& v10 = get<1, 2>(inv);
205  T& v11 = get<1, 3>(inv);
206  T& w00 = get<2, 0>(inv);
207  T& w01 = get<2, 1>(inv);
208  T& w10 = get<3, 0>(inv);
209  T& w11 = get<3, 1>(inv);
210  T& x00 = get<2, 2>(inv);
211  T& x01 = get<2, 3>(inv);
212  T& x10 = get<3, 2>(inv);
213  T& x11 = get<3, 3>(inv);
214 
215  const T det_p = p00 * p11 - p01 * p10;
216  const T one_over_det_p = make_with_value<T>(det_p, 1.0) / det_p;
217  const T inv_p00 = p11 * one_over_det_p;
218  const T inv_p01 = -p01 * one_over_det_p;
219  const T inv_p10 = -p10 * one_over_det_p;
220  const T inv_p11 = p00 * one_over_det_p;
221 
222  const T r_inv_p00 = r00 * inv_p00 + r01 * inv_p10;
223  const T r_inv_p01 = r00 * inv_p01 + r01 * inv_p11;
224  const T r_inv_p10 = r10 * inv_p00 + r11 * inv_p10;
225  const T r_inv_p11 = r10 * inv_p01 + r11 * inv_p11;
226 
227  const T inv_p_q00 = inv_p00 * q00 + inv_p01 * q10;
228  const T inv_p_q01 = inv_p00 * q01 + inv_p01 * q11;
229  const T inv_p_q10 = inv_p10 * q00 + inv_p11 * q10;
230  const T inv_p_q11 = inv_p10 * q01 + inv_p11 * q11;
231 
232  const T inv_x00 = s00 - (r_inv_p00 * q00 + r_inv_p01 * q10);
233  const T inv_x01 = s01 - (r_inv_p00 * q01 + r_inv_p01 * q11);
234  const T inv_x10 = s10 - (r_inv_p10 * q00 + r_inv_p11 * q10);
235  const T inv_x11 = s11 - (r_inv_p10 * q01 + r_inv_p11 * q11);
236 
237  const T det_inv_x = inv_x00 * inv_x11 - inv_x01 * inv_x10;
238  const T one_over_det_inv_x = make_with_value<T>(det_inv_x, 1.0) / det_inv_x;
239  x00 = inv_x11 * one_over_det_inv_x;
240  x01 = -inv_x01 * one_over_det_inv_x;
241  x10 = -inv_x10 * one_over_det_inv_x;
242  x11 = inv_x00 * one_over_det_inv_x;
243 
244  w00 = -x00 * r_inv_p00 - x01 * r_inv_p10;
245  w01 = -x00 * r_inv_p01 - x01 * r_inv_p11;
246  w10 = -x10 * r_inv_p00 - x11 * r_inv_p10;
247  w11 = -x10 * r_inv_p01 - x11 * r_inv_p11;
248 
249  v00 = -inv_p_q00 * x00 - inv_p_q01 * x10;
250  v01 = -inv_p_q00 * x01 - inv_p_q01 * x11;
251  v10 = -inv_p_q10 * x00 - inv_p_q11 * x10;
252  v11 = -inv_p_q10 * x01 - inv_p_q11 * x11;
253 
254  u00 = inv_p00 - (inv_p_q00 * w00 + inv_p_q01 * w10);
255  u01 = inv_p01 - (inv_p_q00 * w01 + inv_p_q01 * w11);
256  u10 = inv_p10 - (inv_p_q10 * w00 + inv_p_q11 * w10);
257  u11 = inv_p11 - (inv_p_q10 * w01 + inv_p_q11 * w11);
258 
259  // det is non-const so that it can be moved into the std::pair:
260  Scalar<T> det{det_p * det_inv_x};
261  return std::make_pair(std::move(det), std::move(inv));
262  }
263 };
264 
265 template <typename Index0>
266 struct DetAndInverseImpl<Symmetry<1, 1>, Index0, Index0,
267  Requires<Index0::dim == 4>> {
268  template <typename T>
269  static determinant_inverse_pair<T, Symmetry<1, 1>, Index0, Index0> apply(
270  const Tensor<T, Symmetry<1, 1>, tmpl::list<Index0, Index0>>&
271  tensor) noexcept {
272  const T& p00 = get<0, 0>(tensor);
273  const T& p01 = get<0, 1>(tensor);
274  const T& p11 = get<1, 1>(tensor);
275  const T& q00 = get<0, 2>(tensor);
276  const T& q01 = get<0, 3>(tensor);
277  const T& q10 = get<1, 2>(tensor);
278  const T& q11 = get<1, 3>(tensor);
279  const T& s00 = get<2, 2>(tensor);
280  const T& s01 = get<2, 3>(tensor);
281  const T& s11 = get<3, 3>(tensor);
282 
283  Tensor<T, Symmetry<1, 1>, inverse_indices<Index0, Index0>> inv{};
284  T& u00 = get<0, 0>(inv);
285  T& u01 = get<0, 1>(inv);
286  T& u11 = get<1, 1>(inv);
287  T& v00 = get<0, 2>(inv);
288  T& v01 = get<0, 3>(inv);
289  T& v10 = get<1, 2>(inv);
290  T& v11 = get<1, 3>(inv);
291  T& x00 = get<2, 2>(inv);
292  T& x01 = get<2, 3>(inv);
293  T& x11 = get<3, 3>(inv);
294 
295  const T det_p = p00 * p11 - p01 * p01;
296  const T one_over_det_p = make_with_value<T>(det_p, 1.0) / det_p;
297  const T inv_p00 = p11 * one_over_det_p;
298  const T inv_p01 = -p01 * one_over_det_p;
299  const T inv_p11 = p00 * one_over_det_p;
300 
301  const T r_inv_p00 = q00 * inv_p00 + q10 * inv_p01;
302  const T r_inv_p01 = q00 * inv_p01 + q10 * inv_p11;
303  const T r_inv_p10 = q01 * inv_p00 + q11 * inv_p01;
304  const T r_inv_p11 = q01 * inv_p01 + q11 * inv_p11;
305 
306  const T inv_x00 = s00 - (r_inv_p00 * q00 + r_inv_p01 * q10);
307  const T inv_x01 = s01 - (r_inv_p00 * q01 + r_inv_p01 * q11);
308  const T inv_x11 = s11 - (r_inv_p10 * q01 + r_inv_p11 * q11);
309 
310  const T det_inv_x = inv_x00 * inv_x11 - inv_x01 * inv_x01;
311  const T one_over_det_inv_x = make_with_value<T>(det_inv_x, 1.0) / det_inv_x;
312  x00 = inv_x11 * one_over_det_inv_x;
313  x01 = -inv_x01 * one_over_det_inv_x;
314  x11 = inv_x00 * one_over_det_inv_x;
315 
316  v00 = -x00 * r_inv_p00 - x01 * r_inv_p10;
317  v10 = -x00 * r_inv_p01 - x01 * r_inv_p11;
318  v01 = -x01 * r_inv_p00 - x11 * r_inv_p10;
319  v11 = -x01 * r_inv_p01 - x11 * r_inv_p11;
320 
321  u00 = inv_p00 - (r_inv_p00 * v00 + r_inv_p10 * v01);
322  u01 = inv_p01 - (r_inv_p00 * v10 + r_inv_p10 * v11);
323  u11 = inv_p11 - (r_inv_p01 * v10 + r_inv_p11 * v11);
324 
325  // det is non-const so that it can be moved into the std::pair:
326  Scalar<T> det{det_p * det_inv_x};
327  return std::make_pair(std::move(det), std::move(inv));
328  }
329 };
330 } // namespace determinant_and_inverse_detail
331 
332 /*!
333  * \ingroup TensorGroup
334  * \brief Computes the determinant and inverse of a rank-2 Tensor.
335  *
336  * Computes the determinant and inverse together, because this leads to fewer
337  * operations compared to computing the determinant independently.
338  *
339  * \param tensor the input rank-2 Tensor.
340  * \return a std::pair that holds the determinant (in `pair.first`) and inverse
341  * (in `pair.second`) of the input tensor.
342  *
343  * \details
344  * Treats the input rank-2 tensor as a matrix. The first (second) index of the
345  * tensor corresponds to the rows (columns) of the matrix. The determinant is a
346  * scalar tensor. The inverse is a rank-2 tensor whose indices are reversed and
347  * of opposite valence relative to the input tensor, i.e. given \f$T_a^b\f$
348  * returns \f$(Tinv)_b^a\f$.
349  *
350  * \note
351  * When inverting a 4x4 spacetime metric, it is typically more efficient to use
352  * the 3+1 decomposition of the 4-metric in terms of lapse, shift, and spatial
353  * 3-metric, in which only the spatial 3-metric needs to be inverted.
354  */
355 template <typename T, typename Symm, typename Index0, typename Index1>
357  const Tensor<T, Symm, tmpl::list<Index0, Index1>>& tensor) noexcept
358  -> std::pair<Scalar<T>, Tensor<T, Symm,
359  tmpl::list<change_index_up_lo<Index1>,
360  change_index_up_lo<Index0>>>> {
361  static_assert(Index0::dim == Index1::dim,
362  "Cannot take the inverse of a Tensor whose Indices are not "
363  "of the same dimensionality.");
364  static_assert(Index0::index_type == Index1::index_type,
365  "Taking the inverse of a mixed Spatial and Spacetime index "
366  "Tensor is not allowed since it's not clear what that means.");
367  static_assert(not std::is_integral<T>::value, "Can't invert a Tensor<int>.");
368  return determinant_and_inverse_detail::DetAndInverseImpl<
369  Symm, Index0, Index1>::apply(tensor);
370 }
Defines the type alias Requires.
constexpr auto apply(F &&f, const DataBox< BoxTags > &box, Args &&... args)
Apply the function f with argument Tags TagsList from DataBox box
Definition: DataBox.hpp:1595
Definition: DeterminantAndInverse.hpp:18
typename detail::SymmetryImpl< std::make_index_sequence< sizeof...(T)>, tmpl::integral_list< std::int32_t, T... > >::type Symmetry
Computes the canonical symmetry from the integers T
Definition: Symmetry.hpp:81
auto determinant_and_inverse(const Tensor< T, Symm, tmpl::list< Index0, Index1 >> &tensor) noexcept -> std::pair< Scalar< T >, Tensor< T, Symm, tmpl::list< change_index_up_lo< Index1 >, change_index_up_lo< Index0 >>>>
Computes the determinant and inverse of a rank-2 Tensor.
Definition: DeterminantAndInverse.hpp:356
Defines classes for Tensor.
Tensor_detail::TensorIndexType< Index::index_type==IndexType::Spatial ? Index::value :Index::value - 1, Index::ul==UpLo::Up ? UpLo::Lo :UpLo::Up, typename Index::Frame, Index::index_type > change_index_up_lo
Change the TensorIndexType to be covariant if it&#39;s contravariant and vice-versa.
Definition: IndexType.hpp:233
Defines macro ASSERT.
typename Requires_detail::requires_impl< B >::template_error_type_failed_to_meet_requirements_on_template_parameters Requires
Express requirements on the template parameters of a function or class, replaces std::enable_if_t ...
Definition: Requires.hpp:67
Tensor< T, Symmetry<>, index_list<> > Scalar
Scalar type.
Definition: TypeAliases.hpp:21
Defines make_with_value.