Jpp
 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 <algorithm>
6 
7 #include "JMath/JMatrix3D.hh"
8 #include "JLang/JException.hh"
9 
10 
11 /**
12  * \author mdejong
13  */
14 
15 namespace JMATH {}
16 namespace JPP { using namespace JMATH; }
17 
18 namespace JMATH {
19 
21 
22 
23  /**
24  * 3 x 3 symmetric matrix
25  */
26  class JMatrix3S :
27  public JMatrix3D
28  {
29  public:
30  /**
31  * Default constructor.
32  */
34  JMatrix3D()
35  {}
36 
37 
38  /**
39  * Contructor.
40  *
41  * \param A matrix
42  */
43  JMatrix3S(const JMatrix3D& A) :
44  JMatrix3D(A)
45  {}
46 
47 
48  /**
49  * Contructor.
50  *
51  * \param __a00 (0,0)
52  * \param __a01 (0,1)
53  * \param __a02 (0,2)
54  * \param __a11 (1,1)
55  * \param __a12 (1,2)
56  * \param __a22 (2,2)
57  */
58  JMatrix3S(const double __a00, const double __a01, const double __a02,
59  /* */ const double __a11, const double __a12,
60  /* */ const double __a22) :
61  JMatrix3D(__a00, __a01, __a02,
62  __a01, __a11, __a12,
63  __a02, __a12, __a22)
64  {}
65 
66 
67  /**
68  * Invert matrix.
69  */
70  void invert()
71  {
72  using std::swap;
73 
74  // LDU decomposition
75 
76  int p0 = 0; // permute row 0
77  int p1 = 0; // permute row 1
78 
79  double val;
80 
81  val = a00;
82 
83  if (fabs(a10) > fabs(val)) {
84  p0 = 1;
85  val = a10;
86  }
87 
88  if (fabs(a20) > fabs(val)) {
89  p0 = 2;
90  val = a20;
91  }
92 
93  switch (p0) {
94 
95  case 1:
96  swap(a00,a10);
97  swap(a01,a11);
98  swap(a02,a12);
99  break;
100 
101  case 2:
102  swap(a00,a20);
103  swap(a01,a21);
104  swap(a02,a22);
105  break;
106  }
107 
108  if (val == 0) {
109  throw JDivisionByZero("LDU decomposition zero pivot");
110  }
111 
112  a10 /= val;
113  a11 -= a10 * a01;
114  a12 -= a10 * a02;
115 
116  a20 /= val;
117  a21 -= a20 * a01;
118  a22 -= a20 * a02;
119 
120  val = a11;
121 
122  if (fabs(a21) > fabs(val)) {
123  p1 = 2;
124  val = a21;
125  }
126 
127  switch (p1) {
128 
129  case 2:
130  swap(a10,a20);
131  swap(a11,a21);
132  swap(a12,a22);
133  break;
134  }
135 
136  if (val == 0) {
137  throw JDivisionByZero("LDU decomposition zero pivot");
138  }
139 
140  a21 /= val;
141  a22 -= a21 * a12;
142 
143  // invert D
144 
145  if (a22 == 0) {
146  throw JDivisionByZero("D matrix not invertable");
147  }
148 
149  a00 = 1.0 / a00;
150  a11 = 1.0 / a11;
151  a22 = 1.0 / a22;
152 
153  // invert U
154 
155  a01 *= -a00;
156  a01 *= a11;
157 
158  a02 *= -a00;
159  a02 -= a01 * a12;
160  a02 *= a22;
161 
162  a12 *= -a11;
163  a12 *= a22;
164 
165  // invert L
166 
167  a21 = -a21;
168  a20 = -a20;
169  a20 -= a21 * a10;
170  a10 = -a10;
171 
172  // U^-1 x L^-1
173 
174  a00 += a01 * a10 + a02 * a20;
175  a01 += a02 * a21;
176  a10 *= a11;
177  a10 += a12 * a20;
178  a11 += a12 * a21;
179  a20 *= a22;
180  a21 *= a22;
181 
182  switch (p1) {
183 
184  case 2:
185  swap(a01,a02);
186  swap(a11,a12);
187  swap(a21,a22);
188  break;
189  }
190 
191  switch (p0) {
192 
193  case 1:
194  swap(a00,a01);
195  swap(a10,a11);
196  swap(a20,a21);
197  break;
198 
199  case 2:
200  swap(a00,a02);
201  swap(a10,a12);
202  swap(a20,a22);
203  break;
204  }
205  }
206  };
207 }
208 
209 #endif
Exceptions.
TPaveText * p1
3 x 3 matrix
void invert()
Invert matrix.
Definition: JMatrix3S.hh:70
JMatrix3S()
Default constructor.
Definition: JMatrix3S.hh:33
3 x 3 symmetric matrix
Definition: JMatrix3S.hh:26
JMatrix3S(const double __a00, const double __a01, const double __a02, const double __a11, const double __a12, const double __a22)
Contructor.
Definition: JMatrix3S.hh:58
Exception for division by zero.
Definition: JException.hh:252
JMatrix3S(const JMatrix3D &A)
Contructor.
Definition: JMatrix3S.hh:43