matext4cgmath/
eigens.rs

1use crate::*;
2
3impl<F: BaseFloat> EigenValues for Matrix2<F> {
4    type EigenValues = [Complex<F>; 2];
5    fn eigenvalues(self) -> [Complex<F>; 2] {
6        solver::solve_quadratic(-self.trace(), self.determinant())
7    }
8}
9
10impl<F: BaseFloat> EigenValues for Matrix3<F> {
11    type EigenValues = [Complex<F>; 3];
12    fn eigenvalues(self) -> [Complex<F>; 3] {
13        let poly2 = self[0][0] * self[1][1] + self[0][0] * self[2][2] + self[1][1] * self[2][2]
14            - self[0][1] * self[1][0]
15            - self[1][2] * self[2][1]
16            - self[0][2] * self[2][0];
17        solver::solve_cubic(-self.trace(), poly2, -self.determinant())
18    }
19}
20
21impl<F: BaseFloat> EigenValues for Matrix4<F> {
22    type EigenValues = [Complex<F>; 4];
23    fn eigenvalues(self) -> [Complex<F>; 4] {
24        let poly2 = self[0][0] * self[1][1]
25            + self[0][0] * self[2][2]
26            + self[0][0] * self[3][3]
27            + self[1][1] * self[2][2]
28            + self[1][1] * self[3][3]
29            + self[2][2] * self[3][3]
30            - self[0][1] * self[1][0]
31            - self[0][2] * self[2][0]
32            - self[0][3] * self[3][0]
33            - self[1][2] * self[2][1]
34            - self[1][3] * self[3][1]
35            - self[2][3] * self[3][2];
36        let poly3 = (self[0][0] * self[1][1] + self[0][0] * self[2][2] + self[1][1] * self[2][2]
37            - self[0][1] * self[1][0]
38            - self[0][2] * self[2][0]
39            - self[1][2] * self[2][1])
40            * self[3][3]
41            + (self[2][0] * self[3][2] - self[2][2] * self[3][0] + self[1][0] * self[3][1]
42                - self[1][1] * self[3][0])
43                * self[0][3]
44            + (self[2][1] * self[3][2] - self[2][2] * self[3][1] + self[3][0] * self[0][1]
45                - self[3][1] * self[0][0])
46                * self[1][3]
47            + (self[3][1] * self[1][2] - self[3][2] * self[1][1] + self[3][0] * self[0][2]
48                - self[3][2] * self[0][0])
49                * self[2][3]
50            + Matrix3::from_cols(self[0].truncate(), self[1].truncate(), self[2].truncate())
51                .determinant();
52        solver::solve_quartic(-self.trace(), poly2, -poly3, self.determinant())
53    }
54}
55
56impl<F: BaseFloat> OperatorNorm for Matrix2<F> {
57    #[inline]
58    fn norm_l1(self) -> F {
59        F::max(
60            F::abs(self[0][0]) + F::abs(self[0][1]),
61            F::abs(self[1][0]) + F::abs(self[1][1]),
62        )
63    }
64    #[inline]
65    fn norm_linf(self) -> F {
66        F::max(
67            F::abs(self[0][0]) + F::abs(self[1][0]),
68            F::abs(self[0][1]) + F::abs(self[1][1]),
69        )
70    }
71    #[inline]
72    fn norm_l2(self) -> F {
73        let eigens = (self.transpose() * self).eigenvalues();
74        F::sqrt(F::max(eigens[0].re, eigens[1].re))
75    }
76}
77
78impl<F: BaseFloat> OperatorNorm for Matrix3<F> {
79    #[inline]
80    fn norm_l1(self) -> F {
81        F::max(
82            F::max(
83                F::abs(self[0][0]) + F::abs(self[0][1]) + F::abs(self[0][2]),
84                F::abs(self[1][0]) + F::abs(self[1][1]) + F::abs(self[1][2]),
85            ),
86            F::abs(self[2][0]) + F::abs(self[2][1]) + F::abs(self[2][2]),
87        )
88    }
89    #[inline]
90    fn norm_linf(self) -> F {
91        F::max(
92            F::max(
93                F::abs(self[0][0]) + F::abs(self[1][0]) + F::abs(self[2][0]),
94                F::abs(self[0][1]) + F::abs(self[1][1]) + F::abs(self[2][1]),
95            ),
96            F::abs(self[0][2]) + F::abs(self[1][2]) + F::abs(self[2][2]),
97        )
98    }
99    #[inline]
100    fn norm_l2(self) -> F {
101        let eigens = (self.transpose() * self).eigenvalues();
102        F::sqrt(F::max(F::max(eigens[0].re, eigens[1].re), eigens[2].re))
103    }
104}
105
106impl<F: BaseFloat> OperatorNorm for Matrix4<F> {
107    #[inline]
108    fn norm_l1(self) -> F {
109        F::max(
110            F::max(
111                F::max(
112                    F::abs(self[0][0])
113                        + F::abs(self[0][1])
114                        + F::abs(self[0][2])
115                        + F::abs(self[0][3]),
116                    F::abs(self[1][0])
117                        + F::abs(self[1][1])
118                        + F::abs(self[1][2])
119                        + F::abs(self[1][3]),
120                ),
121                F::abs(self[2][0]) + F::abs(self[2][1]) + F::abs(self[2][2]) + F::abs(self[2][3]),
122            ),
123            F::abs(self[3][0]) + F::abs(self[3][1]) + F::abs(self[3][2]) + F::abs(self[3][3]),
124        )
125    }
126    #[inline]
127    fn norm_linf(self) -> F {
128        F::max(
129            F::max(
130                F::max(
131                    F::abs(self[0][0])
132                        + F::abs(self[1][0])
133                        + F::abs(self[2][0])
134                        + F::abs(self[3][0]),
135                    F::abs(self[0][1])
136                        + F::abs(self[1][1])
137                        + F::abs(self[2][1])
138                        + F::abs(self[3][1]),
139                ),
140                F::abs(self[0][2]) + F::abs(self[1][2]) + F::abs(self[2][2]) + F::abs(self[3][2]),
141            ),
142            F::abs(self[0][3]) + F::abs(self[1][3]) + F::abs(self[2][3]) + F::abs(self[3][3]),
143        )
144    }
145    #[inline]
146    fn norm_l2(self) -> F {
147        let eigens = (self.transpose() * self).eigenvalues();
148        F::sqrt(F::max(
149            F::max(F::max(eigens[0].re, eigens[1].re), eigens[2].re),
150            eigens[3].re,
151        ))
152    }
153}
matext4cgmath/eigens.rs

matext4cgmath/
eigens.rs
