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...
 
JMatrixNDreset ()
 Set matrix to the null matrix. More...
 
JMatrixNDsetIdentity ()
 Set to identity matrix. More...
 
JMatrixNDtranspose ()
 Transpose. More...
 
JMatrixNDnegate ()
 Negate matrix. More...
 
JMatrixNDadd (const JMatrixND &A)
 Matrix addition. More...
 
JMatrixNDsub (const JMatrixND &A)
 Matrix subtraction. More...
 
JMatrixNDmul (const double factor)
 Scale matrix. More...
 
JMatrixNDmul (const JMatrixND &A, const JMatrixND &B)
 Matrix multiplication. More...
 
JMatrixNDmul (const JSecond_t &object)
 Multiply with object. More...
 
JMatrixNDdiv (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

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

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 
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