conspire/math/tensor/rank_2/inverse/
mod.rs

1#[cfg(test)]
2mod test;
3use crate::math::{Factor, Reference};
4use crate::units::{Dimensionless, UnitInv};
5
6use super::{Quantity, Rank2, Tensor, TensorArray, TensorRank0, TensorRank2, relabel};
7use crate::ABS_TOL;
8
9/// The factors of an LU decomposition, with the permutation applied to the rows.
10type LuFactors<const D: usize, I, J, U> = (
11    TensorRank2<D, I, Factor, U>,
12    TensorRank2<D, Factor, J, U>,
13    Vec<usize>,
14);
15
16/// The unit an inverse carries.
17type Inverse<U> = <U as UnitInv>::Output;
18
19impl<const D: usize, I, J, U> TensorRank2<D, I, J, U> {
20    /// Returns the determinant of the rank-2 tensor.
21    pub fn determinant(&self) -> TensorRank0 {
22        self.canonical().determinant_core()
23    }
24    /// Returns the inverse of the rank-2 tensor.
25    pub fn inverse(&self) -> TensorRank2<D, J, I, Inverse<U>>
26    where
27        U: UnitInv,
28    {
29        relabel(self.canonical().inverse_core())
30    }
31    /// Returns the inverse and determinant of the rank-2 tensor.
32    pub fn inverse_and_determinant(&self) -> (TensorRank2<D, J, I, Inverse<U>>, TensorRank0)
33    where
34        U: UnitInv,
35    {
36        let (inverse, determinant) = self.canonical().inverse_and_determinant_core();
37        (relabel(inverse), determinant)
38    }
39    /// Returns the inverse transpose of the rank-2 tensor.
40    pub fn inverse_transpose(&self) -> TensorRank2<D, I, J, Inverse<U>>
41    where
42        U: UnitInv,
43    {
44        relabel(self.canonical().inverse_transpose_core())
45    }
46    /// Returns the inverse transpose and determinant of the rank-2 tensor.
47    pub fn inverse_transpose_and_determinant(
48        &self,
49    ) -> (TensorRank2<D, I, J, Inverse<U>>, TensorRank0)
50    where
51        U: UnitInv,
52    {
53        let (inverse_transpose, determinant) =
54            self.canonical().inverse_transpose_and_determinant_core();
55        (relabel(inverse_transpose), determinant)
56    }
57    /// Returns the LU decomposition of the rank-2 tensor.
58    pub fn lu_decomposition(&self) -> LuFactors<D, I, J, U> {
59        let (tensor_l, tensor_u, p) = self.canonical().lu_decomposition_core();
60        (relabel(tensor_l), relabel(tensor_u), p)
61    }
62    /// Returns the inverse of the LU decomposition of the rank-2 tensor.
63    pub fn lu_decomposition_inverse(&self) -> LuFactors<D, I, J, U> {
64        let (tensor_l, tensor_u, p) = self.canonical().lu_decomposition_inverse_core();
65        (relabel(tensor_l), relabel(tensor_u), p)
66    }
67}
68
69impl<const D: usize> TensorRank2<D, Reference, Reference, Dimensionless> {
70    fn determinant_core(&self) -> TensorRank0 {
71        if D == 2 {
72            (self[0][0] * self[1][1] - self[0][1] * self[1][0]).value()
73        } else if D == 3 {
74            let c_00 = self[1][1] * self[2][2] - self[1][2] * self[2][1];
75            let c_10 = self[1][2] * self[2][0] - self[1][0] * self[2][2];
76            let c_20 = self[1][0] * self[2][1] - self[1][1] * self[2][0];
77            (self[0][0] * c_00 + self[0][1] * c_10 + self[0][2] * c_20).value()
78        } else if D == 4 {
79            let s0 = self[0][0] * self[1][1] - self[0][1] * self[1][0];
80            let s1 = self[0][0] * self[1][2] - self[0][2] * self[1][0];
81            let s2 = self[0][0] * self[1][3] - self[0][3] * self[1][0];
82            let s3 = self[0][1] * self[1][2] - self[0][2] * self[1][1];
83            let s4 = self[0][1] * self[1][3] - self[0][3] * self[1][1];
84            let s5 = self[0][2] * self[1][3] - self[0][3] * self[1][2];
85            let c5 = self[2][2] * self[3][3] - self[2][3] * self[3][2];
86            let c4 = self[2][1] * self[3][3] - self[2][3] * self[3][1];
87            let c3 = self[2][1] * self[3][2] - self[2][2] * self[3][1];
88            let c2 = self[2][0] * self[3][3] - self[2][3] * self[3][0];
89            let c1 = self[2][0] * self[3][2] - self[2][2] * self[3][0];
90            let c0 = self[2][0] * self[3][1] - self[2][1] * self[3][0];
91            (s0 * c5 - s1 * c4 + s2 * c3 + s3 * c2 - s4 * c1 + s5 * c0).value()
92        } else {
93            let (_, u, p) = self.lu_decomposition_core();
94            let num_swaps = p.iter().enumerate().filter(|(i, p_i)| p_i != &i).count();
95            u.into_iter()
96                .enumerate()
97                .map(|(i, u_i)| u_i[i].value())
98                .product::<TensorRank0>()
99                * if num_swaps % 2 == 0 { 1.0 } else { -1.0 }
100        }
101    }
102    fn inverse_core(&self) -> Self {
103        if D == 2 {
104            let mut adjugate = Self::zero();
105            adjugate[0][0] = self[1][1];
106            adjugate[0][1] = -self[0][1];
107            adjugate[1][0] = -self[1][0];
108            adjugate[1][1] = self[0][0];
109            adjugate / self.determinant_core()
110        } else if D == 3 {
111            let mut adjugate = Self::zero();
112            let c_00 = self[1][1] * self[2][2] - self[1][2] * self[2][1];
113            let c_10 = self[1][2] * self[2][0] - self[1][0] * self[2][2];
114            let c_20 = self[1][0] * self[2][1] - self[1][1] * self[2][0];
115            adjugate[0][0] = c_00;
116            adjugate[0][1] = self[0][2] * self[2][1] - self[0][1] * self[2][2];
117            adjugate[0][2] = self[0][1] * self[1][2] - self[0][2] * self[1][1];
118            adjugate[1][0] = c_10;
119            adjugate[1][1] = self[0][0] * self[2][2] - self[0][2] * self[2][0];
120            adjugate[1][2] = self[0][2] * self[1][0] - self[0][0] * self[1][2];
121            adjugate[2][0] = c_20;
122            adjugate[2][1] = self[0][1] * self[2][0] - self[0][0] * self[2][1];
123            adjugate[2][2] = self[0][0] * self[1][1] - self[0][1] * self[1][0];
124            adjugate / (self[0][0] * c_00 + self[0][1] * c_10 + self[0][2] * c_20)
125        } else if D == 4 {
126            let mut adjugate = Self::zero();
127            let s0 = self[0][0] * self[1][1] - self[0][1] * self[1][0];
128            let s1 = self[0][0] * self[1][2] - self[0][2] * self[1][0];
129            let s2 = self[0][0] * self[1][3] - self[0][3] * self[1][0];
130            let s3 = self[0][1] * self[1][2] - self[0][2] * self[1][1];
131            let s4 = self[0][1] * self[1][3] - self[0][3] * self[1][1];
132            let s5 = self[0][2] * self[1][3] - self[0][3] * self[1][2];
133            let c5 = self[2][2] * self[3][3] - self[2][3] * self[3][2];
134            let c4 = self[2][1] * self[3][3] - self[2][3] * self[3][1];
135            let c3 = self[2][1] * self[3][2] - self[2][2] * self[3][1];
136            let c2 = self[2][0] * self[3][3] - self[2][3] * self[3][0];
137            let c1 = self[2][0] * self[3][2] - self[2][2] * self[3][0];
138            let c0 = self[2][0] * self[3][1] - self[2][1] * self[3][0];
139            adjugate[0][0] = self[1][1] * c5 - self[1][2] * c4 + self[1][3] * c3;
140            adjugate[0][1] = self[0][2] * c4 - self[0][1] * c5 - self[0][3] * c3;
141            adjugate[0][2] = self[3][1] * s5 - self[3][2] * s4 + self[3][3] * s3;
142            adjugate[0][3] = self[2][2] * s4 - self[2][1] * s5 - self[2][3] * s3;
143            adjugate[1][0] = self[1][2] * c2 - self[1][0] * c5 - self[1][3] * c1;
144            adjugate[1][1] = self[0][0] * c5 - self[0][2] * c2 + self[0][3] * c1;
145            adjugate[1][2] = self[3][2] * s2 - self[3][0] * s5 - self[3][3] * s1;
146            adjugate[1][3] = self[2][0] * s5 - self[2][2] * s2 + self[2][3] * s1;
147            adjugate[2][0] = self[1][0] * c4 - self[1][1] * c2 + self[1][3] * c0;
148            adjugate[2][1] = self[0][1] * c2 - self[0][0] * c4 - self[0][3] * c0;
149            adjugate[2][2] = self[3][0] * s4 - self[3][1] * s2 + self[3][3] * s0;
150            adjugate[2][3] = self[2][1] * s2 - self[2][0] * s4 - self[2][3] * s0;
151            adjugate[3][0] = self[1][1] * c1 - self[1][0] * c3 - self[1][2] * c0;
152            adjugate[3][1] = self[0][0] * c3 - self[0][1] * c1 + self[0][2] * c0;
153            adjugate[3][2] = self[3][1] * s1 - self[3][0] * s3 - self[3][2] * s0;
154            adjugate[3][3] = self[2][0] * s3 - self[2][1] * s1 + self[2][2] * s0;
155            adjugate / (s0 * c5 - s1 * c4 + s2 * c3 + s3 * c2 - s4 * c1 + s5 * c0)
156        } else {
157            let (l_inverse, u_inverse, p) = self.lu_decomposition_inverse_core();
158            let mut q = [0; D];
159            p.into_iter().enumerate().for_each(|(i, p_i)| q[p_i] = i);
160            u_inverse
161                .into_iter()
162                .map(|u_inverse_i| {
163                    q.iter()
164                        .map(|&q_j| {
165                            u_inverse_i
166                                .iter()
167                                .zip(l_inverse.iter())
168                                .map(|(u_inverse_ik, l_inverse_k)| u_inverse_ik * l_inverse_k[q_j])
169                                .sum::<Quantity>()
170                        })
171                        .collect()
172                })
173                .collect()
174        }
175    }
176    fn inverse_and_determinant_core(&self) -> (Self, TensorRank0) {
177        if D == 2 {
178            let mut adjugate = Self::zero();
179            adjugate[0][0] = self[1][1];
180            adjugate[0][1] = -self[0][1];
181            adjugate[1][0] = -self[1][0];
182            adjugate[1][1] = self[0][0];
183            let determinant = self.determinant_core();
184            (adjugate / determinant, determinant)
185        } else if D == 3 {
186            let mut adjugate = Self::zero();
187            let c_00 = self[1][1] * self[2][2] - self[1][2] * self[2][1];
188            let c_10 = self[1][2] * self[2][0] - self[1][0] * self[2][2];
189            let c_20 = self[1][0] * self[2][1] - self[1][1] * self[2][0];
190            let determinant = (self[0][0] * c_00 + self[0][1] * c_10 + self[0][2] * c_20).value();
191            adjugate[0][0] = c_00;
192            adjugate[0][1] = self[0][2] * self[2][1] - self[0][1] * self[2][2];
193            adjugate[0][2] = self[0][1] * self[1][2] - self[0][2] * self[1][1];
194            adjugate[1][0] = c_10;
195            adjugate[1][1] = self[0][0] * self[2][2] - self[0][2] * self[2][0];
196            adjugate[1][2] = self[0][2] * self[1][0] - self[0][0] * self[1][2];
197            adjugate[2][0] = c_20;
198            adjugate[2][1] = self[0][1] * self[2][0] - self[0][0] * self[2][1];
199            adjugate[2][2] = self[0][0] * self[1][1] - self[0][1] * self[1][0];
200            (adjugate / determinant, determinant)
201        } else if D == 4 {
202            let mut adjugate = Self::zero();
203            let s0 = self[0][0] * self[1][1] - self[0][1] * self[1][0];
204            let s1 = self[0][0] * self[1][2] - self[0][2] * self[1][0];
205            let s2 = self[0][0] * self[1][3] - self[0][3] * self[1][0];
206            let s3 = self[0][1] * self[1][2] - self[0][2] * self[1][1];
207            let s4 = self[0][1] * self[1][3] - self[0][3] * self[1][1];
208            let s5 = self[0][2] * self[1][3] - self[0][3] * self[1][2];
209            let c5 = self[2][2] * self[3][3] - self[2][3] * self[3][2];
210            let c4 = self[2][1] * self[3][3] - self[2][3] * self[3][1];
211            let c3 = self[2][1] * self[3][2] - self[2][2] * self[3][1];
212            let c2 = self[2][0] * self[3][3] - self[2][3] * self[3][0];
213            let c1 = self[2][0] * self[3][2] - self[2][2] * self[3][0];
214            let c0 = self[2][0] * self[3][1] - self[2][1] * self[3][0];
215            let determinant = (s0 * c5 - s1 * c4 + s2 * c3 + s3 * c2 - s4 * c1 + s5 * c0).value();
216            adjugate[0][0] = self[1][1] * c5 - self[1][2] * c4 + self[1][3] * c3;
217            adjugate[0][1] = self[0][2] * c4 - self[0][1] * c5 - self[0][3] * c3;
218            adjugate[0][2] = self[3][1] * s5 - self[3][2] * s4 + self[3][3] * s3;
219            adjugate[0][3] = self[2][2] * s4 - self[2][1] * s5 - self[2][3] * s3;
220            adjugate[1][0] = self[1][2] * c2 - self[1][0] * c5 - self[1][3] * c1;
221            adjugate[1][1] = self[0][0] * c5 - self[0][2] * c2 + self[0][3] * c1;
222            adjugate[1][2] = self[3][2] * s2 - self[3][0] * s5 - self[3][3] * s1;
223            adjugate[1][3] = self[2][0] * s5 - self[2][2] * s2 + self[2][3] * s1;
224            adjugate[2][0] = self[1][0] * c4 - self[1][1] * c2 + self[1][3] * c0;
225            adjugate[2][1] = self[0][1] * c2 - self[0][0] * c4 - self[0][3] * c0;
226            adjugate[2][2] = self[3][0] * s4 - self[3][1] * s2 + self[3][3] * s0;
227            adjugate[2][3] = self[2][1] * s2 - self[2][0] * s4 - self[2][3] * s0;
228            adjugate[3][0] = self[1][1] * c1 - self[1][0] * c3 - self[1][2] * c0;
229            adjugate[3][1] = self[0][0] * c3 - self[0][1] * c1 + self[0][2] * c0;
230            adjugate[3][2] = self[3][1] * s1 - self[3][0] * s3 - self[3][2] * s0;
231            adjugate[3][3] = self[2][0] * s3 - self[2][1] * s1 + self[2][2] * s0;
232            (adjugate / determinant, determinant)
233        } else {
234            (self.inverse_core(), self.determinant_core())
235        }
236    }
237    fn inverse_transpose_core(&self) -> Self {
238        if D == 2 {
239            let mut adjugate_transpose = Self::zero();
240            adjugate_transpose[0][0] = self[1][1];
241            adjugate_transpose[0][1] = -self[1][0];
242            adjugate_transpose[1][0] = -self[0][1];
243            adjugate_transpose[1][1] = self[0][0];
244            adjugate_transpose / self.determinant_core()
245        } else if D == 3 {
246            let mut adjugate_transpose = Self::zero();
247            let c_00 = self[1][1] * self[2][2] - self[1][2] * self[2][1];
248            let c_10 = self[1][2] * self[2][0] - self[1][0] * self[2][2];
249            let c_20 = self[1][0] * self[2][1] - self[1][1] * self[2][0];
250            adjugate_transpose[0][0] = c_00;
251            adjugate_transpose[1][0] = self[0][2] * self[2][1] - self[0][1] * self[2][2];
252            adjugate_transpose[2][0] = self[0][1] * self[1][2] - self[0][2] * self[1][1];
253            adjugate_transpose[0][1] = c_10;
254            adjugate_transpose[1][1] = self[0][0] * self[2][2] - self[0][2] * self[2][0];
255            adjugate_transpose[2][1] = self[0][2] * self[1][0] - self[0][0] * self[1][2];
256            adjugate_transpose[0][2] = c_20;
257            adjugate_transpose[1][2] = self[0][1] * self[2][0] - self[0][0] * self[2][1];
258            adjugate_transpose[2][2] = self[0][0] * self[1][1] - self[0][1] * self[1][0];
259            adjugate_transpose / (self[0][0] * c_00 + self[0][1] * c_10 + self[0][2] * c_20)
260        } else if D == 4 {
261            let mut adjugate_transpose = Self::zero();
262            let s0 = self[0][0] * self[1][1] - self[0][1] * self[1][0];
263            let s1 = self[0][0] * self[1][2] - self[0][2] * self[1][0];
264            let s2 = self[0][0] * self[1][3] - self[0][3] * self[1][0];
265            let s3 = self[0][1] * self[1][2] - self[0][2] * self[1][1];
266            let s4 = self[0][1] * self[1][3] - self[0][3] * self[1][1];
267            let s5 = self[0][2] * self[1][3] - self[0][3] * self[1][2];
268            let c5 = self[2][2] * self[3][3] - self[2][3] * self[3][2];
269            let c4 = self[2][1] * self[3][3] - self[2][3] * self[3][1];
270            let c3 = self[2][1] * self[3][2] - self[2][2] * self[3][1];
271            let c2 = self[2][0] * self[3][3] - self[2][3] * self[3][0];
272            let c1 = self[2][0] * self[3][2] - self[2][2] * self[3][0];
273            let c0 = self[2][0] * self[3][1] - self[2][1] * self[3][0];
274            adjugate_transpose[0][0] = self[1][1] * c5 - self[1][2] * c4 + self[1][3] * c3;
275            adjugate_transpose[1][0] = self[0][2] * c4 - self[0][1] * c5 - self[0][3] * c3;
276            adjugate_transpose[2][0] = self[3][1] * s5 - self[3][2] * s4 + self[3][3] * s3;
277            adjugate_transpose[3][0] = self[2][2] * s4 - self[2][1] * s5 - self[2][3] * s3;
278            adjugate_transpose[0][1] = self[1][2] * c2 - self[1][0] * c5 - self[1][3] * c1;
279            adjugate_transpose[1][1] = self[0][0] * c5 - self[0][2] * c2 + self[0][3] * c1;
280            adjugate_transpose[2][1] = self[3][2] * s2 - self[3][0] * s5 - self[3][3] * s1;
281            adjugate_transpose[3][1] = self[2][0] * s5 - self[2][2] * s2 + self[2][3] * s1;
282            adjugate_transpose[0][2] = self[1][0] * c4 - self[1][1] * c2 + self[1][3] * c0;
283            adjugate_transpose[1][2] = self[0][1] * c2 - self[0][0] * c4 - self[0][3] * c0;
284            adjugate_transpose[2][2] = self[3][0] * s4 - self[3][1] * s2 + self[3][3] * s0;
285            adjugate_transpose[3][2] = self[2][1] * s2 - self[2][0] * s4 - self[2][3] * s0;
286            adjugate_transpose[0][3] = self[1][1] * c1 - self[1][0] * c3 - self[1][2] * c0;
287            adjugate_transpose[1][3] = self[0][0] * c3 - self[0][1] * c1 + self[0][2] * c0;
288            adjugate_transpose[2][3] = self[3][1] * s1 - self[3][0] * s3 - self[3][2] * s0;
289            adjugate_transpose[3][3] = self[2][0] * s3 - self[2][1] * s1 + self[2][2] * s0;
290            adjugate_transpose / (s0 * c5 - s1 * c4 + s2 * c3 + s3 * c2 - s4 * c1 + s5 * c0)
291        } else {
292            self.inverse_core().transpose()
293        }
294    }
295    fn inverse_transpose_and_determinant_core(&self) -> (Self, TensorRank0) {
296        if D == 2 {
297            let mut adjugate_transpose = Self::zero();
298            adjugate_transpose[0][0] = self[1][1];
299            adjugate_transpose[0][1] = -self[1][0];
300            adjugate_transpose[1][0] = -self[0][1];
301            adjugate_transpose[1][1] = self[0][0];
302            let determinant = self.determinant_core();
303            (adjugate_transpose / determinant, determinant)
304        } else if D == 3 {
305            let mut adjugate_transpose = Self::zero();
306            let c_00 = self[1][1] * self[2][2] - self[1][2] * self[2][1];
307            let c_10 = self[1][2] * self[2][0] - self[1][0] * self[2][2];
308            let c_20 = self[1][0] * self[2][1] - self[1][1] * self[2][0];
309            let determinant = (self[0][0] * c_00 + self[0][1] * c_10 + self[0][2] * c_20).value();
310            adjugate_transpose[0][0] = c_00;
311            adjugate_transpose[1][0] = self[0][2] * self[2][1] - self[0][1] * self[2][2];
312            adjugate_transpose[2][0] = self[0][1] * self[1][2] - self[0][2] * self[1][1];
313            adjugate_transpose[0][1] = c_10;
314            adjugate_transpose[1][1] = self[0][0] * self[2][2] - self[0][2] * self[2][0];
315            adjugate_transpose[2][1] = self[0][2] * self[1][0] - self[0][0] * self[1][2];
316            adjugate_transpose[0][2] = c_20;
317            adjugate_transpose[1][2] = self[0][1] * self[2][0] - self[0][0] * self[2][1];
318            adjugate_transpose[2][2] = self[0][0] * self[1][1] - self[0][1] * self[1][0];
319            (adjugate_transpose / determinant, determinant)
320        } else if D == 4 {
321            let mut adjugate_transpose = Self::zero();
322            let s0 = self[0][0] * self[1][1] - self[0][1] * self[1][0];
323            let s1 = self[0][0] * self[1][2] - self[0][2] * self[1][0];
324            let s2 = self[0][0] * self[1][3] - self[0][3] * self[1][0];
325            let s3 = self[0][1] * self[1][2] - self[0][2] * self[1][1];
326            let s4 = self[0][1] * self[1][3] - self[0][3] * self[1][1];
327            let s5 = self[0][2] * self[1][3] - self[0][3] * self[1][2];
328            let c5 = self[2][2] * self[3][3] - self[2][3] * self[3][2];
329            let c4 = self[2][1] * self[3][3] - self[2][3] * self[3][1];
330            let c3 = self[2][1] * self[3][2] - self[2][2] * self[3][1];
331            let c2 = self[2][0] * self[3][3] - self[2][3] * self[3][0];
332            let c1 = self[2][0] * self[3][2] - self[2][2] * self[3][0];
333            let c0 = self[2][0] * self[3][1] - self[2][1] * self[3][0];
334            let determinant = (s0 * c5 - s1 * c4 + s2 * c3 + s3 * c2 - s4 * c1 + s5 * c0).value();
335            adjugate_transpose[0][0] = self[1][1] * c5 - self[1][2] * c4 + self[1][3] * c3;
336            adjugate_transpose[1][0] = self[0][2] * c4 - self[0][1] * c5 - self[0][3] * c3;
337            adjugate_transpose[2][0] = self[3][1] * s5 - self[3][2] * s4 + self[3][3] * s3;
338            adjugate_transpose[3][0] = self[2][2] * s4 - self[2][1] * s5 - self[2][3] * s3;
339            adjugate_transpose[0][1] = self[1][2] * c2 - self[1][0] * c5 - self[1][3] * c1;
340            adjugate_transpose[1][1] = self[0][0] * c5 - self[0][2] * c2 + self[0][3] * c1;
341            adjugate_transpose[2][1] = self[3][2] * s2 - self[3][0] * s5 - self[3][3] * s1;
342            adjugate_transpose[3][1] = self[2][0] * s5 - self[2][2] * s2 + self[2][3] * s1;
343            adjugate_transpose[0][2] = self[1][0] * c4 - self[1][1] * c2 + self[1][3] * c0;
344            adjugate_transpose[1][2] = self[0][1] * c2 - self[0][0] * c4 - self[0][3] * c0;
345            adjugate_transpose[2][2] = self[3][0] * s4 - self[3][1] * s2 + self[3][3] * s0;
346            adjugate_transpose[3][2] = self[2][1] * s2 - self[2][0] * s4 - self[2][3] * s0;
347            adjugate_transpose[0][3] = self[1][1] * c1 - self[1][0] * c3 - self[1][2] * c0;
348            adjugate_transpose[1][3] = self[0][0] * c3 - self[0][1] * c1 + self[0][2] * c0;
349            adjugate_transpose[2][3] = self[3][1] * s1 - self[3][0] * s3 - self[3][2] * s0;
350            adjugate_transpose[3][3] = self[2][0] * s3 - self[2][1] * s1 + self[2][2] * s0;
351            (adjugate_transpose / determinant, determinant)
352        } else {
353            (self.inverse_transpose_core(), self.determinant_core())
354        }
355    }
356    fn lu_decomposition_core(&self) -> (Self, Self, Vec<usize>) {
357        let n = D;
358        let mut p: Vec<usize> = (0..n).collect();
359        let mut factor;
360        let mut lu = self.clone();
361        let mut max_row;
362        let mut max_val;
363        let mut pivot;
364        for i in 0..n {
365            max_row = i;
366            max_val = lu[max_row][i].abs();
367            for k in i + 1..n {
368                if lu[k][i].abs() > max_val {
369                    max_row = k;
370                    max_val = lu[max_row][i].abs();
371                }
372            }
373            if max_row != i {
374                lu.0.swap(i, max_row);
375                p.swap(i, max_row);
376            }
377            pivot = lu[i][i];
378            if pivot.abs() < ABS_TOL {
379                panic!("LU decomposition failed (zero pivot).")
380            }
381            for j in i + 1..n {
382                if lu[j][i] != 0.0 {
383                    lu[j][i] = lu[j][i] / pivot;
384                    factor = lu[j][i];
385                    for k in i + 1..n {
386                        let update = factor * lu[i][k];
387                        lu[j][k] -= update;
388                    }
389                }
390            }
391        }
392        let mut l = TensorRank2::identity();
393        for i in 0..D {
394            for j in 0..i {
395                l[i][j] = lu[i][j]
396            }
397        }
398        let mut u = TensorRank2::zero();
399        for i in 0..D {
400            for j in i..D {
401                u[i][j] = lu[i][j]
402            }
403        }
404        (l, u, p)
405    }
406    fn lu_decomposition_inverse_core(&self) -> (Self, Self, Vec<usize>) {
407        let (mut tensor_l, mut tensor_u, p) = self.lu_decomposition_core();
408        let mut sum: Quantity;
409        for i in 0..D {
410            tensor_l[i][i] = 1.0 / tensor_l[i][i];
411            for j in 0..i {
412                sum = Quantity::new(0.0);
413                for k in j..i {
414                    sum += tensor_l[i][k] * tensor_l[k][j];
415                }
416                tensor_l[i][j] = -sum * tensor_l[i][i];
417            }
418        }
419        for i in 0..D {
420            tensor_u[i][i] = 1.0 / tensor_u[i][i];
421            for j in 0..i {
422                sum = Quantity::new(0.0);
423                for k in j..i {
424                    sum += tensor_u[j][k] * tensor_u[k][i];
425                }
426                tensor_u[j][i] = -sum * tensor_u[i][i];
427            }
428        }
429        (tensor_l, tensor_u, p)
430    }
431}
conspire/math/tensor/rank_2/inverse/mod.rs

conspire/math/tensor/rank_2/inverse/
mod.rs
