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

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

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