Jpp  18.0.0-rc.3
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JMatrix3S.hh
Go to the documentation of this file.
1 #ifndef __JMATRIX3S__
2 #define __JMATRIX3S__
3 
4 #include <cmath>
5 #include <array>
6 #include <algorithm>
7 
8 #include "JLang/JException.hh"
9 
10 #include "JMath/JConstants.hh"
11 #include "JMath/JMatrix3D.hh"
12 
13 
14 /**
15  * \author mdejong
16  */
17 
18 namespace JMATH {}
19 namespace JPP { using namespace JMATH; }
20 
21 namespace JMATH {
22 
24 
25 
26  /**
27  * 3 x 3 symmetric matrix
28  */
29  class JMatrix3S :
30  public JMatrix3D
31  {
32  public:
33  /**
34  * Type definition of Eigen vakues.
35  */
36  typedef std::array<double, 3> eigen_values;
37 
38 
39  /**
40  * Default constructor.
41  */
43  JMatrix3D()
44  {}
45 
46 
47  /**
48  * Contructor.
49  *
50  * \param A matrix
51  */
52  JMatrix3S(const JMatrix3D& A) :
53  JMatrix3D(A)
54  {}
55 
56 
57  /**
58  * Contructor.\n
59  * The upper triangle is internally set.
60  *
61  * \param __a00 (0,0)
62  * \param __a10 (1,0)
63  * \param __a11 (1,1)
64  * \param __a20 (2,0)
65  * \param __a21 (2,1)
66  * \param __a22 (2,2)
67  */
68  JMatrix3S(const double __a00,
69  const double __a10, const double __a11,
70  const double __a20, const double __a21, const double __a22) :
71  JMatrix3D(__a00, __a10, __a20,
72  __a10, __a11, __a21,
73  __a20, __a21, __a22)
74  {}
75 
76 
77  /**
78  * Invert matrix.
79  */
80  void invert()
81  {
82  using std::swap;
83 
84  // LDU decomposition
85 
86  int p0 = 0; // permute row 0
87  int p1 = 0; // permute row 1
88 
89  double val;
90 
91  val = a00;
92 
93  if (fabs(a10) > fabs(val)) {
94  p0 = 1;
95  val = a10;
96  }
97 
98  if (fabs(a20) > fabs(val)) {
99  p0 = 2;
100  val = a20;
101  }
102 
103  switch (p0) {
104 
105  case 1:
106  swap(a00,a10);
107  swap(a01,a11);
108  swap(a02,a12);
109  break;
110 
111  case 2:
112  swap(a00,a20);
113  swap(a01,a21);
114  swap(a02,a22);
115  break;
116  }
117 
118  if (val == 0) {
119  throw JDivisionByZero("LDU decomposition zero pivot");
120  }
121 
122  a10 /= val;
123  a11 -= a10 * a01;
124  a12 -= a10 * a02;
125 
126  a20 /= val;
127  a21 -= a20 * a01;
128  a22 -= a20 * a02;
129 
130  val = a11;
131 
132  if (fabs(a21) > fabs(val)) {
133  p1 = 2;
134  val = a21;
135  }
136 
137  switch (p1) {
138 
139  case 2:
140  swap(a10,a20);
141  swap(a11,a21);
142  swap(a12,a22);
143  break;
144  }
145 
146  if (val == 0) {
147  throw JDivisionByZero("LDU decomposition zero pivot");
148  }
149 
150  a21 /= val;
151  a22 -= a21 * a12;
152 
153  // invert D
154 
155  if (a22 == 0) {
156  throw JDivisionByZero("D matrix not invertable");
157  }
158 
159  a00 = 1.0 / a00;
160  a11 = 1.0 / a11;
161  a22 = 1.0 / a22;
162 
163  // invert U
164 
165  a01 *= -a00;
166  a01 *= a11;
167 
168  a02 *= -a00;
169  a02 -= a01 * a12;
170  a02 *= a22;
171 
172  a12 *= -a11;
173  a12 *= a22;
174 
175  // invert L
176 
177  a21 = -a21;
178  a20 = -a20;
179  a20 -= a21 * a10;
180  a10 = -a10;
181 
182  // U^-1 x L^-1
183 
184  a00 += a01 * a10 + a02 * a20;
185  a01 += a02 * a21;
186  a10 *= a11;
187  a10 += a12 * a20;
188  a11 += a12 * a21;
189  a20 *= a22;
190  a21 *= a22;
191 
192  switch (p1) {
193 
194  case 2:
195  swap(a01,a02);
196  swap(a11,a12);
197  swap(a21,a22);
198  break;
199  }
200 
201  switch (p0) {
202 
203  case 1:
204  swap(a00,a01);
205  swap(a10,a11);
206  swap(a20,a21);
207  break;
208 
209  case 2:
210  swap(a00,a02);
211  swap(a10,a12);
212  swap(a20,a22);
213  break;
214  }
215  }
216 
217 
218  /**
219  * Get eigen values.\n
220  * The eigen values sorted in decreasing order of absolute values.\n
221  * Algorithm taken from "Eigenvalues of a symmetric 3x3 matrix"\n
222  * by Oliver K. Smith; see <a href="https://dl.acm.org/doi/10.1145/355578.366316">reference</a>.
223  *
224  * \param epsilon precision
225  * \return eigen values
226  */
227  eigen_values getEigenValues(const double epsilon = 1e-6) const
228  {
229  using namespace std;
230 
231  eigen_values eigenvalues;
232 
233  const double p1 = a01*a01 + a02*a02 + a12*a12;
234 
235  if (fabs(p1 - 1.0) > epsilon) {
236 
237  const double q = (a00 + a11 + a22) / 3.0;
238  const double p2 = (a00 - q)*(a00 - q) + (a11 - q)*(a11 - q) + (a22 - q)*(a22 - q) + 2*p1;
239  const double p = sqrt(p2 / 6.0);
240 
241  JMatrix3S B(*this);
242 
243  B.sub(q * JMatrix3D::getIdentity());
244  B.div(p);
245 
246  const double r = B.getDeterminant() / 2.0;
247  const double phi = (r < -1.0 ? PI / 3.0 : (r > 1.0 ? 0.0 : acos(r) / 3.0));
248 
249  eigenvalues[0] = q + 2*p*cos(phi);
250  eigenvalues[2] = q + 2*p*cos(phi + 2*PI/3.0);
251  eigenvalues[1] = 3 * q - eigenvalues[0] - eigenvalues[2];
252 
253  } else {
254 
255  eigenvalues = eigen_values{ a00, a11, a12 };
256  }
257 
258  // sort absolute values in descending order
259 
260  if (fabs(eigenvalues[0]) < fabs(eigenvalues[1])) {
261  swap(eigenvalues[0], eigenvalues[1]);
262  }
263 
264  if (fabs(eigenvalues[1]) < fabs(eigenvalues[2])) {
265 
266  swap(eigenvalues[1], eigenvalues[2]);
267 
268  if (fabs(eigenvalues[0]) < fabs(eigenvalues[1])) {
269  swap(eigenvalues[0], eigenvalues[1]);
270  }
271  }
272 
273  return eigenvalues;
274  }
275  };
276 }
277 
278 #endif
eigen_values getEigenValues(const double epsilon=1e-6) const
Get eigen values.
Definition: JMatrix3S.hh:227
Exceptions.
JMatrix3S(const double __a00, const double __a10, const double __a11, const double __a20, const double __a21, const double __a22)
Contructor.
Definition: JMatrix3S.hh:68
TPaveText * p1
static const JMatrix3D & getIdentity()
Get reference to unique instance of this class object.
JMatrix3D & sub(const JMatrix3D &A)
Matrix subtraction.
3 x 3 matrix
void invert()
Invert matrix.
Definition: JMatrix3S.hh:80
JMatrix3S()
Default constructor.
Definition: JMatrix3S.hh:42
data_type r[M+1]
Definition: JPolint.hh:779
3 x 3 symmetric matrix
Definition: JMatrix3S.hh:29
Mathematical constants.
static const double PI
Mathematical constants.
std::array< double, 3 > eigen_values
Type definition of Eigen vakues.
Definition: JMatrix3S.hh:36
p2
Definition: module-Z:fit.sh:74
JMatrix3D & div(const double factor)
Scale matrix.
Exception for division by zero.
Definition: JException.hh:270
double getDeterminant() const
Get determinant of matrix.
source $JPP_DIR setenv csh $JPP_DIR &dev null eval JShellParser o a A
JMatrix3S(const JMatrix3D &A)
Contructor.
Definition: JMatrix3S.hh:52
const double epsilon
Definition: JQuadrature.cc:21