Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JMatrixNS.hh
Go to the documentation of this file.
1 #ifndef __JMATH__JMATRIXNS__
2 #define __JMATH__JMATRIXNS__
3 
4 #include <vector>
5 
6 #include "JMath/JMatrixND.hh"
7 #include "JLang/JException.hh"
8 
9 
10 /**
11  * \author mdejong
12  */
13 
14 namespace JMATH {}
15 namespace JPP { using namespace JMATH; }
16 
17 namespace JMATH {
18 
21 
22 
23  /**
24  * N x N symmetric matrix
25  */
26  class JMatrixNS :
27  public JMatrixND
28  {
29  public:
30  /**
31  * Default constructor.
32  */
34  JMatrixND()
35  {}
36 
37 
38  /**
39  * Constructor (identity matrix).
40  *
41  * \param size dimension
42  */
43  JMatrixNS(const unsigned int size) :
44  JMatrixND(size)
45  {}
46 
47 
48  /**
49  * Contructor.
50  *
51  * \param A matrix
52  */
53  JMatrixNS(const JMatrixND& A) :
54  JMatrixND(A)
55  {}
56 
57 
58  /**
59  * Invert matrix.
60  */
61  void invert()
62  {
63  // LDU decomposition
64 
65  unsigned int i, j, k;
66  double val, del;
67  row_type row, root;
68  col_type col, cursor;
69 
71 
72  for (root = this->begin(), k = 0; root != this->end(); ++root, ++k) {
73 
74  row_type pivot;
75 
76  val = (*root)[k];
77  pivot = root;
78 
79  for (row = root; ++row != this->end(); ) {
80  if (fabs((*row)[k]) > fabs(val)) {
81  val = (*row)[k];
82  pivot = row;
83  }
84  }
85 
86  if (val == 0) {
87  throw JDivisionByZero("LDU decomposition zero pivot");
88  }
89 
90  if (pivot != root) {
91 
92  std::iter_swap(root, pivot);
93 
94  P.push_back(std::make_pair(std::distance(this->begin(),root),
95  std::distance(this->begin(),pivot)));
96  }
97 
98  // decompose in LDU
99 
100  for (row = root; ++row != this->end(); ) {
101 
102  del = (*row)[k] /= val;
103 
104  for (col = row->begin() + k, cursor = root->begin() + k; ++col != row->end(); ) {
105  (*col) -= del * (*(++cursor));
106  }
107  }
108  }
109 
110 
111  // invert D
112 
113  for (k = 0; k != this->size(); ++k) {
114 
115  double& a = this->at(k,k);
116 
117  if (a == 0) {
118  throw JDivisionByZero("LDU decomposition zero pivot");
119  }
120 
121  a = 1.0/a;
122  }
123 
124 
125  // invert U
126 
127  for (root = this->begin(), k = 0; root != this->end(); ++root, ++k) {
128 
129  col_type pivot = root->begin() + k;
130 
131  for (col = pivot, j = k + 1; ++col != root->end(); ++j) {
132 
133  (*col) *= -(*pivot);
134 
135  for (cursor = pivot, row = root; ++cursor != col; ) {
136  (*col) -= (*cursor) * (*(++row))[j];
137  }
138 
139  (*col) *= this->at(j,j);
140  }
141  }
142 
143 
144  // invert L
145 
146  for (i = size(); i != 0; ) {
147 
148  root = this->begin() + --i;
149 
150  for (j = i; j != 0; ) {
151 
152  double& a = (*root)[--j];
153 
154  a = -a;
155 
156  for (k = j + 1, cursor = root->begin() + k, row = this->begin() + k; k != i; ++cursor, ++row, ++k) {
157  a -= (*cursor) * (*row)[j];
158  }
159  }
160  }
161 
162 
163  // U^-1 x D^-1 x L^-1
164 
165  for (root = this->begin(), i = 0; root != this->end(); ++root, ++i) {
166 
167  val = (*root)[i];
168 
169  for (col = root->begin(), j = 0; j != i; ++col, ++j) {
170 
171  (*col) *= val;
172 
173  for (row = root, cursor = root->begin() + i; ++cursor != root->end(); ) {
174  (*col) += (*cursor) * (*(++row))[j];
175  }
176  }
177 
178  for (j = i, col = root->begin() + j; j != this->size() - 1; ++col, ++j) {
179  for (k = j + 1, row = this->begin() + k, cursor = root->begin() + k; k != this->size(); ++k, ++cursor, ++row) {
180  (*col) += (*(cursor)) * (*(row))[j];
181  }
182  }
183  }
184 
185 
186  for (JPermutationMatrix::reverse_iterator p = P.rbegin(); p != P.rend(); ++p) {
187  for (row = this->begin(); row != this->end(); ++row) {
188  std::swap((*row)[p->first], (*row)[p->second]);
189  }
190  }
191  }
192 
193 
194  /**
195  * Update inverted matrix at given diagonal element.
196  *
197  * If A is the original matrix and A' is such that for each (i,j) != (k,k):
198  * <pre>
199  * A'[i][j] = A[i][j]
200  * A'[k][k] = A[k][k] + extra
201  * </pre>
202  * then JMatrixNS::update(k, extra) will modify the already inverted matrix
203  * of A such that it will be the inverse of A'.
204  *
205  * \param k index of diagonal element
206  * \param extra value to add
207  */
208  void update(const unsigned int k, const double extra)
209  {
210  if (k >= this->size()) {
211  throw JIndexOutOfRange("JMatrixNS::invert(): index out of range.");
212  }
213 
214  unsigned int i, j;
215  row_type row, root;
216  col_type col;
217 
218  const double del = this->at(k,k);
219  const double din = 1.0 / (1.0 + extra * del);
220 
221  root = this->begin() + k;
222 
223  for (i = 0, row = this->begin(); row != this->end(); ++i, ++row) {
224 
225  if (i != k) {
226 
227  const double val = (*row)[k];
228 
229  for (j = 0, col = row->begin(); col != row->end(); ++j, ++col) {
230 
231  if (j != k) {
232  (*col) -= val * (*root)[j] * extra * din;
233  }
234  }
235 
236  (*row)[k] *= din;
237  }
238  }
239 
240  for (j = 0, col = root->begin(); col != root->end(); ++j, ++col) {
241 
242  if (j != k) {
243  (*col) *= din;
244  }
245  }
246 
247  (*root)[k] = del * din;
248  }
249  };
250 }
251 
252 #endif
double at(size_t row, size_t col) const
Get matrix element.
Definition: JMatrixND.hh:120
Exceptions.
std::vector< double >::iterator col_type
Definition: JMatrixND.hh:50
N x N symmetric matrix.
Definition: JMatrixNS.hh:26
JMatrixNS(const JMatrixND &A)
Contructor.
Definition: JMatrixNS.hh:53
void update(const unsigned int k, const double extra)
Update inverted matrix at given diagonal element.
Definition: JMatrixNS.hh:208
JMatrixNS(const unsigned int size)
Constructor (identity matrix).
Definition: JMatrixNS.hh:43
matrix_type::iterator row_type
Definition: JMatrixND.hh:49
N x N matrix.
Definition: JMatrixND.hh:63
Exception for division by zero.
Definition: JException.hh:252
JMatrixNS()
Default constructor.
Definition: JMatrixNS.hh:33
Exception for accessing an index in a collection that is outside of its range.
Definition: JException.hh:90
void invert()
Invert matrix.
Definition: JMatrixNS.hh:61