Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JMatrixND.hh
Go to the documentation of this file.
1 #ifndef __JMATH__JMATRIXND__
2 #define __JMATH__JMATRIXND__
3 
4 #include <limits>
5 #include <cmath>
6 #include <algorithm>
7 #include <ostream>
8 #include <iomanip>
9 #include <vector>
10 #include <utility>
11 
12 #include "JIO/JSerialisable.hh"
13 #include "JLang/JEquals.hh"
14 #include "JLang/JException.hh"
15 #include "JLang/JSingleton.hh"
16 #include "Jeep/JPrint.hh"
17 #include "JMath/JMath.hh"
18 
19 
20 /**
21  * \author mdejong
22  */
23 
24 namespace JMATH {}
25 namespace JPP { using namespace JMATH; }
26 
27 namespace JMATH {
28 
29  using JIO::JReader;
30  using JIO::JWriter;
31  using JLANG::JEquals;
33  using JLANG::JException;
34  using JLANG::JSingleton;
35 
36 
37  /**
38  * Type definition of a NxN matrix.
39  */
40  struct JMatrixND_t :
41  public std::vector< std::vector<double> >
42  {
43 
45 
46  typedef matrix_type::const_iterator const_row_type;
48 
49  typedef matrix_type::iterator row_type;
51  };
52 
53 
54  /**
55  * Type definition of permutation matrix.
56  */
58 
59 
60  /**
61  * N x N matrix.
62  */
63  class JMatrixND :
64  public JMatrixND_t,
65  public JMath <JMatrixND>,
66  public JEquals <JMatrixND>
67  {
68  public:
69 
71 
72 
73  /**
74  * Default constructor.
75  */
77  {}
78 
79 
80  /**
81  * Constructor.
82  *
83  * \param size dimension
84  */
85  JMatrixND(const size_t size)
86  {
87  resize(size);
88  reset();
89  }
90 
91 
92  /**
93  * Resize matrix.
94  * Note that this method resets the matrix.
95  *
96  * \param size dimension
97  */
98  void resize(const size_t size)
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  }
111 
112 
113  /**
114  * Get matrix element.
115  *
116  * \param row row number
117  * \param col column number
118  * \return matrix element at (row,col)
119  */
120  double at(size_t row, size_t col) const
121  {
122  return (*this)[row][col];
123  }
124 
125 
126  /**
127  * Get matrix element.
128  *
129  * \param row row number
130  * \param col column number
131  * \return matrix element at (row,col)
132  */
133  double& at(size_t row, size_t col)
134  {
135  return (*this)[row][col];
136  }
137 
138 
139  /**
140  * Get matrix element.
141  *
142  * \param row row number
143  * \param col column number
144  * \return matrix element at (row,col)
145  */
146  double operator()(size_t row, size_t col) const
147  {
148  return at(row, col);
149  }
150 
151 
152  /**
153  * Get matrix element.
154  *
155  * \param row row number
156  * \param col column number
157  * \return matrix element at (row,col)
158  */
159  double& operator()(size_t row, size_t col)
160  {
161  return at(row, col);
162  }
163 
164 
165  /**
166  * Set matrix to the null matrix.
167  *
168  * \return this matrix
169  */
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  }
180 
181 
182  /**
183  * Set to identity matrix
184  *
185  * \return this matrix
186  */
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  }
200 
201 
202  /**
203  * Transpose.
204  *
205  * \return this matrix
206  */
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  }
217 
218 
219  /**
220  * Negate matrix.
221  *
222  * \return this matrix
223  */
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  }
234 
235 
236  /**
237  * Matrix addition.
238  *
239  * \param A matrix
240  * \return this matrix
241  */
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  }
258 
259 
260  /**
261  * Matrix subtraction.
262  *
263  * \param A matrix
264  * \return this matrix
265  */
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  }
282 
283 
284  /**
285  * Scale matrix.
286  *
287  * \param factor factor
288  * \return this matrix
289  */
290  JMatrixND& mul(const double factor)
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  }
300 
301 
302  /**
303  * Scale matrix.
304  *
305  * \param factor factor
306  * \return this matrix
307  */
308  JMatrixND& div(const double factor)
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  }
318 
319 
320  /**
321  * Matrix multiplication.
322  *
323  * \param A matrix
324  * \param B matrix
325  * \return this matrix
326  */
328  const JMatrixND& B)
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  }
352 
353 
354  /**
355  * Equality.
356  *
357  * \param A matrix
358  * \param eps numerical precision
359  * \return true if matrices identical; else false
360  */
361  bool equals(const JMatrixND& A,
362  const double eps = std::numeric_limits<double>::min()) const
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  }
380 
381 
382  /**
383  * Test identity.
384  *
385  * \param eps numerical precision
386  * \return true if identity matrix; else false
387  */
388  bool isIdentity(const double eps = std::numeric_limits<double>::min()) const
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  }
405 
406 
407  /**
408  * Get determinant of matrix.
409  *
410  * \return determinant of matrix
411  */
412  double getDeterminant() const
413  {
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  }
426 
427 
428  /**
429  * Test positive definiteness.
430  *
431  * \return true if positive definite; else false
432  */
433  bool isPositiveDefinite() const
434  {
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  }
447 
448 
449  /**
450  * Test positive semi-definiteness.
451  *
452  * \return true if positive semi-definite; else false
453  */
454  const bool isPositiveSemiDefinite() const
455  {
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  }
468 
469 
470  /**
471  * Read matrix from input.
472  *
473  * \param in reader
474  * \param matrix matrix
475  * \return reader
476  */
477  friend inline JReader& operator>>(JReader& in, JMatrixND& matrix)
478  {
479  int size;
480 
481  in >> size;
482 
483  matrix.resize(size);
484 
485  for (row_type row = matrix.begin(); row != matrix.end(); ++row) {
486  for (col_type col = row->begin(); col != row->end(); ++col) {
487  in >> *col;
488  }
489  }
490 
491  return in;
492  }
493 
494 
495  /**
496  * Write matrix to output.
497  *
498  * \param out writer
499  * \param matrix matrix
500  * \return writer
501  */
502  friend inline JWriter& operator<<(JWriter& out, const JMatrixND& matrix)
503  {
504  int size = matrix.size();
505 
506  out << size;
507 
508  for (const_row_type row = matrix.begin(); row != matrix.end(); ++row) {
509  for (const_col_type col = row->begin(); col != row->end(); ++col) {
510  out << *col;
511  }
512  }
513 
514  return out;
515  }
516 
517 
518  /**
519  * Print ASCII formatted output.
520  *
521  * \param out output stream
522  * \param A matrix
523  * \return output stream
524  */
525  friend inline std::ostream& operator<<(std::ostream& out, const JMatrixND& A)
526  {
527  using namespace std;
528 
529  JFlags flags(out);
530 
531  for (const_row_type row = A.begin(); row != A.end(); ++row) {
532 
533  for (const_col_type col = row->begin(); col != row->end(); ++col) {
534  out << FIXED(10,2) << *col << ' ';
535  }
536 
537  out << endl;
538  }
539 
540  return out;
541  }
542 
543  protected:
544  /**
545  * Auxiliary class for matrix decomposition.
546  */
547  struct JMatrixHelper :
548  public JMatrixND_t,
549  public JSingleton<JMatrixHelper>
550  {
551  /**
552  * Set matrix helper.
553  *
554  * \param A matrix
555  * \return this matrix helper
556  */
557  const JMatrixHelper& set(const JMatrixND& A)
558  {
559  static_cast<JMatrixND_t&>(*this) = A;
560 
561  return *this;
562  }
563 
564 
565  /**
566  * LU decomposition.
567  *
568  * \param A matrix
569  */
570  void decompose(const JMatrixND& A)
571  {
572  set(A);
573 
574  sign = +1;
575 
576  size_t k;
577  double val, del;
578  row_type row, root;
579  col_type col, cursor;
580 
581  // LU decomposition
582 
583  for (root = this->begin(), k = 0; root != this->end(); ++root, ++k) {
584 
585  row_type pivot;
586 
587  val = (*root)[k];
588  pivot = root;
589 
590  for (row = root; ++row != this->end(); ) {
591  if (fabs((*row)[k]) > fabs(val)) {
592  val = (*row)[k];
593  pivot = row;
594  }
595  }
596 
597  if (val == 0) {
598  throw JDivisionByZero("LU decomposition zero pivot");
599  }
600 
601  if (pivot != root) {
602 
603  iter_swap(root, pivot);
604 
605  sign = -sign;
606  }
607 
608 
609  // decompose in LU
610 
611  for (row = root; ++row != this->end(); ) {
612 
613  del = (*row)[k] /= val;
614 
615  for (col = row->begin() + k, cursor = root->begin() + k; ++col != row->end(); )
616  (*col) -= del * (*(++cursor));
617  }
618  }
619  }
620 
621 
622  /**
623  * Get sign corresponding to permutation matrix.
624  */
625  int getSign() const
626  {
627  return sign;
628  }
629 
630  private:
631  int sign;
632  };
633  };
634 
635 
636  /**
637  * Determine dot product <tt>Y^T x V x Y</tt>, where
638  * - <tt>Y</tt> is a <tt>1xN</tt> matrix; and
639  * - <tt>V</tt> is a <tt>NxN</tt> matrix.
640  *
641  * \param __begin begin of data
642  * \param __end end of data
643  * \param V matrix
644  * \return dot product
645  */
646  template<class T>
647  inline double getDot(T __begin,
648  T __end,
649  const JMatrixND& V)
650  {
651  double dot = 0.0;
652 
653  T yl = __begin;
654 
655  for (JMatrixND::const_row_type row = V.begin(); row != V.end(); ++row, ++yl) {
656 
657  T yr = __begin;
658 
659  for (JMatrixND::const_col_type col = row->begin(); col != row->end(); ++col, ++yr) {
660  dot += (*yl) * (*col) * (*yr);
661  }
662  }
663 
664  return dot;
665  }
666 }
667 
668 #endif
double & at(size_t row, size_t col)
Get matrix element.
Definition: JMatrixND.hh:133
General exception.
Definition: JException.hh:40
double at(size_t row, size_t col) const
Get matrix element.
Definition: JMatrixND.hh:120
Exceptions.
Simple singleton class.
Definition: JSingleton.hh:18
Interface for binary output.
friend JReader & operator>>(JReader &in, JMatrixND &matrix)
Read matrix from input.
Definition: JMatrixND.hh:477
Auxiliary base class for aritmetic operations of derived class types.
Definition: JMath.hh:26
std::vector< double >::iterator col_type
Definition: JMatrixND.hh:50
int getSign() const
Get sign corresponding to permutation matrix.
Definition: JMatrixND.hh:625
double getDot(const JNeutrinoDirection &first, const JNeutrinoDirection &second)
Dot product.
Definition: JAstronomy.hh:409
JMatrixND & transpose()
Transpose.
Definition: JMatrixND.hh:207
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:633
Auxiliary data structure for floating point format specification.
Definition: JPrint.hh:461
const bool isPositiveSemiDefinite() const
Test positive semi-definiteness.
Definition: JMatrixND.hh:454
bool isPositiveDefinite() const
Test positive definiteness.
Definition: JMatrixND.hh:433
const JMatrixHelper & set(const JMatrixND &A)
Set matrix helper.
Definition: JMatrixND.hh:557
void resize(const size_t size)
Resize matrix.
Definition: JMatrixND.hh:98
JMatrixND & add(const JMatrixND &A)
Matrix addition.
Definition: JMatrixND.hh:242
std::vector< std::pair< int, int > > JPermutationMatrix
Type definition of permutation matrix.
Definition: JMatrixND.hh:57
I/O formatting auxiliaries.
matrix_type::const_iterator const_row_type
Definition: JMatrixND.hh:46
JMatrixND & setIdentity()
Set to identity matrix.
Definition: JMatrixND.hh:187
JMatrixND & reset()
Set matrix to the null matrix.
Definition: JMatrixND.hh:170
matrix_type::iterator row_type
Definition: JMatrixND.hh:49
std::vector< double >::const_iterator const_col_type
Definition: JMatrixND.hh:47
Template definition of auxiliary base class for comparison of data structures.
Definition: JEquals.hh:24
N x N matrix.
Definition: JMatrixND.hh:63
friend std::ostream & operator<<(std::ostream &out, const JMatrixND &A)
Print ASCII formatted output.
Definition: JMatrixND.hh:525
friend JWriter & operator<<(JWriter &out, const JMatrixND &matrix)
Write matrix to output.
Definition: JMatrixND.hh:502
Interface for binary input.
JMatrixND(const size_t size)
Constructor.
Definition: JMatrixND.hh:85
JMatrixND & negate()
Negate matrix.
Definition: JMatrixND.hh:224
Auxiliary class for matrix decomposition.
Definition: JMatrixND.hh:547
void decompose(const JMatrixND &A)
LU decomposition.
Definition: JMatrixND.hh:570
bool equals(const JMatrixND &A, const double eps=std::numeric_limits< double >::min()) const
Equality.
Definition: JMatrixND.hh:361
Exception for division by zero.
Definition: JException.hh:252
static data_type & getInstance()
Get unique instance of template class.
Definition: JSingleton.hh:27
Base class for data structures with artithmetic capabilities.
JMatrixND & sub(const JMatrixND &A)
Matrix subtraction.
Definition: JMatrixND.hh:266
double getDeterminant() const
Get determinant of matrix.
Definition: JMatrixND.hh:412
JMatrixND & mul(const JMatrixND &A, const JMatrixND &B)
Matrix multiplication.
Definition: JMatrixND.hh:327
JMatrixND & div(const double factor)
Scale matrix.
Definition: JMatrixND.hh:308
JMatrixND & mul(const double factor)
Scale matrix.
Definition: JMatrixND.hh:290
bool isIdentity(const double eps=std::numeric_limits< double >::min()) const
Test identity.
Definition: JMatrixND.hh:388
Type definition of a NxN matrix.
Definition: JMatrixND.hh:40
Auxiliary class to temporarily modify format specifications.
Definition: JPrint.hh:535
double operator()(size_t row, size_t col) const
Get matrix element.
Definition: JMatrixND.hh:146
JMatrixND()
Default constructor.
Definition: JMatrixND.hh:76
std::vector< std::vector< double > > matrix_type
Definition: JMatrixND.hh:44
double & operator()(size_t row, size_t col)
Get matrix element.
Definition: JMatrixND.hh:159