Jpp
Public Types | Public Member Functions | List of all members
JMATH::JMatrixNS Class Reference

N x N symmetric matrix. More...

#include <JMatrixNS.hh>

Inheritance diagram for JMATH::JMatrixNS:
JMATH::JMatrixND JMATH::JMatrixND_t JMATH::JMath< JMatrixND > JLANG::JEquals< JFirst_t, JSecond_t > std::vector< std::vector< double > > JFIT::JMatrixNZ

Public Types

typedef std::vector< std::vector< double > > matrix_type
 
typedef matrix_type::const_iterator const_row_type
 
typedef std::vector< double >::const_iterator const_col_type
 
typedef matrix_type::iterator row_type
 
typedef std::vector< double >::iterator col_type
 

Public Member Functions

 JMatrixNS ()
 Default constructor. More...
 
 JMatrixNS (const unsigned int size)
 Constructor (identity matrix). More...
 
 JMatrixNS (const JMatrixND &A)
 Contructor. More...
 
void invert ()
 Invert matrix. More...
 
void update (const unsigned int k, const double extra)
 Update inverted matrix at given diagonal element. More...
 
void resize (const size_t size)
 Resize matrix. More...
 
double at (size_t row, size_t col) const
 Get matrix element. More...
 
double & at (size_t row, size_t col)
 Get matrix element. More...
 
double operator() (size_t row, size_t col) const
 Get matrix element. More...
 
double & operator() (size_t row, size_t col)
 Get matrix element. More...
 
JMatrixND & reset ()
 Set matrix to the null matrix. More...
 
JMatrixND & setIdentity ()
 Set to identity matrix. More...
 
JMatrixND & transpose ()
 Transpose. More...
 
JMatrixND & negate ()
 Negate matrix. More...
 
JMatrixND & add (const JMatrixND &A)
 Matrix addition. More...
 
JMatrixND & sub (const JMatrixND &A)
 Matrix subtraction. More...
 
JMatrixND & mul (const double factor)
 Scale matrix. More...
 
JMatrixND & mul (const JMatrixND &A, const JMatrixND &B)
 Matrix multiplication. More...
 
JMatrixND & mul (const JSecond_t &object)
 Multiply with object. More...
 
JMatrixND & div (const double factor)
 Scale matrix. More...
 
bool equals (const JMatrixND &A, const double eps=std::numeric_limits< double >::min()) const
 Equality. More...
 
bool isIdentity (const double eps=std::numeric_limits< double >::min()) const
 Test identity. More...
 
double getDeterminant () const
 Get determinant of matrix. More...
 
bool isPositiveDefinite () const
 Test positive definiteness. More...
 
const bool isPositiveSemiDefinite () const
 Test positive semi-definiteness. More...
 

Detailed Description

N x N symmetric matrix.

Definition at line 26 of file JMatrixNS.hh.

Member Typedef Documentation

◆ matrix_type

typedef std::vector< std::vector<double> > JMATH::JMatrixND_t::matrix_type
inherited

Definition at line 44 of file JMatrixND.hh.

◆ const_row_type

typedef matrix_type::const_iterator JMATH::JMatrixND_t::const_row_type
inherited

Definition at line 46 of file JMatrixND.hh.

◆ const_col_type

typedef std::vector<double>::const_iterator JMATH::JMatrixND_t::const_col_type
inherited

Definition at line 47 of file JMatrixND.hh.

◆ row_type

typedef matrix_type::iterator JMATH::JMatrixND_t::row_type
inherited

Definition at line 49 of file JMatrixND.hh.

◆ col_type

typedef std::vector<double>::iterator JMATH::JMatrixND_t::col_type
inherited

Definition at line 50 of file JMatrixND.hh.

Constructor & Destructor Documentation

◆ JMatrixNS() [1/3]

JMATH::JMatrixNS::JMatrixNS ( )
inline

Default constructor.

Definition at line 33 of file JMatrixNS.hh.

33  :
34  JMatrixND()
35  {}

◆ JMatrixNS() [2/3]

JMATH::JMatrixNS::JMatrixNS ( const unsigned int  size)
inline

Constructor (identity matrix).

Parameters
sizedimension

Definition at line 43 of file JMatrixNS.hh.

43  :
44  JMatrixND(size)
45  {}

◆ JMatrixNS() [3/3]

JMATH::JMatrixNS::JMatrixNS ( const JMatrixND &  A)
inline

Contructor.

Parameters
Amatrix

Definition at line 53 of file JMatrixNS.hh.

53  :
54  JMatrixND(A)
55  {}

Member Function Documentation

◆ invert()

void JMATH::JMatrixNS::invert ( )
inline

Invert matrix.

Definition at line 61 of file JMatrixNS.hh.

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 
70  JPermutationMatrix P;
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  }

◆ update()

void JMATH::JMatrixNS::update ( const unsigned int  k,
const double  extra 
)
inline

Update inverted matrix at given diagonal element.

If A is the original matrix and A' is such that for each (i,j) != (k,k): 

      A'[i][j] = A[i][j]          
      A'[k][k] = A[k][k] + extra

then JMatrixNS::update(k, extra) will modify the already inverted matrix of A such that it will be the inverse of A'.

Parameters
kindex of diagonal element
extravalue to add

Definition at line 208 of file JMatrixNS.hh.

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  }

◆ resize()

void JMATH::JMatrixND::resize ( const size_t  size)
inlineinherited

Resize matrix.

Note that this method resets the matrix.

Parameters
sizedimension

Definition at line 98 of file JMatrixND.hh.

99  {
100  matrix_type::resize(size);
101 
102  for (row_type row = this->begin(); row != this->end(); ++row) {
103 
104  row->resize(size);
105 
106  for (col_type col = row->begin(); col != row->end(); ++col) {
107  *col = 0.0;
108  }
109  }
110  }

◆ at() [1/2]

double JMATH::JMatrixND::at ( size_t  row,
size_t  col 
) const
inlineinherited

Get matrix element.

Parameters
rowrow number
colcolumn number
Returns
matrix element at (row,col)

Definition at line 120 of file JMatrixND.hh.

121  {
122  return (*this)[row][col];
123  }

◆ at() [2/2]

double& JMATH::JMatrixND::at ( size_t  row,
size_t  col 
)
inlineinherited

Get matrix element.

Parameters
rowrow number
colcolumn number
Returns
matrix element at (row,col)

Definition at line 133 of file JMatrixND.hh.

134  {
135  return (*this)[row][col];
136  }

◆ operator()() [1/2]

double JMATH::JMatrixND::operator() ( size_t  row,
size_t  col 
) const
inlineinherited

Get matrix element.

Parameters
rowrow number
colcolumn number
Returns
matrix element at (row,col)

Definition at line 146 of file JMatrixND.hh.

147  {
148  return at(row, col);
149  }

◆ operator()() [2/2]

double& JMATH::JMatrixND::operator() ( size_t  row,
size_t  col 
)
inlineinherited

Get matrix element.

Parameters
rowrow number
colcolumn number
Returns
matrix element at (row,col)

Definition at line 159 of file JMatrixND.hh.

160  {
161  return at(row, col);
162  }

◆ reset()

JMatrixND& JMATH::JMatrixND::reset ( )
inlineinherited

Set matrix to the null matrix.

Returns
this matrix

Definition at line 170 of file JMatrixND.hh.

171  {
172  for (row_type row = this->begin(); row != this->end(); ++row) {
173  for (col_type col = row->begin(); col != row->end(); ++col) {
174  *col = 0.0;
175  }
176  }
177 
178  return *this;
179  }

◆ setIdentity()

JMatrixND& JMATH::JMatrixND::setIdentity ( )
inlineinherited

Set to identity matrix.

Returns
this matrix

Definition at line 187 of file JMatrixND.hh.

188  {
189  for (size_t i = 0; i != this->size(); ++i) {
190 
191  at(i,i) = 1.0;
192 
193  for (size_t j = 0; j != i; ++j) {
194  at(i,j) = at(j,i) = 0.0;
195  }
196  }
197 
198  return *this;
199  }

◆ transpose()

JMatrixND& JMATH::JMatrixND::transpose ( )
inlineinherited

Transpose.

Returns
this matrix

Definition at line 207 of file JMatrixND.hh.

208  {
209  for (size_t i = 0; i != this->size(); ++i) {
210  for (size_t j = 0; j != i; ++j) {
211  std::swap(at(i,j), at(j,i));
212  }
213  }
214 
215  return *this;
216  }

◆ negate()

JMatrixND& JMATH::JMatrixND::negate ( )
inlineinherited

Negate matrix.

Returns
this matrix

Definition at line 224 of file JMatrixND.hh.

225  {
226  for (row_type row = this->begin(); row != this->end(); ++row) {
227  for (col_type col = row->begin(); col != row->end(); ++col) {
228  *col = -*col;
229  }
230  }
231 
232  return *this;
233  }

◆ add()

JMatrixND& JMATH::JMatrixND::add ( const JMatrixND &  A)
inlineinherited

Matrix addition.

Parameters
Amatrix
Returns
this matrix

Definition at line 242 of file JMatrixND.hh.

243  {
244  if (this->size() == A.size()) {
245 
246  for (size_t i = 0; i != this->size(); ++i) {
247  for (size_t j = 0; j != this->size(); ++j) {
248  at(i,j) += A.at(i,j);
249  }
250  }
251 
252  } else {
253  THROW(JException, "JMatrixND::add() inconsistent matrix dimensions " << this->size() << ' ' << A.size());
254  }
255 
256  return *this;
257  }

◆ sub()

JMatrixND& JMATH::JMatrixND::sub ( const JMatrixND &  A)
inlineinherited

Matrix subtraction.

Parameters
Amatrix
Returns
this matrix

Definition at line 266 of file JMatrixND.hh.

267  {
268  if (this->size() == A.size()) {
269 
270  for (size_t i = 0; i != this->size(); ++i) {
271  for (size_t j = 0; j != this->size(); ++j) {
272  at(i,j) -= A.at(i,j);
273  }
274  }
275 
276  } else {
277  THROW(JException, "JMatrixND::sub() inconsistent matrix dimensions " << this->size() << ' ' << A.size());
278  }
279 
280  return *this;
281  }

◆ mul() [1/3]

JMatrixND& JMATH::JMatrixND::mul ( const double  factor)
inlineinherited

Scale matrix.

Parameters
factorfactor
Returns
this matrix

Definition at line 290 of file JMatrixND.hh.

291  {
292  for (row_type row = this->begin(); row != this->end(); ++row) {
293  for (col_type col = row->begin(); col != row->end(); ++col) {
294  *col *= factor;
295  }
296  }
297 
298  return *this;
299  }

◆ mul() [2/3]

JMatrixND& JMATH::JMatrixND::mul ( const JMatrixND &  A,
const JMatrixND &  B 
)
inlineinherited

Matrix multiplication.

Parameters
Amatrix
Bmatrix
Returns
this matrix

Definition at line 327 of file JMatrixND.hh.

329  {
330  if (A.size() == B.size()) {
331 
332  this->resize(A.size());
333 
334  for (size_t i = 0; i != A.size(); ++i) {
335  for (size_t j = 0; j != A.size(); ++j) {
336 
337  this->at(i,j) = 0.0;
338 
339  for (size_t k = 0; k != A.size(); ++k) {
340  this->at(i,j) += A[i][k] * B[k][j];
341  }
342  }
343  }
344 
345  } else {
346 
347  THROW(JException, "JMatrixND::mul() inconsistent matrix dimensions " << A.size() << ' ' << B.size());
348  }
349 
350  return *this;
351  }

◆ mul() [3/3]

JMatrixND & JMATH::JMath< JMatrixND , JSecond_t >::mul ( const JSecond_t &  object)
inlineinherited

Multiply with object.

Parameters
objectobject
Returns
result object

Definition at line 273 of file JMath.hh.

274  {
275  return static_cast<JFirst_t&>(*this) = JCalculator<JFirst_t>::calculator.mul(static_cast<const JFirst_t&>(*this), object);
276  }

◆ div()

JMatrixND& JMATH::JMatrixND::div ( const double  factor)
inlineinherited

Scale matrix.

Parameters
factorfactor
Returns
this matrix

Definition at line 308 of file JMatrixND.hh.

309  {
310  for (row_type row = this->begin(); row != this->end(); ++row) {
311  for (col_type col = row->begin(); col != row->end(); ++col) {
312  *col /= factor;
313  }
314  }
315 
316  return *this;
317  }

◆ equals()

bool JMATH::JMatrixND::equals ( const JMatrixND &  A,
const double  eps = std::numeric_limits<double>::min() 
) const
inlineinherited

Equality.

Parameters
Amatrix
epsnumerical precision
Returns
true if matrices identical; else false

Definition at line 361 of file JMatrixND.hh.

363  {
364  if (this->size() == A.size()) {
365 
366  for (const_row_type row1 = this->begin(), row2 = A.begin(); row1 != this->end(); ++row1, ++row2) {
367  for (const_col_type col1 = row1->begin(), col2 = row2->begin(); col1 != row1->end(); ++col1, ++col2) {
368  if (fabs(*col1 - *col2) > eps) {
369  return false;
370  }
371  }
372  }
373 
374  return true;
375 
376  } else {
377  THROW(JException, "JMatrixND::equals() inconsistent matrix dimensions " << this->size() << ' ' << A.size());
378  }
379  }

◆ isIdentity()

bool JMATH::JMatrixND::isIdentity ( const double  eps = std::numeric_limits<double>::min()) const
inlineinherited

Test identity.

Parameters
epsnumerical precision
Returns
true if identity matrix; else false

Definition at line 388 of file JMatrixND.hh.

389  {
390  for (size_t i = 0; i != this->size(); ++i) {
391 
392  if (fabs(1.0 - at(i,i)) > eps) {
393  return false;
394  };
395 
396  for (size_t j = 0; j != i; ++j) {
397  if (fabs(at(i,j)) > eps || fabs(at(j,i)) > eps) {
398  return false;
399  };
400  }
401  }
402 
403  return true;
404  }

◆ getDeterminant()

double JMATH::JMatrixND::getDeterminant ( ) const
inlineinherited

Get determinant of matrix.

Returns
determinant of matrix

Definition at line 412 of file JMatrixND.hh.

413  {
414  JMatrixHelper& A = JMatrixHelper::getInstance();
415 
416  A.decompose(*this);
417 
418  double det = 1.0;
419 
420  for (size_t i = 0; i != A.size(); ++i) {
421  det *= A[i][i];
422  }
423 
424  return det * A.getSign();
425  }

◆ isPositiveDefinite()

bool JMATH::JMatrixND::isPositiveDefinite ( ) const
inlineinherited

Test positive definiteness.

Returns
true if positive definite; else false

Definition at line 433 of file JMatrixND.hh.

434  {
435  JMatrixHelper& A = JMatrixHelper::getInstance();
436 
437  A.decompose(*this);
438 
439  for (size_t i = 0; i != A.size(); ++i) {
440  if (A[i][i] <= 0.0) {
441  return false;
442  }
443  }
444 
445  return true;
446  }

◆ isPositiveSemiDefinite()

const bool JMATH::JMatrixND::isPositiveSemiDefinite ( ) const
inlineinherited

Test positive semi-definiteness.

Returns
true if positive semi-definite; else false

Definition at line 454 of file JMatrixND.hh.

455  {
456  JMatrixHelper& A = JMatrixHelper::getInstance();
457 
458  A.decompose(*this);
459 
460  for (size_t i = 0; i != A.size(); ++i) {
461  if (A[i][i] < 0.0) {
462  return false;
463  }
464  }
465 
466  return true;
467  }
The documentation for this class was generated from the following file:
JLANG::JSingleton::getInstance
static data_type & getInstance()
Get unique instance of template class.
Definition: JSingleton.hh:27
std::vector
Definition: JSTDTypes.hh:12
JTOOLS::j
int j
Definition: JPolint.hh:634
JMATH::JMatrixND_t::col_type
std::vector< double >::iterator col_type
Definition: JMatrixND.hh:50
distance
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Definition: PhysicsEvent.hh:434
JMATH::JMatrixND_t::row_type
matrix_type::iterator row_type
Definition: JMatrixND.hh:49
THROW
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:669
JMATH::JMatrixND_t::const_row_type
matrix_type::const_iterator const_row_type
Definition: JMatrixND.hh:46
JMATH::JMatrixND_t::const_col_type
std::vector< double >::const_iterator const_col_type
Definition: JMatrixND.hh:47
JMATH::JMatrixND::JMatrixND
JMatrixND()
Default constructor.
Definition: JMatrixND.hh:76
JMATH::JMatrixND::at
double at(size_t row, size_t col) const
Get matrix element.
Definition: JMatrixND.hh:120
JMATH::JCalculator
Auxiliary class for arithmetic operations on objects.
Definition: JCalculator.hh:18
JMATH::JMatrixND::resize
void resize(const size_t size)
Resize matrix.
Definition: JMatrixND.hh:98
Generated by   doxygen 1.8.16